Hello everyone,
I’ve been using festim 2.0b0 (up to date to the latest commit in the fenicsx branch) to run 1D simulations in a single-material slab. I was restarting simulations from checkpoints and using the “read_function_from_file” function to load the initial conditions.
I tried running similar simulations now for the multi-material case and it seems like the checkpointing is not implemented for these cases.
Here is the script that I’m using:
import festim as F
import numpy as np
from ufl import conditional, And, ge, lt
def linear_piecewise_ufl(t, times, values):
assert len(times) == len(values) and len(times) >= 2
expr = values[0]
for i in range(len(times) - 1):
t0, t1 = float(times[i]), float(times[i + 1])
v0, v1 = float(values[i]), float(values[i + 1])
slope = (v1 - v0) / (t1 - t0)
seg = v0 + slope * (t - t0)
expr = conditional(And(ge(t, t0), lt(t, t1)), seg, expr)
expr = conditional(ge(t, float(times[-1])), float(values[-1]), expr)
return expr
thick_W = 2e-3
thick_SS = 5e-3
vertices = np.concatenate([
np.linspace(0, thick_W, 300),
np.linspace(thick_W, thick_W+thick_SS, 400)[1:], # avoid duplicate node at thick_W
])
mesh = F.Mesh1D(vertices)
TUNGSTEN = F.Material(D_0=1.9e-7, E_D=0.2, K_S_0=4.52e21, E_K_S=0.3, name=“W”)
SS316L = F.Material(D_0=1.45e-6, E_D=0.59, K_S_0=1.38e23, E_K_S=0.0883, name=“SS”)
vol_W = F.VolumeSubdomain1D(id=1, borders=[0, thick_W], material=TUNGSTEN)
vol_SS = F.VolumeSubdomain1D(id=2, borders=[thick_W, thick_W + thick_SS], material=SS316L)
surf_in = F.SurfaceSubdomain1D(id=1, x=0)
surf_out = F.SurfaceSubdomain1D(id=2, x=thick_W + thick_SS)
interface_WSS = F.Interface(id=1, subdomains=[vol_W, vol_SS])
my_model = F.HydrogenTransportProblemDiscontinuous()
my_model.mesh = mesh
my_model.subdomains = [vol_W, vol_SS, surf_in, surf_out]
my_model.interfaces = [interface_WSS]
my_model.surface_to_volume = {surf_in: vol_W, surf_out: vol_SS}
mobile_T = F.Species(“T”, subdomains=[vol_W, vol_SS])
trapped_SS = F.Species(“trapped_SS”, mobile=False, subdomains=[vol_SS])
trapped_W = F.Species(“trapped_W”, mobile=False, subdomains=[vol_W])
ss_atom_density = 8.45e28
w_atom_density = 6.3e28
empty_trap_SS = F.ImplicitSpecies(name=“empty_trap_SS”, n=0.0008 * ss_atom_density, others=[trapped_SS])
empty_trap_W = F.ImplicitSpecies(name=“empty_trap_W”, n=1.0e-4 * w_atom_density, others=[trapped_W])
D0E, EDE = 1.45e-6, 0.59
D0W, EDW = 1.9e-7, 0.2
reaction_trap_1 = F.Reaction(
k_0=D0E / (3.5935e-102 * ss_atom_density), E_k=EDE,
p_0=8.53e13, E_p=0.7,
reactant=[mobile_T, empty_trap_SS], product=[trapped_SS], volume=vol_SS,
)
reaction_trap_W1 = F.Reaction(
k_0=D0W / (1.1e-102 * 6 * w_atom_density), E_k=EDW,
p_0=1e13, E_p=0.87,
reactant=[mobile_T, empty_trap_W], product=[trapped_W], volume=vol_W,
)
my_model.species = [mobile_T, trapped_SS, trapped_W]
my_model.reactions = [reaction_trap_1, reaction_trap_W1]
time_points = [0, 100, 200, 300]
power_points = [0, 1, 1, 0]
base_temp = 300.0
peak_temp_W = 1200.0
peak_temp_SS = 800.0
def T_profile(x, t):
P = linear_piecewise_ufl(t, time_points, power_points)
TL_W = base_temp + P * (peak_temp_W - base_temp)
TL_SS = base_temp + P * (peak_temp_SS - base_temp)
TR_W, TR_SS = TL_W, TL_SS
TW = TL_W + (TR_W - TL_W) * (x[0]) / thick_W
TSS = TL_SS + (TR_SS - TL_SS) * (x[0] - thick_W) / thick_SS
return conditional(lt(x[0], thick_W), TW, TSS)
my_model.temperature = T_profile
j0 = 1e20
PTVu, PL, PTVd = 50, 200, 50
RE, RP = 1.0, 4.5e-9
kb = 1.38e-23
eV = 1.602176634e-19
D0W, EDW = 1.9e-7, 0.2
K0W, EKW = 7.94e-17, -2.0
def P_num(t: float) → float:
return float(np.interp(t, time_points, power_points))
def Tsurf(t: float) → float:
return base_temp + P_num(t) * (peak_temp_W - base_temp)
def ion_flux_scalar(t: float) → float:
if t < PTVu:
return j0 * t / PTVu
elif t < PTVu + PL:
return j0
elif t < PTVu + PL + PTVd:
return j0 * (PTVu + PL + PTVd - t) / PTVd
else:
return 0.0
def surface_conc(t: float) → float:
Ts = Tsurf(t)
DE = D0W * np.exp(-EDW * eV / (kb * Ts))
KE = K0W * np.exp(-EKW * eV / (kb * Ts))
phi = max(ion_flux_scalar(t), 0.0)
return RE * phi * RP / DE + np.sqrt(phi / KE)
my_model.boundary_conditions = [
F.FixedConcentrationBC(subdomain=surf_in, species=mobile_T, value=lambda t: surface_conc(t)),
F.FixedConcentrationBC(subdomain=surf_out, species=mobile_T, value=lambda t: 0.0),
]
my_model.settings = F.Settings(
transient=True,
final_time=PTVu + PL + PTVd,
atol=1e11,
rtol=1e-9,
max_iterations=500,
stepsize=F.Stepsize(
1e-12,
growth_factor=1.1,
cutback_factor=0.1,
target_nb_iterations=5,
max_stepsize=1,
milestones=time_points,
),
)
my_model.exports = [
F.VTXSpeciesExport(“results/mobile_conc_W.bp”, field=mobile_T, subdomain=vol_W),
F.VTXSpeciesExport(“results/mobile_conc_SS.bp”, field=mobile_T, subdomain=vol_SS),
F.VTXSpeciesExport(“results/trapped_W.bp”, field=trapped_W, subdomain=vol_W),
F.VTXSpeciesExport(“results/trapped_SS.bp”, field=trapped_SS, subdomain=vol_SS),
F.VTXSpeciesExport(“results/mobile_conc_W.bp”, field=mobile_T, subdomain=vol_W, checkpoint=True, times=time_points),
F.VTXSpeciesExport(“results/mobile_conc_SS.bp”, field=mobile_T, subdomain=vol_SS, checkpoint=True, times=time_points),
F.VTXSpeciesExport(“results/trapped_W.bp”, field=trapped_W, subdomain=vol_W, checkpoint=True, times=time_points),
F.VTXSpeciesExport(“results/trapped_SS.bp”, field=trapped_SS, subdomain=vol_SS, checkpoint=True, times=time_points),
]
my_model.initialise()
my_model.run()
And the full error message:
Traceback (most recent call last):
File “/home/ITER/llealsa/multimaterial_dolfinx10/minimal.py”, line 150, in
my_model.initialise()
File “/home/ITER/llealsa/miniconda3/envs/festim-fenicsx/lib/python3.11/site-packages/festim/hydrogen_transport_problem.py”, line 1209, in initialise
self.initialise_exports()
File “/home/ITER/llealsa/miniconda3/envs/festim-fenicsx/lib/python3.11/site-packages/festim/hydrogen_transport_problem.py”, line 1701, in initialise_exports
raise NotImplementedError(
NotImplementedError: Export type <class ‘festim.exports.vtx.VTXSpeciesExport’> not implemented for mixed-domain approach
When commenting the checkpointing exports, the simulation runs without giving any errors.
I was wondering if there’s an alternative way to export checkpoints and load them as initial conditions in this current version.
Thank you very much for your help!