Simple transient diffusion case

Simple transient diffusion case#

This is a simple transient MMS example. We will only consider diffusion of hydrogen in a unit square domain \(\Omega\) at steady state with an homogeneous diffusion coefficient \(D\). Moreover, a Dirichlet boundary condition will be assumed on the boundaries \(\partial \Omega \).

The problem is therefore:

(5)#\[\begin{split} \begin{align} &\nabla \cdot (D \ \nabla{c}) - \frac{\partial c}{\partial t} = -S \quad \text{on } \Omega ; \ t\geq 0 \\ & c = c_0 \quad \text{on } \partial \Omega ; \ t\geq 0 \\ & c = c_\mathrm{initial} \quad \text{on } \partial \Omega ; \ \text{at } t=0 \end{align} \end{split}\]

The exact solution for mobile concentration is:

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

Note

We use a manufactured solution that varies linearly with time (\(t^1\)), as the backward Euler scheme provides an exact solution in this case.

Injecting (6) in (5), we obtain the expressions of \(S\), \(c_0\), and \(c_\mathrm{initial}\):

(7)#\[\begin{align} & S = 2 + 3 y^2 - (4 + 6t) D \\ & c_0 = c_\mathrm{exact} \\ & c_\mathrm{initial} = c_\mathrm{exact}(t=0) \end{align}\]

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

FESTIM code#

Show code cell content Hide code cell content

import festim as F
import sympy as sp
import fenics as f
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

# Create and mark the mesh
nx = ny = 100
fenics_mesh = f.UnitSquareMesh(nx, ny)


volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
volume_markers.set_all(1)

surface_markers = f.MeshFunction(
    "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1
)
surface_markers.set_all(0)


class Boundary(f.SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary


boundary = Boundary()
boundary.mark(surface_markers, 1)

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

my_model.mesh = F.Mesh(
    fenics_mesh, volume_markers=volume_markers, surface_markers=surface_markers
)

# Variational formulation
exact_solution = 1 + 2 * F.x**2 + 3 * F.t * F.y**2 + 2 * F.t  # exact solution

exact_solution_copy = exact_solution.subs(F.x, F.x)

D = 2

my_model.sources = [
    F.Source(2 + 3 * F.y**2 - (4 + 6 * F.t) * D, volume=1, field="solute"),
]

my_model.boundary_conditions = [
    F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
]

my_model.materials = F.Material(id=1, D_0=D, E_D=0)

my_model.T = F.Temperature(500)  # ignored in this problem

xdmf_file_name = "simple_transient_mobile.xdmf"
my_model.exports = [
    F.XDMFExport(field="solute", filename=xdmf_file_name, checkpoint=True)
]

final_time = 17
slices = 4
slice_size = final_time / slices
milestones = list(np.linspace(slice_size, final_time, slices))

my_model.dt = F.Stepsize(initial_value=1, milestones=milestones)

my_model.settings = F.Settings(
    absolute_tolerance=1e-10,
    relative_tolerance=1e-10,
    final_time=final_time,
)

my_model.initialise()
my_model.run()

Defining initial values
Defining variational problem
Defining source terms
Defining boundary conditions

Time stepping...
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

9 %        1.0e+00 s    Elapsed time so far: 1.2 s
8 %        2.0e+00 s    Elapsed time so far: 1.3 s
6 %        3.0e+00 s    Elapsed time so far: 1.4 s

5 %        4.0e+00 s    Elapsed time so far: 1.5 s
0 %        4.2e+00 s    Elapsed time so far: 1.6 s
5 %        4.5e+00 s    Elapsed time so far: 1.7 s

9 %        4.8e+00 s    Elapsed time so far: 1.8 s
4 %        5.0e+00 s    Elapsed time so far: 1.9 s
9 %        5.2e+00 s    Elapsed time so far: 2.0 s

4 %        5.5e+00 s    Elapsed time so far: 2.1 s
8 %        5.8e+00 s    Elapsed time so far: 2.2 s
3 %        6.0e+00 s    Elapsed time so far: 2.3 s

8 %        6.2e+00 s    Elapsed time so far: 2.4 s
2 %        6.5e+00 s    Elapsed time so far: 2.5 s
7 %        6.8e+00 s    Elapsed time so far: 2.6 s

2 %        7.0e+00 s    Elapsed time so far: 2.6 s
6 %        7.2e+00 s    Elapsed time so far: 2.7 s
1 %        7.5e+00 s    Elapsed time so far: 2.8 s

6 %        7.8e+00 s    Elapsed time so far: 2.9 s
1 %        8.0e+00 s    Elapsed time so far: 3.0 s
5 %        8.2e+00 s    Elapsed time so far: 3.1 s

0 %        8.5e+00 s    Elapsed time so far: 3.2 s
5 %        8.8e+00 s    Elapsed time so far: 3.3 s
9 %        9.0e+00 s    Elapsed time so far: 3.4 s

4 %        9.2e+00 s    Elapsed time so far: 3.5 s
9 %        9.5e+00 s    Elapsed time so far: 3.6 s
4 %        9.8e+00 s    Elapsed time so far: 3.7 s

8 %        1.0e+01 s    Elapsed time so far: 3.8 s
3 %        1.0e+01 s    Elapsed time so far: 3.9 s
8 %        1.0e+01 s    Elapsed time so far: 4.0 s

2 %        1.1e+01 s    Elapsed time so far: 4.1 s
7 %        1.1e+01 s    Elapsed time so far: 4.2 s
2 %        1.1e+01 s    Elapsed time so far: 4.2 s

7 %        1.2e+01 s    Elapsed time so far: 4.4 s
1 %        1.2e+01 s    Elapsed time so far: 4.5 s
6 %        1.2e+01 s    Elapsed time so far: 4.5 s

1 %        1.2e+01 s    Elapsed time so far: 4.6 s
5 %        1.2e+01 s    Elapsed time so far: 4.7 s
0 %        1.3e+01 s    Elapsed time so far: 4.8 s

5 %        1.3e+01 s    Elapsed time so far: 4.9 s
9 %        1.3e+01 s    Elapsed time so far: 5.0 s
4 %        1.4e+01 s    Elapsed time so far: 5.1 s

9 %        1.4e+01 s    Elapsed time so far: 5.2 s
3 %        1.4e+01 s    Elapsed time so far: 5.3 s
8 %        1.4e+01 s    Elapsed time so far: 5.4 s

3 %        1.4e+01 s    Elapsed time so far: 5.5 s
8 %        1.5e+01 s    Elapsed time so far: 5.6 s
2 %        1.5e+01 s    Elapsed time so far: 5.7 s

7 %        1.5e+01 s    Elapsed time so far: 5.8 s
2 %        1.6e+01 s    Elapsed time so far: 5.9 s
7 %        1.6e+01 s    Elapsed time so far: 6.0 s

1 %        1.6e+01 s    Elapsed time so far: 6.1 s
6 %        1.6e+01 s    Elapsed time so far: 6.2 s
1 %        1.6e+01 s    Elapsed time so far: 6.3 s

98.5 %        1.7e+01 s    Elapsed time so far: 6.4 s
100.0 %        1.7e+01 s    Elapsed time so far: 6.5 s

Comparison with exact solution#

Show code cell source Hide code cell source

from fenics import XDMFFile, FunctionSpace, Function, plot


def load_xdmf(mesh, filename, field, element="CG", counter=-1):
    """Loads a XDMF file and store its content to a fenics.Function

    Args:
        mesh (fenics.mesh): the mesh of the function
        filename (str): the XDMF filename
        field (str): the name of the field in the XDMF file
        element (str, optional): Finite element of the function.
            Defaults to "CG".
        counter (int, optional): timestep in the file, -1 is the
            last timestep. Defaults to -1.

    Returns:
        fenics.Function: the content of the XDMF file as a Function
    """

    V = FunctionSpace(mesh, element, 1)
    u = Function(V)

    XDMFFile(filename).read_checkpoint(u, field, counter)
    return u


fig, axs = plt.subplots(
    slices, 3, figsize=(slices * 2.3, slices * 2.5 + 1)
)  # tweak figsize if needed
fig.tight_layout()


def compute_arc_length(xs, ys):
    """Computes the arc length of x,y points based
    on x and y arrays
    """
    points = np.vstack((xs, ys)).T
    distance = np.linalg.norm(points[1:] - points[:-1], axis=1)
    arc_length = np.insert(np.cumsum(distance), 0, [0.0])
    return arc_length


def exists_close(x, list):
    return any(np.isclose(x, t) for t in list)


xdmf_times = F.extract_xdmf_times(xdmf_file_name)
counters = [i for (i, time) in enumerate(xdmf_times) if exists_close(time, milestones)]

for i, counter in enumerate(counters):
    time = xdmf_times[counter]

    c_exact = f.Expression(
        sp.printing.ccode(exact_solution_copy.subs(F.t, time)), degree=4
    )
    c_exact = f.project(c_exact, f.FunctionSpace(my_model.mesh.mesh, "CG", 1))

    computed_solution = load_xdmf(
        fenics_mesh, xdmf_file_name, "mobile_concentration", "CG", counter
    )
    E = f.errornorm(computed_solution, c_exact, "L2")

    # plot exact solution and computed solution
    plt.sca(axs[i, 0])
    if i == 0:
        plt.title(f"Exact")
    plt.annotate(
        f"t={time}s",
        xy=(0.5, 1),
        xytext=(-axs[i, 0].yaxis.labelpad - 3, 0),
        xycoords=axs[i, 0].yaxis.label,
        textcoords="offset points",
        size="large",
        ha="right",
        va="center",
    )
    CS1 = f.plot(c_exact, cmap="inferno")
    plt.sca(axs[i, 1])
    if i == 0:
        plt.title(f"Computed")
    CS2 = f.plot(computed_solution, cmap="inferno")

    plt.colorbar(CS1, ax=[axs[i, 0]], shrink=0.8)
    plt.colorbar(CS2, ax=[axs[i, 1]], shrink=0.8)

    axs[i, 0].sharey(axs[i, 1])
    plt.setp(axs[i, 1].get_yticklabels(), visible=False)

    for CS in [CS1, CS2]:
        CS.set_edgecolor("face")

    # define the profiles
    profiles = [
        {"start": (0.0, 0.0), "end": (1.0, 1.0)},
        {"start": (0.2, 0.8), "end": (0.7, 0.2)},
        {"start": (0.2, 0.6), "end": (0.8, 0.8)},
    ]

    # plot the profiles on the right subplot
    for profile in profiles:
        start_x, start_y = profile["start"]
        end_x, end_y = profile["end"]
        plt.sca(axs[i, 1])
        (l,) = plt.plot([start_x, end_x], [start_y, end_y])

        plt.sca(axs[i, 2])

        points_x_exact = np.linspace(start_x, end_x, num=30)
        points_y_exact = np.linspace(start_y, end_y, num=30)
        arc_length_exact = compute_arc_length(points_x_exact, points_y_exact)
        u_values = [c_exact(x, y) for x, y in zip(points_x_exact, points_y_exact)]

        points_x = np.linspace(start_x, end_x, num=100)
        points_y = np.linspace(start_y, end_y, num=100)
        arc_lengths = sorted(compute_arc_length(points_x, points_y))
        computed_values = [computed_solution(x, y) for x, y in zip(points_x, points_y)]

        (exact_line,) = plt.plot(
            arc_length_exact,
            u_values,
            color=l.get_color(),
            marker="o",
            linestyle="None",
            alpha=0.3,
        )
        (computed_line,) = plt.plot(arc_lengths, computed_values, color=l.get_color())

    plt.sca(axs[i, 2])
    plt.xlabel("Arc length")
    if i == 0:
        legend_marker = mpl.lines.Line2D(
            [],
            [],
            color="black",
            marker=exact_line.get_marker(),
            linestyle="None",
            label="Exact",
        )
        legend_line = mpl.lines.Line2D([], [], color="black", label="Computed")
        plt.legend(
            [legend_marker, legend_line],
            [legend_marker.get_label(), legend_line.get_label()],
        )
    plt.grid(alpha=0.3)
    plt.gca().spines[["right", "top"]].set_visible(False)

plt.show()

../../_images/3d73bb1f9247dc026d9d9fabc73fd903d03f8f71a429079ba4dadd25dd20b28d.png

Compute convergence rates#

It is also possible to compute how the numerical error decreases as we increase the number of cells. By iteratively refining the mesh, we find that the error exhibits a second order convergence rate. This is expected for this particular problem as first order finite elements are used.