Added Solvation with numerical runs + diode some tests

examples/ReproducableCode/DoubleLayerCapacity/Solvation.py

@@ -16,6 +16,7 @@ src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../..',
 sys.path.insert(0, src_path)
 
 from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+from Eq04 import solve_System_4eq
 
 
 # Remove the src directory from sys.path after import
@@ -35,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
@@ -50,32 +51,74 @@ p_R = 0
 
 # Solver settings
 number_cells = 1024
-relax_param = 0.03
-p_right = 0
+relax_param = 0.3
+rtol = 1e-4 # ! Change to 1e-8
+refinement_style = 'hard_log'
+max_iter = 10_000
+return_type = 'Vector'
 
 # phi^L domain
-Vol_start = -1.0
+Vol_start = 0 # ! -1.0
 Volt_end = 1.0
-n_Volts = 5_000
+n_Volts = 8
 
 Phi_Pot_Diff_dim = np.linspace(Vol_start, Volt_end, n_Volts)
 Phi_Pot_Diff_dimless = Phi_Pot_Diff_dim * e0/(k*T)
 
-
+Phi_Pot_Diff_dim_PB = np.linspace(0., 0.25, n_Volts)
+Phi_Pot_Diff_dimless_PB = Phi_Pot_Diff_dim_PB * e0/(k*T)
 
 
 # kappa
-kappa_vec = np.array([0, 5, 10, 15])
-
-Q_DL_dim_ = []
-Q_DL_dimless_ = []
-for kappa in kappa_vec:
-    y_A_R, y_C_R = Molarity / nR_mol, Molarity / nR_mol
-    y_N_R = 1 - y_A_R - y_C_R
-    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, kappa))
-    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, nR_m, e0, LR, kappa))
-    
-C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
+kappa_vec = [0, 5]#, 10, 15])
+kappa_legend = ['PB', 0, 5]#, 10, 15]
+# Solution vectors
+y_A, y_C, y_S, phi, p, x = [], [], [], [], [], []
+
+# Extra run for Poisson-Boltzmann
+phi_PB, y_A_PB, y_C_PB, y_S_PB, p_PB,x_PB = [], [], [], [], [], []
+for i, phi_bcs in enumerate(Phi_Pot_Diff_dimless_PB):
+    y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_R, p_R, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=0, PoissonBoltzmann=True, refinement_style='hard_log', rtol=rtol, max_iter=max_iter, relax_param=0.01, return_type=return_type)
+    y_A_PB.append(y_A_)
+    y_C_PB.append(y_C_)
+    y_S_PB.append(1 - y_A_ - y_C_)
+    phi_PB.append(phi_)
+    p_PB.append(p_)
+    x_PB.append(x_)
+
+for kappa_ in kappa_vec:
+    phi_inner, y_A_inner, y_C_inner, y_S_inner, p_inner,x_inner = [], [], [], [], [], []
+    for i, phi_bcs in enumerate(Phi_Pot_Diff_dimless):
+        y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_R, p_R, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa_, refinement_style='hard_log', rtol=rtol, max_iter=max_iter, return_type=return_type)
+        y_A_inner.append(y_A_)
+        y_C_inner.append(y_C_)
+        y_S_inner.append(1 - y_A_ - y_C_)
+        phi_inner.append(phi_)
+        p_inner.append(p_)
+        x_inner.append(x_)
+    y_A.append(y_A_inner)
+    y_C.append(y_C_inner)
+    y_S.append(y_S_inner)
+    phi.append(phi_inner)
+    p.append(p_inner)
+    x.append(x_inner)
+
+        
+Q = []
+Q_DL_dimless_PB = []
+for j in range(len(Phi_Pot_Diff_dimless_PB)):
+    Q_DL_dimless_PB.append(Q_num_(y_A_PB[j], y_C_PB[j], n(p_PB[j], K), x_PB[j]))
+Q_DL_dimless_PB= np.array(Q_DL_dimless_PB)
+for i in range(len(kappa_vec)):
+    Q_inner = []
+    for j in range(len(Phi_Pot_Diff_dimless)):
+        Q_inner.append(Q_num_(y_A[i][j], y_C[i][j], n(p[i][j], K), x[i][j]))
+    Q.append(Q_inner)
+Q_DL_dimless_ = np.array(Q)
+
+    
+# C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
+C_DL_dimless_PB = C_dl(Q_DL_dimless_PB, Phi_Pot_Diff_dimless_PB)
 C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
 
 
@@ -89,17 +132,19 @@ color_dim = 'tab:red'
 ax = fig.add_subplot(111, label="1")
 ax2 = fig.add_subplot(111, label="2", frame_on=False)
 
-# Plot dimensional data
-for i, q_dl in enumerate(Q_DL_dim_):
-    ax.plot(Phi_Pot_Diff_dim, q_dl, color=colors[i], label=f'$\kappa$: {kappa_vec[i]}')
-ax.grid()
-ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
-ax.set_ylabel('$Q$ [\u03bc$F/cm^3]$', color=color_dim)
-ax.tick_params(axis='x', colors=color_dim)
-ax.tick_params(axis='y', colors=color_dim)
+# # Plot dimensional data
+# for i, q_dl in enumerate(Q_DL_dim_):
+#     ax.plot(Phi_Pot_Diff_dim, q_dl, color=colors[i], label=f'$\kappa$: {kappa_vec[i]}')
+# ax.grid()
+# ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+# ax.set_ylabel('$Q$ [\u03bc$F/cm^3]$', color=color_dim)
+# ax.tick_params(axis='x', colors=color_dim)
+# ax.tick_params(axis='y', colors=color_dim)
+
 
+ax2.plot(Phi_Pot_Diff_dimless_PB, Q_DL_dimless_PB, color=colors[0], label=f'$\kappa$: {kappa_legend[0]}')
 for i, q_dl in enumerate(Q_DL_dimless_):
-    ax2.plot(Phi_Pot_Diff_dimless, q_dl, color=colors[i])
+    ax2.plot(Phi_Pot_Diff_dimless, q_dl, color=colors[i+1], label=f'$\kappa$: {kappa_legend[i+1]}')
 ax2.xaxis.tick_top()
 ax2.yaxis.tick_right()
 ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
@@ -115,6 +160,7 @@ fig.show()
 
 
 # Double Layer Capacity
+Phi_pot_center_array_dimless_PB = Phi_pot_center(Phi_Pot_Diff_dimless_PB)
 Phi_pot_center_array_dimless = Phi_pot_center(Phi_Pot_Diff_dimless)
 Phi_pot_center_array_dim = Phi_pot_center(Phi_Pot_Diff_dim)
 fig = plt.figure()
@@ -123,17 +169,18 @@ color_dim = 'tab:red'
 ax = fig.add_subplot(111, label="1")
 ax2 = fig.add_subplot(111, label="2", frame_on=False)
 
-# Plot dimensional data
-for i, c_dl in enumerate(C_DL_dim):
-    ax.plot(Phi_pot_center_array_dim, c_dl, color=colors[i], label=f'$\kappa$: {kappa_vec[i]}')
-ax.grid()
-ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
-ax.set_ylabel('$C_{dl}$ [\u03bc$F/cm^2]$', color=color_dim)
-ax.tick_params(axis='x', colors=color_dim)
-ax.tick_params(axis='y', colors=color_dim)
+# # Plot dimensional data
+# for i, c_dl in enumerate(C_DL_dim):
+#     ax.plot(Phi_pot_center_array_dim, c_dl, color=colors[i], label=f'$\kappa$: {kappa_vec[i]}')
+# ax.grid()
+# ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+# ax.set_ylabel('$C_{dl}$ [\u03bc$F/cm^2]$', color=color_dim)
+# ax.tick_params(axis='x', colors=color_dim)
+# ax.tick_params(axis='y', colors=color_dim)
 
+ax2.plot(Phi_pot_center_array_dimless_PB, C_DL_dimless_PB, color=colors[0], label=f'$\kappa$: {kappa_legend[0]}')
 for i, c_dl in enumerate(C_DL_dimless):
-    ax2.plot(Phi_pot_center_array_dimless, c_dl, color=colors[i])
+    ax2.plot(Phi_pot_center_array_dimless, c_dl, color=colors[i+1], label=f'$\kappa$: {kappa_legend[i+1]}')
 # ax2.grid()
 ax2.xaxis.tick_top()
 ax2.yaxis.tick_right()
@@ -148,5 +195,5 @@ fig.legend()
 fig.tight_layout()
 fig.show()
 
-# Save the results
-np.savez('../../Data/DoubleLayerCapacity/Solvation.npz', Lambda2=Lambda2, a2=a2, K=K, kappa_vec=kappa_vec, z_A=z_A, z_C=z_C, phi_R=phi_R, p_R=p_R, Vol_start=Vol_start, Volt_end=Volt_end, n_Volts=n_Volts, Phi_pot_center_array_dim=Phi_pot_center_array_dim, Phi_pot_center_array_dimless=Phi_pot_center_array_dimless, Phi_Pot_Diff_dim=Phi_Pot_Diff_dim, Phi_Pot_Diff_dimless=Phi_Pot_Diff_dimless, C_DL_dim=C_DL_dim, C_DL_dimless=C_DL_dimless, Q_DL_dim_=Q_DL_dim_, Q_DL_dimless_=Q_DL_dimless_, Molarity=Molarity)
+# # Save the results
+# np.savez('../../Data/DoubleLayerCapacity/Solvation.npz', Lambda2=Lambda2, a2=a2, K=K, kappa_vec=kappa_vec, z_A=z_A, z_C=z_C, phi_R=phi_R, p_R=p_R, Vol_start=Vol_start, Volt_end=Volt_end, n_Volts=n_Volts, Phi_pot_center_array_dim=Phi_pot_center_array_dim, Phi_pot_center_array_dimless=Phi_pot_center_array_dimless, Phi_Pot_Diff_dim=Phi_Pot_Diff_dim, Phi_Pot_Diff_dimless=Phi_Pot_Diff_dimless, C_DL_dim=C_DL_dim, C_DL_dimless=C_DL_dimless, Q_DL_dim_=Q_DL_dim_, Q_DL_dimless_=Q_DL_dimless_, Molarity=Molarity)

src/ElectrolyticDiode.py

@@ -128,6 +128,8 @@ def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C
 
     # Define Finite Elements
     CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+    # from ufl import FiniteElement
+    # CG1_elem = FiniteElement("Lagrange", msh.ufl_cell(), 1)
 
     # Define Mixed Function Space
     W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
@@ -188,6 +190,8 @@ def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C
 
     # Collect boundary conditions
     bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_top_y_A, bc_bottom_y_C, bc_top_y_C]
+    # bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_bottom_y_C]
+
 
     # def p_bottom_(x):
     #     return np.full_like(x[1], 0)
@@ -206,25 +210,54 @@ def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C
 
     # Variational formulation
     if K == 'incompressible':
-        def nF(y_A, y_C):
-            # n = const = 1
-            return (z_C * y_C + z_A * y_A)
-
-        def J_A(y_A, y_C, phi, p):
+        # def nF(y_A, y_C):
+        #     return (z_C * y_C + z_A * y_A) # n = 1
+    
+        # 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(ln(y_A) + a2 * solvation * p - ln(1-y_A-y_C) + z_A * phi), grad(v_A)) * dx
+        #     + inner(grad(ln(y_C) + a2 * solvation * p- ln(1-y_A-y_C) + z_C * phi), grad(v_C)) * dx
+        # )
+        # def nF(y_A, y_C):
+        #     # n = const = 1
+        #     return (z_C * y_C + z_A * y_A)
+
+        # def J_A(y_A, y_C, phi, p):
         #     g_A = (solvation + 1) * a2 * (p - 1) # g_Aref, but constant and take gradient
         #     mu_A = g_A + ln(y_A)
         #     g_N = a2 * (p - 1) # solvation_N = 0, g_Sref, but constant and take gradient
         #     mu_N = g_N + ln(1 - y_A - y_C)
         #     return grad(mu_A - mu_N + z_A * phi)
-            return grad(ln(y_A) - ln(1 - y_A - y_C) + z_A * phi)
+        #     # return grad(ln(y_A) - ln(1 - y_A - y_C) + z_A * phi)
         
-        def J_C(y_A, y_C, phi, p):
+        # def J_C(y_A, y_C, phi, p):
         #     g_C = (solvation + 1) * a2 * (p - 1) # g_Cref, but constant and take gradient
         #     mu_C = g_C + ln(y_C)
         #     g_N = a2 * (p - 1)
         #     mu_N = g_N + ln(1 - y_A - y_C)
         #     return grad(mu_C - mu_N + z_C * phi)
-            return grad(ln(y_C) - ln(1 - y_A - y_C) + z_C * phi)
+        #     # return grad(ln(y_C) - ln(1 - y_A - y_C) + z_C * phi)
+
+        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 grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
+        
+        def J_C(y_A, y_C, phi, p):
+            return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
+
+        # if PoissonBoltzmann:
+        #     def J_A(y_A, y_C, phi, p):
+        #         return grad(ln(y_A) + z_A * phi)
+        #     def J_C(y_A, y_C, phi, p):
+        #         return grad(ln(y_C) + z_C * phi)
 
         A = (
             inner(grad(phi), grad(v_1)) * dx
@@ -305,23 +338,24 @@ def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C
         raise ValueError('Invalid return_type')
 
 if __name__ == '__main__':
-    phi_bias = 2
-    Bias_type = 'ForwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
-    g_phi = 5
-    y_fixed = 0.01
+    phi_bias = 10#10
+    Bias_type = 'NoBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
+    g_phi = 5 #0.5#5
+    y_fixed = 0.01#0.01
     z_A = -1.0
     z_C = 1.0
     K = 'incompressible'
     Lambda2 = 8.553e-2 # ! Change back to 1e-6
     # g_phi *= np.sqrt(Lambda2) # ! Unsure about this scaling
+    # g_phi *= Lambda2
     a2 = 7.5412e-4
-    number_cells = [20, 100]
+    number_cells = [20,128]#[20, 100]
     Lx = 2
     Ly = 10
     x0 = np.array([0, 0])
     x1 = np.array([Lx, Ly])
     refinement_style = 'uniform'
-    solvation = 0
+    solvation = 3
     PoissonBoltzmann = False
     rtol = 1e-3 # ToDo: Change back to 1e-8, currently just for testing
     relax_param = 0.05