diff --git a/festim/exports/derived_quantities/derived_quantities.py b/festim/exports/derived_quantities/derived_quantities.py index bbdbf5c8e..28a8a3608 100644 --- a/festim/exports/derived_quantities/derived_quantities.py +++ b/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)): diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index 04212000e..d548a26f2 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/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 diff --git a/festim/generic_simulation.py b/festim/generic_simulation.py index ff1607fd2..96d35037a 100644 --- a/festim/generic_simulation.py +++ b/festim/generic_simulation.py @@ -418,6 +418,13 @@ def initialise(self): self.dt.milestones.append(time) self.dt.milestones.sort() + # set Soret to True for SurfaceFlux quantities + 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. diff --git a/test/system/test_misc.py b/test/system/test_misc.py index 26cebb8da..410d02137 100644 --- a/test/system/test_misc.py +++ b/test/system/test_misc.py @@ -788,3 +788,65 @@ def test_catch_bug_804(): for trap in model.traps: assert trap.id is not None + + +@pytest.mark.parametrize( + "coordinates, surface_flux_class", + [ + ("cartesian", F.SurfaceFlux), + ("cylindrical", F.SurfaceFluxCylindrical), + ("spherical", F.SurfaceFluxSpherical), + ], +) +def test_soret_surface_flux_mass_balance(coordinates, surface_flux_class): + """ + 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() + + # mesh doesn't start at zero for non-cartesian coordinates + my_model.mesh = F.MeshFromVertices( + vertices=np.linspace(0.01, 0.05, num=600), type=coordinates + ) + 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 = surface_flux_class(field=0, surface=1) + flux_right = surface_flux_class(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) diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py index c5b20ca54..1d92d3a44 100644 --- a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py +++ b/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)