E-Diode not working, concentrations seem to be uncouupled, fixed 1/Lambda2 in...

E-Diode not working, concentrations seem to be uncouupled, fixed 1/Lambda2 in Var.Form and found out, that Lambda2 and a2 also determine maximum of the concentrations at Neumann Cond

examples/ReproducableCode/ElectrolyticDiode/BackwardBias.py

@@ -30,16 +30,16 @@ y_fixed = 0.01
 z_A = -1.0
 z_C = 1.0
 K = 'incompressible'
-Lambda2 = 8.553e-6
-a2 = 7.5412e-4
-number_cells = [40, 400]
+Lambda2 = 8.553e-2 # ! e-6
+a2 = 7.5412e-2 # ! e-4
+number_cells = [20,100] # ! [40, 400]
 Lx = 2
 Ly = 10
 x0 = np.array([0, 0])
 x1 = np.array([Lx, Ly])
 
 # Solver parameters
-rtol = 1e-8
+rtol = 1e-3 # ! e-8
 relax_param = 0.05
 max_iter = 15_000    
 

examples/ReproducableCode/ElectrolyticDiode/ForwardBias.py

@@ -30,16 +30,16 @@ y_fixed = 0.01
 z_A = -1.0
 z_C = 1.0
 K = 'incompressible'
-Lambda2 = 8.553e-6
-a2 = 7.5412e-4
-number_cells = [40, 400]
+Lambda2 = 8.553e-2 # ! e-6
+a2 = 7.5412e-2 # !e-4
+number_cells = [20, 100] # ! [40, 400]
 Lx = 2
 Ly = 10
 x0 = np.array([0, 0])
 x1 = np.array([Lx, Ly])
 
 # Solver parameters
-rtol = 1e-8
+rtol = 1e-3 # ! e-8
 relax_param = 0.05
 max_iter = 15_000    
 

examples/ReproducableCode/ElectrolyticDiode/NoBias.py

@@ -30,16 +30,16 @@ y_fixed = 0.01
 z_A = -1.0
 z_C = 1.0
 K = 'incompressible'
-Lambda2 = 8.553e-6
-a2 = 7.5412e-4
-number_cells = [40, 400]
+Lambda2 = 8.553e-2 # ! e-6
+a2 = 7.5412e-2 # ! e-4
+number_cells = [20, 100] # ! [40, 400]
 Lx = 2
 Ly = 10
 x0 = np.array([0, 0])
 x1 = np.array([Lx, Ly])
 
 # Solver parameters
-rtol = 1e-8
+rtol = 1e-3 # ! e-8
 relax_param = 0.05
 max_iter = 15_000    
 

examples/Visualizations/VisElectrolyticDiode.py

@@ -54,16 +54,19 @@ legend_width = 8
 levels_colorbar = 20
 levels_contour = 20
 
+vmap_concentrations = np.linspace(0, 1, levels_colorbar)
+vmap_potential = np.linspace(-20, 20, levels_colorbar)
+vmap_pressure = np.linspace(np.min(p_packed), np.max(p_packed), levels_colorbar)
+
 for bias, (y_A, y_C, y_S, phi, p, x, y) in enumerate(zip(y_A_packed, y_C_packed, y_S_packed, phi_packed, p_packed, x_packed, y_packed)): 
-    c = axs[bias,0].tricontourf(x, y, phi, cmap=color_theme_potential, levels=np.linspace(np.min(phi),np.max(phi),levels_colorbar))#, extend="both")
+    c = axs[bias,0].tricontourf(x, y, phi, cmap=color_theme_potential, levels=vmap_potential)#, extend="both")
     # c.cmap.set_under('k')
-    c.set_clim(np.min(phi), np.max(phi))
+    c.set_clim(-20, 20)
     cbar = fig.colorbar(c, ax=axs[bias,0])
     cbar.ax.tick_params(labelsize=labelsize)
     axs[bias,0].tricontour(x, y, phi, colors='black', levels=levels_contour)
     axs[bias,0].tick_params(axis='both', labelsize=labelsize)
 
-    vmap_concentrations = np.linspace(0, 1, levels_colorbar)
     c = axs[bias,1].tricontourf(x, y, y_A, cmap=color_theme_concentration, vmin=np.min(y_A), vmax=np.max(y_A), levels=vmap_concentrations)#, extend='both')
     # cbar = fig.colorbar(c, ax=axs[bias,1], boundaries=np.linspace(0,1,5))
     # cbar.solids.set_edgecolor("face")
@@ -78,12 +81,14 @@ for bias, (y_A, y_C, y_S, phi, p, x, y) in enumerate(zip(y_A_packed, y_C_packed,
     axs[bias,2].tick_params(axis='both', labelsize=labelsize)
 
     c = axs[bias,3].tricontourf(x, y, y_S, cmap=color_theme_concentration, levels=vmap_concentrations)
+    c.set_clim(0, 1)
     cbar = fig.colorbar(c, ax=axs[bias,3])
     cbar.ax.tick_params(labelsize=labelsize)
     axs[bias,3].tricontour(x, y, y_S, colors='black', levels=levels_contour)
     axs[bias,3].tick_params(axis='both', labelsize=labelsize)
 
-    c = axs[bias,4].tricontourf(x, y, p, cmap=color_theme_pressure, levels=np.linspace(np.min(p),np.max(p),levels_colorbar))
+    c = axs[bias,4].tricontourf(x, y, p, cmap=color_theme_pressure, levels=vmap_pressure)
+    c.set_clim(np.min(p_packed), np.max(p_packed))
     cbar = fig.colorbar(c, ax=axs[bias,4])
     cbar.ax.tick_params(labelsize=labelsize)
     axs[bias,4].tricontour(x, y, p, colors='black', levels=levels_contour)

src/ElectrolyticDiode.py

@@ -193,16 +193,32 @@ def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C
     if K == 'incompressible':
         def nF(y_A, y_C):
             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)
+        
+        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)
         
         A = (
             inner(grad(phi), grad(v_1)) * dx
-            - Lambda2 * nF(y_A, y_C) * 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
+            inner(J_A(y_A, y_C, phi, p), grad(v_A)) * dx
+            + inner(J_C(y_A, y_C, phi, p), grad(v_C)) * 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
         )
 
     # Define Neumann boundaries
@@ -274,19 +290,16 @@ def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C
 
 if __name__ == '__main__':
     phi_bias = 10
-    Bias_type = 'ForwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
+    Bias_type = 'BackwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
     g_phi = 5
     y_fixed = 0.01
     z_A = -1.0
     z_C = 1.0
     K = 'incompressible'
-    # Lambda2 = 1e-7#7#1e-3
-    # a2 = 1e-4
-    # Lambda2 = 1e-8#8.33e-8 
-    # a2 = 7.3656e-4
-    Lambda2 = 8.553e-6
-    a2 = 7.5412e-4
-    number_cells = [20, 200]
+    Lambda2 = 8.553e-2 # ! Change back to 1e-6
+    # g_phi *= np.sqrt(Lambda2) # ! Unsure about this scaling
+    a2 = 7.5412e-2
+    number_cells = [20, 100]
     Lx = 2
     Ly = 10
     x0 = np.array([0, 0])
@@ -296,7 +309,7 @@ if __name__ == '__main__':
     PoissonBoltzmann = False
     rtol = 1e-3 # ToDo: Change back to 1e-8, currently just for testing
     relax_param = 0.05
-    max_iter = 15_000    
+    max_iter = 15_000
 
     # Solve the system
     y_A, y_C, phi, p, X = ElectrolyticDiode(Bias_type, phi_bias, g_phi, z_A, z_C, y_fixed, y_fixed, K, Lambda2, a2, number_cells, solvation, PoissonBoltzmann, relax_param, Lx, Ly, rtol, max_iter, return_type='Vector')

src/Eq04.py

@@ -206,10 +206,10 @@ def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
         
         # 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
+            return grad(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
+            return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
         
         # Variational Form
         A = (
@@ -219,8 +219,8 @@ def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
             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
+            inner(J_A(y_A, y_C, phi, p), grad(v_A)) * dx
+            + inner(J_C(y_A, y_C, phi, p), grad(v_C)) * dx
         )
         if PoissonBoltzmann:
             A += (

src/EqN.py

@@ -210,7 +210,7 @@ def solve_System_Neq(phi_left:float, phi_right:float, p_right:float, z_alpha:lis
         def J_alpha(y_alpha, alpha, phi, p):
             mu_alpha = ln(y_alpha[alpha])
             mu_S = ln(1 - sum(y_alpha))
-            return mu_alpha - mu_S + z_alpha[alpha] * phi
+            return grad(mu_alpha - mu_S + z_alpha[alpha] * phi)
         
         # Variational Form
         A = (
@@ -222,7 +222,7 @@ def solve_System_Neq(phi_left:float, phi_right:float, p_right:float, z_alpha:lis
         )
         for alpha in range(len(z_alpha)):
             A += (
-                inner(grad(J_alpha(y_alpha, alpha, phi, p)), grad(v_alpha[alpha])) * dx
+                inner(J_alpha(y_alpha, alpha, phi, p), grad(v_alpha[alpha])) * dx
             )
         if PoissonBoltzmann:
             raise ValueError('Poisson-Boltzmann not implemented for incompressible systems')