Diffusion/recombination problem

variational formulation
1

Hi! I am super new to using this library, so I apologies for the inherent ignorance in the questions! I am trying to solve a simple 1D charge carrier diffusion/recombination problem, and am having a hard time figuring out what’s exactly wrong with my formulation.

The equation I need to solve is:
D \frac{d^2n(x)}{dx^2} = -\alpha n_\gamma e^{-\alpha x} + k_1 n(x) + k_2 n(x)^2+ k_3 n(x)^3
with boundry conditions

S_f n(0) = D \frac{dn(x)}{dx}|_{x=0}; S_b n(d) = -D \frac{dn(x)}{dx}|_{x=d}

This is a nonlinear problem with robin boundry conditions.

I followed the guide, split the variational form for the second derivitive and added the terms for the bc, giving:

F= \int D \frac{dn}{dx} \frac{dv}{dx} dx + \int_{x=0} (S_f n) nv ds - \int_{x=d} (S_b n) nv ds - \int fv dx
where
f = -\alpha n_\gamma e^{-\alpha x} + k_1 n(x) + k_2 n(x)^2+ k_3 n(x)^3

I wrote this code :

from fenics import *
import numpy as np
import sympy as sym

# Parameters
D = 0.007
k1 = 1e5
k2 = 1e-10
k3 = 1e-28
alpha = 1e5
n_g = 1e10
S_f = 10
S_b = 100
d = 500e-9  # domain length

# Create mesh 
nx = 1000
mesh = IntervalMesh(nx, 0, d)
V = FunctionSpace(mesh, 'P', 1)

# Define boundry
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

class BoundaryX0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, self.tol)
class BoundaryXd(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], d, self.tol)
bc_x0 = BoundaryX0()
bc_xd = BoundaryXd()
bc_x0.mark(boundary_markers, 0)
bc_xd.mark(boundary_markers, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# trial and test functions 
n = Function(V)
v = TestFunction(V)

# generation term
x = sym.symbols('x[0]')
G_sym = alpha*n_g*sym.exp(-alpha*x)
G_code = sym.printing.ccode(G_sym)

# assemble terms with robin bc, where s=0
f = k1*n + k2*n**2 + k3*n**3 - Expression(G_code, degree=2)
F = D*dot(grad(n), grad(v))*dx + S_f*n*n*v*ds(0) - S_b*n*n*v*ds(1) - f*v*dx

n = Function(V)
solve(F == 0, n)

# Plot the solution
import matplotlib.pyplot as plt
plot(n)

plt.xlabel('x')
plt.ylabel('n(x)')

plt.show()

And get this error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[23], line 54
     51 F = D*dot(grad(n), grad(v))*dx - f*v*dx
     53 n = Function(V)
---> 54 solve(F == 0, n)
     56 # Plot the solution
     57 import matplotlib.pyplot as plt

File /opt/miniconda3/envs/akash_fenics/lib/python3.12/site-packages/dolfin/fem/solving.py:220, in solve(*args, **kwargs)
    217 # Call variational problem solver if we get an equation (but not a
    218 # tolerance)
    219 elif isinstance(args[0], ufl.classes.Equation):
--> 220     _solve_varproblem(*args, **kwargs)
    222 # Default case, just call the wrapped C++ solve function
    223 else:
    224     if kwargs:

File /opt/miniconda3/envs/akash_fenics/lib/python3.12/site-packages/dolfin/fem/solving.py:266, in _solve_varproblem(*args, **kwargs)
    264 solver = NonlinearVariationalSolver(problem)
    265 solver.parameters.update(solver_parameters)
--> 266 solver.solve()

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 /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1719292225108/work/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  ea4ec9db9a0ab5757694c6a6ff308b135231a4dd
*** -------------------------------------------------------------------------

My questions are:

  1. is there an issue with my formulation
  2. What’s wrong with the code?

Thank you all for the patience over what’s a basic question!

2

You must not overwrite n before solving. Drop the redefinition of n before solve, and you’ll see the code entering in the nonlinear solver.

The bad news is that the nonlinear solver doesn’t converge.

  • some tips are reported in Default absolute tolerance and relative tolerance by @nate
  • out of those you may want to try with a non zero initial guess, if you know a good one
  • seeing that n_g = 1e10, you may also want to make your problem dimensionless to avoid having to deal with big numbers in arithmetic precision.