Adapted Lambda2 from 8.5e-6 to 3.4e-6 for better validation, still proplems...

Adapted Lambda2 from 8.5e-6 to 3.4e-6 for better validation, still proplems with plateau for solvated ions

examples/ReproducableCode/DoubleLayerCapacity/Molarity.py

@@ -29,7 +29,7 @@ T = 293.75 # [K]
 epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
-nR_mol = 55 # [mol/m^3]
+nR_mol = 55 # [mol/m^3] # ? m^3 or dm^3? Should be dm^3 and then we get 55 mol/dm^3 * 1e+3 = 55 mol/m^3
 nR_m = nR_mol * NA * 1/(1e-3)# [1/m^3]
 pR = 1.01325 * 1e+5 # [Pa]
 LR = 20e-9

examples/ReproducableCode/ValidationDreyer.py

@@ -32,7 +32,7 @@ epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
 p_ref = 1.01325 * 1e+5 # [Pa]
-L_ref = 10e-9#20e-9#20e-9 # [m]
+L_ref = 20e-9 # 20e-9[m]
 chi = 80 # [-]
 
 # BCS according to Dreyer
@@ -53,9 +53,14 @@ K = 'incompressible'
 Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (L_ref)**2)
 a2 = (p_ref)/(nR_m * k * T)
 Lambda2, a2 # values not sure, as T, chi and L_ref not given in Dreyer paper
+# ! Lambda2 would be: 8.491106981682818e-06 with these parameter
+# Let's try out smaller values to reproduce the results from the Dreyer paper
+Lambda2 = 3.4e-6
+# Lambda2 = Lambda2 * 1e-1
 
 # Solver parameters
-number_cells = 1028*32
+# ! Refine for PB
+number_cells = 1028# ! *32
 relax_param = 0.005
 refinement_style = 'log'
 max_iter = 100_000
@@ -64,13 +69,14 @@ rtol = 1e-7
 
 
 y_A_dimless, y_C_dimless, y_S_dimless, phi_dimless, p_dimless, x_dimless = [], [], [], [], [], []
-y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_L_dimless, phi_R_dimless, pR, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=0, PoissonBoltzmann=True, relax_param=relax_param, refinement_style=refinement_style, return_type=return_type, max_iter=max_iter, rtol=rtol)
-y_A_dimless.append(y_A_)
-y_C_dimless.append(y_C_)
-y_S_dimless.append(1 - y_A_ - y_C_)
-phi_dimless.append(phi_)
-p_dimless.append(p_)
-x_dimless.append(x_)
+# ! Uncomment
+# y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_L_dimless, phi_R_dimless, pR, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=0, PoissonBoltzmann=True, relax_param=relax_param, refinement_style=refinement_style, return_type=return_type, max_iter=max_iter, rtol=rtol)
+# y_A_dimless.append(y_A_)
+# y_C_dimless.append(y_C_)
+# y_S_dimless.append(1 - y_A_ - y_C_)
+# phi_dimless.append(phi_)
+# p_dimless.append(p_)
+# x_dimless.append(x_)
 
 y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_L_dimless, phi_R_dimless, pR, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells=number_cells, solvation=0, PoissonBoltzmann=False, relax_param=relax_param, refinement_style=refinement_style, return_type=return_type, max_iter=max_iter, rtol=rtol)
 y_A_dimless.append(y_A_)
@@ -101,4 +107,67 @@ n_C_dim = y_C_dim * nRef
 n_S_dim = y_S_dim * nRef
 
 # Store the results
-np.savez('../Data/Validation_Dreyer_BulkSurface.npz', y_A_dimless=y_A_dimless, y_C_dimless=y_C_dimless, y_S_dimless=y_S_dimless, phi_dimless=phi_dimless, p_dimless=p_dimless, x_dimless=x_dimless, y_A_dim=y_A_dim, y_C_dim=y_C_dim, y_S_dim=y_S_dim, n_A_dim=n_A_dim, n_C_dim=n_C_dim, n_S_dim=n_S_dim, phi_dim=phi_dim, p_dim=p_dim, x_dim=x_dim)
\ No newline at end of file
+np.savez('../Data/Validation_Dreyer_BulkSurface.npz', y_A_dimless=y_A_dimless, y_C_dimless=y_C_dimless, y_S_dimless=y_S_dimless, phi_dimless=phi_dimless, p_dimless=p_dimless, x_dimless=x_dimless, y_A_dim=y_A_dim, y_C_dim=y_C_dim, y_S_dim=y_S_dim, n_A_dim=n_A_dim, n_C_dim=n_C_dim, n_S_dim=n_S_dim, phi_dim=phi_dim, p_dim=p_dim, x_dim=x_dim, Lambda2=Lambda2, a2=a2)
+
+# Visualize the results if main
+if __name__ == '__main__':
+    # Options
+    markers = [':', '-.', '-'] #['-', '-.', ':', '--']
+    colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
+    labelsize = 15
+    lw = 2
+    legend_width = 3
+
+    scale_x = 1e-9
+    ticks_x = ticker.FuncFormatter(lambda x, pos: '{0:g}'.format(x/scale_x))
+
+    # Figure
+    fig, axs = plt.subplots(ncols=2, figsize=(13,5))
+    axs[0].plot(x_dim[0], n_A_dim[0], markers[0], color=colors[0], lw=lw, label='Poisson-Boltzmann')
+    axs[0].plot(x_dim[1], n_A_dim[1], markers[1], color=colors[1], lw=lw, label=f'ideal mixture, $\kappa = 0$')
+    # axs[0].plot(x_dim[2], n_A_dim[2], markers[2], color=colors[2], lw=lw, label=f'solvated ions, $\kappa = 8$')
+    axs[0].grid()
+    axs[0].xaxis.set_major_formatter(ticks_x)
+    axs[0].set_xlim(0, 2.5*scale_x)
+    axs[0].set_ylim(0, 60)
+    axs[0].set_xlabel('x [nm]', fontsize=labelsize)
+    axs[0].set_ylabel('$n_A$ [mol/m$^3$]', fontsize=labelsize)
+
+    a = plt.axes([.17, .4, .29, .5])
+    a.plot(x_dim[0], n_C_dim[0], markers[0], color=colors[0], lw=lw)
+    a.plot(x_dim[1], n_C_dim[1], markers[1], color=colors[1], lw=lw)
+    # a.plot(x_dim[2], n_C_dim[2], markers[2], color=colors[2], lw=lw)
+    a.grid()
+    a.xaxis.set_major_formatter(ticks_x)
+    a.set_xlim(0, 4.5*scale_x)
+    a.set_xlabel('x [nm]', fontsize=labelsize)
+    a.set_ylabel('$n_C$ [mol/m$^3$]', fontsize=labelsize)
+
+
+
+    axs[1].plot(x_dim[0], phi_dim[0], markers[0], color=colors[0], lw=lw)
+    axs[1].plot(x_dim[1], phi_dim[1], markers[1], color=colors[1], lw=lw)
+    # axs[1].plot(x_dim[2], phi_dim[2], markers[2], color=colors[2], lw=lw)
+    axs[1].grid()
+    axs[1].xaxis.set_major_formatter(ticks_x)
+    axs[1].set_xlim(0, 2.5*scale_x)
+    axs[1].set_xlabel('x [nm]', fontsize=labelsize)
+    axs[1].set_ylabel('$\\varphi [V]$', fontsize=labelsize)
+
+    a = plt.axes([.18+.5, .4, .29, .5])
+    a.plot(x_dim[0], p_dim[0]*1e-6, markers[0], color=colors[0], lw=lw)
+    a.plot(x_dim[1], p_dim[1]*1e-6, markers[1], color=colors[1], lw=lw)
+    # a.plot(x_dim[2], p_dim[2]*1e-6, markers[2], color=colors[2], lw=lw)
+    a.grid()
+    a.xaxis.set_major_formatter(ticks_x)
+    a.set_xlim(0, 0.9*scale_x)
+    a.set_ylim(0, 2_700)
+    a.set_xlabel('x [nm]', fontsize=labelsize)
+    a.set_ylabel('$p [MPa]$', fontsize=labelsize)
+
+    lgnd = fig.legend(bbox_to_anchor=(0.85, 1.1), fontsize=labelsize, ncol=6)
+    for line in lgnd.get_lines():
+        line.set_linewidth(legend_width)
+
+    fig.tight_layout()
+    fig.show()
\ No newline at end of file

examples/Visualizations/VisValidationDreyer.py

@@ -27,15 +27,17 @@ legend_width = 3
 
 scale_x = 1e-9
 ticks_x = ticker.FuncFormatter(lambda x, pos: '{0:g}'.format(x/scale_x))
+# labels = ['Poisson-Boltzmann', 'ideal mixture, $\kappa = 0$', 'solvated ions, $\kappa = 8$']
+labels = ['ideal mixture, $\kappa = 0$', 'solvated ions, $\kappa = 8$']
 
 # Figure
 fig, axs = plt.subplots(ncols=2, figsize=(13,5))
-axs[0].plot(x[0], n_A[0], markers[0], color=colors[0], lw=lw, label='Poisson-Boltzmann')
-axs[0].plot(x[1], n_A[1], markers[1], color=colors[1], lw=lw, label=f'ideal mixture, $\kappa = 0$')
-axs[0].plot(x[2], n_A[2], markers[2], color=colors[2], lw=lw, label=f'solvated ions, $\kappa = 8$')
+axs[0].plot(x[0], n_A[0], markers[0], color=colors[0], lw=lw, label=labels[0])
+axs[0].plot(x[1], n_A[1], markers[1], color=colors[1], lw=lw, label=labels[1])
+# axs[0].plot(x[2], n_A[2], markers[2], color=colors[2], lw=lw, label=f'solvated ions, $\kappa = 8$')
 axs[0].grid()
 axs[0].xaxis.set_major_formatter(ticks_x)
-axs[0].set_xlim(0, 2.5*scale_x)
+axs[0].set_xlim(0, 1.5*scale_x)
 axs[0].set_ylim(0, 60)
 axs[0].set_xlabel('x [nm]', fontsize=labelsize)
 axs[0].set_ylabel('$n_A$ [mol/m$^3$]', fontsize=labelsize)
@@ -43,10 +45,10 @@ axs[0].set_ylabel('$n_A$ [mol/m$^3$]', fontsize=labelsize)
 a = plt.axes([.17, .4, .29, .5])
 a.plot(x[0], n_C[0], markers[0], color=colors[0], lw=lw)
 a.plot(x[1], n_C[1], markers[1], color=colors[1], lw=lw)
-a.plot(x[2], n_C[2], markers[2], color=colors[2], lw=lw)
+# a.plot(x[2], n_C[2], markers[2], color=colors[2], lw=lw)
 a.grid()
 a.xaxis.set_major_formatter(ticks_x)
-a.set_xlim(0, 4.5*scale_x)
+a.set_xlim(0, 2.5*scale_x)
 a.set_xlabel('x [nm]', fontsize=labelsize)
 a.set_ylabel('$n_C$ [mol/m$^3$]', fontsize=labelsize)
 
@@ -54,7 +56,7 @@ a.set_ylabel('$n_C$ [mol/m$^3$]', fontsize=labelsize)
 
 axs[1].plot(x[0], phi[0], markers[0], color=colors[0], lw=lw)
 axs[1].plot(x[1], phi[1], markers[1], color=colors[1], lw=lw)
-axs[1].plot(x[2], phi[2], markers[2], color=colors[2], lw=lw)
+# axs[1].plot(x[2], phi[2], markers[2], color=colors[2], lw=lw)
 axs[1].grid()
 axs[1].xaxis.set_major_formatter(ticks_x)
 axs[1].set_xlim(0, 2.5*scale_x)
@@ -64,7 +66,7 @@ axs[1].set_ylabel('$\\varphi [V]$', fontsize=labelsize)
 a = plt.axes([.18+.5, .4, .29, .5])
 a.plot(x[0], p[0]*1e-6, markers[0], color=colors[0], lw=lw)
 a.plot(x[1], p[1]*1e-6, markers[1], color=colors[1], lw=lw)
-a.plot(x[2], p[2]*1e-6, markers[2], color=colors[2], lw=lw)
+# a.plot(x[2], p[2]*1e-6, markers[2], color=colors[2], lw=lw)
 a.grid()
 a.xaxis.set_major_formatter(ticks_x)
 a.set_xlim(0, 0.9*scale_x)