Compute spectra of bias operators advected from Lagrangian space

[hersle]

WIP. Current plan:

• Create particle fields for $\delta(\mathbf{q})$ and $s^2(\mathbf{q})$.
  • $1(\mathbf{q})$ and $\delta^2(\mathbf{q})$ does not need to be tracked, but can be computed when needed.
• During ICs, when $\delta(\mathbf{q})$ is computed, also compute $s^2(\mathbf{q})$ and store $\delta(\mathbf{q})$ and $s^2(\mathbf{q})$ for each particle.
  • For simplicity begin with $\delta(\mathbf{q})$ only, then add $s^2(\mathbf{q})$ later.
• Compute power spectra on the fly with arbitrary weight/mass $w_i \in [1(\mathbf{q}_i), \delta(\mathbf{q}_i), \delta^2(\mathbf{q}_i), s^2(\mathbf{q}_i)]$ for each particle.
  • Could temp-abuse particle mass field. Save original masses, set mass to desired weight, density-assign particles to grid, then restore masses.
  • Or modify particles_to_grid with a way to specify how the weight is computed.
• Add new input file options for power-spectrum on-the-fly evaluation.
  • Request which spectra $P_{1\delta}$, $P_{\delta\delta}$, ... to compute. Generalize current behavior, which is equivalent to $P_{11}$ only.
• Eventually generalize to multipoles.
• Eventually generalize to redshift space.

---

Notes to self:

•  Should I create separate fields for $1$, $\delta$, $\delta^2$, etc. and interpolate each of them separately to particle positions? Not create a field only for $\delta$, interpolate $\delta$ to each particle, and then calculate $\delta^2 = \delta \delta$ for each particle, as I am currently doing?
• Subtracting the mean (e.g. $<\delta^2> \neq 0$) before computing spectra should not make a difference since it only affects $k=0$? Or must I subtract it before computing spectra?
• What is advected and should be stored in particles during ICs? Is it $\delta^2(q)$ or $\delta^2(q)-<\delta^2(q)>$ that should be associated with each particle? Or is this equivalent?

---

Alternative: output particle weights to GADGET snapshots, then use Python with nbodykit to load GADGET files, paint particles to grid and compute (auto and cross) power spectra

[hersle]

Run as of current changes with the standard example parameterfile:

./nbody parameterfile.lua

echo '
import matplotlib.pyplot as plt 
import numpy as np

filename = "output/snapshot_TestSim_z0.000/pofk_bias_z0.000.txt"
with open(filename, "r") as f:
    cols = f.readline().strip().split()[1:]
    data = np.loadtxt(f)
    data = {col: data[:, i] for i, col in enumerate(cols)} # dict with column name -> column values
    kcol = cols[0]
    for col in cols[1:]:
        plt.plot(np.log10(data[kcol]), np.log10(data[col]), label = col)
    plt.xlabel(f"log10({kcol})")
    plt.ylabel("log10(P/(Mpc/h)^3)")
    plt.legend()
    plt.savefig("plot.png")
' | python3

...

#=====================================================
#
#                .___  _________   
#                |   | \_   ___ \  
#                |   | /    \  \/  
#                |   | \     \____ 
#                |___|  \______  / 
#                              \/  
#
# Generating initial conditions for N-body
# Order in LPT = 2
# The boxsize is 512 comoving Mpc/h
# The initial redshift zini = 20
# Particle has [Position] x_code = x / Box (24 bytes)
# Particle has [Velocity] v_code = a^2 dxdt / (H0 Box) (24 bytes)
# Particle has [Mass] (8 bytes)
# Particle has [ID] (8 bytes)
# Particle has [1LPT Displacement field] (12 bytes)
# Particle has [2LPT Displacement field] (12 bytes)
# Particle has [Lagrangian position] (24 bytes)
# Particle has [Lagrangian overdensity] (8 bytes)
# Particle has [Lagrangian overdensity squared] (8 bytes)
# Total size of particle is 128 bytes
# We will make 128^3 particles
# Plus a buffer with room for 25%% more particles
# We will allocate ~335.544 MB per task for the particles
#
#=====================================================


#=====================================================
#
#            .___        _____          
#            |   | _____/ ____\____     
#            |   |/    \   __\/  _ \    
#            |   |   |  \  | (  <_> )   
#            |___|___|  /__|  \____/    
#                     \/                
#
# Info about MPIParticles. Dimension of particles is 3
# We have allocated 335.544 MB (mean)
#                   335.544 MB (min)
#                   335.544 MB (max)
# Total particles across all tasks is 2097152
# Task 0 has 2097152 parts. Capacity: 2621440
# The buffer is 80 %% (mean)
#               80 %% (min)
#               80 %% (max) filled across tasks
#
#=====================================================

Storing unique ID in particle
Storing Lagrangian position q in particle
Storing Lagrangian overdensity delta(q) in particle
Minimum Lagrangian overdensity delta(q): -0.586401
Maximum Lagrangian overdensity delta(q): 0.564338
Storing Lagrangian overdensity squared delta^2(q) in particle
Minimum Lagrangian overdensity squared delta^2(q): 1.22858e-14
Maximum Lagrangian overdensity squared delta^2(q): 0.343867
Adding 1LPT position to particle
Adding 1LPT velocity to particle
Storing 1LPT displacment field in particle
Maximum 1LPT displacements: 1.74824 Mpc/h
Maximum 1LPT velocity: 450.956 km/s peculiar
Adding 2LPT position to particle
Adding 2LPT velocity to particle
Storing 2LPT displacment field in particle
Maximum 2LPT displacements: 0.039826 Mpc/h
Maximum 2LPT velocity: 20.5219 km/s peculiar

...