Radioactive decay 2D

Radioactive decay 2D#

🏷 Tags: 2D MMS RadioactiveDecay transient

This MMS case verifies the implementation of radioactive decay in FESTIM. We will only consider diffusion of hydrogen in a unit square domain \(\Omega\) at steady state with a homogeneous diffusion coefficient \(D\). We will impose a radioactive decay (RadioactiveDecay) source over the whole domain with a decay constant \(\lambda\). Moreover, a Dirichlet boundary condition will be assumed on the boundaries \(\partial \Omega \).

The problem is therefore:

(22)#\[\begin{split} \begin{align} &\nabla \cdot (D \ \nabla{c}) -\lambda c = -S \quad \text{on } \Omega \\ & c = c_0 \quad \text{on } \partial \Omega \end{align} \end{split}\]

The manufactured exact solution for mobile concentration is:

(23)#\[ \begin{equation} c_\mathrm{exact} = 1 + 2 x^2 + 3 y^2 \end{equation} \]

Injecting (23) in (22), we obtain the expressions of \(S\) and \(c_0\):

\[\begin{split} \begin{align} & S = \lambda \left(1 + 2 x^2 + 3 y^2 \right) -10 D \\ & c_0 = c_\mathrm{exact} \end{align} \end{split}\]

We can then run a FESTIM model with these values and compare the numerical solution with \(c_\mathrm{exact}\).

import festim as F
import numpy as np

import dolfinx
from mpi4py import MPI
import ufl

decay_constant = 3
D = 2

# Create the FESTIM model
my_model = F.HydrogenTransportProblem()

# Create and mark the mesh
nx = ny = 20
fenics_mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, nx, ny, cell_type=dolfinx.mesh.CellType.quadrilateral
)
my_model.mesh = F.Mesh(fenics_mesh)

my_mat = F.Material(D_0=D, E_D=0)
volume = F.VolumeSubdomain(id=1, material=my_mat)
boundary = F.SurfaceSubdomain(id=1)

my_model.subdomains = [volume, boundary]

H = F.Species("H")
my_model.species = [H]

# Variational formulation
exact_solution = lambda x: 1 + 2 * x[0] ** 2 + 3 * x[1] ** 2  # exact solution

def S(x):
    return decay_constant * exact_solution(x) - 10 * D

my_model.sources = [
    F.ParticleSource(
        value=S, volume=volume, species=H
    ),
]

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=boundary, value=exact_solution, species=H),
]

decay_reaction = F.Reaction(reactant=H, k_0=decay_constant, E_k=0, volume=volume)

my_model.reactions = [decay_reaction]


my_model.temperature = 500.0  # ignored in this problem

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)

my_model.initialise()
my_model.run()

/home/docs/checkouts/readthedocs.org/user_builds/festim-vv-report/conda/festim-2/lib/python3.11/site-packages/festim/coupled_heat_hydrogen_problem.py:1: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  import tqdm.autonotebook

Comparison with exact solution#

First, we compute the \(L^2\)-norm of the error, defined by \(E=\sqrt{\int_\Omega (c-c_\mathrm{exact})^2\mathrm{d} x}\). Secondly, we compute the maximum error at any degree of freedom.

def error_L2(u_computed, u_exact, degree_raise=3):
    # Create higher order function space
    degree = u_computed.function_space.ufl_element().degree
    family = u_computed.function_space.ufl_element().family_name
    mesh = u_computed.function_space.mesh
    W = dolfinx.fem.functionspace(mesh, (family, degree + degree_raise))

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = dolfinx.fem.Function(W)
    if isinstance(u_exact, ufl.core.expr.Expr):
        u_expr = dolfinx.fem.Expression(u_exact, W.element.interpolation_points)
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_exact)

    # Integrate the error
    error = dolfinx.fem.form(
        ufl.inner(u_computed - u_ex_W, u_computed - u_ex_W) * ufl.dx
    )
    error_local = dolfinx.fem.assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)

Hide code cell source

 
computed_solution = H.solution

E_l2 = error_L2(computed_solution, exact_solution)

exact_solution_function = dolfinx.fem.Function(computed_solution.function_space)
exact_solution_function.interpolate(exact_solution)
E_max = np.max(np.abs(exact_solution_function.x.array-computed_solution.x.array))

print(f"L2 error: {E_l2:.2e}")
print(f"Max error: {E_max:.2e}")

L2 error: 2.09e-03
Max error: 2.13e-04

import pyvista
from dolfinx.plot import vtk_mesh

pyvista.start_xvfb()
pyvista.set_jupyter_backend('html')

u_topology, u_cell_types, u_geometry = vtk_mesh(computed_solution.function_space)

u_grid_computed = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid_computed.point_data["c"] = computed_solution.x.array.real
u_grid_computed.set_active_scalars("c")

u_grid_exact = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid_exact.point_data["c_exact"] = exact_solution_function.x.array.real
u_grid_exact.set_active_scalars("c_exact")


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

u_plotter.subplot(0, 0)
u_plotter.add_mesh(u_grid_computed, show_edges=True)
u_plotter.view_xy()

u_plotter.subplot(0, 1)
u_plotter.add_mesh(u_grid_exact, show_edges=True)
u_plotter.view_xy()


if not pyvista.OFF_SCREEN:
    u_plotter.show()
else:
    figure = u_plotter.screenshot("decay.png")

/home/docs/checkouts/readthedocs.org/user_builds/festim-vv-report/conda/festim-2/lib/python3.11/site-packages/pyvista/plotting/utilities/xvfb.py:48: PyVistaDeprecationWarning: This function is deprecated and will be removed in future version of PyVista. Use vtk-osmesa instead.
  warnings.warn(
Contents