festim-dev / festim-dev/FESTIM
[BUG]`HenrysBC` / `SievertsBC` fail when activation energy is exactly `0.0`
- Dominant language
- Python
- Stars
- 135
- Forks
- 45
- Avg merge
- 3d 23h
- Merged PRs (30d)
- 14
Description
**Describe the bug**
When using `HenrysBC` and/or `SievertsBC`, setting the activation energy parameter (`E_H`or `E_S`) to exactly `0.0` causes the simulation to fail.
If the activation energy is set to a very small nonzero value (for example `1e-30`), the simulation runs normally.
Since an activation energy of zero is physically valid (Arrhenius factor = 1), the model is expected to behave identically for `0.0` and very small values.
**To Reproduce**
```
import numpy as np
from mpi4py import MPI
import dolfinx
import festim as F
COMM = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_interval(COMM, 20)
class LeftSurface(F.SurfaceSubdomain):
def locate_boundary_facet_indices(self, mesh):
return dolfinx.mesh.locate_entities_boundary(
mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0)
)
class RightSurface(F.SurfaceSubdomain):
def locate_boundary_facet_indices(self, mesh):
return dolfinx.mesh.locate_entities_boundary(
mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 1.0)
)
left = LeftSurface(id=1)
right = RightSurface(id=2)
model = F.HydrogenTransportProblemDiscontinuous()
model.mesh = F.Mesh(mesh, coordinate_system="cartesian")
mat = F.Material(
D_0=1.0,
E_D=0.0,
K_S_0=1.0,
E_K_S=0.0,
solubility_law="sievert",
)
vol = F.VolumeSubdomain(id=1, material=mat)
model.subdomains = [vol, left, right]
model.surface_to_volume = {left: vol, right: vol}
H = F.Species("H", subdomains=model.volume_subdomains)
model.species = [H]
T = 973.15
model.temperature = T
P_left = 1e5
E_H = 0.0
E_S = 0.0
model.boundary_conditions = [
F.HenrysBC(
subdomain=left,
species=H,
pressure=P_left,
H_0=1.0,
E_H=E_H,
),
# F.SievertsBC(
# subdomain=left,
# species=H,
# pressure=P_left,
# S_0=1.0,
# E_S=E_S,
# ),
F.FixedConcentrationBC(
subdomain=right,
species=H,
value=0.0,
),
]
model.settings = F.Settings(atol=1e-12, rtol=1e-12, transient=False)
model.initialise()
model.run()
```
**Expected behavior**
Setting activation energy to `0.0` should be valid and should not alter solver stability or model behaviour.
**Screenshots**
Contributor guide
No contributing guide indexed for this repository
Assessment
This issue has not been assessed yet.