diff --git a/festim/__init__.py b/festim/__init__.py index 3b199d1ad..1fe1d1397 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -74,14 +74,22 @@ ) from .exports.derived_quantities.hydrogen_flux import HydrogenFlux from .exports.derived_quantities.thermal_flux import ThermalFlux -from .exports.derived_quantities.average_volume import AverageVolume +from .exports.derived_quantities.average_volume import ( + AverageVolume, + AverageVolumeCylindrical, + AverageVolumeSpherical, +) from .exports.derived_quantities.maximum_volume import MaximumVolume from .exports.derived_quantities.minimum_volume import MinimumVolume 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.average_surface import AverageSurface +from .exports.derived_quantities.average_surface import ( + AverageSurface, + AverageSurfaceCylindrical, + AverageSurfaceSpherical, +) from .exports.derived_quantities.point_value import PointValue from .exports.derived_quantities.adsorbed_hydrogen import AdsorbedHydrogen diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index cb7bddad0..c56e0d793 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -30,7 +30,7 @@ def __init__(self, field, surface) -> None: @property def allowed_meshes(self): - return ["cartesian"] + return ["cartesian", "spherical"] @property def export_unit(self): @@ -51,3 +51,59 @@ def compute(self): return f.assemble(self.function * self.ds(self.surface)) / f.assemble( 1 * self.ds(self.surface) ) + + +class AverageSurfaceCylindrical(AverageSurface): + """ + Computes the average value of a field on a given surface + int(f ds) / int (1 * ds) + ds is the surface measure in cylindrical coordinates. + ds = r dr dtheta + + Args: + field (str, int): the field ("solute", 0, 1, "T", "retention") + surface (int): the surface id + + Attributes: + field (str, int): the field ("solute", 0, 1, "T", "retention") + surface (int): the surface id + 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 field + r (ufl.indexed.Indexed): the radius of the cylinder + + Notes: + Units are in H/m3 for hydrogen concentration and K for temperature + """ + + def __init__(self, field, surface) -> None: + super().__init__(field=field, surface=surface) + self.r = None + + @property + def allowed_meshes(self): + return ["cylindrical"] + + def compute(self): + + 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.dx is the same + + avg_surf = f.assemble( + self.function * self.r * self.ds(self.surface) + ) / f.assemble(1 * self.r * self.ds(self.surface)) + + return avg_surf + + +AverageSurfaceSpherical = AverageSurface diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py index a91967e3d..6428c0e5a 100644 --- a/festim/exports/derived_quantities/average_volume.py +++ b/festim/exports/derived_quantities/average_volume.py @@ -19,6 +19,7 @@ class AverageVolume(VolumeQuantity): file function (dolfin.function.function.Function): the solution function of the field + r (ufl.indexed.Indexed): the radius of the cylinder .. note:: Units are in H/m3 for hydrogen concentration and K for temperature @@ -50,3 +51,94 @@ def compute(self): return f.assemble(self.function * self.dx(self.volume)) / f.assemble( 1 * self.dx(self.volume) ) + + +class AverageVolumeCylindrical(AverageVolume): + """ + Computes the average value of a field in a given volume + int(f dx) / int (1 * dx) + dx is the volume measure in cylindrical coordinates. + dx = r dr 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") + volume (int): the volume id + + Attributes: + field (str, int): the field ("solute", 0, 1, "T", "retention") + volume (int): the volume id + 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 field + r (ufl.indexed.Indexed): the radius of the sphere + """ + + def __init__(self, field, volume) -> None: + super().__init__(field=field, volume=volume) + self.r = None + + @property + def allowed_meshes(self): + return ["cylindrical"] + + def compute(self): + + 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 + + avg_vol = f.assemble( + self.function * self.r * self.dx(self.volume) + ) / f.assemble(1 * self.r * self.dx(self.volume)) + + return avg_vol + + +class AverageVolumeSpherical(AverageVolume): + """ + Computes the average value of a field in a given volume + int(f dx) / int (1 * dx) + dx is the volume measure in cylindrical coordinates. + dx = rho 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") + volume (int): the volume id + 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 field + """ + + def __init__(self, field, volume) -> None: + super().__init__(field=field, volume=volume) + self.r = None + + @property + def allowed_meshes(self): + return ["spherical"] + + def compute(self): + + 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 + + avg_vol = f.assemble( + self.function * self.r**2 * self.dx(self.volume) + ) / f.assemble(1 * self.r**2 * self.dx(self.volume)) + + return avg_vol diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py index fdcb3a60d..c5c09aaa8 100644 --- a/test/simulation/test_initialise.py +++ b/test/simulation/test_initialise.py @@ -117,7 +117,6 @@ def test_TXTExport_times_added_to_milestones(tmpdir): 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), diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index 5c7f97a2e..17f47ac8d 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -1,6 +1,14 @@ -from festim import AverageSurface +from festim import ( + AverageSurface, + AverageSurfaceCylindrical, + AverageSurfaceSpherical, + x, + y, +) import fenics as f import pytest +import numpy as np +from sympy.printing import ccode @pytest.mark.parametrize("field, surface", [("solute", 1), ("T", 2)]) @@ -46,3 +54,147 @@ def test_h_average(self): ) computed = self.my_average.compute() assert computed == expected + + +@pytest.mark.parametrize("radius", [1, 4]) +@pytest.mark.parametrize("r0", [3, 5]) +@pytest.mark.parametrize("height", [2, 7]) +def test_compute_cylindrical(r0, radius, height): + """ + Test that AverageSurfaceCylindrical computes the value correctly on a hollow cylinder + + Args: + r0 (float): internal radius + radius (float): cylinder radius + height (float): cylinder height + """ + # creating a mesh with FEniCS + r1 = r0 + radius + z0, z1 = 0, height + + mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10) + + 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 + top_id = 2 + top_surface.mark(surface_markers, top_id) + + my_export = AverageSurfaceCylindrical("solute", top_id) + V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda r: 2 * r + + expr = f.Expression( + ccode(c_fun(x)), + degree=1, + ) + my_export.function = f.interpolate(expr, V) + my_export.ds = ds + + expected_value = 4 / 3 * (r1**3 - r0**3) / (r1**2 - r0**2) + + computed_value = my_export.compute() + + assert np.isclose(computed_value, expected_value) + + +@pytest.mark.parametrize("radius", [2, 4]) +@pytest.mark.parametrize("r0", [3, 5]) +def test_compute_spherical(r0, radius): + """ + Test that AverageSurfaceSpherical computes the average value correctly + on a hollow sphere + + Args: + r0 (float): internal radius + radius (float): sphere radius + """ + # creating a mesh with FEniCS + r1 = r0 + radius + mesh_fenics = f.IntervalMesh(10, r0, r1) + + # marking physical groups (volumes and surfaces) + 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 + right_id = 2 + right_surface.mark(surface_markers, right_id) + ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers) + + my_export = AverageSurfaceSpherical("solute", right_id) + V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda r: 4 * r + expr = f.Expression( + ccode(c_fun(x)), + degree=1, + ) + my_export.function = f.interpolate(expr, V) + + my_export.ds = ds + + expected_value = 4 * r1 + + computed_value = my_export.compute() + + assert np.isclose(computed_value, expected_value) + + +def test_average_surface_cylindrical_title_no_units_solute(): + """A simple test to check that the title is set correctly in + festim.AverageSurfaceCylindrical with a solute field without units""" + + my_export = AverageSurfaceCylindrical("solute", 4) + assert my_export.title == "Average solute surface 4 (H m-3)" + + +def test_average_surface_cylindrical_title_no_units_temperature(): + """A simple test to check that the title is set correctly in + festim.AverageSurfaceCylindrical with a T field without units""" + + my_export = AverageSurfaceCylindrical("T", 5) + assert my_export.title == "Average T surface 5 (K)" + + +def test_average_surface_spherical_title_no_units_solute(): + """A simple test to check that the title is set correctly in + festim.AverageSurfaceSpherical with a solute field without units""" + + my_export = AverageSurfaceSpherical("solute", 6) + assert my_export.title == "Average solute surface 6 (H m-3)" + + +def test_average_surface_spherical_title_no_units_temperature(): + """A simple test to check that the title is set correctly in + festim.AverageSurfaceSpherical with a T field without units""" + + my_export = AverageSurfaceSpherical("T", 9) + assert my_export.title == "Average T surface 9 (K)" + + +def test_avg_surf_cylindrical_allow_meshes(): + """A simple test to check cylindrical meshes are the only + meshes allowed when using AverageSurfaceCylindrical""" + + my_export = AverageSurfaceCylindrical("solute", 2) + + assert my_export.allowed_meshes == ["cylindrical"] + + +def test_avg_surf_spherical_allow_meshes(): + """A simple test to check spherical meshes are one of the + meshes allowed when using AverageSurfaceSpherical""" + + my_export = AverageSurfaceSpherical("T", 6) + + assert "spherical" in my_export.allowed_meshes diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index 14905ca07..fc4f83059 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -1,6 +1,8 @@ -from festim import AverageVolume +from festim import AverageVolume, AverageVolumeCylindrical, AverageVolumeSpherical, x import fenics as f import pytest +import numpy as np +from sympy.printing import ccode @pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)]) @@ -46,3 +48,127 @@ def test_h_average(self): ) computed = self.my_average.compute() assert computed == expected + + +@pytest.mark.parametrize("radius", [2, 4]) +@pytest.mark.parametrize("r0", [3, 5]) +@pytest.mark.parametrize("height", [2, 7]) +def test_compute_cylindrical(r0, radius, height): + """ + Test that AverageVolumeCylindrical computes the value correctly on a hollow cylinder + + Args: + r0 (float): internal radius + radius (float): cylinder radius + height (float): cylinder height + """ + # creating a mesh with FEniCS + r1 = r0 + radius + z0, z1 = 0, height + + mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10) + + volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim()) + volume_markers.set_all(1) + dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers) + + my_export = AverageVolumeCylindrical("solute", 1) + V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda r: 3 * r + expr = f.Expression(ccode(c_fun(x)), degree=1) + my_export.function = f.interpolate(expr, V) + my_export.dx = dx + + expected_value = 2 * (r1**3 - r0**3) / (r1**2 - r0**2) + + computed_value = my_export.compute() + + assert np.isclose(computed_value, expected_value) + + +@pytest.mark.parametrize("radius", [2, 4]) +@pytest.mark.parametrize("r0", [3, 5]) +def test_compute_spherical(r0, radius): + """ + Test that AverageVolumeSpherical computes the average value correctly + on a hollow sphere + + Args: + r0 (float): internal radius + radius (float): sphere radius + """ + # 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) + dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers) + + my_export = AverageVolumeSpherical("solute", 1) + V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda r: 4 * r + expr = f.Expression( + ccode(c_fun(x)), + degree=1, + ) + my_export.function = f.interpolate(expr, V) + + my_export.dx = dx + + expected_value = 3 * (r1**4 - r0**4) / (r1**3 - r0**3) + + computed_value = my_export.compute() + + assert np.isclose(computed_value, expected_value) + + +def test_average_volume_cylindrical_title_no_units_solute(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeCylindrical with a solute field without units""" + + my_export = AverageVolumeCylindrical("solute", 4) + assert my_export.title == "Average solute volume 4 (H m-3)" + + +def test_average_volume_cylindrical_title_no_units_temperature(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeCylindrical with a T field without units""" + + my_export = AverageVolumeCylindrical("T", 5) + assert my_export.title == "Average T volume 5 (K)" + + +def test_average_volume_spherical_title_no_units_solute(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeSpherical with a solute field without units""" + + my_export = AverageVolumeSpherical("solute", 6) + assert my_export.title == "Average solute volume 6 (H m-3)" + + +def test_average_volume_spherical_title_no_units_temperature(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeSpherical with a T field without units""" + + my_export = AverageVolumeSpherical("T", 9) + assert my_export.title == "Average T volume 9 (K)" + + +def test_avg_vol_cylindrical_allow_meshes(): + """A simple test to check cylindrical meshes are the only + meshes allowed when using AverageVolumeCylindrical""" + + my_export = AverageVolumeCylindrical("solute", 3) + + assert my_export.allowed_meshes == ["cylindrical"] + + +def test_avg_vol_spherical_allow_meshes(): + """A simple test to check cylindrical meshes are the only + meshes allowed when using AverageVolumeSpherical""" + + my_export = AverageVolumeSpherical("solute", 5) + + assert my_export.allowed_meshes == ["spherical"]