Merge pull request #1037 from jhdark/diffusion_fix

Only get diffusion coeff if mobile

Author: jhdark

src/festim/hydrogen_transport_problem.py

@@ -814,10 +814,10 @@ def create_formulation(self):
             v = spe.test_function
 
             for vol in self.volume_subdomains:
-                D = vol.material.get_diffusion_coefficient(
-                    self.mesh.mesh, self.temperature_fenics, spe
-                )
                 if spe.mobile:
+                    D = vol.material.get_diffusion_coefficient(
+                        self.mesh.mesh, self.temperature_fenics, spe
+                    )
                     match self.mesh.coordinate_system:
                         case CoordinateSystem.CARTESIAN:
                             self.formulation += ufl.dot(
@@ -1379,13 +1379,13 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
             u_n = spe.subdomain_to_prev_solution[subdomain]
             v = spe.subdomain_to_test_function[subdomain]
 
-            D = subdomain.material.get_diffusion_coefficient(
-                self.mesh.mesh, self.temperature_fenics, spe
-            )
             if self.settings.transient:
                 form += ((u - u_n) / self.dt) * v * self.dx(subdomain.id)
 
             if spe.mobile:
+                D = subdomain.material.get_diffusion_coefficient(
+                    self.mesh.mesh, self.temperature_fenics, spe
+                )
                 match self.mesh.coordinate_system:
                     case CoordinateSystem.CARTESIAN:
                         form += ufl.dot(D * ufl.grad(u), ufl.grad(v)) * self.dx(

test/system_tests/test_misc.py

@@ -245,3 +245,42 @@ def test_del():
     my_model.exports = [F.VTXSpeciesExport(filename="h.bp", field=H)]
 
     my_model.__del__()
+
+
+def test_multispecies_with_immobile():
+    """test to catch bug 1036, that the diffusion coefficent for a non-mobile species
+    is not requested in the dict given to material"""
+
+    my_model = F.HydrogenTransportProblem()
+    my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=50))
+
+    H = F.Species("H")
+    D = F.Species("D")
+    trap = F.Species("trap", mobile=False)
+    trapping_site = F.ImplicitSpecies(n=100, others=[H])
+    my_model.species = [H, trap]
+
+    my_mat = F.Material(
+        D_0={H: 1, D: 1},
+        E_D={H: 0, D: 0},
+    )
+    vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat)
+    my_model.subdomains = [vol]
+
+    my_model.reactions = [
+        F.Reaction(
+            reactant=[H, trapping_site],
+            k_0=1,
+            E_k=0,
+            p_0=1,
+            E_p=0,
+            product=[trap],
+            volume=vol,
+        )
+    ]
+
+    my_model.temperature = 500
+
+    my_model.settings = F.Settings(transient=False, atol=1e-10, rtol=1e-10)
+
+    my_model.initialise()