Soret form for cylidrical and spherical coordinates and MMS tests

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:

test/system/test_cylindrical_coordinates.py

@@ -112,3 +112,107 @@ def test_temperature_MMS():
     error_L2 = fenics.errornorm(expected_solution, produced_solution, "L2")
     print(error_L2)
     assert error_L2 < 1e-7
+
+
+def test_Soret_MMS():
+    """
+    Tests that festim produces the correct concentration field with the Soret flag
+    in cylindrical coordinates
+    """
+
+    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.UnitSquareMesh(250, 250)
+
+    volume_markers = fenics.MeshFunction(
+        "size_t", fenics_mesh, fenics_mesh.topology().dim()
+    )
+    volume_markers.set_all(1)
+
+    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)
+    other_surface.mark(surface_markers, 1)
+
+    # 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 = 4 * festim.k_B
+
+    mms_source = -div(D * grad(exact_solution)) - div(
+        D * Q * exact_solution / (festim.k_B * T**2) * grad(T)
+    )
+
+    my_model.sources = [
+        festim.Source(
+            mms_source,
+            volume=1,
+            field="solute",
+        ),
+    ]
+
+    my_model.boundary_conditions = [
+        festim.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
+    ]
+
+    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 < 5e-5

test/system/test_spherical_coordinates.py

@@ -110,3 +110,83 @@ 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():
+    """
+    Tests that festim produces the correct concentration field with the Soret flag
+    in spherical coordinates
+    """
+
+    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(0, 1, 1000), type="spherical")
+
+    # Variational formulation
+    r = festim.x
+
+    exact_solution = 1 + 7 * r**2  # exact solution
+
+    T = 300 + 20 * r**2
+
+    D = 2
+    Q = 4 * festim.k_B
+
+    mms_source = -div(D * grad(exact_solution)) - div(
+        D * Q * exact_solution / (festim.k_B * T**2) * grad(T)
+    )
+
+    my_model.sources = [
+        festim.Source(
+            mms_source,
+            volume=1,
+            field="solute",
+        ),
+    ]
+
+    my_model.boundary_conditions = [
+        festim.DirichletBC(surfaces=[2], value=exact_solution, field="solute"),
+    ]
+
+    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 < 5e-6

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:
     """