Pre-loaded semi-infinite slab

🏷 Tags: 1D MES transient

This verification case from case Val-1c of TMAP7’s V&V report [3] consists of a semi-infinite slab with no traps under a constant concentration \(c_0\) on the first \(10 \ \mathrm{m}\) of the slab. The concentration at the boundaries is set to zero.

FESTIM Code

import festim as F

from dolfinx.geometry import bb_tree, compute_colliding_cells, compute_collisions_points


class PointValue(F.VolumeQuantity):
    def __init__(self, field, volume, x0, filename=None):
        super().__init__(field, volume, filename)
        self.x0 = x0

    def compute(self):
        u = self.field.solution
        mesh = u.function_space.mesh
        tree = bb_tree(mesh, mesh.geometry.dim)
        cell_candidates = compute_collisions_points(tree, self.x0)
        cell = compute_colliding_cells(mesh, cell_candidates, self.x0).array
        assert len(cell) > 0
        first_cell = cell[0]

        self.value = self.field.solution.eval(self.x0, first_cell)
        self.data.append(self.value)

class Profile(F.Profile1DExport):
    def __init__(self, field, volume, times=None):
        super().__init__(field=field, subdomain=volume, times=times)
        self.data = []   
        self.t = []      
        self.x = None    

    def compute(self, *args, **kwargs):
        super().compute(*args, **kwargs)

        last_profile = self.data[-1].copy()
        self.data[-1] = last_profile     

/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

Show code cell content

Hide code cell content

import festim as F
import numpy as np
import ufl
from matplotlib import pyplot as plt

preloaded_length = 10  # m
C_0 = 1  # atom m^-3
D = 1  # 1 m^2 s^-1

model = F.HydrogenTransportProblem()

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

### Mesh Settings ###
vertices = np.concatenate(
    [
        np.linspace(0, 10, 400),
        np.linspace(10, 100, 1000),
    ]
)

model.mesh = F.Mesh1D(vertices)

material = F.Material(D_0=D, E_D=0)

left_boundary = F.SurfaceSubdomain1D(id=1, x=0)
volume = F.VolumeSubdomain1D(id=1, borders=[0, 100], material=material)

model.subdomains = [left_boundary, volume]

model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H)
]

initial_concentration = lambda x: ufl.conditional(x[0] <= preloaded_length, C_0, 0)
model.initial_conditions = [F.InitialConcentration(value=initial_concentration, species=H, volume=volume)]


model.temperature = 500  # ignored in this problem


test_points = [0.5, preloaded_length, 12]  # m
profile_times = [0.1] + np.linspace(0, 100, num=10).tolist()[1:]


profile_export = F.Profile1DExport(field=H, times=profile_times)

model.exports = [
    PointValue(field=H, volume=volume, x0=np.array([v, 0, 0])) for v in test_points
] + [profile_export
] 


model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=100)
model.settings.stepsize = F.Stepsize(
    initial_value=0.01,
    growth_factor=1.1,
    cutback_factor=0.9,
    target_nb_iterations=10,
    milestones=profile_times,
)

model.initialise()
model.run()

Comparison with exact solution

This is a comparison of the computed concentration profiles at different times with the exact analytical solution (shown in dashed lines).

The analytical solution is given by

\[ c(x, t) = \frac{c_0}{2}\left[ 2\mathrm{erf}\left(\frac{x}{2 \sqrt{Dt}}\right) - \mathrm{erf}\left(\frac{x - h}{2 \sqrt{Dt}}\right) - \mathrm{erf}\left(\frac{x + h}{2 \sqrt{Dt}}\right) \right] \]

where \(h\) is the thickness of the pre-loaded region.

Warning

There is an existing bug with ProfileExport1D class. See (this associated issue for more details)[https://github.com/festim-dev/FESTIM/issues/1046].

Show code cell source

Hide code cell source

from matplotlib import cm
from matplotlib.colors import Normalize
from scipy.special import erf


def exact_solution(x, t):
    sqrt_term = np.sqrt(4 * D * t)
    return (
        C_0
        / 2
        * (
            2 * erf(x / sqrt_term)
            - erf((x - preloaded_length) / sqrt_term)
            - erf((x + preloaded_length) / sqrt_term)
        )
    )


profile_export = model.exports[-1]
times = np.array(profile_export.t)

profiles_1d = [np.ravel(p) for p in profile_export.data]
data = np.stack(profiles_1d, axis=0)

x = profile_export.x
order = np.argsort(x)
x = x[order]
data = data[:, order]

norm = Normalize(vmin=times.min(), vmax=times.max())
cmap = cm.viridis

plt.figure()


for i, t in enumerate(profile_times):
    label = "exact" if i == 0 else ""
    y_name = f"t{t:.2e}s".replace("+", "").replace("-", "").replace(".", "")

    idx = np.argmin(np.abs(times - t))
    t = times[idx]
    y = data[idx, :]

    exact_y = exact_solution(x, t)

    plt.plot(x, exact_y, linestyle="dashed", color="tab:grey", linewidth=3, label=label)
    plt.plot(x, y, color=cmap(norm(t)))


plt.xlabel("x")
plt.ylabel("Concentration (atom / m^3)")

# Add colorbar
sm = cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])  # You can also set this to the range of your data
plt.colorbar(sm, label="Time (s)", ax=plt.gca())

plt.legend()
plt.show()

The results can also be compared with the results obtained by TMAP7.

Show code cell source

Hide code cell source

fig, axs = plt.subplots(1, len(test_points), sharey=True)

for i, x in enumerate(test_points):
    plt.sca(axs[i])

    # plotting computed data
    computed_solution = model.exports[i].data
    t = np.array(model.exports[i].t)
    plt.plot(t, computed_solution, label="FESTIM", linewidth=3)

    # plotting exact solution
    plt.plot(t, exact_solution(x, t), label="Exact", color="orange", linestyle="--")

    # plotting TMAP data
    tmap_data = np.genfromtxt(
        f"./tmap_1c_data/tmap_point_data_{i}.txt", delimiter=" ", names=True
    )
    tmap_t = tmap_data["t"]
    tmap_solution = tmap_data["tmap"]
    plt.scatter(tmap_t, tmap_solution, label="TMAP7", color="purple")

    plt.title(f"x={x} m")
    if i == 0:
        plt.ylabel("Concentration (atom / m$^3$)")
    plt.xlabel("t (s)")

fig.tight_layout()
plt.legend()
plt.show()