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:
- is there an issue with my formulation
- What’s wrong with the code?
Thank you all for the patience over what’s a basic question!