Hi, I am very new to fenics and I have gone through several examples in order to familiarize myself. I am studying what happens to dielectric samples under electron bombardment. I have implemented my system as far as I could and more or less analogous to the advection-diffusion model in the fenics manual. However, when I run my code the initial absolute residual is huge (like 4e+6) which tells me I did something wrong but I cannot figure out what I am doing wrong. My guess is that my variational formulation is not correct, even though I went through it quite carefully. My code is attached below, any help is welcome. If it is necessary to state my entire system (i.e. all boundary conditions etc.) then let me know as well.
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import numpy as np
obj_path = '/home/guus/Documents/python_objects'
import sys
sys.path.insert(1,obj_path)
from const import constants
from Mobility import Mobility
cons = constants()
#Introduce all types of constants
eps_0 = cons.eps_0
kb = cons.k_b
m0 = cons.m_e
eV = -cons.e
q_e = eV
eps_r = 9.8
s = 1.0 #Surface recombination velocity, this will need some tuning
n_i = 1e5 #Intrinsic carrier density, literature required here
E_g = 9.8*eV #energy gap
E_t = 1.2*eV #average trapping energy
rho = 3.95 #Density of alumina [g/cm^3]
E_b = 300 #primary electron energy
N_I = 1e27 #Number of ionized impurities
N_t = N_I #Number of potential trapping sites
C_n = 1e-20 #Capture cross section of electrons
C_p = 1e-15 #Capture cross section of holes
tau_n = 2e-9 #SRH electron lifetime
tau_p = tau_n #SRH hole lifetime
r = 0.99 #fraction of ionized impurities in the valence band, not sure if that makes sense
T = 290 #Temperature, generally just room temperature
#Defining beam parameters
r_spot = 1e-3 #beam spot radius, I consider beam pulse symmetric in all directions (i.e. sigma_x=sigma_y=sigma_z)
I_0 = 70e-9 #beam current
R_p = 93.4*((1e-3*E_b)**1.45/rho**0.91)*1e-9 #Penetration depth
A = np.divide(I_0,((2*np.pi)**(3/2))*q_e*r_spot*r_spot) #make it all a constant so it looks a little neater
Bx = 2e-9 #Where the beam hits (x-coor)
By = Bx #Where the beam hits (y-coor)
#Defining mobilities (and diffusion coeff.)
mobility = Mobility('Al2O3','Al2O3','100')
mu_ipp = mobility.mu_ip_holes(T, r*N_I)
mu_ipn = mobility.mu_ip_electrons(T,(1-r)*N_I)
mu_apn = mobility.mu_ap_electrons(T)
mu_app = mobility.mu_ap_holes(T)
mu_n = (1/mu_ipn + 1/mu_apn)**(-1)
mu_p = (1/mu_ipp + 1/mu_app)**(-1)
D_n = mu_n*(kb*T/eV)
D_p = mu_p*(kb*T/eV)
#A brief check
print('mu_n = ',mu_n)
print('mu_p = ',mu_p)
print('D_n = ',D_n)
print('D_p = ',D_p)
#Numerical setup
t_final = 1e-9 #stop after 1 ms
num_steps = 100 #number of time steps
dt = t_final/num_steps
tol = 1E-14 #General tolerance used for defining subdomains
#Defining domain on which I intend to solve
p0 = Point(0.0, 0.0, 0.0)
p1 = Point(4e-7, 4e-7, 0.0)
nx = 20
ny = 20
mesh = RectangleMesh(p0, p1, nx, ny, 'crossed')
#Defining the subdomains and everything that comes with it
#Subdomain 0
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 2e-7 + tol
#Subdomain 1
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 2e-7 - tol
#Create a cellfunction
materials = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
#Now we can mark the two subdomains
materials.set_all(0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
#Now we are making sure there are two media defined
class Permittivity(UserExpression): # UserExpression instead of Expression
def __init__(self, markers, **kwargs):
super().__init__(**kwargs) # This part is new!
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 0:
values[0] = Constant(eps_0*eps_r)
else:
values[0] = Constant(eps_0)
def value_shape(self):
#return (1,)
return ()
eps = Permittivity(materials, degree=0)
#Time for the boundary conditions
#Boundary at x = 0
class BoundaryX0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)
#Boundary at x = 1
class BoundaryX1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 4e-7, tol)
#Boundary at y = 0
class BoundaryY0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, tol)
#Boundary at y = 1
class BoundaryY1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 4e-7, tol)
#Boundary at left lower domain
class BoundaryY05_L(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 2e-7, tol) # and near(x[0], 0, tol)
#Boundary at right lower domain
class BoundaryY05_R(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 2e-7, tol) # and near(x[0], 1, tol)
#Define the interface, just where the two domains meet
class Interface(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 2e-7, tol)
#Mark the boundaries
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bx0_L = BoundaryX0()
bx0_R = BoundaryX0()
Int = Interface()
bx0.mark(materials, 2)
bx1.mark(materials, 3)
by0.mark(materials, 4)
by1.mark(materials, 5)
bx0_L.mark(materials, 6)
bx0_R.mark(materials, 7)
Int.mark(materials, 8)
#Make all the quantities as one entity
element = VectorElement('P',triangle,3,dim=3)
#Defining the function space in which to solve the problem
V = FunctionSpace(mesh, element)
#Setting initial condition
ic = Expression(('V0','n0','p0'),degree=3,V0=0,n0=1e-10,p0=1e-10)
u0 = project(ic,V)
V_0, n_0, p_0 = split(u0)
#Define test functions
v_1, v_2, v_3 = TestFunctions(V)
# Define functions for electric field and densities
#Densities and electric field
u = Function(V)
#This the variable to which we write the new time to
u_n = Function(V)
#Splitting it so that we have access to all the components, u_1 = potential; u_2 = electrons; u_3 = holes.
u_1, u_2, u_3 = split(u)
#We have u_n so that we can use those for backward euler method in time
u_n1, u_n2, u_n3 = split(u_n)
#Now that all the boundaries are marked I will define the specific boundary conditions
#The potential will be zero at all boundaries
zero = Constant((0.0))
bc0_pot = DirichletBC(V.sub(0), zero, bx0)
bc1_pot = DirichletBC(V.sub(0), zero, bx1)
bc2_pot = DirichletBC(V.sub(0), zero, by0)
bc3_pot = DirichletBC(V.sub(0), zero, by1)
#Now setting boundary conditions for the charge carriers
n_i = Constant((n_i))
bc0_n = DirichletBC(V.sub(1), n_i, by0)
bc0_p = DirichletBC(V.sub(2), n_i, by0)
bc1_n = DirichletBC(V.sub(1), n_i, bx0_L)
bc1_p = DirichletBC(V.sub(2), n_i, bx0_L)
bc2_n = DirichletBC(V.sub(1), n_i, bx0_R)
bc2_p = DirichletBC(V.sub(2), n_i, bx0_R)
bcs = [bc0_pot,bc1_pot,bc2_pot,bc3_pot, bc0_n,bc0_p,bc1_n,bc1_p,bc2_n,bc2_p]
#File("materials.pvd").write(materials)
################################################################################################################################################
"""
#Do some debugging because I clearly defined the boundary conditions wrong :(
for x in mesh.coordinates():
if bx0.inside(x, True): print('%s is on x = 0' % x)
if bx1.inside(x, True): print('%s is on x = 1' % x)
if by0.inside(x, True): print('%s is on y = 0' % x)
if by1.inside(x, True): print('%s is on y = 1' % x)
# Print the Dirichlet conditions
print('Number of Dirichlet conditions:', len(bcs))
if V.ufl_element().degree() == 1: # P1 elements
d2v = dof_to_vertex_map(V)
coor = mesh.coordinates()
for i, bc in enumerate(bcs):
print('Dirichlet condition %d' % i)
boundary_values = bc.get_boundary_values()
print(boundary_values)
"""
################################################################################################################################################
#u_1 = potential; u_2 = electrons; u_3 = holes
#Most of the setup is done and it is time to plug in the system of equations
#Definition of some constants for variational notation as well as the source terms
h = Constant(dt)
q_e = Constant(q_e)
mu_n = Constant(mu_n)
mu_p = Constant(mu_p)
D_n = Constant(D_n)
D_p = Constant(D_p)
tau_n = Constant(tau_n)
tau_p = Constant(tau_p)
#Define all integration regimes and what not
dx = Measure('dx',mesh,subdomain_data=materials)
ds = Measure('ds',mesh,subdomain_data=materials)
#Define source terms
#SRH recombination term for thermal equilibrium
R = (n_i*n_i - u_2*u_3)/(tau_n*(u_2 + n_i) + tau_p*(u_3*n_i))
#Beam source with 2D Gaussian model
S_n = Expression("A*exp(-(pow(x[0] - Bx, 2) + pow(x[1] - (By-0.3*R), 2)) / r)",A=A, Bx=Bx, By=By, R=R_p,r=r_spot*r_spot, degree=2)
S_p = Expression("1e-5*A*exp(-(pow(x[0] - Bx, 2) + pow(x[1] - (By-0.3*R), 2)) / r)",A=A, Bx=Bx, By=By, R=R_p,r=r_spot*r_spot, degree=2)
#Write it all down in the variational form
#I have added a quick uniform source, this will change but I must first fix my code...
F1 = (eps/q_e)*dot(grad(u_1), grad(v_1))*dx - (u_2-u_3)*v_1*dx
F2 = -mu_n*dot(grad(u_1),grad(v_2))*dx(0) + D_n*dot(grad(u_2),grad(v_2))*dx(0) + (S_n - R - (u_n2-u_2)/h)*v_2*dx(0) + (-mu_n*dot(grad(u_1),grad(v_2)) + D_n*dot(grad(u_2),grad(v_2)))*ds(8)
F3 = mu_n*dot(grad(u_1),grad(v_3))*dx(0) + D_n*dot(grad(u_3),grad(v_3))*dx(0) + (S_p - R - (u_n3-u_3)/h)*v_3*dx(0)
F = F1 + F2 + F3
J = derivative(F,u)
t = 0
solve(F == 0, u, bcs, J=J)
"""
for n in range(num_steps):
t += dt
solve(F == 0, u, bcs, J=J)
u_n.assign(u)
"""