Hi Christophe, I tried running your problem, but I get a seg fault with both the PySparse and Trilinos solvers. The Scipy solver seems to hang.
Which solvers are you using? The link to the version of your code that I'm using is below. Could you check that it is as you intended? https://gist.github.com/wd15/5c780ffbbae949cbf04eee5720e5a727 In the code that was in the email there is an indentation issue with the lines beginning "E.constrain" and "J.constrain" so I'm wondering if there are any other issues. Cheers, Daniel On Sat, Apr 1, 2017 at 12:09 PM, Christoph Günther <christophgn...@gmail.com> wrote: > Dear Mailing list, > > I have been trying to solve a Nonlinear Envelope Equation, discriping the > propagation of a laser pulse in a nonlinear optical medium within Fipy. The > Electric field couples to the electron subsystem through the Term posed Rest > in the following and subsequently raises the electron density in the regions > where the electric field is the strongest (within the focus of the > laserbeam). The Property that I am interested in finally is the electron > density. The system of Equations goes as follows: > > dEdz + dEdt = i/2k * d2Ed2r - i *k''/2 dE2dt2 + i*kn''|E|^2 * E + > Rest(|E|^2, |E|, E) > dNdt = Wpi(|E| + N*sigma*|E|^2 > > Basically I would like to reproduce the attached graph from the attached > paper. As I am a FiPy beginner I wonder if such a Problem can actually be > takled with fipy or if speciallized solving methods such as the split step > fourier method (see e.g. https://en.wikipedia.org/wiki/Split-step_method) > are required. > > I tryed to incorporate the Problem into Fipy by spliting the equation into > its real and complex parts and solving for the coupled system of equations. > The problem that I am facing is that the incident pulse which I define with > a time dependent boundary condition is appearing at the boundary, but it > does not propagate through the medium. I am wondering if the E.grad[1] > operator, which I used to express dEdz is actually a correct representation > for this term or if it could be the source for the incorrect behaviour. This > assignement came from the observation that if I pose the problem in the > following way: > > dEdz = i/2k * d2Ed2r - i *k''/2 dE2dt2 + i*kn''|E|^2 * E + R(|E|^2, |E|, > E) - dEd t > > I get the following error message: > > File > "C:\Python27\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", > line 685, in runfile > execfile(filename, namespace) > > File > "C:\Python27\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", > line 71, in execfile > exec(compile(scripttext, filename, 'exec'), glob, loc) > > File "D:/Arbeit/Masterarbeit/Programmiertes/DGL/NLSE_basic.py", line 129, > in <module> > coupledeqn = eqn1 & eqn2 & eqn3 & eqn4 & eqn5 > > File "C:\Python27\lib\site-packages\fipy\terms\term.py", line 446, in > __and__ > elif other == 0: > > File "C:\Python27\lib\site-packages\fipy\variables\variable.py", line > 1375, in __nonzero__ > return bool(self.value) > > ValueError: The truth value of an array with more than one element is > ambiguous. Use a.any() or a.all() > > In the following you can find the full code. I am hopefull that someone can > spot an error that prevents the Pulse from Propagating through the medium. > Thanks in advance for all your efforts. And thanks to the developers for > setting up this very impressive and capable Python module. > > With Best regards, > Christoph > > > from fipy import * > from scipy.special import erf,ellipk,ellipe > from scipy.integrate import quad > > def thornber(E): > e = 1.602e-19 > kb = 1.380e-23 > vs = 2e5 > Eg = 3.3*e > EI = 1e13 > Eph = 0.8e6 > T = 300 > Ekt = (EI*kb*T)/Eg > > #thorn = > numerix.abs(((vs*e*E)/Eg)*numerix.exp((-1*EI)/((E*(1+(E/Eph)))+Ekt))) > > thorn = ((vs*e*E)/Eg)*numerix.exp((-1*EI)/((E*(1+(E/Eph)))+Ekt)) > return numerix.absolute(thorn) > > def integ(j,x,gamma1,gamma2): > return > erf((numerix.pi/2)*numerix.sqrt(((2.0*(int(x+1.0)))-(2*x)+j)/(ellipk(gamma2)*ellipe(gamma2))))*numerix.exp(-1*j*numerix.pi*((ellipk(gamma2)-ellipe(gamma2))/ellipe(gamma1))) > > def Keld(E): > #E = E*100000000.0 > e = 1.602e-19 > hbar = 1.054e-34 > c = 299792458 > wavelength = 1030e-9 > w = 2*numerix.pi*c/wavelength > me = 0.4*9.109e-31 > mh = 3*9.109e-31 > m = (me*mh)/(me+mh) > Eg = 3.3*1.602e-19 > > if E<0.000000001: > return 0 > else: > gamma = ((w*numerix.sqrt(m*Eg))/(e*E)) > gamma1 = (gamma**2)/(1+gamma**2) > gamma2 = 1.0 - gamma1 > x = > (2/numerix.pi)*(Eg/(hbar*w))*(numerix.sqrt(1+gamma**2)/gamma)*ellipe(1/(1+gamma**2)) > Q = numerix.sqrt((numerix.pi)/(2*ellipk(gamma2)))*quad(lambda > j: integ(j,x,gamma1,gamma2),0,10000)[0] > keld = > ((2*w)/(9*numerix.pi))*(((w*m)/(hbar*numerix.sqrt(gamma1)))**(3/2))*Q*numerix.exp(-1*numerix.pi*(int(x+1))*((ellipk(gamma1)-ellipe(gamma1))/ellipe(gamma2))) > #keld = keld*1000000 > return numerix.absolute(keld) > > def multi(E): > tmp=numerix.empty((),) > for i in range(len(E)): > numerix.append(tmp,Keld(E()[i])) > return tmp > > if __name__ == "__main__": > nr = 50 > nz = 50 > dr = 1e-6 > dz = dr > L = dr*nr > mesh = mesh = CylindricalGrid2D(nr=nr, > nz=nz,dr=dr,dz=dz,origin=((-L/2.0,),(0,)))#Grid2D(dx=dx, dy=dy, nx=nx, > ny=ny) > > #create variables > E = CellVariable(name="real part electric Field", mesh=mesh, > hasOld=True, value=0.00000000000000000000000000000000000000001, ) > Estar = CellVariable(name="diff real part electric Field", mesh=mesh, > hasOld=True, value=0.0 ) > J = CellVariable(name="imag part electric Field", mesh=mesh, > hasOld=True, value=0.00000000000000000000000000000000000000001 ) > Jstar = CellVariable(name="diff imag part electric Field", mesh=mesh, > hasOld=True, value=0.0 ) > rho = CellVariable(name="carrier concentration", mesh=mesh, > hasOld=True, value=1.0e12 ) > time = Variable() > > #make constants for eqn system > #physical constants > c = 299792458.0 > > #material Parameters > n = 2.5 > alpha = 0.05 > n2 = 3.5e-16 > vg = 1.0e18 > Beta2 = 1.0e18 > sigma = 1. > tau = 1e-16 > Eg = 3.3*1.602e-19 > #Laser Parameters > wavelength = 1030e-9 #wavelength in nm > w = 2.0*numerix.pi*c/(wavelength) > d = 25e-6 #Focusdepth > wf = 1.0e-6 #Beam radius in Focus > tp = 1e-9 #Puls length > Ein = 100e-3 #Pulsenergy > timeStep = 1e-11 > desiredResidual = 1.0e-8 > totalElapsedTime = 1000*timeStep > > #physical relations > k = n*w/c > > #equation parsmeters > gamma = numerix.array(((1./(2*k), 0.), (0., 0.))) > > #Laserparameter relations > zf = numerix.pi*(wf**2)*n/wavelength > w0 = wf*numerix.sqrt((1.0+((d**2)/(zf**2)))) > f = d+zf**2.0/d > Pin = Ein/(tp*(numerix.pi/2)) > E0 = numerix.sqrt((2.0*Pin)/(numerix.pi*w0)) > > #set upper values > R,Z = mesh.faceCenters > #mask = ((Y > 19.5e-6)) > E.constrain((E0*numerix.cos((k*(R)**2.0)/(2.0*f)))/(numerix.exp((((R)**2)/(w0**2.0))+((time**2.0)/(tp**2.0)))), > where=mesh.facesTop) > J.constrain((E0*(-1)*numerix.sin((k*(R)**2.0)/(2.0*f)))/(numerix.exp((((R)**2.0)/(w0**2.0))+((time**2)/(tp**2)))), > where=mesh.facesTop) > > #setup equation system > eqn1 = TransientTerm(var=E)==Estar > eqn2 = TransientTerm(var=J)==Jstar > eqn3 = E.grad[1] + TransientTerm(coeff=vg, var=E) == > DiffusionTerm(coeff=(-gamma), var=J) \ > +TransientTerm(coeff= Beta2/2.0,var=Jstar)- > (k*n2*((E**2.0))+((J**2.0))*J)+(w0*tau*rho*J)\ > -((sigma*rho*E)/2)-((multi(numerix.sqrt((E**2.0)+(J**2.0)))*E)/(2*((E**2)+(J**2.0)))) > - alpha*E > eqn4 = J.grad[1] + TransientTerm(coeff=vg, var=J) == > DiffusionTerm(coeff=(gamma), var=E) \ > -TransientTerm(coeff= Beta2/2.0,var=Estar)+ > (k*n2*((E**2))+((J**2))*E)-(w0*tau*rho*E)\ > -((sigma*rho*J)/2)-((multi(numerix.sqrt((E**2.0)+(J**2)))*J)/(2.0*((E**2)+(J**2.0)))) > - alpha*J > eqn5 = TransientTerm(var=rho)== > rho*thornber((E**2.0)+(J**2.0))+multi((E**2.0)+(J**2.0)) > coupledeqn = eqn1 & eqn2 & eqn3 & eqn4 & eqn5 > > #solve equation system > while time < totalElapsedTime: > residual = 1e5 > E.updateOld() > Estar.updateOld() > J.updateOld() > Jstar.updateOld() > rho.updateOld() > #viewer = Viewer(vars=J) > #viewer.plot() > time.setValue(time() + timeStep) > while residual > desiredResidual: > residual = coupledeqn.sweep(dt=timeStep) > > viewer = Viewer(vars=J) > viewer.plot() > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Daniel Wheeler _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
