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

[hhy2022]

Describe the bug
When using HenrysBC and/or SievertsBC, setting the activation energy parameter (E_Hor 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

[RemDelaporteMathurin]

I've tracked down the issue to HenryBC.create_new_value_function, especially this line

return lambda T: henrys_law(T, self.H_0, self.E_H, self.pressure)

[RemDelaporteMathurin]

The issue is that if E_H is zero, pressure is not time dependent, and everything is a float, then we try to use a simple float inside an expression... A simple fix in your situation @hhy2022 is to replace HenryBC by

    F.FixedConcentrationBC(
        subdomain=left,
        species=H,
        value=pressure * H_0
    ),

Alternatively we could check if the final value inside FixedConcentrationBC.create_value is a float, then convert it to a fem.Constant

                    self.bc_expr = fem.Expression(
                        self.value(**kwargs),
                        helpers.get_interpolation_points(function_space.element),
                    )

What do you think @jhdark ?