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 ]
