Eq02 now really eq02, before by accident Eq04

examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py

@@ -36,7 +36,7 @@ NA = 6.022e+23 # [1/mol] - Avogadro constant
 nR_mol = 55
 nR_m = nR_mol * NA * 1/(1e-3)# [1/m^3]
 pR = 1.01325 * 1e+5 # [Pa]
-LR = 20e-8
+LR = 20e-9
 chi = 80 # [-]
 
 # Parameter and bcs for the electrolyte
@@ -54,12 +54,13 @@ p_right = 0
 number_cells = 1024
 relax_param = 0.03
 p_right = 0
+rtol = 1e-4 # ! Change back to 1e-8
 
 
 # phi^L domain
-Vol_start = 0
+Vol_start = 0.1 # ! Change back to 0
 Volt_end = 0.75
-n_Volts = 5#0
+n_Volts = 25#0
 
 phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 
@@ -68,7 +69,7 @@ phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 # Solution vectors
 y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
 for i, phi_bcs in enumerate(phi_left):
-    y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-10, max_iter=2500, return_type='Vector', relax_param=relax_param)
+    y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=rtol, max_iter=2500, return_type='Vector', relax_param=relax_param)
     y_S_ = 1 - y_A_ - y_C_
     y_A_num.append(y_A_)
     y_C_num.append(y_C_)
@@ -82,9 +83,9 @@ for j in range(len(phi_left)):
     Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
 Q_num = np.array(Q_num)
 
-# dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
-# C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
-# C_DL_num = np.array(C_DL_num)
+dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
+C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
+C_DL_num = np.array(C_DL_num)
 C_dl_num = C_dl(Q_num, phi_left)
 
 
@@ -115,4 +116,8 @@ plt.legend()
 plt.xlabel('$\delta \\varphi$ [-]')
 plt.ylabel('$C_{dl,num} - C_{dl,ana} [-]$')
 plt.tight_layout()
-plt.show()
\ No newline at end of file
+plt.show()
+
+print('phi_left:', phi_left)
+print('Q_num:', Q_num)
+print('Q_ana:', Q_ana)
\ No newline at end of file

examples/ReproducableCode/ValidationSimplifiedSystem.py

+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This script is used to validate the simplification of the system of four equations to a system of two equations in the one-dimensional equilibrium case.
+'''
+
+# 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
+from Eq02 import solve_System_2eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Further imports
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+# Define the parameters and boundary conditions
+phi_left = 4.0
+phi_right = 0.0
+p_right = 0.0
+y_A_R = 0.01
+y_C_R = 0.01
+z_A = -1.0
+z_C = 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+solvation = 15
+number_cells = 1024#*4
+refinement_style = 'log'
+rtol = 1e-8
+
+# solve the complete system
+y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# solve the simplified system
+y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# Evaluate the difference
+plt.figure()
+# plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+# plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+plt.grid()
+plt.xlim(0, 0.05)
+plt.legend()
+plt.show()
+
+plt.figure()
+# plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+# plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+plt.grid()
+plt.xlim(0, 0.05)
+plt.legend()
+plt.show()
+
+plt.figure()
+# plt.plot(x_4eq, phi_4eq, label='phi_4eq')
+# plt.plot(x_2eq, phi_2eq, label='phi_2eq')
+plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
+plt.grid()
+plt.xlim(0, 0.05)
+plt.legend()
+plt.show()
+
+plt.figure()
+# plt.plot(x_4eq, p_4eq, label='p_4eq')
+# plt.plot(x_2eq, p_2eq, label='p_2eq')
+plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
+plt.grid()
+plt.xlim(0, 0.05)
+plt.legend()
+plt.show()
\ No newline at end of file

src/Eq02.py

@@ -8,7 +8,7 @@ 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, Mesh
+from ufl import TestFunctions, split, dot, grad, dx, inner, Mesh, exp
 from basix.ufl import element, mixed_element
 import matplotlib.pyplot as plt
 
@@ -51,7 +51,7 @@ def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
     msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
     return msh
 
-def solve_System_2eq(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):
+def solve_System_2eq(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, 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):
     '''
     Solve the simplified dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
 
@@ -96,8 +96,6 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
         Number of cells in the mesh
     solvation : float, optional
         solvation number, by default 0
-    PoissonBoltzmann : bool, optional
-        Solve classical Nernst-Planck model with the use of the Poisson-Boltzmann formulation if True, else solve the presented model by Dreyer, Guhlke, Müller, by default False
     relax_param : float, optional
         Relaxation parameter for the Newton solver
         xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
@@ -122,6 +120,8 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
         Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
         If return_type is 'Vector', the solution is returned as numpy arrays
     '''
+    if return_type == 'Scalar':
+        raise NotImplementedError('Scalar return type is not implemented yet')
     # Define boundaries of the domain
     x0 = 0
     x1 = 1
@@ -143,20 +143,18 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
     CG1_elem = element('Lagrange', msh.basix_cell(), 1)
 
     # Define Mixed Function Space
-    W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
+    W_elem = mixed_element([CG1_elem, CG1_elem])#, CG1_elem, CG1_elem])
     W = fem.functionspace(msh, W_elem)
 
     # Define Trial- and Testfunctions
     u = fem.Function(W)
-    y_A, y_C, phi, p = split(u)
-    (v_A, v_C, v_1, v_2) = TestFunctions(W)
+    phi, p = split(u)
+    (v_1, v_2) = TestFunctions(W)
 
     # Collapse function space for bcs
     W0, _ = W.sub(0).collapse()
     W1, _ = W.sub(1).collapse()
-    W2, _ = W.sub(2).collapse()
-    W3, _ = W.sub(3).collapse()
-
+    
     # Define boundary conditions values
     def phi_left_(x):
         return np.full_like(x[0], phi_left)
@@ -164,42 +162,36 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
         return np.full_like(x[0], phi_right)
     def p_right_(x):
         return np.full_like(x[0], p_right)
-    def y_A_right_(x):
-        return np.full_like(x[0], y_A_R)
-    def y_C_right_(x):
-        return np.full_like(x[0], y_C_R)
-
+    
     # Interpolate bcs functions
-    phi_left_bcs = fem.Function(W2)
+    phi_left_bcs = fem.Function(W0)
     phi_left_bcs.interpolate(phi_left_)
-    phi_right_bcs = fem.Function(W2)
+    phi_right_bcs = fem.Function(W0)
     phi_right_bcs.interpolate(phi_right_)
-    p_right_bcs = fem.Function(W3)
+    p_right_bcs = fem.Function(W1)
     p_right_bcs.interpolate(p_right_)
-    y_A_right_bcs = fem.Function(W0)
-    y_A_right_bcs.interpolate(y_A_right_)
-    y_C_right_bcs = fem.Function(W1)
-    y_C_right_bcs.interpolate(y_C_right_)
-
+    
     # Identify dofs for boundary conditions
     # Define boundary conditions
-    facet_left_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Left)
-    facet_right_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Right)
-    bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(2))
-    bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(2))
-
-    facet_right_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Right)
-    bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(3))
-
+    facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
     facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
-    bc_right_y_A = fem.dirichletbc(y_A_right_bcs, facet_right_dofs, W.sub(0))
+    bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
+    bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
 
     facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
-    bc_right_y_C = fem.dirichletbc(y_C_right_bcs, facet_right_dofs, W.sub(1))
+    bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
 
+    
     # Combine boundary conditions into list
-    bcs = [bc_left_phi, bc_right_phi, bc_right_p, bc_right_y_A, bc_right_y_C]
-        
+    bcs = [bc_left_phi, bc_right_phi, bc_right_p]
+
+    def y_A(phi, p):
+        D_A = y_A_R / exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
+        return D_A * exp(-(solvation + 1) * a2 * p - z_A * phi)
+    
+    def y_C(phi, p):
+        D_C = y_C_R / exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
+        return D_C * exp(-(solvation + 1) * a2 * p - z_C * phi)
     
 
     # Define variational problem
@@ -207,30 +199,6 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
         # total free charge density
         def nF(y_A, y_C):
             return (z_C * y_C + z_A * y_A)
-        
-        # Diffusion fluxes for species A and C
-        def J_A(y_A, y_C, phi, p):
-            return ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi
-        
-        def J_C(y_A, y_C, phi, p):
-            return ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi
-        
-        # Variational Form
-        A = (
-            inner(grad(phi), grad(v_1)) * dx
-            - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
-        ) + (
-            inner(grad(p), grad(v_2)) * dx
-            + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
-        ) + (
-            inner(grad(J_A(y_A, y_C, phi, p)), grad(v_A)) * dx
-            + inner(grad(J_C(y_A, y_C, phi, p)), grad(v_C)) * dx
-        )
-        if PoissonBoltzmann:
-            A += (
-                inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_A)) * dx
-                + inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_C)) * dx
-            )
     else: 
         # total number density
         def n(p):
@@ -239,37 +207,16 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
         # total free charge density
         def nF(y_A, y_C, p):
             return (z_C * y_C + z_A * y_A) * n(p)
-        
-        # Diffusion fluxes for species A and C
-        def J_A(y_A, y_C, phi, p):
-            return ln(y_A) + a2 * (solvation + 1) * K * ln(1 + 1/K * (p-1)) + z_A * phi
-        
-        def J_C(y_A, y_C, phi, p):
-            return ln(y_C) + a2 * (solvation + 1)* K * ln(1 + 1/K * (p-1)) + z_C * phi
-
-        A = (
-            inner(grad(phi), grad(v_1)) * dx
-            - 1 / Lambda2 * nF(y_A, y_C, p) * v_1 * dx
-        ) + (
-            inner(grad(p), grad(v_2)) * dx
-            + 1 / a2 * nF(y_A, y_C, p) * dot(grad(phi), grad(v_2)) * dx
-        ) + (
-            inner(grad(J_A(y_A, y_C, phi, p)), grad(v_A)) * dx
-            + inner(grad(J_C(y_A, y_C, phi, p)), grad(v_C)) * dx
-        )
+            # Variational Form
+    A = (
+        inner(grad(phi), grad(v_1)) * dx
+        - 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p)) * v_1 * dx
+    ) + (
+        inner(grad(p), grad(v_2)) * dx
+        + 1 / a2 * nF(y_A(phi, p), y_C(phi, p)) * dot(grad(phi), grad(v_2)) * dx
+    )
     F = A
 
-    # Initialize initial guess for u
-    y_C_init = fem.Function(W1)
-    y_A_init = fem.Function(W0)
-    y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_R))
-    y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_R))
-
-    with u.vector.localForm() as u_loc:
-        u_loc.set(0)
-    u.sub(0).interpolate(y_A_init)
-    u.sub(1).interpolate(y_C_init)
-
     # Define Nonlinear Problem
     problem = NonlinearProblem(F, u, bcs=bcs)
 
@@ -294,20 +241,22 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
     print(f"Number of interations: {n:d}")
 
     # Split the mixed function space into the individual components    
-    y_A, y_C, phi, p = u.split()
+    phi, p = u.split()
     
     # Return the solution
     if return_type=='Vector':
         x_vals = np.array(msh.geometry.x[:,0])
-        y_A_vals = np.array(u.sub(0).collapse().x.array)
-        y_C_vals = np.array(u.sub(1).collapse().x.array)
-        phi_vals = np.array(u.sub(2).collapse().x.array)
-        p_vals = np.array(u.sub(3).collapse().x.array)
+        phi_vals = np.array(u.sub(0).collapse().x.array)
+        p_vals = np.array(u.sub(1).collapse().x.array)
+
+        # Calculate the atomic fractions
+        D_A = y_A_R / np.exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
+        y_A_vals = D_A * np.exp(-(solvation + 1) * a2 * p_vals - z_A * phi_vals)
+    
+        D_C = y_C_R / np.exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
+        y_C_vals = D_C * np.exp(-(solvation + 1) * a2 * p_vals - z_C * phi_vals)
         
         return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
-    elif return_type=='Scalar':
-        return y_A, y_C, phi, p, msh
-    
     
 if __name__ == '__main__':
     # Define the parameters

src/Helper_DoubleLayerCapacity.py

@@ -157,10 +157,15 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
         D_C = y_C_R / (np.exp(-(solvation+1)*a2*p_R-z_C*phi_R))
         D_N = y_N_R / (np.exp(-(solvation+1)*a2*p_R-z_N*phi_R))
         E_p = p_R
-        CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - (solvation+1) * E_p) + D_C * np.exp(-z_C * phi_L - (solvation+1) * E_p) + D_N * np.exp(-z_N * phi_L - (solvation+1) * E_p))
-        CLambda_R = np.log(D_A * np.exp(-z_A * phi_R + E_p * a2) + D_C * np.exp(-z_C * phi_R + E_p * a2) + D_N * np.exp(-z_N * phi_R + E_p * a2))
+
+        CLambda_L = np.log(D_A * np.exp(-z_A * phi_L) + D_C * np.exp(-z_C * phi_L) + D_N * np.exp(-z_N * phi_L))
+        CLambda_R = np.log(D_A * np.exp(-z_A * phi_R) + D_C * np.exp(-z_C * phi_R) + D_N * np.exp(-z_N * phi_R))
+
+        # dx_Phi_L = np.sqrt((CLambda_L - (solvation+1) * a2 * E_p) * 2 / (solvation + 1))
+        # dx_Phi_R = np.sqrt((CLambda_R - (solvation+1) * a2 * E_p) * 2 / (solvation + 1))
+        # Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
         Lambda = np.sqrt(Lambda2)
-        Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(solvation+1)) * (np.sqrt(CLambda_L) - np.sqrt(CLambda_R))
+        Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(solvation+1)) * (np.sqrt(CLambda_L + E_p * a2) - np.sqrt(CLambda_R + E_p * a2))
     else:
         C_A = y_A_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_A*phi_R))
         C_C = y_C_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_C*phi_R))