Skip to content

Only used mixed function spaces #1042

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jhdark merged 9 commits into from
Nov 12, 2025
Merged

Only used mixed function spaces #1042

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
048917f
only use mixed space
jhdark Oct 22, 2025
07ce65a
fix initial condition tests
jhdark Oct 22, 2025
df9691a
fi species test
jhdark Oct 22, 2025
40eeffa
update test_h_transport_problem
jhdark Oct 22, 2025
cb70bf4
fix dc bc tests
jhdark Oct 22, 2025
9d77566
fixes for cov method
jhdark Oct 23, 2025
b07c027
Merge branch 'fenicsx' into test_only_mixed_space
jhdark Oct 27, 2025
a9579da
keep using scifem blockedsolver
jhdark Oct 27, 2025
226fbab
test case with CG traps
jhdark Oct 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
167 changes: 54 additions & 113 deletions src/festim/hydrogen_transport_problem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,6 @@ class HydrogenTransportProblem(problem.ProblemBase):
model
formulation (ufl.form.Form): the formulation of the model
solver (dolfinx.nls.newton.NewtonSolver): the solver of the model
multispecies (bool): True if the model has more than one species.
temperature_fenics (fem.Constant or fem.Function): the
temperature of the model as a fenics object (fem.Constant or
fem.Function).
Expand Down Expand Up @@ -256,10 +255,6 @@ def temperature_time_dependent(self):
else:
return False

@property
def multispecies(self):
return len(self.species) > 1

@property
def species(self) -> list[_species.Species]:
return self._species
Expand Down Expand Up @@ -569,15 +564,13 @@ def define_D_global(self, species):
D.interpolate(D_expr)
return D, D_expr

def define_function_spaces(self, element_degree=1):
"""Creates the function space of the model, creates a mixed element if
model is multispecies. Creates the main solution and previous solution
function u and u_n. Create global DG function spaces of degree 0 and 1
for the global diffusion coefficient.
def define_function_spaces(self, element_degree: int = 1):
"""Creates the function space of the modelw with a mixed element. Creates the
main solution and previous solution function u and u_n. Create global DG
function spaces of degree 0 and 1 for the global diffusion coefficient.

Args:
element_degree (int, optional): Degree order for finite element.
Defaults to 1.
element_degree: Degree order for finite element. Defaults to 1.
"""

element_CG = basix.ufl.element(
Expand All @@ -593,19 +586,16 @@ def define_function_spaces(self, element_degree=1):
basix.LagrangeVariant.equispaced,
)

if not self.multispecies:
element = element_CG
else:
elements = []
for spe in self.species:
if isinstance(spe, _species.Species):
if spe.mobile:
elements.append(element_CG)
elif self._element_for_traps == "DG":
elements.append(element_DG)
else:
elements.append(element_CG)
element = basix.ufl.mixed_element(elements)
elements = []
for spe in self.species:
if isinstance(spe, _species.Species):
if spe.mobile:
elements.append(element_CG)
elif self._element_for_traps == "DG":
elements.append(element_DG)
else:
elements.append(element_CG)
element = basix.ufl.mixed_element(elements)
Copy link
Collaborator

@jorgensd jorgensd Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will there be strong coupling across all these species?

If not, one should probably have a "global" MixedFunctionSpace where each sub-space is its own function space, and this class would keep track of the local->global indexing.

Copy link
Collaborator Author

@jhdark jhdark Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There will be some strong coupling between different combinations of species from the reaction framework.

Copy link
Collaborator

@jorgensd jorgensd Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That there is coupling between them is on for either approach, I’ll try to rephrase:

Given M species, if we make an MxM matrix where M_ij =1 means that there is a coupling term, M_ij=0 means that there is not, how many non-zeros do we usually get?


self.function_space = fem.functionspace(self.mesh.mesh, element)

Expand All @@ -629,31 +619,21 @@ def define_function_spaces(self, element_degree=1):
self.u_n = fem.Function(self.function_space)

def assign_functions_to_species(self):
"""Creates the solution, prev solution, test function and
post-processing solution for each species, if model is multispecies,
created a collapsed function space for each species"""

if not self.multispecies:
sub_solutions = [self.u]
sub_prev_solution = [self.u_n]
sub_test_functions = [ufl.TestFunction(self.function_space)]
self.species[0].sub_function_space = self.function_space
self.species[0].post_processing_solution = self.u
self.species[0].sub_function = self.u
else:
sub_solutions = list(ufl.split(self.u))
sub_prev_solution = list(ufl.split(self.u_n))
sub_test_functions = list(ufl.TestFunctions(self.function_space))

for idx, spe in enumerate(self.species):
spe.sub_function_space = self.function_space.sub(idx)
spe.sub_function = self.u.sub(
idx
) # TODO add this to discontinuous class
spe.post_processing_solution = self.u.sub(idx).collapse()
spe.collapsed_function_space, spe.map_sub_to_main_solution = (
self.function_space.sub(idx).collapse()
)
"""Creates the solution, prev solution, test function and post-processing
solution for each species, as well as a collapsed function space for each
species"""

sub_solutions = list(ufl.split(self.u))
sub_prev_solution = list(ufl.split(self.u_n))
sub_test_functions = list(ufl.TestFunctions(self.function_space))

for idx, spe in enumerate(self.species):
spe.sub_function_space = self.function_space.sub(idx)
spe.sub_function = self.u.sub(idx) # TODO add this to discontinuous class
spe.post_processing_solution = self.u.sub(idx).collapse()
spe.collapsed_function_space, spe.map_sub_to_main_solution = (
self.function_space.sub(idx).collapse()
)
Comment on lines +631 to +636
Copy link
Collaborator

@jorgensd jorgensd Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not just relevant for this PR, but currently, by assigning data in this way, without having a getter/setter structure, users can overwrite spe.sub_function_space in their own code, and invalidate the logic in FESTIM.
In general, I think it is a good idea to do something like:

from mpi4py import MPI
import dolfinx


class Species:
    _sub_function_space: None | dolfinx.fem.function.FunctionSpace

    def __init__(self):
        self._sub_function_space = None

    @property
    def sub_function_space(self) -> dolfinx.fem.FunctionSpace:
        if self._sub_function_space is None:
            raise ValueError("sub_function_space has not been set yet")
        else:
            return self._sub_function_space

    def set_sub_function_space(self, function_space: dolfinx.fem.FunctionSpace) -> None:
        if self._sub_function_space is not None:
            raise ValueError("sub_function_space has already been set")
        self._sub_function_space = function_space

to ensure that data is not set twice (mistakenly).
If you want these to be adaptive and be possible to set multiple times, I would advice having some warnings on it.

Copy link
Collaborator Author

@jhdark jhdark Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, this would be better. I think we often overwrite this or assign it ourselves in a lot of the tests to simplify them, though, so it may make them more complicated


for idx, spe in enumerate(self.species):
spe.solution = sub_solutions[idx]
Expand Down Expand Up @@ -695,19 +675,14 @@ def create_dirichletbc_form(self, bc):
the boundary condition for modifying linear systems.
"""
# create value_fenics
if not self.multispecies:
function_space_value = bc.species.sub_function_space
else:
function_space_value = bc.species.collapsed_function_space

bc.create_value(
temperature=self.temperature_fenics,
function_space=function_space_value,
function_space=bc.species.collapsed_function_space,
t=self.t,
)

# get dofs
if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
if isinstance(bc.value_fenics, (fem.Function)):
function_space_dofs = (
bc.species.sub_function_space,
bc.species.collapsed_function_space,
Expand All @@ -721,16 +696,10 @@ def create_dirichletbc_form(self, bc):
)

# create form
if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
# no need to pass the functionspace since value_fenics is already a Function
function_space_form = None
else:
function_space_form = bc.species.sub_function_space

form = fem.dirichletbc(
value=bc.value_fenics,
dofs=bc_dofs,
V=function_space_form,
V=bc.species.sub_function_space,
)

return form
Expand Down Expand Up @@ -779,10 +748,7 @@ def create_initial_conditions(self):
function_space_value = None
if callable(condition.value):
# if bc.value is a callable then need to provide a functionspace
if not self.multispecies:
function_space_value = condition.species.sub_function_space
else:
function_space_value = condition.species.collapsed_function_space
function_space_value = condition.species.collapsed_function_space

condition.create_expr_fenics(
mesh=self.mesh.mesh,
Expand All @@ -792,11 +758,8 @@ def create_initial_conditions(self):

# assign to previous solution of species
entities = self.volume_meshtags.find(condition.volume.id)
if not self.multispecies:
function_to_interpolate_on = condition.species.prev_solution
else:
idx = self.species.index(condition.species)
function_to_interpolate_on = self.u_n.sub(idx)
idx = self.species.index(condition.species)
function_to_interpolate_on = self.u_n.sub(idx)

function_to_interpolate_on.interpolate(
condition.expr_fenics, cells0=entities
Expand Down Expand Up @@ -961,12 +924,10 @@ def update_time_dependent_values(self):
def update_post_processing_solutions(self):
"""Updates the post-processing solutions of each species"""

# update post-processing for mixed function space
if self.multispecies:
for spe in self.species:
spe.post_processing_solution.x.array[:] = self.u.x.array[
spe.map_sub_to_main_solution
]
for spe in self.species:
spe.post_processing_solution.x.array[:] = self.u.x.array[
spe.map_sub_to_main_solution
]

def post_processing(self):
"""Post processes the model"""
Expand Down Expand Up @@ -1054,26 +1015,16 @@ def post_processing(self):

# computing dofs at each time step is costly so storing it in the export
if export._dofs is None:
# if multispecies then we have a mixed function space
# (mixed element)
if self.multispecies:
index = self.species.index(export.field)
V0, export._dofs = self.u.function_space.sub(index).collapse()
coords = V0.tabulate_dof_coordinates()[:, 0]
export._sort_coords = np.argsort(coords)
x = coords[export._sort_coords]
# if not multispecies then we have a single function space
else:
coords = self.u.function_space.tabulate_dof_coordinates()[:, 0]
export._sort_coords = np.argsort(coords)
x = coords[export._sort_coords]
index = self.species.index(export.field)
V0, export._dofs = self.u.function_space.sub(index).collapse()
coords = V0.tabulate_dof_coordinates()[:, 0]
export._sort_coords = np.argsort(coords)
x = coords[export._sort_coords]

export.x = x

if self.multispecies:
c = self.u.x.array[export._dofs][export._sort_coords]
else:
c = self.u.x.array[export._dofs]
c = self.u.x.array[export._dofs][export._sort_coords]

export.data.append(c)
export.t.append(float(self.t))

Expand Down Expand Up @@ -1666,8 +1617,8 @@ def create_solver(self):
self.solver.convergence_criterion = self.settings.convergence_criterion
self.solver.atol = (
self.settings.atol
if not callable(self.settings.rtol)
else self.settings.rtol(float(self.t))
if not callable(self.settings.atol)
else self.settings.atol(float(self.t))
)
self.solver.rtol = (
self.settings.rtol
Expand Down Expand Up @@ -2070,6 +2021,7 @@ def update_post_processing_solutions(self):
"""Updates the post-processing solutions after each time step"""
# need to compute c = theta * K_S
# this expression is stored in species.dg_expr

for spe in self.species:
if not spe.mobile:
continue
Expand All @@ -2085,12 +2037,6 @@ def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
dolfinx.fem.bcs.DirichletBC: A representation of
the boundary condition for modifying linear systems.
"""
# create value_fenics
if not self.multispecies:
function_space_value = bc.species.sub_function_space
else:
function_space_value = bc.species.collapsed_function_space

# create K_S function
Q0 = fem.functionspace(self.mesh.mesh, ("DG", 0))
K_S0 = fem.Function(Q0)
Expand All @@ -2102,15 +2048,16 @@ def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):

K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))

# create value_fenics
bc.create_value(
temperature=self.temperature_fenics,
function_space=function_space_value,
function_space=bc.species.collapsed_function_space,
t=self.t,
K_S=K_S,
)

# get dofs
if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
if isinstance(bc.value_fenics, fem.Function):
function_space_dofs = (
bc.species.sub_function_space,
bc.species.collapsed_function_space,
Expand All @@ -2124,16 +2071,10 @@ def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
)

# create form
if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
# no need to pass the functionspace since value_fenics is already a Function
function_space_form = None
else:
function_space_form = bc.species.sub_function_space

form = fem.dirichletbc(
value=bc.value_fenics,
dofs=bc_dofs,
V=function_space_form,
V=bc.species.sub_function_space,
)

return form
Expand Down
9 changes: 6 additions & 3 deletions test/test_dirichlet_bc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,11 @@ def value(x, t):

my_model.define_function_spaces()
my_model.define_meshtags_and_measures()
my_model.assign_functions_to_species()

T = fem.Constant(my_model.mesh.mesh, 550.0)
t = fem.Constant(my_model.mesh.mesh, 0.0)
bc.create_value(my_model.function_space, T, t)
bc.create_value(species.collapsed_function_space, T, t)

# check that the value_fenics attribute is set correctly
assert isinstance(bc.value_fenics, fem.Function)
Expand Down Expand Up @@ -116,10 +117,11 @@ def value(x, t, T):

my_model.define_function_spaces()
my_model.define_meshtags_and_measures()
my_model.assign_functions_to_species()

T = fem.Constant(my_model.mesh.mesh, 550.0)
t = fem.Constant(my_model.mesh.mesh, 0.0)
bc.create_value(my_model.function_space, T, t)
bc.create_value(species.collapsed_function_space, T, t)

# check that the value_fenics attribute is set correctly
assert isinstance(bc.value_fenics, fem.Function)
Expand Down Expand Up @@ -205,12 +207,13 @@ def value(x):

my_model.define_function_spaces()
my_model.define_meshtags_and_measures()
my_model.assign_functions_to_species()

T = fem.Constant(my_model.mesh.mesh, 550.0)
t = fem.Constant(my_model.mesh.mesh, 0.0)

# TEST
bc.create_value(my_model.function_space, T, t)
bc.create_value(species.collapsed_function_space, T, t)

# check that the value_fenics attribute is set correctly
assert isinstance(bc.value_fenics, fem.Function)
Expand Down
Loading
Loading