Hi everyone, I’m trying to reproduce the results of the article “Influence of interface conditions on hydrogen transport studies” and I have some questions.
First, what do “OK TMAP” and “OK MUYI” in the source code mean respectively?
Second, I followed the source code, bi_material_flux.py, but it gave an error: AttributeError: module ‘festim.generic_simulation’ has no attribute ‘run’. Can anyone help me with this?
The detailed code is shown as follow:
Hey there
What version of Festim are you running? This may be causing the issue of why it’s not running. Regarding the inline comments, I have no idea, that will have to be a question for @remidm.
For future reference, you can copy code onto discourse using ``` markers before and after the code.
However, I have made a new version of this for the newest version of FESTIM (v1.2.1):
import festim as F
import numpy as np
chemical_pot_conervation = True
if chemical_pot_conervation:
folder = "bi_material/chemical_pot_conservation/"
else:
folder = "bi_material/concentration_continuity/"
my_model = F.Simulation(log_level=20)
# define mesh
size = 4e-03
vertices = np.linspace(0, 4e-03, num=100)
my_model.mesh = F.MeshFromVertices(vertices=vertices)
# define materials
tungsten = F.Material(
D_0=0.8165 * 2.9e-07,
E_D=0.39,
S_0=1.87e24,
E_S=1.04,
borders=[0, size / 2],
id=1,
)
copper = F.Material(
D_0=6.6e-07,
E_D=0.387,
S_0=3.14e24,
E_S=0.572,
borders=[size / 2, size],
id=2,
)
my_model.materials = F.Materials([tungsten, copper])
# define boundary conditions
my_model.boundary_conditions = [
F.DirichletBC(field="solute", value=1e20, surfaces=1),
F.DirichletBC(field="solute", value=0, surfaces=2),
]
# define temperature
my_model.T = F.Temperature(value=500)
# define exports
my_derived_quantities = F.DerivedQuantities(filename=folder + "derived_quantities.csv")
my_derived_quantities.derived_quantities = [
F.SurfaceFlux(field="solute", surface=1),
F.SurfaceFlux(field="solute", surface=2),
]
my_model.exports = F.Exports(
[
F.XDMFExport(
"solute",
label="solute",
folder=folder,
checkpoint=False,
mode=1,
),
my_derived_quantities,
]
)
# define solving parameters
my_model.dt = F.Stepsize(
initial_value=100,
stepsize_change_ratio=1.2,
dt_min=1e-8,
)
my_model.settings = F.Settings(
transient=True,
final_time=1e6,
absolute_tolerance=1e5,
relative_tolerance=1e-10,
maximum_iterations=20,
traps_element_type="DG",
)
# run simulation
my_model.initialise()
my_model.run()
The version I’m using is v1.2.1.
Okay, I see.
It works. Thank you so much!
Yeah this original script was for v0.7.1. There is a way to use that specific version:
pip install festim==0.7.1
However, I would still recommend to use a more recent version
Got it. Thank you very much.
Hi Xin, and welcome!
Oh dear that is taking me back to my PhD years… 
We initially took that problem from a benchmark project where we compared FESTIM to TMAP and also COMSOL.
These comments were from one of my colleagues checking that the parameter was the same as in TMAP or COMSOL (a colleague from China, Muyi, sent us some parameters, hence OK MUYI).
Nothing to worry about! 
Haha, I got it.
That is interesting.
Hi, good morning.
Can I ask that how the 2D geometry model (from “Influence of interface conditions on hydrogen transport studies”, Figure 5) was built in FESTIM?
@xinshen not sure what you are actually referring to:
Thank you very much for the step by step explanation. Your answer solved my question perfectly. 
@xinshen would you mind creating another topic for this please?
Sure.
Sorry. I should’ve done that.
Hi all, I am still reproducing the 2D case and I have a few question.
What does “traps_element_type”: “DG” mean in the line 197 of the code? (iter_2d_case.py, https://github.com/RemDelaporteMathurin/interface_conditions_paper/blob/d40ee0fb3daf0a7622a3e41c84b14014595cbab0/iter_2d_case.py)
Does the “boundary_conditions”: “type”: “dc” mean Dirichlet boundary condition? (line 145 to 147 of the iter_2d_case.py)
a) I am trying to use the new version (v1.2.1) to run the iter_2d_case based on the source code, but no derived_quantities file was successfully exported. I was wondering why? Here is my code.
import festim as F
my_model = F.Simulation(log_level=20)
import meshio
atom_density_W = 6.28e28
atom_density_Cu = 8.43e28
atom_density_CuCrZr = 2.6096e28
id_W = 8
id_Cu = 7
id_CuCrZr = 6
id_W_yUpper = 18
id_W_xlower = 19
id_CuCrZr_Surface = 16
def rhoCp_W(T):
return 5.15356e-6*T**3-8.30703e-2*T**2+5.98312e2*T+2.48160e6
def rhoCp_Cu(T):
return 1.68402e-4*T**3+6.14079e-2*T**2+4.67353e2*T+3.45899e6
def rhoCp_CuCrZr(T):
return -1.79134e-4*T**3+1.51383e-1*T**2+6.22091e2*T+3.46007e6
chemical_pot_conservation = True
if chemical_pot_conservation:
folder = "iter_2d_case/chemical_pot_conservation"
else:
folder = "iter_2d_case/concentration_continuity"
my_model.mesh = F.MeshFromXDMF(volume_file="bi_material/mesh_domains.xdmf", boundary_file="bi_material/mesh_boundaries.xdmf")
W = F.Material(
D_0=2.9e-7*0.8165,
E_D=0.39,
S_0=1.87e24,
E_S=1.04,
rho= rhoCp_W,
id=8,
)
Cu = F.Material(
D_0=6.6e-7,
E_D=0.387,
S_0=3.14e24,
E_S=0.572,
rho= rhoCp_Cu,
id=7,
)
CuCrZr = F.Material(
D_0=3.92e-7,
E_D=0.418,
S_0=4.28e23,
E_S=0.387,
rho= rhoCp_CuCrZr,
id=6,
)
my_model.materials = F.Materials([W, Cu, CuCrZr])
trap_1 = F.Trap(
E_k=0.39,
k_0=2.9e12*0.8165/atom_density_W,
E_p=1.2,
p_0=8.4e12,
density=5e-4*atom_density_W,
materials=W,
)
trap_2 = F.Trap(
E_k=0.387,
k_0=6.6e-7/(3.61e-10**2)/atom_density_Cu,
E_p=0.5,
p_0=7.98e13,
density=5e-5*atom_density_Cu,
materials=Cu,
)
trap_3 = F.Trap(
E_k=0.418,
k_0=3.92e-7/(3.61e-10**2)/atom_density_CuCrZr,
E_p=0.5,
p_0=7.98e13,
density=5e-5*atom_density_CuCrZr,
materials=CuCrZr,
)
trap_4 = F.Trap(
E_k=0.418,
k_0=3.92e-7/(3.61e-10**2)/atom_density_CuCrZr,
E_p=0.83,
p_0=7.98e13,
density=0.04*atom_density_CuCrZr,
materials=CuCrZr,
)
my_model.traps = [trap_1, trap_2, trap_3, trap_4]
my_model.boundary_conditions = [
F.DirichletBC(field="solute", value=1.223634016511513351e23, surfaces=id_W_yUpper),
F.DirichletBC(field="solute", value=0, surfaces=id_W_xlower),
F.RecombinationFlux(Kr_0=2.9e-14, E_Kr=1.92, order=2, surfaces=id_CuCrZr_Surface)
]
stationary_conditions = [
F.DirichletBC(field="temperature", value=1200, surfaces=id_W_yUpper),
F.DirichletBC(field="temperature",value=373, surfaces=id_CuCrZr_Surface)
]
my_model.dt = F.Stepsize(
initial_value=2.4e5,
stepsize_change_ratio=1.2,
t_stop=1e8,
max_stepsize=1e7,
dt_min=1e-8,
)
my_derived_quantities = F.DerivedQuantities(filename=folder + "derived_quantities.csv")
my_derived_quantities.derived_quantities = [
F.TotalVolume(field="retention", volume=6),
F.TotalVolume(field="retention", volume=7),
F.TotalVolume(field="retention", volume=8),
F.SurfaceFlux(field="solute", surface=18),
F.SurfaceFlux(field="solute", surface=16)
]
my_model.exports = F.Exports(
[
F.XDMFExport(
"retention",
label="retention",
folder=folder,
checkpoint=True,
filename="derived_quantities.xdmf"
),
my_derived_quantities,
]
)
all_timesteps = True
my_model.settings = F.Settings(
absolute_tolerance=1e11,
relative_tolerance=1e-10,
maximum_iterations=20,
)
if my_model is not None:
print(my_model.sources)
else:
my_model.initialise()
my_model.run()
This case has no error message, and also no exported files. 
b) If I change the last part of the code from a) to the way as below, the process finished with exit code 1. How can I improve this?
my_model.initialise()
my_model.run()
Error message is : AttributeError: ‘NoneType’ object has no attribute ‘sources’
Thank you for your prompt attention!
This means Discontinuous Galerkin. Using these elements are useful when the trap concentration is discontinuous in space. Otherwise you may observe oscillations (over and undershoots) also known as Gibbs phenomenon. Typically, when running multi-material simulations, I suggest setting the trap element type to “DG”.
Correct
Could you share the error message please?
Can you also share the error message?
Thank you so much!
I have added the error message to the last comment.
Thanks. Would you mind making a MWE (Minimal Working Example)? A MWE can reproduce the bug/error with the minimal amount of code and can be run easily (without any additional file if possible).
In your case, this MWE reproduces the error:
import festim as F
my_model = F.Simulation()
my_model.mesh = F.MeshFromVertices([0, 1, 2, 3])
my_model.materials = F.Material(D_0=1, E_D=0, id=1)
my_model.settings = F.Settings(
absolute_tolerance=1e-10,
relative_tolerance=1e-10,
transient=False,
)
my_model.initialise()
my_model.run()
Produces:
Traceback (most recent call last):
File "/home/remidm/FESTIM/mwe_2.py", line 16, in <module>
my_model.initialise()
File "/home/remidm/FESTIM/festim/generic_simulation.py", line 264, in initialise
self.attribute_source_terms()
File "/home/remidm/FESTIM/festim/generic_simulation.py", line 173, in attribute_source_terms
self.T.sources = []
^^^^^^^^^^^^^^
AttributeError: 'NoneType' object has no attribute 'sources
The problem is simply that the temperature attribute hasn’t been defined.
Adding my_model.T = F.Temperature(500) solves the issue.
For any additional issue, would you mind opening a new topic?
Oh I see, Thank you. I will try making a MWE in the future.
Adding my_model.T = F.Temperature(500) makes the case run, but another error comes out: ‘<’ not supported between instances of ‘int’ and ‘NoneType’. I will try to make a MWE to reproduce this error in a new topic if I can’t figure this problem out.
By the way, if I want to set up different temperature at different places like your code (line 162-175), how can I achieve that in the FESTIM v1.2.1?
Thanks for your suggestion. I will open a new topic if I have additional issues. I thought the question I asked today belongs to this topic, so I just added a comment under this topic. I am sorry for that.
Thank you very much!
Thanks!
The reason why I’m asking for seperate topics is because these issues are not strictly limited to the problem of reproducing the results from the paper but more general. In general, we prefer having small topics 
Thanks!
Haha I got it. No problem!
This is because you forgot to set a final time in the settings.