SIDM i CONCEPT¶

Setup¶

Opdater CONCEPT:
```
./concept -u update master
```

Lav en backup af koden, så I kan tracke jeres ændringer:

mkdir before_sidm
cp -r *.py concept Makefile fft.c utilities tests before_sidm/

Gem følgende parameterfil som params/collision:

params/collision¶

# Input/output
initial_conditions = paths['ics_dir'] + '/' + basename(paths['params']) + '/IC.hdf5'
output_times = {
    'render2D': linspace(0, 1*Gyr, 15),
}
render2D_select = {
    'particles': {'data': False, 'image': False, 'terminal image': True},
}

# Numerical parameters
boxsize = 100*Mpc

# Graphics
render2D_options = {
    'gridsize': {
        'particles': 40,
    },
    'terminal resolution': {
        'particles': 40,
    },
    'interpolation': {
        'particles': 'NGP',
    },
    'extent': {
        'particles': boxsize,
    },
    'enhance': {
        'particles': False,
    },
}

# Debugging options
enable_Hubble = False

Gem følgende script som ICs/collision/generate_IC.py:

ICs/collision/generate_IC.py¶

# This file has to be run in pure Python mode!

# Imports from the CO𝘕CEPT code
from commons import *
from species import Component
from snapshot import save

# Specifications
distance = 0.8*boxsize
collision_time = 0.5*units.Gyr
mass = 1e+18*units.m_sun

# Create particle positions and momenta
posx = [0.5*(boxsize - distance), 0.5*(boxsize + distance)]
posy = [0.5*boxsize]*2
posz = [0.5*boxsize]*2
mom = mass*(0.5*distance)/collision_time
momx = [+mom, -mom]
momy = [0, 0]
momz = [0, 0]

# Create matter particle component from positions and momenta
component = Component('matter', 'matter', N=len(posx), mass=mass)
component.populate(asarray(posx), 'posx')
component.populate(asarray(posy), 'posy')
component.populate(asarray(posz), 'posz')
component.populate(asarray(momx), 'momx')
component.populate(asarray(momy), 'momy')
component.populate(asarray(momz), 'momz')

# Save snapshot
save(component, initial_conditions)

Generer IC (startbetingelses-) snapshot ved at køre ovenstående script gennem CONCEPT:
```
./concept -p params/collision -m ICs/collision/generate_IC.py --pure-python
```

Aside: Investigering af snapshots

Prøv:

./concept -u info ICs/collision

Installer ViTables via

(source concept && $python -m pip install pyqt5==5.14 vitables)

Kør ViTables via

(source concept && $python -m vitables) &

Drag’n’drop HDF5-filer …

Kør CONCEPT med parameterfilen, som starter fra den genererede IC:

./concept -p params/collision

Genspil serien af 2D renders via

./concept -u play

Indfør `collision` kraft¶

Tildel partiklerne en ny kraft, collision, ved at tilføje følgende til parameterfilen:
```
# Physics
select_forces = {
    'particles': {'collision': 'pp'},
}
```
Kør simuleringen på ny og se hvad der sker.

Indfør collision som en ny kraft ved at tiføje følgende til bunden af interactions.py:

# Collision
register('collision', 'pp')
@cython.pheader(
    # Arguments
    method=str,
    receivers=list,
    suppliers=list,
    ᔑdt=dict,
    interaction_type=str,
    printout='bint',
)
def collision(method, receivers, suppliers, ᔑdt, interaction_type, printout):
    if method == 'pp':
        # The particle-particle method
        if printout:
            masterprint(
                'Executing',
                shortrange_progress_message('collision', method, receivers),
                '...',
            )
        # Carry out interaction
        ...
        if printout:
            # Finish progress message
            masterprint('done')
    elif master:
        abort(f'collision() was called with the "{method}" method')

Kør simuleringen igen.

Tip

Kør evt. koden i pure Python mode (--pure-python) for at undgå rekompilering

Når CONCEPT accepterer den nye kraft, kan vi begynde på den egentlige implementering. Erstat ... fra sidste code-snippet med

component_component(
    'collision', receivers, suppliers, collision_pairwise, ᔑdt,
    pairing_level='domain',
)

Her er collision_pairwise en (endnu uskrevet) funktion hvori selve kraftberegningen skal finde sted. Vi laver denne i et nyt modul, collision.py:

collision.py¶

# Import everything from the commons module.
# In the .pyx file, Cython declared variables will also get cimported.
from commons import *

# Cython imports
cimport('from interactions import particle_particle')



# Function implementing pairwise collision
@cython.header(
    # Arguments
    interaction_name=str,
    receiver='Component',
    supplier='Component',
    ᔑdt_rungs=dict,
    rank_supplier='int',
    only_supply='bint',
    pairing_level=str,
    tile_indices_receiver='Py_ssize_t[::1]',
    tile_indices_supplier_paired='Py_ssize_t**',
    tile_indices_supplier_paired_N='Py_ssize_t*',
    extra_args=dict,
    # Locals
    subtiling_r='Tiling',
    returns='void',
)
def collision_pairwise(
    interaction_name, receiver, supplier, ᔑdt_rungs, rank_supplier, only_supply, pairing_level,
    tile_indices_receiver, tile_indices_supplier_paired, tile_indices_supplier_paired_N,
    extra_args,
):
    # Extract momentum buffers
    momx_r = receiver.momx
    momy_r = receiver.momy
    momz_r = receiver.momz
    momx_s = supplier.momx
    momy_s = supplier.momy
    momz_s = supplier.momz
    # Loop over all (receiver, supplier) particle pairs (i, j)
    for i, j, rung_index_i, rung_index_j, x_ji, y_ji, z_ji, apply_to_i, apply_to_j, particle_particle_t_begin, subtiling_r in particle_particle(
        receiver, supplier, pairing_level,
        tile_indices_receiver, tile_indices_supplier_paired, tile_indices_supplier_paired_N,
        rank_supplier, interaction_name, only_supply,
    ):
        d2 = x_ji**2 + y_ji**2 + z_ji**2
        if d2 > (3*units.Mpc)**2:
            continue
        fancyprint('Collision!', fun=terminal.bold_cyan)

Spørgsmål: Hvorfor d2 > (3*units.Mpc)**2 fremfor sqrt(d2) > 3*units.Mpc?
Det nye modul skal registreres i Makefile under pyfiles.
Husk også at give interactions.py adgang til collision_pairwise():
```
cimport('from collision import *')
```
Check at kollisionen detekteres!

Simpel impuls-udveksling¶

Vi indfører her tilfældig impuls-udveksling uden at overveje $\mathrm{LAB}\leftrightarrow\mathrm{CM}$ transformationen (hvorfor går det godt?) og sørger for at alting går som det skal, inden vi indfører den rigtige impuls-udveksling.

Erstat udprintet i collision_pairwise() med følgende:

# Extract momenta
momx_i = momx_r[i]
momy_i = momy_r[i]
momz_i = momz_r[i]
# Compute collision probability
probability = 1
# Determine whether collision takes place
if random() > probability:
    continue
# Collision takes place!
# Draw collision angles.
θ = arccos(2*random() - 1)
ϕ = 2*π*random()
# Construct transformed momentum according to {θ, ϕ}
mom_i = sqrt(momx_i**2 + momy_i**2 + momz_i**2)
momx_i = sin(θ)*cos(ϕ)*mom_i
momy_i = sin(θ)*sin(ϕ)*mom_i
momz_i = cos(θ)       *mom_i
# Get the other momentum from momentum conservation
momx_j = -momx_i
momy_j = -momy_i
momz_j = -momz_i
# Apply momentum changes
momx_r[i] = momx_i
momy_r[i] = momy_i
momz_r[i] = momz_i
momx_s[j] = momx_j
momy_s[j] = momy_j
momz_s[j] = momz_j

Kør simuleringen. I pure Python mode resulterer random_seed = 2 i “bedre” vinkler.
Prøv nu med to processer (-n 2). Det fejler!
- Den nye collision-kraft skal registreres som indeterministisk i kaldet til register():
```
deterministic=False,
```
  Kør igen.
- Interaktionen finder nu kun sted på den ene proces, men den anden proces får ikke besked herom! Fiks:
  - Under # Extract momentum buffers:
    Δmomx_s = supplier.Δmomx Δmomy_s = supplier.Δmomy Δmomz_s = supplier.Δmomz
  - Under # Extract momenta:
    momx_j_ori = momx_s[j] momy_j_ori = momy_s[j] momz_j_ori = momz_s[j]
  - Under # Apply momentum changes:
    if rank_supplier != rank: Δmomx_s[j] += -momx_j_ori + momx_j Δmomy_s[j] += -momy_j_ori + momy_j Δmomz_s[j] += -momz_j_ori + momz_j
- Den ene partikel opfører sig stadig ikke korrekt. Det skyldes at impuls ikke udveksles mellem prosser som standard (f.eks. afhænger tyngdekraft ikke af impuls, men kun af position). Specificer eksplicit at collision afhænger af både position ('pos') og momentum ('mom'), i kaldet til register():
```
dependent=['pos', 'mom'],
```
  Kør igen.
- Nu virker det, men simuleringen sløver ned efter kollisionen. Det skyldes at den ene partikel (fejlagtigt!) går en “rung” op, hvilket skyldes den pludselige acceleration. Da accelerationen i kollisions-tidspunktet er uendelig er det ikke smart at lade disse påvirke tidsskridtene/rungs. Vi kan registrere collision-kraften som værende instantan ved at tilføje
```
instantaneous=True,
```
  i kaldet til register().
  
  Nu skulle det gerne gå godt uanset antallet af processer!

Korrekt impuls-udveksling¶

Slet linjerne efter θ og ϕ er blevet valgt og ned til # Apply momentum changes. Indfør jeres egen kode og omskriv så det passer ind.

Mulig løsning, til inspiration:

# Transform to center-of-mass frame
mass_r = receiver.mass  # Bør rykkes udenfor loopet!
mass_s = supplier.mass  # Bør rykkes udenfor loopet!
boostx_i = (momx_i + momx_j_ori)*mass_r/(mass_r + mass_s)
boosty_i = (momy_i + momy_j_ori)*mass_r/(mass_r + mass_s)
boostz_i = (momz_i + momz_j_ori)*mass_r/(mass_r + mass_s)
momx_i -= boostx_i
momy_i -= boosty_i
momz_i -= boostz_i
# Construct transformed momentum according to {θ, ϕ},
# rotated correctly.
momx2 = momx_i**2
momy2 = momy_i**2
momxy_i = momx_i*momy_i
tmp = (momz_i - sqrt(momx2 + momy2 + momz_i**2))/(machine_ϵ + momx2 + momy2)
sinθ, cosθ = sin(θ), cos(θ)
sinϕ, cosϕ = sin(ϕ), cos(ϕ)
momx_i, momy_i, momz_i = (
    +sinθ*(sinϕ*momxy_i*tmp + cosϕ*(momz_i - momy2*tmp)) + cosθ*momx_i,
    +sinθ*(cosϕ*momxy_i*tmp + sinϕ*(momz_i - momx2*tmp)) + cosθ*momy_i,
    -sinθ*(cosϕ*momx_i      + sinϕ*momy_i              ) + cosθ*momz_i,
)
# Get the other momentum from momentum conservation
momx_j = -momx_i
momy_j = -momy_i
momz_j = -momz_i
# Transform back to the global box frame
momx_i += boostx_i
momy_i += boosty_i
momz_i += boostz_i
momx_j += mass_s/mass_r*boostx_i
momy_j += mass_s/mass_r*boosty_i
momz_j += mass_s/mass_r*boostz_i

Note

Kollisionen i simuleringen er nu ikke magen til før, men nærmere en spejling heraf! Hvorfor? Er det et problem?

Korrekt kollisions-sandsynlighed¶

Vi vil nu indføre korrekt beregning af probability. Vi benytter følgende formalisme (se e.g. side 26 her):

\[ \begin{align}\begin{aligned}P_{ij} &= n_{ij}|\vec{v}_i - \vec{v}_j|\sigma\Delta t\,,\\n_{ij} &\equiv W_i(\vec{r}) * W_j(\vec{r}) \,,\\W_i(\vec{r}) &\equiv r_{\mathrm{s}_i}^{-3} w\biggl(\frac{\vec{r} - \vec{r}_i}{r_{\mathrm{s}_i}}\biggr) \,.\end{aligned}\end{align} \]

Vi vælger (nogenlunde arbitrært) at bruge en Gauss for partikel-formen (“particle shape”/”smoothing kernel”):

\[ \begin{align}\begin{aligned}w(\vec{x}) &= w(|\vec{x}|) = w(x) \equiv (2\pi)^{-3/2}\exp\biggl(-\frac{x^2}{2}\biggr)\\\Rightarrow n_{ij} &= \bigl[2\pi(r_{\mathrm{s}_i}^2 + r_{\mathrm{s}_j}^2)\bigr]^{-3/2} \exp\biggl(-\frac{d^2}{2(r_{\mathrm{s}_i}^2 + r_{\mathrm{s}_j}^2)}\biggr) \,,\end{aligned}\end{align} \]

hvor $d = |\vec{r}_i - \vec{r}_j|$ er afstanden mellem partikel $i$ og $j$.

Tværsnittet er parameteriseret som tværsnit/masse, $(\sigma/m)$. Jeg kan da se to muligheder, alt efter om $m$ skal betyde den samlede masse eller massen af én partikel:

\[\begin{split}\sigma_{ij} = (\sigma/m)\begin{cases} m_i + m_j \\ (m_i + m_j)/2 \end{cases}\end{split}\]

Fra dette lader det til at vi skal bruge $\sigma_{ij} = (\sigma/m)(m_i + m_j)/2$.

Tilbage til koden:

Start med at definere $(r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j})$ og $(\sigma/m)$ ovenfor loopet:

# The sum of squared collision lengths
collision_length_r = 0.1*boxsize
collision_length_s = 0.1*boxsize
collision_length2 = collision_length_r**2 + collision_length_s**2
# Collision cross section per mass
σ_per_mass = 10*units.cm**2/units.g

Redefiner nu probability til

probability = collision_factor*exp(-d2/(2*collision_length2))*vel_relative

Opgave til jer: Kod collision_factor og vel_relative. I får masserne som receiver.mass og supplier.mass. Tidsskridtet $\Delta t$ er lidt specielt; I kan bruge Δt = ᔑdt_rungs['1'][0], for nu.

Mulig løsning nedenfor, skulle det blive nødvendigt:

# Collision cross section
σ_per_mass = 10*units.cm**2/units.g
mass_r = receiver.mass
mass_s = supplier.mass
σ = σ_per_mass*(mass_r + mass_s)/2
# Factor used for the collision probability.
# We have
#   P = factor*exp(-r²/(2*collision_length2))|v⃗ᵢ - v⃗ⱼ|
# with
#   factor = (2π*collision_length2)**(-3/2)*σ*Δt
collision_factor = (2*π*collision_length2)**(-1.5)*σ*ᔑdt_rungs['1'][0]

.
.
.

    vel_relative = sqrt(
        + (momx_i/mass_r - momx_j_ori/mass_s)**2
        + (momy_i/mass_r - momy_j_ori/mass_s)**2
        + (momz_i/mass_r - momz_j_ori/mass_s)**2
    )

Kør en simulering. I skulle gerne få en sandsynlighed til kollisions-tidsskridtet på $18\%$. Jeg fandt at kollison finder sted for random_seed = 6.
Nu hvor collision_length2 er defineret kan vi indføre den i stedet for de hard-coded (3*units.Mpc)**2. Vi kan f.eks. bruge
```
if d2 > 16*collision_length2:
   continue
```
hvilket svarer til $d > 4\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}}$. Allerede ved $3.76\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}}$ har man $3\sigma$ ($99.7\%$) med.

Tværsnit/masse som parameter¶

Lad os forfremme σ_per_mass til en parameter vi kan specificere i parameterfilen. Vi gør det generelt, således at én værdi af $(\sigma/m)$ er knyttet til et par af komponenter. Således kan (baryoner $\leftrightarrow$ baryoner), (baryoner $\leftrightarrow$ mørk stof) og (mørk stof $\leftrightarrow$ mørk stof) alle have forskellige $(\sigma/m)$.

Indfør denne nye parameter til parameterfilen:
```
select_collision_crosssection_per_mass = {
    ('particles', 'particles'): 10*cm**2/g,
}
```
Køres CONCEPT nu vil den fortælle at der er givet en ukendt parameter.

Indfør følgende i bunden af “Physics”-kategorien i parameter-indlæsningen i commons.py for at tilføje select_collision_crosssection_per_mass som en ny parameter:

select_collision_crosssection_per_mass = {}
if user_params.get('select_collision_crosssection_per_mass'):
    if isinstance(user_params['select_collision_crosssection_per_mass'], dict):
        select_collision_crosssection_per_mass = user_params['select_collision_crosssection_per_mass']
        replace_ellipsis(select_collision_crosssection_per_mass)
    else:
        select_collision_crosssection_per_mass = {
            'all'             : user_params['select_collision_crosssection_per_mass'],
            'all combinations': user_params['select_collision_crosssection_per_mass'],
        }
user_params['select_collision_crosssection_per_mass'] = select_collision_crosssection_per_mass

Tilføj også cython-deklarationen

select_collision_crosssection_per_mass=dict,

ovenover.

Vi kunne hive $(\sigma/m)$ ud af select_collision_crosssection_per_mass inde i collision_pairwise(), men det er smartere at gøre det allerede i collision(). Tilføj følgende ovenover # Carry out interaction:

# Get collision cross section per mass
# for each {receiver, supplier} pair.
extra_args = {
    'σ_per_mass': get_pairwise_quantities(
        receivers, suppliers, select_collision_crosssection_per_mass, 'σ_per_mass',
    ),
}

og tilføj

interaction_extra_args=extra_args,

som argument til component_component().

Inde i collision_paiwise() fås σ_per_mass da ved

σ_per_mass = extra_args['σ_per_mass'][receiver.name, supplier.name]

Check om I nu kan styre om der sker kollision eller ej ved at justere på select_collision_crosssection_per_mass i parameterfilen.

Aside: Tracking af ændringer

I kan se præcis hvad I har ændret i en fil ved at køre

file=commons.py; git diff --no-index before_sidm/$file $file

Kollisions-længde som parameter¶

Vi indfører nu kollisions-længden/skalaen $r_{\mathrm{s}}$ som en parameter. Vi gør det på en lidt speciel måde (I ser hvorfor når vi når til tiling).

Indfør parameteren

shortrange_params = {
    'collision': {
        'scale': {
            'particles': 0.1*boxsize,
        },
    },
}

i parameterfilen.

Erstat definitionerne af collision_length_r og collision_length_s i collision_pairwise() med

collision_length_r = get_shortrange_param(receiver, 'collision', 'scale')
collision_length_s = get_shortrange_param(supplier, 'collision', 'scale')

Check at det virker.

Kollisions-rækkevidde som parameter¶

Vi har p.t.

if d2 > 16*collision_length2:
   continue

svarende til en hard-coded rækkevidde på $4\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}}$. Lad os også forfremme denne rækkevidde til en parameter.

Indfør
```
'range': {
    ('particles', 'particles'): '3.763*sqrt(scale_i**2 + scale_j**2)',  # 3σ
},
```
i parameterfilen under shortrange_params['colision'].

Note

De 3.763 fås fra

\[\begin{split}\int_0^R n_{ij}(d) \, 4\pi d^2\, \mathrm{d}d \equiv \mathrm{erf}\biggl(\frac{k}{\sqrt{2}}\biggr) \Rightarrow R = \begin{cases} 1.878\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}} & k = 1 \\ 2.833\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}} & k = 2 \\ 3.763\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}} & k = 3 \end{cases}\end{split}\]

Her svarer $k$ til antallet af $\sigma$ (“standardafvigelser”) man medtager. Med $k = 1$ medtager man altså $\mathrm{erf}(1/\sqrt{2}) = 68.3\%$ af kollisionerne, med $k = 2$ medtager man $\mathrm{erf}(2/\sqrt{2}) = 95.4\%$, med $k = 3$ medtager man $\mathrm{erf}(3/\sqrt{2}) = 99.7\%$, etc.

Indfør collision_range2 i collision_pairwise() som:

collision_range2 = get_shortrange_param((receiver, supplier), 'collision', 'range')**2

Udskfit de hard-codede 16*collision_length2 med collision_range2.
Check at tingene stadig fungerer.

Kosmologisk simulering¶

Nu hvor alting virker og intet er hard-coded kan vi skifte vores simple simulering/parameterfil ud med en fuld kosmologisk simulering/parameterfil.

Gem følgende nye parameterfil som params/collision_cosmo:

params/collision_cosmo¶

# Non-parameter variable used to control the size of the simulation
_size = 20

# Input/output
initial_conditions = [
    {
        'species': 'matter',
        'N'      : _size**3,
    },
]
output_dirs = {
    'render3D': paths['output_dir'] + '/' + basename(paths['params']),
}
output_times = {
    'render3D': [0.02, 0.1, 0.3, 1],
}

# Numerical parameters
boxsize = _size*Mpc
shortrange_params = {
    'collision': {
        'scale': {
            'particles': '0.06*boxsize/cbrt(N)',
        },
        'range': {
            ('particles', 'particles'): '3.763*sqrt(scale_i**2 + scale_j**2)',  # 3σ
        },
    },
}

# Cosmology
H0 = 67*km/(s*Mpc)
Ωcdm = 0.27
Ωb = 0.049
a_begin = 0.02

# Physics
select_forces = {
    'particles': {
        'gravity'  : ('p3m', 2*_size),
        'collision': 'pp',
    },
}
select_collision_crosssection_per_mass = {
    ('particles', 'particles'): 10.0*cm**2/g,
}

# Debugging options
particle_reordering = False

Den store forskel fra før er, at vi nu har mange flere partikler (hvis ICs genereres som en del af selve simuleringen), Hubble-ekspansionen er ikke længere deaktiveret, tyngdekraft er sat til, og kollisions-længden er nu sat til '0.06*boxsize/cbrt(N)', altså $6\%$ af middel-partikel-separationen.

Note

De 0.06 fremkommer således:

På side 26 i denne (og andre steder) påstås det at man bør have $r_{\text{s}} \approx 0.2 L/\sqrt[3]{N}$ (eller større), hvor $L$ er boxsize. Her er $w$ dog ikke en Gauss, men en såkaldt B-spline på formen

\[\begin{split}w(x) = \frac{8}{\pi}\begin{cases} 1 - 6x^2 + 6x^3 & 0 \le x < 1/2 \\ 2(1 - x)^3 & 1/2 \le x < 1 \\ 0 & 1 \le x\,. \end{cases}\end{split}\]

Nedenunfer er plottet vores Gauss $W(r)$ og deres B-spline $W(r)$. Gaussen er plottet to gange; én hvor $r_{\mathrm{s}_{\mathrm{G}}} = r_{\mathrm{s}_{\mathrm{B}}}$ (‘G’ for Gauss, ‘B’ for B-spline), og én hvor $r_{\mathrm{s}_{\mathrm{G}}} = 0.280r_{\mathrm{s}_{\mathrm{B}}}$.

Det ses at med samme $r_{\mathrm{s}}$ er B-splinen langt mere lokaliseret end Gaussen (husk at begge $W$’er er normerede!). Bruges en mindre $r_{\mathrm{s}}$ for Gaussen end for B-splinen — omkring 0.280 gange så stor — kan vi få dem til at ligne hinanden.

De foreslåede $r_{\text{s}} \approx 0.2 L/\sqrt[3]{N}$ for B-splinen svarer da til $r_{\text{s}} \approx 0.2 L/\sqrt[3]{N}\times 0.280 \approx 0.06L/\sqrt[3]{N}$ for Gaussen.

Selvom øjemål er ligeså godt, er den nøjagtige værdi på 0.280 fundet ved at kræve at Gaussen “indeholder $1\sigma$ af partiklen” indenfor samme fysiske radius som B-splinen gør det ved:

\[\begin{split}\mathrm{erf}\biggl(\frac{1}{\sqrt{2}}\biggr) \equiv \int_0^R W(r) \, 4\pi r^2 \, \mathrm{d}r \Rightarrow R = \begin{cases} 0.525r_{\mathrm{s}_{\mathrm{B}}} & \text{B-spline} \\ 1.88r_{\mathrm{s}_{\mathrm{G}}} & \mathrm{Gauss}\,, \end{cases}\end{split}\]

\[R_{\mathrm{G}} \equiv R_{\mathrm{B}} \Rightarrow r_{\mathrm{s}_{\mathrm{G}}} = 0.280r_{\mathrm{s}_{\mathrm{B}}}\]

Køres collision_cosmo i pure Python mode går det meget langsomt, hvorfor vi er nødt til at skifte til kompileret tilstand. Da vi har tilføjet et nyt modul er det sikreste at lave en komplet rekompilering:
```
make clean && ./concept
```
Påbegynd en simulering af collision_cosmo i compiled mode, bare for at se det køre.

Kollisions-counter¶

Vi indfører nu en counter der tæller antallet af kollisioner i hvert tidsskridt.

Indfør en variabel n_interactions i collision_pairwise() som tæller kollisionerne op.
Vi kan ikke bare udskrive n_interactions direkte fra collision_pairwise(), da denne funktion kaldes mange gange per tidsskridt. Tilføj i stedet
```
# Add number of interactions to running total
receiver.n_interactions[supplier.name] += n_interactions
```
efter loopet.

Køres simuleringen nu vil n_interactions blive udskrevet efter de fleste “kicks”, men ikke alle (se e.g. tidsskridt 8). Tilføj følgende ovenfor # Finish progress message i collision() for at rette op på dette:

# Print out number of interactions
masterprint('Interaction count:')
lines = []
pairs = []
for receiver in receivers:
    for supplier in suppliers:
        pair = {receiver, supplier}
        if pair in pairs:
            continue
        pairs.append(pair)
        n_interactions = allreduce(receiver.n_interactions[supplier.name], op=MPI.SUM)
        lines.append(f'{receiver.name} ⟷ {supplier.name}: ${n_interactions}')
masterprint('\n'.join(align_text(lines, indent=4)))

Kør nu simuleringen til ende. Antallet af kollisioner per tidsskridt skulle gerne være 0 i starten og stige med tiden.
Tip

Et hurtigt overblik over antallet af kollisioner kan fås ved
```
grep '⟷\|factor' logs/<ID>
```
Hvis simuleringen stadig er igang er følgende smartere:
```
tail -f -n +1 logs/<ID> | grep '⟷\|factor'
```
Tip

Det totale antal kollisioner i en simulering kan fås ved
```
grep '⟷' logs/<ID> | awk '{s+=$NF}END{print s}'
```
Er man for doven til manuelt at erstatte <ID> i kommandoen ovenfor kan man definere
```
count(){ grep '⟷' logs/$(ls -t logs | head -n 1) | awk '{s+=$NF}END{print s}'; }
```
hvorefter man bare kan køre
```
count
```
for at få det samlede antal kollisioner i den seneste simulering.

Periodicitet¶

Lad en kosmologisk simulering køre til ende og find det totale antal kollisioner.
Spørgsmål: Hvad er den størst mulige afstand mellem to partikler, i en given dimension?

Implementer periodicitet i afstandene x_ji, y_ji og z_ji i collision_pairwise():

# Translate coordinates so they correspond to the nearest image
if x_ji > 0.5*boxsize:
    x_ji -= boxsize
elif x_ji < -0.5*boxsize:
    x_ji += boxsize
if y_ji > 0.5*boxsize:
    y_ji -= boxsize
elif y_ji < -0.5*boxsize:
    y_ji += boxsize
if z_ji > 0.5*boxsize:
    z_ji -= boxsize
elif z_ji < -0.5*boxsize:
    z_ji += boxsize

Kør endnu en simulering og find det totale antal kollisioner. Skal antallet blive større eller mindre nu hvor periodicitet er indbygget?

Skalafaktoren¶

Nu hvor Hubble-ekspansionen er slået til (vi har ikke længere enable_Hubble = False) varierer $a(t)$ gennem simuleringen, hvor den før var $a=1$ hele vejen igennem. Det er automatisk taget højde for idet koden anvender comoving koordinater og impulser. Dog skal hastigheden $|\vec{v}_i - \vec{v}_j|$ i $P_{ij}$ være fysisk og ikke comoving, men vi beregner den ud fra comoving impuls $\vec{q}$. Vi skal dermed have

\[\vec{v}_i = \frac{\vec{q}_i}{m_i} \rightarrow \frac{\vec{q}_i}{am_i} \, .\]

Note

Tænk på det således: Comoving impuls $\vec{q}_i$ er per konstruktion konstant for en partikel (som ikke påvirkes af nogen kræfter). Men den fysiske impuls rødforskydes/mindskes over tid grundet universets udvidelse. Ergo må $\vec{q}_i/a$ være den fysiske impuls.

Den nuværende værdi af skalafaktoren er altid tilgængelig i CONCEPT som universals.a. Altså kan vi bare indføre en faktor 1/universals.a i collision_factor.

MEN da universals.a er værdien af $a(t)$ i starten af tidsskridtet, bliver collision_factor en smule for stor. Det havde været mere fair at bruge værdien som $a(t)$ antog i midten af tidsskridtet. Allerbedst bør man integrere over tidsskridtet. Dvs:

\[a^{-1}\Delta t \rightarrow \int_t^{t + \Delta t} a^{-1}(t)\, \mathrm{d}t \,.\]

Variablen ᔑdt_rungs er i virkeligheden en stor samling af diverse integraler over det nuværende tidsskridt. ᔑdt_rungs['1'] er præcis $\int_t^{t+\Delta t}1\,\mathrm{d}t = \Delta t$. For at indføre $a^{-1}(t)$ korrekt midlet over tidsskridtet skal vi da blot erstatte ᔑdt_rungs['1'] med ᔑdt_rungs['a**(-1)']. (Det efterfølgende [0] skal beholdes. I ser hvad det er i næste sektion).

Indfør ᔑdt_rungs['a**(-1)'][0] i collision_factor og genkør simuleringen. Forventer I nu flere eller færre kollisioner?

Rungs¶

Som I kan se fra simuleringerne bor forskellige partikler på forskellige “rungs” ($0, 1, 2, \dots$). En partikel på rung $n$ har i virkeligheden et tidsskridt på

\[\Delta t_n = \frac{\Delta t}{2^n} \,.\]

Dette giver anledning til forskellige collision_factor for partikler på forskellige rungs. Det ekstra [0] i ᔑdt_rungs['a**(-1)'][0] betyder netop at vi plukker integralet ud svarende til rung 0. En partikel på rung 1 får da fejlagtigt en $\sim$ dobbelt så stor sandsynlighed for at kollidere, som den egentlig skulle, og endnu værre for partikler på højere rungs.

Vi fikser det ved at fjerne [0] og omdøbe collision_factor til collision_factors.

I beregningen af probability skal vi da vælge den korrekte rung/indgang i collision_factors. Dette skal være den største af partikel $i$’s og partikel $j$’s rung (hvorfor)?:

rung_index_max = (rung_index_i if rung_index_i > rung_index_j else rung_index_j)
probability = (
    collision_factors[rung_index_max]
    *exp(-d2/(2*collision_length2))
    *vel_relative
)

Genkør simuleringen. I skulle nu få langt færre kollisioner.

Tiles¶

I dette afsnit ændrer vi algoritmen hvorved partikel-partikel-parringen finder sted, for at opnå bedre performance. Bemærk kørsels-tiden for den seneste simuleringer.

Udskift pairing_level='domain' med pairing_level='tile' i collision().
Tilføj
```
'tilesize': 2*boxsize/_size,
```
under shortrange_params['collision'] i parameterfilen.
Kør simuleringen igen. Det skulle gerne gå hurtigere nu. Havde simuleringen været større ville effekten værre tydeligere.

Note

De 2 i 'tilesize': 2*boxsize/_size er totalt arbitrære. Forskellige værdier leder til forskellig performance. I bør teste (ikke nu) hvilken faktor der får en stor simulering til at køre hurtigt, uden at kræve urimeligt meget ekstra hukommelse.

Definer default-værdier for shortrange_params['collision'] ved at tilføje

'collision': {
    'scale'    : '0.06*boxsize/cbrt(N)',
    'range'    : '3.763*sqrt(scale_i**2 + scale_j**2)',  # 3σ
    'tilesize' : boxsize/3,
    'subtiling': 'automatic',
},

til shortrange_params_defaults i commons.py.

For at få automatisk subtiling (refinement) til at fungere, skal følgende tilføjes som det allersidste i collision_pairwise():

# Add computation time to the running total,
# for use with automatic subtiling refinement.
if j != -1:
    particle_particle_t_final = time()
    subtiling_r.computation_time += particle_particle_t_final - particle_particle_t_begin

og j = -1 skal indsættes ovenfor loopet.

Generaliseret spredningsfordeling¶

Den isotrope spredning (θ = arccos(2*random() - 1), ϕ = 2*π*random()) er hard-coded. Selvom vi pt. ikke har andre spredningsfordelinger i tankerne, er det godt at få en mere generel mekanisme bygget ind.

Tilføj følgende nye parameter til parameterfilen:

select_collision_distribution = {
    ('particles', 'particles'): 'isotropic',
}

Tilføj ligeledes select_collision_distribution til commons.py, på præcis samme måde som vi gjorde det for select_collision_crosssection_per_mass, dog med én tilføjelse:
```
select_collision_distribution.setdefault('all', 'isotropic')
```
Tilføj 'distribution' til extra_args i collision() på samme måde som 'σ_per_mass', og hiv ligeledes distribution ud af extra_args i collision_pairwise() som vi gjordet det for σ_per_mass.

Koden der trækker de to vinkler skal nu byttes ud med

# Draw collision angles
if distribution == 'isotropic':
    θ = arccos(2*random() - 1)
    ϕ = 2*π*random()
else:
    θ = ϕ = 0  # To satisfy the compiler
    abort(
        f'Collision distribution "{distribution}" '
        f'not implemented in collision_pairwise()'
    )

Timestep limiter¶

Vores sandsynligheds-formalisme kræver at probability $P_{ij}$ er additiv, sådan at den totale sandsynlighed for at partikel $i$ kolliderer er

\[P_i = \sum_{j\neq i} P_{ij}\]

Da sandsynligheder kun er (approksimativt) additive når de er små, skal vi sørge for at $\Delta t$ er lille nok til at dette er tilfældet. I den her har de som krav at $P_{ij} < 0.2$.

Jeg har testet og endda bygget en timestep limiter (krav for hvornår en partikel skal hoppe en rung op) ind så tidsskridtet for hver partikel overholder ovenstående, men det viste sig at $P_{ij}$ aldrig kommer tæt på $0.2$ (for rimelige valg af $(\sigma/m)$), så jeg har fjernet det igen.

Konklusion: Vi er heldige at den nye kollisions-kraft foregår på en tidsskala der er større end den partiklerne allerede opdaters med i forvejen, så vi skal ikke tænke yderligere over det. I kan prøve en gang at printe probability ud og se hvad dens største værdi er.

Optimering¶

I dette afsnit indfører vi diverse kode-optimeringer uden at ændre på selve algoritmen, for at opnå bedre performance.

Indfør typer med hjælp fra collision.html.
Tip

Hurtig rekompilering kan fås ved
```
./concept --unsafe
```
Udskift opslag i collision_factors-arrayet med et råt pointer-opslag i
```
collision_factors_ptr = cython.address(collision_factors[:])
```

Indfør caching af mellemresultater via ℝ[]. De vigtigste steder er:

2*π $\rightarrow$ ℝ[2*π].
Periodicitets-beregningen: ℝ[0.5*boxsize] og ℝ[-0.5*boxsize].
Den relative hastighed: /mass_r $\rightarrow$ *ℝ[1/mass_r], samme for s.
Sandsynligheden: -d2/(2*collision_length2) $\rightarrow$ d2*ℝ[-1/(2*collision_length2)].
Boost-definition: ℝ[mass_r/(mass_r + mass_s)].

Transformerede impuls $i$:

tmp = (momz_i - sqrt(ℝ[momx2 + momy2] + momz_i**2))/(machine_ϵ + ℝ[momx2 + momy2])

samt

momx_i, momy_i, momz_i = (
    +sinθ*(sinϕ*ℝ[momxy_i*tmp] + cosϕ*(momz_i - momy2*tmp)) + cosθ*momx_i,
    +sinθ*(cosϕ*ℝ[momxy_i*tmp] + sinϕ*(momz_i - momx2*tmp)) + cosθ*momy_i,
    -sinθ*(cosϕ*momx_i         + sinϕ*momy_i              ) + cosθ*momz_i,
)

Boost: ℝ[mass_s/mass_r].

Unswitch distribution == 'isotropic' og if rank_supplier != rank:.
Fjern particle_reordering = False fra parameterfilen.
Kør simuleringen (kræver rekompilering). Giver optimeringerne en synlig effekt?

Halo finding¶

Sæt _size op til 32 i parameterfilen.
Tilføj snapshot outputs til parameterfilen.
Tilføj
```
snapshot_type = 'gadget2'
```
til parameterfilen
Kør en simulering for at producere nogle snapshots.

Gem nedenstående script i concept-mappen som run_rockstar (og gør eksekverbar med chmod +x run_rockstar):

run_rockstar¶

#!/usr/bin/env bash

# This script installs (if not found) and runs Rockstar on supplied GADGET2 snapshots.
# You might want to correct some constants in rockstar/universal_constants.h:
#   Gc                    : G*(m_sun/Mpc)/(km/s)**2
#   VMAX_CONST            : sqrt(G*m_sun/Mpc)/(km/s)
#   CRITICAL_DENSITY      : 3750/(π*G)*Mpc/m_sun*(km/s)**2
#   HUBBLE_TIME_CONVERSION: 1/(100*km/(s*Mpc)*yr)

set -e
this_file="$(readlink -f "${BASH_SOURCE[0]}")"
this_dir="$(dirname "${this_file}")"

cd "${this_dir}"
source concept

# Get and build Rockstar if not already present
rockstar_dir="${top_dir}/rockstar"
if [ ! -f "${rockstar_dir}/rockstar" ]; then
    rockstar_url="https://bitbucket.org/gfcstanford/rockstar.git"
    git clone "${rockstar_url}" "${rockstar_dir}"
    cd "${rockstar_dir}"
    make
fi
cd "${this_dir}"
rockstar="${rockstar_dir}/rockstar"

# Handle command-line arguments
snapshot=""
softening=0.03
MIN_HALO_OUTPUT_SIZE=20
for arg in $@; do
    if     [ "${arg}" == "-h"   ] || [ "${arg}" == "--h"    ] || [ "${arg}" == "--he" ] \
        || [ "${arg}" == "-hel" ] || [ "${arg}" == "--help" ]; then
        echo "Usage: ./$(basename "${this_file}") [softening=<softening length in simulation in units of mean interparticle distance, defaults to 0.03>] [MIN_HALO_OUTPUT_SIZE=<min number of particles in halos, defaults to 20>] </path/to/snapshot(s)>"
        exit 0
    elif [[ "${arg}" == "softening"* ]]; then
        softening=$("${python}" -B -c "print('${arg}'.split('=')[-1])")
    elif [[ "${arg}" == "MIN_HALO_OUTPUT_SIZE"* ]]; then
        MIN_HALO_OUTPUT_SIZE=$("${python}" -B -c "print('${arg}'.split('=')[-1])")
    elif [ -z "${snapshot}" ]; then
        snapshot="${arg}"
    else
        echo "Incorrect input format! Supply -h for help." >&2
        exit 1
    fi
done
if [ -z "${snapshot}" ]; then
    echo "You must supply a snapshot or a directory of snapshots" >&2
    exit 1
fi
if [ -d "${snapshot}" ]; then
    # Locate GADGET2 snapshots in directory
    for f in "${snapshot}/"*; do
        is_gadget2=$("${python}" -B -c "
import struct
is_gadget2 = False
try:
    with open('${f}', 'rb') as f:
        f.seek(4)
        head = struct.unpack('4s', f.read(struct.calcsize('4s')))
        is_gadget2 = (head[0] == b'HEAD')
except:
    pass
print(is_gadget2)
")
        if [[ "${is_gadget2}" == "True" ]]; then
            "${this_file}" softening=$softening MIN_HALO_OUTPUT_SIZE=$MIN_HALO_OUTPUT_SIZE "${f}"
        fi
    done
    exit 0
elif [ ! -f "${snapshot}" ]; then
    echo "The snapshot \"${snapshot}\" does not exist" >&2
    exit 1
fi
snapshot="$(absolute_path "${snapshot}")"
printf "\n${esc_bold}Halo finding within \"${snapshot}\"${esc_normal}\n"

# Compute the FORCE_RES parameter,
# the force resolution of the simulation in Mpc/h.
info="$("${concept}" -u info ${snapshot} --pure)"
boxsize=$(echo "${info}" | grep "BoxSize" | awk '{print $NF}')  # kpc/h
boxsize="$("${python}" -B -c "print(${boxsize}*1e-3)")"  # Mpc/h
N=$(echo "${info}" | grep "Npart" | awk '{print $3}')
N="${N/,/}"
FORCE_RES=$("${python}" -B -c "print(${softening}*(${boxsize})/${N}**(1/3))")
echo "{softening = ${softening}, boxsize = ${boxsize} Mpc/h, N = ${N}} ⟶ FORCE_RES = ${FORCE_RES}"

# Generate Rockstar parameter file
snapshot_dir="$(dirname "${snapshot}")"
OUTBASE="${rockstar_dir}/output"
rm -rf "${OUTBASE}"
mkdir "${OUTBASE}"
params="${OUTBASE}/params.cfg"
echo "
# Input
FILE_FORMAT = \"GADGET2\"
INBASE = \"${snapshot_dir}\"
FILENAME =  \"$(basename "${snapshot}")\"
GADGET_LENGTH_CONVERSION = 1e-3   # kpc/h        (GADGET) -> Mpc/h  (Rockstar)
GADGET_MASS_CONVERSION   = 1e+10  # 1e+10*Msun/h (GADGET) -> Msun/h (Rockstar)

# Output
OUTBASE = \"${OUTBASE}\"
OUTPUT_FORMAT = \"ASCII\"
DELETE_BINARY_OUTPUT_AFTER_FINISHED = 1

# Numerics
FORCE_RES = ${FORCE_RES}  # Force resolution of simulation in Mpc/h
MIN_HALO_OUTPUT_SIZE = ${MIN_HALO_OUTPUT_SIZE}
RESCALE_PARTICLE_MASS = 1  # Rescale cdm particle masses as to include baryons

# Parallelization (needed for periodic boundaries)
PARALLEL_IO = 1
NUM_READERS = 1
NUM_WRITERS = 8  # number of CPUs, must be 1 or multiple of 8, BUT setting it to 1 disables periodic boundaries
FORK_READERS_FROM_WRITERS = 1
FORK_PROCESSORS_PER_MACHINE = 8  # set equal to NUM_WRITERS
" > "${params}"

# Run Rockstar on the snapshot
cd "${rockstar_dir}"
./$(basename "${rockstar}") -c "${params}" & pid=$!  # master
sleep_time=0.05
sleeps_max=200
sleeps=0
while [ ! -f "${OUTBASE}/auto-rockstar.cfg" ]; do
    sleep ${sleep_time}
    ((sleeps += 1))
    if [ ${sleeps} -gt ${sleeps_max} ]; then
        echo "Could not start writer processes as \"${OUTBASE}/auto-rockstar.cfg\" could not be found" >&2
        kill $pid
        exit 1
    fi
done
./$(basename "${rockstar}") -c "${OUTBASE}/auto-rockstar.cfg"  # writers

# Move result next to input snapshot
halo_file="${snapshot}_halos"
mv "${OUTBASE}/out_0.list" "${halo_file}"
echo "Halo file generated: \"${halo_file}\""

set +e

Kør Rockstar på de producerede snapshots således:
```
./run_rockstar output/collision_cosmo
```
Dette generere en halo-fil per snapshot.

Gem nedenstående (quick and dirty) plotte-script i output/collision_cosmo:

plot_halos.py¶

import glob, re, warnings
import numpy as np
import matplotlib.pyplot as plt

# What to plot?
quantity, unit = 'Rvir', 'kpc/$h$'

# Get file names and the scale factor values
filenames = glob.glob('*_halos')
a_values = [
    float(re.search('a=(.+)_halos', filename).group(1))
    for filename in filenames
]
a_values, filenames = zip(*sorted(zip(a_values, filenames)))

# Read in header
with open(filenames[0]) as f:
    header = f.readline().strip('# ').split()

# Plot quantity for each time
fig, axes = plt.subplots(1, len(a_values), figsize=(3*len(a_values), 3))
for a, filename, ax in zip(a_values, filenames, axes):
    ax.set_title(rf'$a = {a}$')
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')
        data = np.loadtxt(filename, usecols=[header.index(quantity)])
    if data.size == 0:
        continue
    ax.hist(data, 15)
    ax.set_xlabel(unit)

# Save figure
fig.suptitle(quantity)
plt.tight_layout()
plt.savefig('halos.pdf')

Kør plotte-scriptet via
```
(source concept && cd output/collision_cosmo && $python plot_halos.py)
```
Både størrelsen og antallet af haloer skulle gerne vokse over tid.
Som det står øverst i run_rockstar er det værd lige at rette værdierne af diverse konstanter i Rockstar’s universal_constants.h, så de stemmer overens med værdierne brugt af CONCEPT. Som eksempel findes den første værdi (Gc — “gravitational constant”) via
```
./concept -m 'from commons import *; print(G_Newton*(units.m_sun/units.Mpc)/(units.km/units.s)**2)' --pure
```
Når alle fire værdier er rettet skal Rockstar rekompileres:
```
(cd ../rockstar && make clean && make)
```
Genkør Rockstar, omdøb halos.pdf så den ikke overskrives og genkør plotte-scriptet. De to halo-plots skulle gerne være så godt som identiske.

Baryoner og koldt mørkt stof for sig¶

Indtil nu har vores partikler representeret matter, dvs. baryoner og koldt mørt stof, slået sammen.

Separer de to species ud ved at ændre startbetingelserne i parameterfilen til

initial_conditions = [
    {
        'species': 'baryons',
        'N'      : _size**3,
    },
    {
        'species': 'cold dark matter',
        'N'      : _size**3,
    },
]

Vi skal nu opdatere select_forces så det kun er koldt mørkt stof der får kollisionskraften:

select_forces = {
    'baryons': {
        'gravity': ('p3m', 2*_size),
    },
    'cold dark matter': {
        'gravity'  : ('p3m', 2*_size),
        'collision': 'pp',
    },
}

Man burde nok også ændre de mange 'particles' til 'cold dark matter' for at gøre det klart at disse ikke gælder for baryoner, men da baryoner alligevel ikke deltager i kollisioskraften gør det ingen forskel.
Rockstar kan desværre ikke læse snapshots med flere komponenter, så vi er nødt til kun at gemme den ene. Tilføj
```
snapshot_select = {
    'cold dark matter': True,
    'baryons'         : False,
}
```
til parameterfilen.
Kør simuleringen, genkør Rockstar, omdøb halos.pdf så den ikke bliver overskrevet og genkør plotte-scriptet. Plottet skulle gerne være lidt anderledes end før.

Note

Rockstar er blevet sat op til at tage højde for baryonerne ved at skalere mørk stof partiklernes masse op. Det er ikke perfekt, men det er bedre end helt at ignorere baryonernes eksistens.

Collision switch¶

For at se effekten af kollisions-kraften er det selvfølgelig nødvendigt at køre tilsvarende simuleringer uden denne kraft. Her indfører vi en command-line-parameter _collision, som bestemmer om simuleringen skal være med eller uden kollision.

Ændre select_forces i parameterfilen til

select_forces = {
    'particles': {
        'gravity': ('p3m', 2*_size),
    },
}
if _collision:
    select_forces.update(
        {
            'cold dark matter': {
                'collision': 'pp',
            },
        }
    )

Opdater også output_dirs så output fra simuleringer med og uden kollision ikke overskriver hinanden:
```
paths['output_dir'] + '/' + basename(paths['params']) + f'/collision={_collision}',
```

Nu skal simuleringen køres enten som

./concept -p params/collision_cosmo -c "_collision = True"

eller

./concept -p params/collision_cosmo -c "_collision = False"

Kør en fuld simulering både med og uden kollision. Hvad er forskellen i beregningstid?
Kør Rockstar på de nye snapshots.

Erstat plotte-scriptet med dette:

plot_halos.py¶

import glob, re, warnings
import numpy as np
import matplotlib.pyplot as plt

# What to plot?
quantity, unit = 'Rvir', 'kpc/$h$'

# Get file names and theu scale factor values
filenames = glob.glob('collision=False/*_halos')
a_values = [
    float(re.search('a=(.+)_halos', filename).group(1))
    for filename in filenames
]
a_values, filenames = zip(*sorted(zip(a_values, filenames)))

# Read in header
with open(filenames[0]) as f:
    header = f.readline().strip('# ').split()

# Plot quantity for each time
fig, axes = plt.subplots(1, len(a_values), figsize=(3*len(a_values), 3))
for a, filename, ax in zip(a_values, filenames, axes):
    ax.set_title(rf'$a = {a}$')
    for collision, color in {False: 'C0', True: 'C1'}.items():
        with warnings.catch_warnings():
            warnings.simplefilter('ignore')
            data = np.loadtxt(
                filename.replace('collision=False', f'collision={collision}'),
                usecols=[header.index(quantity)],
            )
        if data.size == 0:
            continue
        ax.hist(data, 15, color=color, alpha=0.4,
            label=('with' + 'out'*(not collision) + ' collisions'),
        )
        ax.set_xlabel(unit)
ax.legend()
fig.suptitle(quantity)
plt.tight_layout()
plt.savefig('halos.pdf')

og kør plotte-scriptet.

Tip

Man kan køre begge simuleringer (uden og med kollisioner) samt Rockstar og plotte-sciptet i én kommando således:

for collision in False True; do
    ./concept -p params/collision_cosmo -c "_collision = $collision"
    ./run_rockstar output/collision_cosmo/collision=$collision
done && (source concept && cd output/collision_cosmo && $python plot_halos.py)

Ting der kan/bør testses¶

Science¶

Hvordan ændres power spectra og halo-plots med valget af $(\sigma/m)$?
- Hvilke halo-egenskaber (masse, radius, eccentricitet, …) er interessante at kigge på?
- Kan vi reproducere eksisterende resultater i literaturen?
Undersøg andre spredningsfordelinger end den isotrope.
Undersøg hastigheds-afhængigt tværsnit.

Numerics¶

Vi ignorerer partikelpar med en afstand $d > 3.763\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}}$ idet vi forventer at dette kun negligerer $0.3\%$ af alle potentielle kollisioner.
- Passer dette? Får vi $\sim 0.3\%$ flere kollisioner hvis vi bruger e.g. $10\sqrt{r^2_{\mathrm{s}_i} + r^2_{\mathrm{s}_j}}$?
- Er de 3.763, svarende til $3\sigma$, rimelige? Kunne man nøjes med $2\sigma$, og ville dette øge computeringshastigheden mærktbart?
Valget $r_{\mathrm{s}} = 0.06L/\sqrt[3]{N}$. Hvordan ser halo-plots (og power spectra) ud med $r_{\mathrm{s}} \in \{0.03, 0.04, \dots , 1.00\}L/\sqrt[3]{N}$? Da dette er en ufysisk parameter skulle resultaterne helst være ufølsomme overfor det præcise valg af denne!
Der bør findes en (nogenlunde) optimal måde at vælge 'tilesize' på. Her kan potentielt være meget beregningstid at spare.
Hvad sker der hvis vi har forskelligt antal partikler for baryoner og koldt mørkt stof? Lige nu har vi lige mange af hver species, hvilket leder til at baryon-partiklerne er $\sim 5$ gange lettere, hvilket potentielt ikke er godt?
Kan man nøjes med bare at have én komponent, matter, men sænke $(\sigma/m)$ således at det effektivt svarer til at det kun er “koldt mørkt stof”-delen af partiklerne der besidder kollisions-kraften?
Hvor stokastisk er simuleringerne? Dvs. hvor stor forskel er der i resultaterne, hvis man ændrer random_seed? (Her bør man nok starte fra et snapshot, så valget af random_seed kun påvirker tidsuviklingen (via kollisioner), ikke startbetingelserne).
For fastlåst boxsize, hvor mange partikler (hvor stor _size) skal der til for at opnå konvergens i power spectra og halo-plots (for en given $k_{\mathrm{max}}$ eller tilsvarende $r_{\mathrm{min}}$)?

SIDM i CONCEPT¶

Setup¶

Aside: Investigering af snapshots

Indfør collision kraft¶

Simpel impuls-udveksling¶

Korrekt impuls-udveksling¶

Korrekt kollisions-sandsynlighed¶

Tværsnit/masse som parameter¶

Aside: Tracking af ændringer

Kollisions-længde som parameter¶

Kollisions-rækkevidde som parameter¶

Kosmologisk simulering¶

Kollisions-counter¶

Periodicitet¶

Skalafaktoren¶

Rungs¶

Tiles¶

Generaliseret spredningsfordeling¶

Timestep limiter¶

Optimering¶

Halo finding¶

Baryoner og koldt mørkt stof for sig¶

Collision switch¶

Ting der kan/bør testses¶

Science¶

Numerics¶

Indfør `collision` kraft¶