Hello researchers
I’m using VScode connected to Ubuntu WSL and have installed FEniCSx using conda.
I’m trying to implement a P2D model of a solid-state battery cathode, which is divided into macroscopic and microscopic scales, linked by the BV equations. The BV equation requires both the potential \phi_s, \ \phi_e at the macroscopic scale and the concentration c(r=r_p) at the microscopic circular boundary. Each point at the macroscopic scale corresponds to a circle at the microscopic scale, so I used 6*6 circles as representatives.
Macroscopic scale: a square with side length L
\nabla \cdot ( - \sigma_e \nabla \phi_e) = i_{BV} \frac{3}{r_p}
\phi_e = 0 \quad at \quad x=L
\nabla \cdot ( - \sigma_s \nabla \phi_s) = - i_{BV} \frac{3}{r_p}
- \sigma_s \nabla \phi_s \cdot \boldsymbol{n} = - I_{appl} \quad at \quad x=0
Microscopic scale: a circle with raidus r_p (in this code using 2D circle rather than 1D line)
\frac{\partial c}{\partial t} - D \nabla^2 c = 0
- D \nabla c \cdot \boldsymbol{n} = - \frac{i_{BV}}{F} \quad at \quad r = r_p
BV equation:
i_{BV} = i_0 \left( \exp(\frac{F \eta}{2 R T}) - \exp(-\frac{F \eta}{2 R T}) \right)
\eta = \phi_s - \phi_e - U_{OCP} (\frac{c(r=r_p)}{c_{max}})
To ensure the uniqueness of the \phi_s solution, I solved \phi_s and \phi_e together, and rewrote the BV equations as functions of \phi_s and \phi_e.
However, during the solution process, my potential doesn’t seem to evolve. The following value is always displayed.
phi_e range: [-1.4609e-02, 0.0000e+00] V
phi_s range: [4.1637e+00, 4.1637e+00] V
Is this a problem with my code? Or is there a problem with the way I ensured the uniqueness of the solution?
# Simplified code for questions
import numpy as np
from mpi4py import MPI
import gmsh
import ufl
from dolfinx import fem, io, mesh
from dolfinx.fem import (FunctionSpace, Function, form, assemble_scalar,
locate_dofs_topological, dirichletbc, Constant,
Expression)
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmsh import model_to_mesh
import dolfinx.fem.petsc
import time
import os
from scipy.interpolate import interp1d
import dolfinx.nls as nls
from dolfinx.fem.petsc import NonlinearProblem
import pyvista
from dolfinx import plot
os.environ["FENICSX_JIT_TIMEOUT"] = "300"
os.environ["FENICSX_CACHE_DIR"] = "/tmp/fenics_cache"
# ============================================================================
# 1. Parameter Settings
# ============================================================================
L0 = 1.0e-6 # Characteristic length
L_macro = 40.0e-6 / L0 # Macroscopic square side length [m]
R_micro = 2.0e-6 / L0 # Microscopic circle radius [m]
h_macro = 4.0e-6 / L0 # Macroscopic mesh size [m]
h_micro = 0.2e-6 / L0 # Microscopic mesh size [m]
N_macro = 6 # Number of points per direction (6×6=36 points)
# Material properties
sigma_e = 1.0e-2 / (L0 ** 2) # Solid electrolyte conductivity [S/m]
sigma_s = 100.0 / (L0 ** 2) # Active material conductivity [S/m]
D_s = 1.0e-15 / (L0 ** 2) # Active material diffusion coefficient [m²/s]
c_s_max = 50000.0 # Maximum lithium concentration [mol/m³]
# Physical constants
F = 96485.0 # Faraday constant [C/mol]
R = 8.314 # Gas constant [J/(mol·K)]
T = 298.0 # Temperature [K]
alpha = 0.5 # Charge transfer coefficient
# Solver parameters
max_iter = 1000 # Maximum iterations per time step
tolerance = 1e-6 # Convergence tolerance
# Variable time steps
time_intervals = [
(0.0, 0.1, 0.01), # 0:0.01:0.1
(0.2, 1.0, 0.1), # 0.2:0.1:1
(1.0, 5.0, 1.0) # 1:1:5
]
# Boundary conditions
I_app = -10.0 / L0 # Applied current density [A/m²] (negative for discharge)
# ============================================================================
# 2. Generate Time Step Sequence
# ============================================================================
def generate_time_steps():
time_steps = []
for start, end, step in time_intervals:
num_steps = int((end - start) / step) + 1
interval_steps = np.linspace(start, end, num_steps)
if not time_steps:
time_steps.extend(interval_steps)
else:
time_steps.extend(interval_steps[1:])
return np.array(time_steps)
# ============================================================================
# 3. Create Meshes
# ============================================================================
def create_macro_mesh():
print("Creating macroscopic geometry...")
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
# Create square
gmsh.model.add("macro_domain")
p1 = gmsh.model.geo.addPoint(0, 0, 0, h_macro)
p2 = gmsh.model.geo.addPoint(L_macro, 0, 0, h_macro)
p3 = gmsh.model.geo.addPoint(L_macro, L_macro, 0, h_macro)
p4 = gmsh.model.geo.addPoint(0, L_macro, 0, h_macro)
l1 = gmsh.model.geo.addLine(p1, p2) # Bottom boundary
l2 = gmsh.model.geo.addLine(p2, p3) # Right boundary
l3 = gmsh.model.geo.addLine(p3, p4) # Top boundary
l4 = gmsh.model.geo.addLine(p4, p1) # Left boundary
square_loop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
square_surface = gmsh.model.geo.addPlaneSurface([square_loop])
gmsh.model.geo.synchronize()
# Mark boundaries
gmsh.model.addPhysicalGroup(2, [square_surface], 1, name="macro_domain")
gmsh.model.addPhysicalGroup(1, [l1], 1, name="bottom")
gmsh.model.addPhysicalGroup(1, [l2], 2, name="right")
gmsh.model.addPhysicalGroup(1, [l3], 3, name="top")
gmsh.model.addPhysicalGroup(1, [l4], 4, name="left")
# Generate mesh
gmsh.model.mesh.generate(2)
# Convert to dolfinx mesh
mesh_comm = MPI.COMM_WORLD
model_rank = 0
result = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
macro_mesh, cell_tags, facet_tags = result.mesh, result.cell_tags, result.facet_tags
gmsh.finalize()
print(f"Macro mesh: {macro_mesh.topology.index_map(2).size_local} elements")
return macro_mesh, facet_tags
def create_micro_mesh():
print("Creating microscopic geometry...")
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
# Create circle
gmsh.model.add("micro_domain")
center = gmsh.model.geo.addPoint(0, 0, 0, h_micro)
p1 = gmsh.model.geo.addPoint(R_micro, 0, 0, h_micro)
p2 = gmsh.model.geo.addPoint(0, R_micro, 0, h_micro)
p3 = gmsh.model.geo.addPoint(-R_micro, 0, 0, h_micro)
p4 = gmsh.model.geo.addPoint(0, -R_micro, 0, h_micro)
# Create arcs
arc1 = gmsh.model.geo.addCircleArc(p1, center, p2)
arc2 = gmsh.model.geo.addCircleArc(p2, center, p3)
arc3 = gmsh.model.geo.addCircleArc(p3, center, p4)
arc4 = gmsh.model.geo.addCircleArc(p4, center, p1)
circle_loop = gmsh.model.geo.addCurveLoop([arc1, arc2, arc3, arc4])
circle_surface = gmsh.model.geo.addPlaneSurface([circle_loop])
gmsh.model.geo.synchronize()
# Mark boundaries
gmsh.model.addPhysicalGroup(2, [circle_surface], 1, name="micro_domain")
gmsh.model.addPhysicalGroup(1, [arc1, arc2, arc3, arc4], 1, name="boundary")
# Generate mesh
gmsh.model.mesh.generate(2)
# Convert to dolfinx mesh
mesh_comm = MPI.COMM_WORLD
model_rank = 0
result = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
micro_mesh, micro_cell_tags, micro_facet_tags = result.mesh, result.cell_tags, result.facet_tags
gmsh.finalize()
print(f"Micro mesh: {micro_mesh.topology.index_map(2).size_local} elements")
return micro_mesh, micro_facet_tags
# ============================================================================
# 4. Load OCP Data
# ============================================================================
def load_ocp_data(filename="OCP.txt"):
data = np.loadtxt(filename)
soc_data = data[:, 0] # SOC (0-1)
ocp_data = data[:, 1] # Voltage
# Create interpolation function
ocp_interp = interp1d(soc_data, ocp_data, kind='linear', fill_value='extrapolate')
return ocp_interp
# ============================================================================
# 5. Physical Model Functions
# ============================================================================
def BV_current_density(phi_s, phi_e, c_s_surf, ocp_func):
soc = c_s_surf / c_s_max
U_ocp = float(ocp_func(soc))
eta = phi_s - phi_e - U_ocp
i0 = 1.0 / L0
if abs(eta) < 1e-12:
j = 0.0
else:
arg = alpha * F * eta / (R * T)
arg = np.clip(arg, -100.0, 100.0)
j = i0 * np.sinh(arg)
return j
def solve_micro_diffusion(c_s_old, j, dt, V_micro, dx_micro, ds_micro):
u = ufl.TrialFunction(V_micro)
v = ufl.TestFunction(V_micro)
# Weak form (implicit Euler)
a = u * v * dx_micro + dt * D_s * ufl.dot(ufl.grad(u), ufl.grad(v)) * dx_micro
L = c_s_old * v * dx_micro - dt * (j / F) * v * ds_micro
c_s_new = Function(V_micro)
problem = fem.petsc.LinearProblem(
a, L, u=c_s_new,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
petsc_options_prefix="c_s_"
)
problem.solve()
return c_s_new
def get_surface_concentration(c_s_func, R_micro):
V_micro = c_s_func.function_space
dof_coords = V_micro.tabulate_dof_coordinates()
boundary_indices = []
for i, coord in enumerate(dof_coords):
distance = np.linalg.norm(coord[:2])
if abs(distance - R_micro) < 1e-6 / L0:
boundary_indices.append(i)
if boundary_indices:
values = c_s_func.x.array[boundary_indices]
return np.mean(values)
else:
return 0.0
# ============================================================================
# 6. Spatial Interpolator
# ============================================================================
def create_spatial_interpolator(macro_mesh, N_macro):
# Create interpolation grid
x_coords = np.linspace(0, L_macro, N_macro)
y_coords = np.linspace(0, L_macro, N_macro)
x_grid, y_grid = np.meshgrid(x_coords, y_coords, indexing='ij')
points = np.column_stack([x_grid.ravel(), y_grid.ravel()])
# Create scalar function space for interpolation
V_scalar = fem.functionspace(macro_mesh, ("Lagrange", 1))
return V_scalar, points, x_coords, y_coords
def discrete_to_spatial(discrete_values, V_scalar, x_coords, y_coords, N_macro):
spatial_func = Function(V_scalar)
dof_coords = V_scalar.tabulate_dof_coordinates()
if discrete_values.ndim == 1:
values_grid = discrete_values.reshape(N_macro, N_macro)
else:
values_grid = discrete_values
for i in range(dof_coords.shape[0]):
x, y = dof_coords[i, :2]
x_idx = np.argmin(np.abs(x_coords - x))
y_idx = np.argmin(np.abs(y_coords - y))
value = values_grid[x_idx, y_idx]
spatial_func.x.array[i] = float(value)
return spatial_func
# ============================================================================
# 7. Main Solver Logic
# ============================================================================
# Create results directory
if not os.path.exists("results"):
os.makedirs("results")
# Create geometries
print("=" * 60)
print("P2D Model for All-Solid-State Battery Cathode")
print("=" * 60)
macro_mesh, facet_tags = create_macro_mesh()
micro_mesh, micro_facet_tags = create_micro_mesh()
# Load OCP data
ocp_func = load_ocp_data()
# Generate time step sequence
time_steps = generate_time_steps()
n_steps = len(time_steps) - 1
print(f"Time step sequence: {time_steps[:10]}... ({len(time_steps)} time points)")
print(f"Number of time steps: {n_steps}")
# Set up function spaces
V_macro_coupled = fem.functionspace(macro_mesh, ("Lagrange", 1, (2,)))
V_macro_scalar = fem.functionspace(macro_mesh, ("Lagrange", 1))
V_micro = fem.functionspace(micro_mesh, ("Lagrange", 1))
# Define measures
dx_macro = ufl.Measure("dx", domain=macro_mesh)
ds_macro = ufl.Measure("ds", domain=macro_mesh, subdomain_data=facet_tags)
dx_micro = ufl.Measure("dx", domain=micro_mesh)
ds_micro = ufl.Measure("ds", domain=micro_mesh, subdomain_data=micro_facet_tags)
# Initialize functions
phi_coupled = Function(V_macro_coupled, name="phi_coupled")
phi_coupled_old = Function(V_macro_coupled, name="phi_coupled_old")
initial_phi_s = float(ocp_func(1.0e-3) + 0.0 * 1.0e-4)
phi_coupled_old.x.array[0::2] = 0.0
phi_coupled_old.x.array[1::2] = initial_phi_s
phi_coupled.x.array[:] = phi_coupled_old.x.array
# Create spatial interpolator
V_scalar_interp, discrete_points, x_coords, y_coords = create_spatial_interpolator(macro_mesh, N_macro)
# Initialize microscopic functions
micro_functions = []
c_s_surf_discrete = np.zeros(N_macro * N_macro)
initial_c_s = 1.0e-3 * c_s_max
for i in range(N_macro * N_macro):
c_s = Function(V_micro, name=f"c_s_{i}")
c_s.x.array[:] = initial_c_s
micro_functions.append(c_s)
c_s_surf_discrete[i] = initial_c_s
# Initial surface concentration function
c_surf_func = discrete_to_spatial(c_s_surf_discrete, V_scalar_interp, x_coords, y_coords, N_macro)
# Initial OCP function
ocp_discrete = np.zeros(N_macro * N_macro)
for i in range(N_macro * N_macro):
soc = c_s_surf_discrete[i] / c_s_max
ocp_discrete[i] = float(ocp_func(soc))
U_ocp_func = discrete_to_spatial(ocp_discrete, V_scalar_interp, x_coords, y_coords, N_macro)
# Set boundary conditions
def right_boundary(x):
return np.isclose(x[0], L_macro)
fdim = macro_mesh.topology.dim - 1
right_facets = locate_entities_boundary(macro_mesh, fdim, right_boundary)
V_phi_e = V_macro_coupled.sub(0)
dofs_right_phi_e = locate_dofs_topological(V_phi_e, fdim, right_facets)
bc_phi_e = dirichletbc(Constant(macro_mesh, 0.0), dofs_right_phi_e, V_phi_e)
bcs = [bc_phi_e]
# ============================================================================
# 8. Time Loop
# ============================================================================
print("\nStarting time loop...")
print(f"Total time steps: {n_steps}")
print(f"Total simulation time: {time_steps[-1]:.2f}s")
start_time = time.time()
convergence_stats = []
for step in range(n_steps):
t_current = time_steps[step]
t_next = time_steps[step + 1]
dt = t_next - t_current
if step % 10 == 0:
print(f"\nTime step {step}/{n_steps}, t = {t_current:.3f}s, dt = {dt:.3f}s")
# Step 1: Solve macroscopic potential equations
print(" Step 1: Solving macroscopic potential equations...")
# Update OCP function
for i in range(N_macro * N_macro):
soc = c_s_surf_discrete[i] / c_s_max
ocp_discrete[i] = float(ocp_func(soc))
U_ocp_func = discrete_to_spatial(ocp_discrete, V_scalar_interp, x_coords, y_coords, N_macro)
phi = Function(V_macro_coupled)
phi.x.array[:] = phi_coupled.x.array
v = ufl.TestFunction(V_macro_coupled)
phi_e, phi_s = ufl.split(phi)
v_e, v_s = ufl.split(v)
U_ocp = U_ocp_func
i0 = Constant(macro_mesh, 1.0 / L0)
eta = phi_s - phi_e - U_ocp
arg = alpha * F * eta / (R * T)
j_BV = i0 * ufl.sinh(arg)
F_eq = (sigma_e * ufl.dot(ufl.grad(phi_e), ufl.grad(v_e)) +
sigma_s * ufl.dot(ufl.grad(phi_s), ufl.grad(v_s))) * dx_macro
F_eq += -j_BV * (3 / R_micro) * (v_e - v_s) * dx_macro
F_eq += -I_app * v_s * ds_macro(4)
petsc_options = {
"snes_type": "newtonls",
"snes_linesearch_type": "basic",
"snes_max_it": 200,
"snes_rtol": 1e-5,
"snes_atol": 1e-5,
"ksp_type": "gmres",
"ksp_max_it": 100,
"ksp_rtol": 1e-6,
"pc_type": "ilu",
}
# Create nonlinear problem
macro_problem = NonlinearProblem(F_eq, phi,
petsc_options_prefix="phi_bv_",
bcs=bcs,petsc_options=petsc_options)
try:
macro_problem.solve()
phi_coupled.x.array[:] = phi.x.array
converged = True
except Exception as e:
print(f"Nonlinear solver failed: {e}")
converged = False
if not converged:
print(" Warning: Macroscopic potential solution did not fully converge, continuing...")
# Step 2: Solve microscopic diffusion equations
print(" Step 2: Solving microscopic diffusion equations...")
# Extract potentials from coupled solution
phi_e_func = Function(V_macro_scalar)
phi_s_func = Function(V_macro_scalar)
phi_e_func.x.array[:] = phi_coupled.x.array[0::2]
phi_s_func.x.array[:] = phi_coupled.x.array[1::2]
# Solve at discrete points
for i, point in enumerate(discrete_points):
point_np = np.array(point, dtype=np.float64).reshape(1, -1)
phi_e_expr = Expression(phi_e_func, point_np)
phi_s_expr = Expression(phi_s_func, point_np)
U_ocp_expr = Expression(U_ocp_func, point_np)
phi_e_val = phi_e_expr.eval(macro_mesh, point_np)[0,0]
phi_s_val = phi_s_expr.eval(macro_mesh, point_np)[0,0]
U_ocp_val = U_ocp_expr.eval(macro_mesh, point_np)[0,0]
c_surf_val = c_s_surf_discrete[i]
eta = phi_s_val - phi_e_val - U_ocp_val
i0_val = 1.0 / L0
if abs(eta) < 1e-12:
j_val = 0.0
else:
arg = alpha * F * eta / (R * T)
arg = np.clip(arg, -100.0, 100.0)
j_val = i0_val * np.sinh(arg)
# Solve microscopic diffusion
c_s_new = solve_micro_diffusion(
micro_functions[i], j_val, dt,
V_micro, dx_micro, ds_micro
)
micro_functions[i].x.array[:] = c_s_new.x.array
c_s_surf_discrete[i] = get_surface_concentration(c_s_new, R_micro)
# Update surface concentration function
c_surf_func = discrete_to_spatial(c_s_surf_discrete, V_scalar_interp, x_coords, y_coords, N_macro)
# Save old solution
phi_coupled_old.x.array[:] = phi_coupled.x.array
# Debug information
if step % 2 == 0:
print(f"\nStep {step} debug information (t={t_current:.3f}s, dt={dt:.3f}s):")
print("-" * 40)
phi_e_array = phi_coupled.x.array[0::2]
phi_s_array = phi_coupled.x.array[1::2]
print(f"phi_e range: [{np.min(phi_e_array):.4e}, {np.max(phi_e_array):.4e}] V")
print(f"phi_s range: [{np.min(phi_s_array):.4e}, {np.max(phi_s_array):.4e}] V")
c_surf_array = c_surf_func.x.array
print(f"Surface concentration range: [{np.min(c_surf_array):.2f}, {np.max(c_surf_array):.2f}] mol/m³")
print(f"SOC range: [{np.min(c_surf_array)/c_s_max:.6f}, {np.max(c_surf_array)/c_s_max:.6f}]")
U_ocp_array = U_ocp_func.x.array
eta_array = phi_s_array - phi_e_array - U_ocp_array
print(f"Overpotential (eta) range: [{np.min(eta_array):.4e}, {np.max(eta_array):.4e}] V")
print(f"Average |eta|: {np.mean(np.abs(eta_array)):.4e} V")
print(f"Maximum |eta|: {np.max(np.abs(eta_array)):.4e} V")
print(f"Applied current density: {I_app} A/m²")
print("-" * 40)
convergence_stats.append(converged)
end_time = time.time()
total_time = end_time - start_time
successful_steps = sum(convergence_stats)
print(f"\nCalculation completed!")
print(f"Total computation time: {total_time:.2f} s")
print(f"Average time per step: {total_time/n_steps:.3f} s")
print(f"Successful steps: {successful_steps}/{n_steps} ({successful_steps/n_steps*100:.1f}%)")
# ============================================================================
# 9. Results Output
# ============================================================================
print("\nCalculation results statistics:")
print("=" * 50)
c_surf_array = c_surf_func.x.array
avg_surf_conc = np.mean(c_surf_array)
print(f"Average surface concentration: {avg_surf_conc:.2f} mol/m³")
print(f"Surface concentration range: {np.min(c_surf_array):.2f} - {np.max(c_surf_array):.2f} mol/m³")
print(f"Average SOC: {avg_surf_conc/c_s_max:.6f}")
phi_e_min = np.min(phi_coupled.x.array[0::2])
phi_e_max = np.max(phi_coupled.x.array[0::2])
phi_s_min = np.min(phi_coupled.x.array[1::2])
phi_s_max = np.max(phi_coupled.x.array[1::2])
print(f"Solid electrolyte potential range: {phi_e_min:.4e} - {phi_e_max:.4e} V")
print(f"Active material potential range: {phi_s_min:.4e} - {phi_s_max:.4e} V")
print(f"Potential difference range: {(phi_s_min-phi_e_min):.4e} - {(phi_s_max-phi_e_max):.4e} V")
print("\nSaving final results...")
final_time = time_steps[-1]
phi_e = Function(V_macro_scalar)
phi_s = Function(V_macro_scalar)
phi_e.x.array[:] = phi_coupled.x.array[0::2]
phi_s.x.array[:] = phi_coupled.x.array[1::2]
with io.XDMFFile(macro_mesh.comm,
f"results/final_results_t{final_time:.1f}s.xdmf", "w",
encoding=io.XDMFFile.Encoding.ASCII) as xdmf:
xdmf.write_mesh(macro_mesh)
xdmf.write_function(phi_e)
xdmf.write_function(phi_s)
xdmf.write_function(c_surf_func)
xdmf.write_function(U_ocp_func)
print(f"Final results saved to results/final_results_t{final_time:.1f}s.xdmf")
print("=" * 50)
print("\nProgram execution completed!")
OCP data is
| 0 | 4.184617547806524 |
|---|---|
| 0.02052834550915783 | 4.149555461411479 |
| 0.04262486142300803 | 4.123168889603085 |
| 0.06472137733685823 | 4.0982688852205165 |
| 0.08681789325070821 | 4.0749483587278865 |
| 0.10891440916455841 | 4.052278205484055 |
| 0.1310109250784086 | 4.030351335953136 |
| 0.1531074409922586 | 4.008981929206901 |
| 0.1752039569061088 | 3.988355806173579 |
| 0.19730047281995877 | 3.967636772676143 |
| 0.20834873077688387 | 3.95750953208771 |
| 0.23044524669073407 | 3.937255050910844 |
| 0.2525417626045843 | 3.9169076592698637 |
| 0.27463827851843425 | 3.8969319094853403 |
| 0.29673479443228445 | 3.877606532949615 |
| 0.31883131034613466 | 3.8590244401268023 |
| 0.34092782625998463 | 3.8415572728733585 |
| 0.36302434217383484 | 3.825019210261055 |
| 0.38512085808768504 | 3.8100606255386906 |
| 0.40721737400153524 | 3.797331891955064 |
| 0.429313889915385 | 3.7853464420843497 |
| 0.4514104058292352 | 3.7748475596394604 |
| 0.4845551797000103 | 3.7612826318788075 |
| 0.5176999535707858 | 3.749297182008093 |
| 0.5508447274415609 | 3.738612478634976 |
| 0.6060860172261862 | 3.7218885950944443 |
| 0.6502790490538866 | 3.7086953091902473 |
| 0.7165685967954367 | 3.6874188129081267 |
| 0.771809886580062 | 3.6653990329130934 |
| 0.8049546604508371 | 3.6496042540137026 |
| 0.8491476922785375 | 3.62210275663594 |
| 0.8822924661493126 | 3.593114691832352 |
| 0.924978917346523 | 3.520818736943596 |
| 0.9325118204989722 | 3.4987176602924634 |
| 0.9445644655428902 | 3.4508873533665434 |
| 0.9531017557823325 | 3.404919385076865 |
| 0.959630271847788 | 3.356861963683111 |
| 0.9626434331087679 | 3.333832556644705 |
| 0.968526271761156 | 3.2820504579784666 |
| 0.9711807233482097 | 3.252923027478708 |
| 0.973189497522196 | 3.2306941989394184 |
| 0.9752987104048816 | 3.2039685039370074 |
| 0.9772070458701689 | 3.174227864374096 |
| 0.9792158200441552 | 3.150466013176924 |
| 0.9809232780920434 | 3.1250689378113448 |
| 0.9821285425964357 | 3.0989053511168247 |
| 0.9842377554791213 | 3.066575606620601 |
| 0.9863469683618069 | 3.0367668327173387 |
| 0.9872509167401007 | 3.0094279286517756 |
| 0.9881548651183949 | 2.984439659328298 |
| 0.9900488521967246 | 2.9535456027180276 |
| 0.9912684650880736 | 2.925878193797204 |
| 0.9922728521750668 | 2.9033938614816 |
| 0.9932772392620599 | 2.8829535593765065 |
| 0.994281626349053 | 2.8574031817451386 |
| 0.9949512177403816 | 2.8270834002892493 |
| 0.9962904005230393 | 2.7931865659649686 |
| 0.99808394889267 | 2.7554206744565093 |
| 0.9989687660883546 | 2.7173871123252447 |
| 1 | 2.7005238630885424 |