Coverage for fxdgm/Eq04.py: 96%
113 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 script implements the fenics solver for a ternary electrolyte (N=3, A,C,S)
6'''
8import numpy as np
9from mpi4py import MPI
10from dolfinx import mesh, fem, log
11from dolfinx.fem.petsc import NonlinearProblem
12from dolfinx.nls.petsc import NewtonSolver
13from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
14from basix.ufl import element, mixed_element
15import matplotlib.pyplot as plt
16from fxdgm.RefinedMesh1D import create_refined_mesh
19def solve_System_4eq(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):
20 '''
21 Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
23 System of equations:
24 λ²Δ φ =−L²n^F
26 a²∇p=−n^F∇ φ
28 div(J_A)=0
30 div(J_C)=0
32 with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.
34 ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
36 Parameters
37 ----------
38 phi_left : float
39 Value of φ at the left boundary
40 phi_right : float
41 Value of φ at the right boundary
42 p_right : float
43 Value of p at the right boundary
44 z_A : float
45 Charge number of species A
46 z_C : float
47 Charge number of species C
48 y_A_R : float
49 Atomic fractions of species A at right boundary
50 y_C_R : float
51 Atomic fractions of species C at right boundary
52 K : float | str
53 Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
54 Lambda2 : float
55 Dimensionless parameter
56 a2 : float
57 Dimensionless parameter
58 number_cells : int
59 Number of cells in the mesh
60 solvation : float, optional
61 solvation number, by default 0
62 PoissonBoltzmann : bool, optional
63 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
64 relax_param : float, optional
65 Relaxation parameter for the Newton solver
66 xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
67 , by default None -> Determined automatically
68 x0 : float, optional
69 Left boundary of the domain, by default 0
70 x1 : float, optional
71 Right boundary of the domain, by default 1
72 refinement_style : str, optional
73 Specify for refinement towards zero
74 Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
75 return_type : str, optional
76 'Vector' or 'Scalar', 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Scalar'
77 rtol : float, optional
78 Relative tolerance for Newton solver, by default 1e-8
79 max_iter : float, optional
80 Maximum number of Newton iterations, by default 500
82 Returns
83 -------
84 y_A, y_C, phi, p, msh
85 Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
86 If return_type is 'Vector', the solution is returned as numpy arrays
87 '''
88 # Define boundaries of the domain
89 x0 = 0
90 x1 = 1
92 # Define boundaries for the boundary conditions
93 def Left(x):
94 return np.isclose(x[0], x0)
96 def Right(x):
97 return np.isclose(x[0], x1)
99 # Create mesh
100 if refinement_style == 'uniform':
101 msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
102 else:
103 msh = create_refined_mesh(refinement_style, number_cells)
105 # Define Finite Elements
106 CG1_elem = element('Lagrange', msh.basix_cell(), 1)
108 # Define Mixed Function Space
109 W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
110 W = fem.functionspace(msh, W_elem)
112 # Define Trial- and Testfunctions
113 u = fem.Function(W)
114 y_A, y_C, phi, p = split(u)
115 (v_A, v_C, v_1, v_2) = TestFunctions(W)
117 # Collapse function space for bcs
118 W0, _ = W.sub(0).collapse()
119 W1, _ = W.sub(1).collapse()
120 W2, _ = W.sub(2).collapse()
121 W3, _ = W.sub(3).collapse()
123 # Define boundary conditions values
124 def phi_left_(x):
125 return np.full_like(x[0], phi_left)
126 def phi_right_(x):
127 return np.full_like(x[0], phi_right)
128 def p_right_(x):
129 return np.full_like(x[0], p_right)
130 def y_A_right_(x):
131 return np.full_like(x[0], y_A_R)
132 def y_C_right_(x):
133 return np.full_like(x[0], y_C_R)
135 # Interpolate bcs functions
136 phi_left_bcs = fem.Function(W2)
137 phi_left_bcs.interpolate(phi_left_)
138 phi_right_bcs = fem.Function(W2)
139 phi_right_bcs.interpolate(phi_right_)
140 p_right_bcs = fem.Function(W3)
141 p_right_bcs.interpolate(p_right_)
142 y_A_right_bcs = fem.Function(W0)
143 y_A_right_bcs.interpolate(y_A_right_)
144 y_C_right_bcs = fem.Function(W1)
145 y_C_right_bcs.interpolate(y_C_right_)
147 # Identify dofs for boundary conditions
148 # Define boundary conditions
149 facet_left_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Left)
150 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Right)
151 bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(2))
152 bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(2))
154 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Right)
155 bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(3))
157 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
158 bc_right_y_A = fem.dirichletbc(y_A_right_bcs, facet_right_dofs, W.sub(0))
160 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
161 bc_right_y_C = fem.dirichletbc(y_C_right_bcs, facet_right_dofs, W.sub(1))
163 # Combine boundary conditions into list
164 bcs = [bc_left_phi, bc_right_phi, bc_right_p, bc_right_y_A, bc_right_y_C]
168 # Define variational problem
169 if K == 'incompressible':
170 # total free charge density
171 def nF(y_A, y_C):
172 return (z_C * y_C + z_A * y_A)
174 # Diffusion fluxes for species A and C
175 def J_A(y_A, y_C, phi, p):
176 return grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
177 # return grad(ln(y_A) + a2 * (p - 1) - solvation * ln(1-y_A-y_C) + z_A * phi)
179 def J_C(y_A, y_C, phi, p):
180 return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
181 # return grad(ln(y_C) + a2 * (p - 1) - solvation * ln(1-y_A-y_C) + z_C * phi)
183 # Variational Form
184 A = (
185 inner(grad(phi), grad(v_1)) * dx
186 - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
187 ) + (
188 inner(grad(p), grad(v_2)) * dx
189 + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
190 ) + (
191 inner(J_A(y_A, y_C, phi, p), grad(v_A)) * dx
192 + inner(J_C(y_A, y_C, phi, p), grad(v_C)) * dx
193 )
194 if PoissonBoltzmann:
195 A += (
196 inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_A)) * dx
197 + inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_C)) * dx
198 )
199 else:
200 # total number density
201 def n(p):
202 return (p-1)/K + 1
204 # total free charge density
205 def nF(y_A, y_C, p):
206 return (z_C * y_C + z_A * y_A) * n(p)
208 # Diffusion fluxes for species A and C
209 def J_A(y_A, y_C, phi, p):
210 return ln(y_A) + a2 * (solvation + 1) * K * ln(1 + 1/K * (p-1)) + z_A * phi
212 def J_C(y_A, y_C, phi, p):
213 return ln(y_C) + a2 * (solvation + 1)* K * ln(1 + 1/K * (p-1)) + z_C * phi
215 A = (
216 inner(grad(phi), grad(v_1)) * dx
217 - 1 / Lambda2 * nF(y_A, y_C, p) * v_1 * dx
218 ) + (
219 inner(grad(p), grad(v_2)) * dx
220 + 1 / a2 * nF(y_A, y_C, p) * dot(grad(phi), grad(v_2)) * dx
221 ) + (
222 inner(grad(J_A(y_A, y_C, phi, p)), grad(v_A)) * dx
223 + inner(grad(J_C(y_A, y_C, phi, p)), grad(v_C)) * dx
224 )
225 F = A
227 # Initialize initial guess for u
228 y_C_init = fem.Function(W1)
229 y_A_init = fem.Function(W0)
230 y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_R))
231 y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_R))
233 with u.vector.localForm() as u_loc:
234 u_loc.set(0)
235 u.sub(0).interpolate(y_A_init)
236 u.sub(1).interpolate(y_C_init)
238 # Define Nonlinear Problem
239 problem = NonlinearProblem(F, u, bcs=bcs)
241 # Define Newton Solver and solver settings
242 solver = NewtonSolver(MPI.COMM_WORLD, problem)
243 solver.convergence_criterion = "incremental"
244 solver.rtol = rtol
245 if relax_param != None:
246 solver.relaxation_parameter = relax_param
247 else:
248 if phi_right == phi_left:
249 solver.relaxation_parameter = 1.0
250 else:
251 solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
252 solver.max_it = max_iter
253 solver.report = True
255 # Solve the problem
256 log.set_log_level(log.LogLevel.INFO)
257 n, converged = solver.solve(u)
258 assert (converged)
259 print(f"Number of interations: {n:d}")
261 # Split the mixed function space into the individual components
262 y_A, y_C, phi, p = u.split()
264 # Return the solution
265 if return_type=='Vector':
266 x_vals = np.array(msh.geometry.x[:,0])
267 y_A_vals = np.array(u.sub(0).collapse().x.array)
268 y_C_vals = np.array(u.sub(1).collapse().x.array)
269 phi_vals = np.array(u.sub(2).collapse().x.array)
270 p_vals = np.array(u.sub(3).collapse().x.array)
272 return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
273 elif return_type=='Scalar':
274 return y_A, y_C, phi, p, msh
277if __name__ == '__main__': # pragma: no cover # dont cover main in coverage
278 # Define the parameters
279 phi_left = 10.0
280 phi_right = 0.0
281 p_right = 0.0
282 y_A_R = 0.01#1/3
283 y_C_R = 0.01#1/3
284 z_A = -1.0
285 z_C = 1.0
286 K = 'incompressible'
287 Lambda2 = 8.553e-6
288 a2 = 7.5412e-4
289 solvation = 0
290 number_cells = 1024
291 relax_param = .1
292 rtol = 1e-4
293 max_iter = 500
295 # Solve the system
296 y_A, y_C, phi, p, x = 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=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
298 # Plot the solution
299 plt.plot(x, phi)
300 plt.xlim(0,0.05)
301 plt.grid()
302 plt.xlabel('x [-]')
303 plt.ylabel('$\\varphi$ [-]')
304 plt.show()
306 plt.plot(x, y_A, '--', color='tab:blue', label='$y_A$')
307 plt.plot(x, y_C, '-', color='tab:blue', label='$y_C$')
308 plt.plot(x, 1 - y_A - y_C, ':', color='tab:blue', label='$y_S$')
309 plt.xlim(0,0.05)
310 plt.legend()
311 plt.grid()
312 plt.xlabel('x [-]')
313 plt.ylabel('$y_\\alpha$ [-]')
314 plt.show()
316 plt.plot(x, p)
317 plt.xlim(0,0.05)
318 plt.grid()
319 plt.xlabel('x [-]')
320 plt.ylabel('$p$ [-]')
321 plt.show()