Coverage for fxdgm/ElectrolyticDiode.py: 94%
146 statements
« prev ^ index » next coverage.py v7.6.12, created at 2025-02-13 13:26 +0000
1'''
2Jan Habscheid
3Jan.Habscheid@rwth-aachen.de
5This module solves the dimensionless system of equations for the example of an electric diode
7# ! Disclaimer: The results conform with [A Numerical Strategy for Nernst–Planck Systems with Solvation Effect, J. Fuhrmann, DOI:10.1002/fuce.201500215]. As the size of the domain, the boundary conditions, and the parameters are chosen differently, it can not be expected to match the same results quantitatively but qualitatively. In [Fuhrmann], only the potential and the cations are visualized. The potential behaves the same, and the cations are pushed to the outside of the domain for the forward bias and to the center of the domain for the backward bias. Furthermore, a rectification effect in the concentrations of the ions cannot be seen, as the range of the concentrations is the same along the different biases. On the other hand, this rectification process can be observed in [Entropy and convergence analysis for two finite volume schemes for a Nernst-Planck-Poisson system with ion volume constraints, B. Gaudeaul, J. Fuhrmann, DOI:10.1007/s00211-022-01279-y]. The electric potential can be verified to match the behavior from [Gaudeaul & Fuhrmann]. Furthermore, the ion concentrations can also be validated qualitatively. However, in [Gaudeaul & Fuhrmann], a rectification effect can be observed, reducing the concentrations of anions and cations for the reverse bias and no bias compared to the forward bias.
8'''
10import numpy as np
11from mpi4py import MPI
12from dolfinx import mesh, fem, log, io
13from dolfinx.fem.petsc import NonlinearProblem
14from dolfinx.nls.petsc import NewtonSolver
15from ufl import TestFunctions, split, dot, grad, dx, inner, ln
16from basix.ufl import element, mixed_element
17import matplotlib.pyplot as plt
18from ufl import Measure
19from dolfinx.mesh import locate_entities, meshtags
21def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C:float, y_A_bath:float, y_C_bath:float, K:float|str, Lambda2:float, a2:float, number_cells:list, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, Lx:float=2, Ly:float=10, rtol:float=1e-8, max_iter:float=500, return_type:str='Vector'):
22 '''
23 Solves the system of equations for the example of an electric diode
25 Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
27 ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
29 Parameters
30 ----------
31 Bias_type : str
32 ForwardBias, NoBias, BackwardBias
33 phi_bias : float
34 Bias in φ
35 g_phi : float
36 Neumann boundary condition for φ
37 z_A : float
38 Charge number of species A
39 z_C : float
40 Charge number of species C
41 y_A_bath : float
42 Bath concentration of anions
43 y_C_bath : float
44 Bath concentration of cations
45 K : float | str
46 Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
47 ! Only implemented for incompressible mixtures, yet
48 Lambda2 : float
49 Dimensionless parameter
50 a2 : float
51 Dimensionless parameter
52 number_cells : int
53 Number of cells in the mesh
54 solvation : float, optional
55 solvation number, by default 0
56 PoissonBoltzmann : bool, optional
57 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
58 relax_param : float, optional
59 Relaxation parameter for the Newton solver
60 xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
61 , by default None -> Determined automatically
62 Lx : float, optional
63 Length of domain in x-direction, by default 2
64 Ly : float, optional
65 Length of domain in y-direction, by default 10
66 rtol : float, optional
67 Relative tolerance for Newton solver, by default 1e-8
68 max_iter : float, optional
69 Maximum number of Newton iterations, by default 500
70 return_type : str, optional
71 Type of return value, by default 'Vector'
72 'Vector' -> Returns the solution as a numpy array for triangle elements
73 'Extended' -> Returns the solution as both, a numpy array for triangle elements and the fenicsx function u
76 Returns
77 -------
78 y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
79 Returns atomic fractions for species A and C, electric potential, pressure, and the mesh as numpy arrays for triangle elements. x_vector is a list of both dimensions [x, y]
80 If return_type is 'Extended', also returns the fenicsx function u
81 '''
82 if K != 'incompressible':
83 raise NotImplementedError('Only incompressible electrolytes are implemented yet')
84 x0 = np.array([0, 0])
85 x1 = np.array([Lx, Ly])
86 match Bias_type:
87 case 'ForwardBias': phi_bias = phi_bias
88 case 'NoBias': phi_bias = 0
89 case 'BackwardBias': phi_bias = -phi_bias
90 case _: raise ValueError('Invalid Bias_type')
92 # Define boundaries
93 geom_tol = 1E-4 # ! Geometric tolerance, may need to be adjusted
94 def Left(x):
95 return np.isclose(x[0], x0[0], geom_tol, geom_tol)
97 def Right(x):
98 return np.isclose(x[0], x1[0], geom_tol, geom_tol)
100 def Top(x):
101 return np.isclose(x[1], x1[1], geom_tol, geom_tol)
103 def Bottom(x):
104 return np.isclose(x[1], x0[1], geom_tol, geom_tol)
106 def Right_Top(x):
107 NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
108 RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
109 return np.logical_and(RightPart, NeumannPart)
111 def Right_Bottom(x):
112 NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
113 RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
114 return np.logical_and(RightPart, NeumannPart)
116 def Left_Top(x):
117 NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
118 LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
119 return np.logical_and(LeftPart, NeumannPart)
121 def Left_Bottom(x):
122 NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
123 LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
124 return np.logical_and(LeftPart, NeumannPart)
126 # Create the mesh
127 msh = mesh.create_rectangle(MPI.COMM_WORLD, [x0, x1], number_cells)
129 # Define Finite Elements
130 CG1_elem = element('Lagrange', msh.basix_cell(), 1)
131 # from ufl import FiniteElement
132 # CG1_elem = FiniteElement("Lagrange", msh.ufl_cell(), 1)
134 # Define Mixed Function Space
135 W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
136 W = fem.functionspace(msh, W_elem)
138 # Define Trial- and Testfunctions
139 u = fem.Function(W)
140 y_A, y_C, phi, p = split(u)
141 (v_A, v_C, v_1, v_2) = TestFunctions(W)
143 # Define boundary conditions
144 W0, _ = W.sub(0).collapse()
145 W1, _ = W.sub(1).collapse()
146 W2, _ = W.sub(2).collapse()
147 W3, _ = W.sub(3).collapse()
149 def phi_bottom_(x):
150 return np.full_like(x[1], 0)
151 def phi_top_(x):
152 return np.full_like(x[1], phi_bias)
153 def y_A_bottom_(x):
154 return np.full_like(x[1], y_A_bath)
155 def y_A_top_(x):
156 return np.full_like(x[1], y_A_bath)
157 def y_C_bottom_(x):
158 return np.full_like(x[1], y_C_bath)
159 def y_C_top_(x):
160 return np.full_like(x[1], y_C_bath)
162 phi_bottom_bcs = fem.Function(W2)
163 phi_bottom_bcs.interpolate(phi_bottom_)
164 phi_top_bcs = fem.Function(W2)
165 phi_top_bcs.interpolate(phi_top_)
166 y_A_bottom_bcs = fem.Function(W0)
167 y_A_bottom_bcs.interpolate(y_A_bottom_)
168 y_A_top_bcs = fem.Function(W0)
169 y_A_top_bcs.interpolate(y_A_top_)
170 y_C_bottom_bcs = fem.Function(W1)
171 y_C_bottom_bcs.interpolate(y_C_bottom_)
172 y_C_top_bcs = fem.Function(W1)
173 y_C_top_bcs.interpolate(y_C_top_)
175 # Identify face tags and create dirichlet conditions
176 facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Bottom)
177 facet_top_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Top)
178 bc_bottom_phi = fem.dirichletbc(phi_bottom_bcs, facet_bottom_dofs, W.sub(2))
179 bc_top_phi = fem.dirichletbc(phi_top_bcs, facet_top_dofs, W.sub(2))
181 facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Bottom)
182 bc_bottom_y_A = fem.dirichletbc(y_A_bottom_bcs, facet_bottom_dofs, W.sub(0))
183 facet_top_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Top)
184 bc_top_y_A = fem.dirichletbc(y_A_top_bcs, facet_top_dofs, W.sub(0))
186 facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Bottom)
187 bc_bottom_y_C = fem.dirichletbc(y_C_bottom_bcs, facet_bottom_dofs, W.sub(1))
188 facet_top_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Top)
189 bc_top_y_C = fem.dirichletbc(y_C_top_bcs, facet_top_dofs, W.sub(1))
191 # Collect boundary conditions
192 bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_top_y_A, bc_bottom_y_C, bc_top_y_C]
193 # bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_bottom_y_C]
196 # def p_bottom_(x):
197 # return np.full_like(x[1], 0)
198 # def p_top_(x):
199 # return np.full_like(x[1], 0)
200 # p_bottom_bcs = fem.Function(W3)
201 # p_bottom_bcs.interpolate(p_bottom_)
202 # p_top_bcs = fem.Function(W3)
203 # p_top_bcs.interpolate(p_top_)
204 # facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Bottom)
205 # facet_top_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Top)
206 # bc_bottom_p = fem.dirichletbc(p_bottom_bcs, facet_bottom_dofs, W.sub(3))
207 # bc_top_p = fem.dirichletbc(p_top_bcs, facet_top_dofs, W.sub(3))
208 # bcs.append(bc_bottom_p)
209 # bcs.append(bc_top_p)
211 # Variational formulation
212 if K == 'incompressible':
213 # def nF(y_A, y_C):
214 # return (z_C * y_C + z_A * y_A) # n = 1
216 # A = (
217 # inner(grad(phi), grad(v_1)) * dx
218 # - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
219 # ) + (
220 # inner(grad(p), grad(v_2)) * dx
221 # + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
222 # ) + (
223 # inner(grad(ln(y_A) + a2 * solvation * p - ln(1-y_A-y_C) + z_A * phi), grad(v_A)) * dx
224 # + inner(grad(ln(y_C) + a2 * solvation * p- ln(1-y_A-y_C) + z_C * phi), grad(v_C)) * dx
225 # )
226 # def nF(y_A, y_C):
227 # # n = const = 1
228 # return (z_C * y_C + z_A * y_A)
230 # def J_A(y_A, y_C, phi, p):
231 # g_A = (solvation + 1) * a2 * (p - 1) # g_Aref, but constant and take gradient
232 # mu_A = g_A + ln(y_A)
233 # g_N = a2 * (p - 1) # solvation_N = 0, g_Sref, but constant and take gradient
234 # mu_N = g_N + ln(1 - y_A - y_C)
235 # return grad(mu_A - mu_N + z_A * phi)
236 # # return grad(ln(y_A) - ln(1 - y_A - y_C) + z_A * phi)
238 # def J_C(y_A, y_C, phi, p):
239 # g_C = (solvation + 1) * a2 * (p - 1) # g_Cref, but constant and take gradient
240 # mu_C = g_C + ln(y_C)
241 # g_N = a2 * (p - 1)
242 # mu_N = g_N + ln(1 - y_A - y_C)
243 # return grad(mu_C - mu_N + z_C * phi)
244 # # return grad(ln(y_C) - ln(1 - y_A - y_C) + z_C * phi)
246 def nF(y_A, y_C):
247 return (z_C * y_C + z_A * y_A)
249 # Diffusion fluxes for species A and C
250 def J_A(y_A, y_C, phi, p):
251 return grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
253 def J_C(y_A, y_C, phi, p):
254 return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
256 # if PoissonBoltzmann:
257 # def J_A(y_A, y_C, phi, p):
258 # return grad(ln(y_A) + z_A * phi)
259 # def J_C(y_A, y_C, phi, p):
260 # return grad(ln(y_C) + z_C * phi)
262 A = (
263 inner(grad(phi), grad(v_1)) * dx
264 - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
265 ) + (
266 inner(grad(p), grad(v_2)) * dx
267 + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
268 ) + (
269 inner(J_A(y_A, y_C, phi, p), grad(v_A)) * dx
270 + inner(J_C(y_A, y_C, phi, p), grad(v_C)) * dx
271 )
273 # Define Neumann boundaries
274 boundaries = [(0, Right_Top), (1, Right_Bottom), (2, Left_Top), (3, Left_Bottom)]
275 facet_indices, facet_markers = [], []
276 fdim = msh.topology.dim - 1
277 for (marker, locator) in boundaries:
278 facets = locate_entities(msh, fdim, locator)
279 facet_indices.append(facets)
280 facet_markers.append(np.full_like(facets, marker))
281 facet_indices = np.hstack(facet_indices).astype(np.int32)
282 facet_markers = np.hstack(facet_markers).astype(np.int32)
283 sorted_facets = np.argsort(facet_indices)
284 facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
285 ds = Measure("ds", domain=msh, subdomain_data=facet_tag)
286 L = 0
287 L += (
288 inner(g_phi, v_1) * ds(0)
289 + inner(-g_phi, v_1) * ds(1)
290 # Uncomment, if you want to apply the Neumann bcs on both sides/ left side
291 # + inner(g_phi, v_1) * ds(2)
292 # + inner(-g_phi, v_1) * ds(3)
293 )
295 F = A + L
297 # Initialize initial guess for u
298 y_C_init = fem.Function(W1)
299 y_A_init = fem.Function(W0)
300 y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_bath))
301 y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_bath))
303 with u.vector.localForm() as u_loc:
304 u_loc.set(0)
305 u.sub(0).interpolate(y_A_init)
306 u.sub(1).interpolate(y_C_init)
308 # Define Nonlinear Problem
309 problem = NonlinearProblem(F, u, bcs=bcs)
311 solver = NewtonSolver(MPI.COMM_WORLD, problem)
312 solver.convergence_criterion = 'residual' #"incremental"
313 solver.rtol = rtol
314 solver.relaxation_parameter = relax_param
315 solver.max_it = max_iter
316 solver.report = True
318 # Solve the problem
319 log.set_log_level(log.LogLevel.INFO)
320 n, converged = solver.solve(u)
321 assert (converged)
322 print(f"Number of interations: {n:d}")
324 # Extract solution
325 y_A_vector = u.sub(0).collapse().x.array
326 y_C_vector = u.sub(1).collapse().x.array
327 phi_vector = u.sub(2).collapse().x.array
328 p_vector = u.sub(3).collapse().x.array
330 X = W.sub(0).collapse()[0].tabulate_dof_coordinates()
331 x_vector = [X[:, 0], X[:, 1]]
333 if return_type == 'Vector':
334 return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
335 elif return_type == 'Extended':
336 return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector, u
337 else:
338 raise ValueError('Invalid return_type')
340if __name__ == '__main__': # pragma: no cover # dont cover main in coverage
341 phi_bias = 10#10
342 Bias_type = 'ForwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
343 g_phi = 350#5
344 y_fixed = 0.01#0.01
345 z_A = -1.0
346 z_C = 1.0
347 K = 'incompressible'
348 Lambda2 = 8.553e-6 # ! Change back to 1e-6
349 # g_phi *= np.sqrt(Lambda2) # ! Unsure about this scaling
350 # g_phi *= Lambda2
351 # g_phi = g_phi / np.sqrt(Lambda2)
352 a2 = 7.5412e-4
353 number_cells = [20, 128]#[20, 100]
354 Lx = 0.02
355 Ly = 0.1
356 x0 = np.array([0, 0])
357 x1 = np.array([Lx, Ly])
358 refinement_style = 'uniform'
359 solvation = 5
360 PoissonBoltzmann = False
361 rtol = 1e-3 # ToDo: Change back to 1e-8, currently just for testing
362 relax_param = 0.15 # 0.1
363 max_iter = 15_000
365 # Solve the system
366 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')
367 x, y = X[0], X[1]
368 y_S = 1 - y_A - y_C
370 # Plot results
371 levelsf = 10
372 levels = 10
373 fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(8,10))
375 c = axs[0,0].tricontourf(x, y, y_A, levelsf)
376 fig.colorbar(c)
377 axs[0,0].tricontour(x, y, y_A, levels, colors='black')
378 axs[0,0].set_title('$y_A$')
380 c = axs[0,1].tricontourf(x, y, y_C, levelsf)
381 fig.colorbar(c)
382 axs[0,1].tricontour(x, y, y_C, levels, colors='black')
383 axs[0,1].set_title('$y_C$')
385 c = axs[0,2].tricontourf(x, y, y_S, levelsf)
386 fig.colorbar(c)
387 axs[0,2].tricontour(x, y, y_S, levels, colors='black')
388 axs[0,2].set_title('$y_S$')
390 c = axs[1,0].tricontourf(x, y, phi, levelsf)
391 fig.colorbar(c)
392 axs[1,0].tricontour(x, y, phi, levels, colors='black')
393 axs[1,0].set_title('$\\varphi$')
395 c = axs[1,1].tricontourf(x, y, p, levelsf)
396 fig.colorbar(c)
397 axs[1,1].tricontour(x, y, p, levels, colors='black')
398 axs[1,1].set_title('$p$')
400 fig.tight_layout()
401 fig.show()