diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 8b3a074de..a74fb89d6 100644 --- a/src/festim/__init__.py +++ b/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 diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 31a287244..513ad2b92 100644 --- a/src/festim/exports/__init__.py +++ b/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", ] diff --git a/src/festim/exports/surface_temperature.py b/src/festim/exports/surface_temperature.py new file mode 100644 index 000000000..e15de834c --- /dev/null +++ b/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) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 133b9fd78..0727d8e33 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/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 diff --git a/test/test_surface_temperature.py b/test/test_surface_temperature.py new file mode 100644 index 000000000..373ea91d3 --- /dev/null +++ b/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()