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 SurfaceTemperature
 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",
+    "SurfaceTemperature"
 ]
 
 from .average_surface import AverageSurface
 from .average_volume import AverageVolume
 from .surface_flux import SurfaceFlux
 from .surface_quantity import SurfaceQuantity
+from .surface_temperature import SurfaceTemperature
 from .total_surface import TotalSurface
 from .total_volume import TotalVolume
 from .volume_quantity import VolumeQuantity

src/festim/exports/surface_temperature.py

@@ -0,0 +1,102 @@
+import csv
+from dolfinx import fem
+import ufl
+import festim as F
+
+class SurfaceTemperature:
+    """Exports the average temperature on a given surface.
+
+    Args:
+        temperature_field (fem.Constant or fem.Function): the temperature field to be computed
+        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, temperature_field, surface, filename: str = None) -> None:
+        self.temperature_field = temperature_field
+        self.surface = surface
+        self.filename = filename
+
+        self.t = []
+        self.data = []
+        self._first_time_export = True
+
+    @property
+    def filename(self):
+        return self._filename
+
+    @filename.setter
+    def filename(self, value):
+        if value is None:
+            self._filename = None
+        elif not isinstance(value, str):
+            raise TypeError("filename must be of type str")
+        elif not value.endswith(".csv") and not value.endswith(".txt"):
+            raise ValueError("filename must end with .csv or .txt")
+        self._filename = value
+
+    @property
+    def surface(self):
+        return self._surface
+
+    @surface.setter
+    def surface(self, value):
+        if not isinstance(value, (int, F.SurfaceSubdomain)) or isinstance(value, bool):
+            raise TypeError("surface should be an int or F.SurfaceSubdomain")
+        self._surface = value
+
+    @property
+    def temperature_field(self):
+        return self._temperature_field
+
+    @temperature_field.setter
+    def temperature_field(self, value):
+        # check that temperature field is float, int, fem.Constant, fem.Function, or fem.Expression
+        if not isinstance(value, (fem.Constant, fem.Function, fem.Expression, int, float)):
+            raise TypeError("field must be of type float, int, fem.Constant, fem.Function, or fem.Expression")
+
+        self._temperature_field = value
+
+    @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)
+
+    def write(self, t):
+        """Writes the time and temperature value to the file.
+
+        Args:
+            t (float): current time value
+        """
+        if self.filename is not None:
+            if self._first_time_export:
+                header = ["t(s)", f"{self.title}"]
+                with open(self.filename, mode="w+", newline="") as file:
+                    writer = csv.writer(file)
+                    writer.writerow(header)
+                self._first_time_export = False
+            with open(self.filename, mode="a", newline="") as file:
+                writer = csv.writer(file)
+                writer.writerow([t, self.value])

src/festim/hydrogen_transport_problem.py

@@ -369,6 +369,10 @@
         a string, find species object in self.species"""
 
         for export in self.exports:
+
+            if isinstance(export, exports.SurfaceTemperature):
+                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,8 +419,11 @@
                 export.D = D
                 export.D_expr = D_expr
 
+            elif isinstance(export, exports.SurfaceTemperature):
+                export.temperature_field = self.temperature_fenics
+
             # reset the data and time for SurfaceQuantity and VolumeQuantity
-            if isinstance(export, (exports.SurfaceQuantity, exports.VolumeQuantity)):
+            if isinstance(export, (exports.SurfaceQuantity, exports.VolumeQuantity, exports.SurfaceTemperature)):
                 export.t = []
                 export.data = []
 
@@ -786,6 +793,8 @@
             if source.temperature_dependent:
                 source.update(t=t)
 
+    surface_temp_processed = False
+
     def post_processing(self):
         """Post processes the model"""
 
@@ -818,6 +827,17 @@
                 # if filename given write export data to file
                 if export.filename is not None:
                     export.write(t=float(self.t))
+
+            elif isinstance(export, exports.SurfaceTemperature):
+                export.compute(self.ds)  # compute surface temp
+
+                export.t.append(float(self.t))  # update export time
+
+                # 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,200 @@
+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.SurfaceTemperature(
+        temperature_field=my_model.temperature_fenics, surface=dummy_surface
+    )
+    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.SurfaceTemperature(
+        temperature_field=temp, 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.SurfaceTemperature(
+        filename=os.path.join(tmp_path, f"{value}"),
+        temperature_field=my_model.temperature_fenics,
+        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.SurfaceTemperature(
+        filename=filename,
+        temperature_field=my_model.temperature_fenics,
+        surface=F.SurfaceSubdomain1D(id=35, x=1),
+    )
+    my_export.value = 2.0
+    my_export.write(0)
+    my_export.write(1)
+
+    my_export2 = F.SurfaceTemperature(
+        filename=filename,
+        temperature_field=my_model.temperature_fenics,
+        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.SurfaceTemperature(
+            filename=1,
+            temperature_field=my_model.temperature_fenics,
+            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.SurfaceTemperature(
+            filename=os.path.join(tmp_path, "my_export.xdmf"),
+            temperature_field=my_model.temperature_fenics,
+            surface=F.SurfaceSubdomain1D(id=1, x=1),
+        )
+
+
+def test_field_setter_raises_TypeError():
+    """Test that a TypeError is raised when the field is not an int, float, fem.Constant, fem.Expression, or fem.Function"""
+
+    with pytest.raises(TypeError):
+
+        F.SurfaceTemperature(
+            temperature_field="str",
+            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.SurfaceTemperature(
+        filename=os.path.join(tmp_path, f"{value}"),
+        temperature_field=my_model.temperature_fenics,
+        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.SurfaceTemperature(
+            temperature_field=my_model.temperature_fenics,
+            surface="1",
+        )