PETSc problem when solving mixed 3D nonlinear variational problem

General
1

Hello,
We are using fenics for solving the equations below (which describes the electric field inside a pn junction):

\mathrm{\nabla}^2U=\frac{-e_0}{\varepsilon\varepsilon_0}\left(N_D-N_A+p-n\right)
-\frac{k_BT}{e_0}\nabla^2n+\nabla n\cdot \nabla U = 0
\frac{k_BT}{e_0}\nabla^2p+\nabla p\cdot \nabla U = 0

But we encountered an error which said:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.736e+07 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
…
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside ./dolfin/la/PETScMatrix.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Here is the minimal executable code:

import fenics

e0 = 1.60217733e-19
perm_mat = 11.7
perm0 = 8.854187817e-12   #F/m  
kboltz=8.617385e-5 #eV/K
u_T = kboltz * 300/1 # u_T = kT/q
n_i = 6.95e-3 #Silicon, um^-3

mesh3D = fenics.BoxMesh(fenics.Point(0, 0, 0),
                        fenics.Point(1300, 1300, 55), 
                        5, 5, 2000)
P1 = fenics.FiniteElement('P', fenics.tetrahedron, 1)
element = fenics.MixedElement([P1, P1, P1])
V = fenics.FunctionSpace(mesh3D, element)

funcs = fenics.Function(V)
u, p, n = fenics.split(funcs) # Electric Potential, Electron Density, Hole Density
v_u, v_p, v_n = fenics.TestFunctions(V)

Neff = fenics.Expression("x[2] < 1 - tol? 1e6 : (x[2] < 2 - tol ? -5 : (x[2] < 3 + tol ? -1.917e4 : (x[2] < 53 + tol ? -5 : -1e6)))", degree = 0, tol = 1e-14)
Neff_0 = eval("1e6 if (z < 1) else (-5 if (z < 2) else (-1.917e4 if (z < 3) else -5 if (z < 53) else -1e6))".replace("z","0"))
Neff_end = eval("1e6 if (z < 1) else (-5 if (z < 2) else (-1.917e4 if (z < 3) else -5 if (z < 53) else -1e6))".replace("z","55"))

u_D = fenics.Expression('x[2]<tol ? p_1:p_2',
                        degree = 2, tol = 1E-14,
                        p_1 = 200, p_2 = 0)

p_D = fenics.Expression('x[2]<tol ? p_1:p_2',
                        degree = 2, tol = 1E-14,
                        p_1 = 0.5*(-Neff_0 + (Neff_0**2 + 4*n_i**2)**0.5), 
                        p_2 = 0.5*(-Neff_end + (Neff_end**2 + 4*n_i**2)**0.5))
  
n_D = fenics.Expression('x[2]<tol ? p_1:p_2',
                        degree = 2, tol = 1E-14,
                        p_1 = 0.5*(Neff_0 + (Neff_0**2 + 4*n_i**2)**0.5), 
                        p_2 = 0.5*(Neff_end + (Neff_end**2 + 4*n_i**2)**0.5))

def boundary(x, on_boundary):
    return abs(x[2])<1e-14 or abs(x[2]-55)<1e-14
  
u_bc = fenics.DirichletBC(V.sub(0), u_D, boundary)
p_bc = fenics.DirichletBC(V.sub(1), p_D, boundary)
n_bc = fenics.DirichletBC(V.sub(2), n_D, boundary)

g_expression = fenics.Constant(0)

a_u = fenics.dot(fenics.grad(u), fenics.grad(v_u))*fenics.dx
a_p = fenics.dot(fenics.grad(p), fenics.grad(v_p))*u_T*fenics.dx
a_n = fenics.dot(fenics.grad(n), fenics.grad(v_n))*(-u_T)*fenics.dx

L_u = ((Neff+n-p)*1e6*e0/perm0/perm_mat)*v_u*fenics.dx
L_p = (fenics.inner(fenics.grad(u),fenics.grad(p))+g_expression)*v_p*fenics.dx
L_n = (fenics.inner(fenics.grad(u),fenics.grad(n))+g_expression)*v_n*fenics.dx

funcs_ = fenics.Function(V)
fenics.solve(a_u+a_p+a_n - L_u+L_p+L_n == 0 , funcs_, [u_bc, p_bc, n_bc])
2

You are solving a non-linear problem wrt the variable funcs

And thus should call:

fenics.solve(a_u+a_p+a_n - L_u+L_p+L_n == 0 , funcs, [u_bc, p_bc, n_bc])

Because to solve a non-linear problem, one has to compute the jacobian of your residual, and one has to know what variable should be differented wrt.

3

Thank you, the Newton iterator works normally after I applied your code by changing the last two lines into

fenics.solve(a_u+a_p+a_n-L_u-L_p-L_n == 0 , funcs, [u_bc, p_bc, n_bc])

(the equation was wrong before), but a new error appeared:

UMFPACK V5.7.9 (Oct 20, 2019): ERROR: out of memory

*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside ./dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown

Is this a bug or a math problem that requires me to adjust the conditions of the solver?

4

Change the underlying solver of your problem (as UMFPACK, the default LU solver is running out of memory). See for instance: UMFPACK error: out of memory despite system having free memory - #2 by plugged

5

Thank you, and… the solution indeed failed to converge after I changed the solver.

*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = -1.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.

6

Note that non-linear solvers should have suitable initial conditions to ensure convergence. Having a Zero initial guess might be too far away from a physical solution.