Merge branch 'fenicsx' into milestones_from_exports

Author: jhdark

examples/multi_material_2d_bis.py

@@ -129,7 +129,7 @@ def half(x):
 
 my_model.temperature = lambda x: 400 + 10 * x[1] + 100 * x[0]
 
-my_model.settings = F.Settings(atol=None, rtol=1e-5, transient=False)
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-5, transient=False)
 
 my_model.exports = [
     F.VTXSpeciesExport(f"c_{subdomain.id}.bp", field=H, subdomain=subdomain)
@@ -149,6 +149,7 @@ def half(x):
 entity_maps = {sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains}
 entity_maps[mesh] = bottom_domain.submesh_to_mesh
 
+ds = ufl.Measure("ds", domain=mesh, subdomain_data=my_model.facet_meshtags)
 ds_b = ufl.Measure("ds", domain=bottom_domain.submesh, subdomain_data=bottom_domain.ft)
 ds_t = ufl.Measure("ds", domain=top_domain.submesh, subdomain_data=top_domain.ft)
 dx_b = ufl.Measure("dx", domain=bottom_domain.submesh)
@@ -191,3 +192,16 @@ def half(x):
     entity_maps={mesh: top_domain.submesh_to_mesh},
 )
 print(dolfinx.fem.assemble_scalar(form))
+
+
+# using submesh function
+u = H.subdomain_to_post_processing_solution[bottom_domain]
+form = dolfinx.fem.form(
+    ufl.dot(
+        D * ufl.grad(u),
+        n_b,
+    )
+    * ds(bottom_surface.id),
+    entity_maps={sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains},
+)
+print(dolfinx.fem.assemble_scalar(form))

src/festim/exports/average_volume.py

@@ -19,12 +19,12 @@ class AverageVolume(VolumeQuantity):
     def title(self):
         return f"Average {self.field.name} volume {self.volume.id}"
 
-    def compute(self, dx):
+    def compute(self, u, dx, entity_maps=None):
         """
         Computes the average value of solution function within the defined volume
         subdomain, and appends it to the data list
         """
         self.value = fem.assemble_scalar(
-            fem.form(self.field.solution * dx(self.volume.id))
+            fem.form(u * dx(self.volume.id), entity_maps=entity_maps)
         ) / fem.assemble_scalar(fem.form(1 * dx(self.volume.id)))
         self.data.append(self.value)

src/festim/exports/surface_flux.py

@@ -20,26 +20,27 @@ class SurfaceFlux(SurfaceQuantity):
     def title(self):
         return f"{self.field.name} flux surface {self.surface.id}"
 
-    def compute(self, ds):
+    def compute(self, u, ds: ufl.Measure, entity_maps=None):
         """Computes the value of the flux at the surface
 
         Args:
-            ds (ufl.Measure): surface measure of the model
+            u: field for which the flux is computed
+            ds: surface measure of the model
+            entity_maps: entity maps relating parent mesh and submesh
         """
 
         # obtain mesh normal from field
         # if case multispecies, solution is an index, use sub_function_space
-        if isinstance(self.field.solution, ufl.indexed.Indexed):
+        if isinstance(u, ufl.indexed.Indexed):
             mesh = self.field.sub_function_space.mesh
         else:
-            mesh = self.field.solution.function_space.mesh
+            mesh = u.function_space.mesh
         n = ufl.FacetNormal(mesh)
 
         self.value = fem.assemble_scalar(
             fem.form(
-                -self.D
-                * ufl.dot(ufl.grad(self.field.solution), n)
-                * ds(self.surface.id)
+                -self.D * ufl.dot(ufl.grad(u), n) * ds(self.surface.id),
+                entity_maps=entity_maps,
             )
         )
         self.data.append(self.value)

src/festim/exports/total_volume.py

@@ -20,15 +20,17 @@ class TotalVolume(VolumeQuantity):
     def title(self):
         return f"Total {self.field.name} volume {self.volume.id}"
 
-    def compute(self, dx: ufl.Measure):
+    def compute(self, u, dx: ufl.Measure, entity_maps=None):
         """
         Computes the value of the total volume of the field in the volume subdomain
         and appends it to the data list
 
         Args:
-            dx (ufl.Measure): volume measure of the model
+            u: field for which the total volume is computed
+            dx: volume measure of the model
+            entity_maps: entity maps relating parent mesh and submesh
         """
         self.value = fem.assemble_scalar(
-            fem.form(self.field.solution * dx(self.volume.id))
+            fem.form(u * dx(self.volume.id), entity_maps=entity_maps)
         )
         self.data.append(self.value)

src/festim/hydrogen_transport_problem.py

@@ -947,9 +947,7 @@ def post_processing(self):
                             "Advection terms are not currently accounted for in the "
                             "evaluation of surface flux values"
                         )
-                    export.compute(
-                        self.ds,
-                    )
+                    export.compute(export.field.solution, self.ds)
                 else:
                     export.compute()
                 # update export data
@@ -960,7 +958,7 @@ def post_processing(self):
                     export.write(t=float(self.t))
             elif isinstance(export, exports.VolumeQuantity):
                 if isinstance(export, exports.TotalVolume | exports.AverageVolume):
-                    export.compute(self.dx)
+                    export.compute(u=export.field.solution, dx=self.dx)
                 else:
                     export.compute()
                 # update export data
@@ -1064,6 +1062,21 @@ def initialise(self):
 
         for subdomain in self.volume_subdomains:
             self.define_function_spaces(subdomain)
+        # create global DG function spaces of degree 0 and 1
+        element_DG0 = basix.ufl.element(
+            "DG",
+            self.mesh.mesh.basix_cell(),
+            0,
+            basix.LagrangeVariant.equispaced,
+        )
+        element_DG1 = basix.ufl.element(
+            "DG",
+            self.mesh.mesh.basix_cell(),
+            1,
+            basix.LagrangeVariant.equispaced,
+        )
+        self.V_DG_0 = fem.functionspace(self.mesh.mesh, element_DG0)
+        self.V_DG_1 = fem.functionspace(self.mesh.mesh, element_DG1)
 
         self.define_temperature()
         self.convert_source_input_values_to_fenics_objects()
@@ -1503,8 +1516,32 @@ def initialise_exports(self):
                         f"Export type {type(export)} not implemented for "
                         f"mixed-domain approach"
                     )
-            else:
-                raise NotImplementedError(f"Export type {type(export)} not implemented")
+
+        # compute diffusivity function for surface fluxes
+
+        spe_to_D_global = {}  # links species to global D function
+        spe_to_D_global_expr = {}  # links species to D expression
+
+        for export in self.exports:
+            if isinstance(export, exports.SurfaceQuantity):
+                if export.field in spe_to_D_global:
+                    # if already computed then use the same D
+                    D = spe_to_D_global[export.field]
+                    D_expr = spe_to_D_global_expr[export.field]
+                else:
+                    # compute D and add it to the dict
+                    D, D_expr = self.define_D_global(export.field)
+                    spe_to_D_global[export.field] = D
+                    spe_to_D_global_expr[export.field] = D_expr
+
+                # add the global D to the export
+                export.D = D
+                export.D_expr = D_expr
+
+            # reset the data and time for SurfaceQuantity and VolumeQuantity
+            if isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity):
+                export.t = []
+                export.data = []
 
     def post_processing(self):
         # update post-processing solutions (for each species in each subdomain)
@@ -1536,6 +1573,56 @@ def post_processing(self):
                 if export.is_it_time_to_export(float(self.t)):
                     export.writer.write(float(self.t))
 
+            # handle derived quantities
+            if isinstance(export, exports.SurfaceQuantity):
+                if isinstance(
+                    export,
+                    exports.SurfaceFlux | exports.TotalSurface | exports.AverageSurface,
+                ):
+                    if len(self.advection_terms) > 0:
+                        warnings.warn(
+                            "Advection terms are not currently accounted for in the "
+                            "evaluation of surface flux values"
+                        )
+                    export_surf = export.surface
+                    export_vol = self.surface_to_volume[export_surf]
+                    submesh_function = (
+                        export.field.subdomain_to_post_processing_solution[export_vol]
+                    )
+                    export.compute(
+                        u=submesh_function,
+                        ds=self.ds,
+                        entity_maps={
+                            sd.submesh: sd.parent_to_submesh
+                            for sd in self.volume_subdomains
+                        },
+                    )
+                else:
+                    export.compute()
+
+            elif isinstance(export, exports.VolumeQuantity):
+                if isinstance(export, exports.TotalVolume | exports.AverageVolume):
+                    export.compute(
+                        u=export.field.subdomain_to_post_processing_solution[
+                            export_vol
+                        ],
+                        dx=self.dx,
+                        entity_maps={
+                            sd.submesh: sd.parent_to_submesh
+                            for sd in self.volume_subdomains
+                        },
+                    )
+                else:
+                    export.compute()
+
+            if isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity):
+                # update export data
+                export.t.append(float(self.t))
+
+                # if filename given write export data to file
+                if export.filename is not None:
+                    export.write(t=float(self.t))
+
     def iterate(self):
         """Iterates the model for a given time step"""
         if self.show_progress_bar:

test/system_tests/test_multi_mat_penalty.py

@@ -161,3 +161,74 @@ def c_exact_bot_np(x):
 
     assert L2_error_top < 1e-3
     assert L2_error_bot < 1e-3
+
+
+def test_derived_quantities_multi_mat():
+    K_S_top = 1.0
+    K_S_bot = 1.0
+    D_top = 2.0
+    D_bot = 1.0
+
+    mesh, mt, ct = generate_mesh(20)
+
+    my_model = F.HydrogenTransportProblemDiscontinuous()
+    my_model.mesh = F.Mesh(mesh)
+    my_model.volume_meshtags = ct
+    my_model.facet_meshtags = mt
+
+    material_top = F.Material(D_0=D_top, E_D=0, K_S_0=K_S_top, E_K_S=0)
+    material_bottom = F.Material(D_0=D_bot, E_D=0, K_S_0=K_S_bot, E_K_S=0)
+
+    material_top.solubility_law = "sievert"
+    material_bottom.solubility_law = "sievert"
+
+    top_domain = F.VolumeSubdomain(4, material=material_top)
+    bottom_domain = F.VolumeSubdomain(3, material=material_bottom)
+
+    top_surface = F.SurfaceSubdomain(id=1)
+    bottom_surface = F.SurfaceSubdomain(id=2)
+    my_model.subdomains = [
+        bottom_domain,
+        top_domain,
+        top_surface,
+        bottom_surface,
+    ]
+
+    my_model.interfaces = [
+        F.Interface(5, (bottom_domain, top_domain), penalty_term=1),
+    ]
+    my_model.surface_to_volume = {
+        top_surface: top_domain,
+        bottom_surface: bottom_domain,
+    }
+
+    H = F.Species("H", mobile=True)
+
+    my_model.species = [H]
+
+    for species in my_model.species:
+        species.subdomains = [bottom_domain, top_domain]
+
+    my_model.boundary_conditions = [
+        F.FixedConcentrationBC(top_surface, value=1, species=H),
+        F.FixedConcentrationBC(bottom_surface, value=0, species=H),
+    ]
+
+    my_model.temperature = 500.0
+
+    my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+    my_model.exports = [
+        F.SurfaceFlux(field=H, surface=top_surface),
+        F.SurfaceFlux(field=H, surface=bottom_surface),
+        F.AverageVolume(field=H, volume=bottom_domain),
+        F.AverageVolume(field=H, volume=top_domain),
+        F.TotalVolume(field=H, volume=bottom_domain),
+        F.TotalVolume(field=H, volume=top_domain),
+    ]
+
+    my_model.initialise()
+    my_model.run()
+
+    print("Top surface flux:", my_model.exports[0].data)
+    print("Bottom surface flux:", my_model.exports[1].data)
+    assert np.isclose(my_model.exports[0].data[0], -my_model.exports[1].data[0])

test/system_tests/test_surface_reactions.py

@@ -14,7 +14,7 @@ def __init__(self, reaction: F.SurfaceReactionBC):
         )
         self.reaction = reaction.flux_bcs[0]
 
-    def compute(self, ds):
+    def compute(self, u, ds, entity_maps=None):
         self.value = fem.assemble_scalar(
             fem.form(self.reaction.value_fenics * ds(self.surface.id))
         )

test/test_average_volume.py

@@ -30,7 +30,7 @@ def test_average_volume_compute_1D():
     my_export = F.AverageVolume(field=my_species, volume=dummy_volume)
 
     # RUN
-    my_export.compute(dx)
+    my_export.compute(my_species.solution, dx)
 
     # TEST
     expected_value = 3.5

test/test_surface_quantity.py

@@ -45,7 +45,7 @@ def test_surface_flux_export_compute():
     my_export.D = D
 
     # RUN
-    my_export.compute(ds=ds)
+    my_export.compute(my_species.solution, ds=ds)
 
     # TEST
     # flux = -D grad(c)_ \cdot n = -D dc/dx = -D * 4 * x

test/test_total_volume.py

@@ -36,7 +36,7 @@ def test_total_volume_export_compute():
     )
 
     # RUN
-    my_export.compute(dx=dx)
+    my_export.compute(my_species.solution, dx=dx)
 
     # TEST
     # total = int[0,L] c dx = 2/3 * L^3 + L