diff --git a/src/festim/__init__.py b/src/festim/__init__.py index b2095f861..4a7b5244a 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -42,6 +42,7 @@ CustomFieldExport, ExportBaseClass, ReactionRateExport, + VTXInterfaceResidualExport, VTXSpeciesExport, VTXTemperatureExport, ) diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index ad3988734..597ba7c09 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -16,6 +16,7 @@ CustomFieldExport, ExportBaseClass, ReactionRateExport, + VTXInterfaceResidualExport, VTXSpeciesExport, VTXTemperatureExport, ) @@ -38,6 +39,7 @@ "SurfaceQuantity", "TotalSurface", "TotalVolume", + "VTXInterfaceResidualExport", "VTXSpeciesExport", "VTXTemperatureExport", "VolumeQuantity", diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 26e67f1f0..1d5ee3db7 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -4,6 +4,8 @@ from pathlib import Path from typing import Union +import dolfinx +import numpy as np import ufl from dolfinx import fem, io @@ -11,6 +13,7 @@ from festim.helpers import get_interpolation_points from festim.reaction import Reaction from festim.species import ImplicitSpecies, Species +from festim.subdomain.interface import Interface, interface_condition_term from festim.subdomain.volume_subdomain import VolumeSubdomain @@ -345,6 +348,209 @@ def check_valid_inputs(self, kwargs: dict): ) +class VTXInterfaceResidualExport(ExportBaseClass): + """Export the interface condition residual to a VTX file. + + This quantity measures how well the penalty/Nitsche interface condition is + satisfied at the interface. It is the absolute difference between the two sides' + terms, ``|right - left|``, so it is always non-negative, zero when the condition + holds exactly; the lower the value, the better. + + Each side's term depends on the solubility laws of the two subdomains: + + - **Same law on both sides** (Henry-Henry or Sievert-Sievert): + ``residual = |c_1/K_S_1 - c_0/K_S_0|`` + - **Henry (side 0) - Sievert (side 1)**: + ``residual = |(c_1/K_S_1)^2 - c_0/K_S_0|`` + - **Sievert (side 0) - Henry (side 1)**: + ``residual = |c_1/K_S_1 - (c_0/K_S_0)^2|`` + + Args: + field: The species whose interface residual is exported. + filename: The name of the output file. + interface: The interface between the two subdomains. + times: if provided, the field will be exported at these timesteps. + Otherwise exports at all timesteps. Defaults to None. + + Attributes: + field: The species to export. + filename: The name of the output file. + interface: The interface between the two subdomains. + function: Residual function on the interface submesh. Set by + ``initialise``. + writer: VTXWriter used to write the output file. Set by ``initialise``. + """ + + field: Species + interface: Interface + times: list[float] | list[int] | None + + function: fem.Function + writer: io.VTXWriter + + def __init__( + self, + field: Species, + filename: str | Path, + interface: Interface, + times: list[float | int] | None = None, + ): + super().__init__(filename, ".bp", times) + self.field = field + self.interface = interface + + def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None: + """Create the interface submesh, interpolation data, and VTX writer. + + Called by the problem during ``initialise_exports``. Builds a CG1 + function space on the interface submesh, pre-computes + ``create_interpolation_data`` for both subdomain concentrations and + (if space-dependent) temperature, and opens the VTXWriter. + + Args: + temperature_fenics: Temperature field on the parent mesh. Either a + ``fem.Constant`` (uniform) or a ``fem.Function`` (spatially + varying). + """ + parent_mesh = self.interface.parent_mesh + fdim = parent_mesh.topology.dim - 1 + interface_facets = self.interface.mt.find(self.interface.id) + + interface_submesh, _, _, _ = dolfinx.mesh.create_submesh( + parent_mesh, fdim, interface_facets + ) + V_interface = fem.functionspace(interface_submesh, ("CG", 1)) + + self._c_0_interface = fem.Function( + V_interface, name=f"{self.field.name}_interface_c0" + ) + self._c_1_interface = fem.Function( + V_interface, name=f"{self.field.name}_interface_c1" + ) + + self._f_0_interface = fem.Function( + V_interface, name=f"{self.field.name}_interface_f0" + ) + self._f_1_interface = fem.Function( + V_interface, name=f"{self.field.name}_interface_f1" + ) + self.function = fem.Function(V_interface) + self.function.name = f"{self.field.name}_interface_residual" + + imap = interface_submesh.topology.index_map(interface_submesh.topology.dim) + self._interface_cells = np.arange( + imap.size_local + imap.num_ghosts, dtype=np.int32 + ) + + subdomain_0, subdomain_1 = self.interface.subdomains + V_0 = self.field.subdomain_to_post_processing_solution[ + subdomain_0 + ].function_space + V_1 = self.field.subdomain_to_post_processing_solution[ + subdomain_1 + ].function_space + + self._interp_data_0 = fem.create_interpolation_data( + V_interface, V_0, self._interface_cells, padding=1e-11 + ) + self._interp_data_1 = fem.create_interpolation_data( + V_interface, V_1, self._interface_cells, padding=1e-11 + ) + + self._temperature_fenics = temperature_fenics + if isinstance(temperature_fenics, fem.Constant): + self._T_func = None + else: + self._T_func = fem.Function(V_interface) + self._T_interp_data = fem.create_interpolation_data( + V_interface, + temperature_fenics.function_space, + self._interface_cells, + padding=1e-11, + ) + + self.writer = io.VTXWriter( + comm=interface_submesh.comm, + filename=self.filename, + output=[self.function, self._f_0_interface, self._f_1_interface], + engine="BP5", + ) + + def set_dolfinx_expression(self) -> None: + """Compute the interface condition residual. + + Interpolates the concentration from each subdomain onto the interface + submesh, evaluates ``K_S = K_S_0 * exp(-E_K_S / (k_B * T))`` at the + current temperature, then computes ``|right - left|`` via + :func:`interface_condition_term`: + + - **Same law or Henry (side 0)**: ``left = c_0 / K_S_0`` + - **Sievert (side 0, mixed only)**: ``left = (c_0 / K_S_0)^2`` + + and symmetrically for ``right``. + """ + subdomain_0, subdomain_1 = self.interface.subdomains + + # Get the concentration from each subdomain on the interface submesh + c_0 = self.field.subdomain_to_post_processing_solution[subdomain_0] + c_1 = self.field.subdomain_to_post_processing_solution[subdomain_1] + + # FIXME: once we support dolfinx 0.11 we should be able to mix parent and + # sub-meshes in the same expression + + # NOTE: we need to do this multiple times to update the residual + self._c_0_interface.interpolate_nonmatching( + c_0, self._interface_cells, interpolation_data=self._interp_data_0 + ) + self._c_1_interface.interpolate_nonmatching( + c_1, self._interface_cells, interpolation_data=self._interp_data_1 + ) + + # Get the temperature on the interface submesh + if self._T_func is not None: + self._T_func.interpolate_nonmatching( + self._temperature_fenics, + self._interface_cells, + interpolation_data=self._T_interp_data, + ) + T = self._T_func + else: + T = self._temperature_fenics + + # Compute f_i = (c_i / K_S_i) ^ {1, 2} on the interface submesh + K_S_0 = subdomain_0.material.get_K_S_0(self.field) * ufl.exp( + -subdomain_0.material.get_E_K_S(self.field) / (_k_B * T) + ) + f0 = interface_condition_term( + self._c_0_interface, + K_S_0, + subdomain_0.material.solubility_law, + subdomain_1.material.solubility_law, + ) + + K_S_1 = subdomain_1.material.get_K_S_0(self.field) * ufl.exp( + -subdomain_1.material.get_E_K_S(self.field) / (_k_B * T) + ) + f1 = interface_condition_term( + self._c_1_interface, + K_S_1, + subdomain_1.material.solubility_law, + subdomain_0.material.solubility_law, + ) + + self.f_0_expr = fem.Expression( + f0, get_interpolation_points(self.function.function_space.element) + ) + self.f_1_expr = fem.Expression( + f1, get_interpolation_points(self.function.function_space.element) + ) + + self.residual_expr = fem.Expression( + abs(f1 - f0), + get_interpolation_points(self.function.function_space.element), + ) + + class ReactionRateExport(CustomFieldExport): """Export a reaction rate to a VTX file @@ -370,7 +576,6 @@ def __init__( subdomain: VolumeSubdomain | None = None, checkpoint: bool = False, ): - reactant_names = [reactant.name for reactant in reaction.reactant] if isinstance(reaction.product, list): product_names = [product.name for product in reaction.product] diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 1e6682ccb..68705043e 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1738,6 +1738,9 @@ def initialise_exports(self): self.temperature_fenics, engine="BP5", ) + elif isinstance(export, exports.VTXInterfaceResidualExport): + export.initialise(self.temperature_fenics) + elif isinstance(export, exports.CustomFieldExport): # need to find an appropriate function space on the right submesh V = self.subdomain_to_V_CG1[export.subdomain] @@ -1853,6 +1856,14 @@ def post_processing(self): export.writer.write(float(self.t)) elif isinstance(export, exports.VTXTemperatureExport): export.writer.write(float(self.t)) + elif isinstance(export, exports.VTXInterfaceResidualExport): + # FIXME: don't need to remake the whole expression everytime just + # need to update "sub" functions + export.set_dolfinx_expression() + export.function.interpolate(export.residual_expr) + export._f_0_interface.interpolate(export.f_0_expr) + export._f_1_interface.interpolate(export.f_1_expr) + export.writer.write(float(self.t)) else: raise NotImplementedError( f"Export type {type(export)} not implemented" diff --git a/src/festim/subdomain/interface.py b/src/festim/subdomain/interface.py index 8b568d40e..088e56a40 100644 --- a/src/festim/subdomain/interface.py +++ b/src/festim/subdomain/interface.py @@ -1,8 +1,10 @@ +from abc import ABC, abstractmethod from enum import Enum from typing import TYPE_CHECKING import dolfinx import numpy as np +import numpy.typing as npt import ufl from dolfinx.cpp.fem import compute_integration_domains from packaging.version import Version @@ -13,7 +15,38 @@ if TYPE_CHECKING: from festim.species import Species -from abc import ABC, abstractmethod + +def interface_condition_term( + c: ufl.core.expr.Expr | npt.NDArray[np.floating], + K_S: ufl.core.expr.Expr | npt.NDArray[np.floating], + law: SolubilityLaw, + other_law: SolubilityLaw, +) -> ufl.core.expr.Expr | npt.NDArray[np.floating]: + """Return the interface condition expression for one side. + + When both sides share the same solubility law the condition simplifies to + ``c/K_S`` on each side. When the laws differ Henry is expressed as + ``c/K_S`` and Sievert as ``(c/K_S)^2``. + + Args: + c: Concentration on this side. + K_S: Solubility coefficient on this side. + law: Solubility law for this side. + other_law: Solubility law for the other side. + + Returns: + The interface condition term for this side. + + Raises: + ValueError: If ``law`` is not ``HENRY`` or ``SIEVERT`` in the mixed-law + case. + """ + if law == other_law or law == SolubilityLaw.HENRY: + return c / K_S + elif law == SolubilityLaw.SIEVERT: + return (c / K_S) ** 2 + else: + raise ValueError(f"Unsupported solubility law: {law}") class InterfaceMethod(Enum): @@ -393,32 +426,12 @@ def penalty_method(self, dS, species, temperature): u_0, u_1 = self.us(species) v_0, v_1 = self.vs(species) K_0, K_1 = self.Ks(species, temperature) - if subdomain_0.material.solubility_law == subdomain_1.material.solubility_law: - left = u_0 / K_0 - right = u_1 / K_1 - else: - match subdomain_0.material.solubility_law: - case SolubilityLaw.HENRY: - left = u_0 / K_0 - case SolubilityLaw.SIEVERT: - left = (u_0 / K_0) ** 2 - case _: - raise ValueError( - "Unsupported material law " - + f"{subdomain_0.material.solubility_law}" - ) - - match subdomain_1.material.solubility_law: - case SolubilityLaw.HENRY: - right = u_1 / K_1 - case SolubilityLaw.SIEVERT: - right = (u_1 / K_1) ** 2 - case _: - raise ValueError( - f"Unsupported material law " - f"{subdomain_1.material.solubility_law}" - ) - + law_0, law_1 = ( + subdomain_0.material.solubility_law, + subdomain_1.material.solubility_law, + ) + left = interface_condition_term(u_0, K_0, law_0, law_1) + right = interface_condition_term(u_1, K_1, law_1, law_0) equality = right - left F_0 = self.penalty_term * ufl.inner(equality, v_0) * dS(self.id) diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py index d01b4dca5..32c155270 100644 --- a/test/system_tests/test_multi_mat_penalty.py +++ b/test/system_tests/test_multi_mat_penalty.py @@ -1,8 +1,8 @@ from mpi4py import MPI import dolfinx -import dolfinx.fem.petsc import numpy as np +import pytest import ufl import festim as F @@ -283,3 +283,73 @@ def test_penalty_multispecies(): my_model.initialise() my_model.run() + + +@pytest.mark.parametrize("T_value", [500, lambda x: 300 + 10 * x[1] + 100 * x[0]]) +def test_interface_residual_export_below_tolerance(tmpdir, T_value): + """A high penalty term drives the exported interface residual below 1e-3.""" + K_S_top, K_S_bot, D_top, D_bot = 3.0, 6.0, 2.0, 5.0 + + def c_exact_top_ufl(x): + return 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0]) + + def c_exact_bot_ufl(x): + return K_S_bot / K_S_top**2 * c_exact_top_ufl(x) ** 2 + + mesh, mt, ct = generate_mesh(20) + + my_model = F.HydrogenTransportProblemDiscontinuous() + my_model.mesh = F.Mesh(mesh) + my_model.volume_meshtags = ct + my_model.facet_meshtags = mt + + material_top = F.Material(D_0=D_top, E_D=0, K_S_0=K_S_top, E_K_S=0) + material_bottom = F.Material(D_0=D_bot, E_D=0, K_S_0=K_S_bot, E_K_S=0) + material_top.solubility_law = "sievert" + material_bottom.solubility_law = "henry" + + top_domain = F.VolumeSubdomain(4, material=material_top) + bottom_domain = F.VolumeSubdomain(3, material=material_bottom) + top_surface = F.SurfaceSubdomain(id=1) + bottom_surface = F.SurfaceSubdomain(id=2) + my_model.subdomains = [bottom_domain, top_domain, top_surface, bottom_surface] + + interface = F.Interface(5, (bottom_domain, top_domain), penalty_term=10000) + my_model.interfaces = [interface] + + H = F.Species("H", mobile=True) + my_model.species = [H] + H.subdomains = [bottom_domain, top_domain] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(top_surface, value=c_exact_top_ufl, species=H), + F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H), + ] + + x = ufl.SpatialCoordinate(mesh) + my_model.sources = [ + F.ParticleSource( + volume=top_domain, + species=H, + value=-ufl.div(D_top * ufl.grad(c_exact_top_ufl(x))), + ), + F.ParticleSource( + volume=bottom_domain, + species=H, + value=-ufl.div(D_bot * ufl.grad(c_exact_bot_ufl(x))), + ), + ] + my_model.temperature = T_value + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + residual_export = F.VTXInterfaceResidualExport( + field=H, + filename=str(tmpdir.join("interface_residual.bp")), + interface=interface, + ) + my_model.exports = [residual_export] + + my_model.initialise() + my_model.run() + + assert np.max(np.abs(residual_export.function.x.array)) < 1e-3 diff --git a/test/test_interface.py b/test/test_interface.py new file mode 100644 index 000000000..5dfbf260b --- /dev/null +++ b/test/test_interface.py @@ -0,0 +1,35 @@ +import pytest + +from festim.material import SolubilityLaw +from festim.subdomain.interface import interface_condition_term + + +@pytest.mark.parametrize( + "law", + [SolubilityLaw.SIEVERT, SolubilityLaw.HENRY, SolubilityLaw.NONE], +) +def test_interface_condition_term_same_law_returns_ratio(law): + """Identical solubility laws on both sides give the linear term c / K_S.""" + assert interface_condition_term(6.0, 3.0, law, law) == 2.0 + + +def test_interface_condition_term_henry_side_of_mixed_interface_returns_ratio(): + """A Henry side facing a different law gives the linear term c / K_S.""" + result = interface_condition_term( + 6.0, 3.0, SolubilityLaw.HENRY, SolubilityLaw.SIEVERT + ) + assert result == 2.0 + + +def test_interface_condition_term_sievert_side_of_mixed_interface_returns_squared(): + """A Sievert side facing a different law gives the squared term (c / K_S) ** 2.""" + result = interface_condition_term( + 6.0, 3.0, SolubilityLaw.SIEVERT, SolubilityLaw.HENRY + ) + assert result == 4.0 + + +def test_interface_condition_term_raises_for_unsupported_mixed_law(): + """An unsupported law (NONE) facing a different law raises a ValueError.""" + with pytest.raises(ValueError, match="Unsupported solubility law"): + interface_condition_term(6.0, 3.0, SolubilityLaw.NONE, SolubilityLaw.HENRY) diff --git a/test/test_vtx.py b/test/test_vtx.py index 200810dfb..4764d3a28 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -538,3 +538,39 @@ def test_export_base_class_times_and_extension(tmp_path): def test_export_base_class_no_times(tmp_path): export = F.ExportBaseClass(filename=tmp_path / "correct.bp", ext=".bp", times=None) assert export.times is None + + +def test_vtx_interface_residual_export_stores_field_and_interface(tmpdir): + """Test that the interface residual export stores its field and interface.""" + mat = F.Material(D_0=1, E_D=0) + interface = F.Interface( + id=5, + subdomains=( + F.VolumeSubdomain(id=1, material=mat), + F.VolumeSubdomain(id=2, material=mat), + ), + ) + species = F.Species("H") + filename = str(tmpdir.join("residual.bp")) + my_export = F.VTXInterfaceResidualExport( + field=species, filename=filename, interface=interface + ) + assert my_export.field is species + assert my_export.interface is interface + + +def test_vtx_interface_residual_export_suffix_converter(tmpdir): + """Test that a non-.bp filename is corrected to a .bp suffix.""" + mat = F.Material(D_0=1, E_D=0) + interface = F.Interface( + id=5, + subdomains=( + F.VolumeSubdomain(id=1, material=mat), + F.VolumeSubdomain(id=2, material=mat), + ), + ) + filename = str(tmpdir.join("residual.txt")) + my_export = F.VTXInterfaceResidualExport( + field=F.Species("H"), filename=filename, interface=interface + ) + assert my_export.filename.suffix == ".bp"