In FESTIM verification, surface kinetic model gives the code like that,
F.SurfaceKinetics(
k_sb=k_sb,
k_bs=k_bs,
lambda_IS=lambda_IS,
n_surf=n_surf,
n_IS=n_IS,
J_vs=J_vs,
surfaces=1,
initial_condition=exact_solution_cs(t=0),
t=F.t,
),
but this setting is only available in festim1.4, but not found in the festim2 api.
why kinetic surface model is absent for festim new version?
Can alll boundary conditions be set as the plasma dirven permeation case?
Hi @williamhayes123 and welcome to the FESTIM discourse 
Indeed, this surface kinetics model is not implemented in FESTIM2. The reason for that is that FESTIM1 used RealElements in legacy-fenics, which until very recently were not available in FEniCSx.
However, with the (very) recent developments coming up to v2.2 (currently available as a pre-release v2.2rc.2), we introduced manifolds which will allow to simulate surface kinetics.
You can find the relevant tutorials there: Manifold subdomains — FESTIM tutorial
Depending on your actual problem, you may not need the full kinetic surface model if the rates are fast enough and you could get away with dissociation/recombination fluxes instead.
Happy to discuss more
Thankyou for your reply. The updated tutorials:Manifold subdomains — FESTIM tutorial is important for solving the surface kinetics.
In kinetic surface model shown in Kinetic surface model — FESTIM V&V book, I try to rebuilt the model with FESTIM2.2 version with a manifold introduced. Although manifolds can define a surface subvolume with dim=1 for 2D simulation, subsurface storage term λ_IS *(∂c_m)/∂t is missing in the simulation.
To solve the problem, i read the code and try to understand what happen in solution process, finally found there is a simple way to reproduce the surfacekinetics in FESTIM2.2.
a newer api named surfacestorage should be added to the festim, surfacestroage act as a newer boundary conditon to give the subsurface storage term λ_IS *(∂c_m)/∂t
in folder festim/boundary_conditions
add surface_storage.py; (in folder festim/boundary_conditions)
add surfacestorage weak term integrals in the hydrogen_transport_problem.py, (after the bolck if isinstance(bc, boundary_conditions.ParticleFluxBC)
if isinstance(bc, boundary_conditions.SurfaceStorageBC):
target, measure, restriction = self.flux_bc_target(bc)
if subdomain == target:
v = bc.species.subdomain_to_test_function[subdomain]
c = bc.species.subdomain_to_solution[subdomain]
c_n = bc.species.subdomain_to_prev_solution[subdomain]
dt = self.subdomain_dt(subdomain)
dc_dt = (c - c_n) / dt
form_coupling += (
bc.coefficient
* self.restrict(dc_dt, restriction)
* self.restrict(v, restriction)
* measure
)
After finishing those preparing works, kinetic surface model can be rebuilt simply, the code is attached in the 8_1_KineticSurfaceModel.py
Kinetic Suface Model.zip (43.2 KB)
Besides, i also use another validation case Isothermal absorption of protium in titanium — FESTIM V&V book to examine the surface_storage.py.
I found if pressure varys during the simulation, present F.Enclousure can not directly transfer data to manifold because Enclosure needs a surface plane in present definition.
so I simply rebuilt the V&V case code.
It seems work well.
Your new SurfaceStorageBC term suggests a specific conservation check before relying on the rebuilt pressure-dependent case.
For constant lambda_IS and a fixed boundary Gamma, the kinetic equations imply the combined inventory
N = integral_Omega(c_m dx)
+ integral_Gamma((lambda_IS * trace(c_m) + c_s) ds).
The bulk loses lambda_IS*d(trace(c_m))/dt + J_bs - J_sb through that surface, while the adsorbed inventory gains J_bs - J_sb + J_vs. Adding the trace-storage inventory therefore cancels the exchange terms: dN/dt equals the volume source plus J_vs, minus flux leaving through other boundaries.
That tests the added weak term specifically, rather than just checking that the solve converges. For a closed regression, set the external input and remaining boundary fluxes to zero; the rear Dirichlet condition in the manufactured-solution example is not a closed boundary. Compare lambda_IS=0, positive lambda with storage included, and the same positive lambda with storage deliberately omitted as a negative control. Refine mesh and time step before choosing a balance tolerance.
If the enclosure is dynamically coupled, include its inventory too and distinguish internal exchange from imposed external input. Also keep units consistent: lambda_IS is a length when c_m is volume concentration, so lambda_IS*c_m has the same inventory-per-area units as c_s.
I derived this check from the linked kinetic equations; I have not run your attachments or verified an Enclosure/manifold API fix.
@williamhayes123 you could re-arrange the kinetic surface model equation λ_IS *(∂c_m)/∂t = f(…) into ∂c_m/∂t = f(…) / λ_IS
No?
EDIT: I misunderstood something. No this is not possible and yes a new BC object would be required I guess
@VVKulagin happy to have your opinion on this since you were the one to implement the kinetic surface model in FESTIM1
@williamhayes123 inspecting changes in a zip file is not ideal, could you open a PR?
@williamhayes123 another thought: could you try running the same validation case without this lambda term in the flux? I’m pretty sure you’d get the same result as lambda_IS is very small.
I can confirm that the lambda_IS term can be neglected. It wasn’t even present in the Pick and Sonnenberg model to be fair. It was introduced by Hodille or Schmidt.
Here’s the same V&V case in FESTIM2 with a custom problem class to add (or not) the lambda_IS term:
"""Isothermal absorption of H in Ti (Hirooka et al. 1981), FESTIM2 port of the
FESTIM V&V case validation/surface_model/hirooka.
Kinetic surface model on a point manifold at x=0 of a 1D half sample:
dc_s/dt = J_vs + J_bs - J_sb on Gamma (x=0)
-D dc_m/dx . n = lambda_IS dc_m/dt + J_bs - J_sb on Gamma
J_bs = K_bs lambda_abs c_m (1 - c_s/n_surf), lambda_abs = n_surf/n_IS
J_sb = K_sb c_s (1 - c_m/n_IS)
J_vs = J_ads(P_H2) - J_des
The gas pressure follows the absorbed inventory explicitly (lagged by one step),
as in the FESTIM1 V&V script. The lambda_IS storage term is optional so that
lambda_IS = 0 can be compared with lambda_IS = n_surf/n_IS.
"""
import argparse
import json
import numpy as np
from dolfinx import fem
import festim as F
# ------------------------- parameters (as in the V&V) -------------------------
N_A = 6.022e23
e = 1.602e-19
k_B = F.k_B # eV/K
M_H2 = 2.016e-3 / N_A # kg
A = 1e-2 * 1.3e-2 # m2, one face of the strip
L = 1e-3 / 2 # m, half thickness
V = A * L
X_max = 3
n_b = 9.4e4 * N_A # Ti atoms m-3
n_surf = X_max * 2.16e-5 * N_A # m-2
n_IS = X_max * n_b # m-3
N_Ti = n_b * V # Ti atoms in the half sample
lambda_Ti = n_surf / n_IS # m
lambda_abs = n_surf / n_IS # m, surface-subsurface distance (FESTIM 1.4)
# Wille & Davis (h-transport-materials, author="wille")
D0 = 9.0e-7 # m2/s
E_D = 12400 * 4.184 / 96485.33212 # eV (12400 cal/mol)
k_sb = 5.1413505e9 # s-1
E_sb = 1.0087995 # eV
k_bs = 1.0012013e10 # s-1, multiplies lambda_abs * c_m as in FESTIM 1.4
E_bs = 1.0518870 # eV
k_des = 3.4115671e-11 # m2/s
E_des = 5.6197034e-01 # eV
V_ch = 2.95e-3 # m3
P0 = 1.3e4 # Pa
def S0(T):
return 0.0143 * np.exp(1.99e3 / 96485.33212 / k_B / T)
class TiAbsorption(F.HydrogenTransportProblemDiscontinuous):
"""Adds the optional lambda_IS storage term and the explicit pressure update."""
def __init__(self, lambda_IS, exchange_bc, bulk, mobile, X, exports_map, **kw):
super().__init__(**kw)
self.lambda_IS = lambda_IS
self.exchange_bc = exchange_bc
self.bulk = bulk
self.mobile = mobile
self.X = X
self.exports_map = exports_map
def create_subdomain_formulation(self, subdomain):
super().create_subdomain_formulation(subdomain)
if self.lambda_IS == 0 or subdomain is not self.bulk:
return
target, measure, restriction = self.flux_bc_target(self.exchange_bc)
c = self.mobile.subdomain_to_solution[subdomain]
c_n = self.mobile.subdomain_to_prev_solution[subdomain]
v = self.mobile.subdomain_to_test_function[subdomain]
dc_dt = (c - c_n) / self.subdomain_dt(subdomain)
subdomain.F += (
self.lambda_IS
* self.restrict(dc_dt, restriction)
* self.restrict(v, restriction)
* measure
)
def inventory(self):
"""H per unit area in the half sample: bulk + adsorbed + subsurface storage."""
ex = self.exports_map
return (
ex["bulk"].value + ex["surface"].value + self.lambda_IS * ex["trace"].value
)
def iterate(self):
super().iterate()
# number of H atoms absorbed by the whole strip (two faces)
self.X.value = 2 * A * self.inventory()
def run(T0, lambda_IS, nx=1000):
K_bs = k_bs * np.exp(-E_bs / k_B / T0)
K_sb = k_sb * np.exp(-E_sb / k_B / T0)
K_des = k_des * np.exp(-E_des / k_B / T0)
gas_factor = 2 * S0(T0) / np.sqrt(2 * np.pi * M_H2 * k_B * T0 * e)
mesh = F.Mesh1D(np.linspace(0, L, nx))
X = fem.Constant(mesh.mesh, 0.0) # absorbed H atoms (whole strip)
def P_H2(X):
return k_B * T0 * e / V_ch * (P0 * V_ch / (k_B * T0 * e) - X / 2)
def J_bs(c_m, c_s):
return K_bs * lambda_abs * c_m * (1 - c_s / n_surf)
def J_sb(c_m, c_s):
return K_sb * c_s * (1 - c_m / n_IS)
def J_vs(c_s):
J_ads = gas_factor * (1 - c_s / n_surf) ** 2 * P_H2(X)
J_des = 2 * K_des * c_s**2
return J_ads - J_des
material = F.Material(D_0=D0, E_D=E_D)
bulk = F.VolumeSubdomain1D(id=1, borders=[0, L], material=material)
gamma = F.VolumeSubdomain(
id=2, dim=0, material=material, locator=lambda x: np.isclose(x[0], 0)
)
H_b = F.Species("H_bulk", subdomains=[bulk])
H_s = F.Species("H_surface", mobile=False, subdomains=[gamma])
dependencies = {"c_m": H_b, "c_s": H_s}
exchange = F.ParticleFluxBC(
gamma,
species=H_b,
value=lambda c_m, c_s: J_sb(c_m, c_s) - J_bs(c_m, c_s),
species_dependent_value=dependencies,
)
surface_source = F.ParticleSource(
volume=gamma,
species=H_s,
value=lambda c_m, c_s: J_vs(c_s) + J_bs(c_m, c_s) - J_sb(c_m, c_s),
species_dependent_value=dependencies,
)
exports_map = {
"bulk": F.TotalVolume(field=H_b, volume=bulk),
"surface": F.TotalVolume(field=H_s, volume=gamma),
"trace": F.TotalSurface(field=H_b, surface=gamma),
}
model = TiAbsorption(
lambda_IS=lambda_IS,
exchange_bc=exchange,
bulk=bulk,
mobile=H_b,
X=X,
exports_map=exports_map,
mesh=mesh,
subdomains=[bulk, gamma],
species=[H_b, H_s],
boundary_conditions=[exchange],
sources=[surface_source],
temperature=T0,
exports=list(exports_map.values()),
settings=F.Settings(
atol=1e6,
rtol=1e-8,
max_iterations=50,
final_time=25 * 60,
stepsize=F.Stepsize(
initial_value=1e-3,
growth_factor=1.1,
cutback_factor=0.9,
target_nb_iterations=4,
max_stepsize=5,
),
),
)
model.show_progress_bar = False
model.initialise()
model.run()
ex = exports_map
t = np.array(ex["bulk"].t)
bulk_inv = np.array(ex["bulk"].data)
surf_inv = np.array(ex["surface"].data)
trace = np.array(ex["trace"].data)
return {
"t": t.tolist(),
"H_Ti_bulk_surface": ((bulk_inv + surf_inv) / (n_b * L)).tolist(),
"H_Ti_storage": (lambda_IS * trace / (n_b * L)).tolist(),
"c_s": surf_inv.tolist(),
"c_m0": trace.tolist(),
}
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--T", type=float, required=True, help="temperature in K")
parser.add_argument("--lam", type=float, required=True, help="lambda_IS in m")
parser.add_argument("--nx", type=int, default=1000)
parser.add_argument("--out", required=True)
args = parser.parse_args()
result = run(args.T, args.lam, args.nx)
result.update(T=args.T, lambda_IS=args.lam, nx=args.nx)
with open(args.out, "w") as f:
json.dump(result, f)
"""Compare lambda_IS = 0 with lambda_IS = n_surf/n_IS and with Hirooka's data."""
import json
import sys
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
here = Path(__file__).parent
data_dir = Path(sys.argv[1]) # V-V-report/.../hirooka/content_data
temperatures = [723, 773, 823, 873, 923]
def load(T, tag):
r = json.loads((here / "results" / f"T{T}_lam{tag}.json").read_text())
t = np.array(r["t"]) / 60
inventory = np.array(r["H_Ti_bulk_surface"]) + np.array(r["H_Ti_storage"])
return t, inventory, np.array(r["H_Ti_storage"])
fig, (ax, ax_diff) = plt.subplots(
1, 2, figsize=(12, 5), gridspec_kw={"width_ratios": [3, 2]}
)
colors = plt.cm.viridis(np.linspace(0, 0.9, len(temperatures)))
print(
f"{'T (C)':>6} {'H/Ti(25min)':>12} {'max|dH/Ti|':>11} {'max rel':>9} "
f"{'max storage':>12} {'RMSE lam=0':>11} {'RMSE lam_Ti':>12}"
)
for T, color in zip(temperatures, colors):
t0, x0, _ = load(T, "0")
t1, x1, storage = load(T, "Ti")
x1_on_t0 = np.interp(t0, t1, x1)
diff = x1_on_t0 - x0
mask = t0 > 0.5 # relative difference once there is something absorbed
rel = np.max(np.abs(diff[mask]) / x0[mask])
exp = np.loadtxt(data_dir / f"{T - 273}.csv", delimiter=",", skiprows=1)
rmse0 = np.sqrt(np.mean((np.interp(exp[:, 0], t0, x0) - exp[:, 1]) ** 2))
rmse1 = np.sqrt(np.mean((np.interp(exp[:, 0], t1, x1) - exp[:, 1]) ** 2))
print(
f"{T - 273:>6} {x0[-1]:>12.4f} {np.max(np.abs(diff)):>11.2e} {rel:>9.2e} "
f"{storage.max():>12.2e} {rmse0:>11.4f} {rmse1:>12.4f}"
)
ax.plot(t0, x0, color=color, lw=2.5, label=f"{T - 273} °C, λ_IS = 0")
ax.plot(t1, x1, color="k", lw=1, ls="--")
ax.plot(exp[:, 0], exp[:, 1], "o", color=color, mfc="none", ms=7)
ax_diff.plot(t0, np.abs(diff), color=color, label=f"{T - 273} °C")
ax.plot([], [], color="k", lw=1, ls="--", label="λ_IS = n_surf/n_IS")
ax.plot([], [], "o", color="gray", mfc="none", label="Hirooka et al. (exp.)")
ax.set_xlim(0, 25)
ax.set_ylim(0, 0.8)
ax.set_xlabel("Time (min)")
ax.set_ylabel("Content (H/Ti)")
ax.legend(fontsize=8, ncol=2)
ax.set_title("FESTIM2, kinetic surface model on a point manifold")
ax_diff.set_yscale("log")
ax_diff.set_xlim(0, 25)
ax_diff.set_xlabel("Time (min)")
ax_diff.set_ylabel("|H/Ti(λ_IS = n_surf/n_IS) − H/Ti(λ_IS = 0)|")
ax_diff.set_title("Effect of the storage term")
ax_diff.legend(fontsize=8)
fig.tight_layout()
fig.savefig(here / "hirooka_lambda_comparison.png", dpi=150)
print("saved", here / "hirooka_lambda_comparison.png")
[Sorry these scripts aren’t super reproducible because they use local data on my workspace, but this is just to give you an idea of how to achieve this]
@williamhayes123 @VVKulagin if you still want to use the lambda_IS term then you can use the existing machinery:
"""Hirooka H-Ti case with stock FESTIM2 classes only, keeping lambda_IS.
The only non-standard step is the closed-chamber pressure, updated in a
hand-written time loop instead of a subclassed iterate().
"""
import json
from pathlib import Path
import numpy as np
from dolfinx import fem
import festim as F
import hirooka_HTi as H # parameters only, from the same folder
T0 = 873.0
lam = H.lambda_Ti # subsurface storage length, kept this time
K_bs = H.k_bs * np.exp(-H.E_bs / H.k_B / T0)
K_sb = H.k_sb * np.exp(-H.E_sb / H.k_B / T0)
K_des = H.k_des * np.exp(-H.E_des / H.k_B / T0)
gas = 2 * H.S0(T0) / np.sqrt(2 * np.pi * H.M_H2 * H.k_B * T0 * H.e)
mesh = F.Mesh1D(np.linspace(0, H.L, 1000))
X = fem.Constant(mesh.mesh, 0.0) # H atoms absorbed by the whole strip
def J_bs(c_m, c_s):
return K_bs * H.lambda_abs * c_m * (1 - c_s / H.n_surf)
def J_sb(c_m, c_s):
return K_sb * c_s * (1 - c_m / H.n_IS)
def J_vs(c_s):
P = H.k_B * T0 * H.e / H.V_ch * (H.P0 * H.V_ch / (H.k_B * T0 * H.e) - X / 2)
return gas * (1 - c_s / H.n_surf) ** 2 * P - 2 * K_des * c_s**2
mat = F.Material(D_0=H.D0, E_D=H.E_D)
bulk = F.VolumeSubdomain1D(id=1, borders=[0, H.L], material=mat)
gamma = F.VolumeSubdomain(
id=2, dim=0, material=mat, locator=lambda x: np.isclose(x[0], 0)
)
H_b = F.Species("H_bulk", subdomains=[bulk])
H_s = F.Species("H_surface", mobile=False, subdomains=[gamma])
deps = {"c_m": H_b, "c_s": H_s}
def exchange_value(c_m, c_s):
# evaluated during model.initialise(), when dt and the previous solution exist
c_n = H_b.subdomain_to_prev_solution[bulk]
return J_sb(c_m, c_s) - J_bs(c_m, c_s) - lam * (c_m - c_n) / model.dt
N_bulk = F.TotalVolume(field=H_b, volume=bulk)
N_surf = F.TotalVolume(field=H_s, volume=gamma)
c_m0 = F.TotalSurface(field=H_b, surface=gamma) # c_m at x=0
model = F.HydrogenTransportProblemDiscontinuous(
mesh=mesh,
subdomains=[bulk, gamma],
species=[H_b, H_s],
boundary_conditions=[
F.ParticleFluxBC(
gamma,
species=H_b,
value=exchange_value,
species_dependent_value=deps,
)
],
sources=[
F.ParticleSource(
volume=gamma,
species=H_s,
value=lambda c_m, c_s: J_vs(c_s) + J_bs(c_m, c_s) - J_sb(c_m, c_s),
species_dependent_value=deps,
)
],
temperature=T0,
exports=[N_bulk, N_surf, c_m0],
settings=F.Settings(
atol=1e6,
rtol=1e-8,
max_iterations=50,
final_time=25 * 60,
stepsize=F.Stepsize(
initial_value=1e-3,
growth_factor=1.1,
cutback_factor=0.9,
target_nb_iterations=4,
max_stepsize=5,
),
),
)
model.show_progress_bar = False
model.initialise()
while model.t.value < model.settings.final_time:
model.iterate()
X.value = 2 * H.A * (N_bulk.value + N_surf.value + lam * c_m0.value)
t = np.array(N_bulk.t)
H_Ti = (np.array(N_bulk.data) + np.array(N_surf.data) + lam * np.array(c_m0.data)) / (H.n_b * H.L)
ref = json.loads((Path(__file__).parent / "results/T873_lamTi.json").read_text())
t_ref = np.array(ref["t"])
H_ref = np.array(ref["H_Ti_bulk_surface"]) + np.array(ref["H_Ti_storage"])
print("same time steps:", np.array_equal(t, t_ref), len(t), len(t_ref))
print("H/Ti at 25 min: no class %.6f, subclass %.6f" % (H_Ti[-1], H_ref[-1]))
print("max |difference|:", np.max(np.abs(np.interp(t_ref, t, H_Ti) - H_ref)))
# plot both runs and the measurements on the same graph
import matplotlib # noqa: E402
matplotlib.use("Agg")
import matplotlib.pyplot as plt # noqa: E402
here = Path(__file__).parent
exp = np.loadtxt(here / f"exp_data/{T0 - 273:.0f}.csv", delimiter=",", skiprows=1)
fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(t / 60, H_Ti, lw=3, alpha=0.6, label="stock classes + closure")
ax.plot(t_ref / 60, H_ref, "k--", lw=1.2, label="custom subclass")
ax.plot(exp[:, 0], exp[:, 1], "o", mfc="none", color="tab:red", label="Hirooka et al.")
ax.set_xlim(0, 25)
ax.set_ylim(0, 0.8)
ax.set_xlabel("Time (min)")
ax.set_ylabel("Content (H/Ti)")
ax.set_title(f"T = {T0 - 273:.0f} °C, λ_IS = n_surf/n_IS")
ax.legend()
fig.tight_layout()
fig.savefig(here / "no_class_lambda.png", dpi=150)
print("saved", here / "no_class_lambda.png")
The trick is to add the lambda_IS dc/dt term in the model with a custom particle flux BC:
def exchange_value(c_m, c_s):
# evaluated during model.initialise(), when dt and the previous solution exist
c_n = H_b.subdomain_to_prev_solution[bulk]
return J_sb(c_m, c_s) - J_bs(c_m, c_s) - lam * (c_m - c_n) / model.dt
deps = {"c_m": H_b, "c_s": H_s}
F.ParticleFluxBC(
gamma,
species=H_b,
value=exchange_value,
species_dependent_value=deps,
)
Glad to discuss the problem.
I check the result between lambda_IS term on and off.
In isothermal absorption of protium in titanium case, the lambda_IS term can be negelected.
In Kinetic surface model — FESTIM V&V book case,
a) positive lambda with storage included
b) positive lamdda with storage neglated:
c) lambda_IS=0
In Kinetic surface model — FESTIM V&V book, I thought that lambda_IS means a sub-surface layer with the physical assumpution of limited exchange velocity between subsurface layer and surface layer which results in the accumlation in the subsurface layer.
In boundary, FESTIM particlefluxbc adds particle flux into the integrals to solve the bulk iteration. It forces -D(grad c dot n)=J_bs-J_sb+lambda_IS*(∂c_m)/∂t, if lambda_IS=0,.no subsurface exits in the bulk, and FESTIM calculate the -D(grad c dot n)=J_bs-J_sb which is exactly the exchange term between bulk and surface. but when lamda_IS is not 0, the surface accumulation at boundary will exist, thus why a bulit-in surfacestorage.py can give the right result. The larger the lamda_IS is,the surfacestorage are more important in the equations.
In fact, I think everything should be approached from a physics perspective，the J_bs-J_sb physcially describes the surface and bulk exchange, and FESTIM can put particle flux coupled into the bulk equation. The physics is described directly. There is no need to add a subsurface to the equation.
@remidm, @williamhayes123, @nickharris808, Hello!
I would like to contribute to the discussion.
A more general idea is to distinguish three hydrogen populations: adsorbed hydrogen at the surface, hydrogen in subsurface sites (which may differ from bulk sites), and hydrogen in the bulk. Such models have been investigated in several recent papers (Ono et al. (2019), Hamamoto et al. (2020), Shimohata et al. (2021) ).
The governing equations are:
where q_{ss} is the subsurface inventory per unit area.
If we assume that the subsurface concentration equals the bulk concentration at the interface in 1D-0D coupling (c_{ss}=c_m(0,t) and q_{ss}=\lambda_{\mathrm{IS}}c_{ss}), the subsurface balance reduces to
This is precisely the dynamic boundary condition used in FESTIM 1.4.
To illustrate the role of these populations, here is a small test based on Remi’s solution above (though in 2D geometry, since 1D+0D coupling doesn’t work for me with FESTIM2.2rc2), comparing three formulations:
import matplotlib.pyplot as plt
import numpy as np
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import festim as F
D = 0.25
n_surf = 1.0
n_IS = 3.0
k_sb = 1.4
k_bs = 0.65
J_vs = 0.4
final_time = 4.0
nx = 60
dx = 1 / nx
lambda_values = (1e-3, 1)
modes = ("full", "storage", "no_storage")
def j_sb(c_s, c_ss):
return k_sb * c_s * (1 - c_ss / n_IS)
def j_bs(c_s, c_ss):
return k_bs * c_ss * (1 - c_s / n_surf)
def run(mode, lambda_IS):
mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, nx, 2))
material = F.Material(D_0=D, E_D=0.0)
bulk = F.VolumeSubdomain(
id=1,
material=material,
locator=lambda x: np.ones(x.shape[1], dtype=bool),
)
surface = F.VolumeSubdomain(
id=2,
dim=1,
material=material,
locator=lambda x: np.isclose(x[0], 0.0),
)
c_m = F.Species("c_m", subdomains=[bulk])
c_s = F.Species("c_s", mobile=False, subdomains=[surface])
species = [c_m, c_s]
if mode == "full":
q_ss = F.Species("q_ss", mobile=False, subdomains=[surface])
species.append(q_ss)
def j_dif(c_bulk, q):
return D / dx * (q / lambda_IS - c_bulk)
bcs = [
F.ParticleFluxBC(
subdomain=surface,
species=c_m,
value=lambda c_bulk, q: j_dif(c_bulk, q),
species_dependent_value={"c_bulk": c_m, "q": q_ss},
)
]
sources = [
F.ParticleSource(
volume=surface,
species=c_s,
value=lambda c_surface, q: (
J_vs
+ j_bs(c_surface, q / lambda_IS)
- j_sb(c_surface, q / lambda_IS)
),
species_dependent_value={"c_surface": c_s, "q": q_ss},
),
F.ParticleSource(
volume=surface,
species=q_ss,
value=lambda c_bulk, c_surface, q: (
j_sb(c_surface, q / lambda_IS)
- j_bs(c_surface, q / lambda_IS)
- j_dif(c_bulk, q)
),
species_dependent_value={"c_bulk": c_m, "c_surface": c_s, "q": q_ss},
),
]
else:
def surface_flux(c_bulk, c_surface):
value = j_sb(c_surface, c_bulk) - j_bs(c_surface, c_bulk)
if mode == "storage":
old_bulk = c_m.subdomain_to_prev_solution[bulk]
value -= lambda_IS * (c_bulk - old_bulk) / model.dt
return value
bcs = [
F.ParticleFluxBC(
subdomain=surface,
species=c_m,
value=surface_flux,
species_dependent_value={"c_bulk": c_m, "c_surface": c_s},
)
]
sources = [
F.ParticleSource(
volume=surface,
species=c_s,
value=lambda c_bulk, c_surface: (
J_vs + j_bs(c_surface, c_bulk) - j_sb(c_surface, c_bulk)
),
species_dependent_value={"c_bulk": c_m, "c_surface": c_s},
)
]
N_bulk = F.TotalVolume(field=c_m, volume=bulk)
N_surf = F.TotalVolume(field=c_s, volume=surface)
c_m0 = F.TotalSurface(field=c_m, surface=surface)
exports = [N_bulk, N_surf, c_m0]
if mode == "full":
N_ss = F.TotalVolume(field=q_ss, volume=surface)
exports.append(N_ss)
model = F.HydrogenTransportProblemDiscontinuous(
mesh=mesh,
subdomains=[bulk, surface],
species=species,
boundary_conditions=bcs,
sources=sources,
temperature=300,
exports=exports,
settings=F.Settings(
atol=1e-9,
rtol=1e-9,
max_iterations=40,
final_time=final_time,
stepsize=F.Stepsize(
initial_value=1e-3,
growth_factor=1.1,
cutback_factor=0.5,
target_nb_iterations=4,
max_stepsize=0.05,
),
),
)
model.show_progress_bar = False
model.initialise()
model.run()
t = np.asarray(N_bulk.t)
cs = np.asarray(N_surf.data)
cm0 = np.asarray(c_m0.data)
q = (
np.asarray(N_ss.data)
if mode == "full"
else lambda_IS * cm0
if mode == "storage"
else np.zeros_like(t)
)
nm = np.asarray(N_bulk.data)
total = nm + cs + q
err = np.max(np.abs(total - J_vs * t))
print(f"lambda_IS={lambda_IS:g}, {mode}: mass-balance error = {err:.3e}")
return t, cs, cm0, q, nm
fig, axes = plt.subplots(2, 4, figsize=(14, 6))
labels = (
"Surface inventory",
"Bulk boundary concentration",
"Subsurface inventory",
"Bulk inventory",
)
for row, lambda_IS in enumerate(lambda_values):
results = {mode: run(mode, lambda_IS) for mode in modes}
t, cs, cm0, _, _ = results["storage"]
t0, cs0, cm00, _, _ = results["no_storage"]
print(
f"lambda_IS={lambda_IS:g}: "
f"max |Delta c_s|={np.max(np.abs(cs - np.interp(t, t0, cs0))):.3e}, "
f"max |Delta c_m(0)|={np.max(np.abs(cm0 - np.interp(t, t0, cm00))):.3e}"
)
for mode, (time, *curves) in results.items():
for ax, values in zip(axes[row], curves):
ax.plot(time, values, label=mode)
for ax, label in zip(axes[row], labels):
ax.set(xlabel="Time", ylabel=label)
axes[row, 0].set_title(f"lambda_IS = {lambda_IS:g}")
axes[0, 0].legend()
fig.tight_layout()
fig.savefig("surface_kinetics_three_models.png", dpi=150)
plt.show()
The results show that the full model and the reduced formulation with subsurface storage give nearly identical solutions in this test. Omitting the storage term changes the distribution of hydrogen between the surface and the bulk, while the total inventory remains conserved. Of course, the magnitude of the effect depends on the problem parameters.
I hope you find this helpful.
@williamhayes123 in the verification case of the V&V book, the lambda_IS term is purposefully exagerated.
Either way, in my previous comment I showed how one can still include this term in FESTIM 2 to have the full expression. Let me know if that doesn’t solve your problem!
It works. the result are equal to the result with surfacestorage included Thanks a lot.
