From c690db2952347b51f8f80995c1e64322e3d34c0a Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Sat, 10 Aug 2024 15:00:15 -0400
Subject: [PATCH 01/10] simple cylindrical case

---
 report/verification/mms/simple_cylindrical.md | 282 ++++++++++++++++++
 1 file changed, 282 insertions(+)
 create mode 100644 report/verification/mms/simple_cylindrical.md

diff --git a/report/verification/mms/simple_cylindrical.md b/report/verification/mms/simple_cylindrical.md
new file mode 100644
index 0000000..f9a6757
--- /dev/null
+++ b/report/verification/mms/simple_cylindrical.md
@@ -0,0 +1,282 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.4
+kernelspec:
+  display_name: vv-festim-report-env
+  language: python
+  name: python3
+---
+
+# Simple diffusion case
+
+```{tags} 2D, MMS, steady state, cylindrical
+```
+
+This is a simple MMS example.
+We will only consider diffusion of hydrogen in a unit square domain $\Omega$ at steady state with an homogeneous diffusion coefficient $D$.
+Moreover, a Dirichlet boundary condition will be assumed on the boundaries $\partial \Omega $.
+
+The problem is therefore:
+
+$$
+\begin{align}
+    & \nabla \cdot (D \ \nabla{c}) = -S  \quad \text{on }  \Omega  \\
+    & c = c_0 \quad \text{on }  \partial \Omega
+\end{align}
+$$(problem_simple_cylindrical)
+
+The exact solution for mobile concentration is:
+
+$$
+\begin{equation}
+    c_\mathrm{exact} = 1 + r^2
+\end{equation}
+$$(c_exact_simple_cylindrical)
+
+Injecting {eq}`c_exact_simple_cylindrical` in {eq}`problem_simple_cylindrical`, we obtain the expressions of $S$ and $c_0$:
+
+$$
+\begin{align}
+    & S = -4D \\
+    & c_0 = c_\mathrm{exact}
+\end{align}
+$$
+
+We can then run a FESTIM model with these values and compare the numerical solution with $c_\mathrm{exact}$.
+
++++
+
+## FESTIM code
+
+```{code-cell} ipython3
+import festim as F
+import sympy as sp
+import fenics as f
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Create and mark the mesh
+nx = ny = 100
+fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(1, np.pi / 2), nx, ny)
+
+volume_markers = f.MeshFunction(
+    "size_t", fenics_mesh, fenics_mesh.topology().dim()
+)
+volume_markers.set_all(1)
+
+surface_markers = f.MeshFunction(
+    "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1
+)
+surface_markers.set_all(0)
+
+        
+class Boundary(f.SubDomain):
+    def inside(self, x, on_boundary):
+        return on_boundary
+
+
+boundary = Boundary()
+boundary.mark(surface_markers, 1)
+
+# Create the FESTIM model
+my_model = F.Simulation()
+
+my_model.mesh = F.Mesh(
+    fenics_mesh,
+    volume_markers=volume_markers,
+    surface_markers=surface_markers,
+    type="cylindrical",
+)
+
+r = F.x
+exact_solution = r**2
+D = 2
+
+my_model.sources = [
+    F.Source(- 4 * D, volume=1, field="solute"),
+]
+
+my_model.boundary_conditions = [
+    F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
+]
+
+my_model.materials = F.Material(id=1, D_0=D, E_D=0)
+        
+my_model.T = F.Temperature(500)  # ignored in this problem
+
+my_model.settings = F.Settings(
+    absolute_tolerance=1e-10,
+    relative_tolerance=1e-10,
+    transient=False,
+)
+
+my_model.initialise()
+my_model.run()
+```
+
+## Comparison with exact solution
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+c_exact = f.Expression(sp.printing.ccode(exact_solution), degree=4)
+c_exact = f.project(c_exact, f.FunctionSpace(fenics_mesh, "CG", 1))
+
+computed_solution = my_model.h_transport_problem.mobile.post_processing_solution
+# turn computed solution to cartesian coordinates
+
+
+E = f.errornorm(computed_solution, c_exact, "L2")
+print(f"L2 error: {E:.2e}")
+
+# plot exact solution and computed solution
+fig, axs = plt.subplots(1, 3, figsize=(15, 5))
+plt.sca(axs[0])
+plt.title("Exact solution")
+plt.xlabel("r")
+plt.ylabel("$\\theta$")
+CS1 = f.plot(c_exact, cmap="inferno")
+plt.sca(axs[1])
+plt.xlabel("r")
+plt.title("Computed solution")
+CS2 = f.plot(computed_solution, cmap="inferno")
+
+plt.colorbar(CS1, ax=[axs[0]], shrink=0.8)
+plt.colorbar(CS2, ax=[axs[1]], shrink=0.8)
+
+axs[0].sharey(axs[1])
+plt.setp(axs[1].get_yticklabels(), visible=False)
+
+for CS in [CS1, CS2]:
+    CS.set_edgecolor("face")
+
+
+def compute_arc_length(xs, ys):
+    """Computes the arc length of x,y points based
+    on x and y arrays
+    """
+    points = np.vstack((xs, ys)).T
+    distance = np.linalg.norm(points[1:] - points[:-1], axis=1)
+    arc_length = np.insert(np.cumsum(distance), 0, [0.0])
+    return arc_length
+
+
+# define the profiles
+profiles = [
+    {"start": (0.0, 0.0), "end": (1.0, np.pi/2)},
+    {"start": (0.2, np.pi/6), "end": (0.7, np.pi/4)},
+    {"start": (0.5, 0.1), "end": (0.8, 1.5)},
+]
+
+# plot the profiles on the right subplot
+for i, profile in enumerate(profiles):
+    start_x, start_y = profile["start"]
+    end_x, end_y = profile["end"]
+    plt.sca(axs[1])
+    (l,) = plt.plot([start_x, end_x], [start_y, end_y])
+
+    plt.sca(axs[2])
+
+    points_x_exact = np.linspace(start_x, end_x, num=30)
+    points_y_exact = np.linspace(start_y, end_y, num=30)
+    arc_length_exact = compute_arc_length(points_x_exact, points_y_exact)
+    u_values = [c_exact(x, y) for x, y in zip(points_x_exact, points_y_exact)]
+
+    points_x = np.linspace(start_x, end_x, num=100)
+    points_y = np.linspace(start_y, end_y, num=100)
+    arc_lengths = compute_arc_length(points_x, points_y)
+    computed_values = [computed_solution(x, y) for x, y in zip(points_x, points_y)]
+
+    (exact_line,) = plt.plot(
+        arc_length_exact,
+        u_values,
+        color=l.get_color(),
+        marker="o",
+        linestyle="None",
+        alpha=0.3,
+    )
+    (computed_line,) = plt.plot(arc_lengths, computed_values, color=l.get_color())
+
+plt.sca(axs[2])
+plt.xlabel("Arc length")
+plt.ylabel("Solution")
+
+legend_marker = mpl.lines.Line2D(
+    [],
+    [],
+    color="black",
+    marker=exact_line.get_marker(),
+    linestyle="None",
+    label="Exact",
+)
+legend_line = mpl.lines.Line2D([], [], color="black", label="Computed")
+plt.legend(
+    [legend_marker, legend_line], [legend_marker.get_label(), legend_line.get_label()]
+)
+
+plt.grid(alpha=0.3)
+plt.gca().spines[["right", "top"]].set_visible(False)
+plt.show()
+```
+
+## Compute convergence rates
+
+It is also possible to compute how the numerical error decreases as we increase the number of cells.
+By iteratively refining the mesh, we find that the error exhibits a second order convergence rate.
+This is expected for this particular problem as first order finite elements are used.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+errors = []
+ns = [5, 10, 20, 30, 50, 100, 150]
+
+for n in ns:
+    nx = ny = n
+    fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(1, np.pi / 2), nx, ny)
+
+    volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
+    volume_markers.set_all(1)
+
+    surface_markers = f.MeshFunction(
+        "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1
+    )
+    surface_markers.set_all(0)
+
+    class Boundary(f.SubDomain):
+        def inside(self, x, on_boundary):
+            return on_boundary
+
+    boundary = Boundary()
+    boundary.mark(surface_markers, 1)
+
+    my_model.mesh = F.Mesh(
+        fenics_mesh, volume_markers=volume_markers, surface_markers=surface_markers, type="cylindrical"
+    )
+
+    my_model.initialise()
+    my_model.run()
+
+    computed_solution = my_model.h_transport_problem.mobile.post_processing_solution
+    errors.append(f.errornorm(computed_solution, c_exact, "L2"))
+
+h = 1 / np.array(ns)
+
+plt.loglog(h, errors, marker="o")
+plt.xlabel("Element size")
+plt.ylabel("L2 error")
+
+plt.loglog(h, 2 * h**2, linestyle="--", color="black")
+plt.annotate(
+    "2nd order", (h[0], 2 * h[0] ** 2), textcoords="offset points", xytext=(10, 0)
+)
+
+plt.grid(alpha=0.3)
+plt.gca().spines[["right", "top"]].set_visible(False)
+```

From d1ecd7d9306ee0f293f9f17f389bfbb12fcd7943 Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Sat, 10 Aug 2024 21:11:30 -0400
Subject: [PATCH 02/10] changed to unit disk + added spherical case (not
 working atm)

---
 ...cylindrical.md => simple_non_cartesian.md} | 177 +++++++-----------
 1 file changed, 70 insertions(+), 107 deletions(-)
 rename report/verification/mms/{simple_cylindrical.md => simple_non_cartesian.md} (55%)

diff --git a/report/verification/mms/simple_cylindrical.md b/report/verification/mms/simple_non_cartesian.md
similarity index 55%
rename from report/verification/mms/simple_cylindrical.md
rename to report/verification/mms/simple_non_cartesian.md
index f9a6757..118451f 100644
--- a/report/verification/mms/simple_cylindrical.md
+++ b/report/verification/mms/simple_non_cartesian.md
@@ -12,13 +12,13 @@ kernelspec:
   name: python3
 ---
 
-# Simple diffusion case
+# Simple non-cartesian diffusion cases
 
-```{tags} 2D, MMS, steady state, cylindrical
+```{tags} 2D, MMS, steady state, cylindrical, spherical
 ```
 
-This is a simple MMS example.
-We will only consider diffusion of hydrogen in a unit square domain $\Omega$ at steady state with an homogeneous diffusion coefficient $D$.
+This is a simple MMS example on both cylindrical and spherical meshes.
+We will only consider diffusion of hydrogen in a unit disk domain $\Omega$ at steady state with an homogeneous diffusion coefficient $D$.
 Moreover, a Dirichlet boundary condition will be assumed on the boundaries $\partial \Omega $.
 
 The problem is therefore:
@@ -62,12 +62,11 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 # Create and mark the mesh
+R = 1.0
 nx = ny = 100
-fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(1, np.pi / 2), nx, ny)
+fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(R, 2 * np.pi), nx, ny)
 
-volume_markers = f.MeshFunction(
-    "size_t", fenics_mesh, fenics_mesh.topology().dim()
-)
+volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
 volume_markers.set_all(1)
 
 surface_markers = f.MeshFunction(
@@ -75,7 +74,6 @@ surface_markers = f.MeshFunction(
 )
 surface_markers.set_all(0)
 
-        
 class Boundary(f.SubDomain):
     def inside(self, x, on_boundary):
         return on_boundary
@@ -84,40 +82,45 @@ class Boundary(f.SubDomain):
 boundary = Boundary()
 boundary.mark(surface_markers, 1)
 
-# Create the FESTIM model
-my_model = F.Simulation()
-
-my_model.mesh = F.Mesh(
-    fenics_mesh,
-    volume_markers=volume_markers,
-    surface_markers=surface_markers,
-    type="cylindrical",
-)
-
 r = F.x
 exact_solution = r**2
 D = 2
 
-my_model.sources = [
-    F.Source(- 4 * D, volume=1, field="solute"),
-]
 
-my_model.boundary_conditions = [
-    F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
-]
+def run_model(type: str):
 
-my_model.materials = F.Material(id=1, D_0=D, E_D=0)
-        
-my_model.T = F.Temperature(500)  # ignored in this problem
+    # Create the FESTIM model
+    my_model = F.Simulation()
 
-my_model.settings = F.Settings(
-    absolute_tolerance=1e-10,
-    relative_tolerance=1e-10,
-    transient=False,
-)
+    my_model.mesh = F.Mesh(
+        fenics_mesh,
+        volume_markers=volume_markers,
+        surface_markers=surface_markers,
+        type=type,
+    )
+
+    my_model.sources = [
+        F.Source(-4 * D, volume=1, field="solute"),
+    ]
+
+    my_model.boundary_conditions = [
+        F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
+    ]
+
+    my_model.materials = F.Material(id=1, D_0=D, E_D=0)
+
+    my_model.T = F.Temperature(500)  # ignored in this problem
 
-my_model.initialise()
-my_model.run()
+    my_model.settings = F.Settings(
+        absolute_tolerance=1e-10,
+        relative_tolerance=1e-10,
+        transient=False,
+    )
+
+    my_model.initialise()
+    my_model.run()
+
+    return my_model.h_transport_problem.mobile.post_processing_solution
 ```
 
 ## Comparison with exact solution
@@ -128,32 +131,48 @@ my_model.run()
 c_exact = f.Expression(sp.printing.ccode(exact_solution), degree=4)
 c_exact = f.project(c_exact, f.FunctionSpace(fenics_mesh, "CG", 1))
 
-computed_solution = my_model.h_transport_problem.mobile.post_processing_solution
-# turn computed solution to cartesian coordinates
-
+cylindrical_solution = run_model("cylindrical")
+spherical_solution = run_model("spherical")
 
-E = f.errornorm(computed_solution, c_exact, "L2")
-print(f"L2 error: {E:.2e}")
+E_cyl = f.errornorm(cylindrical_solution, c_exact, "L2")
+E_sph = f.errornorm(spherical_solution, c_exact, "L2")
+print(f"L2 error, cylindrical: {E_cyl:.2e}")
+print(f"L2 error, spherical: {E_sph:.2e}")
 
 # plot exact solution and computed solution
-fig, axs = plt.subplots(1, 3, figsize=(15, 5))
+fig, axs = plt.subplots(
+    1,
+    4,
+    figsize=(15, 5),
+    sharey=True,
+)
+
 plt.sca(axs[0])
 plt.title("Exact solution")
 plt.xlabel("r")
 plt.ylabel("$\\theta$")
 CS1 = f.plot(c_exact, cmap="inferno")
+plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
 plt.sca(axs[1])
 plt.xlabel("r")
-plt.title("Computed solution")
-CS2 = f.plot(computed_solution, cmap="inferno")
+plt.title("Cylindrical solution")
+CS2 = f.plot(cylindrical_solution, cmap="inferno")
+plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
+plt.sca(axs[2])
+plt.xlabel("r")
+plt.title("Spherical solution")
+CS3 = f.plot(spherical_solution, cmap="inferno")
+plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
 
 plt.colorbar(CS1, ax=[axs[0]], shrink=0.8)
 plt.colorbar(CS2, ax=[axs[1]], shrink=0.8)
+plt.colorbar(CS3, ax=[axs[2]], shrink=0.8)
 
-axs[0].sharey(axs[1])
+# axs[0].sharey(axs[1])
 plt.setp(axs[1].get_yticklabels(), visible=False)
+plt.setp(axs[2].get_yticklabels(), visible=False)
 
-for CS in [CS1, CS2]:
+for CS in [CS1, CS2, CS3]:
     CS.set_edgecolor("face")
 
 
@@ -169,9 +188,9 @@ def compute_arc_length(xs, ys):
 
 # define the profiles
 profiles = [
-    {"start": (0.0, 0.0), "end": (1.0, np.pi/2)},
-    {"start": (0.2, np.pi/6), "end": (0.7, np.pi/4)},
-    {"start": (0.5, 0.1), "end": (0.8, 1.5)},
+    {"start": (0.0, 0.0), "end": (1.0, 2 * np.pi)},
+    {"start": (0.2, 2 * np.pi / 3), "end": (0.7, np.pi)},
+    {"start": (0.5, 1), "end": (0.8, 4)},
 ]
 
 # plot the profiles on the right subplot
@@ -181,7 +200,7 @@ for i, profile in enumerate(profiles):
     plt.sca(axs[1])
     (l,) = plt.plot([start_x, end_x], [start_y, end_y])
 
-    plt.sca(axs[2])
+    plt.sca(axs[-1])
 
     points_x_exact = np.linspace(start_x, end_x, num=30)
     points_y_exact = np.linspace(start_y, end_y, num=30)
@@ -191,7 +210,7 @@ for i, profile in enumerate(profiles):
     points_x = np.linspace(start_x, end_x, num=100)
     points_y = np.linspace(start_y, end_y, num=100)
     arc_lengths = compute_arc_length(points_x, points_y)
-    computed_values = [computed_solution(x, y) for x, y in zip(points_x, points_y)]
+    computed_values = [cylindrical_solution(x, y) for x, y in zip(points_x, points_y)]
 
     (exact_line,) = plt.plot(
         arc_length_exact,
@@ -203,7 +222,7 @@ for i, profile in enumerate(profiles):
     )
     (computed_line,) = plt.plot(arc_lengths, computed_values, color=l.get_color())
 
-plt.sca(axs[2])
+plt.sca(axs[-1])
 plt.xlabel("Arc length")
 plt.ylabel("Solution")
 
@@ -224,59 +243,3 @@ plt.grid(alpha=0.3)
 plt.gca().spines[["right", "top"]].set_visible(False)
 plt.show()
 ```
-
-## Compute convergence rates
-
-It is also possible to compute how the numerical error decreases as we increase the number of cells.
-By iteratively refining the mesh, we find that the error exhibits a second order convergence rate.
-This is expected for this particular problem as first order finite elements are used.
-
-```{code-cell} ipython3
-:tags: [hide-input]
-
-errors = []
-ns = [5, 10, 20, 30, 50, 100, 150]
-
-for n in ns:
-    nx = ny = n
-    fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(1, np.pi / 2), nx, ny)
-
-    volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
-    volume_markers.set_all(1)
-
-    surface_markers = f.MeshFunction(
-        "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1
-    )
-    surface_markers.set_all(0)
-
-    class Boundary(f.SubDomain):
-        def inside(self, x, on_boundary):
-            return on_boundary
-
-    boundary = Boundary()
-    boundary.mark(surface_markers, 1)
-
-    my_model.mesh = F.Mesh(
-        fenics_mesh, volume_markers=volume_markers, surface_markers=surface_markers, type="cylindrical"
-    )
-
-    my_model.initialise()
-    my_model.run()
-
-    computed_solution = my_model.h_transport_problem.mobile.post_processing_solution
-    errors.append(f.errornorm(computed_solution, c_exact, "L2"))
-
-h = 1 / np.array(ns)
-
-plt.loglog(h, errors, marker="o")
-plt.xlabel("Element size")
-plt.ylabel("L2 error")
-
-plt.loglog(h, 2 * h**2, linestyle="--", color="black")
-plt.annotate(
-    "2nd order", (h[0], 2 * h[0] ** 2), textcoords="offset points", xytext=(10, 0)
-)
-
-plt.grid(alpha=0.3)
-plt.gca().spines[["right", "top"]].set_visible(False)
-```

From 2c3f9c7884d2001b694d199967c89192117c59da Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Tue, 13 Aug 2024 13:29:23 -0400
Subject: [PATCH 03/10] started soret

---
 report/verification/mms/soret.md              |   3 +-
 .../verification/mms/soret_non_cartesian.md   | 245 ++++++++++++++++++
 2 files changed, 246 insertions(+), 2 deletions(-)
 create mode 100644 report/verification/mms/soret_non_cartesian.md

diff --git a/report/verification/mms/soret.md b/report/verification/mms/soret.md
index 9710701..a1a2298 100644
--- a/report/verification/mms/soret.md
+++ b/report/verification/mms/soret.md
@@ -14,7 +14,7 @@ kernelspec:
 
 # Soret Effect
 
-```{tags} 2D, MMS, soret, steady state
+```{tags} 2D, MMS, steady state, soret
 ```
 
 This MMS case verifies the implementation of the Soret effect in FESTIM.
@@ -47,7 +47,6 @@ For this problem, we choose:
     & D = 2
 \end{align}
 
-
 Injecting {eq}`c_exact_soret` in {eq}`problem_soret`, we obtain the expressions of $S$ and $c_0$:
 
 $$
diff --git a/report/verification/mms/soret_non_cartesian.md b/report/verification/mms/soret_non_cartesian.md
new file mode 100644
index 0000000..e1b08b0
--- /dev/null
+++ b/report/verification/mms/soret_non_cartesian.md
@@ -0,0 +1,245 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.4
+kernelspec:
+  display_name: vv-festim-report-env
+  language: python
+  name: python3
+---
+
+# Simple non-cartesian diffusion cases
+
+```{tags} 2D, MMS, steady state, cylindrical, spherical, soret
+```
+
+This MMS case verifies the implementation of the Soret effect in FESTIM on both cylindrical and spherical meshes.
+We will only consider diffusion of hydrogen in a unit disk domain $\Omega$ at steady state with a homogeneous diffusion coefficient $D$.
+Moreover, a Dirichlet boundary condition will be assumed on the boundaries $\partial \Omega $.
+
+The problem is therefore:
+
+$$
+\begin{align}
+    & \nabla \cdot (D \ \nabla{c}) = -S  \quad \text{on }  \Omega  \\
+    & c = c_0 \quad \text{on }  \partial \Omega
+\end{align}
+$$(problem_simple_cylindrical)
+
+The exact solution for mobile concentration is:
+
+$$
+\begin{equation}
+    c_\mathrm{exact} = 1 + r^2
+\end{equation}
+$$(c_exact_simple_cylindrical)
+
+Injecting {eq}`c_exact_simple_cylindrical` in {eq}`problem_simple_cylindrical`, we obtain the expressions of $S$ and $c_0$:
+
+$$
+\begin{align}
+    & S = -4D \\
+    & c_0 = c_\mathrm{exact}
+\end{align}
+$$
+
+We can then run a FESTIM model with these values and compare the numerical solution with $c_\mathrm{exact}$.
+
++++
+
+## FESTIM code
+
+```{code-cell}
+import festim as F
+import sympy as sp
+import fenics as f
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Create and mark the mesh
+R = 1.0
+nx = ny = 100
+fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(R, 2 * np.pi), nx, ny)
+
+volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
+volume_markers.set_all(1)
+
+surface_markers = f.MeshFunction(
+    "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1
+)
+surface_markers.set_all(0)
+
+class Boundary(f.SubDomain):
+    def inside(self, x, on_boundary):
+        return on_boundary
+
+
+boundary = Boundary()
+boundary.mark(surface_markers, 1)
+
+r = F.x
+exact_solution = r**2
+D = 2
+
+
+def run_model(type: str):
+
+    # Create the FESTIM model
+    my_model = F.Simulation()
+
+    my_model.mesh = F.Mesh(
+        fenics_mesh,
+        volume_markers=volume_markers,
+        surface_markers=surface_markers,
+        type=type,
+    )
+
+    my_model.sources = [
+        F.Source(-4 * D, volume=1, field="solute"),
+    ]
+
+    my_model.boundary_conditions = [
+        F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
+    ]
+
+    my_model.materials = F.Material(id=1, D_0=D, E_D=0)
+
+    my_model.T = F.Temperature(500)  # ignored in this problem
+
+    my_model.settings = F.Settings(
+        absolute_tolerance=1e-10,
+        relative_tolerance=1e-10,
+        transient=False,
+    )
+
+    my_model.initialise()
+    my_model.run()
+
+    return my_model.h_transport_problem.mobile.post_processing_solution
+```
+
+## Comparison with exact solution
+
+```{code-cell}
+:tags: [hide-input]
+
+c_exact = f.Expression(sp.printing.ccode(exact_solution), degree=4)
+c_exact = f.project(c_exact, f.FunctionSpace(fenics_mesh, "CG", 1))
+
+cylindrical_solution = run_model("cylindrical")
+spherical_solution = run_model("spherical")
+
+E_cyl = f.errornorm(cylindrical_solution, c_exact, "L2")
+E_sph = f.errornorm(spherical_solution, c_exact, "L2")
+print(f"L2 error, cylindrical: {E_cyl:.2e}")
+print(f"L2 error, spherical: {E_sph:.2e}")
+
+# plot exact solution and computed solution
+fig, axs = plt.subplots(
+    1,
+    4,
+    figsize=(15, 5),
+    sharey=True,
+)
+
+plt.sca(axs[0])
+plt.title("Exact solution")
+plt.xlabel("r")
+plt.ylabel("$\\theta$")
+CS1 = f.plot(c_exact, cmap="inferno")
+plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
+plt.sca(axs[1])
+plt.xlabel("r")
+plt.title("Cylindrical solution")
+CS2 = f.plot(cylindrical_solution, cmap="inferno")
+plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
+plt.sca(axs[2])
+plt.xlabel("r")
+plt.title("Spherical solution")
+CS3 = f.plot(spherical_solution, cmap="inferno")
+plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
+
+plt.colorbar(CS1, ax=[axs[0]], shrink=0.8)
+plt.colorbar(CS2, ax=[axs[1]], shrink=0.8)
+plt.colorbar(CS3, ax=[axs[2]], shrink=0.8)
+
+# axs[0].sharey(axs[1])
+plt.setp(axs[1].get_yticklabels(), visible=False)
+plt.setp(axs[2].get_yticklabels(), visible=False)
+
+for CS in [CS1, CS2, CS3]:
+    CS.set_edgecolor("face")
+
+
+def compute_arc_length(xs, ys):
+    """Computes the arc length of x,y points based
+    on x and y arrays
+    """
+    points = np.vstack((xs, ys)).T
+    distance = np.linalg.norm(points[1:] - points[:-1], axis=1)
+    arc_length = np.insert(np.cumsum(distance), 0, [0.0])
+    return arc_length
+
+
+# define the profiles
+profiles = [
+    {"start": (0.0, 0.0), "end": (1.0, 2 * np.pi)},
+    {"start": (0.2, 2 * np.pi / 3), "end": (0.7, np.pi)},
+    {"start": (0.5, 1), "end": (0.8, 4)},
+]
+
+# plot the profiles on the right subplot
+for i, profile in enumerate(profiles):
+    start_x, start_y = profile["start"]
+    end_x, end_y = profile["end"]
+    plt.sca(axs[1])
+    (l,) = plt.plot([start_x, end_x], [start_y, end_y])
+
+    plt.sca(axs[-1])
+
+    points_x_exact = np.linspace(start_x, end_x, num=30)
+    points_y_exact = np.linspace(start_y, end_y, num=30)
+    arc_length_exact = compute_arc_length(points_x_exact, points_y_exact)
+    u_values = [c_exact(x, y) for x, y in zip(points_x_exact, points_y_exact)]
+
+    points_x = np.linspace(start_x, end_x, num=100)
+    points_y = np.linspace(start_y, end_y, num=100)
+    arc_lengths = compute_arc_length(points_x, points_y)
+    computed_values = [cylindrical_solution(x, y) for x, y in zip(points_x, points_y)]
+
+    (exact_line,) = plt.plot(
+        arc_length_exact,
+        u_values,
+        color=l.get_color(),
+        marker="o",
+        linestyle="None",
+        alpha=0.3,
+    )
+    (computed_line,) = plt.plot(arc_lengths, computed_values, color=l.get_color())
+
+plt.sca(axs[-1])
+plt.xlabel("Arc length")
+plt.ylabel("Solution")
+
+legend_marker = mpl.lines.Line2D(
+    [],
+    [],
+    color="black",
+    marker=exact_line.get_marker(),
+    linestyle="None",
+    label="Exact",
+)
+legend_line = mpl.lines.Line2D([], [], color="black", label="Computed")
+plt.legend(
+    [legend_marker, legend_line], [legend_marker.get_label(), legend_line.get_label()]
+)
+
+plt.grid(alpha=0.3)
+plt.gca().spines[["right", "top"]].set_visible(False)
+plt.show()
+```

From be4b9f8e1afc7b0357ea35b13bee9990db36cf67 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Wed, 14 Aug 2024 10:00:48 +0200
Subject: [PATCH 04/10] changed requirement on FESTIM

---
 environment.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/environment.yml b/environment.yml
index 3399ca1..b7cdac8 100644
--- a/environment.yml
+++ b/environment.yml
@@ -22,4 +22,4 @@ dependencies:
   - scipy
   - pip
   - pip:
-    - festim==1.3
+    - festim~=1.3

From a77719d7092bed1d5cad9b05678ae03f902555d5 Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Tue, 20 Aug 2024 14:52:40 -0400
Subject: [PATCH 05/10] cylindrical only + rename + add to ToC

---
 report/_toc.yml                               |  1 +
 ...non_cartesian.md => simple_cylindrical.md} | 34 ++++++-------------
 2 files changed, 12 insertions(+), 23 deletions(-)
 rename report/verification/mms/{simple_non_cartesian.md => simple_cylindrical.md} (83%)

diff --git a/report/_toc.yml b/report/_toc.yml
index e0f0539..90cc405 100644
--- a/report/_toc.yml
+++ b/report/_toc.yml
@@ -17,6 +17,7 @@ parts:
           - file: verification/mms/simple
           - file: verification/mms/simple_transient
           - file: verification/mms/simple_temperature
+          - file: verification/mms/simple_cylindrical
       - file: verification/mms/fluxes
       - file: verification/mms/heat_transfer
       - file: verification/mms/discontinuity
diff --git a/report/verification/mms/simple_non_cartesian.md b/report/verification/mms/simple_cylindrical.md
similarity index 83%
rename from report/verification/mms/simple_non_cartesian.md
rename to report/verification/mms/simple_cylindrical.md
index 118451f..1e6c104 100644
--- a/report/verification/mms/simple_non_cartesian.md
+++ b/report/verification/mms/simple_cylindrical.md
@@ -12,13 +12,13 @@ kernelspec:
   name: python3
 ---
 
-# Simple non-cartesian diffusion cases
+# Simple cylindrical case
 
-```{tags} 2D, MMS, steady state, cylindrical, spherical
+```{tags} 2D, MMS, steady state, cylindrical
 ```
 
-This is a simple MMS example on both cylindrical and spherical meshes.
-We will only consider diffusion of hydrogen in a unit disk domain $\Omega$ at steady state with an homogeneous diffusion coefficient $D$.
+This is a simple MMS example on a cylindrical mesh.
+We will only consider diffusion of hydrogen in a unit square domain $\Omega$ at steady state with an homogeneous diffusion coefficient $D$.
 Moreover, a Dirichlet boundary condition will be assumed on the boundaries $\partial \Omega $.
 
 The problem is therefore:
@@ -62,9 +62,8 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 # Create and mark the mesh
-R = 1.0
 nx = ny = 100
-fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(R, 2 * np.pi), nx, ny)
+fenics_mesh = f.UnitSquareMesh(nx, ny)
 
 volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
 volume_markers.set_all(1)
@@ -132,17 +131,14 @@ c_exact = f.Expression(sp.printing.ccode(exact_solution), degree=4)
 c_exact = f.project(c_exact, f.FunctionSpace(fenics_mesh, "CG", 1))
 
 cylindrical_solution = run_model("cylindrical")
-spherical_solution = run_model("spherical")
 
 E_cyl = f.errornorm(cylindrical_solution, c_exact, "L2")
-E_sph = f.errornorm(spherical_solution, c_exact, "L2")
 print(f"L2 error, cylindrical: {E_cyl:.2e}")
-print(f"L2 error, spherical: {E_sph:.2e}")
 
 # plot exact solution and computed solution
 fig, axs = plt.subplots(
     1,
-    4,
+    3,
     figsize=(15, 5),
     sharey=True,
 )
@@ -150,29 +146,21 @@ fig, axs = plt.subplots(
 plt.sca(axs[0])
 plt.title("Exact solution")
 plt.xlabel("r")
-plt.ylabel("$\\theta$")
+plt.ylabel("z")
 CS1 = f.plot(c_exact, cmap="inferno")
-plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
 plt.sca(axs[1])
 plt.xlabel("r")
 plt.title("Cylindrical solution")
 CS2 = f.plot(cylindrical_solution, cmap="inferno")
-plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
-plt.sca(axs[2])
-plt.xlabel("r")
-plt.title("Spherical solution")
-CS3 = f.plot(spherical_solution, cmap="inferno")
-plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
 
 plt.colorbar(CS1, ax=[axs[0]], shrink=0.8)
 plt.colorbar(CS2, ax=[axs[1]], shrink=0.8)
-plt.colorbar(CS3, ax=[axs[2]], shrink=0.8)
 
 # axs[0].sharey(axs[1])
 plt.setp(axs[1].get_yticklabels(), visible=False)
 plt.setp(axs[2].get_yticklabels(), visible=False)
 
-for CS in [CS1, CS2, CS3]:
+for CS in [CS1, CS2]:
     CS.set_edgecolor("face")
 
 
@@ -188,9 +176,9 @@ def compute_arc_length(xs, ys):
 
 # define the profiles
 profiles = [
-    {"start": (0.0, 0.0), "end": (1.0, 2 * np.pi)},
-    {"start": (0.2, 2 * np.pi / 3), "end": (0.7, np.pi)},
-    {"start": (0.5, 1), "end": (0.8, 4)},
+    {"start": (0.0, 0.0), "end": (1.0, 1.0)},
+    {"start": (0.2, 0.8), "end": (0.7, 0.2)},
+    {"start": (0.2, 0.6), "end": (0.8, 0.8)},
 ]
 
 # plot the profiles on the right subplot

From 4c5e72acd4396993f66c5d0031039d5549762b4e Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Fri, 23 Aug 2024 13:38:29 -0400
Subject: [PATCH 06/10] changed equation  name to avoid conflict w prev non
 cartesian case

---
 report/verification/mms/soret_non_cartesian.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/report/verification/mms/soret_non_cartesian.md b/report/verification/mms/soret_non_cartesian.md
index e1b08b0..80aef74 100644
--- a/report/verification/mms/soret_non_cartesian.md
+++ b/report/verification/mms/soret_non_cartesian.md
@@ -28,7 +28,7 @@ $$
     & \nabla \cdot (D \ \nabla{c}) = -S  \quad \text{on }  \Omega  \\
     & c = c_0 \quad \text{on }  \partial \Omega
 \end{align}
-$$(problem_simple_cylindrical)
+$$(problem_simple_soret_cylindrical)
 
 The exact solution for mobile concentration is:
 
@@ -36,9 +36,9 @@ $$
 \begin{equation}
     c_\mathrm{exact} = 1 + r^2
 \end{equation}
-$$(c_exact_simple_cylindrical)
+$$(c_exact_simple_soret_cylindrical)
 
-Injecting {eq}`c_exact_simple_cylindrical` in {eq}`problem_simple_cylindrical`, we obtain the expressions of $S$ and $c_0$:
+Injecting {eq}`c_exact_simple_soret_cylindrical` in {eq}`problem_simple_soret_cylindrical`, we obtain the expressions of $S$ and $c_0$:
 
 $$
 \begin{align}

From 8f1b2d0b71e18c86aa69dba943dfbb71d8c70fba Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Fri, 23 Aug 2024 13:54:56 -0400
Subject: [PATCH 07/10] fixed exact_solution in code

---
 report/verification/mms/simple_cylindrical.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/report/verification/mms/simple_cylindrical.md b/report/verification/mms/simple_cylindrical.md
index 1e6c104..53d3626 100644
--- a/report/verification/mms/simple_cylindrical.md
+++ b/report/verification/mms/simple_cylindrical.md
@@ -82,7 +82,7 @@ boundary = Boundary()
 boundary.mark(surface_markers, 1)
 
 r = F.x
-exact_solution = r**2
+exact_solution = 1 + r**2
 D = 2
 
 

From 1df33abc47f3beedb8b151eb248358f3b1265e92 Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Fri, 23 Aug 2024 14:19:02 -0400
Subject: [PATCH 08/10] starts at r = 1

---
 report/verification/mms/simple_cylindrical.md | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/report/verification/mms/simple_cylindrical.md b/report/verification/mms/simple_cylindrical.md
index 53d3626..588c041 100644
--- a/report/verification/mms/simple_cylindrical.md
+++ b/report/verification/mms/simple_cylindrical.md
@@ -18,7 +18,7 @@ kernelspec:
 ```
 
 This is a simple MMS example on a cylindrical mesh.
-We will only consider diffusion of hydrogen in a unit square domain $\Omega$ at steady state with an homogeneous diffusion coefficient $D$.
+We will only consider diffusion of hydrogen in a unit square domain $\Omega$ s.t. $r \in [1, 2]$ and $z \in [0, 1]$ at steady state with a homogeneous diffusion coefficient $D$.
 Moreover, a Dirichlet boundary condition will be assumed on the boundaries $\partial \Omega $.
 
 The problem is therefore:
@@ -63,7 +63,7 @@ import numpy as np
 
 # Create and mark the mesh
 nx = ny = 100
-fenics_mesh = f.UnitSquareMesh(nx, ny)
+fenics_mesh = f.RectangleMesh(f.Point(1, 0), f.Point(2, 1), nx, ny)
 
 volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
 volume_markers.set_all(1)
@@ -140,7 +140,6 @@ fig, axs = plt.subplots(
     1,
     3,
     figsize=(15, 5),
-    sharey=True,
 )
 
 plt.sca(axs[0])
@@ -176,9 +175,9 @@ def compute_arc_length(xs, ys):
 
 # define the profiles
 profiles = [
-    {"start": (0.0, 0.0), "end": (1.0, 1.0)},
-    {"start": (0.2, 0.8), "end": (0.7, 0.2)},
-    {"start": (0.2, 0.6), "end": (0.8, 0.8)},
+    {"start": (1.0, 0.0), "end": (2.0, 1.0)},
+    {"start": (1.2, 0.8), "end": (1.7, 0.2)},
+    {"start": (1.2, 0.6), "end": (1.8, 0.8)},
 ]
 
 # plot the profiles on the right subplot

From 51bfe7da0cfb824e7891fb63acdc78e69d23d2f1 Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Thu, 29 Aug 2024 12:38:41 -0400
Subject: [PATCH 09/10] file for cylindrical soret

---
 .../mms/{soret_non_cartesian.md => soret_cylindrical.md}          | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename report/verification/mms/{soret_non_cartesian.md => soret_cylindrical.md} (100%)

diff --git a/report/verification/mms/soret_non_cartesian.md b/report/verification/mms/soret_cylindrical.md
similarity index 100%
rename from report/verification/mms/soret_non_cartesian.md
rename to report/verification/mms/soret_cylindrical.md

From c61a01238147599fc2c7670f331254eaaebfb28c Mon Sep 17 00:00:00 2001
From: rekomodo <jairrs@mit.edu>
Date: Thu, 29 Aug 2024 14:15:32 -0400
Subject: [PATCH 10/10] cylindrical soret + added to ToC + section for soret
 cases

---
 report/_toc.yml                              |   5 +-
 report/verification/mms/soret_cylindrical.md | 172 +++++++++++--------
 report/verification/soret.md                 |   1 +
 3 files changed, 103 insertions(+), 75 deletions(-)
 create mode 100644 report/verification/soret.md

diff --git a/report/_toc.yml b/report/_toc.yml
index 90cc405..ca1216b 100644
--- a/report/_toc.yml
+++ b/report/_toc.yml
@@ -26,7 +26,10 @@ parts:
         sections:
           - file: verification/mes/radioactive_decay
           - file: verification/mms/radioactive_decay
-      - file: verification/mms/soret
+      - file: verification/soret
+        sections:
+          - file: verification/mms/soret
+          - file: verification/mms/soret_cylindrical
       - file: verification/tmap
         sections:
           - file: verification/mes/depleting_source
diff --git a/report/verification/mms/soret_cylindrical.md b/report/verification/mms/soret_cylindrical.md
index 80aef74..1bbb433 100644
--- a/report/verification/mms/soret_cylindrical.md
+++ b/report/verification/mms/soret_cylindrical.md
@@ -12,37 +12,48 @@ kernelspec:
   name: python3
 ---
 
-# Simple non-cartesian diffusion cases
+# Soret Effect - Cylindrical
 
 ```{tags} 2D, MMS, steady state, cylindrical, spherical, soret
 ```
 
-This MMS case verifies the implementation of the Soret effect in FESTIM on both cylindrical and spherical meshes.
-We will only consider diffusion of hydrogen in a unit disk domain $\Omega$ at steady state with a homogeneous diffusion coefficient $D$.
+This MMS case verifies the implementation of the Soret effect in FESTIM on cylindrical meshes.
+We will only consider diffusion of hydrogen in a unit square domain $\Omega$ s.t. $r \in [1, 2]$ and $z \in [0, 1]$ at steady state with a homogeneous diffusion coefficient $D$.
 Moreover, a Dirichlet boundary condition will be assumed on the boundaries $\partial \Omega $.
 
 The problem is therefore:
 
 $$
 \begin{align}
-    & \nabla \cdot (D \ \nabla{c}) = -S  \quad \text{on }  \Omega  \\
+    & - \nabla\cdot\vec{\varphi} = -S  \quad \text{on }  \Omega  \\
+    & \vec{\varphi} = -D \ \nabla{c} - D\frac{Q^* c}{k_B T^2} \ \nabla{T} \\
     & c = c_0 \quad \text{on }  \partial \Omega
 \end{align}
 $$(problem_simple_soret_cylindrical)
 
-The exact solution for mobile concentration is:
+The manufactured exact solution for mobile concentration is:
 
 $$
 \begin{equation}
-    c_\mathrm{exact} = 1 + r^2
+    c_\mathrm{exact} = 1 + 3r + 5z
 \end{equation}
 $$(c_exact_simple_soret_cylindrical)
 
+For this problem, we choose:
+
+$$
+\begin{align}
+    & T = 300 + 30r^2 + 50z^2 \\
+    & Q^* = 4 \\
+    & D = 2
+\end{align}
+$$
+
 Injecting {eq}`c_exact_simple_soret_cylindrical` in {eq}`problem_simple_soret_cylindrical`, we obtain the expressions of $S$ and $c_0$:
 
 $$
 \begin{align}
-    & S = -4D \\
+    & S = - D\nabla \cdot \left(\frac{Q^* c_\mathrm{exact}}{k_B T^2} \ \nabla{T} \right) -12 D \\
     & c_0 = c_\mathrm{exact}
 \end{align}
 $$
@@ -53,7 +64,7 @@ We can then run a FESTIM model with these values and compare the numerical solut
 
 ## FESTIM code
 
-```{code-cell}
+```{code-cell} ipython3
 import festim as F
 import sympy as sp
 import fenics as f
@@ -62,9 +73,8 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 # Create and mark the mesh
-R = 1.0
 nx = ny = 100
-fenics_mesh = f.RectangleMesh(f.Point(0, 0), f.Point(R, 2 * np.pi), nx, ny)
+fenics_mesh = f.RectangleMesh(f.Point(1.0, 0), f.Point(2.0, 1.0), nx, ny)
 
 volume_markers = f.MeshFunction("size_t", fenics_mesh, fenics_mesh.topology().dim())
 volume_markers.set_all(1)
@@ -72,8 +82,10 @@ volume_markers.set_all(1)
 surface_markers = f.MeshFunction(
     "size_t", fenics_mesh, fenics_mesh.topology().dim() - 1
 )
+
 surface_markers.set_all(0)
 
+
 class Boundary(f.SubDomain):
     def inside(self, x, on_boundary):
         return on_boundary
@@ -82,97 +94,108 @@ class Boundary(f.SubDomain):
 boundary = Boundary()
 boundary.mark(surface_markers, 1)
 
+# Create the FESTIM model
+my_model = F.Simulation()
+
+my_model.mesh = F.Mesh(
+    fenics_mesh, volume_markers=volume_markers, surface_markers=surface_markers, type="cylindrical"
+)
+
+# Variational formulation
 r = F.x
-exact_solution = r**2
+z = F.y
+
+exact_solution = 1 + 3 * r + 5 * z # exact solution
+
+T = 300 + 30*r**2 + 50*z**2
+
 D = 2
+Q = 4
+
+def grad_cylindrical(u):
+    """Computes the gradient of a function u in cylindrical coordinates.
 
+    Args:
+        u (sympy.Expr): a sympy function
 
-def run_model(type: str):
+    Returns:
+        sympy.Matrix: the gradient of u
+    """
 
-    # Create the FESTIM model
-    my_model = F.Simulation()
+    return sp.Matrix([sp.diff(u, r), sp.diff(u, z)])
 
-    my_model.mesh = F.Mesh(
-        fenics_mesh,
-        volume_markers=volume_markers,
-        surface_markers=surface_markers,
-        type=type,
-    )
 
-    my_model.sources = [
-        F.Source(-4 * D, volume=1, field="solute"),
-    ]
+def div_cylindrical(u):
+    """Computes the divergence of a vector field u in cylindrical coordinates.
 
-    my_model.boundary_conditions = [
-        F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
-    ]
+    Args:
+        u (sympy.Matrix): a sympy vector field
 
-    my_model.materials = F.Material(id=1, D_0=D, E_D=0)
+    Returns:
+        sympy.Expr: the divergence of u
+    """
 
-    my_model.T = F.Temperature(500)  # ignored in this problem
+    return 1/r * sp.diff(u[0] * r, r) + sp.diff(u[1], z)
 
-    my_model.settings = F.Settings(
-        absolute_tolerance=1e-10,
-        relative_tolerance=1e-10,
-        transient=False,
-    )
+diffusion_term = -D * grad_cylindrical(exact_solution)
+soret_term = -D * (Q * exact_solution) / (F.k_B * T**2) * grad_cylindrical(T)
 
-    my_model.initialise()
-    my_model.run()
+mms_source = div_cylindrical( diffusion_term + soret_term )
 
-    return my_model.h_transport_problem.mobile.post_processing_solution
+my_model.sources = [
+    F.Source(
+        mms_source,
+        volume=1,
+        field="solute",
+    ),
+]
+
+my_model.boundary_conditions = [
+    F.DirichletBC(surfaces=[1], value=exact_solution, field="solute"),
+]
+
+my_model.materials = F.Material(id=1, D_0=D, E_D=0, Q=Q)
+
+my_model.T = F.Temperature(T)
+
+my_model.settings = F.Settings(
+    absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False, soret=True
+)
+
+my_model.initialise()
+my_model.run()
 ```
 
 ## Comparison with exact solution
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 c_exact = f.Expression(sp.printing.ccode(exact_solution), degree=4)
-c_exact = f.project(c_exact, f.FunctionSpace(fenics_mesh, "CG", 1))
-
-cylindrical_solution = run_model("cylindrical")
-spherical_solution = run_model("spherical")
+c_exact = f.project(c_exact, f.FunctionSpace(my_model.mesh.mesh, "CG", 1))
 
-E_cyl = f.errornorm(cylindrical_solution, c_exact, "L2")
-E_sph = f.errornorm(spherical_solution, c_exact, "L2")
-print(f"L2 error, cylindrical: {E_cyl:.2e}")
-print(f"L2 error, spherical: {E_sph:.2e}")
+computed_solution = my_model.h_transport_problem.mobile.post_processing_solution
+E = f.errornorm(computed_solution, c_exact, "L2")
+print(f"L2 error: {E:.2e}")
 
 # plot exact solution and computed solution
-fig, axs = plt.subplots(
-    1,
-    4,
-    figsize=(15, 5),
-    sharey=True,
-)
-
+fig, axs = plt.subplots(1, 3, figsize=(15, 5))
 plt.sca(axs[0])
 plt.title("Exact solution")
 plt.xlabel("r")
-plt.ylabel("$\\theta$")
+plt.ylabel("z")
 CS1 = f.plot(c_exact, cmap="inferno")
-plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
 plt.sca(axs[1])
 plt.xlabel("r")
-plt.title("Cylindrical solution")
-CS2 = f.plot(cylindrical_solution, cmap="inferno")
-plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
-plt.sca(axs[2])
-plt.xlabel("r")
-plt.title("Spherical solution")
-CS3 = f.plot(spherical_solution, cmap="inferno")
-plt.gca().set_aspect(1 / (2 * np.pi)) # TODO: change this
+plt.title("Computed solution")
+CS2 = f.plot(computed_solution, cmap="inferno")
 
-plt.colorbar(CS1, ax=[axs[0]], shrink=0.8)
-plt.colorbar(CS2, ax=[axs[1]], shrink=0.8)
-plt.colorbar(CS3, ax=[axs[2]], shrink=0.8)
+plt.colorbar(CS2, ax=[axs[0], axs[1]], shrink=0.8)
 
-# axs[0].sharey(axs[1])
+axs[0].sharey(axs[1])
 plt.setp(axs[1].get_yticklabels(), visible=False)
-plt.setp(axs[2].get_yticklabels(), visible=False)
 
-for CS in [CS1, CS2, CS3]:
+for CS in [CS1, CS2]:
     CS.set_edgecolor("face")
 
 
@@ -188,9 +211,9 @@ def compute_arc_length(xs, ys):
 
 # define the profiles
 profiles = [
-    {"start": (0.0, 0.0), "end": (1.0, 2 * np.pi)},
-    {"start": (0.2, 2 * np.pi / 3), "end": (0.7, np.pi)},
-    {"start": (0.5, 1), "end": (0.8, 4)},
+    {"start": (1.0, 0.0), "end": (2.0, 1.0)},
+    {"start": (1.2, 0.8), "end": (1.7, 0.2)},
+    {"start": (1.2, 0.6), "end": (1.8, 0.8)},
 ]
 
 # plot the profiles on the right subplot
@@ -200,7 +223,7 @@ for i, profile in enumerate(profiles):
     plt.sca(axs[1])
     (l,) = plt.plot([start_x, end_x], [start_y, end_y])
 
-    plt.sca(axs[-1])
+    plt.sca(axs[2])
 
     points_x_exact = np.linspace(start_x, end_x, num=30)
     points_y_exact = np.linspace(start_y, end_y, num=30)
@@ -210,7 +233,7 @@ for i, profile in enumerate(profiles):
     points_x = np.linspace(start_x, end_x, num=100)
     points_y = np.linspace(start_y, end_y, num=100)
     arc_lengths = compute_arc_length(points_x, points_y)
-    computed_values = [cylindrical_solution(x, y) for x, y in zip(points_x, points_y)]
+    computed_values = [computed_solution(x, y) for x, y in zip(points_x, points_y)]
 
     (exact_line,) = plt.plot(
         arc_length_exact,
@@ -222,7 +245,7 @@ for i, profile in enumerate(profiles):
     )
     (computed_line,) = plt.plot(arc_lengths, computed_values, color=l.get_color())
 
-plt.sca(axs[-1])
+plt.sca(axs[2])
 plt.xlabel("Arc length")
 plt.ylabel("Solution")
 
@@ -234,6 +257,7 @@ legend_marker = mpl.lines.Line2D(
     linestyle="None",
     label="Exact",
 )
+
 legend_line = mpl.lines.Line2D([], [], color="black", label="Computed")
 plt.legend(
     [legend_marker, legend_line], [legend_marker.get_label(), legend_line.get_label()]
diff --git a/report/verification/soret.md b/report/verification/soret.md
new file mode 100644
index 0000000..184aeba
--- /dev/null
+++ b/report/verification/soret.md
@@ -0,0 +1 @@
+# Soret effect