Export average surface temperature

src/festim/__init__.py

@@ -32,6 +32,7 @@
 from .exports.profile_1d import Profile1DExport
 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

@@ -7,6 +7,7 @@
 from .profile_1d import Profile1DExport
 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
@@ -30,4 +31,5 @@
     "VTXTemperatureExport",
     "VolumeQuantity",
     "XDMFExport",
+    "AverageSurfaceTemperature",
 ]

src/festim/exports/surface_temperature.py

@@ -0,0 +1,55 @@
+import csv
+from dolfinx import fem
+import ufl
+from .surface_quantity import SurfaceQuantity
+from festim.subdomain.surface_subdomain import SurfaceSubdomain
+
+
+class AverageSurfaceTemperature(SurfaceQuantity):
+    """Exports the average temperature on a given surface.
+
+    Args:
+        surface: the surface subdomain
+        filename: name of the file to which the average surface temperature is exported
+
+    Attributes:
+        temperature_field: 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
+    """
+
+    surface: int | SurfaceSubdomain
+    filename: str | None
+
+    temperature_field: fem.Constant | fem.Function
+
+    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
+
+        self.value = fem.assemble_scalar(
+            fem.form(temperature_field * ds(self.surface.id))
+        ) / fem.assemble_scalar(
+            fem.form(1 * ds(self.surface.id))
+        )  # integral over surface / surface area
+
+        self.data.append(self.value)

src/festim/hydrogen_transport_problem.py

@@ -411,6 +411,8 @@ def initialise_exports(self):
         a string, find species object in self.species"""
 
         for export in self.exports:
+            if isinstance(export, exports.AverageSurfaceTemperature):
+                continue
             if isinstance(export, exports.ExportBaseClass):
                 if export.times:
                     for time in export.times:
@@ -478,6 +480,8 @@ def initialise_exports(self):
 
         for export in self.exports:
             if isinstance(export, exports.SurfaceQuantity):
+                if isinstance(export, exports.AverageSurfaceTemperature):
+                    continue
                 if export.field in spe_to_D_global:
                     # if already computed then use the same D
                     D = spe_to_D_global[export.field]
@@ -491,6 +495,10 @@ def initialise_exports(self):
                 # add the global D to the export
                 export.D = D
                 export.D_expr = D_expr
+
+            if isinstance(export, exports.AverageSurfaceTemperature):
+                export.temperature_field = self.temperature_fenics
+
             if isinstance(
                 export,
                 exports.MaximumVolume
@@ -1026,6 +1034,8 @@ def post_processing(self):
                             "evaluation of surface flux values"
                         )
                     export.compute(export.field.solution, self.ds)
+                elif isinstance(export, exports.AverageSurfaceTemperature):
+                    export.compute(self.ds)
                 else:
                     export.compute()
                 # update export data
@@ -1034,6 +1044,7 @@ def post_processing(self):
                 # 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(u=export.field.solution, dx=self.dx)
@@ -1708,6 +1719,11 @@ def initialise_exports(self):
         # HydrogenTransportProblem
         for export in self.exports:
             if isinstance(export, exports.SurfaceQuantity):
+                if isinstance(export, exports.AverageSurfaceTemperature):
+                    raise NotImplementedError(
+                        f"Export type {type(export)} not implemented for "
+                        f"mixed-domain approach"
+                    )
                 volume = self.surface_to_volume[export.surface]
                 D = volume.material.get_diffusion_coefficient(
                     self.mesh.mesh, self.temperature_fenics, export.field

test/test_surface_temperature.py

@@ -0,0 +1,127 @@
+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=L)
+    dummy_volume = F.VolumeSubdomain1D(
+        id=1, borders=(0.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)
+
+    dt = F.Stepsize(initial_value=1)
+    settings = F.Settings(atol=1e05, rtol=1e-10, stepsize=dt, final_time=10)
+    my_model = F.HydrogenTransportProblem(
+        mesh=my_mesh, temperature=T_function, settings=settings
+    )
+
+    my_model.species = [F.Species("H")]
+    my_model.subdomains = [
+        dummy_surface,
+        dummy_volume,
+    ]
+
+    my_model.t = fem.Constant(my_mesh.mesh, 0.0)
+    dt = fem.Constant(my_mesh.mesh, 1.0)
+
+    my_model.define_temperature()
+
+    my_export = F.AverageSurfaceTemperature(surface=dummy_surface)
+    my_export.temperature_field = my_model.temperature_fenics
+    my_model.exports = [my_export]
+
+    # RUN
+    for _ in range(3):
+        my_model.t.value += dt.value
+        my_model.update_time_dependent_values()
+
+    my_model.initialise()
+    my_model.run()
+
+    # 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_not_implemented_error_raised_with_initialise_exports_multiple_subdomains():
+    """Test that NotImplementedError is raised for problems with multiple volume domains."""
+
+    # BUILD
+    L = 1.0
+    test_mesh = F.Mesh1D(vertices=np.linspace(0, L, num=101))
+
+    my_mat = F.Material(D_0=1, E_D=1)
+    dummy_surface = F.SurfaceSubdomain1D(id=1, x=L)
+    vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=my_mat)
+    vol_2 = F.VolumeSubdomain1D(id=1, borders=[0.5, L], material=my_mat)
+
+    H = F.Species("H")
+    my_model = F.HydrogenTransportProblemDiscontinuous(
+        mesh=test_mesh,
+        temperature=10,
+        subdomains=[dummy_surface, vol_1, vol_2],
+        species=[F.Species("H")],
+    )
+
+    my_model.exports = [F.AverageSurfaceTemperature(surface=dummy_surface)]
+
+    # TEST
+    with pytest.raises(NotImplementedError):
+        my_model.initialise_exports()