---

Introduction to Numerical Cosmology

---

Overview

• Objective: Introduce numerical methods in cosmology.
• Agenda:
  1. Numerical Cosmology
  2. Numerical Methods
  3. Python in Numerical Cosmology
  4. Hands-On Activity

---

1) Numerical Cosmology

---

Why Do We Need Numerical Cosmology?

1.1) Complexity of Cosmological Equations

• Many problems in cosmology (and not only) are nonlinear and coupled.
• Analytical solutions are often infeasible or impossible.
• Example:
  • Many-body problem

Dynamic Range of Scales

• Cosmology spans from subatomic particles to the entire universe.
• Numerical methods allow (in principle) simulating multiple scales simultaneously.
• Time-dependent phenomena require evolving solutions.

Discussion: Think of cosmological problems that cannot be solved analytically and why.

Why Now? The Role of Modern Computing

1.2) Advancements in Computational Power

• Modern Laptops vs. Italian Supercomputers at Cineca:

  • Y-MP8/464 (1992): ~1.33 GFlop/s.
  • SP Power3 375 MHz 16-way (2000): ~100 GFlop/s.
  • Today's Laptops: ~10 GFlop/s.

• Impact on Research:

  • Complex simulations are now accessible.
  • Encourages innovation and the development of new models.

Are you curious to benchmark your PC?

1. Get the latest version of HPL: Download HPL

   • On Linux, installation is straightforward:

     ./configure --prefix=$HOME && make && make install
     

   • The xhpl executable will be located in your $HOME/bin directory. It requires a specific parameter file.

2. Create a benchmark parameter file:
   Use this website to generate one: Advanced Clustering HPL.dat Tuning Tool

3. Compare your results:
   Check your benchmark outcome against Cineca's historical data on Top500.

Visualization

How is this computational power distributed

1.3) Modern HPC design

Key idea: a supercomputer is not one very large computer. It is a large number of compute nodes connected by a fast network.

cluster → nodes → sockets / CPUs → cores → threads
           ↘ GPUs → streaming units → many lightweight threads

For cosmology codes, the important question is usually: where is the data, and which processing unit owns it?

A typical modern HPC machine

                 ┌────────────────────────────┐
                 │ Login / front-end nodes    │
                 │ edit, compile, submit jobs │
                 └──────────────┬─────────────┘
                                │
                         job scheduler
                                │
        ┌───────────────────────┴───────────────────────┐
        │                  fast network                 │
        └───────┬───────────────────┬───────────────────┘
                │                   │
       ─────────▼─────────┐ ┌───────▼──────────┐
      │ Compute node 0    │ │ Compute node 1   │   ...
      │ CPUs + RAM        │ │ CPUs + RAM       │
      │ GPUs + GPU memory │ │ GPUs + GPU memory│
      └───────────────────┘ └──────────────────┘

• CPU cores are optimized for flexible control flow and relatively complex tasks.
• GPUs are optimized for applying the same operation to very many data elements.
• The network connects nodes, but communication is much slower than local arithmetic.

Memory view: why topology matters

Shared memory inside a node

• CPU threads on the same node can usually access the same RAM.
• This is convenient, but race conditions appear if two threads update the same variable.

Distributed memory across nodes

• Each node has its own memory.
• A process cannot directly read the memory of another node.
• Data must be explicitly exchanged.

GPU memory

• GPUs often have separate high-bandwidth memory.
• Moving data between CPU memory and GPU memory can dominate the runtime if done too often.

Memory view: why topology matters

Shared memory inside a node

• CPU threads on the same node can usually access the same RAM.
• This is convenient, but race conditions appear if two threads update the same variable.

Distributed memory across nodes

• Each node has its own memory.
• A process cannot directly read the memory of another node.
• Data must be explicitly exchanged.

GPU memory

• GPUs often have separate high-bandwidth memory.
• Moving data between CPU memory and GPU memory can dominate the runtime if done too often.

OpenMP vs MPI

Feature	OpenMP	MPI
Main idea	Multiple threads inside one process	Multiple processes exchanging messages
Memory model	Shared memory	Distributed memory
Typical scale	One node, sometimes one shared-memory domain	Many nodes
Code style	Add compiler directives around loops/regions	Explicitly split data and communicate
Common use	Parallel loops, reductions, task parallelism	Domain decomposition, distributed arrays, large simulations

Practical rule: OpenMP is often easier to add to an existing loop. MPI requires more design, but it scales naturally to many nodes.

Example 1: integral with OpenMP only

Goal: approximate

$$ I = \int_a^b f(x)\,dx $$

using a parallel loop over integration bins.

function integrate_openmp(a, b, N):
    dx  = (b - a) / N
    sum = 0

    OpenMP parallel for reduction(+:sum)
    for i = 0 ... N-1:
        x_mid = a + (i + 0.5) * dx
        sum = sum + f(x_mid)

    I = dx * sum
    return I

Interpretation: all threads share the same memory, but each thread keeps a private temporary contribution to sum. The reduction combines them safely at the end.

Example 2: integral with MPI only

Goal: approximate the same integral, but distribute the bins across independent processes.

function integrate_mpi(a, b, N):
    MPI_Init()
    rank = MPI_Comm_rank()
    size = MPI_Comm_size()

    dx = (b - a) / N

    i_start, i_end = split_indices(N, rank, size)
    local_sum = 0

    for i = i_start ... i_end-1:
        x_mid = a + (i + 0.5) * dx
        local_sum = local_sum + f(x_mid)

    global_sum = MPI_Reduce(local_sum, operation=sum, root=0)
    if rank == 0:
        I = dx * global_sum
        print(I)
    MPI_Finalize()

Interpretation: each process owns only part of the interval. The final sum is obtained by explicit communication.

Why GPU porting is not automatic

A GPU is powerful only when the algorithm exposes enough parallel work and uses memory efficiently.

Main challenges:

• Data movement: CPU–GPU transfers can be more expensive than the computation.
• Memory layout: contiguous and coalesced accesses are much faster than irregular accesses.
• Branch divergence: different execution paths inside the same GPU group reduce efficiency.
• Small kernels: launching many tiny GPU tasks can waste time.
• Precision choices: double precision may be much slower or less available than single precision on some GPUs.
• Legacy code structure: global state, pointer-heavy data structures, and complex dependencies are hard to port.
• Debugging and reproducibility: parallel reductions can change floating-point ordering and therefore the last digits.

Why GPU porting is not automatic

A GPU is powerful only when the algorithm exposes enough parallel work and uses memory efficiently.

Main challenges:

• Data movement: CPU–GPU transfers can be more expensive than the computation.
• Memory layout: contiguous and coalesced accesses are much faster than irregular accesses.
• Branch divergence: different execution paths inside the same GPU group reduce efficiency.
• Small kernels: launching many tiny GPU tasks can waste time.
• Precision choices: double precision may be much slower or less available than single precision on some GPUs.
• Legacy code structure: global state, pointer-heavy data structures, and complex dependencies are hard to port.
• Debugging and reproducibility: parallel reductions can change floating-point ordering and therefore the last digits.

Extra Care Needed in Numerical Solutions

1.4) Numerical Stability and Error Analysis

• Numerical Stability:
  • Importance of choosing appropriate algorithms to prevent instability.
  • Avoid pitfalls like rounding errors and divergence.
• Error Types:
  • Truncation Errors: Due to approximations in numerical methods.
  • Round-off Errors: Due to finite precision in computations.

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

# Define the variable
x = sp.symbols('x')

# Define the function
f = sp.ln(1 + x)

# Degree of expansion (adjust as needed)
degree = 15

# Perform the Taylor expansion around x=0
taylor_expansion = sp.series(f, x, 0, degree + 1).removeO()

# Convert the symbolic expression to a numerical function for plotting
taylor_func = sp.lambdify(x, taylor_expansion, 'numpy')
original_func = sp.lambdify(x, f, 'numpy')

# Generate x values for the plot
x_vals = np.linspace(-0.9, 1.5, 500)  # Avoid x=-1 to prevent a singularity
# Correct the formatting issue by ensuring the LaTeX string is properly prepared
expansion_expr = sp.latex(taylor_expansion)

# Plot the original function and the Taylor expansion with the updated title
plt.figure(figsize=(12, 6))
plt.plot(x_vals, original_func(x_vals), label="ln(1+x)", linewidth=2, linestyle='--')
plt.plot(x_vals, taylor_func(x_vals), label=f"Taylor Expansion (degree {degree})", linewidth=2)
plt.axhline(0, color='black', linewidth=0.5, linestyle='--')
plt.axvline(0, color='black', linewidth=0.5, linestyle='--')
plt.legend()
plt.title(f"Taylor Expansion of ln(1+x): ${expansion_expr}$", fontsize=14)
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.ylim([-2.5, 2.5])
plt.tight_layout()
plt.savefig("taylor.png")
plt.show()

from scipy.interpolate import pade

# Define the variable
x = sp.symbols('x')
# Define the function
f = sp.ln(1 + x)
# Fixed order of Taylor expansion
order = 15
# Extract Taylor series coefficients using SymPy and convert to numerical values
taylor_series = f.series(x, 0, order).removeO()
coefficients = [float(taylor_series.coeff(x, i)) for i in range(order)]
# Define degrees for the Padé approximant
m, n = 3, 3  # Numerator and denominator degrees
# Compute the Padé approximant using SciPy
p, q = pade(coefficients, m)
# Define the original function
original_func = sp.lambdify(x, f, 'numpy')
# Generate x values for plotting
x_vals = np.linspace(-0.9, 2.5, 500)  # Avoid x = -1 to prevent singularity
# Evaluate the Padé approximant and the original function
pade_values = p(x_vals) / q(x_vals)
original_values = original_func(x_vals)

# Plot the original function and the Padé approximant
plt.figure(figsize=(10, 6))
plt.plot(x_vals, original_values, label="ln(1+x)", linewidth=2, linestyle='--')
plt.plot(x_vals, pade_values, label=f"Padé Approximant [{m}/{n}]", linewidth=2)
plt.axhline(0, color='black', linewidth=0.5, linestyle='--')
plt.axvline(0, color='black', linewidth=0.5, linestyle='--')
plt.legend()
plt.title(f"Padé Approximant of ln(1+x): Degrees [{m}/{n}]")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.savefig("pade.png")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Set the random seed for reproducibility
np.random.seed(1234)
# Generate an array 'x' of 1,000,000 random numbers between 0 and 1
# Convert the numbers to 32-bit floating point format to simulate single precision
x = np.random.rand(1_000_000).astype(np.float32)
# Initialize an empty list to store the computed means
sums = []
# Repeat the process 3,000 times
for i in range(3_000):
    # Randomly shuffle the elements of 'x' in place
    # This changes the order of the numbers in 'x'
    np.random.shuffle(x)
    # Compute the sum of 'x' using single-precision arithmetic
    # The 'dtype' parameter ensures that the sum is computed using 32-bit floats
    sum_i = np.sum(x, dtype=np.float32) / x.size  # Divide by the number of elements to get the mean
    # Append the computed mean to the list 'sums'
    sums.append(sum_i)
# Plot a histogram of the computed means
plt.hist(sums, bins=30, edgecolor='black')
plt.xlabel('Mean Value')
plt.ylabel('Frequency')
plt.title('Distribution of Computed Means (Single Precision)')
plt.savefig("roundoff.png")
plt.show()

Validation and Verification

• Validation:
  • Compare numerical results with known solutions.
• Verification:
  • Compare with independent implementations.

---

2) Numerical Methods in Cosmology

---

Exact vs. Approximate Numerical Methods

• Objective: Understand different numerical approaches in cosmology.
• Types:
  • Exact Numerical Solutions
  • Approximate Methods

Exact Numerical Solutions

2.1) Definition and Examples

• Definition:
  • Solutions limited only by computational precision.
• Examples:
  • Boltzmann Codes: e.g., CAMB, CLASS.
  • N-body Simulations: Simulate gravitational interactions of particles.

Considerations for Exact Methods

• Computational Resources:
  • High demand on CPU/GPU time and memory.
• Precision vs. Cost:
  • Higher resolution increases accuracy but also computational cost.

Discussion: How does increasing resolution affect simulations?

import numpy as np
import matplotlib.pyplot as plt
# Define the original function and its derivative
def f(x):
    return np.sin(x)
def f_prime(x):
    return np.cos(x)
# Set up the range and step size
x_exact = np.linspace(0, 2 * np.pi, 400)
y_exact = f(x_exact)
# Euler's method parameters
h = np.pi / 10  # Step size
x_approx = np.arange(0, 2 * np.pi + h, h)
y_approx = np.zeros(len(x_approx))
y_approx[0] = f(0)  # Initial condition
# Apply Euler's method
for i in range(1, len(x_approx)):
    y_approx[i] = y_approx[i - 1] + h * f_prime(x_approx[i - 1])

# Plotting
plt.figure(figsize=(10, 6))
# Plot the exact function as a filled curve
plt.fill_between(x_exact, y_exact, color='lightblue', alpha=0.5, label='Exact Function $f(x) = \sin(x)$')
# Plot the approximate function
plt.plot(x_approx, y_approx, 'ro-', label='Approximate Function via Euler\'s Method')
# Plot tangent lines at each step
for i in range(len(x_approx) - 1):
    slope = f_prime(x_approx[i])
    x_tangent = np.linspace(x_approx[i], x_approx[i+1], 10)
    y_tangent = y_approx[i] + slope * (x_tangent - x_approx[i])
    plt.plot(x_tangent, y_tangent, 'r--', alpha=0.6)
# Labels and legend
plt.title('Numerical Integration Using Euler\'s Method')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend()
plt.grid(True)

import numpy as np
import matplotlib.pyplot as plt
# Define the original function and its derivatives
def f(x):
    return np.sin(10*x)
def f_prime(x):
    return 10*np.cos(10*x)
def f_double_prime(x):
    return -100*np.sin(10*x)

# Initial conditions
x0 = 0
xf = 2 * np.pi
y0 = f(x0)
# Parameters for adaptive timestep
k = 0.05            # Proportionality constant
epsilon = 1e-9      # Small number to prevent division by zero
h_min = 0.0         # Minimum timestep
h_max = (xf-x0)/100 # Maximum timestep
# Arrays to store the approximated solution
x_approx = [x0]
y_approx = [y0]
h_values = []
x = x0
y = y0

# Adaptive Euler's method
while x < xf:
    # Compute second derivative
    f_dd = f_double_prime(x)   
    # Compute timestep h
    h = k * np.sqrt(np.pi / (np.abs(f_dd) + epsilon))
    h = np.clip(h, h_min, h_max)    
    # Ensure we don't overshoot the final value
    if x + h > xf:
        h = xf - x    
    # Store timestep
    h_values.append(h)    
    # Update x and y using Euler's method
    y += h * f_prime(x)
    x += h    
    # Store the results
    x_approx.append(x)
    y_approx.append(y)
# Exact solution for comparison
x_exact = np.linspace(x0, xf, 400)
y_exact = f(x_exact)

# Plotting
plt.figure(figsize=(12, 6))
# Plot the exact function as a filled curve
plt.fill_between(x_exact, y_exact, color='lightblue', alpha=0.5, label='Exact Function $f(x) = \sin(x)$')
# Plot the approximate function
plt.plot(x_approx, y_approx, 'ro-', label='Approximate Function with Adaptive Timestep')
# Plot tangent lines at each step
for i in range(len(x_approx) - 1):
    slope = f_prime(x_approx[i])
    x_tangent = np.linspace(x_approx[i], x_approx[i+1], 10)
    y_tangent = y_approx[i] + slope * (x_tangent - x_approx[i])
    plt.plot(x_tangent, y_tangent, 'r--', alpha=0.6)
# Labels and legend
plt.title('Numerical Integration with Adaptive Timestep')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend()
plt.grid(True)

Approximate Methods

2.2) Definition and Examples

• Definition:
  • Use approximations to simplify computations.
• Examples:
  • Emulators: Machine learning to interpolate simulation results.
  • Fitting Functions: Empirical formulas derived from data.

Advantages and Limitations

• Advantages:
  • Faster computations.
  • Useful for parameter studies and real-time analysis.
• Limitations:
  • Valid only within certain parameter ranges.
  • Less precise than exact methods.

---

3) Python in Numerical Cosmology

---

3.1) Advantages of Python

• Ease of Use:
  • Simple syntax and readable code.
• Scientific Libraries:
  • NumPy: Efficient array operations.
  • SciPy: Advanced algorithms and functions.
  • Matplotlib: Plotting and visualization.
• Community Support:
  • Active development and support.
  • Cosmology-specific packages (e.g., Astropy, Healpy).

3.2) Performance Considerations

• Interpreted Language:
  • Slower execution compared to compiled languages (C/C++).
• Optimization Techniques:
  • Use efficient libraries (e.g., NumPy, SciPy).
  • Numba: Just-In-Time (JIT) compilation.
  • Cython: Compile Python code to C for speed.
  • Parallel Computing: Utilize multi-threading and multi-processing.

---

When Performance Matters

• Computationally Intensive Tasks:
  • Profile code to find bottlenecks.
  • Integrate Python with faster languages if necessary.

---

4) Hands-On Activity in Python

---

Review of Matter Power Spectrum and Growth Function

4.1) Key Concepts

• Matter Power Spectrum:
  • Describes the distribution of matter density fluctuations as a function of scale.
• Growth Function ( D(a) ):
  • Quantifies the growth of density perturbations over time.
• Key Equation: $$ \ddot{\delta} + 2H\dot{\delta} - \frac{3}{2} H^2 \Omega_m \delta = 0 $$
  • ( $\dot{\delta}$ ): Time derivative of ( $\delta$ ).
  • ( H ): Hubble parameter.
  • ( $\Omega_m$ ): Matter density constribution.

Solving the Growth Function Numerically

4.2) Setup and Implementation

• Define Cosmological Parameters:
  • Let's assume a Flat-$\Lambda$CDM unverse.
  • Is the cosmic time variable the best choice for writing the ODE equations?
  • Which validation and verification strategy we are going to use?
• Initial Conditions:
  • Which initial conditions are apropriate?
• Numerical Solver:
  • Use solve_ivp from SciPy.

Coding Session

• Implementation Steps:
  1. Define the system of equations.
  2. Set initial conditions.
  3. Choose a range for the scale factor ( a ).
  4. Implement the numerical solver.
• Live Coding:
  • Write code step by step.
  • Encourage troubleshooting and debugging.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Define constants
H0 = 70.0  # Hubble parameter at present time in km/s/Mpc
Omega_m = 0.3
Omega_Lambda = 1.0 - Omega_m

# Define the E(a)
def E(a):
    return np.sqrt(Omega_m / a**3 + Omega_Lambda)
# Define the Hubble parameter H(a)
def H(a):
    return H0 * E(a)
# Define Oma(a)
def Om(a):
    return Omega_m / a**3 / E(a)**2
# Define the derivative of H wrt lna
def dlnH_dlna (a):
    return -3/2 / E(a) **2 * Omega_m / a**3 

# Define the system of equations
def growth_function(lna, y):
    a = np.exp(lna) # Get the scale factor
    dydt = [y[1],
            -y[1] * (dlnH_dlna(a)+2) + 3/2 * Om(a) * y[0]]
    return dydt

from astropy.cosmology import FlatLambdaCDM
cosmo = FlatLambdaCDM(H0=H0, Om0=Omega_m, Tcmb0=0)
z = 1
print(cosmo.Om(z)/Om(1/(1+z)), cosmo.H(z).value/H(1/(1+z)))

# Initial conditions and time span
a0 = 1e-2
y0 = [a0, a0]  # Initial values for D and D'
t_span = [np.log(a0), np.log(1)]  # From a ~ 0 to a = 1

# Solve the ODE
sol = solve_ivp(growth_function, t_span, y0, method='RK45', dense_output=True, rtol=1e-6, atol=1e-9)

# Extract results
a = np.linspace(0.1, 1, 100)
D = sol.sol

# Plotting the results
plt.loglog(a, D(np.log(a))[0], label='Numerical Solution')
plt.loglog(a, a, label='EdS')
plt.xlabel('Scale factor a')
plt.ylabel('Growth function D(a)')
plt.legend()
plt.savefig("growth.png")