Coupled solution of variational form semiconductor equations

Hi everyone,
I am new to Fenics/dolfinx and I am struggling to solve a set of coupled equations. Here is a mathematical description of the problem which is the steady state Poisson-Nernst-Planck set of equations with two charged species, with variables to be solved \phi, c_1 and c_2.

(1) \quad -\nabla ^2 \phi = \frac{(c_1 - c_2)}{\varepsilon}
(2) \quad \nabla \cdot \left[ D_1 \nabla c_1 + \frac{D_1z_1q}{kT}c_1\nabla\phi\right] = 0
(3) \quad \nabla \cdot \left[ D_2 \nabla c_2 + \frac{D_2z_2q}{kT}c_2\nabla\phi\right] = 0

Then scaling the variables and simplifying I get,

(1) \quad -\nabla ^2 \bar{\phi} = (\bar{c}_1-\bar{c}_2)
(2) \quad \nabla^2\bar{c}_1 + \nabla\bar{c}_1\cdot \nabla\bar{\phi} +\bar{c}_1\nabla^2\bar{\phi} = 0
(2) \quad \nabla^2\bar{c}_2 + \nabla\bar{c}_2\cdot \nabla\bar{\phi} +\bar{c}_2\nabla^2\bar{\phi} = 0
where the bar represents the adimensional variable.

My first guess at the variational problem formulation is as follows
(1) \quad \int_\Omega \nabla\bar\phi \cdot \nabla v \, dx - \int_{\partial\Omega}\frac{\partial\bar\phi}{\partial n} v \, ds - \int_\Omega (c_1 - c_2) = 0
(2) \quad -\int_\Omega \nabla\bar{c}_1\cdot d_1 \, dx + \int_{\partial\Omega}\frac{\partial\bar{c}_1}{\partial n}d_1 \, ds + \int_\Omega \nabla \bar{c}_1\cdot \nabla\bar\phi \cdot d_1\, dx+ \int_\Omega \bar{c}_1\nabla^2\bar\phi \cdot d_1\, dx = 0
(3) \quad ....

However I am not sure if this is the correct formulation since when I implement it with a mixed element space and a newton solver it reaches the maximum number of iterations without converging.
The problem could still be elsewhere (I upload also the code below) but I thought I started asking for the most basic problem first.

#%% Import statements

from mpi4py import MPI
from petsc4py import PETSc
from import gmshio, XDMFFile
from dolfinx import mesh, fem
from dolfinx.fem import functionspace
import numpy as np
import ufl
from dolfinx import default_scalar_type, default_real_type
from dolfinx import geometry as gm
from dolfinx.fem.petsc import LinearProblem
import pyvista
from dolfinx import plot
    import gmsh
except ImportError:
    print("This demo requires gmsh to be installed")

import matplotlib.pyplot as plt

import dolfinx.fem.petsc
import basix.ufl
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

#%% Generate mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

#%% Define model
# Create connectivity
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)

element_family = basix.ElementFamily.P
element_degree = 1
variant = basix.LagrangeVariant.equispaced

# Define the element using Basix
el1 = basix.ufl.element(element_family, domain.topology.cell_name(), element_degree, variant)
el2 = basix.ufl.element(element_family, domain.topology.cell_name(), element_degree, variant)
el3 = basix.ufl.element(element_family, domain.topology.cell_name(), element_degree, variant)
mel = basix.ufl.mixed_element(([el1, el2, el3]))
Z = fem.functionspace(domain, mel)
V = fem.functionspace(domain, ("Lagrange", 1))
C1 = fem.functionspace(domain, ("Lagrange", 1))
C2 = fem.functionspace(domain, ("Lagrange", 1))

# Define test functions
(v, d1, d2) = ufl.TestFunction(Z)

u = fem.Function(Z) 

# Split mixed functions
phi, c1, c2 = ufl.split(u)
phi0, c10, c20 = ufl.split(u0)

# Zero u
u.x.array[:] = 0.0

# Variational formulation
F1 = ufl.inner(ufl.grad(phi), ufl.grad(v)) * ufl.dx - ufl.inner(c1-c2, v) * ufl.dx
F2 = ufl.inner(ufl.grad(c1), ufl.grad(d1)) * ufl.dx - ufl.inner(ufl.inner(ufl.grad(c1), ufl.grad(phi)), d1) * ufl.dx - ufl.inner(c1*ufl.div(ufl.grad(phi)), d1)*ufl.dx
F3 = ufl.inner(ufl.grad(c2), ufl.grad(d2)) * ufl.dx - ufl.inner(ufl.inner(ufl.grad(c2), ufl.grad(phi)), d2) * ufl.dx - ufl.inner(c2*ufl.div(ufl.grad(phi)), d2)*ufl.dx
F = F1 + F2 + F3

def bulk(x):
    return np.isclose(x[0], 0)

def surface(x):
    return np.isclose(x[0], 1)

def sides(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))

def zerofunction(x):
    return 0*x[0]

def onefunction(x):
    return 1+0*x[0]

boundaries = [(1, bulk),
              (2, surface),
              (3, sides)]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(domain, fdim, locator)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values, space, testfunction):
        self._type = type
        if type == "Dirichlet":
            u_D = fem.Function(V)
            facets = facet_tag.find(marker)
            dofs = fem.locate_dofs_topological(space, fdim, facets)
            self._bc = fem.dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = ufl.inner(values, testfunction) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * ufl.inner(u-values[1], testfunction)* ds(marker)
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    def bc(self):
        return self._bc

    def type(self):
        return self._type

# Define the Dirichlet conditions
zero_bulk_potential = zerofunction
applied_potential = onefunction
bulk_concentration = onefunction

# Define the Neumann conditions
noflux1 = fem.Function(C1)
noflux2 = fem.Function(C2)
no_charge = fem.Function(V)

# Define boundary conditions
boundary_conditions = [BoundaryCondition("Dirichlet", 1, zero_bulk_potential, V, v),
                       BoundaryCondition("Dirichlet", 2, applied_potential, V, v),
                       BoundaryCondition("Dirichlet", 1, bulk_concentration, C1, d1),
                       BoundaryCondition("Dirichlet", 1, bulk_concentration, C2, d2),
                       BoundaryCondition("Neumann", 3, no_charge, V, v),
                       BoundaryCondition("Neumann", 3, noflux1, C1, d1),
                       BoundaryCondition("Neumann", 3, noflux2, C2, d2),
                       BoundaryCondition("Neumann", 2, noflux1, C1, d1),
                       BoundaryCondition("Neumann", 2, noflux2, C2, d2)

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        F -= condition.bc

# %%Create nonlinear problem and Newton solver

problem = NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2
solver.max_it = 200

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
sys = PETSc.Sys()  # type: ignore
# For factorisation prefer superlu_dist, then MUMPS, then default
if sys.hasExternalPackage("superlu_dist"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
elif sys.hasExternalPackage("mumps"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"


Thank you in advance

I don’t have a direct answer but are you aware of simudo? It uses legacy fenics to solve the same set of equations. The associated publication has some details about the implementation choices they made with regards to the elements type and the nonlinear solver.

Thank you for your response and for the link, I did not know about it, it can be useful for more advanced solver settings.
However they solve for quasi-fermi levels and currents and not for concentrations, that is what I am interested in, so the problem is a little different in its base formulation.
In their algorithm I see a loop of Newton iterations, that I am not able to reproduce in fenics/dolfinx code. If my mathematical formulation is correct that may be the issue, but I don’t know how to do it