Coverage for fxdgm/EqN.py: 94%
110 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 the generic system of equations for N species
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
18def solve_System_Neq(phi_left:float, phi_right:float, p_right:float, z_alpha:list, y_R:list, 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='Vector', rtol:float=1e-8, max_iter:float=500):
19 '''
20 Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
22 System of equations:
23 λ²Δ φ =−L²n^F
25 a²∇p=−n^F∇ φ
27 div(J_α)=0 α∈ {1,...,N−1}
29 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.
31 ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
33 Parameters
34 ----------
35 phi_left : float
36 Value of φ at the left boundary
37 phi_right : float
38 Value of φ at the right boundary
39 p_right : float
40 Value of p at the right boundary
41 z_alpha : list
42 Charge numbers for species α = 1,...,N-1
43 y_R : list
44 Atomic fractions at right boundary for species α = 1,...,N-1
45 K : float | str
46 Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
47 Lambda2 : float
48 Dimensionless parameter
49 a2 : float
50 Dimensionless parameter
51 number_cells : int
52 Number of cells in the mesh
53 solvation : float, optional
54 solvation number, not implemented yet, by default 0
55 PoissonBoltzmann : bool, optional
56 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, Not implemented yet, by default False
57 relax_param : float, optional
58 Relaxation parameter for the Newton solver
59 xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
60 , by default None -> Determined automatically
61 x0 : float, optional
62 Left boundary of the domain, by default 0
63 x1 : float, optional
64 Right boundary of the domain, by default 1
65 refinement_style : str, optional
66 Specify for refinement towards zero
67 Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
68 return_type : str, optional
69 'Vector' or 'Scalar' (not implemented yet, should be implemented in a later version), 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Vector'
70 rtol : float, optional
71 Relative tolerance for Newton solver, by default 1e-8
72 max_iter : float, optional
73 Maximum number of Newton iterations, by default 500
75 Returns
76 -------
77 y_A, y_C, phi, p, msh
78 Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
79 If return_type is 'Vector', the solution is returned as numpy arrays
80 Only return_type 'Vector' is implemented yet
81 '''
82 if solvation != 0:
83 raise NotImplementedError('Solvation number not implemented yet')
84 if PoissonBoltzmann:
85 raise NotImplementedError('Poisson-Boltzmann not implemented yet')
86 # Define boundaries of the domain
87 x0 = 0
88 x1 = 1
90 # Define boundaries
91 def Left(x):
92 return np.isclose(x[0], x0)
94 def Right(x):
95 return np.isclose(x[0], x1)
97 # Create mesh
98 if refinement_style == 'uniform':
99 msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
100 else:
101 msh = create_refined_mesh(refinement_style, number_cells)
103 # Define Finite Elements
104 CG1_elem = element('Lagrange', msh.basix_cell(), 1)
106 # Define Mixed Function Space
107 Elem_list = [CG1_elem, CG1_elem]
108 [Elem_list.append(CG1_elem) for _ in range(len(z_alpha))]
109 W_elem = mixed_element(Elem_list)
110 W = fem.functionspace(msh, W_elem)
112 # Define Trial- and Testfunctions
113 u = fem.Function(W)
114 my_TrialFunctions = split(u)
115 my_TestFunctions = TestFunctions(W)
117 phi, p = my_TrialFunctions[0], my_TrialFunctions[1]
118 v_1, v_2 = my_TestFunctions[0], my_TestFunctions[1]
119 y_alpha = my_TrialFunctions[2:]
120 v_alpha = my_TestFunctions[2:]
122 # Collapse function space for bcs
123 W_ = []
124 [W_.append(W.sub(i).collapse()[0]) for i in range(len(z_alpha)+2)]
126 # Define boundary conditions values
127 def phi_left_(x):
128 return np.full_like(x[0], phi_left)
129 def phi_right_(x):
130 return np.full_like(x[0], phi_right)
131 def p_right_(x):
132 return np.full_like(x[0], p_right)
134 # Interpolate bcs functions
135 phi_left_bcs = fem.Function(W_[0])
136 phi_left_bcs.interpolate(phi_left_)
137 phi_right_bcs = fem.Function(W_[0])
138 phi_right_bcs.interpolate(phi_right_)
139 p_right_bcs = fem.Function(W_[1])
140 p_right_bcs.interpolate(p_right_)
142 # Identify dofs for boundary conditions
143 facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
144 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
145 bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
146 bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
148 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
149 bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
151 # Combine boundary conditions for electric potential and pressure into list
152 bcs = [bc_left_phi, bc_right_phi, bc_right_p]
154 # Repeat the same for the boundary conditoins for the atomic fractions
155 for i in range(len(z_alpha)):
156 y_right_bcs = fem.Function(W_[i+2])
157 def y_right_(x):
158 return np.full_like(x[0], y_R[i])
159 y_right_bcs.interpolate(y_right_)
160 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(i+2), W.sub(i+2).collapse()[0]), Right)
161 bc_right_y = fem.dirichletbc(y_right_bcs, facet_right_dofs, W.sub(i+2))
162 bcs.append(bc_right_y)
164 # Define variational problem
165 if K == 'incompressible':
166 # total free charge density
167 def nF(y_alpha):
168 nF = 0
169 for i in range(len(z_alpha)):
170 nF += z_alpha[i] * y_alpha[i]
171 return nF
173 # Diffusion fluxes for species A and C
174 def J_alpha(y_alpha, alpha, phi, p):
175 mu_alpha = ln(y_alpha[alpha])
176 mu_S = ln(1 - sum(y_alpha))
177 return grad(mu_alpha - mu_S + z_alpha[alpha] * phi)
179 # Variational Form
180 A = (
181 inner(grad(phi), grad(v_1)) * dx
182 - 1 / Lambda2 * nF(y_alpha) * v_1 * dx
183 ) + (
184 inner(grad(p), grad(v_2)) * dx
185 + 1 / a2 * nF(y_alpha) * dot(grad(phi), grad(v_2)) * dx
186 )
187 for alpha in range(len(z_alpha)):
188 A += (
189 inner(J_alpha(y_alpha, alpha, phi, p), grad(v_alpha[alpha])) * dx
190 )
191 if PoissonBoltzmann:
192 raise ValueError('Poisson-Boltzmann not implemented for incompressible systems')
193 else:
194 raise ValueError('Only incompressible systems are implemented')
195 F = A
197 # Initialize initial guess for u
198 with u.vector.localForm() as u_loc:
199 u_loc.set(0)
201 # Initialize initial guess for u
202 for alpha in range(len(z_alpha)):
203 y_alpha_init = fem.Function(W_[alpha+2])
204 y_alpha_init.interpolate(lambda x: np.full_like(x[0], y_R[alpha]))
205 u.sub(alpha+2).interpolate(y_alpha_init)
207 # Define Nonlinear Problem
208 problem = NonlinearProblem(F, u, bcs=bcs)
210 # Define Newton Solver
211 solver = NewtonSolver(MPI.COMM_WORLD, problem)
212 solver.convergence_criterion = "incremental"
213 solver.rtol = rtol
214 if relax_param != None:
215 solver.relaxation_parameter = relax_param
216 else:
217 if phi_right == phi_left:
218 solver.relaxation_parameter = 1.0
219 else:
220 solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
221 solver.max_it = max_iter
222 solver.report = True
225 log.set_log_level(log.LogLevel.INFO)
226 n, converged = solver.solve(u)
227 assert (converged)
228 print(f"Number of (interations: {n:d}")
230 # Return the solution
231 if return_type=='Vector':
232 x = np.array(msh.geometry.x[:,0])
233 phi = np.array(u.sub(0).collapse().x.array)
234 p = np.array(u.sub(1).collapse().x.array)
235 y = []
236 [y.append(u.sub(i+2).collapse().x.array) for i in range(len(z_alpha))]
237 y = np.array(y)
238 return y, phi, p, x
241if __name__ == '__main__': # pragma: no cover # dont cover main in coverage
242 # Define the parameters
243 phi_left = 8.0
244 phi_right = 0.0
245 p_right = 0.0
246 y_R = [3/6, 1/6, 1/6]
247 z_alpha = [-1.0, 1.0, 2.0]
248 K = 'incompressible'
249 Lambda2 = 8.553e-6
250 a2 = 7.5412e-4
251 number_cells = 1024
252 relax_param = .05
253 rtol = 1e-4
254 max_iter = 2_500
255 refinement_style = 'hard_log'
256 return_type = 'Vector'
258 # Solve the system
259 y, phi, p, x = solve_System_Neq(phi_left, phi_right, p_right, z_alpha, y_R, K, Lambda2, a2, number_cells, relax_param=relax_param, refinement_style=refinement_style, return_type=return_type, max_iter=max_iter, rtol=rtol)
261 # Plot the solution
262 plt.figure()
263 plt.plot(x, phi)
264 plt.xlim(0,0.05)
265 plt.grid()
266 plt.xlabel('x [-]')
267 plt.ylabel('$\\varphi$ [-]')
268 plt.show()
270 plt.figure()
271 plt.xlim(0,0.05)
272 plt.plot(x, p)
273 plt.grid()
274 plt.xlabel('x [-]')
275 plt.ylabel('$p$ [-]')
276 plt.show()
278 plt.figure()
279 for i in range(len(z_alpha)):
280 plt.plot(x, y[i], label=f'$y_{i}$')
281 plt.xlim(0,0.05)
282 plt.legend()
283 plt.grid()
284 plt.xlabel('x [-]')
285 plt.ylabel('$y_i$ [-]')
286 plt.show()