Error with HenrysBC and zero activation energy

User Support
1

Hi all!

I tried running a test problem using HenrysBC using a zero Henrys constant activation energy (E_H=0). Here is MWE:

MWE 
import festim as F
import numpy as np

model = F.HydrogenTransportProblem()

vertices = np.linspace(0, 10, 50)

model.mesh = F.Mesh1D(vertices)

material = F.Material(D_0=1, E_D=0)

left_surface = F.SurfaceSubdomain1D(id=1, x=0)
volume = F.VolumeSubdomain1D(id=2, borders=[0, 10], material=material)
model.subdomains = [left_surface, volume]

H = F.Species("H")
model.species = [H]

model.boundary_conditions = [
    F.HenrysBC(subdomain=left_surface, H_0=1, E_H=0, pressure=1, species=H)
]

model.temperature = 10

model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=30)

model.settings.stepsize = F.Stepsize(0.05)

model.initialise()
model.run()

which raises the following error

Error message 
Could not extract MPI communicator for Expression. Maybe you need to pass a communicator?

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[12], line 29
     25 model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=30)
     27 model.settings.stepsize = F.Stepsize(0.05)
---> 29 model.initialise()
     30 model.run()

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/festim/hydrogen_transport_problem.py:288, in HydrogenTransportProblem.initialise(self)
    285 self.create_implicit_species_value_fenics()
    287 self.define_temperature()
--> 288 self.define_boundary_conditions()
    289 self.convert_source_input_values_to_fenics_objects()
    290 self.convert_advection_term_to_fenics_objects()

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/festim/hydrogen_transport_problem.py:577, in HydrogenTransportProblem.define_boundary_conditions(self)
    570     if isinstance(bc, boundary_conditions.ParticleFluxBC):
    571         bc.create_value_fenics(
    572             mesh=self.mesh.mesh,
    573             temperature=self.temperature_fenics,
    574             t=self.t,
    575         )
--> 577 super().define_boundary_conditions()

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/festim/problem.py:116, in ProblemBase.define_boundary_conditions(self)
    114 for bc in self.boundary_conditions:
    115     if isinstance(bc, F.DirichletBCBase):
--> 116         form = self.create_dirichletbc_form(bc)
    117         self.bc_forms.append(form)

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/festim/hydrogen_transport_problem.py:595, in HydrogenTransportProblem.create_dirichletbc_form(self, bc)
    592 else:
    593     function_space_value = bc.species.collapsed_function_space
--> 595 bc.create_value(
    596     temperature=self.temperature_fenics,
    597     function_space=function_space_value,
    598     t=self.t,
    599 )
    601 # get dofs
    602 if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/festim/boundary_conditions/dirichlet_bc.py:246, in FixedConcentrationBC.create_value(self, function_space, temperature, t, K_S)
    242             kwargs["T"] = temperature
    244         # store the expression of the boundary condition
    245         # to update the value_fenics later
--> 246         self.bc_expr = fem.Expression(
    247             self.value(**kwargs),
    248             helpers.get_interpolation_points(function_space.element),
    249         )
    250         self.value_fenics.interpolate(self.bc_expr)
    252 # if K_S is provided, divide the value by K_S (change of variable method)

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/dolfinx/fem/function.py:131, in Expression.__init__(self, e, X, comm, form_compiler_options, jit_options, dtype)
    129 if comm is None:
    130     try:
--> 131         mesh = ufl.domain.extract_unique_domain(e).ufl_cargo()
    132         comm = mesh.comm
    133     except AttributeError:

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/ufl/domain.py:242, in extract_unique_domain(expr)
    240 def extract_unique_domain(expr):
    241     """Return the single unique domain expression is defined on or throw an error."""
--> 242     domains = extract_domains(expr)
    243     if len(domains) == 1:
    244         return domains[0]

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/ufl/domain.py:232, in extract_domains(expr)
    230 """Return all domains expression is defined on."""
    231 domainlist = []
--> 232 for t in traverse_unique_terminals(expr):
    233     domainlist.extend(t.ufl_domains())
    234 return sorted(
    235     join_domains(domainlist),
    236     key=lambda D: (D.topological_dimension(), D.ufl_cell(), D.ufl_id()),
    237 )

File ~/anaconda3/envs/vv-festim-report-env/lib/python3.13/site-packages/ufl/corealg/traversal.py:141, in traverse_unique_terminals(expr, visited)
    139 """Traverse unique terminals."""
    140 for op in unique_pre_traversal(expr, visited=visited):
--> 141     if op._ufl_is_terminal_:
    142         yield op

AttributeError: 'float' object has no attribute '_ufl_is_terminal_'
2

What version of FESTIM is this using?

3

Right, sorry

It is the 2.0a2 version.

4

Do you have the error with SievertsBC?

5

Yes, the same problem.

6

I made a FixedConcentrationBC version giving the same error:

import festim as F
import numpy as np
import ufl

model = F.HydrogenTransportProblem()

vertices = np.linspace(0, 10, 50)

model.mesh = F.Mesh1D(vertices)

material = F.Material(D_0=1, E_D=0)

left_surface = F.SurfaceSubdomain1D(id=1, x=0)
volume = F.VolumeSubdomain1D(id=2, borders=[0, 10], material=material)
model.subdomains = [left_surface, volume]

H = F.Species("H")
model.species = [H]


def arrhenius(T, x0, Ex):
    return x0 * ufl.exp(-Ex / F.k_B / T)


model.boundary_conditions = [
    F.FixedConcentrationBC(
        subdomain=left_surface, value=lambda T: arrhenius(T, 1, 0), species=H
    ),
]

model.temperature = 1000

model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=30)

model.settings.stepsize = F.Stepsize(0.05)

model.initialise()
model.run()
7

And here’s a pure dolfinx version:

import dolfinx
from mpi4py import MPI
import ufl

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("CG", 1))

T = dolfinx.fem.Function(V)


def arrhenius(T, x0, Ex):
    return x0 * ufl.exp(-Ex / T)


bc_expr = dolfinx.fem.Expression(
    arrhenius(T=T, x0=1, Ex=0), V.element.interpolation_points()
)
8 
bc_expr = dolfinx.fem.Expression(
    dolfinx.fem.Constant(mesh, 1.0),
    V.element.interpolation_points(),
)

This works though, so I guess you could use a FixedConcentrationBC instead and don’t bother with the arrhenius expression of K_H

9

@VVKulagin the workaround for now is:

solubility = 1.0
pressure = 1.0

model.boundary_conditions = [
    F.FixedConcentrationBC(
        subdomain=left_surface, value=solubility * pressure, species=H
    ),
]
1 Like