Merge pull request #1013 from festim-dev/profile_1d

Profile 1d export

Author: RemDelaporteMathurin

src/festim/__init__.py

@@ -29,6 +29,7 @@
 from .exports.maximum_volume import MaximumVolume
 from .exports.minimum_surface import MinimumSurface
 from .exports.minimum_volume import MinimumVolume
+from .exports.profile_1d import Profile1DExport
 from .exports.surface_flux import SurfaceFlux
 from .exports.surface_quantity import SurfaceQuantity
 from .exports.total_surface import TotalSurface

src/festim/exports/__init__.py

@@ -4,6 +4,7 @@
 from .maximum_volume import MaximumVolume
 from .minimum_surface import MinimumSurface
 from .minimum_volume import MinimumVolume
+from .profile_1d import Profile1DExport
 from .surface_flux import SurfaceFlux
 from .surface_quantity import SurfaceQuantity
 from .total_surface import TotalSurface
@@ -20,6 +21,7 @@
     "MaximumVolume",
     "MinimumSurface",
     "MinimumVolume",
+    "Profile1DExport",
     "SurfaceFlux",
     "SurfaceQuantity",
     "TotalSurface",

src/festim/exports/profile_1d.py

@@ -0,0 +1,36 @@
+import numpy as np
+from festim.species import Species
+from festim.subdomain import VolumeSubdomain
+
+
+class Profile1DExport:
+    """Class to export 1D profiles of a field in a simulation.
+
+    Args:
+        field: the species for which the profile is computed.
+        subdomain: the volume subdomain to compute the profile on. If None, the profile
+            is computed over the entire domain.
+
+    Attributes:
+        field: the species for which the profile is computed.
+        subdomain: the volume subdomain to compute the profile on. If None, the profile
+            is computed over the entire domain.
+        x: the coordinates along which the profile is computed.
+        data: the computed profile data.
+    """
+
+    x: np.ndarray
+    data: list
+    field: Species
+    subdomain: VolumeSubdomain | None
+    _dofs: np.ndarray
+    _sort_coords: np.ndarray
+
+    def __init__(self, field: Species, subdomain: VolumeSubdomain = None):
+        self.field = field
+        self.data = []
+        self.x = None
+        self.subdomain = subdomain
+
+        self._dofs = None
+        self._sort_coords = None

src/festim/hydrogen_transport_problem.py

@@ -12,6 +12,7 @@
 import ufl
 from dolfinx import fem
 from scifem import BlockedNewtonSolver
+import numpy as np
 
 import festim.boundary_conditions
 import festim.problem
@@ -971,12 +972,25 @@ def post_processing(self):
                     export.write(t=float(self.t))
             if isinstance(export, exports.XDMFExport):
                 export.write(float(self.t))
+            if isinstance(export, exports.Profile1DExport):
+                if export._dofs is None:
+                    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
+
+                c = self.u.x.array[export._dofs][export._sort_coords]
+
+                export.data.append(c)
 
 
 class HydrogenTransportProblemDiscontinuous(HydrogenTransportProblem):
     interfaces: list[_subdomain.Interface]
     surface_to_volume: dict
     method_interface: str = "penalty"
+    subdomain_to_species: dict
 
     def __init__(
         self,
@@ -1031,6 +1045,7 @@ def __init__(
         )
         self.interfaces = interfaces or []
         self.surface_to_volume = surface_to_volume or {}
+        self.subdomain_to_species = {}  # maps subdomain to species defined in it
 
     def initialise(self):
         # check that all species have a list of F.VolumeSubdomain as this is
@@ -1162,13 +1177,13 @@ def define_function_spaces(
                 Defaults to 1.
         """
         # get number of species defined in the subdomain
-        all_species = [
+        self.subdomain_to_species[subdomain] = [
             species for species in self.species if subdomain in species.subdomains
         ]
 
         # instead of using the set function we use a list to keep the order
         unique_species = []
-        for species in all_species:
+        for species in self.subdomain_to_species[subdomain]:
             if species not in unique_species:
                 unique_species.append(species)
         nb_species = len(unique_species)
@@ -1364,7 +1379,7 @@ def mixed_term(u, v, n):
             h_1 = 2 * cr(res[1])
 
             all_mobile_species = [spe for spe in self.species if spe.mobile]
-
+            # TODO only do this if the species in defined in both domains of the interface?
             for H in all_mobile_species:
                 v_b = H.subdomain_to_test_function[subdomain_0](res[0])
                 v_t = H.subdomain_to_test_function[subdomain_1](res[1])
@@ -1426,17 +1441,13 @@ def mixed_term(u, v, n):
                         interface.id
                     ) - 0.5 * mixed_term(
                         v_b, (u_b / K_b - u_t / K_t), n_0
-                    ) * dInterface(
-                        interface.id
-                    )
+                    ) * dInterface(interface.id)
 
                     F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface(
                         interface.id
                     ) - 0.5 * mixed_term(
                         v_t, (u_b / K_b - u_t / K_t), n_0
-                    ) * dInterface(
-                        interface.id
-                    )
+                    ) * dInterface(interface.id)
                     F_0 += (
                         2
                         * gamma
@@ -1631,6 +1642,25 @@ def post_processing(self):
                 if export.filename is not None:
                     export.write(t=float(self.t))
 
+            elif isinstance(export, exports.Profile1DExport):
+                assert export.subdomain, (
+                    "Profile1DExport requires a subdomain to be set"
+                )
+                u = export.subdomain.u
+                if export._dofs is None:
+                    index = self.subdomain_to_species[export.subdomain].index(
+                        export.field
+                    )
+                    V0, export._dofs = 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
+
+                c = u.x.array[export._dofs][export._sort_coords]
+
+                export.data.append(c)
+
     def iterate(self):
         """Iterates the model for a given time step"""
         if self.show_progress_bar:

test/test_profile.py

@@ -0,0 +1,127 @@
+import festim as F
+import numpy as np
+
+
+def test_profile():
+    my_model = F.HydrogenTransportProblem()
+
+    protium = F.Species("H")
+    deuterium = F.Species("D")
+    tritium = F.Species("T")
+    my_model.species = [protium, deuterium, tritium]
+
+    my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100))
+
+    left_surf = F.SurfaceSubdomain1D(id=1, x=0)
+    right_surf = F.SurfaceSubdomain1D(id=2, x=1)
+
+    # assumes the same diffusivity for all species
+    material = F.Material(D_0=1, E_D=0)
+
+    vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material)
+
+    my_model.subdomains = [vol, left_surf, right_surf]
+
+    my_model.boundary_conditions = [
+        # Protium BCs
+        F.FixedConcentrationBC(left_surf, value=10, species=protium),
+        F.FixedConcentrationBC(right_surf, value=0, species=protium),
+        # Deuterium BCs
+        F.FixedConcentrationBC(left_surf, value=5, species=deuterium),
+        F.FixedConcentrationBC(right_surf, value=0, species=deuterium),
+        # Tritium BCs
+        F.FixedConcentrationBC(left_surf, value=0, species=tritium),
+        F.FixedConcentrationBC(right_surf, value=2, species=tritium),
+    ]
+
+    my_model.temperature = 300
+
+    my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5)
+
+    my_model.settings.stepsize = F.Stepsize(1)
+
+    my_model.exports = [
+        F.Profile1DExport(protium),
+        F.Profile1DExport(deuterium),
+    ]
+
+    my_model.initialise()
+    my_model.run()
+
+    assert my_model.exports[0].x is not None
+    assert my_model.exports[1].x is not None
+    assert len(my_model.exports[0].data) > 0
+    assert len(my_model.exports[1].data) > 0
+
+
+def test_profile_discontinuous():
+    my_model = F.HydrogenTransportProblemDiscontinuous()
+
+    protium = F.Species("H")
+    deuterium = F.Species("D")
+    tritium = F.Species("T")
+    my_model.species = [protium, deuterium, tritium]
+
+    vertices_left = np.linspace(0, 0.5, 50)
+    vertices_right = np.linspace(0.5, 1, 50)
+    vertices = np.concatenate((vertices_left, vertices_right))
+
+    my_model.mesh = F.Mesh1D(vertices)
+
+    left_surf = F.SurfaceSubdomain1D(id=1, x=0)
+    right_surf = F.SurfaceSubdomain1D(id=2, x=1)
+
+    # assumes the same diffusivity for all species
+    material_left = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0)
+    material_right = F.Material(D_0=1, E_D=0, K_S_0=2, E_K_S=0)
+
+    vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=material_left)
+    vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=material_right)
+
+    my_model.interfaces = [
+        F.Interface(id=3, subdomains=[vol1, vol2], penalty_term=1000)
+    ]
+
+    my_model.subdomains = [vol1, vol2, left_surf, right_surf]
+
+    for spe in my_model.species:
+        spe.subdomains = [vol1, vol2]
+
+    my_model.surface_to_volume = {
+        left_surf: vol1,
+        right_surf: vol2,
+    }
+
+    my_model.boundary_conditions = [
+        # Protium BCs
+        F.FixedConcentrationBC(left_surf, value=10, species=protium),
+        F.FixedConcentrationBC(right_surf, value=0, species=protium),
+        # Deuterium BCs
+        F.FixedConcentrationBC(left_surf, value=5, species=deuterium),
+        F.FixedConcentrationBC(right_surf, value=0, species=deuterium),
+        # Tritium BCs
+        F.FixedConcentrationBC(left_surf, value=0, species=tritium),
+        F.FixedConcentrationBC(right_surf, value=2, species=tritium),
+    ]
+
+    my_model.temperature = 300
+
+    my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=0.5)
+
+    my_model.settings.stepsize = F.Stepsize(0.1)
+
+    my_model.exports = [
+        F.Profile1DExport(field=protium, subdomain=vol1),
+        F.Profile1DExport(field=deuterium, subdomain=vol1),
+        F.Profile1DExport(field=tritium, subdomain=vol1),
+        F.Profile1DExport(field=protium, subdomain=vol2),
+        F.Profile1DExport(field=deuterium, subdomain=vol2),
+        F.Profile1DExport(field=tritium, subdomain=vol2),
+    ]
+
+    my_model.initialise()
+    my_model.run()
+
+    for export in my_model.exports:
+        assert export.x is not None
+        assert len(export.data) > 0