From 53350cbeb8e543afd1cc0e2dd2ddf6dfefbaf555 Mon Sep 17 00:00:00 2001
From: JanHab <Jan.Habscheid@web.de>
Date: Fri, 16 Aug 2024 16:29:19 +0200
Subject: [PATCH] Added 4eq system for oneD + refined mesh for oneD

---
 .vscode/settings.json                         |   6 +-
 environment.yml                               |   2 +-
 src/Eq04.py                                   | 311 ++++++++++++++++++
 .../oneD_refined_mesh.cpython-312.pyc         | Bin 0 -> 2838 bytes
 src/oneD_refined_mesh.py                      |  49 +++
 5 files changed, 365 insertions(+), 3 deletions(-)
 create mode 100644 src/Eq04.py
 create mode 100644 src/__pycache__/oneD_refined_mesh.cpython-312.pyc
 create mode 100644 src/oneD_refined_mesh.py

diff --git a/.vscode/settings.json b/.vscode/settings.json
index cbbbea7..7f757a8 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -1,5 +1,7 @@
 {
     "cSpell.enableFiletypes": [
-        "!markdown"
-    ]
+        "!markdown",
+        "!python"
+    ],
+    "python.analysis.typeCheckingMode": "off"
 }
\ No newline at end of file
diff --git a/environment.yml b/environment.yml
index ca883de..46c7db4 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,4 +1,4 @@
-name: fenicsx
+name: fenicsx-env
 channels:
   - conda-forge
   - defaults
diff --git a/src/Eq04.py b/src/Eq04.py
new file mode 100644
index 0000000..bed5462
--- /dev/null
+++ b/src/Eq04.py
@@ -0,0 +1,311 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+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, Cell
+from basix.ufl import element, mixed_element
+from petsc4py import PETSc
+import matplotlib.pyplot as plt
+
+from oneD_refined_mesh import create_refined_mesh
+
+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_α)=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.
+
+    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 ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi
+        
+        def J_C(y_A, y_C, phi, p):
+            return ln(y_C) + a2 * (p - 1) * (solvation + 1) + 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(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
+        )
+        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 = 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_4eq(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/src/__pycache__/oneD_refined_mesh.cpython-312.pyc b/src/__pycache__/oneD_refined_mesh.cpython-312.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..8bd99bacda254a65624f8611c0fba32b9248ffc7
GIT binary patch
literal 2838
zcmbVOeP|Td7N42@-t32KYA_pH#}At3CAw&|UnxdSjL{mCn6%3F$uOCF_amL1+1#0}
zMrP|icz@VZsD*+_d7y&eiuofCgr<LW-#?E)TgYT>y9|-{N1-4To~?z)llR9vcQzAO
zu^;xr?78ROd(QdYbI-?q5e^3#82=GJPaX6z%y*P=hIe&xe-I{@7{nwQ#3D{+lN?Z&
z?8>>5ZkDdOWlzqV^xAWe?92I+ewcf40Qqnb`LRFcMS;J%lTFiHNAO%%@TjQpgW`lH
zWpEZzZs!`epPbP%J46w-l${6%?@|!QJ6Ov#G;#!pJBPIlr9(iibp)rf3eI6g7c_lZ
z##Trv<R&l?BrMCC<t11z5JlovuR=)xJ{NrVm*MvkW55O8et}`gR-L8WbWAp$uP-r(
zLoNfZeGd4WR2~L!;6CI&*TgVI?jozWK%2+bR{U*kx~Hy%VL-3<9B?;R5!9_9-<<!Z
zyMCt85(?0>MHdP@r43u|Q5#RXNgWHlXVKkY!3hTkG*H83xKkYB=7Kj@JyztDzkm;#
zY#$g6(!b!*KkS3hst-a>Y2)|$AWZfDQ6K0u#_VyPbNpTO7@qno7QGFa?eud_*J;CU
z=+Pa@8D2-@ub+cx_15Fad7IQjEI#Ba`cb4902<Vzjuk!T%uv*TJ+y8v=0u=lwb+zq
ziXnsRXFd!|T$+)%)6YPpgpJT-vm^bG8sRkSNPfh9(Bg*k#W2h08=NzfEzayA)i*lx
z$+$CnNDZG6GMXG*#fT9Z1FfvB`3%#sUOWqF6y#y_5iIIh<3(Oo@D7xPI-zA%MU;73
zMffvWJp+678iIhPVUrhiUe91&#wnejPzwqY$#l?GHzE=-hjmOeo87TCZ1(X|#J|dG
zI<dDXmO=H5t=G1!Wz>R<_zCP7-^mZ>b+8Ut3Fg~nHQnCDw`W9xgvR;Wf~~f()Pc3E
zqOVyyZdXBGP4RXiguON6^}C#~3*#dfv4v?_UDGoB1g~`gWD!mCDV6XA4O7hnhm>sJ
z*jZc5asggKv;8hDW_M{OSbhh0`l~PL_@Q?IO1K{@D6azQ^e3})eAY|rx_CCLSz%jX
zP3nI)eq0`!$oHn_HYHKdzy958`9TbmcSPCp9~09=RMNPEB~XSL(BaWCVjg!yEKbQ=
ze!!#UMG0Fjkg=L1RgoYl6!-uuR;-Yu%7vT){X>-Aw!HAhpp4MkZ27Qk_ZC>q=K(FD
z)fiORVI|O)B(xJ)!M@?4A>rikiDUg8Q9=tPfo95B5lK}E%2HdJ09c}r^1B{i*lo8j
z;T=r$tb{cy0581S$^7;^QYG#*1TvwyW_b}#6O4hXIWem|u1>H#wn594(=r;(avF#Y
z<H1{roT4YB#2Xp}bz)LfGGd}%l?t?>YKaj7PbubdStXrFXY~X$>>c0{Nhhj2t>c`E
zu&gCCA|)Wb`i0dzvi-6%KW#Oy`l)f9bODMmeA@RA9wp`;&-gda>|OTpEAdU2J7+tu
zh~-%M*h1{a@oKzhK{VsNm0f0Q|ND_zpmip>9B5x@d-ihAY|oXca%`@<Y%KJdiQa1G
z{z}a3Jh*u1cH4=iv2n9)eCf@TX8hz*^0XO0eIZ(lKYMNc^|q^R<<augt>*bZSGJjZ
z4^?*!RHn>bgVilZ7PpvNUSAwF`QuBY<0e0D#>X#2S6bUHKR5eadGkVI>4l?9qbF{+
zj=g6C5taj8PYFaY;{pY;`P%O5Jy(0mQ}eM~-SbAJ&wS}{wR><eW_BN~ZhL+4fVpjG
zamw5}vNU$e+<MB4pL#No3z3!PmW%IJo1ZTotF7BKGi9#hOGj3su@BF`cm5i8J#aNp
z?x^y6ZcTie{d=~;R$f^g{YOvr#o?vRBh~0=X|UE5D)rYwEu}-XaHKS_5{`Y`KJ%w5
zxI8+KZm;jDhF>cU)LJ%{hHJa`f8#z--cS<0^|YCuwwX6(4u9?0UUeT>2}ZsR#?4^-
zZ{3&o&F(9|JTFvI)!<O6uhtYP9sTc59>#N!)o4QelW#q|-Q5wge1d>fNr2dd;LcAI
zXdV#AxqNmHFm^WikSfE!y0i6ca}`o@JsVw6P6_Q6b`afQ?!%NJo$&L({~82iLB{*Y
zUSMhA(RRbnJr~QeKXNYCeLu*s&3BlVJIuyA%rp1=+u3Nj`5psgq1z^vJwH>9>;bX=
E1G3uNxBvhE

literal 0
HcmV?d00001

diff --git a/src/oneD_refined_mesh.py b/src/oneD_refined_mesh.py
new file mode 100644
index 0000000..a335934
--- /dev/null
+++ b/src/oneD_refined_mesh.py
@@ -0,0 +1,49 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh
+from ufl import Mesh
+
+
+# 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
\ No newline at end of file
-- 
GitLab