Hi Christoph,
I have it running now and it seems to be running fine with no issues.
I don't get a "ValueError". Is this the code that give the
"ValueError" or is that with a modification of this code?
Anyway, I have the entire environment in a Dockerfile that you're
welcome to use, I'm using the following
https://github.com/wd15/extremefill2D-dockerize/blob/master/base/Dockerfile
if you want to check how your environment might differ from mine. To
run your problem I use
$ docker pull docker.io/wd15/extremefill2d_base
$ docker run -i -t wd15/extremefill2d_base:latest /bin/bash
# cd ..
# wget
https://gist.githubusercontent.com/wd15/5c780ffbbae949cbf04eee5720e5a727/raw/27dfe73dc0e2370c3db04fa46333c406bc5606cd/christoph.py
# python christoph.py
Obviously, you'll need to install Docker to use this.
Cheers,
Daniel
On Thu, Apr 13, 2017 at 4:18 AM, Christoph Günther
<[email protected]> wrote:
> Hey Daniel,
>
> thanks for looking into this. I am sorry for the identation error. I
> thoroughly checked the code you put up on Github and it is now in the way
> that I intended it to be.
> I also downloaded the code and do not get any error messages when running
> the code on my PC. I use the pysparse solver to run the Program:
> print DefaultSolver
> <class 'fipy.solvers.pysparse.linearPCGSolver.LinearPCGSolver'>
>
> I use the Python Version:
> sys.version: 2.7.10 (default, May 23 2015, 09:40:32) [MSC v.1500 32 bit
> (Intel)]
>
> and the Package versions:
> pysparse.__version__ Out[35]: '1.3'
> fipy.__version__ Out[33]: '3.1.3'
> scipy.__version__ Out[31]: '0.15.1'
>
> Does the program give a seg error right away or does it loop through the
> solution a couple of times? If you uncomment the lines 132 and 133 you
> should get a picture with each loop step that is similar to the one I
> attach.
>
> Thanks for your efforts,
> Christoph
>
>
>
> Am 12.04.2017 um 16:24 schrieb Daniel Wheeler:
>>
>> 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
>> <[email protected]> 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
>>> [email protected]
>>> http://www.ctcms.nist.gov/fipy
>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>>
>>
>>
>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
