From 450b50eb963591b7022c94f34725fd9a63abeb2a Mon Sep 17 00:00:00 2001 From: rekomodo Date: Thu, 29 Aug 2024 16:49:53 -0400 Subject: [PATCH 1/9] initial implementation of TotalVolumeCylindrical --- .../derived_quantities/total_volume.py | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/festim/exports/derived_quantities/total_volume.py b/festim/exports/derived_quantities/total_volume.py index fd770dde1..b82dbc78e 100644 --- a/festim/exports/derived_quantities/total_volume.py +++ b/festim/exports/derived_quantities/total_volume.py @@ -1,5 +1,6 @@ from festim import VolumeQuantity import fenics as f +import numpy as np class TotalVolume(VolumeQuantity): @@ -58,3 +59,54 @@ def title(self): def compute(self): return f.assemble(self.function * self.dx(self.volume)) + + +class TotalVolumeCylindrical(TotalVolume): + """ + Computes the total value of a field in a given volume + int(f r dx) + + Args: + field (str, int): the field ("solute", 0, 1, "T", "retention") + volume (int): the volume 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). + + Attributes: + field (str, int): the field ("solute", 0, 1, "T", "retention") + volume (int): the volume id + export_unit (str): the unit of the derived quantity for exporting + title (str): the title of the derived quantity + show_units (bool): show the units in the title in the derived quantities + file + function (dolfin.function.function.Function): the solution function of + the hydrogen solute field + + .. note:: + units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen + concentration and K m in 1D, K m2 in 2D and K m3 in 3D domains for temperature + + """ + + def __init__(self, field, volume, azimuth_range = (0, 2*np.pi)): + super().__init__(field=field, volume=volume) + 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 + + @property + def allowed_meshes(self): + return ["cylindrical"] + + def compute(self): + theta_0, theta_1 = self.azimuth_range + return (theta_1 - theta_0) * f.assemble(self.function * self.r * self.dx(self.volume)) \ No newline at end of file From d8748c4bb9fdbfa7137f4bde3f4f07a12bc83b8b Mon Sep 17 00:00:00 2001 From: rekomodo Date: Thu, 29 Aug 2024 16:52:32 -0400 Subject: [PATCH 2/9] added TotalVolumeCylindrical to imports in __init__.py --- festim/__init__.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/festim/__init__.py b/festim/__init__.py index 10ce0105c..4026984a9 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -81,7 +81,10 @@ from .exports.derived_quantities.minimum_surface import MinimumSurface from .exports.derived_quantities.maximum_surface import MaximumSurface from .exports.derived_quantities.total_surface import TotalSurface -from .exports.derived_quantities.total_volume import TotalVolume +from .exports.derived_quantities.total_volume import ( + TotalVolume, + TotalVolumeCylindrical +) from .exports.derived_quantities.average_surface import AverageSurface from .exports.derived_quantities.point_value import PointValue from .exports.derived_quantities.adsorbed_hydrogen import AdsorbedHydrogen From e48463f6375c4c543ea45424b2ef9e6a3cc4d4e6 Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 13:36:33 -0400 Subject: [PATCH 3/9] tests for cylindrical total volume --- .../test_total_volume.py | 60 ++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py index 5fd3a8cb3..156f2efd7 100644 --- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_total_volume.py @@ -1,8 +1,9 @@ -from festim import TotalVolume +from festim import TotalVolume, TotalVolumeCylindrical import fenics as f import pytest from .tools import c_1D, c_2D, c_3D import pytest +import numpy as np @pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)]) @@ -46,6 +47,63 @@ def test_minimum(self): assert produced == expected +@pytest.mark.parametrize("radius", [2, 3]) +@pytest.mark.parametrize("r0", [0, 2]) +@pytest.mark.parametrize("height", [2, 3]) +@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi/2)]) +def test_compute_cylindrical(r0, radius, height, azimuth_range): + """ + Test that TotalVolumeCylindrical computes the volume correctly on a cylinder + + Args: + r0 (float): internal radius + radius (float): cylinder radius + height (float): cylinder height + azimuth_range (tuple): range of the azimuthal angle + soret (bool): if True, the Soret effect will be set + """ + # 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) + + volume = 1 + my_total = TotalVolumeCylindrical("solute", volume, azimuth_range) + + dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers) + + V = f.FunctionSpace(mesh_fenics, "P", 1) + r = f.interpolate(f.Expression("x[0]", degree=1), V) + c = 2*r + + my_total.dx = dx + my_total.function = c + my_total.r = f.Expression("x[0]", degree=1) + + az0, az1 = azimuth_range + + expected_value = f.assemble((az1 - az0) * c * r * dx(volume)) + computed_value = my_total.compute() + + assert np.isclose(expected_value, computed_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): + TotalVolumeCylindrical("solute", 1, azimuth_range=azimuth_range) + + @pytest.mark.parametrize( "function, field, expected_title", [ From 2989b2a1bbff466378d74f778e7d467349056610 Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 13:46:27 -0400 Subject: [PATCH 4/9] initial impl of TotalVolumeSpherical --- .../derived_quantities/total_volume.py | 68 ++++++++++++++++++- 1 file changed, 67 insertions(+), 1 deletion(-) diff --git a/festim/exports/derived_quantities/total_volume.py b/festim/exports/derived_quantities/total_volume.py index b82dbc78e..a9299216b 100644 --- a/festim/exports/derived_quantities/total_volume.py +++ b/festim/exports/derived_quantities/total_volume.py @@ -109,4 +109,70 @@ def allowed_meshes(self): def compute(self): theta_0, theta_1 = self.azimuth_range - return (theta_1 - theta_0) * f.assemble(self.function * self.r * self.dx(self.volume)) \ No newline at end of file + return (theta_1 - theta_0) * f.assemble(self.function * self.r * self.dx(self.volume)) + + +class TotalVolumeSpherical(TotalVolume): + """ + Computes the total value of a field in a given volume + int(f r dx) + + Args: + field (str, int): the field ("solute", 0, 1, "T", "retention") + volume (int): the volume 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). + + Attributes: + field (str, int): the field ("solute", 0, 1, "T", "retention") + volume (int): the volume id + export_unit (str): the unit of the derived quantity for exporting + title (str): the title of the derived quantity + show_units (bool): show the units in the title in the derived quantities + file + function (dolfin.function.function.Function): the solution function of + the hydrogen solute field + + .. note:: + units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen + concentration and K m in 1D, K m2 in 2D and K m3 in 3D domains for temperature + + """ + + def __init__(self, field, volume, azimuth_range = (0, np.pi), polar_range = (-np.pi, np.pi)): + super().__init__(field=field, volume=volume) + self.r = None + self.azimuth_range = azimuth_range + self.polar_range = polar_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] > np.pi: + raise ValueError("Azimuthal range must be between 0 and pi") + self._azimuth_range = value + + @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 allowed_meshes(self): + return ["spherical"] + + def compute(self): + phi_0, phi_1 = self.azimuth_range + theta_0, theta_1 = self.polar_range + + return (phi_1 - phi_0) * (theta_1 - theta_0) * f.assemble(self.function * self.r**2 * self.dx(self.volume)) \ No newline at end of file From 79f5b1660fd7af53ae81a8dfd5d62e7c7feb4dec Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 13:46:46 -0400 Subject: [PATCH 5/9] ran black on total_volume.py --- .../derived_quantities/total_volume.py | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/festim/exports/derived_quantities/total_volume.py b/festim/exports/derived_quantities/total_volume.py index a9299216b..fda231981 100644 --- a/festim/exports/derived_quantities/total_volume.py +++ b/festim/exports/derived_quantities/total_volume.py @@ -88,7 +88,7 @@ class TotalVolumeCylindrical(TotalVolume): """ - def __init__(self, field, volume, azimuth_range = (0, 2*np.pi)): + def __init__(self, field, volume, azimuth_range=(0, 2 * np.pi)): super().__init__(field=field, volume=volume) self.r = None self.azimuth_range = azimuth_range @@ -106,11 +106,13 @@ def azimuth_range(self, value): @property def allowed_meshes(self): return ["cylindrical"] - + def compute(self): theta_0, theta_1 = self.azimuth_range - return (theta_1 - theta_0) * f.assemble(self.function * self.r * self.dx(self.volume)) - + return (theta_1 - theta_0) * f.assemble( + self.function * self.r * self.dx(self.volume) + ) + class TotalVolumeSpherical(TotalVolume): """ @@ -141,7 +143,9 @@ class TotalVolumeSpherical(TotalVolume): """ - def __init__(self, field, volume, azimuth_range = (0, np.pi), polar_range = (-np.pi, np.pi)): + def __init__( + self, field, volume, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) + ): super().__init__(field=field, volume=volume) self.r = None self.azimuth_range = azimuth_range @@ -170,9 +174,13 @@ def polar_range(self, value): @property def allowed_meshes(self): return ["spherical"] - + def compute(self): phi_0, phi_1 = self.azimuth_range theta_0, theta_1 = self.polar_range - return (phi_1 - phi_0) * (theta_1 - theta_0) * f.assemble(self.function * self.r**2 * self.dx(self.volume)) \ No newline at end of file + return ( + (phi_1 - phi_0) + * (theta_1 - theta_0) + * f.assemble(self.function * self.r**2 * self.dx(self.volume)) + ) From 7d9aa8822b48966e4e7d95766b8771ed9aef6090 Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 13:49:16 -0400 Subject: [PATCH 6/9] fix cylinder test docstring --- .../test_exports/test_derived_quantities/test_total_volume.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py index 156f2efd7..98174d5ca 100644 --- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_total_volume.py @@ -60,7 +60,6 @@ def test_compute_cylindrical(r0, radius, height, azimuth_range): radius (float): cylinder radius height (float): cylinder height azimuth_range (tuple): range of the azimuthal angle - soret (bool): if True, the Soret effect will be set """ # creating a mesh with FEniCS r1 = r0 + radius From a53c22e8f7e34e726b21aeff4f53ecfb53366f67 Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 14:00:04 -0400 Subject: [PATCH 7/9] added to exports + made tests for TotalVolumeSpherical --- festim/__init__.py | 3 +- .../test_total_volume.py | 67 ++++++++++++++++++- 2 files changed, 68 insertions(+), 2 deletions(-) diff --git a/festim/__init__.py b/festim/__init__.py index 4026984a9..9e397e764 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -83,7 +83,8 @@ from .exports.derived_quantities.total_surface import TotalSurface from .exports.derived_quantities.total_volume import ( TotalVolume, - TotalVolumeCylindrical + TotalVolumeCylindrical, + TotalVolumeSpherical ) from .exports.derived_quantities.average_surface import AverageSurface from .exports.derived_quantities.point_value import PointValue diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py index 98174d5ca..d5f6003a4 100644 --- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_total_volume.py @@ -1,4 +1,4 @@ -from festim import TotalVolume, TotalVolumeCylindrical +from festim import TotalVolume, TotalVolumeCylindrical, TotalVolumeSpherical import fenics as f import pytest from .tools import c_1D, c_2D, c_3D @@ -103,6 +103,71 @@ def test_azimuthal_range_cylindrical(azimuth_range): TotalVolumeCylindrical("solute", 1, azimuth_range=azimuth_range) +@pytest.mark.parametrize("radius", [2, 3]) +@pytest.mark.parametrize("r0", [0, 2]) +@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi/2)]) +@pytest.mark.parametrize("polar_range", [(-np.pi, np.pi), (np.pi/2, -np.pi/2)]) +def test_compute_spherical(r0, radius, azimuth_range, polar_range): + """ + Test that TotalVolumeSpherical computes the volume correctly on a spherical mesh + + Args: + r0 (float): internal radius + radius (float): sphere radius + azimuth_range (tuple): range of the azimuthal angle + """ + # creating a mesh with FEniCS + r1 = r0 + radius + mesh_fenics = f.IntervalMesh(100, 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) + + volume = 1 + my_total = TotalVolumeSpherical("solute", volume, azimuth_range, polar_range) + + dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers) + + V = f.FunctionSpace(mesh_fenics, "P", 1) + r = f.interpolate(f.Expression("x[0]", degree=1), V) + c = 2*r + + my_total.dx = dx + my_total.function = c + my_total.r = f.Expression("x[0]", degree=1) + + az0, az1 = azimuth_range + th0, th1 = polar_range + + expected_value = f.assemble((az1 - az0) * (th1 - th0) * c * r**2 * dx(volume)) + computed_value = my_total.compute() + + assert np.isclose(expected_value, computed_value) + + +@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): + TotalVolumeSpherical("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): + TotalVolumeSpherical("solute", 1, polar_range=polar_range) + + @pytest.mark.parametrize( "function, field, expected_title", [ From aadb1efbb3ed06b1644fe597addb43ae00fd910f Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 14:00:58 -0400 Subject: [PATCH 8/9] rank black on tests --- .../test_derived_quantities/test_total_volume.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py index d5f6003a4..87eec1602 100644 --- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_total_volume.py @@ -50,7 +50,7 @@ def test_minimum(self): @pytest.mark.parametrize("radius", [2, 3]) @pytest.mark.parametrize("r0", [0, 2]) @pytest.mark.parametrize("height", [2, 3]) -@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi/2)]) +@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi / 2)]) def test_compute_cylindrical(r0, radius, height, azimuth_range): """ Test that TotalVolumeCylindrical computes the volume correctly on a cylinder @@ -78,7 +78,7 @@ def test_compute_cylindrical(r0, radius, height, azimuth_range): V = f.FunctionSpace(mesh_fenics, "P", 1) r = f.interpolate(f.Expression("x[0]", degree=1), V) - c = 2*r + c = 2 * r my_total.dx = dx my_total.function = c @@ -105,8 +105,8 @@ def test_azimuthal_range_cylindrical(azimuth_range): @pytest.mark.parametrize("radius", [2, 3]) @pytest.mark.parametrize("r0", [0, 2]) -@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi/2)]) -@pytest.mark.parametrize("polar_range", [(-np.pi, np.pi), (np.pi/2, -np.pi/2)]) +@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi / 2)]) +@pytest.mark.parametrize("polar_range", [(-np.pi, np.pi), (np.pi / 2, -np.pi / 2)]) def test_compute_spherical(r0, radius, azimuth_range, polar_range): """ Test that TotalVolumeSpherical computes the volume correctly on a spherical mesh @@ -131,7 +131,7 @@ def test_compute_spherical(r0, radius, azimuth_range, polar_range): V = f.FunctionSpace(mesh_fenics, "P", 1) r = f.interpolate(f.Expression("x[0]", degree=1), V) - c = 2*r + c = 2 * r my_total.dx = dx my_total.function = c From 9dfb219ded55fc2617abe0704e94f2a6db71fbf9 Mon Sep 17 00:00:00 2001 From: rekomodo Date: Sun, 1 Sep 2024 14:06:43 -0400 Subject: [PATCH 9/9] ran black on __init__.py --- festim/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/__init__.py b/festim/__init__.py index 9e397e764..422b6f7bf 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -84,7 +84,7 @@ from .exports.derived_quantities.total_volume import ( TotalVolume, TotalVolumeCylindrical, - TotalVolumeSpherical + TotalVolumeSpherical, ) from .exports.derived_quantities.average_surface import AverageSurface from .exports.derived_quantities.point_value import PointValue