diff --git a/examples/OneDimensionalEquilibrium-AnInstructiveExample.py b/examples/OneDimensionalEquilibrium-AnInstructiveExample.py
new file mode 100644
index 0000000000000000000000000000000000000000..aebf81f89794ac843f5bfd72d9f95b8b85c3ac9e
--- /dev/null
+++ b/examples/OneDimensionalEquilibrium-AnInstructiveExample.py
@@ -0,0 +1,129 @@
+'''
+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
diff --git a/examples/ParameterAnalysis/Lambda.py b/examples/ParameterAnalysis/Lambda.py
new file mode 100644
index 0000000000000000000000000000000000000000..b5fc5f31cc21d67d8c3243f59c584da7b95f884d
--- /dev/null
+++ b/examples/ParameterAnalysis/Lambda.py
@@ -0,0 +1,102 @@
+'''
+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
diff --git a/examples/ParameterAnalysis/ScalingPhi.py b/examples/ParameterAnalysis/ScalingPhi.py
new file mode 100644
index 0000000000000000000000000000000000000000..95da151770d1b49a02208487a8fb7359778b4d23
--- /dev/null
+++ b/examples/ParameterAnalysis/ScalingPhi.py
@@ -0,0 +1,72 @@
+'''
+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
diff --git a/examples/ParameterAnalysis/a.py b/examples/ParameterAnalysis/a.py
new file mode 100644
index 0000000000000000000000000000000000000000..042e3b48f91e467a598b8f23d4443925d8886292
--- /dev/null
+++ b/examples/ParameterAnalysis/a.py
@@ -0,0 +1,68 @@
+'''
+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
diff --git a/src/Eq04.py b/src/Eq04.py
index bed54623f96042dd885e597cb7af934b230bf87d..37c292f29d4fd86cc332aebb57468cad0f463fa1 100644
--- a/src/Eq04.py
+++ b/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
diff --git a/src/__pycache__/Eq04.cpython-312.pyc b/src/__pycache__/Eq04.cpython-312.pyc
index d115d2ed6ed42485613169199be3432cb7a86dc5..e94b79d647eeb82bdf59ecca6b780f96647b96e5 100644
Binary files a/src/__pycache__/Eq04.cpython-312.pyc and b/src/__pycache__/Eq04.cpython-312.pyc differ
diff --git a/src/oneD_refined_mesh.py b/src/oneD_refined_mesh.py
deleted file mode 100644
index d34bca8c0ee31d48e0ef9a40e3c0ee31b8f1e405..0000000000000000000000000000000000000000
--- a/src/oneD_refined_mesh.py
+++ /dev/null
@@ -1,50 +0,0 @@
-'''
-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