Merge pull request #688 from RemDelaporteMathurin/warn-about-SurfaceF…

…lux-in-non-cartesian

Warn about derived quantities in non cartesian

festim/exports/derived_quantities/hydrogen_flux.py

@@ -2,5 +2,7 @@
 
 
 class HydrogenFlux(SurfaceFlux):
+    """Equivalent to SurfaceFlux("solute", ...)"""
+
     def __init__(self, surface) -> None:
         super().__init__(field="solute", surface=surface)

festim/exports/derived_quantities/surface_flux.py

@@ -4,6 +4,21 @@
 
 class SurfaceFlux(SurfaceQuantity):
     def __init__(self, field, surface) -> None:
+        """
+        Object to compute the flux J of a field u through a surface
+        J = integral(-prop * grad(u) . n ds)
+        where prop is the property of the field (D, thermal conductivity, etc)
+        u is the field
+        n is the normal vector of the surface
+        ds is the surface measure.
+
+        Note: this will probably won't work correctly for axisymmetric meshes.
+        Use only with cartesian coordinates.
+
+        Args:
+            field (str, int):  the field ("solute", 0, 1, "T", "retention")
+            surface (int): the surface id
+        """
         super().__init__(field=field, surface=surface)
         self.title = "Flux surface {}: {}".format(self.surface, self.field)
 

festim/generic_simulation.py

@@ -297,6 +297,26 @@ def initialise(self):
 
         self.h_transport_problem.initialise(self.mesh, self.materials, self.dt)
 
+        # raise warning if SurfaceFlux is used with non-cartesian meshes
+        for export in self.exports:
+            if isinstance(export, festim.DerivedQuantities):
+                allowed_quantities = (
+                    festim.MaximumSurface,
+                    festim.MinimumSurface,
+                    festim.MaximumVolume,
+                    festim.MinimumVolume,
+                    festim.PointValue,
+                )
+                if any(
+                    [
+                        not isinstance(q, allowed_quantities)
+                        for q in export.derived_quantities
+                    ]
+                ):
+                    if self.mesh.type != "cartesian":
+                        warnings.warn(
+                            "Some derived quantities may not work as intended for non-cartesian meshes"
+                        )
         self.exports.initialise_derived_quantities(
             self.mesh.dx, self.mesh.ds, self.materials
         )

test/simulation/test_initialise.py

@@ -1,5 +1,6 @@
 import festim as F
 from pathlib import Path
+import pytest
 
 
 def test_initialise_changes_nb_of_sources():
@@ -103,3 +104,38 @@ def test_TXTExport_times_added_to_milestones(tmpdir):
 
     # test
     assert my_model.dt.milestones == txt_export.times
+
+
+@pytest.mark.parametrize(
+    "quantity",
+    [
+        F.SurfaceFlux(field="solute", surface=1),
+        F.TotalVolume(field="solute", volume=1),
+        F.TotalSurface(field="solute", surface=1),
+        F.AverageSurface(field="solute", surface=1),
+        F.AverageVolume(field="solute", volume=1),
+        F.HydrogenFlux(surface=1),
+        F.ThermalFlux(surface=1),
+    ],
+)
+@pytest.mark.parametrize("sys", ["cylindrical", "spherical"])
+def test_cartesian_and_surface_flux_warning(quantity, sys):
+    """Creates a Simulation object and checks that, if either a cylindrical
+    or spherical meshes are given with a SurfaceFlux, a warning is raised.
+    """
+    # build
+    my_model = F.Simulation()
+    my_model.mesh = F.MeshFromVertices([1, 2, 3], type=sys)
+    my_model.materials = F.Material(id=1, D_0=1, E_D=0)
+    my_model.T = F.Temperature(100)
+    my_model.dt = F.Stepsize(initial_value=3)
+    my_model.settings = F.Settings(
+        absolute_tolerance=1e-10, relative_tolerance=1e-10, final_time=4
+    )
+
+    derived_quantities = F.DerivedQuantities([quantity])
+    my_model.exports = [derived_quantities]
+
+    # test
+    with pytest.warns(UserWarning, match="Some derived quantities .* non-cartesian"):
+        my_model.initialise()