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 dolfinx.io 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
try:
import gmsh
except ImportError:
print("This demo requires gmsh to be installed")
exit(0)
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
# %% SET BOUNDARY CONDITIONS
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_indices.append(facets)
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)
u_D.interpolate(values)
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)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
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)
noflux1.interpolate(zerofunction)
noflux2 = fem.Function(C2)
noflux2.interpolate(zerofunction)
no_charge = fem.Function(V)
no_charge.interpolate(zerofunction)
# 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":
bcs.append(condition.bc)
else:
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"
ksp.setFromOptions()
solver.solve(u)
Thank you in advance