festim-dev / festim-dev/FESTIM

[BUG]`HenrysBC` / `SievertsBC` fail when activation energy is exactly `0.0`

Open
#1,064 2 comments 0 reactions 0 assignees View on GitHub
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**

Image

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.