diff --git a/.gitignore b/.gitignore
index c249ef59dcd0e06a9fc899cc27473529d465e55f..40ab7b5596ab2be7842145809f9f9ff8b1992a19 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,4 +10,9 @@
# Ignore pip folder
*egg-info/*
# Ignore venv files
-venv
\ No newline at end of file
+venv
+# Ignore .eggs
+*.eggs
+FENICSxDGM.egg-info
+# Ignore build files
+build
\ No newline at end of file
diff --git a/build/lib/FENICSxDGM/ElectrolyticDiode.py b/build/lib/FENICSxDGM/ElectrolyticDiode.py
deleted file mode 100644
index 836b48c19f548142f22d74f01577dde207448d40..0000000000000000000000000000000000000000
--- a/build/lib/FENICSxDGM/ElectrolyticDiode.py
+++ /dev/null
@@ -1,401 +0,0 @@
-'''
-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
diff --git a/build/lib/FENICSxDGM/Eq02.py b/build/lib/FENICSxDGM/Eq02.py
deleted file mode 100644
index 0ec56e7491732d7519b83f8a79fdd3ffc8a34bcf..0000000000000000000000000000000000000000
--- a/build/lib/FENICSxDGM/Eq02.py
+++ /dev/null
@@ -1,309 +0,0 @@
-'''
-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()
diff --git a/build/lib/FENICSxDGM/Eq04.py b/build/lib/FENICSxDGM/Eq04.py
deleted file mode 100644
index 4856d848e3847f7229ac43754add17cc33cb8f61..0000000000000000000000000000000000000000
--- a/build/lib/FENICSxDGM/Eq04.py
+++ /dev/null
@@ -1,358 +0,0 @@
-'''
-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()
diff --git a/build/lib/FENICSxDGM/EqN.py b/build/lib/FENICSxDGM/EqN.py
deleted file mode 100644
index 9b48ed187e5cf869fa7338df76a185aa1c19b42a..0000000000000000000000000000000000000000
--- a/build/lib/FENICSxDGM/EqN.py
+++ /dev/null
@@ -1,324 +0,0 @@
-'''
-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
diff --git a/build/lib/FENICSxDGM/Helper_DoubleLayerCapacity.py b/build/lib/FENICSxDGM/Helper_DoubleLayerCapacity.py
deleted file mode 100644
index e47b37efa121e9c1b5ae89707ed25366f1a01fcc..0000000000000000000000000000000000000000
--- a/build/lib/FENICSxDGM/Helper_DoubleLayerCapacity.py
+++ /dev/null
@@ -1,411 +0,0 @@
-'''
-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
diff --git a/build/lib/FENICSxDGM/__init__.py b/build/lib/FENICSxDGM/__init__.py
deleted file mode 100644
index 18742b0a46d34b5481c083238d2930ed976c323e..0000000000000000000000000000000000000000
--- a/build/lib/FENICSxDGM/__init__.py
+++ /dev/null
@@ -1,5 +0,0 @@
-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
diff --git a/dist/FENICSxDGM-1.0-py3-none-any.whl b/dist/FENICSxDGM-1.0-py3-none-any.whl
deleted file mode 100644
index de6fbfd32a614c32ba7a13917281e11d70c97f32..0000000000000000000000000000000000000000
Binary files a/dist/FENICSxDGM-1.0-py3-none-any.whl and /dev/null differ
diff --git a/dist/fenicsxdgm-1.0.linux-x86_64.tar.gz b/dist/fenicsxdgm-1.0.linux-x86_64.tar.gz
deleted file mode 100644
index 8e3f717f7cefacf05580c1564c869962d5fdfbb6..0000000000000000000000000000000000000000
Binary files a/dist/fenicsxdgm-1.0.linux-x86_64.tar.gz and /dev/null differ
diff --git a/dist/fenicsxdgm-1.0.tar.gz b/dist/fenicsxdgm-1.0.tar.gz
deleted file mode 100644
index 9e2fb5d7abd0441e5741e6b79817dc81ccab1327..0000000000000000000000000000000000000000
Binary files a/dist/fenicsxdgm-1.0.tar.gz and /dev/null differ
diff --git a/setup.py b/setup.py
index 8ebf6ca98835a76dc39f383911eebc8f50f22a0a..84f340168dc9950b55c37199af54badf28696190 100644
--- a/setup.py
+++ b/setup.py
@@ -8,11 +8,8 @@ setup(
author="Jan Habscheid, Lambert Theisen, Manuel Torrilhon",
author_email="Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de",
packages=["FENICSxDGM"],
- # packages=find_packages("src"),
package_dir={"FENICSxDGM":"src"},
- # package_dir={"": "src"},
package_data={"FENICSxDGM": ["tests/*"]},
- # package_data={"": ["tests/*"]},
python_requires=">=3.12.3",
install_requires=[
"pyvista == 0.43.10",
diff --git a/src/FENICSxDGM.egg-info/PKG-INFO b/src/FENICSxDGM.egg-info/PKG-INFO
deleted file mode 100644
index 7e72ec7f367d3f0270c19fc46ddca6be7e681709..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/PKG-INFO
+++ /dev/null
@@ -1,117 +0,0 @@
-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)
-
-[](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/pipelines)
-[](https://janhab.pages.rwth-aachen.de/bsc-electrolytemodels/)
-[](https://doi.org/10.5281/zenodo.13645296)
-[](https://git.rwth-aachen.de/jan.habscheid/bsc-electrolytemodels/-/tags)
-[](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
diff --git a/src/FENICSxDGM.egg-info/SOURCES.txt b/src/FENICSxDGM.egg-info/SOURCES.txt
deleted file mode 100644
index d9be7a858f3a0ecafa61865c948a3291f417312e..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/SOURCES.txt
+++ /dev/null
@@ -1,23 +0,0 @@
-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
diff --git a/src/FENICSxDGM.egg-info/dependency_links.txt b/src/FENICSxDGM.egg-info/dependency_links.txt
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/dependency_links.txt
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/src/FENICSxDGM.egg-info/entry_points.txt b/src/FENICSxDGM.egg-info/entry_points.txt
deleted file mode 100644
index 7bcc5d3ffd0275679d5e3d0d72b5a5c1dabcd4d6..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/entry_points.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-[console_scripts]
-my-example-utility = example.example_module:main
diff --git a/src/FENICSxDGM.egg-info/requires.txt b/src/FENICSxDGM.egg-info/requires.txt
deleted file mode 100644
index ef44023b95a3d14015bfa0e082666fa7b1190c0b..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/requires.txt
+++ /dev/null
@@ -1,23 +0,0 @@
-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
diff --git a/src/FENICSxDGM.egg-info/top_level.txt b/src/FENICSxDGM.egg-info/top_level.txt
deleted file mode 100644
index a0ebb6eb232b72eaaab423f11a02573a18408250..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/top_level.txt
+++ /dev/null
@@ -1 +0,0 @@
-FENICSxDGM
diff --git a/src/FENICSxDGM.egg-info/zip-safe b/src/FENICSxDGM.egg-info/zip-safe
deleted file mode 100644
index 8b137891791fe96927ad78e64b0aad7bded08bdc..0000000000000000000000000000000000000000
--- a/src/FENICSxDGM.egg-info/zip-safe
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/test_.py b/test_.py
deleted file mode 100644
index 36bb1e0c403e0ccb623857adb520edd824d76dbf..0000000000000000000000000000000000000000
--- a/test_.py
+++ /dev/null
@@ -1,40 +0,0 @@
-'''
-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