### Poisson+internal Neumann: Why I am getting NaN and Inf over all cells ?

109
views
0
12 weeks ago by
Hi there

I am trying to solve a nonlinear Poisson equation by Picard iteration,  However I am getting no errors but nonsense values from the first round, I cant see whether I have made a mistake or should I play with the solver options,

I have tried to summarize the problem into a short format, I hope it is clear enough

Cheers


from dolfin import *
import numpy as np

#-----------------------------------------------
eps0 = 8.85E-12
q = 1.60217657E-19
KB = 1.3806488E-23
T =300
REV = 0
Q_interface =  2.12E-2

#-----------------------------------------------
width_A = 2E-9
width_B = 2E-9
x_A   = 0
x_B = width_A +x_A
x_end   = width_B + x_B

mesh = IntervalMesh( int(abs(x_end-x_A)/.5E-10) , x_A , x_end)
V = FunctionSpace( mesh, "DG" , 0)
u , v = TrialFunction(V) , TestFunction(V)
#-----------------------------------------------
tol = 1E-14  ;
class Subdomain_1(SubDomain):
def inside(self , x , on_boundary):
return    (x[0] >= x_A - tol )   and  (x[0] <= x_B + tol)

class Subdomain_2(SubDomain):
def inside(self , x , on_boundary):
return     (x[0] >= x_B - tol )   and   (x[0] <= x_end + tol)

subdomains = CellFunction("size_t" , mesh)
subdomains.set_all(0)
Subdomain_1().mark(subdomains,1)
Subdomain_2().mark(subdomains,2)

boundary_markers = FacetFunction("size_t" , mesh)
boundary_markers.set_all(0)
Neuman_boundary= CompiledSubDomain(  "fabs(x[0]-"+ str(x_B) +")<DOLFIN_EPS" )
Neuman_boundary.mark(boundary_markers,1)
Dirich_boundary= CompiledSubDomain(  "fabs(x[0]-"+ str(x_end) +")<DOLFIN_EPS" )
Dirich_boundary.mark(boundary_markers,2)
bc = DirichletBC(V , Constant(REV) , boundary_markers , 2 )

dx =Measure('dx',domain=mesh,subdomain_data=subdomains)
ds =Measure('ds',domain=mesh,subdomain_data=boundary_markers)
dS =Measure('dS',domain=mesh,subdomain_data=boundary_markers)

#-----------------------------------------------
dielectric_v = eps0* np.array([10.4,10.34 ] )
eps  = Function(V)
for cell_no in range (len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
eps.vector()[cell_no] =dielectric_v[subdomain_no-1]

nominal_doping_profile = np.array([1E14 ,2E14 ] )
Nd  = Function(V)
for cell_no in range (len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
Nd.vector()[cell_no] =nominal_doping_profile[subdomain_no-1]
#-----------------------------------------------

L =  q*(Nd )*v * dx(1) \
+q*(Nd )*v * dx(2) \
-avg(Q_interface)*v('+')*dS(1)

u= Function(V)
solve(a==L,u,bc)

​
Community: FEniCS Project

1
12 weeks ago by
You define DG0 function space, and then your energy involves a gradient. Did you try with CG1 for instance?
0
12 weeks ago by
Hi Martin,
As you already mentioned,I had made a mistake using DG0 and using CG1 resolved that specific problem.