Skip to content
Snippets Groups Projects
Commit 8d9097ee authored by Jan Habscheid's avatar Jan Habscheid
Browse files

Cleaning up from failed merges

parents 511cbc14 ba16f88f
No related branches found
No related tags found
1 merge request!8Package
Pipeline #1521022 passed
Stage: prepare
Stage: build
Hide whitespace changes
Inline Side-by-side
Showing
with 6 additions and 2020 deletions
.gitignore
+ 6
− 1
View file @ 8d9097ee
...@@ -10,4 +10,9 @@ ...@@ -10,4 +10,9 @@
# Ignore pip folder # Ignore pip folder
*egg-info/* *egg-info/*
# Ignore venv files # Ignore venv files
venv venv
\ No newline at end of file # Ignore .eggs
*.eggs
FENICSxDGM.egg-info
# Ignore build files
build
\ No newline at end of file
build/lib/FENICSxDGM/ElectrolyticDiode.py deleted 100644 → 0
+ 0
− 401
View file @ 511cbc14
'''
Jan Habscheid
Jan.Habscheid@rwth-aachen.de
This module solves the dimensionless system of equations for the example of an electric diode
# ! 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.
'''
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, log, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import TestFunctions, split, dot, grad, dx, inner, ln
from basix.ufl import element, mixed_element
import matplotlib.pyplot as plt
from ufl import Measure
from dolfinx.mesh import locate_entities, meshtags
def 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'):
'''
Solves the system of equations for the example of an electric diode
Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
! If the Newton solver diverges, you may try to reduce the relaxation parameter.
Parameters
----------
Bias_type : str
ForwardBias, NoBias, BackwardBias
phi_bias : float
Bias in φ
g_phi : float
Neumann boundary condition for φ
z_A : float
Charge number of species A
z_C : float
Charge number of species C
y_A_bath : float
Bath concentration of anions
y_C_bath : float
Bath concentration of cations
K : float | str
Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
! Only implemented for incompressible mixtures, yet
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
number_cells : int
Number of cells in the mesh
solvation : float, optional
solvation number, by default 0
PoissonBoltzmann : bool, optional
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
relax_param : float, optional
Relaxation parameter for the Newton solver
xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
, by default None -> Determined automatically
Lx : float, optional
Length of domain in x-direction, by default 2
Ly : float, optional
Length of domain in y-direction, by default 10
rtol : float, optional
Relative tolerance for Newton solver, by default 1e-8
max_iter : float, optional
Maximum number of Newton iterations, by default 500
return_type : str, optional
Type of return value, by default 'Vector'
'Vector' -> Returns the solution as a numpy array for triangle elements
'Extended' -> Returns the solution as both, a numpy array for triangle elements and the fenicsx function u
Returns
-------
y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
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]
If return_type is 'Extended', also returns the fenicsx function u
'''
if K != 'incompressible':
raise NotImplementedError('Only incompressible electrolytes are implemented yet')
x0 = np.array([0, 0])
x1 = np.array([Lx, Ly])
match Bias_type:
case 'ForwardBias': phi_bias = phi_bias
case 'NoBias': phi_bias = 0
case 'BackwardBias': phi_bias = -phi_bias
case _: raise ValueError('Invalid Bias_type')
# Define boundaries
geom_tol = 1E-4 # ! Geometric tolerance, may need to be adjusted
def Left(x):
return np.isclose(x[0], x0[0], geom_tol, geom_tol)
def Right(x):
return np.isclose(x[0], x1[0], geom_tol, geom_tol)
def Top(x):
return np.isclose(x[1], x1[1], geom_tol, geom_tol)
def Bottom(x):
return np.isclose(x[1], x0[1], geom_tol, geom_tol)
def Right_Top(x):
NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
return np.logical_and(RightPart, NeumannPart)
def Right_Bottom(x):
NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
return np.logical_and(RightPart, NeumannPart)
def Left_Top(x):
NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
return np.logical_and(LeftPart, NeumannPart)
def Left_Bottom(x):
NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
return np.logical_and(LeftPart, NeumannPart)
# Create the mesh
msh = mesh.create_rectangle(MPI.COMM_WORLD, [x0, x1], number_cells)
# Define Finite Elements
CG1_elem = element('Lagrange', msh.basix_cell(), 1)
# from ufl import FiniteElement
# CG1_elem = FiniteElement("Lagrange", msh.ufl_cell(), 1)
# Define Mixed Function Space
W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
W = fem.functionspace(msh, W_elem)
# Define Trial- and Testfunctions
u = fem.Function(W)
y_A, y_C, phi, p = split(u)
(v_A, v_C, v_1, v_2) = TestFunctions(W)
# Define boundary conditions
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()
W2, _ = W.sub(2).collapse()
W3, _ = W.sub(3).collapse()
def phi_bottom_(x):
return np.full_like(x[1], 0)
def phi_top_(x):
return np.full_like(x[1], phi_bias)
def y_A_bottom_(x):
return np.full_like(x[1], y_A_bath)
def y_A_top_(x):
return np.full_like(x[1], y_A_bath)
def y_C_bottom_(x):
return np.full_like(x[1], y_C_bath)
def y_C_top_(x):
return np.full_like(x[1], y_C_bath)
phi_bottom_bcs = fem.Function(W2)
phi_bottom_bcs.interpolate(phi_bottom_)
phi_top_bcs = fem.Function(W2)
phi_top_bcs.interpolate(phi_top_)
y_A_bottom_bcs = fem.Function(W0)
y_A_bottom_bcs.interpolate(y_A_bottom_)
y_A_top_bcs = fem.Function(W0)
y_A_top_bcs.interpolate(y_A_top_)
y_C_bottom_bcs = fem.Function(W1)
y_C_bottom_bcs.interpolate(y_C_bottom_)
y_C_top_bcs = fem.Function(W1)
y_C_top_bcs.interpolate(y_C_top_)
# Identify face tags and create dirichlet conditions
facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Bottom)
facet_top_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Top)
bc_bottom_phi = fem.dirichletbc(phi_bottom_bcs, facet_bottom_dofs, W.sub(2))
bc_top_phi = fem.dirichletbc(phi_top_bcs, facet_top_dofs, W.sub(2))
facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Bottom)
bc_bottom_y_A = fem.dirichletbc(y_A_bottom_bcs, facet_bottom_dofs, W.sub(0))
facet_top_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Top)
bc_top_y_A = fem.dirichletbc(y_A_top_bcs, facet_top_dofs, W.sub(0))
facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Bottom)
bc_bottom_y_C = fem.dirichletbc(y_C_bottom_bcs, facet_bottom_dofs, W.sub(1))
facet_top_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Top)
bc_top_y_C = fem.dirichletbc(y_C_top_bcs, facet_top_dofs, W.sub(1))
# Collect boundary conditions
bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_top_y_A, bc_bottom_y_C, bc_top_y_C]
# bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_bottom_y_C]
# def p_bottom_(x):
# return np.full_like(x[1], 0)
# def p_top_(x):
# return np.full_like(x[1], 0)
# p_bottom_bcs = fem.Function(W3)
# p_bottom_bcs.interpolate(p_bottom_)
# p_top_bcs = fem.Function(W3)
# p_top_bcs.interpolate(p_top_)
# facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Bottom)
# facet_top_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Top)
# bc_bottom_p = fem.dirichletbc(p_bottom_bcs, facet_bottom_dofs, W.sub(3))
# bc_top_p = fem.dirichletbc(p_top_bcs, facet_top_dofs, W.sub(3))
# bcs.append(bc_bottom_p)
# bcs.append(bc_top_p)
# Variational formulation
if K == 'incompressible':
# def nF(y_A, y_C):
# return (z_C * y_C + z_A * y_A) # n = 1
# A = (
# inner(grad(phi), grad(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
# )
# def nF(y_A, y_C):
# # n = const = 1
# 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)
# # return grad(ln(y_A) - ln(1 - y_A - y_C) + 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)
# # return grad(ln(y_C) - ln(1 - y_A - y_C) + z_C * phi)
def nF(y_A, y_C):
return (z_C * y_C + z_A * y_A)
# Diffusion fluxes for species A and C
def J_A(y_A, y_C, phi, p):
return grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
def J_C(y_A, y_C, phi, p):
return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
# if PoissonBoltzmann:
# def J_A(y_A, y_C, phi, p):
# return grad(ln(y_A) + z_A * phi)
# def J_C(y_A, y_C, phi, p):
# return grad(ln(y_C) + z_C * phi)
A = (
inner(grad(phi), grad(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(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
)
# Define Neumann boundaries
boundaries = [(0, Right_Top), (1, Right_Bottom), (2, Left_Top), (3, Left_Bottom)]
facet_indices, facet_markers = [], []
fdim = msh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(msh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=msh, subdomain_data=facet_tag)
L = 0
L += (
inner(g_phi, v_1) * ds(0)
+ inner(-g_phi, v_1) * ds(1)
# Uncomment, if you want to apply the Neumann bcs on both sides/ left side
# + inner(g_phi, v_1) * ds(2)
# + inner(-g_phi, v_1) * ds(3)
)
F = A + L
# Initialize initial guess for u
y_C_init = fem.Function(W1)
y_A_init = fem.Function(W0)
y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_bath))
y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_bath))
with u.vector.localForm() as u_loc:
u_loc.set(0)
u.sub(0).interpolate(y_A_init)
u.sub(1).interpolate(y_C_init)
# Define Nonlinear Problem
problem = NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = 'residual' #"incremental"
solver.rtol = rtol
solver.relaxation_parameter = relax_param
solver.max_it = max_iter
solver.report = True
# Solve the problem
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(u)
assert (converged)
print(f"Number of interations: {n:d}")
# Extract solution
y_A_vector = u.sub(0).collapse().x.array
y_C_vector = u.sub(1).collapse().x.array
phi_vector = u.sub(2).collapse().x.array
p_vector = u.sub(3).collapse().x.array
X = W.sub(0).collapse()[0].tabulate_dof_coordinates()
x_vector = [X[:, 0], X[:, 1]]
if return_type == 'Vector':
return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
elif return_type == 'Extended':
return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector, u
else:
raise ValueError('Invalid return_type')
if __name__ == '__main__':
phi_bias = 10#10
Bias_type = 'ForwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
g_phi = 350#5
y_fixed = 0.01#0.01
z_A = -1.0
z_C = 1.0
K = 'incompressible'
Lambda2 = 8.553e-6 # ! Change back to 1e-6
# g_phi *= np.sqrt(Lambda2) # ! Unsure about this scaling
# g_phi *= Lambda2
# g_phi = g_phi / np.sqrt(Lambda2)
a2 = 7.5412e-4
number_cells = [20, 128]#[20, 100]
Lx = 0.02
Ly = 0.1
x0 = np.array([0, 0])
x1 = np.array([Lx, Ly])
refinement_style = 'uniform'
solvation = 5
PoissonBoltzmann = False
rtol = 1e-3 # ToDo: Change back to 1e-8, currently just for testing
relax_param = 0.15 # 0.1
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')
x, y = X[0], X[1]
y_S = 1 - y_A - y_C
# Plot results
levelsf = 10
levels = 10
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(8,10))
c = axs[0,0].tricontourf(x, y, y_A, levelsf)
fig.colorbar(c)
axs[0,0].tricontour(x, y, y_A, levels, colors='black')
axs[0,0].set_title('$y_A$')
c = axs[0,1].tricontourf(x, y, y_C, levelsf)
fig.colorbar(c)
axs[0,1].tricontour(x, y, y_C, levels, colors='black')
axs[0,1].set_title('$y_C$')
c = axs[0,2].tricontourf(x, y, y_S, levelsf)
fig.colorbar(c)
axs[0,2].tricontour(x, y, y_S, levels, colors='black')
axs[0,2].set_title('$y_S$')
c = axs[1,0].tricontourf(x, y, phi, levelsf)
fig.colorbar(c)
axs[1,0].tricontour(x, y, phi, levels, colors='black')
axs[1,0].set_title('$\\varphi$')
c = axs[1,1].tricontourf(x, y, p, levelsf)
fig.colorbar(c)
axs[1,1].tricontour(x, y, p, levels, colors='black')
axs[1,1].set_title('$p$')
fig.tight_layout()
fig.show()
\ No newline at end of file
build/lib/FENICSxDGM/Eq02.py deleted 100644 → 0
+ 0
− 309
View file @ 511cbc14
'''
Jan Habscheid
Jan.Habscheid@rwth-aachen.de
This file implements the fenics solver for the reduced systme of equations to two equations for the electric potential and the pressure.
'''
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import TestFunctions, split, dot, grad, dx, inner, Mesh, exp
from basix.ufl import element, mixed_element
import matplotlib.pyplot as plt
# Define mesh
def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
'''
Creates a one-dimensional mesh with a refined region at the left boundary
Parameters
----------
refinement_style : str
How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
number_cells : int
Number of cells in the mesh
Returns
-------
Mesh
One-dimensional mesh, ready for use in FEniCSx
'''
if refinement_style == 'log':
coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
elif refinement_style == 'hard_log':
coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
elif refinement_style == 'hard_hard_log':
coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
num_vertices = len(coordinates_np)
num_cells = num_vertices - 1
cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
gdim = 1
shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
degree = 1
domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
coordinates_np_ = []
[coordinates_np_.append([coord]) for coord in coordinates_np]
msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
return msh
def 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, 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):
'''
Solve the simplified dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
System of equations:
λ²Δ φ =−L²n^F
a²∇p=−n^F∇ φ
with the space charge
n^F = z_A y_A(φ, p) + z_C y_C(φ ,p)
if the mixture is compressible:
y_alpha = C_alpha * (K+p−1)^(−κ+1)a²K exp(−z_α φ)
if the mixture is incompressible:
y_alpha = D_alpha * exp(−(κ+1)a²p−z_α φ)
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.
! If the Newton solver diverges, you may try to reduce the relaxation parameter.
Parameters
----------
phi_left : float
Value of φ at the left boundary
phi_right : float
Value of φ at the right boundary
p_right : float
Value of p at the right boundary
z_A : float
Charge number of species A
z_C : float
Charge number of species C
y_A_R : float
Atomic fractions of species A at right boundary
y_C_R : float
Atomic fractions of species C at right boundary
K : float | str
Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
number_cells : int
Number of cells in the mesh
solvation : float, optional
solvation number, by default 0
relax_param : float, optional
Relaxation parameter for the Newton solver
xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
, by default None -> Determined automatically
x0 : float, optional
Left boundary of the domain, by default 0
x1 : float, optional
Right boundary of the domain, by default 1
refinement_style : str, optional
Specify for refinement towards zero
Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
return_type : str, optional
'Vector' or 'Scalar', 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Scalar'
rtol : float, optional
Relative tolerance for Newton solver, by default 1e-8
max_iter : float, optional
Maximum number of Newton iterations, by default 500
Returns
-------
y_A, y_C, phi, p, msh
Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
If return_type is 'Vector', the solution is returned as numpy arrays
'''
if return_type == 'Scalar':
raise NotImplementedError('Scalar return type is not implemented yet')
# Define boundaries of the domain
x0 = 0
x1 = 1
# Define boundaries for the boundary conditions
def Left(x):
return np.isclose(x[0], x0)
def Right(x):
return np.isclose(x[0], x1)
# Create mesh
if refinement_style == 'uniform':
msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
else:
msh = create_refined_mesh(refinement_style, number_cells)
# Define Finite Elements
CG1_elem = element('Lagrange', msh.basix_cell(), 1)
# Define Mixed Function Space
W_elem = mixed_element([CG1_elem, CG1_elem])#, CG1_elem, CG1_elem])
W = fem.functionspace(msh, W_elem)
# Define Trial- and Testfunctions
u = fem.Function(W)
phi, p = split(u)
(v_1, v_2) = TestFunctions(W)
# Collapse function space for bcs
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()
# Define boundary conditions values
def phi_left_(x):
return np.full_like(x[0], phi_left)
def phi_right_(x):
return np.full_like(x[0], phi_right)
def p_right_(x):
return np.full_like(x[0], p_right)
# Interpolate bcs functions
phi_left_bcs = fem.Function(W0)
phi_left_bcs.interpolate(phi_left_)
phi_right_bcs = fem.Function(W0)
phi_right_bcs.interpolate(phi_right_)
p_right_bcs = fem.Function(W1)
p_right_bcs.interpolate(p_right_)
# Identify dofs for boundary conditions
# Define boundary conditions
facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
# Combine boundary conditions into list
bcs = [bc_left_phi, bc_right_phi, bc_right_p]
def y_A(phi, p):
D_A = y_A_R / exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
return D_A * exp(-(solvation + 1) * a2 * p - z_A * phi)
def y_C(phi, p):
D_C = y_C_R / exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
return D_C * exp(-(solvation + 1) * a2 * p - z_C * phi)
# Define variational problem
if K == 'incompressible':
# total free charge density
def nF(y_A, y_C, p):
return (z_C * y_C + z_A * y_A)
else:
# total number density
def n(p):
return (p-1)/K + 1
# total free charge density
def nF(y_A, y_C, p):
return (z_C * y_C + z_A * y_A) * n(p)
# Variational Form
A = (
inner(grad(phi), grad(v_1)) * dx
- 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p), p) * v_1 * dx
) + (
inner(grad(p), grad(v_2)) * dx
+ 1 / a2 * nF(y_A(phi, p), y_C(phi, p), p) * dot(grad(phi), grad(v_2)) * dx
)
F = A
# Define Nonlinear Problem
problem = NonlinearProblem(F, u, bcs=bcs)
# Define Newton Solver and solver settings
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = rtol
if relax_param != None:
solver.relaxation_parameter = relax_param
else:
if phi_right == phi_left:
solver.relaxation_parameter = 1.0
else:
solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
solver.max_it = max_iter
solver.report = True
# Solve the problem
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(u)
assert (converged)
print(f"Number of interations: {n:d}")
# Split the mixed function space into the individual components
phi, p = u.split()
# Return the solution
if return_type=='Vector':
x_vals = np.array(msh.geometry.x[:,0])
phi_vals = np.array(u.sub(0).collapse().x.array)
p_vals = np.array(u.sub(1).collapse().x.array)
# Calculate the atomic fractions
D_A = y_A_R / np.exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
y_A_vals = D_A * np.exp(-(solvation + 1) * a2 * p_vals - z_A * phi_vals)
D_C = y_C_R / np.exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
y_C_vals = D_C * np.exp(-(solvation + 1) * a2 * p_vals - z_C * phi_vals)
return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
if __name__ == '__main__':
# Define the parameters
phi_left = 5.0
phi_right = 0.0
p_right = 0.0
y_A_R = 1/3
y_C_R = 1/3
z_A = -1.0
z_C = 1.0
K = 'incompressible'
Lambda2 = 8.553e-6
a2 = 7.5412e-4
number_cells = 1024
relax_param = .1
rtol = 1e-4
max_iter = 500
# Solve the system
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)
# Plot the solution
plt.plot(x, phi)
plt.xlim(0,0.05)
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$\\varphi$ [-]')
plt.show()
plt.plot(x, y_A, '--', color='tab:blue', label='$y_A$')
plt.plot(x, y_C, '-', color='tab:blue', label='$y_C$')
plt.plot(x, 1 - y_A - y_C, ':', color='tab:blue', label='$y_S$')
plt.xlim(0,0.05)
plt.legend()
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$y_\\alpha$ [-]')
plt.show()
plt.plot(x, p)
plt.xlim(0,0.05)
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$p$ [-]')
plt.show()
build/lib/FENICSxDGM/Eq04.py deleted 100644 → 0
+ 0
− 358
View file @ 511cbc14
'''
Jan Habscheid
Jan.Habscheid@rwth-aachen.de
This script implements the fenics solver for a ternary electrolyte (N=3, A,C,S)
'''
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
from basix.ufl import element, mixed_element
import matplotlib.pyplot as plt
# Define mesh
def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
'''
Creates a one-dimensional mesh with a refined region at the left boundary
Parameters
----------
refinement_style : str
How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
number_cells : int
Number of cells in the mesh
Returns
-------
Mesh
One-dimensional mesh, ready for use in FEniCSx
'''
if refinement_style == 'log':
coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
elif refinement_style == 'hard_log':
coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
elif refinement_style == 'hard_hard_log':
coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
num_vertices = len(coordinates_np)
num_cells = num_vertices - 1
cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
gdim = 1
shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
degree = 1
domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
coordinates_np_ = []
[coordinates_np_.append([coord]) for coord in coordinates_np]
msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
return msh
def 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):
'''
Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
System of equations:
λ²Δ φ =−L²n^F
a²∇p=−n^F∇ φ
div(J_A)=0
div(J_C)=0
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.
! If the Newton solver diverges, you may try to reduce the relaxation parameter.
Parameters
----------
phi_left : float
Value of φ at the left boundary
phi_right : float
Value of φ at the right boundary
p_right : float
Value of p at the right boundary
z_A : float
Charge number of species A
z_C : float
Charge number of species C
y_A_R : float
Atomic fractions of species A at right boundary
y_C_R : float
Atomic fractions of species C at right boundary
K : float | str
Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
number_cells : int
Number of cells in the mesh
solvation : float, optional
solvation number, by default 0
PoissonBoltzmann : bool, optional
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
relax_param : float, optional
Relaxation parameter for the Newton solver
xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
, by default None -> Determined automatically
x0 : float, optional
Left boundary of the domain, by default 0
x1 : float, optional
Right boundary of the domain, by default 1
refinement_style : str, optional
Specify for refinement towards zero
Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
return_type : str, optional
'Vector' or 'Scalar', 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Scalar'
rtol : float, optional
Relative tolerance for Newton solver, by default 1e-8
max_iter : float, optional
Maximum number of Newton iterations, by default 500
Returns
-------
y_A, y_C, phi, p, msh
Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
If return_type is 'Vector', the solution is returned as numpy arrays
'''
# Define boundaries of the domain
x0 = 0
x1 = 1
# Define boundaries for the boundary conditions
def Left(x):
return np.isclose(x[0], x0)
def Right(x):
return np.isclose(x[0], x1)
# Create mesh
if refinement_style == 'uniform':
msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
else:
msh = create_refined_mesh(refinement_style, number_cells)
# Define Finite Elements
CG1_elem = element('Lagrange', msh.basix_cell(), 1)
# Define Mixed Function Space
W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
W = fem.functionspace(msh, W_elem)
# Define Trial- and Testfunctions
u = fem.Function(W)
y_A, y_C, phi, p = split(u)
(v_A, v_C, v_1, v_2) = TestFunctions(W)
# Collapse function space for bcs
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()
W2, _ = W.sub(2).collapse()
W3, _ = W.sub(3).collapse()
# Define boundary conditions values
def phi_left_(x):
return np.full_like(x[0], phi_left)
def phi_right_(x):
return np.full_like(x[0], phi_right)
def p_right_(x):
return np.full_like(x[0], p_right)
def y_A_right_(x):
return np.full_like(x[0], y_A_R)
def y_C_right_(x):
return np.full_like(x[0], y_C_R)
# Interpolate bcs functions
phi_left_bcs = fem.Function(W2)
phi_left_bcs.interpolate(phi_left_)
phi_right_bcs = fem.Function(W2)
phi_right_bcs.interpolate(phi_right_)
p_right_bcs = fem.Function(W3)
p_right_bcs.interpolate(p_right_)
y_A_right_bcs = fem.Function(W0)
y_A_right_bcs.interpolate(y_A_right_)
y_C_right_bcs = fem.Function(W1)
y_C_right_bcs.interpolate(y_C_right_)
# Identify dofs for boundary conditions
# Define boundary conditions
facet_left_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Left)
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Right)
bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(2))
bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(2))
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Right)
bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(3))
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
bc_right_y_A = fem.dirichletbc(y_A_right_bcs, facet_right_dofs, W.sub(0))
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
bc_right_y_C = fem.dirichletbc(y_C_right_bcs, facet_right_dofs, W.sub(1))
# Combine boundary conditions into list
bcs = [bc_left_phi, bc_right_phi, bc_right_p, bc_right_y_A, bc_right_y_C]
# Define variational problem
if K == 'incompressible':
# total free charge density
def nF(y_A, y_C):
return (z_C * y_C + z_A * y_A)
# Diffusion fluxes for species A and C
def J_A(y_A, y_C, phi, p):
return grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
# return grad(ln(y_A) + a2 * (p - 1) - solvation * ln(1-y_A-y_C) + z_A * phi)
def J_C(y_A, y_C, phi, p):
return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
# return grad(ln(y_C) + a2 * (p - 1) - solvation * ln(1-y_A-y_C) + z_C * phi)
# Variational Form
A = (
inner(grad(phi), grad(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(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 += (
inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_A)) * dx
+ inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_C)) * dx
)
else:
# total number density
def n(p):
return (p-1)/K + 1
# total free charge density
def nF(y_A, y_C, p):
return (z_C * y_C + z_A * y_A) * n(p)
# Diffusion fluxes for species A and C
def J_A(y_A, y_C, phi, p):
return ln(y_A) + a2 * (solvation + 1) * K * ln(1 + 1/K * (p-1)) + z_A * phi
def J_C(y_A, y_C, phi, p):
return ln(y_C) + a2 * (solvation + 1)* K * ln(1 + 1/K * (p-1)) + z_C * phi
A = (
inner(grad(phi), grad(v_1)) * dx
- 1 / Lambda2 * nF(y_A, y_C, p) * v_1 * dx
) + (
inner(grad(p), grad(v_2)) * dx
+ 1 / a2 * nF(y_A, y_C, p) * 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
)
F = A
# Initialize initial guess for u
y_C_init = fem.Function(W1)
y_A_init = fem.Function(W0)
y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_R))
y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_R))
with u.vector.localForm() as u_loc:
u_loc.set(0)
u.sub(0).interpolate(y_A_init)
u.sub(1).interpolate(y_C_init)
# Define Nonlinear Problem
problem = NonlinearProblem(F, u, bcs=bcs)
# Define Newton Solver and solver settings
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = rtol
if relax_param != None:
solver.relaxation_parameter = relax_param
else:
if phi_right == phi_left:
solver.relaxation_parameter = 1.0
else:
solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
solver.max_it = max_iter
solver.report = True
# Solve the problem
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(u)
assert (converged)
print(f"Number of interations: {n:d}")
# Split the mixed function space into the individual components
y_A, y_C, phi, p = u.split()
# Return the solution
if return_type=='Vector':
x_vals = np.array(msh.geometry.x[:,0])
y_A_vals = np.array(u.sub(0).collapse().x.array)
y_C_vals = np.array(u.sub(1).collapse().x.array)
phi_vals = np.array(u.sub(2).collapse().x.array)
p_vals = np.array(u.sub(3).collapse().x.array)
return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
elif return_type=='Scalar':
return y_A, y_C, phi, p, msh
if __name__ == '__main__':
# Define the parameters
phi_left = 10.0
phi_right = 0.0
p_right = 0.0
y_A_R = 0.01#1/3
y_C_R = 0.01#1/3
z_A = -1.0
z_C = 1.0
K = 'incompressible'
Lambda2 = 8.553e-6
a2 = 7.5412e-4
solvation = 0
number_cells = 1024
relax_param = .1
rtol = 1e-4
max_iter = 500
# Solve the system
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)
# Plot the solution
plt.plot(x, phi)
plt.xlim(0,0.05)
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$\\varphi$ [-]')
plt.show()
plt.plot(x, y_A, '--', color='tab:blue', label='$y_A$')
plt.plot(x, y_C, '-', color='tab:blue', label='$y_C$')
plt.plot(x, 1 - y_A - y_C, ':', color='tab:blue', label='$y_S$')
plt.xlim(0,0.05)
plt.legend()
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$y_\\alpha$ [-]')
plt.show()
plt.plot(x, p)
plt.xlim(0,0.05)
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$p$ [-]')
plt.show()
build/lib/FENICSxDGM/EqN.py deleted 100644 → 0
+ 0
− 324
View file @ 511cbc14
'''
Jan Habscheid
Jan.Habscheid@rwth-aachen.de
This script implements the fenics solver for the generic system of equations for N species
'''
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
from basix.ufl import element, mixed_element
import matplotlib.pyplot as plt
# Define mesh
def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
'''
Creates a one-dimensional mesh with a refined region at the left boundary
Parameters
----------
refinement_style : str
How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
number_cells : int
Number of cells in the mesh
Returns
-------
Mesh
One-dimensional mesh, ready for use in FEniCSx
'''
if refinement_style == 'log':
coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
elif refinement_style == 'hard_log':
coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
elif refinement_style == 'hard_hard_log':
coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
num_vertices = len(coordinates_np)
num_cells = num_vertices - 1
cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
gdim = 1
shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
degree = 1
domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
coordinates_np_ = []
[coordinates_np_.append([coord]) for coord in coordinates_np]
msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
return msh
def 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):
'''
Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
System of equations:
λ²Δ φ =−L²n^F
a²∇p=−n^F∇ φ
div(J_α)=0 α∈ {1,...,N−1}
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.
! If the Newton solver diverges, you may try to reduce the relaxation parameter.
Parameters
----------
phi_left : float
Value of φ at the left boundary
phi_right : float
Value of φ at the right boundary
p_right : float
Value of p at the right boundary
z_alpha : list
Charge numbers for species α = 1,...,N-1
y_R : list
Atomic fractions at right boundary for species α = 1,...,N-1
K : float | str
Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
number_cells : int
Number of cells in the mesh
solvation : float, optional
solvation number, not implemented yet, by default 0
PoissonBoltzmann : bool, optional
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
relax_param : float, optional
Relaxation parameter for the Newton solver
xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
, by default None -> Determined automatically
x0 : float, optional
Left boundary of the domain, by default 0
x1 : float, optional
Right boundary of the domain, by default 1
refinement_style : str, optional
Specify for refinement towards zero
Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
return_type : str, optional
'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'
rtol : float, optional
Relative tolerance for Newton solver, by default 1e-8
max_iter : float, optional
Maximum number of Newton iterations, by default 500
Returns
-------
y_A, y_C, phi, p, msh
Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
If return_type is 'Vector', the solution is returned as numpy arrays
Only return_type 'Vector' is implemented yet
'''
if solvation != 0:
raise NotImplementedError('Solvation number not implemented yet')
if PoissonBoltzmann:
raise NotImplementedError('Poisson-Boltzmann not implemented yet')
# Define boundaries of the domain
x0 = 0
x1 = 1
# Define boundaries
def Left(x):
return np.isclose(x[0], x0)
def Right(x):
return np.isclose(x[0], x1)
# Create mesh
if refinement_style == 'uniform':
msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
else:
msh = create_refined_mesh(refinement_style, number_cells)
# Define Finite Elements
CG1_elem = element('Lagrange', msh.basix_cell(), 1)
# Define Mixed Function Space
Elem_list = [CG1_elem, CG1_elem]
[Elem_list.append(CG1_elem) for _ in range(len(z_alpha))]
W_elem = mixed_element(Elem_list)
W = fem.functionspace(msh, W_elem)
# Define Trial- and Testfunctions
u = fem.Function(W)
my_TrialFunctions = split(u)
my_TestFunctions = TestFunctions(W)
phi, p = my_TrialFunctions[0], my_TrialFunctions[1]
v_1, v_2 = my_TestFunctions[0], my_TestFunctions[1]
y_alpha = my_TrialFunctions[2:]
v_alpha = my_TestFunctions[2:]
# Collapse function space for bcs
W_ = []
[W_.append(W.sub(i).collapse()[0]) for i in range(len(z_alpha)+2)]
# Define boundary conditions values
def phi_left_(x):
return np.full_like(x[0], phi_left)
def phi_right_(x):
return np.full_like(x[0], phi_right)
def p_right_(x):
return np.full_like(x[0], p_right)
# Interpolate bcs functions
phi_left_bcs = fem.Function(W_[0])
phi_left_bcs.interpolate(phi_left_)
phi_right_bcs = fem.Function(W_[0])
phi_right_bcs.interpolate(phi_right_)
p_right_bcs = fem.Function(W_[1])
p_right_bcs.interpolate(p_right_)
# Identify dofs for boundary conditions
facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
# Combine boundary conditions for electric potential and pressure into list
bcs = [bc_left_phi, bc_right_phi, bc_right_p]
# Repeat the same for the boundary conditoins for the atomic fractions
for i in range(len(z_alpha)):
y_right_bcs = fem.Function(W_[i+2])
def y_right_(x):
return np.full_like(x[0], y_R[i])
y_right_bcs.interpolate(y_right_)
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(i+2), W.sub(i+2).collapse()[0]), Right)
bc_right_y = fem.dirichletbc(y_right_bcs, facet_right_dofs, W.sub(i+2))
bcs.append(bc_right_y)
# Define variational problem
if K == 'incompressible':
# total free charge density
def nF(y_alpha):
nF = 0
for i in range(len(z_alpha)):
nF += z_alpha[i] * y_alpha[i]
return nF
# Diffusion fluxes for species A and C
def J_alpha(y_alpha, alpha, phi, p):
mu_alpha = ln(y_alpha[alpha])
mu_S = ln(1 - sum(y_alpha))
return grad(mu_alpha - mu_S + z_alpha[alpha] * phi)
# Variational Form
A = (
inner(grad(phi), grad(v_1)) * dx
- 1 / Lambda2 * nF(y_alpha) * v_1 * dx
) + (
inner(grad(p), grad(v_2)) * dx
+ 1 / a2 * nF(y_alpha) * dot(grad(phi), grad(v_2)) * dx
)
for alpha in range(len(z_alpha)):
A += (
inner(J_alpha(y_alpha, alpha, phi, p), grad(v_alpha[alpha])) * dx
)
if PoissonBoltzmann:
raise ValueError('Poisson-Boltzmann not implemented for incompressible systems')
else:
raise ValueError('Only incompressible systems are implemented')
F = A
# Initialize initial guess for u
with u.vector.localForm() as u_loc:
u_loc.set(0)
# Initialize initial guess for u
for alpha in range(len(z_alpha)):
y_alpha_init = fem.Function(W_[alpha+2])
y_alpha_init.interpolate(lambda x: np.full_like(x[0], y_R[alpha]))
u.sub(alpha+2).interpolate(y_alpha_init)
# Define Nonlinear Problem
problem = NonlinearProblem(F, u, bcs=bcs)
# Define Newton Solver
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = rtol
if relax_param != None:
solver.relaxation_parameter = relax_param
else:
if phi_right == phi_left:
solver.relaxation_parameter = 1.0
else:
solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
solver.max_it = max_iter
solver.report = True
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(u)
assert (converged)
print(f"Number of (interations: {n:d}")
# Return the solution
if return_type=='Vector':
x = np.array(msh.geometry.x[:,0])
phi = np.array(u.sub(0).collapse().x.array)
p = np.array(u.sub(1).collapse().x.array)
y = []
[y.append(u.sub(i+2).collapse().x.array) for i in range(len(z_alpha))]
y = np.array(y)
return y, phi, p, x
if __name__ == '__main__':
# Define the parameters
phi_left = 8.0
phi_right = 0.0
p_right = 0.0
y_R = [3/6, 1/6, 1/6]
z_alpha = [-1.0, 1.0, 2.0]
K = 'incompressible'
Lambda2 = 8.553e-6
a2 = 7.5412e-4
number_cells = 1024
relax_param = .05
rtol = 1e-4
max_iter = 2_500
refinement_style = 'hard_log'
return_type = 'Vector'
# Solve the system
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)
# Plot the solution
plt.figure()
plt.plot(x, phi)
plt.xlim(0,0.05)
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$\\varphi$ [-]')
plt.show()
plt.figure()
plt.xlim(0,0.05)
plt.plot(x, p)
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$p$ [-]')
plt.show()
plt.figure()
for i in range(len(z_alpha)):
plt.plot(x, y[i], label=f'$y_{i}$')
plt.xlim(0,0.05)
plt.legend()
plt.grid()
plt.xlabel('x [-]')
plt.ylabel('$y_i$ [-]')
plt.show()
\ No newline at end of file
build/lib/FENICSxDGM/Helper_DoubleLayerCapacity.py deleted 100644 → 0
+ 0
− 411
View file @ 511cbc14
'''
Jan Habscheid
Jan.Habscheid@rwth-aachen.de
'''
import numpy as np
from scipy.optimize import fixed_point, fsolve
# Helper functions
def Phi_pot_center(Phi_pot:np.ndarray) -> np.ndarray:
'''
Returns vector with the center of the electric potential values.
Parameters
----------
Phi_pot : np.ndarray
Input vector with electric potential values.
Returns
-------
np.ndarray
Vector with the center of the electric potential values.
'''
return (Phi_pot[1:] + Phi_pot[:-1]) / 2
def dx(Phi_pot:np.ndarray) -> float:
'''
Returns the difference between the first two electric potential values.
Assumes that the electric potential values are equally spaced.
Parameters
----------
Phi_pot : np.ndarray
Input vector with electric potential values.
Returns
-------
float
Difference between the first two electric potential values.
'''
return Phi_pot[1] - Phi_pot[0]
def C_dl(Q_DL:np.ndarray, Phi_pot:np.ndarray) -> np.ndarray:
'''
Double Layer Capacity
Parameters
----------
Q_DL : np.ndarray
Charge of the system
Phi_pot : np.ndarray
Electric potential values
Returns
-------
np.ndarray
Double Layer Capacity
'''
return (Q_DL[1:] - Q_DL[:-1]) / dx(Phi_pot)
def n(p:np.ndarray, K:str|float) -> np.ndarray:
'''
Calculates the total number density
Parameters
----------
p : np.ndarray
Pressure
K : str | float
Bulk modulus, use 'incompressible' for an incompressible mixture
Returns
-------
np.ndarray
Total number density
'''
if K == 'incompressible':
return np.ones_like(p)
n_Dimensionless = (p-1) / K + 1
return n_Dimensionless
def Q_num_(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float=-1.0, z_C:float=1.0) -> float:
'''
Calculates the charge of the system
Q = ∫_Ω n^F dΩ
Parameters
----------
y_A : np.ndarray
Anion fraction
y_C : np.ndarray
Cation fraction
n : np.ndarray
Total number density
x : np.ndarray
Spatial discretization
z_A : float, optional
Charge number of anions, by default -1.0
z_C : float, optional
Charge number of cations, by default 1.0
Returns
-------
float
Charge of the system
'''
nF_dimensionless = (z_C * y_C + z_A * y_A) * n
nF_int = -np.trapz(nF_dimensionless, x)
return nF_int
def Q_num_dim(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float, z_C:float, nR_m:float, e0:float, LR:float) -> float:
'''
Calculates the charge of the system
Q = ∫_Ω n^F dΩ
Parameters
----------
y_A : np.ndarray
Anion fraction
y_C : np.ndarray
Cation fraction
n : np.ndarray
Total number density
x : np.ndarray
Spatial discretization
z_A : float, optional
Charge number of anions, by default -1.0
z_C : float, optional
Charge number of cations, by default 1.0
nR_m : float
Reference number density in 1/m^3
e0 : float
Dielectric constant
LR : float
Reference length in m
Returns
-------
float
Charge of the system
'''
Q_DL = Q_num_(y_A, y_C, n, x, z_A, z_C)
Q_DL *= nR_m * e0 * LR
Q_DL *= 1e+6
Q_DL *= 1/(1e+4)
return Q_DL
def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
'''
Calculates charge of the system using the analytical method in dimensionless units
Q = λ^2(∂ₓ φᴸ − ∂ₓ φᴿ)
Parameters
----------
y_A_R : float
Value of anion fraction at the right boundary
y_C_R : float
Value of cation fraction at the right boundary
y_N_R : float
Value of neutral fraction at the right boundary
z_A : float
Charge number of anions
z_C : float
Charge number of cations
phi_L : float
Electric potential at the left boundary [V]
phi_R : float
Electric potential at the right boundary
p_R : float
Pressure at the right boundary
K : str | float
Bulk modulus, use 'incompressible' for an incompressible mixture
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
solvation : float
Solvation number
Returns
-------
float
Charge of the system in dimensionless units
'''
if solvation != 0:
raise ValueError('Solvation number must be 0 for the analytical method')
z_N = 0
E_p = p_R
if K == 'incompressible':
# Assume E_p = p_R = 0
D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
D_N = y_N_R / (np.exp(-a2*p_R))
# Analytical method
E_p = p_R
CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - E_p) + D_C * np.exp(-z_C * phi_L - E_p) + D_N * np.exp(-z_N * phi_L - E_p))
Lambda = np.sqrt(Lambda2)
Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2) * (np.sqrt(CLambda_L))
return Q_DL
else:
C_A = y_A_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
C_C = y_C_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
C_N = y_N_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
CLambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N
Left = np.sqrt((1/CLambda_L)**(-1/(a2*K)) + 1 - K - E_p)
Lambda = np.sqrt(Lambda2)
a = np.sqrt(a2)
Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * a * np.sqrt(2) * (Left)# - Right)
return Q_DL
raise ValueError('Invalid input for K')
def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, solvation:float) -> float:
'''
Calculates charge of the system using the analytical method in dimensionless units
Q = λ^2(∂ₓ φᴸ − ∂ₓ φᴿ)
Parameters
----------
y_A_R : float
Value of anion fraction at the right boundary
y_C_R : float
Value of cation fraction at the right boundary
y_N_R : float
Value of neutral fraction at the right boundary
z_A : float
Charge number of anions
z_C : float
Charge number of cations
phi_L : float
Value of electric potential at the left boundary (dimensionless units)
phi_R : float
Value of electric potential at the right
p_R : float
Value of pressure at the right boundary
K : str | float
Bulk modulus, use 'incompressible' for an incompress
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
nR_m : float
Reference number density in 1/m^3
e0 : float
Dielectric constant
LR : float
Reference length in m
solvation : float
Solvation number
Returns
-------
float
Charge of the system in µAs/cm³
'''
Q_DL = Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
Q_DL *= nR_m * e0 * LR
Q_DL *= 1e+6
Q_DL *= 1/(1e+4)
return Q_DL
def C_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
'''
Calculates the double layer capacity of the system using the analytical method in dimensionless units
C_dl = ∂Q/∂ φᴸ
Parameters
----------
y_A_R : float
Value of anion fraction at the right boundary
y_C_R : float
Value of cation fraction at the right boundary
y_N_R : float
Value of neutral fraction at the right boundary
z_A : float
Charge number of anions
z_C : float
Charge number of cations
phi_L : float
Electric potential at the left boundary [V]
phi_R : float
Electric potential at the right boundary
p_R : float
Pressure at the right boundary
K : str | float
Bulk modulus, use 'incompressible' for an incompressible mixture
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
solvation : float
Solvation number
Returns
-------
float
Double layer capacity of the system in dimensionless units
'''
if solvation != 0:
raise ValueError('Solvation number must be 0 for the analytical method')
z_N = 0
E_p = p_R
if K == 'incompressible':
# Assume E_p = p_R = 0
D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
D_N = y_N_R / (np.exp(-a2*p_R))
CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - E_p * a2) + D_C * np.exp(-z_C * phi_L - E_p * a2) + D_N * np.exp(-z_N * phi_L - E_p * a2))
xi_CLambda_L = D_A * np.exp(-z_A * phi_L - E_p * a2) + D_C * np.exp(-z_C * phi_L - E_p * a2) + D_N * np.exp(-z_N * phi_L - E_p * a2)
dxi_CLambda_L = - D_A * z_A * np.exp(-z_A * phi_L - E_p * a2) - D_C * z_C * np.exp(-z_C * phi_L - E_p * a2) - D_N * z_N * np.exp(-z_N * phi_L - E_p * a2)
Lambda = np.sqrt(Lambda2)
PreTerm = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2)
DerivedTerm = dxi_CLambda_L / (2 * xi_CLambda_L * np.sqrt(np.log(xi_CLambda_L)))
return PreTerm * DerivedTerm
else:
C_A = y_A_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
C_C = y_C_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
C_N = y_N_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
z = 1 - K - E_p
y = 1/(a2*K)
CLambda = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
dClambda = -C_A * z_A * np.exp(-z_A*phi_L) - C_C * z_C * np.exp(-z_C*phi_L) - C_N * z_N * np.exp(-z_N*phi_L)
Lambda = np.sqrt(Lambda2)
a = np.sqrt(a2)
DerivedTerm = y * CLambda**(y-1) * dClambda / (2 * np.sqrt(CLambda**y + z))
PreTerm = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * a * np.sqrt(2)
return PreTerm * DerivedTerm
raise ValueError('Invalid input for K')
def C_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, k:float, T:float, solvation:float) -> float:
'''
Calculates the double layer capacity of the system using the analytical method in dimensionless units
C_dl = ∂Q/∂ φᴸ
Parameters
----------
y_A_R : float
Value of anion fraction at the right boundary
y_C_R : float
Value of cation fraction at the right boundary
y_N_R : float
Value of neutral fraction at the right boundary
z_A : float
Charge number of anions
z_C : float
Charge number of cations
phi_L : float
Value of electric potential at the left boundary (dimensionless units)
phi_R : float
Value of electric potential at the right
p_R : float
Value of pressure at the right boundary
K : str | float
Bulk modulus, use 'incompressible' for an incompress
Lambda2 : float
Dimensionless parameter
a2 : float
Dimensionless parameter
nR_m : float
Reference number density in 1/m^3
e0 : float
Dielectric constant
LR : float
Reference length in m
k : float
Boltzmann constant
T : float
Temperature
solvation : float
Solvation number
Returns
-------
float
Double layer capacity of the system in µAs/cm²
'''
C_DL = C_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
C_DL *= nR_m * e0 * LR
# ! ToDo
C_DL *= 1e+6
C_DL *= 1/(1e+4)
C_DL *= 1 /(k*T/e0)
return C_DL
\ No newline at end of file
build/lib/FENICSxDGM/__init__.py deleted 100644 → 0
+ 0
− 5
View file @ 511cbc14
from .ElectrolyticDiode import ElectrolyticDiode
from .Eq02 import solve_System_2eq, create_refined_mesh
from .Eq04 import solve_System_4eq
from .EqN import solve_System_Neq
from .Helper_DoubleLayerCapacity import *
\ No newline at end of file
dist/FENICSxDGM-1.0-py3-none-any.whl deleted 100644 → 0
+ 0
− 0
View file @ 511cbc14
File deleted
dist/fenicsxdgm-1.0.linux-x86_64.tar.gz deleted 100644 → 0
+ 0
− 0
View file @ 511cbc14
File deleted
dist/fenicsxdgm-1.0.tar.gz deleted 100644 → 0
+ 0
− 0
View file @ 511cbc14
File deleted
setup.py
+ 0
− 3
View file @ 8d9097ee
...@@ -8,11 +8,8 @@ setup( ...@@ -8,11 +8,8 @@ setup(
author="Jan Habscheid, Lambert Theisen, Manuel Torrilhon", author="Jan Habscheid, Lambert Theisen, Manuel Torrilhon",
author_email="Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de", author_email="Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de",
packages=["FENICSxDGM"], packages=["FENICSxDGM"],
# packages=find_packages("src"),
package_dir={"FENICSxDGM":"src"}, package_dir={"FENICSxDGM":"src"},
# package_dir={"": "src"},
package_data={"FENICSxDGM": ["tests/*"]}, package_data={"FENICSxDGM": ["tests/*"]},
# package_data={"": ["tests/*"]},
python_requires=">=3.12.3", python_requires=">=3.12.3",
install_requires=[ install_requires=[
"pyvista == 0.43.10", "pyvista == 0.43.10",
......
src/FENICSxDGM.egg-info/PKG-INFO deleted 100644 → 0
+ 0
− 117
View file @ 511cbc14
Metadata-Version: 2.1
Name: FENICSxDGM
Version: 1.0
Summary: A description is yet to follow
Home-page: https://git.rwth-aachen.de/JanHab/Bsc-ElectrolyteModels
Author: Jan Habscheid, Lambert Theisen, Manuel Torrilhon
Author-email: Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de
License: GNU General Public License v3.0 or later
Keywords: fenics-project,Battery Simulation,finite element method
Classifier: License :: OSI Approved :: GNU General Public License v3 (GPLv3)
Classifier: Programming Language :: Python :: 3
Requires-Python: >=3.12.3
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: pyvista==0.43.10
Requires-Dist: numpy==1.26.4
Requires-Dist: scipy==1.14.0
Provides-Extra: notebook
Requires-Dist: jupyter>=1.0.0; extra == "notebook"
Provides-Extra: tests
Requires-Dist: pytest==8.3.3; extra == "tests"
Provides-Extra: docs
Requires-Dist: sphinx==7.3.7; extra == "docs"
Requires-Dist: myst-parser==4.0.0; extra == "docs"
Requires-Dist: sphinx-copybotton==0.5.2; extra == "docs"
Requires-Dist: sphinx_rtd_theme==3.0.1; extra == "docs"
Provides-Extra: build
Requires-Dist: build>=0.7.0; extra == "build"
Provides-Extra: dev
Requires-Dist: FENICSxDGM[tests]; extra == "dev"
Requires-Dist: FENICSxDGM[docs]; extra == "dev"
Requires-Dist: FENICSxDGM[build]; extra == "dev"
# Reproducibility Repository for Numerical Treatment of a Thermodynamically Consistent Electrolyte Model (B.Sc. Thesis - Jan Habscheid)
[![Pipeline Status](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/badges/main/pipeline.svg)](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/pipelines)
[![Documentation](https://img.shields.io/badge/docs-latest-blue)](https://janhab.pages.rwth-aachen.de/bsc-electrolytemodels/)
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.13645296.svg)](https://doi.org/10.5281/zenodo.13645296)
[![GitLab Version](https://img.shields.io/badge/version-1.0-blue.svg)](https://git.rwth-aachen.de/jan.habscheid/bsc-electrolytemodels/-/tags)
[![License](https://img.shields.io/badge/license-GPLv3-blue)](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/-/blob/main/LICENSE?ref_type=heads)
## Thesis
This repository contains the code to reproduce the results presented in the bachelor thesis: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model
Find the thesis at [https://doi.org/10.18154/RWTH-2024-09837](https://doi.org/10.18154/RWTH-2024-09837)
### Abstract
Batteries play a crucial role in the energy transition. The production of green energy depends on external factors. Storing energy in batteries is necessary to access green energy at any time.
Better optimized batteries are essential for the future. Lifetime, loading time, and energy loss are just some aspects that must be improved to prepare for a greener future. Numerical simulations are crucial to understanding and optimizing batteries' behavior. Those simulations enable researchers to test many different materials without considerable additional expenses to, for example, find the best combination of anions and cations.
The classical Nernst-Planck model for the ion transport in an electrolyte fails to predict the correct concentration in the boundaries of the electrolyte. This work will present and analyze a thermodynamically consistent electrolyte model with dimensionless units under isothermal conditions. A simplified version of the system for the one-dimensional equilibrium of an ideal mixture and the incompressible limit will be considered. The numerical implementation of the model with the open-source software FEniCSx will be discussed. Furthermore, the influence of different boundary conditions, material parameters, solvation, and compressibility on the electric potential, pressure, and ion concentration will be investigated, and the model will be compared with the classical Nernst-Planck model. Examples of the double layer capacity and electrolytic diode will be considered.
## Installation
As a numerical solver, mainly FEniCSx was used and installed via conda.
All the calculations were performed on a Linux machine. According to the documentation, everything should work well on macOS, but this was not tested. FEniCSx offers some beta versions for Windows support, but it is recommended to use WSL2 instead.
```
conda create --name fenicsx-env python=3.12.3 -y
conda activate fenicsx-env
conda install -c conda-forge fenics-dolfinx=0.8.0 mpich=4.2.1 pyvista=0.43.10 gcc=12.4.0 matplotlib=3.8.4 numpy=1.26.4 scipy=1.14.0 pytest==8.3.3 -y
```
### Alternative installation
Use the "environment.yml" file to install all necessary environments
```
conda env create -f environment.yml
```
### macOS installation using Docker
```
docker compose build
docker compose run solver
```
### Testing
Use pytest with
```
python -m pytest
```
to verify that everything was installed correctly.
## Usage
Find the visualizations from the thesis and some extra calculations in the "examples" folder.
In the subfolder "ReproducableCode" is the code, to execute the calculations with some first visualizations.
The subfolder "Data" stores the data for all the simulations in a *.npz file, which can be read with numpy `np.load(file.npz)`.
"Visualizations" creates the necessary figures from the thesis and stores them in *.svg format in "Figures".
In "src" there are the generic FEniCSx implementations, that were used to calculate the examples.
## Contact
**Author**
- Jan Habscheid
- Jan.Habscheid@rwth-aachen.de
**Supervisor**
- Dr. Lambert Theissen
- ACoM - Applied and Computational Mathematics
- RWTH Aachen University
- theisen@acom.rwth-aachen.de
**Supervisor**
- Prof. Dr. Manuel Torrilhon
- ACoM - Applied and Computational Mathematics
- RWTH Aachen University
- mt@acom.rwth-aachen.de
src/FENICSxDGM.egg-info/SOURCES.txt deleted 100644 → 0
+ 0
− 23
View file @ 511cbc14
LICENSE
README.md
pyproject.toml
setup.cfg
setup.py
src/ElectrolyticDiode.py
src/Eq02.py
src/Eq04.py
src/EqN.py
src/Helper_DoubleLayerCapacity.py
src/__init__.py
src/FENICSxDGM.egg-info/PKG-INFO
src/FENICSxDGM.egg-info/SOURCES.txt
src/FENICSxDGM.egg-info/dependency_links.txt
src/FENICSxDGM.egg-info/entry_points.txt
src/FENICSxDGM.egg-info/requires.txt
src/FENICSxDGM.egg-info/top_level.txt
src/FENICSxDGM.egg-info/zip-safe
tests/test_ElectrolyticDiode.py
tests/test_Eq02.py
tests/test_Eq04.py
tests/test_EqN.py
tests/test_Helper_DoubleLayerCapacity.py
\ No newline at end of file
src/FENICSxDGM.egg-info/dependency_links.txt deleted 100644 → 0
+ 0
− 1
View file @ 511cbc14
src/FENICSxDGM.egg-info/entry_points.txt deleted 100644 → 0
+ 0
− 2
View file @ 511cbc14
[console_scripts]
my-example-utility = example.example_module:main
src/FENICSxDGM.egg-info/requires.txt deleted 100644 → 0
+ 0
− 23
View file @ 511cbc14
pyvista==0.43.10
numpy==1.26.4
scipy==1.14.0
[build]
build>=0.7.0
[dev]
FENICSxDGM[tests]
FENICSxDGM[docs]
FENICSxDGM[build]
[docs]
sphinx==7.3.7
myst-parser==4.0.0
sphinx-copybotton==0.5.2
sphinx_rtd_theme==3.0.1
[notebook]
jupyter>=1.0.0
[tests]
pytest==8.3.3
src/FENICSxDGM.egg-info/top_level.txt deleted 100644 → 0
+ 0
− 1
View file @ 511cbc14
FENICSxDGM
src/FENICSxDGM.egg-info/zip-safe deleted 100644 → 0
+ 0
− 1
View file @ 511cbc14
test_.py deleted 100644 → 0
+ 0
− 40
View file @ 511cbc14
'''
Tests the Eq02 implementation in src.Eq02.py
'''
from src.Eq02 import solve_System_2eq
import numpy as np
# Define the testing tolerance
rtol = 1e-10
atol = 1e-10
# Define parameter to use in the test
phi_left = 10.0
phi_right = 0.0
p_right = 0.0
y_A_R = 0.01
y_C_R = 0.01
z_A = -1.0
z_C = 1.0
K = 'incompressible'
Lambda2 = 8.553e-6
a2 = 7.5412e-4
solvation = 0
number_cells = 32
relax_param = .1
rtol = 1e-4
max_iter = 500
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, solvation=solvation, relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
y_A_NP, y_C_NP, phi_NP, p_NP, x_NP = solve_System_2eq(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, PoissonBoltzmann=True, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
solvation = 5
y_A_solvation, y_C_solvation, phi_solvation, p_solvation, x_solvation = solve_System_2eq(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='log', return_type='Vector', max_iter=max_iter, rtol=rtol)
solvation = 0
K = 5_000
y_A_compressible, y_C_compressible, phi_compressible, p_compressible, x_compressible = solve_System_2eq(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='hard_log', return_type='Vector', max_iter=max_iter, rtol=rtol)
np.savez('tests/TestData/Eq02.npz', y_A=y_A, y_C=y_C, phi=phi, p=p, x=x, y_A_NP=y_A_NP, y_C_NP=y_C_NP, phi_NP=phi_NP, p_NP=p_NP, x_NP=x_NP, y_A_solvation=y_A_solvation, y_C_solvation=y_C_solvation, phi_solvation=phi_solvation, p_solvation=p_solvation, x_solvation=x_solvation, y_A_compressible=y_A_compressible, y_C_compressible=y_C_compressible, phi_compressible=phi_compressible, p_compressible=p_compressible, x_compressible=x_compressible)
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment