Cosmic Structure Formation

Source code in C (Click to expand)

#include <math.h>
#include <stdio.h>
#include <hdf5.h>

#include "grav_sim.h"

#define INIT_CONDITION_FILE "../../ics_swift.hdf5"
#define A_FINAL 1.0 // Scale factor at the end of simulation
#define NUM_STEPS 1000

ErrorStatus read_init_condition(
    CosmologicalSystem *__restrict system,
    int *__restrict pm_grid_size
);

int main(void)
{
    CosmologicalSystem system;
    int pm_grid_size;
    ErrorStatus error_status = WRAP_TRACEBACK(read_init_condition(
        &system,
        &pm_grid_size
    ));
    if (error_status.return_code != GRAV_SUCCESS)
    {
        goto error;
    }

    /* Output parameters */
    OutputParam output_param = get_new_output_param();
    output_param.method = OUTPUT_METHOD_HDF5;
    output_param.output_dir = "snapshots/";
    output_param.output_interval = (A_FINAL - system.scale_factor) / 100.0;
    output_param.output_initial = true;
    output_param.coordinate_output_dtype = OUTPUT_DTYPE_FLOAT;
    output_param.velocity_output_dtype = OUTPUT_DTYPE_FLOAT;
    output_param.mass_output_dtype = OUTPUT_DTYPE_FLOAT;

    /* Simulation status */
    SimulationStatus simulation_status;

    /* Settings */
    Settings settings = get_new_settings();
    settings.verbose = GRAV_VERBOSITY_NORMAL;

    int return_code = launch_cosmological_simulation(
        &system,
        &output_param,
        &simulation_status,
        &settings,
        A_FINAL,
        NUM_STEPS,
        pm_grid_size
    );
    if (return_code != 0)
    {
        free_cosmological_system(&system);
        error_status = WRAP_RAISE_ERROR(
            GRAV_FAILURE,
            "Simulation failed."
        );
        goto error;
    }

    /* Free memory */
    free_cosmological_system(&system);

    return 0;

error:
    print_and_free_traceback(&error_status);
    return 1;
}

ErrorStatus read_init_condition(
    CosmologicalSystem *__restrict system,
    int *__restrict pm_grid_size
)
{
    ErrorStatus error_status;

    /* Open file */
    hid_t file = H5Fopen(INIT_CONDITION_FILE, H5F_ACC_RDONLY, H5P_DEFAULT);
    if (file == H5I_INVALID_HID)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to open initial condition file");
        goto err_open_file;
    }

    /* Open groups */
    hid_t header = H5Gopen(file, "/Header", H5P_DEFAULT);
    hid_t units = H5Gopen(file, "/Units", H5P_DEFAULT);
    hid_t part_type_1 = H5Gopen(file, "/PartType1", H5P_DEFAULT);
    if (header == H5I_INVALID_HID
        || units == H5I_INVALID_HID
        || part_type_1 == H5I_INVALID_HID
    )
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to open groups in initial condition file");
        goto err_open_group;
    }

    /* Read header attributes */
    hid_t num_particles_attr = H5Aopen(header, "NumPart_ThisFile", H5P_DEFAULT);
    hid_t box_size_attr = H5Aopen(header, "BoxSize", H5P_DEFAULT);
    hid_t grid_size_attr = H5Aopen(header, "suggested_pmgrid", H5P_DEFAULT);
    hid_t omega_m_attr = H5Aopen(header, "Omega0", H5P_DEFAULT);
    hid_t omega_lambda_attr = H5Aopen(header, "OmegaLambda", H5P_DEFAULT);
    hid_t redshift_attr = H5Aopen(header, "Redshift", H5P_DEFAULT);
    hid_t h_attr = H5Aopen(header, "HubbleParam", H5P_DEFAULT);
    hid_t mass_table_attr = H5Aopen(header, "MassTable", H5P_DEFAULT);
    if (
        num_particles_attr == H5I_INVALID_HID
        || box_size_attr == H5I_INVALID_HID
        || grid_size_attr == H5I_INVALID_HID
        || omega_m_attr == H5I_INVALID_HID
        || omega_lambda_attr == H5I_INVALID_HID
        || redshift_attr == H5I_INVALID_HID
        || h_attr == H5I_INVALID_HID
        || mass_table_attr == H5I_INVALID_HID
    )
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to open attributes in header");
        goto err_header_attr;
    }

    uint32_t temp_num_particles_arr[6];
    if (H5Aread(num_particles_attr, H5T_NATIVE_UINT32, &temp_num_particles_arr) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read number of particles from header");
        goto err_header_attr;
    }
    double box_size;
    if (H5Aread(box_size_attr, H5T_NATIVE_DOUBLE, &box_size) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read box size from header");
        goto err_header_attr;
    }
    uint32_t temp_grid_size;
    if (H5Aread(grid_size_attr, H5T_NATIVE_UINT32, &temp_grid_size) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read grid size from header");
        goto err_header_attr;
    }
    *pm_grid_size = temp_grid_size;

    double omega_m;
    if (H5Aread(omega_m_attr, H5T_NATIVE_DOUBLE, &omega_m) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read omega_m from header");
        goto err_header_attr;
    }
    double omega_lambda;
    if (H5Aread(omega_lambda_attr, H5T_NATIVE_DOUBLE, &omega_lambda) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read omega_lambda from header");
        goto err_header_attr;
    }
    double redshift;
    if (H5Aread(redshift_attr, H5T_NATIVE_DOUBLE, &redshift) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read redshift from header");
        goto err_header_attr;
    }
    double h;
    if (H5Aread(h_attr, H5T_NATIVE_DOUBLE, &h) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read Hubble parameter from header");
        goto err_header_attr;
    }

    error_status = WRAP_TRACEBACK(get_initialized_cosmological_system(
        system,
        temp_num_particles_arr[1]
    ));
    if (error_status.return_code != GRAV_SUCCESS)
    {
        goto err_initialize_system;
    }
    system->box_width = box_size / 2.0;
    system->box_center[0] = system->box_width;
    system->box_center[1] = system->box_width;
    system->box_center[2] = system->box_width;
    system->scale_factor = 1.0 / (redshift + 1.0);
    system->omega_m = omega_m;
    system->omega_lambda = omega_lambda;
    system->h = h;

    /* Read units attributes */
    hid_t unit_length_attr = H5Aopen(units, "Unit length in cgs (U_L)", H5P_DEFAULT);
    hid_t unit_mass_attr = H5Aopen(units, "Unit mass in cgs (U_M)", H5P_DEFAULT);
    hid_t unit_time_attr = H5Aopen(units, "Unit time in cgs (U_t)", H5P_DEFAULT);

    if (
        unit_length_attr == H5I_INVALID_HID
        || unit_mass_attr == H5I_INVALID_HID
        || unit_time_attr == H5I_INVALID_HID
    )
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to open attributes in units");
        goto err_units_attr;
    }

    if (H5Aread(unit_length_attr, H5T_NATIVE_DOUBLE, &(system->unit_length)) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read unit length from units");
        goto err_units_attr;
    }
    if (H5Aread(unit_mass_attr, H5T_NATIVE_DOUBLE, &(system->unit_mass)) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read unit mass from units");
        goto err_units_attr;
    }
    if (H5Aread(unit_time_attr, H5T_NATIVE_DOUBLE, &(system->unit_time)) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read unit time from units");
        goto err_units_attr;
    }

    /* Open particle datasets */
    hid_t part_type_1_masses = H5Dopen(part_type_1, "Masses", H5P_DEFAULT);
    hid_t part_type_1_coordinates = H5Dopen(part_type_1, "Coordinates", H5P_DEFAULT);
    hid_t part_type_1_velocities = H5Dopen(part_type_1, "Velocities", H5P_DEFAULT);
    if (
        part_type_1_masses == H5I_INVALID_HID
        || part_type_1_coordinates == H5I_INVALID_HID
        || part_type_1_velocities == H5I_INVALID_HID
    )
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to open datasets in PartType1");
        goto err_open_particle_dataset;
    }

    /* Read particle data */
    if (H5Dread(part_type_1_masses, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, system->m) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read masses from PartType1");
        goto err_read_particle_data;
    }
    if (H5Dread(part_type_1_coordinates, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, system->x) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read coordinates from PartType1");
        goto err_read_particle_data;
    }
    if (H5Dread(part_type_1_velocities, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, system->v) < 0)
    {
        error_status = WRAP_RAISE_ERROR(GRAV_OS_ERROR, "Failed to read velocities from PartType1");
        goto err_read_particle_data;
    }

    /* Close HDF5 objects */
    H5Dclose(part_type_1_masses);
    H5Dclose(part_type_1_coordinates);
    H5Dclose(part_type_1_velocities);

    H5Aclose(unit_length_attr);
    H5Aclose(unit_mass_attr);
    H5Aclose(unit_time_attr);

    H5Aclose(num_particles_attr);
    H5Aclose(box_size_attr);
    H5Aclose(grid_size_attr);
    H5Aclose(omega_m_attr);
    H5Aclose(omega_lambda_attr);
    H5Aclose(redshift_attr);
    H5Aclose(h_attr);
    H5Aclose(mass_table_attr);

    H5Gclose(header);
    H5Gclose(units);
    H5Gclose(part_type_1);

    H5Fclose(file);

    return make_success_error_status();

err_read_particle_data:
err_open_particle_dataset:
    H5Dclose(part_type_1_masses);
    H5Dclose(part_type_1_coordinates);
    H5Dclose(part_type_1_velocities);
err_units_attr:
    H5Aclose(unit_length_attr);
    H5Aclose(unit_mass_attr);
    H5Aclose(unit_time_attr);
    free_cosmological_system(system);
err_initialize_system:
err_header_attr:
    H5Aclose(num_particles_attr);
    H5Aclose(box_size_attr);
    H5Aclose(grid_size_attr);
    H5Aclose(omega_m_attr);
    H5Aclose(omega_lambda_attr);
    H5Aclose(redshift_attr);
    H5Aclose(h_attr);
    H5Aclose(mass_table_attr);
err_open_group:
    H5Gclose(header);
    H5Gclose(units);
    H5Gclose(part_type_1);
err_open_file:
    H5Fclose(file);

    return error_status;
}

Source code in Python (Click to expand)

from pathlib import Path

import h5py
import numpy as np
from grav_sim import GravitySimulatorAPI

INIT_CONDITION_FILE = Path(__file__).parent.parent / "ics_swift.hdf5"
A_FINAL = 1.0  # Scale factor at the end of simulation
NUM_STEPS = 1000


def main() -> None:
    gs = GravitySimulatorAPI()

    system = gs.get_new_cosmological_system()
    pm_grid_size = read_init_condition(system)

    _, _, output_param, settings = gs.get_new_parameters()
    output_param.method = "hdf5"
    output_param.output_dir = Path(__file__).parent / "snapshots/"
    output_param.output_interval = (A_FINAL - system.scale_factor) / 100
    output_param.output_initial = True
    output_param.coordinate_output_dtype = "float"
    output_param.velocity_output_dtype = "float"
    output_param.mass_output_dtype = "float"

    settings.verbose = "normal"

    gs.launch_cosmological_simulation(
        system,
        output_param,
        settings,
        A_FINAL,
        NUM_STEPS,
        pm_grid_size,
    )


def read_init_condition(system) -> int:
    with h5py.File(INIT_CONDITION_FILE, "r") as f:
        header = f["/Header"]
        units = f["/Units"]
        part_type_1 = f["/PartType1"]

        # Read header attributes
        num_particles = int(header.attrs["NumPart_ThisFile"][1])
        box_size = float(header.attrs["BoxSize"])
        pm_grid_size = int(header.attrs["suggested_pmgrid"])
        omega_m = float(header.attrs["Omega0"])
        omega_lambda = float(header.attrs["OmegaLambda"])
        redshift = float(header.attrs["Redshift"])
        h_param = float(header.attrs["HubbleParam"])

        system.box_width = box_size / 2.0
        system.box_center = np.array(
            [system.box_width, system.box_width, system.box_width], dtype=np.float64
        )
        system.scale_factor = 1.0 / (redshift + 1.0)
        system.omega_m = omega_m
        system.omega_lambda = omega_lambda
        system.h = h_param

        # Read unit attributes
        system.unit_length_in_cgs = float(units.attrs["Unit length in cgs (U_L)"])
        system.unit_mass_in_cgs = float(units.attrs["Unit mass in cgs (U_M)"])
        system.unit_time_in_cgs = float(units.attrs["Unit time in cgs (U_t)"])

        # Read particle datasets
        system.add(
            part_type_1["Coordinates"][:],
            part_type_1["Velocities"][:],
            part_type_1["Masses"][:],
        )

        return pm_grid_size


if __name__ == "__main__":
    main()

Related topics:

• Comoving coordinates
• Particle mesh algorithm

In this example, we will simulate the formation of cosmological structures in a universe with \(\Lambda\)CDM model.

The initial conditions were generated with MUSIC v2¹ at \(z = 50\), with a resolution of \(128^3 \sim 2\) million particles in a 30 Mpc \(/ h\) box. The MUSIC configuration file and initial conditions are provided in the repository, but with a smaller particles count of \(64^3\) to reduce the repository size.

To run the simulation, we used the LeapFrog integrator in comoving coordinates, with a particle mesh algorithm with mesh size = \(256^3 \sim 16\) million cells (following the suggestions from the MUSIC initial conditions).

Gallery

The simulation was done on Macbook Air M1. The simulation ran surprisingly fast and is finished in less than 10 minutes. The visualization was done with gadgetviewer as it is directly compatible with our HDF5 output format.

---

1. Oliver Hahn and Tom Abel. Multi-scale initial conditions for cosmological simulations. Mon. Not. Roy. Astron. Soc., 415:2101–2121, 2011. arXiv:1103.6031, doi:10.1111/j.1365-2966.2011.18820.x. ↩

Comments