Coverage for fxdgm/Eq02.py: 96%
101 statements
« prev ^ index » next coverage.py v7.6.12, created at 2025-02-13 12:39 +0000
1'''
2Jan Habscheid
3Jan.Habscheid@rwth-aachen.de
5This file implements the fenics solver for the reduced systme of equations to two equations for the electric potential and the pressure.
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, Mesh, exp
14from basix.ufl import element, mixed_element
15import matplotlib.pyplot as plt
16from fxdgm.RefinedMesh1D import create_refined_mesh
19def solve_System_2eq(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 simplified 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 with the space charge
29 n^F = z_A y_A(φ, p) + z_C y_C(φ ,p)
31 if the mixture is compressible:
32 y_alpha = C_alpha * (K+p−1)^(−κ+1)a²K exp(−z_α φ)
34 if the mixture is incompressible:
35 y_alpha = D_alpha * exp(−(κ+1)a²p−z_α φ)
37 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.
39 ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
41 Parameters
42 ----------
43 phi_left : float
44 Value of φ at the left boundary
45 phi_right : float
46 Value of φ at the right boundary
47 p_right : float
48 Value of p at the right boundary
49 z_A : float
50 Charge number of species A
51 z_C : float
52 Charge number of species C
53 y_A_R : float
54 Atomic fractions of species A at right boundary
55 y_C_R : float
56 Atomic fractions of species C at right boundary
57 K : float | str
58 Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
59 Lambda2 : float
60 Dimensionless parameter
61 a2 : float
62 Dimensionless parameter
63 number_cells : int
64 Number of cells in the mesh
65 solvation : float, optional
66 solvation number, by default 0
67 PoissonBoltzmann : bool, optional
68 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
69 relax_param : float, optional
70 Relaxation parameter for the Newton solver
71 xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
72 , by default None -> Determined automatically
73 x0 : float, optional
74 Left boundary of the domain, by default 0
75 x1 : float, optional
76 Right boundary of the domain, by default 1
77 refinement_style : str, optional
78 Specify for refinement towards zero
79 Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
80 return_type : str, optional
81 'Vector' or 'Scalar', 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Scalar'
82 rtol : float, optional
83 Relative tolerance for Newton solver, by default 1e-8
84 max_iter : float, optional
85 Maximum number of Newton iterations, by default 500
87 Returns
88 -------
89 y_A, y_C, phi, p, msh
90 Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
91 If return_type is 'Vector', the solution is returned as numpy arrays
92 '''
93 if return_type == 'Scalar':
94 raise NotImplementedError('Scalar return type is not implemented yet')
95 # Define boundaries of the domain
96 x0 = 0
97 x1 = 1
99 # Define boundaries for the boundary conditions
100 def Left(x):
101 return np.isclose(x[0], x0)
103 def Right(x):
104 return np.isclose(x[0], x1)
106 # Create mesh
107 if refinement_style == 'uniform':
108 msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
109 else:
110 msh = create_refined_mesh(refinement_style, number_cells)
112 # Define Finite Elements
113 CG1_elem = element('Lagrange', msh.basix_cell(), 1)
115 # Define Mixed Function Space
116 W_elem = mixed_element([CG1_elem, CG1_elem])#, CG1_elem, CG1_elem])
117 W = fem.functionspace(msh, W_elem)
119 # Define Trial- and Testfunctions
120 u = fem.Function(W)
121 phi, p = split(u)
122 (v_1, v_2) = TestFunctions(W)
124 # Collapse function space for bcs
125 W0, _ = W.sub(0).collapse()
126 W1, _ = W.sub(1).collapse()
128 # Define boundary conditions values
129 def phi_left_(x):
130 return np.full_like(x[0], phi_left)
131 def phi_right_(x):
132 return np.full_like(x[0], phi_right)
133 def p_right_(x):
134 return np.full_like(x[0], p_right)
136 # Interpolate bcs functions
137 phi_left_bcs = fem.Function(W0)
138 phi_left_bcs.interpolate(phi_left_)
139 phi_right_bcs = fem.Function(W0)
140 phi_right_bcs.interpolate(phi_right_)
141 p_right_bcs = fem.Function(W1)
142 p_right_bcs.interpolate(p_right_)
144 # Identify dofs for boundary conditions
145 # Define boundary conditions
146 facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
147 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
148 bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
149 bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
151 facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
152 bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
155 # Combine boundary conditions into list
156 bcs = [bc_left_phi, bc_right_phi, bc_right_p]
158 def y_A(phi, p):
159 if PoissonBoltzmann == False:
160 D_A = y_A_R / exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
161 return D_A * exp(-(solvation + 1) * a2 * p - z_A * phi)
162 elif PoissonBoltzmann == True:
163 D_A = y_A_R / exp(- z_A * phi_right)
164 return D_A * exp(- z_A * phi)
166 def y_C(phi, p):
167 if PoissonBoltzmann == False:
168 D_C = y_C_R / exp(-(solvation + 1) * a2 * - z_C * phi_right)
169 return D_C * exp(-(solvation + 1) * a2 * p - z_C * phi)
170 elif PoissonBoltzmann == True:
171 D_C = y_C_R / exp(- z_C * phi_right)
172 return D_C * exp(- z_C * phi)
175 # Define variational problem
176 if K == 'incompressible':
177 # total free charge density
178 def nF(y_A, y_C, p):
179 return (z_C * y_C + z_A * y_A)
180 else:
181 # total number density
182 def n(p):
183 return (p-1)/K + 1
185 # total free charge density
186 def nF(y_A, y_C, p):
187 return (z_C * y_C + z_A * y_A) * n(p)
188 # Variational Form
189 A = (
190 inner(grad(phi), grad(v_1)) * dx
191 - 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p), p) * v_1 * dx
192 ) + (
193 inner(grad(p), grad(v_2)) * dx
194 + 1 / a2 * nF(y_A(phi, p), y_C(phi, p), p) * dot(grad(phi), grad(v_2)) * dx
195 )
196 F = A
198 # Define Nonlinear Problem
199 problem = NonlinearProblem(F, u, bcs=bcs)
201 # Define Newton Solver and solver settings
202 solver = NewtonSolver(MPI.COMM_WORLD, problem)
203 solver.convergence_criterion = "incremental"
204 solver.rtol = rtol
205 if relax_param != None:
206 solver.relaxation_parameter = relax_param
207 else:
208 if phi_right == phi_left:
209 solver.relaxation_parameter = 1.0
210 else:
211 solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
212 solver.max_it = max_iter
213 solver.report = True
215 # Solve the problem
216 log.set_log_level(log.LogLevel.INFO)
217 n, converged = solver.solve(u)
218 assert (converged)
219 print(f"Number of interations: {n:d}")
221 # Split the mixed function space into the individual components
222 phi, p = u.split()
224 # Return the solution
225 if return_type=='Vector':
226 x_vals = np.array(msh.geometry.x[:,0])
227 phi_vals = np.array(u.sub(0).collapse().x.array)
228 p_vals = np.array(u.sub(1).collapse().x.array)
230 # Calculate the atomic fractions
231 D_A = y_A_R / np.exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
232 y_A_vals = D_A * np.exp(-(solvation + 1) * a2 * p_vals - z_A * phi_vals)
234 D_C = y_C_R / np.exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
235 y_C_vals = D_C * np.exp(-(solvation + 1) * a2 * p_vals - z_C * phi_vals)
237 if PoissonBoltzmann:
238 D_A = y_A_R / np.exp(- z_A * phi_right)
239 y_A_vals = D_A * np.exp(- z_A * phi_vals)
240 D_C = y_C_R / np.exp(- z_C * phi_right)
241 y_C_vals = D_C * np.exp(- z_C * phi_vals)
243 return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
245if __name__ == '__main__': # pragma: no cover # dont cover main in coverage
246 # Define the parameters
247 phi_left = 5.0
248 phi_right = 0.0
249 p_right = 0.0
250 y_A_R = 1/3
251 y_C_R = 1/3
252 z_A = -1.0
253 z_C = 1.0
254 K = 'incompressible'
255 Lambda2 = 8.553e-6
256 a2 = 7.5412e-4
257 number_cells = 1024
258 relax_param = .1
259 rtol = 1e-4
260 max_iter = 500
262 # Solve the system
263 y_A, y_C, phi, p, x = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
265 # Plot the solution
266 plt.plot(x, phi)
267 plt.xlim(0,0.05)
268 plt.grid()
269 plt.xlabel('x [-]')
270 plt.ylabel('$\\varphi$ [-]')
271 plt.show()
273 plt.plot(x, y_A, '--', color='tab:blue', label='$y_A$')
274 plt.plot(x, y_C, '-', color='tab:blue', label='$y_C$')
275 plt.plot(x, 1 - y_A - y_C, ':', color='tab:blue', label='$y_S$')
276 plt.xlim(0,0.05)
277 plt.legend()
278 plt.grid()
279 plt.xlabel('x [-]')
280 plt.ylabel('$y_\\alpha$ [-]')
281 plt.show()
283 plt.plot(x, p)
284 plt.xlim(0,0.05)
285 plt.grid()
286 plt.xlabel('x [-]')
287 plt.ylabel('$p$ [-]')
288 plt.show()