Steady state problem

Steady state problem#

Let’s consider a unit square \(\Omega = [0, 1] \times [0, 1]\).

\(c_m\) is the concentration of mobile hydrogen, and \(c_t\) is the concentration of trapped hydrogen.

We want to solve the following steady state problem:

(21)#\[\begin{align} \nabla \cdot (D \nabla c_m) - R &= 0\quad \text{on } \Omega\\ R &= 0\quad \text{on } \Omega\\ R &= k c_m (n - c_t) - p \ c_t \\ c_m &= 1 \quad \text{on } \Gamma_\mathrm{inlet} \\ c_m &= 0\quad \text{on } \Gamma_\mathrm{outlet} \end{align}\]

Note

In steady state, this problem is equivalent to:

(22)#\[\begin{align} \nabla \cdot (D \nabla c_m) &= 0\quad \text{on } \Omega\\ R &= 0\quad \text{on } \Omega \end{align}\]

Meaning that for the governing equation for mobile concentration is the same as a pure diffusive (no trap) case.

Furthermore,

\[ k c_m (n - c_t) = p \ c_t \]

Leading to the direct expression for \(c_t\):

\[ c_t = n \left(1 + \frac{p}{k\ c_m} \right)^{-1} \]

This is sometimes referred as Oriani’s equilibrium. We won’t solve it the “direct” way here as we are building up to a transient (non-equilibrium) case.

Setting up the problem#

We start by importing the required modules and libraries:

import dolfinx
import numpy as np
import ufl
from dolfinx.fem.petsc import NonlinearProblem
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import basix

We then create a mesh using create_unit_square

nx = ny = 96

domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.quadrilateral)

Now that we have a mesh, we need to define a function space and appropriate functions and test functions:

# we create a mixed element with two components, both continuous galerkin degree 1
cg_element = basix.ufl.element("Lagrange", domain.basix_cell(), degree=1)

mixed_element = basix.ufl.mixed_element([cg_element, cg_element])

# then we make a functionspace from the mixed element
V = fem.functionspace(domain, mixed_element)

# we create a "main" function u which is a vector of the two components
u = fem.Function(V)

# to use the components in variational forms, we use ufl.split
# the first will be the mobile concentration cm, the second the trapped concentration ct
cm, ct = ufl.split(u)

# we create test functions for both components
v_cm, v_ct = ufl.TestFunctions(V)

Every finite element problem needs boundary conditions. Here we define three Dirichlet boundary conditions:

def inlet(x):
    return np.logical_and(np.isclose(x[0], 0), x[1] <= 0.5)

def outlet(x):
    return np.logical_and(np.isclose(x[0], 1), x[1] >= 0.5)

V0, submap = V.sub(0).collapse()

# the trick here was to pass both the subspace and the collapsed space to locate_dofs_geometrical
# in FESTIM we don't need this since we use meshtags for everything
# https://fenicsproject.discourse.group/t/dolfinx-dirichlet-bcs-for-mixed-function-spaces/7844/2

dofs_outlet = fem.locate_dofs_geometrical((V.sub(0), V0), outlet)
dofs_inlet = fem.locate_dofs_geometrical((V.sub(0), V0), inlet)

c_inlet = fem.Constant(domain, 1.0)
c_outlet = fem.Constant(domain, 0.0)

bc_outlet = fem.dirichletbc(c_outlet, dofs_outlet[0], V.sub(0))
bc_inlet = fem.dirichletbc(c_inlet, dofs_inlet[0], V.sub(0))

We then define the variational formulation (weak form)

# Problem parameters
k = 0.1  # trapping rate
p = 0.1  # detrapping rate
n = 1  # total trapping sites
D = 2.0 # diffusion coefficient

trapping = k * cm * (n - ct)
detrapping = p * ct

# NOTE everything is bundled in one variational form F
# the difference between the different equations is made with the test functions v_cm and v_ct
F_mobile = (
    D*ufl.dot(ufl.grad(cm), ufl.grad(v_cm)) * ufl.dx
    - trapping * v_cm * ufl.dx
    + detrapping * v_cm * ufl.dx
)
F_trapped = +trapping * v_ct * ufl.dx - detrapping * v_ct * ufl.dx

F = F_mobile + F_trapped

Now we can create a Newton solver:

# taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28
use_superlu = PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"

petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 1e-10,
    "snes_rtol": 1e-10,
    "snes_max_it": 100,
    "snes_divergence_tolerance": "PETSC_UNLIMITED",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
}

problem = NonlinearProblem(
    F,
    u,
    bcs=[bc_outlet, bc_inlet],
    petsc_options=petsc_options,
    petsc_options_prefix="Poisson",
)

Solving#

problem.solve()

converged = problem.solver.getConvergedReason()
num_iter = problem.solver.getIterationNumber()

assert converged > 0, f"Solver did not converge, got {converged}."
print(
    f"Solver converged after {num_iter} iterations with converged reason {converged}."
)

Solver converged after 2 iterations with converged reason 2.

Post processing#

# we first split the main solution u into its components with .split()
cm_post, ct_post = u.split()  # NOTE this is different from ufl.split(u)

# for postprocessing, it's easier to work with collapsed functions
cm_post = cm_post.collapse()
ct_post = ct_post.collapse()

Visualise the results with `pyvista`#

import pyvista
from dolfinx import plot

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(tdim, tdim)

# We can do this once since both components share the same function space
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(cm_post.function_space)


plotter = pyvista.Plotter(shape=(1, 2))


plotter.subplot(0, 0)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["cm"] = cm_post.x.array.real
u_grid.set_active_scalars("cm")
plotter.add_mesh(u_grid, show_edges=False)
plotter.view_xy()

plotter.subplot(0, 1)
ct_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
ct_grid.point_data["ct"] = ct_post.x.array.real
ct_grid.set_active_scalars("ct")
plotter.add_mesh(ct_grid, show_edges=False)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()

2026-05-22 15:22:33.314 (   0.514s) [    761CEDACA740]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

Export the results to a file#

from dolfinx.io import VTXWriter

writer = VTXWriter(MPI.COMM_WORLD, "solution.bp", cm_post)

writer.write(t=0.0)

Compute derived quantities#

Inventories:

from scifem import assemble_scalar
inventory_cm = assemble_scalar(cm_post * ufl.dx)
inventory_ct = assemble_scalar(ct_post * ufl.dx)

print(f"Total inventory of mobile species: {inventory_cm:.4f}")
print(f"Total inventory of trapped species: {inventory_ct:.4f}")

Total inventory of mobile species: 0.5000
Total inventory of trapped species: 0.3179

Outgassing fluxes:

inlet_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, inlet)
outlet_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, outlet)

inlet_values = np.full_like(inlet_facets, 1.0)
outlet_values = np.full_like(outlet_facets, 2.0)

entities = np.concatenate([inlet_facets, outlet_facets])
values = np.concatenate([inlet_values, outlet_values])

facet_meshtags = dolfinx.mesh.meshtags(domain, fdim, entities=entities, values=values)

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_meshtags)

n = ufl.FacetNormal(domain)

flux_form = -ufl.dot(D * ufl.grad(cm_post), n)

flux_inlet = -assemble_scalar(flux_form * ds(1))
flux_outlet = assemble_scalar(flux_form * ds(2))

print(f"Inlet flux: {flux_inlet:.4f}")
print(f"Outlet flux: {flux_outlet:.4f}")

print(f"Rel difference: {(flux_inlet - flux_outlet)/flux_inlet:.2%}")

Inlet flux: 1.2679
Outlet flux: 1.2679
Rel difference: -0.00%