Skip to content

Subdomain temperature fields #1060

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

Open
jhdark wants to merge 5 commits into
base: fenicsx
Choose a base branch
from
Open

Subdomain temperature fields #1060

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
50b265c
move decoraters
jhdark Jan 29, 2026
aa2fe1d
working fix
jhdark Jan 29, 2026
d0a587d
removed typoed attribute
jhdark Jan 30, 2026
94b2715
use dolfinx blocked solver instead of scifem
jhdark Jan 30, 2026
0a0422d
fix iterate to work with 0.10 solver
jhdark Jan 30, 2026
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
50 changes: 11 additions & 39 deletions src/festim/boundary_conditions/dirichlet_bc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,47 +243,19 @@ def create_value(
if "T" in arguments:
kwargs["T"] = temperature

try:
self.value_fenics = fem.Function(function_space)

# store the expression of the boundary condition
# to update the value_fenics later
self.bc_expr = fem.Expression(
self.value(**kwargs),
helpers.get_interpolation_points(function_space.element),
)
self.value_fenics.interpolate(self.bc_expr)
except RuntimeError:
# if this fails, it is probably because the temperature is a Function
# from the parent mesh and this is used in a mixed domain problem.
# In this case, we need to interpolate the temperature on the submesh

submesh = mesh
mesh_cell_map = submesh.topology.index_map(submesh.topology.dim)
num_cells_on_proc = (
mesh_cell_map.size_local + mesh_cell_map.num_ghosts
)
cells = np.arange(num_cells_on_proc, dtype=np.int32)
parent_functionspace = temperature.function_space
interpolation_data = fem.create_interpolation_data(
function_space, parent_functionspace, cells
)

temperature_sub = fem.Function(function_space)
temperature_sub.interpolate_nonmatching(
temperature,
cells=cells,
interpolation_data=interpolation_data,
)
self.value_fenics = fem.Function(function_space)

# override the kwargs with the temperature_sub
kwargs["T"] = temperature_sub
# store the expression of the boundary condition
# to update the value_fenics later
assert isinstance(self.value(**kwargs), ufl.core.expr.Expr), (
f"{type(self.value(**kwargs))}"
)
self.bc_expr = fem.Expression(
self.value(**kwargs),
helpers.get_interpolation_points(function_space.element),
)

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

# if K_S is provided, divide the value by K_S (change of variable method)
if K_S is not None:
Expand Down
157 changes: 124 additions & 33 deletions src/festim/hydrogen_transport_problem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
from enum import Enum

from mpi4py import MPI
from petsc4py import PETSc

import adios4dolfinx
import basix
Expand All @@ -12,6 +13,8 @@
import tqdm.autonotebook
import ufl
from dolfinx import fem
from dolfinx.nls.petsc import NewtonSolver
from packaging.version import Version
from scifem import BlockedNewtonSolver

import festim.boundary_conditions
Expand Down Expand Up @@ -206,7 +209,6 @@ def __init__(
self.temperature_fenics = None

self._element_for_traps = "DG"
self.petcs_options = petsc_options

self._temperature_as_function = None

Expand Down Expand Up @@ -1090,6 +1092,19 @@ def __init__(
self.surface_to_volume = surface_to_volume or {}
self.subdomain_to_species = {} # maps subdomain to species defined in it

@property
def method_interface(self):
return self._method_interface

@method_interface.setter
def method_interface(self, value):
if isinstance(value, InterfaceMethod):
self._method_interface = value
elif isinstance(value, str):
self._method_interface = InterfaceMethod.from_string(value)
else:
raise TypeError("method_interface must be of type str or InterfaceMethod")

def initialise(self):
# check that all species have a list of F.VolumeSubdomain as this is
# different from F.HydrogenTransportProblem
Expand Down Expand Up @@ -1159,18 +1174,25 @@ def initialise(self):
self.create_solver()
self.initialise_exports()

@property
def method_interface(self):
return self._method_interface
def define_temperature(self):
super().define_temperature()

@method_interface.setter
def method_interface(self, value):
if isinstance(value, InterfaceMethod):
self._method_interface = value
elif isinstance(value, str):
self._method_interface = InterfaceMethod.from_string(value)
else:
raise TypeError("method_interface must be of type str or InterfaceMethod")
# pass temperature function to each subdomain
if isinstance(self.temperature_fenics, fem.Function):
for subdomain in self.volume_subdomains:
element_CG = basix.ufl.element(
basix.ElementFamily.P,
subdomain.submesh.basix_cell(),
1, # could expose?
basix.LagrangeVariant.equispaced,
)
V = dolfinx.fem.functionspace(subdomain.submesh, element_CG)
sub_T = dolfinx.fem.Function(V)
from festim.helpers import nmm_interpolate

nmm_interpolate(f_out=sub_T, f_in=self.temperature_fenics)

subdomain.sub_T = sub_T

def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC):
"""
Expand All @@ -1190,8 +1212,16 @@ def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC):
sub_V = bc.species.subdomain_to_function_space[volume_subdomain]
collapsed_V, _ = sub_V.collapse()

# in the discontinuous case, if the temperature is given as a function
# then we can't use the temperature on the parent mesh
# see issue #1007
if isinstance(self.temperature_fenics, fem.Function):
temp = volume_subdomain.sub_T
else:
temp = self.temperature_fenics

bc.create_value(
temperature=self.temperature_fenics,
temperature=temp,
function_space=collapsed_V,
t=self.t,
)
Expand Down Expand Up @@ -1609,25 +1639,78 @@ def mixed_term(u, v, n):
)

def create_solver(self):
self.solver = BlockedNewtonSolver(
self.forms,
[subdomain.u for subdomain in self.volume_subdomains],
J=self.J,
bcs=self.bc_forms,
petsc_options=self.petsc_options,
)
self.solver.max_iterations = self.settings.max_iterations
self.solver.convergence_criterion = self.settings.convergence_criterion
self.solver.atol = (
self.settings.atol
if not callable(self.settings.atol)
else self.settings.atol(float(self.t))
)
self.solver.rtol = (
self.settings.rtol
if not callable(self.settings.rtol)
else self.settings.rtol(float(self.t))
)
if Version(dolfinx.__version__) == Version("0.9.0"):
self.solver = BlockedNewtonSolver(
self.forms,
[subdomain.u for subdomain in self.volume_subdomains],
J=self.J,
bcs=self.bc_forms,
petsc_options=self.petsc_options,
)
self.solver.max_iterations = self.settings.max_iterations
self.solver.convergence_criterion = self.settings.convergence_criterion
self.solver.atol = (
self.settings.atol
if not callable(self.settings.atol)
else self.settings.atol(float(self.t))
)
self.solver.rtol = (
self.settings.rtol
if not callable(self.settings.rtol)
else self.settings.rtol(float(self.t))
)

elif Version(dolfinx.__version__) > Version("0.9.0"):
from dolfinx.fem.petsc import NonlinearProblem

if self.petsc_options is None:
# taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28
use_superlu = (
PETSc.IntType == np.int64
) # or PETSc.ScalarType == np.complex64
sys = PETSc.Sys() # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
linear_solver = "superlu_dist"
else:
linear_solver = "petsc"

petsc_options = {
"snes_type": "newtonls",
"snes_linesearch_type": "none",
"snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps)
* 1e-2,
# TODO : make atol and rtol callable
"snes_atol": self.settings.atol,
"snes_rtol": self.settings.rtol,
"snes_max_it": self.settings.max_iterations,
"snes_divergence_tolerance": "PETSC_UNLIMITED",
"ksp_type": "preonly",
"pc_type": "lu",
"snes_monitor": None,
"ksp_monitor": None,
"pc_factor_mat_solver_type": linear_solver,
}
else:
petsc_options = self.petsc_options

self.solver = NonlinearProblem(
self.forms,
[subdomain.u for subdomain in self.volume_subdomains],
bcs=self.bc_forms,
J=self.J,
petsc_options=petsc_options,
petsc_options_prefix="festim_solver",
)

# Delete PETSc options post setting them, ref:
# https://gitlab.com/petsc/petsc/-/issues/1201
snes = self.solver.solver
prefix = snes.getOptionsPrefix()
opts = PETSc.Options()
for k in petsc_options.keys():
del opts[f"{prefix}{k}"]

def create_flux_values_fenics(self):
"""For each particle flux create the ``value_fenics`` attribute"""
Expand Down Expand Up @@ -1803,7 +1886,15 @@ def iterate(self):
self.update_time_dependent_values()

# Solve main problem
nb_its, converged = self.solver.solve()
if Version(dolfinx.__version__) == Version("0.9.0"):
nb_its, converged = self.solver.solve(self.u)
elif Version(dolfinx.__version__) > Version("0.9.0"):
_ = self.solver.solve()
converged_reason = self.solver.solver.getConvergedReason()
assert converged_reason > 0, (
f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information."
)
nb_its = self.solver.solver.getIterationNumber()

# post processing
self.post_processing()
Expand Down
18 changes: 16 additions & 2 deletions src/festim/subdomain/volume_subdomain.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

import dolfinx
import numpy as np
from dolfinx import fem
from dolfinx.mesh import Mesh, locate_entities
from numpy import typing as npt
from scifem.mesh import transfer_meshtags_to_submesh
Expand All @@ -23,8 +24,20 @@ class VolumeSubdomain:
Volume subdomain class

Args:
id (int): the id of the volume subdomain
material (festim.Material): the material assigned to the subdomain
id: the id of the volume subdomain
submesh: the submesh of the volume subdomain
cell_map: the cell map of the volume subdomain
parent_mesh: the parent mesh of the volume subdomain
parent_to_submesh: the parent to submesh map of the volume subdomain
v_map: the vertex map of the volume subdomain
n_map: the normal map of the volume subdomain
facet_to_parent: the facet to parent map of the volume subdomain
ft: the facet meshtags of the volume subdomain
padded: whether the subdomain is padded (for 0.9 compatibility)
u: the solution function of the subdomain
u_n: the previous solution function of the subdomain
material: the material assigned to the subdomain
sub_T: the sub temperature field in the subdomain
"""

id: int
Expand All @@ -40,6 +53,7 @@ class VolumeSubdomain:
u: dolfinx.fem.Function
u_n: dolfinx.fem.Function
material: Material
sub_T: fem.Function | float

def __init__(self, id, material, locator: Callable | None = None):
self.id = id
Expand Down
Loading