Python file for running Christoph Gunther's problem

christoph.py

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()