bug fix for min/max values not working on discontinuous problem

.gitignore

@@ -113,10 +113,15 @@ venv.bak/
 *.xdmf
 *.h5
 *.txt
+*.png
+*.csv
 
 # Docs files
 docs/source/_build/
 # docs/source/_static/
 
+# patch files
+*.patch
+
 # Visual Studio Code settings
 .vscode/

src/festim/exports/maximum_surface.py

@@ -4,6 +4,7 @@
 import numpy as np
 
 from festim.exports.surface_quantity import SurfaceQuantity
+from festim.subdomain.volume_subdomain import VolumeSubdomain
 
 
 class MaximumSurface(SurfaceQuantity):
@@ -17,28 +18,59 @@ class MaximumSurface(SurfaceQuantity):
 
     Attributes:
         see `festim.SurfaceQuantity`
+        facet_meshtags: the facet meshtags of the parent mesh
+        volume: the volume subdomain the surface bounds. Set by the problem;
+            ``None`` outside `festim.HydrogenTransportProblemDiscontinuous`
     """
 
-    facet_meshtags: dolfinx.mesh.MeshTags
+    facet_meshtags: dolfinx.mesh.MeshTags | None = None
+    volume: VolumeSubdomain | None = None
 
     @property
     def title(self):
         return f"Maximum {self.field.name} surface {self.surface.id}"
 
+    @property
+    def is_submesh(self) -> bool:
+        """Whether the field's solution lives on the submesh of ``volume``, as in
+        ``festim.HydrogenTransportProblemDiscontinuous``. See issue #1191."""
+        return (
+            self.volume is not None
+            and self.volume in self.field.subdomain_to_post_processing_solution
+        )
+
+    @property
+    def solution(self):
+        if self.is_submesh:
+            return self.field.subdomain_to_post_processing_solution[self.volume]
+        return self.field.post_processing_solution
+
+    @property
+    def meshtags(self):
+        """Facet meshtags of whichever mesh ``solution`` lives on."""
+        return self.volume.ft if self.is_submesh else self.facet_meshtags
+
     def compute(self):
         """Computes the maximum value of the field on the defined surface subdomain, and
-        appends it to the data list."""
-        solution = self.field.post_processing_solution
-        entities = self.facet_meshtags.find(self.surface.id)
-        if isinstance(solution, dolfinx.fem.Function):
-            V = solution.function_space
-        else:
-            V = self.field.sub_function_space
+        appends it to the data list.
+        """
+        assert self.meshtags is not None, (
+            "facet meshtags must be set before computing the max surface value"
+        )
+        solution = self.solution
+        V = (
+            solution.function_space
+            if isinstance(solution, dolfinx.fem.Function)
+            else self.field.sub_function_space
+        )
         mesh = V.mesh
-
+        fdim = mesh.topology.dim - 1
+        mesh.topology.create_connectivity(fdim, mesh.topology.dim)
         dofs = dolfinx.fem.locate_dofs_topological(
-            V=V, entity_dim=mesh.topology.dim - 1, entities=entities
+            V=V, entity_dim=fdim, entities=self.meshtags.find(self.surface.id)
         )
+        values = solution.x.array[dofs]
 
-        self.value = mesh.comm.allreduce(np.max(solution.x.array[dofs]), op=MPI.MAX)
+        local_max = np.max(values) if values.size > 0 else -np.inf
+        self.value = mesh.comm.allreduce(local_max, op=MPI.MAX)
         self.data.append(self.value)

src/festim/exports/maximum_volume.py

@@ -17,29 +17,58 @@ class MaximumVolume(VolumeQuantity):
 
     Attributes:
         see `festim.VolumeQuantity`
+        volume_meshtags: the cell meshtags of the mesh the field is defined on. Is
+            ``None`` when the field is defined on a submesh that already coincides
+            with the volume subdomain (``festim.HydrogenTransportProblemDiscontinuous``)
     """
 
-    volume_meshtags: dolfinx.mesh.MeshTags
+    volume_meshtags: dolfinx.mesh.MeshTags | None = None
 
     @property
     def title(self):
         return f"Maximum {self.field.name} volume {self.volume.id}"
 
+    @property
+    def is_submesh(self) -> bool:
+        """Whether the field's solution lives on a submesh of ``volume`` itself, as in
+        ``festim.HydrogenTransportProblemDiscontinuous``. See issue #1191."""
+        return self.volume in self.field.subdomain_to_post_processing_solution
+
+    @property
+    def solution(self):
+        if self.is_submesh:
+            return self.field.subdomain_to_post_processing_solution[self.volume]
+        return self.field.post_processing_solution
+
     def compute(self):
-        """Computes the maximum value of solution function within the defined volume
-        subdomain, and appends it to the data list."""
-        solution = self.field.post_processing_solution
-        entities = self.volume_meshtags.find(self.volume.id)
+        """
+        Computes the maximum value of solution function within the defined volume
+        subdomain, and appends it to the data list.
+        """
+        solution = self.solution
+        V = (
+            solution.function_space
+            if isinstance(solution, dolfinx.fem.Function)
+            else self.field.sub_function_space
+        )
+        mesh = V.mesh
 
-        if isinstance(solution, dolfinx.fem.Function):
-            V = solution.function_space
+        if self.is_submesh:
+            # the mesh of the field is the volume subdomain itself (submesh case)
+            values = solution.x.array
         else:
-            V = self.field.sub_function_space
-        mesh = V.mesh
-        mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
-        dofs = dolfinx.fem.locate_dofs_topological(
-            V=V, entity_dim=mesh.topology.dim, entities=entities
-        )
+            assert self.volume_meshtags is not None, (
+                "volume_meshtags must be set before computing the max volume value"
+            )
+            entities = self.volume_meshtags.find(self.volume.id)
+            mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
+            dofs = dolfinx.fem.locate_dofs_topological(
+                V=V, entity_dim=mesh.topology.dim, entities=entities
+            )
+            values = solution.x.array[dofs]
+
+        # a process may hold no dof of the volume at all, np.max would then raise
+        local_max = np.max(values) if values.size > 0 else -np.inf
 
-        self.value = mesh.comm.allreduce(np.max(solution.x.array[dofs]), op=MPI.MAX)
+        self.value = mesh.comm.allreduce(local_max, op=MPI.MAX)
         self.data.append(self.value)

src/festim/exports/minimum_surface.py

@@ -4,6 +4,7 @@
 import numpy as np
 
 from festim.exports.surface_quantity import SurfaceQuantity
+from festim.subdomain.volume_subdomain import VolumeSubdomain
 
 
 class MinimumSurface(SurfaceQuantity):
@@ -17,27 +18,59 @@ class MinimumSurface(SurfaceQuantity):
 
     Attributes:
         see `festim.SurfaceQuantity`
+        facet_meshtags: the facet meshtags of the parent mesh
+        volume: the volume subdomain the surface bounds. Set by the problem;
+            ``None`` outside `festim.HydrogenTransportProblemDiscontinuous`
     """
 
-    facet_meshtags: dolfinx.mesh.MeshTags
+    facet_meshtags: dolfinx.mesh.MeshTags | None = None
+    volume: VolumeSubdomain | None = None
 
     @property
     def title(self):
         return f"Minimum {self.field.name} surface {self.surface.id}"
 
+    @property
+    def is_submesh(self) -> bool:
+        """Whether the field's solution lives on the submesh of ``volume``, as in
+        ``festim.HydrogenTransportProblemDiscontinuous``. See issue #1191."""
+        return (
+            self.volume is not None
+            and self.volume in self.field.subdomain_to_post_processing_solution
+        )
+
+    @property
+    def solution(self):
+        if self.is_submesh:
+            return self.field.subdomain_to_post_processing_solution[self.volume]
+        return self.field.post_processing_solution
+
+    @property
+    def meshtags(self):
+        """Facet meshtags of whichever mesh ``solution`` lives on."""
+        return self.volume.ft if self.is_submesh else self.facet_meshtags
+
     def compute(self):
         """Computes the minimum value of the field on the defined surface subdomain, and
-        appends it to the data list."""
-        solution = self.field.post_processing_solution
-        entities = self.facet_meshtags.find(self.surface.id)
-        if isinstance(solution, dolfinx.fem.Function):
-            V = solution.function_space
-        else:
-            V = self.field.sub_function_space
+        appends it to the data list.
+        """
+        assert self.meshtags is not None, (
+            "facet meshtags must be set before computing the min surface value"
+        )
+        solution = self.solution
+        V = (
+            solution.function_space
+            if isinstance(solution, dolfinx.fem.Function)
+            else self.field.sub_function_space
+        )
         mesh = V.mesh
+        fdim = mesh.topology.dim - 1
+        mesh.topology.create_connectivity(fdim, mesh.topology.dim)
         dofs = dolfinx.fem.locate_dofs_topological(
-            V=V, entity_dim=mesh.topology.dim - 1, entities=entities
+            V=V, entity_dim=fdim, entities=self.meshtags.find(self.surface.id)
         )
+        values = solution.x.array[dofs]
 
-        self.value = mesh.comm.allreduce(np.min(solution.x.array[dofs]), op=MPI.MIN)
+        local_min = np.min(values) if values.size > 0 else np.inf
+        self.value = mesh.comm.allreduce(local_min, op=MPI.MIN)
         self.data.append(self.value)

src/festim/exports/minimum_volume.py

@@ -17,29 +17,58 @@ class MinimumVolume(VolumeQuantity):
 
     Attributes:
         see `festim.VolumeQuantity`
+        volume_meshtags: the cell meshtags of the mesh the field is defined on. Is
+            ``None`` when the field is defined on a submesh that already coincides
+            with the volume subdomain (``festim.HydrogenTransportProblemDiscontinuous``)
     """
 
-    volume_meshtags: dolfinx.mesh.MeshTags
+    volume_meshtags: dolfinx.mesh.MeshTags | None = None
 
     @property
     def title(self):
         return f"Minimum {self.field.name} volume {self.volume.id}"
 
+    @property
+    def is_submesh(self) -> bool:
+        """Whether the field's solution lives on a submesh of ``volume`` itself, as in
+        ``festim.HydrogenTransportProblemDiscontinuous``. See issue #1191."""
+        return self.volume in self.field.subdomain_to_post_processing_solution
+
+    @property
+    def solution(self):
+        if self.is_submesh:
+            return self.field.subdomain_to_post_processing_solution[self.volume]
+        return self.field.post_processing_solution
+
     def compute(self):
-        """Computes the minimum value of solution function within the defined volume
-        subdomain, and appends it to the data list."""
-        solution = self.field.post_processing_solution
-        entities = self.volume_meshtags.find(self.volume.id)
+        """
+        Computes the minimum value of solution function within the defined volume
+        subdomain, and appends it to the data list.
+        """
+        solution = self.solution
+        V = (
+            solution.function_space
+            if isinstance(solution, dolfinx.fem.Function)
+            else self.field.sub_function_space
+        )
+        mesh = V.mesh
 
-        if isinstance(solution, dolfinx.fem.Function):
-            V = solution.function_space
+        if self.is_submesh:
+            # the mesh of the field is the volume subdomain itself (submesh case)
+            values = solution.x.array
         else:
-            V = self.field.sub_function_space
-        mesh = V.mesh
-        mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
-        dofs = dolfinx.fem.locate_dofs_topological(
-            V=V, entity_dim=mesh.topology.dim, entities=entities
-        )
+            assert self.volume_meshtags is not None, (
+                "volume meshtags must be set before computing the min volume value"
+            )
+            entities = self.volume_meshtags.find(self.volume.id)
+            mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
+            dofs = dolfinx.fem.locate_dofs_topological(
+                V=V, entity_dim=mesh.topology.dim, entities=entities
+            )
+            values = solution.x.array[dofs]
+
+        # a process may hold no dof of the volume at all, np.min would then raise
+        local_min = np.min(values) if values.size > 0 else np.inf
 
-        self.value = mesh.comm.allreduce(np.min(solution.x.array[dofs]), op=MPI.MIN)
+        self.value = mesh.comm.allreduce(local_min, op=MPI.MIN)
         self.data.append(self.value)

src/festim/hydrogen_transport_problem.py

@@ -3241,6 +3241,34 @@ def initialise_exports(self):
                         mesh=mesh,
                     )
 
+            # the extrema exports read the solution on the submesh of the volume
+            # subdomain their location belongs to, so give them that volume and the
+            # meshtags of the parent mesh, then check the species is defined there
+            is_extremum = isinstance(
+                export,
+                exports.MaximumVolume
+                | exports.MinimumVolume
+                | exports.MaximumSurface
+                | exports.MinimumSurface,
+            )
+            if is_extremum:
+                if isinstance(export, exports.SurfaceQuantity):
+                    export.facet_meshtags = self.facet_meshtags
+                    export.volume = self.surface_to_volume[export.surface]
+                    location = f"surface {export.surface.id}"
+                else:
+                    export.volume_meshtags = self.volume_meshtags
+                    location = f"volume {export.volume.id}"
+
+                # a field that is not a Species (a bare name) is skipped: it fails
+                # earlier, for unrelated reasons
+                if isinstance(export.field, _species.Species) and not export.is_submesh:
+                    raise ValueError(
+                        f"Cannot compute {export.title}: species "
+                        f"{export.field.name} is not defined in the volume subdomain "
+                        f"{export.volume.id} that {location} belongs to"
+                    )
+
             # reset the data and time for SurfaceQuantity and VolumeQuantity
             if isinstance(export, exports.DerivedQuantity):
                 export.t = []