Merge pull request #918 from festim-dev/penalty

Penalty method for interfaces

src/festim/__init__.py

@@ -40,6 +40,7 @@
 from .hydrogen_transport_problem import (
     HTransportProblemDiscontinuous,
     HydrogenTransportProblem,
+    HTransportProblemPenalty,
 )
 from .initial_condition import InitialCondition, InitialTemperature
 from .material import Material

src/festim/hydrogen_transport_problem.py

@@ -1203,7 +1203,7 @@ def create_solver(self):
             self.J,
             [subdomain.u for subdomain in self.volume_subdomains],
             bcs=self.bc_forms,
-            max_iterations=10,
+            max_iterations=self.settings.max_iterations,
             petsc_options=self.petsc_options,
         )
 
@@ -1299,3 +1299,107 @@ def run(self):
     def __del__(self):
         for vtxfile in self._vtxfiles:
             vtxfile.close()
+
+
+class HTransportProblemPenalty(HTransportProblemDiscontinuous):
+    def create_formulation(self):
+        """
+        Takes all the formulations for each subdomain and adds the interface conditions.
+
+        Finally compute the jacobian matrix and store it in the ``J`` attribute,
+        adds the ``entity_maps`` to the forms and store them in the ``forms`` attribute
+        """
+        mesh = self.mesh.mesh
+        mt = self.facet_meshtags
+
+        for interface in self.interfaces:
+            interface.mesh = mesh
+            interface.mt = mt
+
+        integral_data = [
+            interface.compute_mapped_interior_facet_data(mesh)
+            for interface in self.interfaces
+        ]
+        [interface.pad_parent_maps() for interface in self.interfaces]
+        dInterface = ufl.Measure("dS", domain=mesh, subdomain_data=integral_data)
+
+        entity_maps = {
+            sd.submesh: sd.parent_to_submesh for sd in self.volume_subdomains
+        }
+        for interface in self.interfaces:
+            subdomain_0, subdomain_1 = interface.subdomains
+            res = interface.restriction
+
+            all_mobile_species = [spe for spe in self.species if spe.mobile]
+            if len(all_mobile_species) > 1:
+                raise NotImplementedError("Multiple mobile species not implemented")
+            H = all_mobile_species[0]
+            v_b = H.subdomain_to_test_function[subdomain_0](res[0])
+            v_t = H.subdomain_to_test_function[subdomain_1](res[1])
+
+            u_b = H.subdomain_to_solution[subdomain_0](res[0])
+            u_t = H.subdomain_to_solution[subdomain_1](res[1])
+
+            K_b = subdomain_0.material.get_solubility_coefficient(
+                self.mesh.mesh, self.temperature_fenics(res[0]), H
+            )
+            K_t = subdomain_1.material.get_solubility_coefficient(
+                self.mesh.mesh, self.temperature_fenics(res[1]), H
+            )
+
+            if (
+                subdomain_0.material.solubility_law
+                == subdomain_1.material.solubility_law
+            ):
+                left = u_b / K_b
+                right = u_t / K_t
+            else:
+                if subdomain_0.material.solubility_law == "henry":
+                    left = u_b / K_b
+                elif subdomain_0.material.solubility_law == "sievert":
+                    left = (u_b / K_b) ** 2
+                else:
+                    raise ValueError(
+                        f"Unknown material law {subdomain_0.material.solubility_law}"
+                    )
+
+                if subdomain_1.material.solubility_law == "henry":
+                    right = u_t / K_t
+                elif subdomain_1.material.solubility_law == "sievert":
+                    right = (u_t / K_t) ** 2
+                else:
+                    raise ValueError(
+                        f"Unknown material law {subdomain_1.material.solubility_law}"
+                    )
+
+            equality = right - left
+
+            F_0 = (
+                interface.penalty_term
+                * ufl.inner(equality, v_b)
+                * dInterface(interface.id)
+            )
+            F_1 = (
+                -interface.penalty_term
+                * ufl.inner(equality, v_t)
+                * dInterface(interface.id)
+            )
+
+            subdomain_0.F += F_0
+            subdomain_1.F += F_1
+
+        J = []
+        # this is the symbolic differentiation of the Jacobian
+        for subdomain1 in self.volume_subdomains:
+            jac = []
+            for subdomain2 in self.volume_subdomains:
+                jac.append(
+                    ufl.derivative(subdomain1.F, subdomain2.u),
+                )
+            J.append(jac)
+        # compile jacobian (J) and residual (F)
+        self.forms = dolfinx.fem.form(
+            [subdomain.F for subdomain in self.volume_subdomains],
+            entity_maps=entity_maps,
+        )
+        self.J = dolfinx.fem.form(J, entity_maps=entity_maps)

src/festim/material.py

@@ -20,6 +20,7 @@ class Material:
         thermal_conductivity (float, callable): the thermal conductivity of the material (W/m/K)
         density (float, callable): the density of the material (kg/m3)
         heat_capacity (float, callable): the heat capacity of the material (J/kg/K)
+        solubility_law (str): the solubility law of the material ("sievert" or "henry")
 
     Attributes:
         D_0 (float or dict): the pre-exponential factor of the
@@ -55,6 +56,7 @@ def __init__(
         density=None,
         heat_capacity=None,
         name=None,
+        solubility_law=None,
     ) -> None:
         self.D_0 = D_0
         self.E_D = E_D
@@ -65,6 +67,7 @@ def __init__(
         self.density = density
         self.heat_capacity = heat_capacity
         self.name = name
+        self.solubility_law = solubility_law
 
     def get_D_0(self, species=None):
         """Returns the pre-exponential factor of the diffusion coefficient

test/system_tests/test_multi_mat_penalty.py

@@ -0,0 +1,161 @@
+from mpi4py import MPI
+
+import dolfinx
+import dolfinx.fem.petsc
+import numpy as np
+import ufl
+
+import festim as F
+
+from .tools import error_L2
+
+
+def generate_mesh(n=20):
+    def bottom_boundary(x):
+        return np.isclose(x[1], 0.0)
+
+    def top_boundary(x):
+        return np.isclose(x[1], 1.0)
+
+    def half(x):
+        return x[1] <= 0.5 + 1e-14
+
+    mesh = dolfinx.mesh.create_unit_square(
+        MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.triangle
+    )
+
+    # Split domain in half and set an interface tag of 5
+    gdim = mesh.geometry.dim
+    tdim = mesh.topology.dim
+    fdim = tdim - 1
+    top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
+    bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary)
+    num_facets_local = (
+        mesh.topology.index_map(fdim).size_local
+        + mesh.topology.index_map(fdim).num_ghosts
+    )
+    facets = np.arange(num_facets_local, dtype=np.int32)
+    values = np.full_like(facets, 0, dtype=np.int32)
+    values[top_facets] = 1
+    values[bottom_facets] = 2
+
+    bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
+    num_cells_local = (
+        mesh.topology.index_map(tdim).size_local
+        + mesh.topology.index_map(tdim).num_ghosts
+    )
+    cells = np.full(num_cells_local, 4, dtype=np.int32)
+    cells[bottom_cells] = 3
+    ct = dolfinx.mesh.meshtags(
+        mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
+    )
+    all_b_facets = dolfinx.mesh.compute_incident_entities(
+        mesh.topology, ct.find(3), tdim, fdim
+    )
+    all_t_facets = dolfinx.mesh.compute_incident_entities(
+        mesh.topology, ct.find(4), tdim, fdim
+    )
+    interface = np.intersect1d(all_b_facets, all_t_facets)
+    values[interface] = 5
+
+    mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)
+    return mesh, mt, ct
+
+
+def test_2_materials_2d_mms():
+    """
+    MMS case for a 2D problem with 2 materials
+    adapted from https://festim-vv-report.readthedocs.io/en/v1.0/verification/mms/discontinuity.html
+    """
+    K_S_top = 3.0
+    K_S_bot = 6.0
+    D_top = 2.0
+    D_bot = 5.0
+    c_exact_top_ufl = (
+        lambda x: 3
+        + ufl.sin(ufl.pi * (2 * x[1] + 0.5))
+        + 0.1 * ufl.cos(2 * ufl.pi * x[0])
+    )
+    c_exact_top_np = (
+        lambda x: 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0])
+    )
+
+    def c_exact_bot_ufl(x):
+        return K_S_bot / K_S_top**2 * c_exact_top_ufl(x) ** 2
+
+    def c_exact_bot_np(x):
+        return K_S_bot / K_S_top**2 * c_exact_top_np(x) ** 2
+
+    mesh, mt, ct = generate_mesh(100)
+
+    my_model = F.HTransportProblemPenalty()
+    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,
+    ]
+
+    my_model.interfaces = [
+        F.Interface(5, (bottom_domain, top_domain), penalty_term=1),
+    ]
+    my_model.surface_to_volume = {
+        top_surface: top_domain,
+        bottom_surface: bottom_domain,
+    }
+
+    H = F.Species("H", mobile=True)
+
+    my_model.species = [H]
+
+    for species in my_model.species:
+        species.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)
+
+    source_top_val = -ufl.div(D_top * ufl.grad(c_exact_top_ufl(x)))
+    source_bottom_val = -ufl.div(D_bot * ufl.grad(c_exact_bot_ufl(x)))
+    my_model.sources = [
+        F.ParticleSource(volume=top_domain, species=H, value=source_top_val),
+        F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val),
+    ]
+
+    my_model.temperature = 500.0  # lambda x: 300 + 10 * x[1] + 100 * x[0]
+
+    my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+    my_model.exports = [
+        F.VTXSpeciesExport(f"u_{subdomain.id}.bp", field=H, subdomain=subdomain)
+        for subdomain in my_model.volume_subdomains
+    ]
+
+    my_model.initialise()
+    my_model.run()
+
+    c_top_computed = H.subdomain_to_post_processing_solution[top_domain]
+    c_bot_computed = H.subdomain_to_post_processing_solution[bottom_domain]
+
+    L2_error_top = error_L2(c_top_computed, c_exact_top_np)
+    L2_error_bot = error_L2(c_bot_computed, c_exact_bot_np)
+
+    assert L2_error_top < 1e-3
+    assert L2_error_bot < 1e-3