Export average surface temperature

src/festim/__init__.py

@@ -30,6 +30,7 @@
 from .exports.minimum_volume import MinimumVolume
 from .exports.surface_flux import SurfaceFlux
 from .exports.surface_quantity import SurfaceQuantity
+from .exports.surface_temperature import AverageSurfaceTemperature
 from .exports.total_surface import TotalSurface
 from .exports.total_volume import TotalVolume
 from .exports.volume_quantity import VolumeQuantity

src/festim/exports/__init__.py

@@ -9,12 +9,14 @@
     "AverageSurface",
     "AverageVolume",
     "TotalVolume",
+    "AverageSurfaceTemperature",
 ]
 
 from .average_surface import AverageSurface
 from .average_volume import AverageVolume
 from .surface_flux import SurfaceFlux
 from .surface_quantity import SurfaceQuantity
+from .surface_temperature import AverageSurfaceTemperature
 from .total_surface import TotalSurface
 from .total_volume import TotalVolume
 from .volume_quantity import VolumeQuantity

src/festim/exports/surface_temperature.py

@@ -0,0 +1,51 @@
+import csv
+from dolfinx import fem
+import ufl
+from .surface_quantity import SurfaceQuantity
+
+
+class AverageSurfaceTemperature(SurfaceQuantity):
+    """Exports the average temperature on a given surface.
+
+    Args:
+        surface (int or festim.SurfaceSubdomain): the surface subdomain
+        filename (str, optional): name of the file to which the average surface temperature is exported
+
+    Attributes:
+        temperature_field (fem.Constant or fem.Function): the temperature field
+        surface (int or festim.SurfaceSubdomain): the surface subdomain
+        filename (str): name of the file to which the surface temperature is exported
+        t (list): list of time values
+        data (list): list of average temperature values on the surface
+    """
+
+    def __init__(self, surface, filename: str = None) -> None:
+        self.surface = surface
+        self.filename = filename
+
+        self.temperature_field = None
+        self.t = []
+        self.data = []
+        self._first_time_export = True
+
+    @property
+    def title(self):
+        return f"Temperature surface {self.surface.id}"
+
+    def compute(self, ds):
+        """Computes the average temperature on the surface.
+
+        Args:
+            ds (ufl.Measure): surface measure of the model
+        """
+        temperature_field = self.temperature_field
+
+        surface_integral = fem.assemble_scalar(
+            fem.form(temperature_field * ds(self.surface.id))
+        )  # integral over surface
+
+        surface_area = fem.assemble_scalar(fem.form(1 * ds(self.surface.id)))
+
+        self.value = surface_integral / surface_area  # avg temp
+
+        self.data.append(self.value)

src/festim/hydrogen_transport_problem.py

@@ -369,6 +369,9 @@
         a string, find species object in self.species"""
 
         for export in self.exports:
+            if isinstance(export, exports.AverageSurfaceTemperature):
+                continue
+
             # if name of species is given then replace with species object
             if isinstance(export.field, list):
                 for idx, field in enumerate(export.field):
@@ -415,6 +418,9 @@
                 export.D = D
                 export.D_expr = D_expr
 
+            if isinstance(export, exports.AverageSurfaceTemperature):
+                export.temperature_field = self.temperature_fenics
+
             # reset the data and time for SurfaceQuantity and VolumeQuantity
             if isinstance(export, (exports.SurfaceQuantity, exports.VolumeQuantity)):
                 export.t = []
@@ -805,7 +811,12 @@
             if isinstance(export, exports.SurfaceQuantity):
                 if isinstance(
                     export,
-                    (exports.SurfaceFlux, exports.TotalSurface, exports.AverageSurface),
+                    (
+                        exports.SurfaceFlux,
+                        exports.TotalSurface,
+                        exports.AverageSurface,
+                        exports.AverageSurfaceTemperature,
+                    ),
                 ):
                     export.compute(
                         self.ds,
@@ -818,6 +829,7 @@
                 # if filename given write export data to file
                 if export.filename is not None:
                     export.write(t=float(self.t))
+
             elif isinstance(export, exports.VolumeQuantity):
                 if isinstance(export, (exports.TotalVolume, exports.AverageVolume)):
                     export.compute(self.dx)

test/test_surface_temperature.py

@@ -0,0 +1,191 @@
+import numpy as np
+import ufl
+from dolfinx import fem
+import pytest
+import os
+
+import festim as F
+
+
+@pytest.mark.parametrize(
+    "T_function, expected_values",
+    [
+        (3, 3),
+        (lambda t: t, 3.0),
+        (lambda x, t: 1.0 + x[0] + t, 10.0),
+    ],
+)
+def test_surface_temperature_compute_1D(T_function, expected_values):
+    """Test that the average surface temperature export computes the correct value."""
+
+    # BUILD
+    L = 6.0
+    my_mesh = F.Mesh1D(np.linspace(0, L, 10000))
+    dummy_surface = F.SurfaceSubdomain1D(id=1, x=6)
+    dummy_volume = F.VolumeSubdomain1D(
+        id=1, borders=[0, L], material=F.Material(D_0=1, E_D=1, name="dummy")
+    )
+    facet_meshtags, temp = my_mesh.define_meshtags(
+        surface_subdomains=[dummy_surface], volume_subdomains=[dummy_volume]
+    )
+
+    ds = ufl.Measure("ds", domain=my_mesh.mesh, subdomain_data=facet_meshtags)
+
+    my_model = F.HydrogenTransportProblem(
+        mesh=my_mesh,
+        temperature=T_function,
+    )
+    my_model.t = fem.Constant(my_model.mesh.mesh, 0.0)
+    dt = fem.Constant(my_mesh.mesh, 1.0)
+
+    my_model.define_temperature()
+    my_model.initialise_exports()
+
+    # RUN
+    for i in range(3):
+        my_model.t.value += dt.value
+        my_model.update_time_dependent_values()
+
+    my_export = F.AverageSurfaceTemperature(surface=dummy_surface)
+    my_export.temperature_field = my_model.temperature_fenics
+    my_export.compute(ds)
+
+    # TEST
+    assert np.isclose(my_export.value, expected_values)
+
+
+def test_title(tmp_path):
+    surf_1 = F.SurfaceSubdomain(id=1)
+    results = "test.csv"
+    temp = 400
+    surface_temp = F.AverageSurfaceTemperature(surface=surf_1, filename=results)
+
+    my_model = F.HydrogenTransportProblem(
+        temperature=temp,
+    )
+    surface_temp.filename = os.path.join(tmp_path, "test.csv")
+    surface_temp.value = 1
+
+    assert surface_temp.title == "Temperature surface 1"
+
+
+@pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"])
+def test_title_generation(tmp_path, value):
+    """Test that the title is made to be written to the header in a csv or txt file"""
+    my_model = F.HydrogenTransportProblem(
+        mesh=F.Mesh1D(np.linspace(0, 6.0, 10000)), temperature=500
+    )
+    my_model.define_temperature()
+
+    my_export = F.AverageSurfaceTemperature(
+        filename=os.path.join(tmp_path, f"{value}"),
+        surface=F.SurfaceSubdomain1D(id=35, x=1),
+    )
+    my_export.value = 2.0
+    my_export.write(0)
+    title = np.genfromtxt(my_export.filename, delimiter=",", max_rows=1, dtype=str)
+
+    expected_title = "Temperature surface 35"
+
+    assert title[1] == expected_title
+
+
+def test_write_overwrite(tmp_path):
+    """Test that the write method overwrites the file if it already exists"""
+    filename = os.path.join(tmp_path, "my_export.csv")
+    my_model = F.HydrogenTransportProblem(
+        mesh=F.Mesh1D(np.linspace(0, 6.0, 10000)), temperature=500
+    )
+    my_model.define_temperature()
+
+    my_export = F.AverageSurfaceTemperature(
+        filename=filename,
+        surface=F.SurfaceSubdomain1D(id=35, x=1),
+    )
+    my_export.value = 2.0
+    my_export.write(0)
+    my_export.write(1)
+
+    my_export2 = F.AverageSurfaceTemperature(
+        filename=filename,
+        surface=F.SurfaceSubdomain1D(id=1, x=1),
+    )
+    my_export2.value = 3.0
+    my_export2.write(1)
+    my_export2.write(2)
+    my_export2.write(3)
+
+    data = np.genfromtxt(filename, delimiter=",", names=True)
+    file_length = data.size
+    expected_length = 3
+
+    assert file_length == expected_length
+
+
+def test_filename_setter_raises_TypeError():
+    """Test that a TypeError is raised when the filename is not a string"""
+
+    with pytest.raises(TypeError, match="filename must be of type str"):
+        my_model = F.HydrogenTransportProblem(
+            mesh=F.Mesh1D(np.linspace(0, 6.0, 10000)), temperature=500
+        )
+        my_model.define_temperature()
+
+        F.AverageSurfaceTemperature(
+            filename=1,
+            surface=F.SurfaceSubdomain1D(id=1, x=1),
+        )
+
+
+def test_filename_setter_raises_ValueError(tmp_path):
+    """Test that a ValueError is raised when the filename does not end with .csv or .txt"""
+
+    with pytest.raises(ValueError):
+        my_model = F.HydrogenTransportProblem(
+            mesh=F.Mesh1D(np.linspace(0, 6.0, 10000)), temperature=500
+        )
+        my_model.define_temperature()
+
+        F.AverageSurfaceTemperature(
+            filename=os.path.join(tmp_path, "my_export.xdmf"),
+            surface=F.SurfaceSubdomain1D(id=1, x=1),
+        )
+
+
+@pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"])
+def test_writer(tmp_path, value):
+    """Test that the writes values at each timestep to either a csv or txt file"""
+    my_model = F.HydrogenTransportProblem(
+        mesh=F.Mesh1D(np.linspace(0, 6.0, 10000)), temperature=500
+    )
+    my_model.define_temperature()
+
+    my_export = F.AverageSurfaceTemperature(
+        filename=os.path.join(tmp_path, f"{value}"),
+        surface=F.SurfaceSubdomain1D(id=1, x=0),
+    )
+    my_export.value = 2.0
+
+    for i in range(10):
+        my_export.write(i)
+        file_length = len(np.genfromtxt(my_export.filename, delimiter=","))
+
+        expected_length = i + 2
+
+        assert file_length == expected_length
+
+
+def test_surface_setter_raises_TypeError():
+    """Test that a TypeError is raised when the surface is not a
+    F.SurfaceSubdomain"""
+
+    with pytest.raises(
+        TypeError, match="surface should be an int or F.SurfaceSubdomain"
+    ):
+        my_model = F.HydrogenTransportProblem(
+            mesh=F.Mesh1D(np.linspace(0, 6.0, 10000)), temperature=500
+        )
+        my_model.define_temperature()
+        F.AverageSurfaceTemperature(
+            surface="1",
+        )