From 8633936cdd60b3388a3ce81335540c1b7fee53bc Mon Sep 17 00:00:00 2001
From: JanHab <Jan.Habscheid@web.de>
Date: Sat, 17 Aug 2024 16:36:35 +0200
Subject: [PATCH] Added PoissonBoltzmann + reduced system of equations with
y(phi,p)
---
examples/Compressibility.py | 132 +++++++
.../ValidationAnalyticalMethod.py | 22 ++
examples/PoissonBoltzmann.py | 4 +-
src/Eq02.py | 355 ++++++++++++++++++
4 files changed, 511 insertions(+), 2 deletions(-)
create mode 100644 examples/Compressibility.py
create mode 100644 examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py
create mode 100644 src/Eq02.py
diff --git a/examples/Compressibility.py b/examples/Compressibility.py
new file mode 100644
index 0000000..d6b8176
--- /dev/null
+++ b/examples/Compressibility.py
@@ -0,0 +1,132 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Further imports
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Define the parameters and boundary conditions
+phi_left = 6.0
+phi_right = 0.0
+p_right = 0.0
+y_A_L, y_C_L = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+number_cells = 1024*8
+K_vec = ['incompressible', 15_000, 5_000, 1_500, 500]
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+refinement_style = 'hard_hard_log'
+rtol = 1e-7
+
+relax_param = 0.08
+max_iter = 10_000
+
+# Calculate the total number density, based on the pressure
+def n_expr(p, K):
+ if K == 'incompressible':
+ return np.ones_like(p)
+ return (p-1) / K + 1
+
+# Solve the system
+y_A, y_C, y_S, phi, p, n, x = [], [], [], [], [], [], []
+for K in K_vec:
+ y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_L, y_C_L, K, Lambda2, a2, number_cells, relax_param = relax_param, max_iter=max_iter, refinement_style=refinement_style, return_type='Vector', rtol=rtol)
+ y_A.append(y_A_)
+ y_C.append(y_C_)
+ y_S.append(1 - y_A_ - y_C_)
+ phi.append(phi_)
+ p.append(p_)
+ n.append(n_expr(p_, K))
+ x.append(x_)
+
+
+# Plot the results
+fig, axs = plt.subplots(2, 2, figsize=(30, 20))
+# Define plotting parameter
+labelsize = 30
+lw = 4
+legend_width = 8
+xlim = 0.05
+markers = ['-.', '--', '-', ':']
+colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
+
+# Legend for compressibility scaling
+axs[0,0].plot(0, 0, label='Incompressible', color=colors[0])
+[axs[0,0].plot(0, 0, label=f'$\kappa$ = {K_vec[i]}', color=colors[i]) for i in range(1, len(K_vec))]
+
+# Electric potential
+[axs[0,0].plot(x[i], phi[i], lw=lw, color=colors[i]) for i in range(len(K_vec))]
+axs[0,0].set_xlim(0,xlim)
+axs[0,0].grid()
+axs[0,0].set_xlabel('x [-]', fontsize=labelsize)
+axs[0,0].set_ylabel('$\\varphi$ [-]', fontsize=labelsize)
+axs[0,0].tick_params(axis='both', labelsize=labelsize)
+
+# Concentrations
+for i in range(len(K_vec)):
+ clr = colors[i]
+ axs[0,1].plot(x[i], y_A[i], markers[0], color=clr, lw=lw)
+ axs[0,1].plot(x[i], y_C[i], markers[1], color=clr, lw=lw)
+axs[0,1].plot(0, 0.1, color='grey', linestyle='--', label='Anions')
+axs[0,1].plot(0, 0.1, color='grey', linestyle=':', label='Cations')
+axs[0,1].set_xlim(0,xlim)
+axs[0,1].grid()
+axs[0,1].set_xlabel('x [-]', fontsize=labelsize)
+axs[0,1].set_ylabel('$y_\\alpha$ [-]', fontsize=labelsize)
+axs[0,1].tick_params(axis='both', labelsize=labelsize)
+
+a = plt.axes([.75, .77, .2, .2])
+for i in range(len(K_vec)):
+ clr = colors[i]
+ a.plot(x[i], y_S[i], markers[0], color=clr, lw=lw)
+axs[0,1].plot(0, 0.1, color='grey', linestyle='--', label='Solvent')
+a.set_xlim(0,xlim)
+a.grid()
+a.set_xlabel('$\delta \\varphi$ [-]', fontsize=labelsize)
+a.set_ylabel('$y_S$ [-]', fontsize=labelsize)
+a.tick_params(axis='both', labelsize=labelsize)
+
+# Number densities
+for i in range(len(K_vec)):
+ clr = colors[i]
+ axs[1,1].plot(x[i], y_A[i] * n[i], markers[0], color=clr, lw=lw)
+ axs[1,1].plot(x[i], y_C[i] * n[i], markers[1], color=clr, lw=lw)
+axs[1,1].set_xscale('log')
+axs[1,1].set_yscale('log')
+axs[1,1].set_xlim(0,xlim)
+axs[1,1].grid()
+axs[1,1].set_xlabel('log(x) [-]', fontsize=labelsize)
+axs[1,1].set_ylabel('log($n_\\alpha$) [-]', fontsize=labelsize)
+axs[1,1].tick_params(axis='both', labelsize=labelsize)
+
+# Pressure
+[axs[1,0].plot(x[i], p[i], lw=lw, color=colors[i]) for i in range(len(K_vec))]
+axs[1,0].set_yscale('log')
+axs[1,0].set_xlim(0,xlim)
+axs[1,0].set_ylim(1e-9, np.max(p))
+axs[1,0].grid()
+axs[1,0].set_xlabel('x [-]', fontsize=labelsize)
+axs[1,0].set_ylabel('log($p$) [-]', fontsize=labelsize)
+axs[1,0].tick_params(axis='both', labelsize=labelsize)
+
+
+lgnd = fig.legend(bbox_to_anchor=(0.98, 1.05), fontsize=labelsize, ncol=7, markerscale=60)
+for line in lgnd.get_lines():
+ line.set_linewidth(legend_width)
+fig.tight_layout()
+fig.show()
\ No newline at end of file
diff --git a/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py b/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py
new file mode 100644
index 0000000..d0923e4
--- /dev/null
+++ b/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py
@@ -0,0 +1,22 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+
+# Define Parameter and buondary conditions
diff --git a/examples/PoissonBoltzmann.py b/examples/PoissonBoltzmann.py
index c3a428a..a887676 100644
--- a/examples/PoissonBoltzmann.py
+++ b/examples/PoissonBoltzmann.py
@@ -168,7 +168,7 @@ plt.plot(phi_left_vec, y_C_error_L2, 'o-', label='$y_C$')
plt.plot(phi_left_vec, y_S_error_L2, 'o-', label='$y_S$')
plt.yscale('log')
plt.legend()
-plt.xlabel('$\Delta \\varphi$')
+plt.xlabel('$\delta \\varphi$')
plt.ylabel('log($L_2$)')
plt.grid()
plt.tight_layout()
@@ -191,7 +191,7 @@ plt.plot(phi_left_vec, y_C_error_inf, 'o-', label='$y_C$')
plt.plot(phi_left_vec, y_S_error_inf, 'o-', label='$y_S$')
plt.yscale('log')
plt.legend()
-plt.xlabel('$\Delta \\varphi$')
+plt.xlabel('$\delta \\varphi$')
plt.ylabel('log($L_\infty$)')
plt.grid()
plt.tight_layout()
diff --git a/src/Eq02.py b/src/Eq02.py
new file mode 100644
index 0000000..45901a1
--- /dev/null
+++ b/src/Eq02.py
@@ -0,0 +1,355 @@
+'''
+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, 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_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
+ '''
+ 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:
+ n^F = z_A y_A(φ, p) + z_C y_C(φ ,p)
+ if compressible:
+ y_alpha = C_alpha * (K+p−1)â½âˆ’κ+1)a²K exp(−z_α φ)
+ if 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
+ 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_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()
--
GitLab