diff --git a/festim/concentration/mobile.py b/festim/concentration/mobile.py index b8d80baa0..ec2f666f1 100644 --- a/festim/concentration/mobile.py +++ b/festim/concentration/mobile.py @@ -56,9 +56,6 @@ def create_diffusion_form( soret (bool, optional): If True, Soret effect is assumed. Defaults to False. """ - if soret and mesh.type in ["cylindrical", "spherical"]: - msg = "Soret effect not implemented in {} coordinates".format(mesh.type) - raise ValueError(msg) F = 0 for material in materials: @@ -96,6 +93,19 @@ def create_diffusion_form( r = SpatialCoordinate(mesh.mesh)[0] F += r * dot(D * grad(c_0), grad(self.test_function / r)) * dx + if soret: + Q = material.Q + if callable(Q): + Q = Q(T.T) + F += ( + r + * dot( + D * Q * c_0 / (k_B * T.T**2) * grad(T.T), + grad(self.test_function / r), + ) + * dx + ) + elif mesh.type == "spherical": r = SpatialCoordinate(mesh.mesh)[0] F += ( @@ -106,6 +116,21 @@ def create_diffusion_form( * dx ) + if soret: + Q = material.Q + if callable(Q): + Q = Q(T.T) + F += ( + D + * r + * r + * dot( + Q * c_0 / (k_B * T.T**2) * grad(T.T), + grad(self.test_function / r / r), + ) + * dx + ) + # add the trapping terms F_trapping = 0 if traps is not None: diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py index c81577053..04212000e 100644 --- a/festim/exports/derived_quantities/surface_flux.py +++ b/festim/exports/derived_quantities/surface_flux.py @@ -137,10 +137,6 @@ def azimuth_range(self, value): self._azimuth_range = value def compute(self, soret=False): - if soret: - raise NotImplementedError( - "Soret effect not implemented for cylindrical coordinates" - ) if self.r is None: mesh = ( @@ -159,6 +155,16 @@ 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"]: + flux += f.assemble( + self.prop + * self.r + * self.function + * self.Q + / (k_B * self.T**2) + * f.dot(f.grad(self.T), self.n) + * self.ds(self.surface) + ) flux *= self.azimuth_range[1] - self.azimuth_range[0] return flux @@ -229,10 +235,6 @@ def azimuth_range(self, value): self._azimuth_range = value def compute(self, soret=False): - if soret: - raise NotImplementedError( - "Soret effect not implemented for spherical coordinates" - ) if self.r is None: mesh = ( @@ -250,6 +252,16 @@ 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"]: + flux += f.assemble( + self.prop + * self.r**2 + * self.function + * self.Q + / (k_B * self.T**2) + * f.dot(f.grad(self.T), self.n) + * self.ds(self.surface) + ) flux *= (self.polar_range[1] - self.polar_range[0]) * ( -np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0]) ) diff --git a/test/system/test_cylindrical_coordinates.py b/test/system/test_cylindrical_coordinates.py index e29014f30..1dabcd9b4 100644 --- a/test/system/test_cylindrical_coordinates.py +++ b/test/system/test_cylindrical_coordinates.py @@ -112,3 +112,114 @@ def test_temperature_MMS(): error_L2 = fenics.errornorm(expected_solution, produced_solution, "L2") print(error_L2) assert error_L2 < 1e-7 + + +def test_MMS_Soret_cylindrical(): + """ + Tests that festim produces the correct concentration field with the Soret flag + in cylindrical coordinates with FluxBC on the right surface and DirichletBC on others + """ + + def grad(u): + """Computes the gradient of a function u. + + Args: + u (sympy.Expr): a sympy function + + Returns: + sympy.Matrix: the gradient of u + """ + return sp.Matrix([sp.diff(u, r), sp.diff(u, z)]) + + def div(u): + """Computes the divergence of a vector field u. + + Args: + u (sympy.Matrix): a sympy vector field + + Returns: + sympy.Expr: the divergence of u + """ + return sp.simplify(sp.diff(r * u[0], r) / r + sp.diff(u[1], z)) + + # Create and mark the mesh + fenics_mesh = fenics.RectangleMesh(fenics.Point(1, 0), fenics.Point(2, 1), 100, 100) + + volume_markers = fenics.MeshFunction( + "size_t", fenics_mesh, fenics_mesh.topology().dim() + ) + volume_markers.set_all(1) + + right_surface = fenics.CompiledSubDomain("near(x[0], 2.0)") + other_surface = fenics.CompiledSubDomain( + "near(x[0], 1.0) || near(x[1], 0.0) || near(x[1], 1.0)" + ) + + surface_markers = fenics.MeshFunction( + "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1 + ) + + surface_markers.set_all(0) + left_surface_id = 1 + other_surface_id = 2 + + right_surface.mark(surface_markers, left_surface_id) + other_surface.mark(surface_markers, other_surface_id) + + # Create the FESTIM model + my_model = festim.Simulation() + + my_model.mesh = festim.Mesh( + fenics_mesh, + volume_markers=volume_markers, + surface_markers=surface_markers, + type="cylindrical", + ) + + # Variational formulation + r = festim.x + z = festim.y + exact_solution = 1 + 7 * r**2 + 3 * z**2 # exact solution + + T = 300 + 30 * r**2 + 40 * z + + D = 2 + Q = lambda T: 4 * festim.k_B * T + + flux = -D * ( + grad(exact_solution) + Q(T) * exact_solution / (festim.k_B * T**2) * grad(T) + ) + mms_source = div(flux) + + my_model.sources = [ + festim.Source( + mms_source, + volume=1, + field="solute", + ), + ] + + my_model.boundary_conditions = [ + festim.FluxBC(surfaces=left_surface_id, value=-flux[0], field=0), + festim.DirichletBC(surfaces=other_surface_id, value=exact_solution, field=0), + ] + + my_model.materials = festim.Material(id=1, D_0=D, E_D=0, Q=Q) + + my_model.T = festim.Temperature(T) + + my_model.settings = festim.Settings( + absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False, soret=True + ) + + my_model.initialise() + my_model.run() + + expected_solution = fenics.Expression(sp.printing.ccode(exact_solution), degree=4) + expected_solution = fenics.project( + expected_solution, fenics.FunctionSpace(my_model.mesh.mesh, "CG", 1) + ) + + produced_solution = my_model.h_transport_problem.mobile.post_processing_solution + error_L2 = fenics.errornorm(expected_solution, produced_solution, "L2") + assert error_L2 < 2e-4 diff --git a/test/system/test_spherical_coordinates.py b/test/system/test_spherical_coordinates.py index 283ca1b61..5e2e0ea13 100644 --- a/test/system/test_spherical_coordinates.py +++ b/test/system/test_spherical_coordinates.py @@ -110,3 +110,85 @@ def test_MMS_temperature(): expected_solution = fenics.interpolate(u_expr, my_sim.h_transport_problem.V) error_L2 = fenics.errornorm(expected_solution, produced_solution, "L2") assert error_L2 < 1e-7 + + +def test_MMS_Soret_spherical(): + """ + Tests that festim produces the correct concentration field with the Soret flag + in spherical coordinates with DirichletBC on the left surface and FluxBC on the right surface + """ + + def grad(u): + """Computes the gradient of a function u. + + Args: + u (sympy.Expr): a sympy function + + Returns: + sympy.Matrix: the gradient of u + """ + return sp.diff(u, r) + + def div(u): + """Computes the divergence of a vector field u. + + Args: + u (sympy.Matrix): a sympy vector field + + Returns: + sympy.Expr: the divergence of u + """ + return sp.simplify(sp.diff(r**2 * u, r) / r**2) + + # Create the FESTIM model + my_model = festim.Simulation() + + my_model.mesh = festim.MeshFromVertices(np.linspace(1, 2, 100), type="spherical") + + # Variational formulation + r = festim.x + + exact_solution = 1 + 7 * r**2 # exact solution + + T = 300 + 20 * r**2 + + D = 2 + Q = lambda T: 4 * festim.k_B * T + + flux = -D * ( + grad(exact_solution) + Q(T) * exact_solution / (festim.k_B * T**2) * grad(T) + ) + mms_source = div(flux) + + my_model.sources = [ + festim.Source( + mms_source, + volume=1, + field="solute", + ), + ] + + my_model.boundary_conditions = [ + festim.DirichletBC(surfaces=1, value=exact_solution, field=0), + festim.FluxBC(surfaces=2, value=-flux, field=0), + ] + + my_model.materials = festim.Material(id=1, D_0=D, E_D=0, Q=Q) + + my_model.T = festim.Temperature(T) + + my_model.settings = festim.Settings( + absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False, soret=True + ) + + my_model.initialise() + my_model.run() + + expected_solution = fenics.Expression(sp.printing.ccode(exact_solution), degree=4) + expected_solution = fenics.project( + expected_solution, fenics.FunctionSpace(my_model.mesh.mesh, "CG", 1) + ) + + produced_solution = my_model.h_transport_problem.mobile.post_processing_solution + error_L2 = fenics.errornorm(expected_solution, produced_solution, "L2") + assert error_L2 < 1e-4 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 b6d3323be..c5b20ca54 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 @@ -3,8 +3,9 @@ import math import numpy as np import pytest -from festim import SurfaceFlux, k_B +from festim import SurfaceFlux, k_B, x, y from .tools import c_1D, c_2D, c_3D +from sympy.printing import ccode @pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)]) @@ -95,7 +96,8 @@ def test_h_flux_with_soret(self): @pytest.mark.parametrize("height", [2, 3]) @pytest.mark.parametrize("c_top", [2, 3]) @pytest.mark.parametrize("c_bottom", [2, 3]) -def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): +@pytest.mark.parametrize("soret", [False, True]) +def test_compute_cylindrical(r0, radius, height, c_top, c_bottom, soret): """ Test that SurfaceFluxCylindrical computes the flux correctly on a hollow cylinder @@ -105,6 +107,7 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): height (float): cylinder height c_top (float): concentration top c_bottom (float): concentration bottom + soret (bool): if True, the Soret effect will be set """ # creating a mesh with FEniCS r1 = r0 + radius @@ -137,8 +140,9 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): my_flux = SurfaceFluxCylindrical("solute", top_id) my_flux.D = f.Constant(2) V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda z: c_bottom + (c_top - c_bottom) / (z1 - z0) * z expr = f.Expression( - f"{c_bottom} + {(c_top - c_bottom)/(z1-z0)} * x[1]", + ccode(c_fun(y)), degree=1, ) my_flux.function = f.interpolate(expr, V) @@ -146,7 +150,6 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): my_flux.n = f.FacetNormal(mesh_fenics) my_flux.ds = ds - computed_value = my_flux.compute() expected_value = ( -2 * math.pi @@ -155,6 +158,27 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): / height * (0.5 * r1**2 - 0.5 * r0**2) ) + + if soret: + my_flux.Q = f.Constant(2 * k_B) + growth_rate = 7 + T = lambda z: 3 + growth_rate * z + T_expr = f.Expression(ccode(T(y)), degree=1) + my_flux.T = f.interpolate(T_expr, V) + expected_value += ( + 2 + * math.pi + * float(my_flux.D) + * float(my_flux.Q) + * c_fun(z1) + / k_B + / T(z1) ** 2 + * growth_rate + * (0.5 * r1**2 - 0.5 * r0**2) + ) + + computed_value = my_flux.compute(soret=soret) + assert np.isclose(computed_value, expected_value) @@ -162,7 +186,8 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): @pytest.mark.parametrize("radius", [1, 2]) @pytest.mark.parametrize("c_left", [3, 4]) @pytest.mark.parametrize("c_right", [1, 2]) -def test_compute_spherical(r0, radius, c_left, c_right): +@pytest.mark.parametrize("soret", [False, True]) +def test_compute_spherical(r0, radius, c_left, c_right, soret): """ Test that SurfaceFluxSpherical computes the flux correctly on a hollow sphere @@ -171,6 +196,7 @@ def test_compute_spherical(r0, radius, c_left, c_right): radius (float): sphere radius c_left (float): concentration left c_right (float): concentration right + soret (bool): if True, the Soret effect will be set """ # creating a mesh with FEniCS r1 = r0 + radius @@ -201,8 +227,9 @@ def test_compute_spherical(r0, radius, c_left, c_right): my_flux = SurfaceFluxSpherical("solute", right_id) my_flux.D = f.Constant(2) V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda r: c_left + (c_right - c_left) / (r1 - r0) * r expr = f.Expression( - f"{c_left} + {(c_right - c_left)/(r1-r0)} * x[0]", + ccode(c_fun(x)), degree=1, ) my_flux.function = f.interpolate(expr, V) @@ -210,10 +237,30 @@ def test_compute_spherical(r0, radius, c_left, c_right): my_flux.n = f.FacetNormal(mesh_fenics) my_flux.ds = ds - computed_value = my_flux.compute() # expected value is the integral of the flux over the surface flux_value = float(my_flux.D) * (c_left - c_right) / (r1 - r0) expected_value = -4 * math.pi * flux_value * r1**2 + + if soret: + growth_rate = 3 + T = lambda r: 10 + growth_rate * r + T_expr = f.Expression(ccode(T(x)), degree=1) + my_flux.Q = f.Constant(5 * k_B) + my_flux.T = f.interpolate(T_expr, V) + expected_value += ( + 4 + * math.pi + * float(my_flux.D) + * float(my_flux.Q) + * c_fun(r1) + / k_B + / T(r1) ** 2 + * growth_rate + * r1**2 + ) + + computed_value = my_flux.compute(soret=soret) + assert np.isclose(computed_value, expected_value) @@ -228,15 +275,6 @@ def test_azimuthal_range_cylindrical(azimuth_range): SurfaceFluxCylindrical("solute", 1, azimuth_range=azimuth_range) -def test_soret_raises_error(): - """ - Tests that an error is raised when the Soret effect is used with SurfaceFluxCylindrical - """ - my_flux = SurfaceFluxCylindrical("T", 1) - with pytest.raises(NotImplementedError): - my_flux.compute(soret=True) - - @pytest.mark.parametrize( "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)] ) @@ -259,15 +297,6 @@ def test_polar_range_spherical(polar_range): SurfaceFluxSpherical("solute", 1, polar_range=polar_range) -def test_soret_raises_error_spherical(): - """ - Tests that an error is raised when the Soret effect is used with SurfaceFluxSpherical - """ - my_flux = SurfaceFluxSpherical("T", 1) - with pytest.raises(NotImplementedError): - my_flux.compute(soret=True) - - @pytest.mark.parametrize( "function, field, expected_title", [ diff --git a/test/unit/test_mobile.py b/test/unit/test_mobile.py index adbb0d709..17a1d4c13 100644 --- a/test/unit/test_mobile.py +++ b/test/unit/test_mobile.py @@ -258,20 +258,6 @@ def test_with_trap_conglo_transient(self): print(my_mobile.F) assert my_mobile.F.equals(expected_form) - def test_error_soret_cylindrical_spherical(self): - """Tests that the appropriate error is raised when trying to use Soret - with cylindrical or spherical system - """ - my_mobile = festim.Mobile() - - for system in ["cylindrical", "spherical"]: - mesh = festim.Mesh(type=system) - expected_error_msg = "not implemented " + "in {} coordinates".format(system) - with pytest.raises(ValueError, match=expected_error_msg): - my_mobile.create_diffusion_form( - materials=None, mesh=mesh, T=None, soret=True - ) - class TestInitialise: """