Fix SurfaceFlux with Soret

festim/exports/derived_quantities/derived_quantities.py

@@ -154,7 +154,6 @@ def assign_properties_to_quantities(self, materials):
             quantity.Q = materials.Q
 
     def compute(self, t):
-        # TODO need to support for soret flag in surface flux
         row = [t]
         for quantity in self:
             if isinstance(quantity, (MaximumVolume, MinimumVolume)):

festim/exports/derived_quantities/surface_flux.py

@@ -35,6 +35,7 @@ class SurfaceFlux(SurfaceQuantity):
 
     def __init__(self, field, surface) -> None:
         super().__init__(field=field, surface=surface)
+        self.soret = None
 
     @property
     def export_unit(self):
@@ -70,11 +71,11 @@ def prop(self):
         }
         return field_to_prop[self.field]
 
-    def compute(self, soret=False):
+    def compute(self):
         flux = f.assemble(
             self.prop * f.dot(f.grad(self.function), self.n) * self.ds(self.surface)
         )
-        if soret and self.field in [0, "0", "solute"]:
+        if self.soret and self.field in [0, "0", "solute"]:
             flux += f.assemble(
                 self.prop
                 * self.function
@@ -136,7 +137,7 @@ def azimuth_range(self, value):
             raise ValueError("Azimuthal range must be between 0 and pi")
         self._azimuth_range = value
 
-    def compute(self, soret=False):
+    def compute(self):
 
         if self.r is None:
             mesh = (
@@ -155,7 +156,7 @@ def compute(self, soret=False):
             * f.dot(f.grad(self.function), self.n)
             * self.ds(self.surface)
         )
-        if soret and self.field in [0, "0", "solute"]:
+        if self.soret and self.field in [0, "0", "solute"]:
             flux += f.assemble(
                 self.prop
                 * self.r
@@ -234,7 +235,7 @@ def azimuth_range(self, value):
             raise ValueError("Azimuthal range must be between 0 and pi")
         self._azimuth_range = value
 
-    def compute(self, soret=False):
+    def compute(self):
 
         if self.r is None:
             mesh = (
@@ -252,7 +253,7 @@ def compute(self, soret=False):
             * f.dot(f.grad(self.function), self.n)
             * self.ds(self.surface)
         )
-        if soret and self.field in [0, "0", "solute"]:
+        if self.soret and self.field in [0, "0", "solute"]:
             flux += f.assemble(
                 self.prop
                 * self.r**2

festim/generic_simulation.py

@@ -418,6 +418,14 @@ def initialise(self):
                         self.dt.milestones.append(time)
                 self.dt.milestones.sort()
 
+        # set Soret to True
+        for export in self.exports:
+            if isinstance(export, festim.DerivedQuantities):
+                for q in export:
+                    if isinstance(q, festim.SurfaceFlux):
+                        q.soret = self.settings.soret
+                        q.T = self.T.T
+
     def run(self, completion_tone=False):
         """Runs the model.
 

test/system/test_misc.py

@@ -788,3 +788,54 @@ def test_catch_bug_804():
 
     for trap in model.traps:
         assert trap.id is not None
+
+
+def test_soret_surface_flux_mass_balance():
+    """
+    Test to catch bug 830
+
+    1D steady state simulation with Soret on.
+    Compute the surface flux on both surfaces.
+    The surface fluxes should be equal in magnitude
+    but opposite in sign.
+
+    The case is made so that the concentration profile
+    is flat (ie. diffusive flux is zero) but there's a
+    high T gradient (ie. Soret flux is not zero).
+    """
+    my_model = F.Simulation()
+
+    my_model.mesh = F.MeshFromVertices(vertices=np.linspace(0, 0.05, num=100))
+    Q = 2
+    D = 3
+    c = 2
+    my_model.materials = F.Material(id=1, D_0=D, E_D=0, Q=Q)
+    grad_T = 1000
+
+    T_val = lambda x: 700 + grad_T * x
+    my_model.T = F.Temperature(value=T_val(F.x))
+
+    my_model.boundary_conditions = [
+        F.DirichletBC(surfaces=[1], value=c, field=0),
+        F.DirichletBC(surfaces=[2], value=c, field=0),
+    ]
+    my_model.settings = F.Settings(
+        absolute_tolerance=1e-10,
+        relative_tolerance=1e-10,
+        transient=False,
+        soret=True,
+    )
+
+    flux_left = F.SurfaceFlux(field=0, surface=1)
+    flux_right = F.SurfaceFlux(field=0, surface=2)
+    my_model.exports = [F.DerivedQuantities([flux_left, flux_right])]
+
+    my_model.initialise()
+    my_model.run()
+
+    # these two should be equal in steady state
+    print(flux_left.data[0])
+    print(flux_right.data[0])
+
+    assert not np.isclose(flux_left.data[0], 0)
+    assert np.isclose(np.abs(flux_left.data[0]), np.abs(flux_right.data[0]), rtol=1e-2)

test/unit/test_exports/test_derived_quantities/test_surface_flux.py

@@ -53,6 +53,7 @@ class TestCompute:
     my_h_flux.ds = ds
     my_h_flux.T = T
     my_h_flux.Q = Q
+    my_h_flux.soret = True
 
     my_heat_flux = SurfaceFlux("T", surface)
     my_heat_flux.D = D
@@ -87,7 +88,7 @@ def test_h_flux_with_soret(self):
             * f.dot(f.grad(self.T), self.n)
             * self.ds(self.surface)
         )
-        flux = self.my_h_flux.compute(soret=True)
+        flux = self.my_h_flux.compute()
         assert flux == expected_flux
 
 
@@ -149,6 +150,7 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom, soret):
 
     my_flux.n = f.FacetNormal(mesh_fenics)
     my_flux.ds = ds
+    my_flux.soret = soret
 
     expected_value = (
         -2
@@ -177,7 +179,7 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom, soret):
             * (0.5 * r1**2 - 0.5 * r0**2)
         )
 
-    computed_value = my_flux.compute(soret=soret)
+    computed_value = my_flux.compute()
 
     assert np.isclose(computed_value, expected_value)
 
@@ -236,6 +238,7 @@ def test_compute_spherical(r0, radius, c_left, c_right, soret):
 
     my_flux.n = f.FacetNormal(mesh_fenics)
     my_flux.ds = ds
+    my_flux.soret = soret
 
     # expected value is the integral of the flux over the surface
     flux_value = float(my_flux.D) * (c_left - c_right) / (r1 - r0)
@@ -259,7 +262,7 @@ def test_compute_spherical(r0, radius, c_left, c_right, soret):
             * r1**2
         )
 
-    computed_value = my_flux.compute(soret=soret)
+    computed_value = my_flux.compute()
 
     assert np.isclose(computed_value, expected_value)