- edited description
higher numerical precision
Please take a look at the code below. For drift-diffusion type problems my solution u is sometimes constant in certain regions of the mesh resulting in a zero gradient. This is problematic because I need to calculate
J = n*u.dx(0)
where n has very big values and u.dx(0) has very small values. This results in J=0 which is unphysical, instead of the more accurate small number I want to calculate. I tried rescaling the problem and such but that did not help. Please see the code below that illustrates the problem. If Doping is too high (1e5) spikes start to appear in the current density Jn. Reduce Doping to 1e3 and the spikes disappear.
I am looking for suggestions on how to implement float128 precision or use a higher than double precision linear algebra library with dolfin. I know PETSc can be compiled with float128. Is it possible to use Epetra with float128? Any suggestion or experience with this type of problems is welcome.
"Author: Chaffra Affouda, copyright 2014"
from dolfin import *
from dolfin.cpp.la import NewtonSolver
parameters['linear_algebra_backend'] = 'PETSc'
class DDProblem(NonlinearProblem):
def __init__(self, F, J, bcs=None, form_compiler_parameters=None, ident_zeros=False):
NonlinearProblem.__init__(self)
self._F = F
self._J = J
self.bcs = bcs
self.form_compiler_parameters = form_compiler_parameters
self.ident_zeros = ident_zeros
def F(self, b, x):
assemble(self._F, tensor=b,
form_compiler_parameters=self.form_compiler_parameters)
#Apply boundary conditions
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self._J, tensor=A,
form_compiler_parameters=self.form_compiler_parameters)
#Apply boundary conditions
for bc in self.bcs:
bc.apply(A)
if self.ident_zeros:
A.ident_zeros()
class DDSolver(NewtonSolver):
def __init__(self, du_max=1.1):
self.du_max = du_max
NewtonSolver.__init__(self)
def update_solution(self, x, dx):
self.damp_potentials(x, dx)
NewtonSolver.update_solution(self, x, dx)
def damp_potentials(self, x, dx):
"""This will strongly influence the convergence"""
du = dx.array()
dumax = self.du_max
i = (du > dumax)
du[i] = dumax#*(1 + np.log(du[i])/dumax)
i = (du < -dumax)
du[i] = -dumax#*(1 + np.log(-du[i])/dumax)
dx[:] = du
t1 = 1
t2 = 1
t3 = 0.1
x0 = 0.0
x1 = x0 + t1
x2 = x1+t2
x3 = x2+t3
mesh = IntervalMesh(100,x0,x2)
hmax = mesh.hmax()
def inside_region1(x, on_boundary):
return between(x[0],(x0,x1+hmax))
def inside_region2(x, on_boundary):
return between(x[0],(x1,x2+hmax))
cell_domains = CellFunction('size_t',mesh)
region1 = AutoSubDomain(inside_region1)
region2 = AutoSubDomain(inside_region2)
region1.mark(cell_domains,1)
region2.mark(cell_domains,2)
Q = FunctionSpace(mesh,'CG',1)
DG1 = FunctionSpace(mesh,'DG',1)
W = MixedFunctionSpace([Q,Q,Q])
class Doping(Expression):
def eval(self, v, x):
if inside_region1(x,False):
v[0] = -1e5
elif inside_region2(x,False):
v[0] = 1e5
class Ni(Expression):
def eval(self, v, x):
if inside_region1(x,False):
v[0] = 1e-5
elif inside_region2(x,False):
v[0] = 1e-5
u_mixed = Function(W)
v_mixed = TestFunction(W)
du_mixed = TrialFunction(W)
u, un, up = split(u_mixed)
v, vn, vp = split(v_mixed)
du, dn, dp = split(du_mixed)
epsilon = Constant(1.)
T = Constant(1.)
ni = Ni() #Constant(1.0) #Ni()
d = Doping()
n = ni*exp((u-un)/T)
p = ni*exp((up-u)/T)
q = 1.0
rho = (p+d-n)*q
mun = Constant(100.0)
mup = Constant(100.0)
tau_n = Constant(1.0)
tau_p = Constant(1.0)
R = (n*p - ni**2) + (n*p-ni**2)/(tau_n*(p+ni**2)+tau_p*(n+ni**2))
dx = Measure("dx")[cell_domains]
dS = Measure("dS")
L1 = inner(epsilon*grad(u),grad(v))*dx(1) - rho*v*dx(1)
L1 += inner(epsilon*grad(u),grad(v))*dx(2) - rho*v*dx(2)
L2 = inner(mun*n*grad(un), grad(vn))*dx(1) - R*vn*dx(1)
L2 += inner(mun*n*grad(un), grad(vn))*dx(2) - R*vn*dx(2)
L3 = inner(mup*p*grad(up), grad(vp))*dx(1) + R*vp*dx(1)
L3 += inner(mup*p*grad(up), grad(vp))*dx(2) + R*vp*dx(2)
F = L1+L2+L3
def inside_top_contact(x, on_boundary):
return near(x[0],x2,)
def inside_bottom_contact(x, on_boundary):
return near(x[0],x0,)
top_contact = AutoSubDomain(inside_top_contact)
bottom_contact = AutoSubDomain(inside_bottom_contact)
bias = Constant(1.0)
J = derivative(F,u_mixed,du_mixed)
F0 = L1
J0 = derivative(F0,u,du)
F = Form(F)
J = Form(J)
solver = DDSolver()
solver.parameters['maximum_iterations'] = 100
solver.parameters['absolute_tolerance'] = 1e-10
solver.parameters['relative_tolerance'] = 1e-30
solver.parameters['relaxation_parameter'] = 1.0
solver.parameters['error_on_nonconvergence'] = False
solver.parameters['convergence_criterion'] = 'incremental'
solver.parameters['linear_solver'] = 'lu'
u_mixed.vector()[:] = 0.0
problem_phi = DDProblem(F0, J0, bcs=[], ident_zeros=True)
solver.solve(problem_phi, u_mixed.vector())
plot(u, title='equilibrium potential',
interactive=False)
plt = plot(ln(abs(n))/ln(10), title="equilibrium electron density", interactive=False, lc='1,0,0')
#bc_bottom_phi = DirichletBC(W.sub(0),u0,bottom_contact)
bc_bottom_n = DirichletBC(W.sub(1),0.0,bottom_contact)
bc_bottom_p = DirichletBC(W.sub(2),0.0,bottom_contact)
#bc_top_phi = DirichletBC(W.sub(0),u0+bias,top_contact)
bc_top_n = DirichletBC(W.sub(1),bias,top_contact)
bc_top_p = DirichletBC(W.sub(2),bias,top_contact)
bcs_phi = []
bcs_n = [bc_bottom_n, bc_top_n]
bcs_p = [bc_bottom_p, bc_top_p]
problem_dd = DDProblem(F,J,bcs=bcs_phi+bcs_n+bcs_p, ident_zeros=False)
solver.parameters['relaxation_parameter'] = 1.0
solver.du_max = 1.0
solver.solve(problem_dd, u_mixed.vector())
Jn = (-q*mun*n*un.dx(0))
#Jp = (-q*mup*p*up.dx(0))
#f = project(Jn+Jp,DG1)
#plot(f,title='total current')
plt= plot(ln(abs(Jn)+1e-10)/ln(10),title='Jn', lc='1,0,0')
plt = plot(un, title='un',lc='1,0,0')
Comments (4)
-
reporter -
reporter - edited description
-
You could try setting
parameters['form_compiler']['representation'] = 'quadrature'
as the 'tensor' representation in ffc in some cases generate numerically unstable code.
-
- changed status to wontfix
Thread suggests how to increase FFC precision. DOLFIN won't support 128 bit floats any time soon,
- Log in to comment