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..66aa741a3 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,212 @@ 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 zero when the condition holds exactly; the lower + the value, the better. + + The residual is ``right - left`` where 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``. + + Args: + t: Current simulation time. + """ + 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( + f0 - f1, + get_interpolation_points(self.function.function_space.element), + ) + + class ReactionRateExport(CustomFieldExport): """Export a reaction rate to a VTX file @@ -370,7 +579,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)