From 973504e735e6723634ad2e85729d02ee621d1c59 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 24 Jan 2024 17:40:51 -0500 Subject: [PATCH 1/9] added tests --- .../test_surface_flux.py | 132 +++++++++++++++++- 1 file changed, 131 insertions(+), 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py index 582f9c3f6..815db9f67 100644 --- a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py +++ b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py @@ -1,5 +1,8 @@ -from festim import SurfaceFlux, R +from festim import SurfaceFlux, R, SurfaceFluxCylindrical, SurfaceFluxSpherical import fenics as f +import math +import numpy as np +import pytest def test_title_H(): @@ -80,3 +83,130 @@ def test_h_flux_with_soret(self): ) flux = self.my_h_flux.compute(soret=True) assert flux == expected_flux + + +@pytest.mark.parametrize("radius", [2, 3]) +@pytest.mark.parametrize("r0", [0, 2]) +@pytest.mark.parametrize("height", [2, 3]) +@pytest.mark.parametrize("c_top", [2, 3]) +@pytest.mark.parametrize("c_bottom", [2, 3]) +def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): + """ + Test that SurfaceFluxCylindrical computes the flux correctly on a hollow cylinder + + Args: + r0 (float): internal radius + radius (float): cylinder radius + height (float): cylinder height + c_top (float): concentration top + c_bottom (float): concentration bottom + """ + # creating a mesh with FEniCS + r1 = r0 + radius + z0 = 0 + z1 = z0 + height + mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10) + + # marking physical groups (volumes and surfaces) + volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim()) + volume_markers.set_all(1) + + bot_surface = f.CompiledSubDomain( + f"on_boundary && near(x[1], {z0}, tol)", tol=1e-14 + ) + top_surface = f.CompiledSubDomain( + f"on_boundary && near(x[1], {z1}, tol)", tol=1e-14 + ) + + surface_markers = f.MeshFunction( + "size_t", mesh_fenics, mesh_fenics.topology().dim() - 1 + ) + surface_markers.set_all(0) + ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers) + # Surface ids + bottom_id = 1 + top_id = 2 + bot_surface.mark(surface_markers, bottom_id) + top_surface.mark(surface_markers, top_id) + + my_flux = SurfaceFluxCylindrical("solute", top_id) + my_flux.D = f.Constant(2) + V = f.FunctionSpace(mesh_fenics, "P", 1) + expr = f.Expression( + f"{c_bottom} + {(c_top - c_bottom)/(z1-z0)} * x[1]", + degree=1, + ) + my_flux.function = f.interpolate(expr, V) + + my_flux.n = f.FacetNormal(mesh_fenics) + my_flux.ds = ds + + computed_value = my_flux.compute() + expected_value = ( + -2 + * math.pi + * float(my_flux.D) + * (c_bottom - c_top) + / height + * (0.5 * r1**2 - 0.5 * r0**2) + ) + assert np.isclose(computed_value, expected_value) + + +@pytest.mark.parametrize("r0", [1, 2]) +@pytest.mark.parametrize("radius", [1, 2]) +@pytest.mark.parametrize("c_left", [3, 4]) +@pytest.mark.parametrize("c_right", [1, 2]) +def test_compute_spherical(r0, radius, c_left, c_right): + """ + Test that SurfaceFluxSpherical computes the flux correctly on a hollow sphere + + Args: + r0 (float): internal radius + radius (float): sphere radius + c_left (float): concentration left + c_right (float): concentration right + """ + # creating a mesh with FEniCS + r1 = r0 + radius + mesh_fenics = f.IntervalMesh(10, r0, r1) + + # marking physical groups (volumes and surfaces) + volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim()) + volume_markers.set_all(1) + + left_surface = f.CompiledSubDomain( + f"on_boundary && near(x[0], {r0}, tol)", tol=1e-14 + ) + right_surface = f.CompiledSubDomain( + f"on_boundary && near(x[0], {r1}, tol)", tol=1e-14 + ) + + surface_markers = f.MeshFunction( + "size_t", mesh_fenics, mesh_fenics.topology().dim() - 1 + ) + surface_markers.set_all(0) + # Surface ids + left_id = 1 + right_id = 2 + left_surface.mark(surface_markers, left_id) + right_surface.mark(surface_markers, right_id) + ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers) + + my_flux = SurfaceFluxSpherical("solute", right_id) + my_flux.D = f.Constant(2) + V = f.FunctionSpace(mesh_fenics, "P", 1) + expr = f.Expression( + f"{c_left} + {(c_right - c_left)/(r1-r0)} * x[0]", + degree=1, + ) + my_flux.function = f.interpolate(expr, V) + + my_flux.n = f.FacetNormal(mesh_fenics) + my_flux.ds = ds + + computed_value = my_flux.compute() + # expected value is the integral of the flux over the surface + flux_value = float(my_flux.D) * (c_left - c_right) / (r1 - r0) + expected_value = -4 * math.pi * flux_value * r1**2 + assert np.isclose(computed_value, expected_value) From 0e2c9711da8d5acbeee054c85f7b111c6ad6cd51 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 24 Jan 2024 17:41:04 -0500 Subject: [PATCH 2/9] implemented Cylindrical and spherical SurfaceFlux --- festim/__init__.py | 6 +- .../derived_quantities/surface_flux.py | 71 ++++++++++++++++++- 2 files changed, 74 insertions(+), 3 deletions(-) diff --git a/festim/__init__.py b/festim/__init__.py index 2d1134f65..9d90142ab 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -67,7 +67,11 @@ VolumeQuantity, SurfaceQuantity, ) -from .exports.derived_quantities.surface_flux import SurfaceFlux +from .exports.derived_quantities.surface_flux import ( + SurfaceFlux, + SurfaceFluxCylindrical, + SurfaceFluxSpherical, +) from .exports.derived_quantities.hydrogen_flux import HydrogenFlux from .exports.derived_quantities.thermal_flux import ThermalFlux from .exports.derived_quantities.average_volume import AverageVolume diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index 38d66e419..bb0da7171 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/festim/exports/derived_quantities/surface_flux.py @@ -1,5 +1,6 @@ from festim import SurfaceQuantity, R import fenics as f +import numpy as np class SurfaceFlux(SurfaceQuantity): @@ -7,14 +8,17 @@ def __init__(self, field, surface) -> None: super().__init__(field=field, surface=surface) self.title = "Flux surface {}: {}".format(self.surface, self.field) - def compute(self, soret=False): + @property + def prop(self): field_to_prop = { "0": self.D, "solute": self.D, 0: self.D, "T": self.thermal_cond, } - self.prop = field_to_prop[self.field] + return field_to_prop[self.field] + + def compute(self, soret=False): flux = f.assemble( self.prop * f.dot(f.grad(self.function), self.n) * self.ds(self.surface) ) @@ -28,3 +32,66 @@ def compute(self, soret=False): * self.ds(self.surface) ) return flux + + +class SurfaceFluxCylindrical(SurfaceFlux): + def __init__(self, field, surface) -> None: + super().__init__(field, surface) + self.r = None + + def compute(self, soret=False): + if soret: + raise NotImplementedError( + "Soret effect not implemented for cylindrical coordinates" + ) + + if self.r is None: + mesh = ( + self.function.function_space().mesh() + ) # get the mesh from the function + rthetaz = f.SpatialCoordinate(mesh) # get the coordinates from the mesh + self.r = rthetaz[0] # only care about r here + + # dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr + # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz + # in both cases the expression with self.ds is the same + # we assume full cylinder theta = 2 pi + flux = f.assemble( + self.prop + * self.r + * f.dot(f.grad(self.function), self.n) + * self.ds(self.surface) + ) + theta = 2 * np.pi + flux *= theta + return flux + + +class SurfaceFluxSpherical(SurfaceFlux): + def __init__(self, field, surface) -> None: + super().__init__(field, surface) + self.r = None + + def compute(self, soret=False): + if soret: + raise NotImplementedError( + "Soret effect not implemented for spherical coordinates" + ) + + if self.r is None: + mesh = ( + self.function.function_space().mesh() + ) # get the mesh from the function + rthetaphi = f.SpatialCoordinate(mesh) # get the coordinates from the mesh + self.r = rthetaphi[0] # only care about r here + + # dS_r = r^2 sin(theta) dtheta dphi , assuming central symmetry dS_r = r^2 + flux = f.assemble( + self.prop + * self.r**2 + * f.dot(f.grad(self.function), self.n) + * self.ds(self.surface) + ) + # we assume a full sphere so multiply by 4 pi + flux *= 4 * np.pi + return flux From db1115032fa113db2c9905849a077db4723db644 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Sat, 2 Mar 2024 19:44:57 -0500 Subject: [PATCH 3/9] Add azimuth and polar range constraints to SurfaceFluxCylindrical and SurfaceFluxSpherical classes --- .../derived_quantities/surface_flux.py | 51 ++++++++++++++++--- 1 file changed, 44 insertions(+), 7 deletions(-) diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index 452837d31..32a8c5699 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/festim/exports/derived_quantities/surface_flux.py @@ -50,9 +50,20 @@ def compute(self, soret=False): class SurfaceFluxCylindrical(SurfaceFlux): - def __init__(self, field, surface) -> None: + def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None: super().__init__(field, surface) self.r = None + self.azimuth_range = azimuth_range + + @property + def azimuth_range(self): + return self._azimuth_range + + @azimuth_range.setter + def azimuth_range(self, value): + if value[0] < 0 or value[1] > 2 * np.pi: + raise ValueError("Azimuthal range must be between 0 and pi") + self._azimuth_range = value def compute(self, soret=False): if soret: @@ -77,15 +88,38 @@ def compute(self, soret=False): * f.dot(f.grad(self.function), self.n) * self.ds(self.surface) ) - theta = 2 * np.pi - flux *= theta + flux *= self.azimuth_cov return flux class SurfaceFluxSpherical(SurfaceFlux): - def __init__(self, field, surface) -> None: + def __init__( + self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) + ) -> None: super().__init__(field, surface) self.r = None + self.polar_range = polar_range + self.azimuth_range = azimuth_range + + @property + def polar_range(self): + return self._polar_range + + @polar_range.setter + def polar_range(self, value): + if value[0] < -np.pi or value[1] > np.pi: + raise ValueError("Polar range must be between - pi and pi") + self._polar_range = value + + @property + def azimuth_range(self): + return self._azimuth_range + + @azimuth_range.setter + def azimuth_range(self, value): + if value[0] < 0 or value[1] > np.pi: + raise ValueError("Azimuthal range must be between 0 and pi") + self._azimuth_range = value def compute(self, soret=False): if soret: @@ -100,13 +134,16 @@ def compute(self, soret=False): rthetaphi = f.SpatialCoordinate(mesh) # get the coordinates from the mesh self.r = rthetaphi[0] # only care about r here - # dS_r = r^2 sin(theta) dtheta dphi , assuming central symmetry dS_r = r^2 + # dS_r = r^2 sin(theta) dtheta dphi + # integral(f dS_r) = integral(f r^2 sin(theta) dtheta dphi) + # = (phi2 - phi1) * (-cos(theta2) + cos(theta1)) * f r^2 flux = f.assemble( self.prop * self.r**2 * f.dot(f.grad(self.function), self.n) * self.ds(self.surface) ) - # we assume a full sphere so multiply by 4 pi - flux *= 4 * np.pi + flux *= (self.polar_range[1] - self.polar_range[0]) * ( + -np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0]) + ) return flux From 4118209c2b19d1245bd8a346df3a97bb9eab46de Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Sat, 2 Mar 2024 20:06:32 -0500 Subject: [PATCH 4/9] fixed bug --- festim/exports/derived_quantities/surface_flux.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index 32a8c5699..94be6c076 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/festim/exports/derived_quantities/surface_flux.py @@ -88,7 +88,7 @@ def compute(self, soret=False): * f.dot(f.grad(self.function), self.n) * self.ds(self.surface) ) - flux *= self.azimuth_cov + flux *= self.azimuth_range[1] - self.azimuth_range[0] return flux From c1ff7dcfa3f94feba073a9a9ac1ea531a85e6a04 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Sat, 2 Mar 2024 20:06:44 -0500 Subject: [PATCH 5/9] modified the test to warn users --- festim/generic_simulation.py | 39 ++++++++++++++++++------------ test/simulation/test_initialise.py | 2 +- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py index fe149271f..e61db7b54 100644 --- a/festim/generic_simulation.py +++ b/festim/generic_simulation.py @@ -297,25 +297,32 @@ def initialise(self): self.h_transport_problem.initialise(self.mesh, self.materials, self.dt) - # raise warning if SurfaceFlux is used with non-cartesian meshes + # raise warning if the derived quantities don't match the type of mesh + # eg. SurfaceFlux is used with cylindrical mesh + all_types_quantities = [ + festim.MaximumSurface, + festim.MinimumSurface, + festim.MaximumVolume, + festim.MinimumVolume, + festim.PointValue, + ] # these quantities can be used with any mesh + allowed_quantities = { + "cartesian": [ + festim.SurfaceFlux, + festim.AverageSurface, + festim.AverageVolume, + ] + + all_types_quantities, + "cylindrical": [festim.SurfaceFluxCylindrical] + all_types_quantities, + "spherical": [festim.SurfaceFluxSpherical] + all_types_quantities, + } + 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": + for q in export.derived_quantities: + if not isinstance(q, tuple(allowed_quantities[self.mesh.type])): warnings.warn( - "Some derived quantities may not work as intended for non-cartesian meshes" + f"{type(q)} may not work as intended for {self.mesh.type} meshes" ) self.exports.initialise_derived_quantities( self.mesh.dx, self.mesh.ds, self.materials diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py index 6f6ad95bc..da78696f0 100644 --- a/test/simulation/test_initialise.py +++ b/test/simulation/test_initialise.py @@ -137,5 +137,5 @@ def test_cartesian_and_surface_flux_warning(quantity, sys): my_model.exports = [derived_quantities] # test - with pytest.warns(UserWarning, match="Some derived quantities .* non-cartesian"): + with pytest.warns(UserWarning, match=f"may not work as intended for {sys} meshes"): my_model.initialise() From 8fb7972d17ea61d53bc0150bb8200f55ab43bbde Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Sat, 2 Mar 2024 20:10:25 -0500 Subject: [PATCH 6/9] added docstrings --- .../derived_quantities/surface_flux.py | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index 94be6c076..c9f8df0f9 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/festim/exports/derived_quantities/surface_flux.py @@ -50,6 +50,22 @@ def compute(self, soret=False): class SurfaceFluxCylindrical(SurfaceFlux): + """ + 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 in cylindrical coordinates. + ds = r dr dtheta or ds = r dz dtheta + + Args: + field (str, int): the field ("solute", 0, 1, "T", "retention") + surface (int): the surface id + azimuth_range (tuple, optional): Range of the azimuthal angle + (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi). + """ + def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None: super().__init__(field, surface) self.r = None @@ -93,6 +109,24 @@ def compute(self, soret=False): class SurfaceFluxSpherical(SurfaceFlux): + """ + 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 in spherical coordinates. + ds = r^2 sin(theta) dtheta dphi + + Args: + field (str, int): the field ("solute", 0, 1, "T", "retention") + surface (int): the surface id + azimuth_range (tuple, optional): Range of the azimuthal angle + (phi) needs to be between 0 and pi. Defaults to (0, np.pi). + polar_range (tuple, optional): Range of the polar angle + (theta) needs to be between - pi and pi. Defaults to (-np.pi, np.pi). + """ + def __init__( self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) ) -> None: From cab76834f3071139264aa018f9673672d6095771 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Sat, 2 Mar 2024 20:12:18 -0500 Subject: [PATCH 7/9] added class to documentation --- docs/source/api/exports.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/source/api/exports.rst b/docs/source/api/exports.rst index a064f1fd6..92ee969c4 100644 --- a/docs/source/api/exports.rst +++ b/docs/source/api/exports.rst @@ -47,6 +47,14 @@ Exports :members: :show-inheritance: +.. autoclass:: SurfaceFluxCylindrical + :members: + :show-inheritance: + +.. autoclass:: SurfaceFluxSpherical + :members: + :show-inheritance: + .. autoclass:: HydrogenFlux :members: :show-inheritance: From f178287f7c4afe65ba1076f8c2ddea327a536048 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Sat, 2 Mar 2024 20:14:45 -0500 Subject: [PATCH 8/9] added a note in docstrings --- festim/exports/derived_quantities/surface_flux.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index c9f8df0f9..f4ea91daa 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/festim/exports/derived_quantities/surface_flux.py @@ -59,6 +59,8 @@ class SurfaceFluxCylindrical(SurfaceFlux): ds is the surface measure in cylindrical coordinates. ds = r dr dtheta or ds = r dz dtheta + Note: for particle fluxes J is given in H/s, for heat fluxes J is given in W + Args: field (str, int): the field ("solute", 0, 1, "T", "retention") surface (int): the surface id @@ -118,6 +120,8 @@ class SurfaceFluxSpherical(SurfaceFlux): ds is the surface measure in spherical coordinates. ds = r^2 sin(theta) dtheta dphi + Note: for particle fluxes J is given in H/s, for heat fluxes J is given in W + Args: field (str, int): the field ("solute", 0, 1, "T", "retention") surface (int): the surface id From 7d841f069d571103d80802ae11dca1afa7903f11 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 20 Mar 2024 08:48:20 -0400 Subject: [PATCH 9/9] more tests --- .../test_surface_flux.py | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py index 52f5b74f4..668d805e6 100644 --- a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py +++ b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py @@ -210,3 +210,54 @@ def test_compute_spherical(r0, radius, c_left, c_right): flux_value = float(my_flux.D) * (c_left - c_right) / (r1 - r0) expected_value = -4 * math.pi * flux_value * r1**2 assert np.isclose(computed_value, expected_value) + + +@pytest.mark.parametrize( + "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)] +) +def test_azimuthal_range_cylindrical(azimuth_range): + """ + Tests that an error is raised when the azimuthal range is out of bounds + """ + with pytest.raises(ValueError): + SurfaceFluxCylindrical("solute", 1, azimuth_range=azimuth_range) + + +def test_soret_raises_error(): + """ + Tests that an error is raised when the Soret effect is used with SurfaceFluxCylindrical + """ + my_flux = SurfaceFluxCylindrical("T", 1) + with pytest.raises(NotImplementedError): + my_flux.compute(soret=True) + + +@pytest.mark.parametrize( + "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)] +) +def test_azimuthal_range_spherical(azimuth_range): + """ + Tests that an error is raised when the azimuthal range is out of bounds + """ + with pytest.raises(ValueError): + SurfaceFluxSpherical("solute", 1, azimuth_range=azimuth_range) + + +@pytest.mark.parametrize( + "polar_range", [(0, 2 * np.pi), (-2 * np.pi, 0), (-2 * np.pi, 3 * np.pi)] +) +def test_polar_range_spherical(polar_range): + """ + Tests that an error is raised when the polar range is out of bounds + """ + with pytest.raises(ValueError): + SurfaceFluxSpherical("solute", 1, polar_range=polar_range) + + +def test_soret_raises_error_spherical(): + """ + Tests that an error is raised when the Soret effect is used with SurfaceFluxSpherical + """ + my_flux = SurfaceFluxSpherical("T", 1) + with pytest.raises(NotImplementedError): + my_flux.compute(soret=True)