Added ParameterAnalysis + InstructiveExample

examples/OneDimensionalEquilibrium-AnInstructiveExample.py

+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+import numpy as np
+from scipy.integrate import odeint
+import matplotlib.pyplot as plt
+
+# Define the parameters and boundary conditions
+Lambda2 = 8.553e-6
+z_A = -1
+z_C = 1
+y_L = 1/3
+
+# Define the equilibrium concentrations
+def y_C(phi, phi_prime):
+    return y_L * np.exp(-z_C * phi - 1/2 * Lambda2 * phi_prime**2)
+
+def y_A(phi, phi_prime):
+    return y_L * np.exp(-z_A * phi - 1/2 * Lambda2 * phi_prime**2)
+
+# Define the ODEs as a system of first-order ODEs
+def func(PHI, x):
+    phi, phi_prime = PHI
+    f = -1/Lambda2 * (z_C * y_C(phi, phi_prime) + z_A * y_A(phi, phi_prime))
+    return [phi_prime, f]
+
+# Initial conditions
+Y0 = [0, 1e-9]
+x_end = np.array([-0.95, -1.0, -1.05, -1.1])*1e-1
+
+# Storage for the results
+y_C_end, y_A_end, y_S_end = [], [], []
+phi_end, x_vec = [], []
+
+# Solve the ODEs for different values of x_end
+for x_ in x_end:
+    x = np.linspace(x_, 0, 1024)
+
+    Y = odeint(func, Y0, x)
+
+    Y_mirror = np.flip(Y, axis=0)
+    phi_end.append(Y_mirror[:, 0])
+    y_C_end.append(y_C(Y_mirror[:, 0], Y_mirror[:, 1]))
+    y_A_end.append(y_A(Y_mirror[:, 0], Y_mirror[:, 1]))
+    y_S_end.append(1-y_C_end[-1]-y_A_end[-1])
+    x_vec.append(x)
+
+# Visualize the results
+# create 2x1 subplots
+fig, axs = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(14-14/3, 12), constrained_layout=True)
+
+# clear subplots
+for ax in axs:
+    ax.remove()
+
+# add subfigure per subplot
+gridspec = axs[0].get_subplotspec().get_gridspec()
+subfigs = [fig.add_subfigure(gs) for gs in gridspec]
+
+for row, subfig in enumerate(subfigs):
+    subfig.suptitle(f'$L_x = $ {abs(round(x_end[row], 4))}', fontsize=22)
+
+    # create 1x2 subplots per subfig
+    axs = subfig.subplots(nrows=1, ncols=2)
+
+    for col, ax in enumerate(axs):
+        labelsize = 20
+        ax.set_ylim(-0.05, 1.05)
+        match col:
+            case 0: 
+                ax.plot(x_vec[row], phi_end[row])
+                ax.set_ylim(-0.5, phi_end[-1][0]+0.5)
+                ax.set_ylabel('$\\varphi$', fontsize=labelsize)
+            case 1: 
+                ax.plot(x_vec[row], y_C_end[row], label='$y_C$')
+                ax.plot(x_vec[row], y_A_end[row], label='$y_A$')
+                ax.plot(x_vec[row], y_S_end[row], label='$y_S$')
+                ax.legend()
+                ax.set_ylabel('$y_\\alpha$', fontsize=labelsize)
+        ax.set_xlim(min(x_end), 0)
+        ax.set_xlabel('x [-]', fontsize=labelsize)
+        ax.grid(True)
+        ax.tick_params(axis='x', labelrotation=45)
+fig.show()
+
+
+
+
+# Streamplot
+phi_range = np.arange(-0.1, 0.1, 0.01)  # Range for x-axis
+dphi_range = np.linspace(-5e-3, 5e-3, 10)  # Range for y-axis
+
+# Create a meshgrid with different sizes for x and y axes
+phi_values, dphi_values = np.meshgrid(phi_range, dphi_range)
+initial_conditions = np.array([phi_values.ravel(), dphi_values.ravel()]).T
+
+# Points at which the solution is requested
+x_domain = np.linspace(0, 1.0e-1, 1024)
+
+# Prepare storage for the results
+results = np.zeros((len(initial_conditions), len(x_domain), 2))
+
+# Numerical integration for each initial condition
+for i, Y0 in enumerate(initial_conditions):
+    Y = odeint(func, Y0, x_domain)
+    results[i, :, :] = Y
+
+# Reshape results for plotting
+Y_reshaped = results[:, :, 0].reshape(len(dphi_range), len(phi_range), len(x_domain))
+
+# Select a specific slice for plotting, e.g., the last time point -> last point of phi
+Y_slice = Y_reshaped[:, :, -1]
+
+# Compute derivatives for the streamplot
+U, V = np.gradient(Y_slice, dphi_range, phi_range)
+
+# Plotting
+plt.subplots(tight_layout=True)
+# lw = 5*Y_slice/Y_slice.max()
+strm = plt.streamplot(phi_values, dphi_values, U, V, color=U, linewidth=2, cmap='autumn', density=10)
+plt.colorbar(strm.lines)
+plt.xlim(-0.01, 0.01)
+plt.ylim(dphi_values.min(), dphi_values.max())
+plt.xlabel('$\\varphi$')
+plt.xticks(rotation='vertical')
+plt.ylabel("$\phi$")
+plt.show()
\ No newline at end of file

examples/ParameterAnalysis/Lambda.py

+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+
+# Define Parameter and buondary conditions
+phi_right = 0.0
+p_right = 0.0
+y_A_R, y_C_R = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+K = 'incompressible'
+a2 = 7.5412e-4
+number_cells = 1024*4
+refinement_style = 'hard_log'
+rtol = 1e-8
+max_iter = 10_000
+
+# Define left values of electric potential and values for lambda^2
+phi_left_vec = [-1.0, 1.0, 4.0, 8.0]
+Lambda2_vec = [8.553e-5, 8.553e-6, 8.553e-7]
+y_A, y_C, y_S, phi, p, x = [], [], [], [], [], []
+
+# Solve the system for different values of phi_left
+for phi_left_ in phi_left_vec:
+    y_A_bcs, y_C_bcs, y_S_bcs, phi_bcs, p_bcs, x_bcs = [], [], [], [], [], []
+    for l2 in Lambda2_vec:
+        y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_left_, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, l2, a2, number_cells, relax_param=0.05, x0=0, x1=1, refinement_style=refinement_style, return_type='Vector', max_iter=max_iter, rtol=rtol)
+        y_A_bcs.append(y_A_)
+        y_C_bcs.append(y_C_)
+        y_S_bcs.append(1 - y_A_ - y_C_)
+        phi_bcs.append(phi_)
+        p_bcs.append(p_)
+        x_bcs.append(x_)
+    y_A.append(y_A_bcs)
+    y_C.append(y_C_bcs)
+    y_S.append(y_S_bcs)
+    phi.append(phi_bcs)
+    p.append(p_bcs)
+    x.append(x_bcs)
+    
+# Visualize the results
+fig, axs = plt.subplots(2, 2)
+xlim = 0.08
+markers = ['--', '-', ':']
+colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
+for i in range(len(phi_left_vec)):
+    # Plot phi
+
+    match i:
+        case 0: 
+            i_ = 0
+            j_ = 0
+        case 1: 
+            i_ = 0
+            j_ = 1
+        case 2: 
+            i_ = 1
+            j_ = 0
+        case 3: 
+            i_ = 1
+            j_ = 1
+    axs[i_,j_].set_title(f'$\delta \\varphi$ = {phi_left_vec[i]}')#, fontsize=labelsize)
+    for j in range(len(Lambda2_vec)):
+        clr = colors[j]
+        axs[i_,j_].plot(x[i][j], y_A[i][j], markers[0], color=clr)#, lw=lw)
+        axs[i_,j_].plot(x[i][j], y_C[i][j], markers[1], color=clr)#, lw=lw)
+        axs[i_,j_].plot(x[i][j], y_S[i][j], markers[2], color=clr)#, lw=lw)
+    axs[i_,j_].grid()
+    axs[i_,j_].set_xlim(0,xlim)
+    axs[i_,j_].set_xlabel('$x$ [-]')#, fontsize=labelsize)
+    axs[i_,j_].set_ylabel('$y_\\alpha$ [-]')#, fontsize=labelsize)
+    axs[i_,j_].tick_params(axis='both')#, labelsize=labelsize)
+
+dummy, = axs[i_,j_].plot(2, y_A_R, color='grey', linestyle='--', label='$y_A$')
+dummy, = axs[i_,j_].plot(2, y_A_R, color='grey', linestyle='-', label='$y_C$')
+dummy, = axs[i_,j_].plot(2, y_A_R, color='grey', linestyle=':', label='$y_S$')
+for index, Lambda2_cur in enumerate(Lambda2_vec):
+    dummy, = axs[i_,j_].plot(10, y_A_R, color=colors[index], linestyle='-', label=f'$\lambda^2$ = {Lambda2_cur}')
+
+lines_labels = [ax.get_legend_handles_labels() for ax in fig.axes]
+lines, labels = [sum(lol, []) for lol in zip(*lines_labels)]
+
+# Finally, the legend (that maybe you'll customize differently)
+lgnd = fig.legend(lines, labels, bbox_to_anchor=(0.88,1.15), ncol=3)
+fig.tight_layout()
+fig.show()
\ No newline at end of file

examples/ParameterAnalysis/ScalingPhi.py

+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+
+# Define Parameter and buondary conditions
+phi_right = 0.0
+p_right = 0.0
+y_A_R, y_C_R = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+number_cells = 1024*4
+refinement_style = 'hard_log'
+rtol = 1e-8
+max_iter = 10_000
+
+phi_left_vector = [1,2,4,8]
+
+phi, x = [], []
+# Solve the system for different values of phi_left
+for phi_L in phi_left_vector:
+    y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_L, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=0.05, x0=0, x1=1, refinement_style=refinement_style, return_type='Vector', rtol=rtol, max_iter=max_iter)
+    phi.append(phi_)
+    x.append(x_)
+    
+# Normalize the electric potentials
+phi_scaled = []
+for i, phi_L in enumerate(phi_left_vector):
+    phi_scaled.append(phi[i]/phi_L)
+    
+# Visualize the results
+# Non-normalized electric potential
+for i, phi_L in enumerate(phi_left_vector):
+    plt.plot(x[i], phi[i], label='$\delta \\varphi$ = ' + str(phi_L))
+plt.legend()
+plt.xlabel('$x$ [-]')
+plt.ylabel('$\\varphi$ [-]')
+plt.ylim(0, max(phi_left_vector) + 0.1)
+plt.xlim(0,0.05)
+plt.grid()
+plt.tight_layout()
+plt.show()
+
+# Normalized electric potential
+for i, phi_L in enumerate(phi_left_vector):
+    plt.plot(x[i], phi_scaled[i], label='$\delta \\varphi$ = ' + str(phi_L))
+plt.legend()
+plt.xlabel('$x$ [-]')
+plt.ylabel('$\\varphi/\delta \\varphi$ [-]')
+plt.ylim(0, 1.05)
+plt.xlim(0,0.05)
+plt.grid()
+plt.tight_layout()
+plt.show()
\ No newline at end of file

examples/ParameterAnalysis/a.py

+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Define Parameter and buondary conditions
+phi_right = 0.0
+p_right = 0.0
+y_A_R, y_C_R = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+number_cells = 1024*4
+refinement_style = 'hard_log'
+rtol = 1e-8
+max_iter = 10_000
+
+phi_left_vec = [0.1, 1.0, 4.0, 8.0]
+a2_vec = [7.5412e-3, 7.5412e-4, 7.5412e-5]
+y_A, y_C, y_S, phi, p, x = [], [], [], [], [], []
+for phi_left_ in phi_left_vec:
+    y_A_bcs, y_C_bcs, y_S_bcs, phi_bcs, p_bcs, x_bcs = [], [], [], [], [], []
+    for a2 in a2_vec:
+        y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_left_, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=0.05, x0=0, x1=1, refinement_style=refinement_style, return_type='Vector', max_iter=max_iter, rtol=rtol)
+        y_A_bcs.append(y_A_)
+        y_C_bcs.append(y_C_)
+        y_S_bcs.append(1 - y_A_ - y_C_)
+        phi_bcs.append(phi_)
+        p_bcs.append(p_)
+        x_bcs.append(x_)
+    y_A.append(y_A_bcs)
+    y_C.append(y_C_bcs)
+    y_S.append(y_S_bcs)
+    phi.append(phi_bcs)
+    p.append(p_bcs)
+    x.append(x_bcs)
+
+p = np.array(p)
+
+plt.figure()
+for j in range(len(a2_vec)):
+    plt.plot(x[-1][j], p[-1][j], label='$a^2$ = {}'.format(a2_vec[j]))
+plt.grid()
+plt.ylim(0.0, np.max(p[-1][0]))
+plt.xlim(0,0.08)
+plt.xlabel('$x$ [-]')
+plt.ylabel('$p$ [-]')
+plt.tick_params(axis='both')
+plt.legend(loc='best')
+plt.tight_layout()
+plt.show()
\ No newline at end of file

src/Eq04.py

@@ -8,12 +8,48 @@ from mpi4py import MPI
 from dolfinx import mesh, fem, log
 from dolfinx.fem.petsc import NonlinearProblem
 from dolfinx.nls.petsc import NewtonSolver
-from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Cell
+from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
 from basix.ufl import element, mixed_element
-from petsc4py import PETSc
 import matplotlib.pyplot as plt
 
-from oneD_refined_mesh import create_refined_mesh
+# Define mesh
+def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
+    '''
+    Creates a one-dimensional mesh with a refined region at the left boundary
+
+    Parameters
+    ----------
+    refinement_style : str
+        How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
+    number_cells : int
+        Number of cells in the mesh
+
+    Returns
+    -------
+    Mesh
+        One-dimensional mesh, ready for use in FEniCSx
+    '''
+    if refinement_style == 'log':
+        coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
+    elif refinement_style == 'hard_log':
+        coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
+        coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
+        coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+    elif refinement_style == 'hard_hard_log':
+        coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
+        coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
+        coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+    num_vertices = len(coordinates_np)
+    num_cells = num_vertices - 1
+    cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
+    gdim = 1
+    shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
+    degree = 1
+    domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
+    coordinates_np_ = []
+    [coordinates_np_.append([coord]) for coord in coordinates_np]
+    msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
+    return msh
 
 def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
     '''
@@ -26,6 +62,8 @@ def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
 
     with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.
 
+    ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
     Parameters
     ----------
     phi_left : float

src/oneD_refined_mesh.py

-'''
-Jan Habscheid
-Jan.Habscheid@rwth-aachen.de
-'''
-
-import numpy as np
-from mpi4py import MPI
-from dolfinx import mesh
-from basix.ufl import element
-from ufl import Mesh
-
-
-# Define mesh
-def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
-    '''
-    Creates a one-dimensional mesh with a refined region at the left boundary
-
-    Parameters
-    ----------
-    refinement_style : str
-        How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
-    number_cells : int
-        Number of cells in the mesh
-
-    Returns
-    -------
-    Mesh
-        One-dimensional mesh, ready for use in FEniCSx
-    '''
-    if refinement_style == 'log':
-        coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
-    elif refinement_style == 'hard_log':
-        coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
-        coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
-        coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
-    elif refinement_style == 'hard_hard_log':
-        coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
-        coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
-        coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
-    num_vertices = len(coordinates_np)
-    num_cells = num_vertices - 1
-    cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
-    gdim = 1
-    shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
-    degree = 1
-    domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
-    coordinates_np_ = []
-    [coordinates_np_.append([coord]) for coord in coordinates_np]
-    msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
-    return msh
\ No newline at end of file