Merge pull request #823 from KulaginVladimir/Soret_in_cylindrical_sph…

…erical

Soret effect for cylindrical and spherical systems

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:

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])
         )

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

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