Merge remote-tracking branch 'jhdark/surface_flux_export' into reacti…

…ons-in-formulation

festim/exports/surface_flux.py

@@ -1,7 +1,6 @@
 from dolfinx import fem
 import festim as F
 import ufl
-import csv
 
 
 class SurfaceFlux(F.SurfaceQuantity):
@@ -16,23 +15,23 @@ class SurfaceFlux(F.SurfaceQuantity):
         field (festim.Species): species for which the surface flux is computed
         surface (festim.SurfaceSubdomain1D): surface subdomain
         filename (str): name of the file to which the surface flux is exported
-        t (list): list of time values
-        data (list): list of surface flux values
-        title (str): title of the exported data
     """
 
     def __init__(
         self,
-        field,
+        field: F.Species,
         surface: F.SurfaceSubdomain1D,
         filename: str = None,
     ) -> None:
         super().__init__(field, surface, filename)
 
-        self.t = []
-        self.data = []
-
     def compute(self, n, ds):
+        """Computes the value of the surface flux at the surface
+
+        Args:
+            n (ufl.geometry.FacetNormal): normal vector to the surface
+            ds (ufl.Measure): surface measure of the model
+        """
         self.value = fem.assemble_scalar(
             fem.form(
                 -self.D
@@ -41,11 +40,3 @@ def compute(self, n, ds):
             )
         )
         self.data.append(self.value)
-
-    def initialise_export(self):
-        self.title = "Flux surface {}: {}".format(self.surface.id, self.field.name)
-
-        if self.filename is not None:
-            with open(self.filename, mode="w", newline="") as file:
-                writer = csv.writer(file)
-                writer.writerow(["t(s)", f"{self.title}"])

festim/exports/surface_quantity.py

@@ -1,5 +1,6 @@
 import festim as F
 import csv
+import os
 
 
 class SurfaceQuantity:
@@ -14,13 +15,18 @@ class SurfaceQuantity:
         field (festim.Species): species for which the surface flux is computed
         surface (festim.SurfaceSubdomain1D): surface subdomain
         filename (str): name of the file to which the surface flux is exported
+        t (list): list of time values
+        data (list): list of values of the surface quantity
     """
 
     def __init__(self, field, surface, filename: str = None) -> None:
         self.field = field
         self.surface = surface
         self.filename = filename
 
+        self.t = []
+        self.data = []
+
     @property
     def filename(self):
         return self._filename
@@ -31,8 +37,8 @@ def filename(self, value):
             self._filename = None
         elif not isinstance(value, str):
             raise TypeError("filename must be of type str")
-        elif not value.endswith(".csv"):
-            raise ValueError("filename must end with .csv")
+        elif not value.endswith(".csv") and not value.endswith(".txt"):
+            raise ValueError("filename must end with .csv or .txt")
         self._filename = value
 
     @property
@@ -60,6 +66,17 @@ def field(self, value):
         self._field = value
 
     def write(self, t):
+        """If the filename doesnt exist yet, create it and write the header,
+        then append the time and value to the file"""
+
+        if not os.path.isfile(self.filename):
+            title = "Flux surface {}: {}".format(self.surface.id, self.field.name)
+
+            if self.filename is not None:
+                with open(self.filename, mode="w", newline="") as file:
+                    writer = csv.writer(file)
+                    writer.writerow(["t(s)", f"{title}"])
+
         with open(self.filename, mode="a", newline="") as file:
             writer = csv.writer(file)
             writer.writerow([t, self.value])

festim/hydrogen_transport_problem.py

@@ -271,8 +271,6 @@ def initialise_exports(self):
                 export.define_writer(MPI.COMM_WORLD)
                 if isinstance(export, F.XDMFExport):
                     export.writer.write_mesh(self.mesh.mesh)
-            elif isinstance(export, F.SurfaceQuantity):
-                export.initialise_export()
 
         # compute diffusivity function for surface fluxes
 

festim/material.py

@@ -59,7 +59,7 @@ def get_D_0(self, species=None):
                 raise ValueError(f"{species} is not in D_0 keys")
 
         else:
-            raise ValueError("D_0 must be either a float or a dict")
+            raise TypeError("D_0 must be either a float, int or a dict")
 
     def get_E_D(self, species=None):
         """Returns the activation energy of the diffusion coefficient
@@ -85,7 +85,7 @@ def get_E_D(self, species=None):
                 raise ValueError(f"{species} is not in E_D keys")
 
         else:
-            raise ValueError("E_D must be either a float or a dict")
+            raise TypeError("E_D must be either a float, int or a dict")
 
     def get_diffusion_coefficient(self, mesh, temperature, species=None):
         """Defines the diffusion coefficient

test/test_h_transport_problem.py

@@ -229,6 +229,8 @@ def test_initialise_exports_find_species_with_one_field():
 
 
 def test_define_D_global_different_temperatures():
+    """Test that the D_global object is correctly defined when the temperature
+    is different in the volume subdomains"""
     D_0, E_D = 1.5, 0.1
     my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat")
     surf = F.SurfaceSubdomain1D(id=1, x=0)
@@ -242,12 +244,6 @@ def test_define_D_global_different_temperatures():
         ],
         species=[H],
         temperature=lambda x: 100.0 * x[0] + 50,
-        exports=[
-            F.SurfaceFlux(
-                field=H,
-                surface=surf,
-            ),
-        ],
     )
 
     my_model.define_function_spaces()
@@ -267,11 +263,12 @@ def test_define_D_global_different_temperatures():
 
 
 def test_define_D_global_different_materials():
+    """Test that the D_global object is correctly defined when the material
+    is different in the volume subdomains"""
     D_0_left, E_D_left = 1.0, 0.1
     D_0_right, E_D_right = 2.0, 0.2
     my_mat_L = F.Material(D_0=D_0_left, E_D=E_D_left, name="my_mat_L")
     my_mat_R = F.Material(D_0=D_0_right, E_D=E_D_right, name="my_mat_R")
-    surf = F.SurfaceSubdomain1D(id=1, x=0)
     H = F.Species("H")
 
     my_model = F.HydrogenTransportProblem(
@@ -280,14 +277,8 @@ def test_define_D_global_different_materials():
             F.VolumeSubdomain1D(id=1, borders=[0, 2], material=my_mat_L),
             F.VolumeSubdomain1D(id=2, borders=[2, 4], material=my_mat_R),
         ],
-        species=[F.Species("H"), F.Species("D")],
+        species=[H],
         temperature=500,
-        exports=[
-            F.SurfaceFlux(
-                field=H,
-                surface=surf,
-            ),
-        ],
     )
 
     my_model.define_function_spaces()
@@ -298,15 +289,19 @@ def test_define_D_global_different_materials():
 
     computed_values = [D_computed.x.array[0], D_computed.x.array[-1]]
 
-    D_analytical_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature))
-    D_analytical_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature))
+    D_expected_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature))
+    D_expected_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature))
 
-    expected_values = [D_analytical_left, D_analytical_right]
+    expected_values = [D_expected_left, D_expected_right]
 
     assert np.isclose(computed_values, expected_values).all()
 
 
 def test_initialise_exports_multiple_exports_same_species():
+    """Test that the diffusion coefficient within the D_global object function is the same
+    for multiple exports of the same species, and that D_global object is only
+    created once per species"""
+
     D_0, E_D = 1.5, 0.1
     my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat")
     surf_1 = F.SurfaceSubdomain1D(id=1, x=0)
@@ -338,8 +333,93 @@ def test_initialise_exports_multiple_exports_same_species():
     my_model.define_temperature()
     my_model.initialise_exports()
 
-    Ds = []
-    for export in my_model.exports:
-        Ds.append(export.D)
+    Ds = [export.D for export in my_model.exports]
 
-    assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0])
+    assert Ds[0].x.array[0] == Ds[1].x.array[0]
+
+
+def test_define_D_global_multispecies():
+    """Test that the D_global object is correctly defined when there are multiple
+    species in one subdomain"""
+    A = F.Species("A")
+    B = F.Species("B")
+
+    D_0_A, D_0_B = 1.0, 2.0
+    E_D_A, E_D_B = 0.1, 0.2
+
+    my_mat = F.Material(
+        D_0={A: D_0_A, B: D_0_B}, E_D={A: E_D_A, B: E_D_B}, name="my_mat"
+    )
+    surf = F.SurfaceSubdomain1D(id=1, x=1)
+
+    my_model = F.HydrogenTransportProblem(
+        mesh=F.Mesh1D(np.linspace(0, 1, num=101)),
+        subdomains=[
+            F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat),
+        ],
+        species=[F.Species("A"), F.Species("B")],
+        temperature=500,
+    )
+
+    my_model.define_function_spaces()
+    my_model.define_markers_and_measures()
+    my_model.define_temperature()
+
+    D_A_computed, D_A_expr = my_model.define_D_global(A)
+    D_B_computed, D_B_expr = my_model.define_D_global(B)
+
+    computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]]
+
+    D_expected_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature))
+    D_expected_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature))
+
+    expected_values = [D_expected_A, D_expected_B]
+
+    assert np.isclose(computed_values, expected_values).all()
+
+
+def test_post_processing_update_D_global():
+    """Test that the D_global object is updated at each time
+    step when temperture is time dependent"""
+    my_mesh = F.Mesh1D(np.linspace(0, 1, num=11))
+    my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat")
+    surf = F.SurfaceSubdomain1D(id=1, x=1)
+
+    # create species and interpolate solution
+    H = F.Species("H")
+    V = fem.FunctionSpace(my_mesh.mesh, ("CG", 1))
+    u = fem.Function(V)
+    u.interpolate(lambda x: 2 * x[0] ** 2 + 1)
+    H.solution = u
+
+    my_export = F.SurfaceFlux(
+        field=H,
+        surface=surf,
+    )
+
+    # Build the model
+    my_model = F.HydrogenTransportProblem(
+        mesh=my_mesh,
+        subdomains=[F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat), surf],
+        species=[H],
+        temperature=lambda t: 500 * t,
+        exports=[my_export],
+    )
+
+    my_model.define_function_spaces()
+    my_model.define_markers_and_measures()
+    my_model.t = fem.Constant(my_model.mesh.mesh, 1.0)
+    my_model.define_temperature()
+    my_model.initialise_exports()
+
+    # RUN
+    my_model.post_processing()
+    value_t_1 = my_export.D.x.array[-1]
+
+    my_model.t = fem.Constant(my_model.mesh.mesh, 2.0)
+    my_model.update_time_dependent_values()
+    my_model.post_processing()
+    value_t_2 = my_export.D.x.array[-1]
+
+    # TEST
+    assert value_t_1 != value_t_2

test/test_material.py

@@ -98,6 +98,7 @@ def test_multispecies_dict_different_keys():
 def test_D_0_type_raises_error():
     """Test that a value error is raised in the get_diffusion_coefficient
     function"""
+    # TODO remove this when material class is updated
     A = F.Species("A")
     my_mat = F.Material(D_0=[1, 1], E_D=0.1)
 
@@ -128,3 +129,71 @@ def test_error_raised_when_species_not_not_in_D_0_dict():
 
     with pytest.raises(ValueError, match="J is not in D_0 keys"):
         my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=J)
+
+
+def test_D_0_raises_ValueError_if_species_not_provided_in_dict():
+    """Test that a value error is raised in the get_diffusion_coefficient
+    function"""
+    # TODO remove this when material class is updated
+    A = F.Species("A")
+    B = F.Species("B")
+    my_mat = F.Material(D_0={A: 1, B: 2}, E_D=1)
+
+    with pytest.raises(ValueError, match="species must be provided if D_0 is a dict"):
+        my_mat.get_D_0()
+
+
+def test_D_0_raises_ValueError_if_species_given_not_in_dict_keys():
+    """Test that a value error is raised in the get_diffusion_coefficient
+    function"""
+    # TODO remove this when material class is updated
+    A = F.Species("A")
+    B = F.Species("B")
+    J = F.Species("J")
+    my_mat = F.Material(D_0={A: 1, B: 2}, E_D=1)
+
+    with pytest.raises(ValueError, match="J is not in D_0 keys"):
+        my_mat.get_D_0(species=J)
+
+
+def test_raises_TypeError_when_D_0_is_not_correct_type():
+    """Test that a TypeError is raised when D_0 is not a float or a dict"""
+
+    my_mat = F.Material(D_0=[1, 2], E_D=1)
+
+    with pytest.raises(TypeError, match="D_0 must be either a float, int or a dict"):
+        my_mat.get_D_0()
+
+
+def test_E_D_raises_ValueError_if_species_not_provided_in_dict():
+    """Test that a value error is raised in the get_diffusion_coefficient
+    function"""
+    # TODO remove this when material class is updated
+    A = F.Species("A")
+    B = F.Species("B")
+    my_mat = F.Material(D_0=1, E_D={A: 1, B: 2})
+
+    with pytest.raises(ValueError, match="species must be provided if E_D is a dict"):
+        my_mat.get_E_D()
+
+
+def test_E_D_raises_ValueError_if_species_given_not_in_dict_keys():
+    """Test that a value error is raised in the get_diffusion_coefficient
+    function"""
+    # TODO remove this when material class is updated
+    A = F.Species("A")
+    B = F.Species("B")
+    J = F.Species("J")
+    my_mat = F.Material(D_0=1, E_D={A: 1, B: 2})
+
+    with pytest.raises(ValueError, match="J is not in E_D keys"):
+        my_mat.get_E_D(species=J)
+
+
+def test_raises_TypeError_when_E_D_is_not_correct_type():
+    """Test that a TypeError is raised when E_D is not a float or a dict"""
+
+    my_mat = F.Material(D_0=1, E_D=[1, 2])
+
+    with pytest.raises(TypeError, match="E_D must be either a float, int or a dict"):
+        my_mat.get_E_D()

test/test_permeation_problem.py

@@ -62,6 +62,7 @@ def test_permeation_problem(mesh_size=1001):
         ),
     ]
     outgassing_flux = F.SurfaceFlux(
+        filename="outgassing_flux.txt",
         field=mobile_H,
         surface=right_surface,
     )