BUG: c=0 trick for steady state problems

[RemDelaporteMathurin]

In HydrogenTransportProblem we want apply c = 0 to the formulation in steady state in subdomains where immobile species are not involved in a reaction. We do this here by looking in which subdomains each immobile species isn't defined:

FESTIM/src/festim/hydrogen_transport_problem.py

Lines 683 to 699 in 013937a

	# check if each species is defined in all volumes
	if not self.settings.transient:
	for spe in self.species:
	# if species mobile, already defined in diffusion term
	if not spe.mobile:
	not_defined_in_volume = self.volume_subdomains.copy()
	for vol in self.volume_subdomains:
	# check reactions
	for reaction in self.reactions:
	if vol == reaction.volume:
	not_defined_in_volume.remove(vol)
	# add c = 0 to formulation where needed
	for vol in not_defined_in_volume:
	self.formulation += (
	spe.solution * spe.test_function * self.dx(vol.id)
	)

However, this will apply $c = 0$ in all the subdomains where there is a reaction, even if this reaction doesn't involve this specific species.

These lines are the issue:

FESTIM/src/festim/hydrogen_transport_problem.py

Lines 689 to 693 in 013937a

	for vol in self.volume_subdomains:
	# check reactions
	for reaction in self.reactions:
	if vol == reaction.volume:
	not_defined_in_volume.remove(vol)

It should check if the species is in the product of the reaction:

 for vol in self.volume_subdomains: 
     # check reactions 
     for reaction in self.reactions: 
         if spe in reaction.product:  # <-- check if the species is in the reaction
             if vol == reaction.volume: 
                 not_defined_in_volume.remove(vol) 

To reproduce the bug, run this on fenicsx:

import festim as F
import numpy as np

my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh1D(
    np.concatenate([np.linspace(0, 0.5, 10), np.linspace(0.5, 1, 10)])
)

mat = F.Material(D_0=1, E_D=0.1)

vol1 = F.VolumeSubdomain1D(id=1, material=mat, borders=[0, 0.5])
vol2 = F.VolumeSubdomain1D(id=2, material=mat, borders=[0.5, 1])

surface1 = F.SurfaceSubdomain1D(id=4, x=0)
surface2 = F.SurfaceSubdomain1D(id=5, x=1)

my_model.subdomains = [vol1, vol2, surface1, surface2]

mobile = F.Species(name="H", mobile=True)
trapped = F.Species(name="H_trapped", mobile=False)
trapped2 = F.Species(name="H_trapped", mobile=False)
empty = F.ImplicitSpecies(n=0.5, others=[trapped])
empty2 = F.ImplicitSpecies(n=0.5, others=[trapped2])

my_model.species = [mobile, trapped, trapped2]

my_model.reactions = [
    F.Reaction(
        reactant=[mobile, empty],
        product=[trapped],
        k_0=1,
        E_k=mat.E_D,
        p_0=0.1,
        E_p=0.87,
        volume=vol1,
    ),
    F.Reaction(
        reactant=[mobile, empty2],
        product=[trapped2],
        k_0=1,
        E_k=mat.E_D,
        p_0=0.1,
        E_p=0.87,
        volume=vol2,
    ),
]

my_model.temperature = 600

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=surface1, species=mobile, value=0),
    F.FixedConcentrationBC(subdomain=surface2, species=mobile, value=1),
]

my_model.settings = F.Settings(
    atol=1e-10,
    rtol=1e-10,
    transient=False,
)


my_model.initialise()
my_model.run()

Produces:

Traceback (most recent call last):
  File "/home/remidm/FESTIM/mwe_trick.py", line 76, in <module>
    my_model.run()
  File "/home/remidm/FESTIM/src/festim/problem.py", line 150, in run
    self.solver.solve(self.u)
  File "/home/remidm/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/nls/petsc.py", line 51, in solve
    n, converged = super().solve(u.x.petsc_vec)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 76, Error in external library