Hello, everyone. Currently, I’m facing an issue with the mixed simulation of the DG and BDM methods for Equation (1). My current results do not converge. For the electric potential \psi, electron concentration n, and hole concentration p, I’m using DG elements. For the electric field \mathbf{E}, electron current density \mathbf{J}_{n}, and hole current density \mathbf{J}_{p}, I’m using BDM elements. I suspect there might be a problem with my variational derivation, especially regarding the boundary conditions and the interior terms.
For the Dirichlet (D - type) boundaries, I use the Nitsche’s method for implementation. The code is as follows
def nitsche(u, v, g, coeff, n_facet, h, alpha, ds_k):
"""Nitsche method for weakly enforcing Dirichlet boundary conditions"""
return (
-inner(coeff * grad(u), n_facet * v) * ds_k
- inner(coeff * grad(v), n_facet * u) * ds_k
+ alpha / h * inner(u, v) * ds_k
+ inner(coeff * grad(v), n_facet * g) * ds_k
- alpha / h * inner(g, v) * ds_k
)
F1 += nitsche(u_psi, v_psi, psi_anode, lam2, n, h, alpha_psi, ds(1))
F1 += nitsche(u_psi, v_psi, psi_cathode, lam2, n, h, alpha_psi, ds(2))
For the continuity of the internal interfaces, it is directly defined as
def standard_ip(u, v, n, h_avg, gamma):
"""Standard interior penalty method"""
F = -inner(jump(v, n), ufl.avg(grad(u))) * dS
F += -inner(ufl.avg(grad(v)), jump(u, n)) * dS
F += +gamma / h_avg * inner(jump(v, n), jump(u, n)) * dS
return F
F2 += standard_ip(u_psi, v_psi, n, h_avg, alpha_n)
So was for electron concentration n, and hole concentration p
Details :
System of Equations
The system of equations is given by:
\begin{cases}
\lambda^{2}\nabla\cdot\mathbf{E}+n - p - N = 0\\
-\nabla\cdot\mathbf{J}_{n}+R_{n}(n,p)=0\\
\nabla\cdot\mathbf{J}_{p}+R_{p}(n,p)=0\\
\mathbf{E}+\nabla\psi = 0\\
\mathbf{J}_{n}-\mu_{n}n\mathbf{E}-D_{n}\nabla n = 0\\
\mathbf{J}_{p}-\mu_{p}p\mathbf{E}+D_{p}\nabla p = 0
\end{cases} \tag{1}
- First equation: \lambda^{2}\nabla\cdot\mathbf{E}+n - p - N = 0. Here, \lambda is likely a characteristic length scale, \nabla\cdot\mathbf{E} represents the divergence of the electric field \mathbf{E}. The terms n and p denote the electron and hole concentrations respectively, and N is the doping concentration. This equation is related to the charge - neutrality condition in a semiconductor, adjusted by the divergence of the electric field term which accounts for electrostatic effects.
- Second equation: -\nabla\cdot\mathbf{J}_{n}+R_{n}(n,p)=0. \mathbf{J}_{n} is the electron current density. The divergence term \nabla\cdot\mathbf{J}_{n} describes the spatial variation of the electron current. R_{n}(n,p) is a recombination - generation term for electrons, which depends on the electron and hole concentrations n and p. This equation enforces the conservation of electrons, accounting for the changes in electron current density and the recombination - generation processes.
- Third equation: \nabla\cdot\mathbf{J}_{p}+R_{p}(n,p)=0. \mathbf{J}_{p} is the hole current density. Similar to the electron - related equation, the divergence term \nabla\cdot\mathbf{J}_{p} describes the spatial variation of the hole current, and R_{p}(n,p) is the recombination - generation term for holes. This equation enforces the conservation of holes.
- Fourth equation: \mathbf{E}+\nabla\psi = 0. This is the electrostatic relation where \mathbf{E} is the electric field and \psi is the electric potential. It shows that the electric field is the negative gradient of the electric potential, which is a fundamental relationship in electrostatics.
- Fifth equation: \mathbf{J}_{n}-\mu_{n}n\mathbf{E}-D_{n}\nabla n = 0. Here, \mu_{n} is the electron mobility, which describes how easily electrons can move in the presence of an electric field, and D_{n} is the electron diffusion coefficient. This equation combines the drift (due to the electric field \mathbf{E}) and diffusion (due to the concentration gradient \nabla n) components of the electron current density.
- Sixth equation: \mathbf{J}_{p}-\mu_{p}p\mathbf{E}+D_{p}\nabla p = 0. \mu_{p} is the hole mobility and D_{p} is the hole diffusion coefficient. Similar to the electron - related equation, it combines the drift and diffusion components of the hole current density.
Boundary Conditions
Dirichlet Boundary Conditions
For Dirichlet boundaries (under Maxwell - Boltzmann statistics), the boundary conditions are given by:
\begin{align*}
\left.\psi\right|_{\partial\Omega_{\text{Dirichlet}}}&=V_{\text{applied}}+V_{T}\ln\left(\frac{n}{n_{\text{ie}}}\right)\\
\left.n\right|_{\partial\Omega_{\text{Dirichlet}}}&=\frac{1}{2}\left(N+\sqrt{N^{2}+4n_{\text{ie}}^{2}}\right)\\
\left.p\right|_{\partial\Omega_{\text{Dirichlet}}}&=\frac{1}{2}\left(-N+\sqrt{N^{2}+4n_{\text{ie}}^{2}}\right)
\end{align*}
where V_{\text{applied}} is the applied external voltage, which may vary across different electrodes, and V_{T} is the thermal voltage. n_{\text{ie}} is an intrinsic - like electron concentration. These conditions specify the values of the electric potential \psi, electron concentration n, and hole concentration p at the Dirichlet boundary.
Neumann Boundary Conditions
For Neumann boundaries, no - flux conditions are assumed:
\begin{align*}
\left.\nabla\psi\right|_{\partial\Omega_{\text{Neumann}}}\cdot\mathbf{n}&=0\\
\left.\mathbf{J}_{n}\right|_{\partial\Omega_{\text{Neumann}}}\cdot\mathbf{n}&=0\\
\left.\mathbf{J}_{p}\right|_{\partial\Omega_{\text{Neumann}}}\cdot\mathbf{n}&=0
\end{align*}
where \mathbf{n} is the unit outer - normal vector of \partial\Omega_{\text{Neumann}}. These conditions imply that there is no flux of the electric - potential gradient, electron current density, and hole current density across the Neumann boundary, i.e., there is no flow of these quantities through the boundary.
import numpy as np
import matplotlib.pyplot as plt
from petsc4py import PETSc
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities, locate_entities_boundary, meshtags
import dolfinx
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, TestFunctions,
grad, inner, div, exp, ln, conditional, lt, gt, jump)
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import basix
# Use small value instead of zero to avoid numerical issues
EPS = np.finfo(float).eps
##############################
### 1. Parameter Setup - PIN Diode ###
##############################
# Geometric parameters (extracted from Figure 14 in referenced paper)
L_DIODE = 120e-4 # Total diode length [cm]
H_DIODE = 10e-4 # Diode height [cm]
L_P = 4e-4 # P-region length [cm]
L_N = 4e-4 # N-region length [cm]
L_I = L_DIODE - L_P - L_N # Intrinsic region length [cm]
# Voltage settings
V_ANODE = 0.8 # Anode voltage [V] - forward bias
V_CATHODE = 0.0 # Cathode voltage [V] - grounded
# Material parameters
q = 1.60217662e-19 # Elementary charge [C]
eps0 = 8.854e-14 # Vacuum permittivity [F/cm]
epsr_Si = 11.68 # Relative permittivity of silicon
eps_Si = eps0 * epsr_Si # Silicon permittivity [F/cm]
kB = 1.38064852e-23 # Boltzmann constant [J/K]
T = 300 # Temperature [K]
Vth = kB * T / q # Thermal voltage [V]
# Doping concentrations (based on Figure 14 in the paper)
ni = 1.08738184e10 # Intrinsic carrier concentration [cm^-3]
Na_p = 8e17 # Acceptor doping concentration in P-region [cm^-3]
Nd_n = 2e19 # Donor doping concentration in N-region [cm^-3]
N_i = 8e13 # Background doping concentration in intrinsic region [cm^-3]
# Carrier mobility and diffusion coefficients
mun = 1417.0 # Electron mobility [cm^2/V/s]
mup = 470.5 # Hole mobility [cm^2/V/s]
Dn = Vth * mun # Electron diffusion coefficient [cm^2/s]
Dp = Vth * mup # Hole diffusion coefficient [cm^2/s]
# Recombination parameters
tau_n = 1.0e-3 # Electron lifetime [s]
tau_p = 3.0e-4 # Hole lifetime [s]
##############################
### 2. Scaling and Dimensionless Parameters ###
##############################
# Scaling parameters
L_scale = H_DIODE # Length scaling factor [cm]
V_scale = Vth # Potential scaling factor [V]
C_scale = max(Na_p, Nd_n) # Concentration scaling factor [cm^-3]
mu_scale = max(mun, mup) # Mobility scaling factor [cm^2/V/s]
# Other scaling factors
t_scale = L_scale ** 2 / (V_scale * mu_scale) # Time scaling factor [s]
J_scale = q * V_scale * mu_scale * C_scale / L_scale # Current density scaling factor [A/cm^2]
R_scale = C_scale / t_scale # Recombination rate scaling factor
# Dimensionless parameters
# Poisson equation coefficient (relates electric field to charge density)
lmda = eps_Si * V_scale / (L_scale ** 2 * q * C_scale)
mun_norm = mun / mu_scale # Dimensionless electron mobility
mup_norm = mup / mu_scale # Dimensionless hole mobility
Dn_norm = Dn / (mu_scale * V_scale) # Dimensionless electron diffusion coefficient
Dp_norm = Dp / (mu_scale * V_scale) # Dimensionless hole diffusion coefficient
tau_n_norm = tau_n / t_scale # Dimensionless electron lifetime
tau_p_norm = tau_p / t_scale # Dimensionless hole lifetime
ni_norm = ni / C_scale # Dimensionless intrinsic carrier concentration
Na_p_norm = Na_p / C_scale # Dimensionless P-region doping concentration
Nd_n_norm = Nd_n / C_scale # Dimensionless N-region doping concentration
N_i_norm = N_i / C_scale # Dimensionless intrinsic region doping concentration
# Normalized voltage values
V_anode_norm = V_ANODE / V_scale
V_cathode_norm = V_CATHODE / V_scale
# Normalized geometric parameters
L_diode_norm = L_DIODE / L_scale
H_diode_norm = H_DIODE / L_scale
L_p_norm = L_P / L_scale
L_n_norm = L_N / L_scale
L_i_norm = L_I / L_scale
##############################
### 3. Mesh Generation ###
##############################
# Mesh resolution
nx = 120 # Number of elements in x direction (ensure ~10 elements per micron)
ny = 20 # Number of elements in y direction (ensure ~20 elements per micron)
# Create rectangular mesh
comm = MPI.COMM_WORLD
msh = create_rectangle(
comm,
points=((0.0, 0.0), (L_diode_norm, H_diode_norm)), # Normalized to [0,1]x[0,H/L]
n=(nx, ny),
cell_type=CellType.triangle,
)
# Get mesh coordinates for further processing
x = SpatialCoordinate(msh)
# Define subdomains
def p_region(x):
"""P-region: left side with length L_p"""
return np.logical_and(x[0] <= L_p_norm, x[1] >= 0)
def i_region(x):
"""Intrinsic region: middle with length L_i"""
return np.logical_and.reduce((
x[0] > L_p_norm,
x[0] < (L_p_norm + L_i_norm),
x[1] >= 0
))
def n_region(x):
"""N-region: right side with length L_n"""
return np.logical_and(x[0] >= (L_p_norm + L_i_norm), x[1] >= 0)
# Mark subdomains
cells = msh.topology.index_map(msh.topology.dim).size_local
cell_markers = np.zeros(cells, dtype=np.int32)
p_cells = dolfinx.mesh.locate_entities(msh, msh.topology.dim, p_region)
i_cells = dolfinx.mesh.locate_entities(msh, msh.topology.dim, i_region)
n_cells = dolfinx.mesh.locate_entities(msh, msh.topology.dim, n_region)
cell_markers[p_cells] = 1 # P-region
cell_markers[i_cells] = 2 # I-region
cell_markers[n_cells] = 3 # N-region
cell_tags = dolfinx.mesh.meshtags(msh, msh.topology.dim, np.arange(cells), cell_markers)
# Create measures based on subdomains for integration
dx = ufl.Measure('dx', domain=msh, subdomain_data=cell_tags)
# Mark boundaries
def anode_boundary(x):
"""Anode boundary (left side)"""
return np.isclose(x[0], 0.0)
def cathode_boundary(x):
"""Cathode boundary (right side)"""
return np.isclose(x[0], L_diode_norm)
def top_boundary(x):
"""Top boundary"""
return np.isclose(x[1], H_diode_norm)
def bottom_boundary(x):
"""Bottom boundary"""
return np.isclose(x[1], 0.0)
# Mark boundaries
msh.topology.create_entities(msh.topology.dim - 1)
facets = msh.topology.index_map(msh.topology.dim - 1).size_local
facet_markers = np.zeros(facets, dtype=np.int32)
anode_facets = dolfinx.mesh.locate_entities_boundary(msh, msh.topology.dim - 1, anode_boundary)
cathode_facets = dolfinx.mesh.locate_entities_boundary(msh, msh.topology.dim - 1, cathode_boundary)
top_facets = dolfinx.mesh.locate_entities_boundary(msh, msh.topology.dim - 1, top_boundary)
bottom_facets = dolfinx.mesh.locate_entities_boundary(msh, msh.topology.dim - 1, bottom_boundary)
# Boundary markers: 1-anode, 2-cathode, 3-top, 4-bottom
facet_markers[anode_facets] = 1
facet_markers[cathode_facets] = 2
facet_markers[top_facets] = 3
facet_markers[bottom_facets] = 4
facet_tags = dolfinx.mesh.meshtags(msh, msh.topology.dim - 1, np.arange(facets), facet_markers)
# Create measures based on boundaries for integration
ds = ufl.Measure('ds', domain=msh, subdomain_data=facet_tags)
dS = ufl.Measure('dS', domain=msh, subdomain_data=facet_tags)
##############################
### 4. Doping Distribution and Initialization ###
##############################
# Create Discontinuous Galerkin (DG) function space
elem_dg1 = basix.ufl.element("DG", msh.ufl_cell().cellname(), 2)
elem_bdm = basix.ufl.element("BDM", msh.ufl_cell().cellname(), 1)
V_dg = dolfinx.fem.functionspace(msh, elem_dg1)
# Doping distribution function
def doping_function(x):
"""Define doping distribution Nd-Na"""
n_doping = np.zeros_like(x[0])
n_mask = n_region(x)
n_doping[n_mask] = Nd_n_norm
i_mask = i_region(x)
n_doping[i_mask] = N_i_norm
p_doping = np.zeros_like(x[0])
p_mask = p_region(x)
p_doping[p_mask] = Na_p_norm
return n_doping - p_doping # Net doping (Nd-Na)
# Set up doping function
doping = Function(V_dg, name="doping")
doping.interpolate(doping_function)
# Define mixed element for the HDG method
# Using two spaces: DG space for local computations, trace space for boundaries and interfaces
elem_mix = basix.ufl.mixed_element([elem_dg1, elem_dg1, elem_dg1, elem_bdm, elem_bdm, elem_bdm])
V = dolfinx.fem.functionspace(msh, elem_mix)
# Create test functions for the variational form
(v_psi, v_n, v_p, v_e, v_jn, v_jp) = TestFunctions(V)
# Solution function
u_all = Function(V)
# Split solution function into components
(u_psi, u_n, u_p, u_e, u_jn, u_jp) = ufl.split(u_all)
##############################
### 5. Define Carrier Recombination Rate Functions ###
##############################
# Shockley-Read-Hall recombination rate
def R_SRH(n_val, p_val):
"""Shockley-Read-Hall recombination rate"""
return (n_val * p_val - ni_norm ** 2) / (tau_p_norm * (n_val + ni_norm) + tau_n_norm * (p_val + ni_norm) + EPS)
##############################
### 6. Compute Boundary Conditions ###
##############################
# Calculate equilibrium carrier concentrations at ohmic contacts in heavily doped regions
def calc_equilibrium_carriers(N_doping):
"""Calculate equilibrium carrier concentrations for a given doping level"""
n_eq = 0.5 * (N_doping + np.sqrt(N_doping ** 2 + 4 * ni_norm ** 2))
p_eq = 0.5 * (-N_doping + np.sqrt(N_doping ** 2 + 4 * ni_norm ** 2))
return n_eq, p_eq
# Calculate anode boundary conditions (P-region)
n_anode, p_anode = calc_equilibrium_carriers(-Na_p_norm)
psi_anode = V_anode_norm + np.log(n_anode / ni_norm)
# Calculate cathode boundary conditions (N-region)
n_cathode, p_cathode = calc_equilibrium_carriers(Nd_n_norm)
psi_cathode = V_cathode_norm + np.log(n_cathode / ni_norm)
# Normalize potential to avoid numerical issues
min_V = np.min([psi_anode, psi_cathode])
psi_anode = psi_anode - min_V
psi_cathode = psi_cathode - min_V
print(f"Anode boundary conditions: n={n_anode}, p={p_anode}, psi={psi_anode}")
print(f"Cathode boundary conditions: n={n_cathode}, p={p_cathode}, psi={psi_cathode}")
##############################
### 7. Initialize Variables ###
##############################
def init_carrier_concentrations():
"""Initialize carrier concentrations and potential (using linear interpolation)"""
# Create first-order continuous Galerkin element function space
elem_cg1 = basix.ufl.element("CG", msh.ufl_cell().cellname(), 1)
V_cg = dolfinx.fem.functionspace(msh, elem_cg1)
# Create function space for doping function
doping_cg = dolfinx.fem.functionspace(msh, elem_cg1)
doping = Function(doping_cg, name="doping")
doping.interpolate(doping_function)
# Create coordinate function for linear interpolation
x = SpatialCoordinate(msh)
# Use linear interpolation directly, from anode to cathode
# Assumes x[0] is the coordinate along the device length
L = L_diode_norm
# Create function for initialization
psi_init = Function(V_cg)
psi_init.interpolate(lambda x: psi_anode + (psi_cathode - psi_anode) * x[0] / L)
# Interpolate CG solution to DG space
u_all.sub(0).interpolate(psi_init)
# Create carrier concentration functions
n_init = Function(V_cg)
p_init = Function(V_cg)
doping_values = doping.x.array[:]
# Calculate equilibrium carrier concentrations
n_eq_values = np.zeros_like(doping_values)
p_eq_values = np.zeros_like(doping_values)
for i in range(len(doping_values)):
n_eq_values[i] = 0.5 * (doping_values[i] + np.sqrt(doping_values[i] ** 2 + 4 * ni_norm ** 2))
p_eq_values[i] = 0.5 * (-doping_values[i] + np.sqrt(doping_values[i] ** 2 + 4 * ni_norm ** 2))
# Set values to functions
n_init.x.array[:] = n_eq_values
p_init.x.array[:] = p_eq_values
# Set to final solution functions
u_all.sub(1).interpolate(n_init)
u_all.sub(2).interpolate(p_init)
print("Variables initialized successfully")
# Execute initialization
init_carrier_concentrations()
# Uncomment to plot data: plot_mosfet_data(n, p, psi, doping)
##############################
### 8. Full-Coupled HDG Method Implementation ###
##############################
# Define normal vector and cell size for numerical stability
n = ufl.FacetNormal(msh)
h = ufl.CellDiameter(msh)
h_avg = ufl.avg(h)
# The semiconductor equations system we're solving:
# λ²∇·E + n - p - N = 0 (Poisson equation)
# -∇·Jn + R(n,p) = 0 (Electron continuity)
# ∇·Jp + R(n,p) = 0 (Hole continuity)
# E + ∇ψ = 0 (Electric field definition)
# Jn - μn·n·E - Dn∇n = 0 (Electron current density)
# Jp - μp·p·E + Dp∇p = 0 (Hole current density)
# Define Nitsche method function for Dirichlet boundary conditions
def nitsche(u, v, g, coeff, n_facet, h, alpha, ds_k):
"""Nitsche method for weakly enforcing Dirichlet boundary conditions"""
return (
-inner(coeff * grad(u), n_facet * v) * ds_k
- inner(coeff * grad(v), n_facet * u) * ds_k
+ alpha / h * inner(u, v) * ds_k
+ inner(coeff * grad(v), n_facet * g) * ds_k
- alpha / h * inner(g, v) * ds_k
)
# Define SIPG method for interior interface handling
def sipg(u, v, coeff, n_facet, h_avg, alpha, dS):
"""Symmetric Interior Penalty Galerkin method"""
return (
-inner(ufl.avg(coeff * grad(u)), jump(v, n_facet)) * dS
- inner(ufl.avg(coeff * grad(v)), jump(u, n_facet)) * dS
+ alpha / h_avg * inner(jump(u, n_facet), jump(v, n_facet)) * dS
)
# Standard interior penalty method
def standard_ip(u, v, n, h_avg, gamma):
"""Standard interior penalty method"""
F = -inner(jump(v, n), ufl.avg(grad(u))) * dS
F += -inner(ufl.avg(grad(v)), jump(u, n)) * dS
F += +gamma / h_avg * inner(jump(v, n), jump(u, n)) * dS
return F
# Set penalty parameters
alpha_psi = Constant(msh, ScalarType(20.0)) # Penalty parameter for potential
alpha_n = Constant(msh, ScalarType(20.0)) # Penalty parameter for electron concentration
alpha_p = Constant(msh, ScalarType(20.0)) # Penalty parameter for hole concentration
# Constant coefficients
lam2 = Constant(msh, ScalarType(lmda)) # Poisson equation coefficient
Dn_c = Constant(msh, ScalarType(Dn_norm)) # Normalized electron diffusion coefficient
Dp_c = Constant(msh, ScalarType(Dp_norm)) # Normalized hole diffusion coefficient
mun_c = Constant(msh, ScalarType(mun_norm)) # Normalized electron mobility
mup_c = Constant(msh, ScalarType(mup_norm)) # Normalized hole mobility
# Build variational form for the six-variable equation system
# 1. Poisson equation: λ²∇·E + n - p - N = 0
# Integration by parts: ∫λ²∇·E v_psi dx = -∫λ²E·∇v_psi dx + ∫λ²(E·n)v_psi ds
F1 = -lam2 * inner(u_e, grad(v_psi)) * dx + inner((u_n - u_p - doping), v_psi) * dx
# Handle interior interface jumps
F1 += standard_ip(u_psi, v_psi, n, h_avg, alpha_psi)
# Handle boundary conditions
F1 += nitsche(u_psi, v_psi, psi_anode, lam2, n, h, alpha_psi, ds(1))
F1 += nitsche(u_psi, v_psi, psi_cathode, lam2, n, h, alpha_psi, ds(2))
# 2. Electron continuity equation: -∇·Jn + R = 0
# Integration by parts: ∫-∇·Jn v_n dx = ∫Jn·∇v_n dx - ∫(Jn·n)v_n ds
F2 = inner(u_jn, grad(v_n)) * dx + inner(R_SRH(u_n, u_p), v_n) * dx
# Handle interior interface jumps
F2 += standard_ip(u_n, v_n, n, h_avg, alpha_n)
# Handle boundary conditions
F2 += nitsche(u_n, v_n, n_anode, Dn_c, n, h, alpha_n, ds(1)) # Anode
F2 += nitsche(u_n, v_n, n_cathode, Dn_c, n, h, alpha_n, ds(2)) # Cathode
# 3. Hole continuity equation: ∇·Jp + R = 0
# Integration by parts: ∫∇·Jp v_p dx = -∫Jp·∇v_p dx + ∫(Jp·n)v_p ds
F3 = -inner(u_jp, grad(v_p)) * dx + inner(R_SRH(u_n, u_p), v_p) * dx
# Handle interior interface jumps
F3 += standard_ip(u_p, v_p, n, h_avg, alpha_p)
# Handle boundary conditions
F3 += nitsche(u_p, v_p, p_anode, Dp_c, n, h, alpha_p, ds(1))
F3 += nitsche(u_p, v_p, p_cathode, Dp_c, n, h, alpha_p, ds(2))
# 4. Electric field equation (algebraic): E + ∇ψ = 0
F4 = inner(u_e + grad(u_psi), v_e) * dx
# 5. Electron current density equation (algebraic): Jn - μn·n·E - Dn·∇n = 0
F5 = inner(u_jn - mun_c * u_n * u_e - Dn_c * grad(u_n), v_jn) * dx
# 6. Hole current density equation (algebraic): Jp - μp·p·E + Dp·∇p = 0
F6 = inner(u_jp - mup_c * u_p * u_e + Dp_c * grad(u_p), v_jp) * dx
# Neumann boundary conditions are naturally satisfied
# Assemble complete variational form
F = F1 + F2 + F3 + F4 + F5 + F6
# Calculate Jacobian matrix
J = ufl.derivative(F, u_all)
# Define problem and solver
problem = NonlinearProblem(F, u_all, J=J)
solver = dolfinx.nls.petsc.NewtonSolver(msh.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.atol = 1e-12
solver.report = True
# Solve the problem
solver.solve(u_all)
Currently, I suspect that there are errors in my variational derivations for the boundary and internal continuity, especially in the coefficient terms of the Nitsche’s method. I’d be extremely grateful for any relevant tips.