proto.h File Reference

#include <bigfile.h>
#include <stdint.h>
#include "power.h"
#include "allvars.h"
#include "libgadget/petapm.h"

Include dependency graph for proto.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	IDGenerator

Typedefs
typedef struct IDGenerator	IDGenerator

Functions
void	idgen_init (IDGenerator *idgen, PetaPM *pm, int Ngrid, double BoxSize)
uint64_t	idgen_create_id_from_index (IDGenerator *idgen, int index)
void	idgen_create_pos_from_index (IDGenerator *idgen, int index, double pos[3])
void	displacement_fields (PetaPM *pm, enum TransferType Type, struct ic_part_data *dispICP, const int NumPart, Cosmology *CP, const struct genic_config GenicConfig)
int	setup_grid (IDGenerator *idgen, double shift, double mass, struct ic_part_data *ICP)
int	setup_glass (IDGenerator *idgen, PetaPM *pm, double shift, int seed, double mass, struct ic_part_data *ICP, const double UnitLength_in_cm, const char *OutputDir)
void	glass_evolve (PetaPM *pm, int nsteps, const char *pkoutname, struct ic_part_data *ICP, const int NumPart, const double UnitLength_in_cm, const char *OutputDir)
void	saveheader (BigFile *bf, int64_t TotNumPartCDM, int64_t TotNumPartGas, int64_t TotNuPart, double nufrac, const double BoxSize, Cosmology *CP, const struct genic_config GenicConfig)
void	compute_mass (double *mass, int64_t TotNumPartCDM, int64_t TotNumPartGas, int64_t TotNuPart, double nufrac, const double BoxSize, Cosmology *CP, const struct genic_config GenicConfig)
void	write_particle_data (IDGenerator *idgen, const int Type, BigFile *bf, const uint64_t FirstID, const int SavePrePos, int NumFiles, int NumWriters, struct ic_part_data *curICP)
void	read_parameterfile (char *fname, struct genic_config *GenicConfig, int *ShowBacktrace, double *MaxMemSizePerNode, Cosmology *CP)
void	_bigfile_utils_create_block_from_c_array (BigFile *bf, void *baseptr, const char *name, const char *dtype, size_t dims[], ptrdiff_t elsize, int NumFiles, int NumWriters, MPI_Comm comm)

Typedef Documentation

IDGenerator

typedef struct IDGenerator IDGenerator

Function Documentation

_bigfile_utils_create_block_from_c_array()

void _bigfile_utils_create_block_from_c_array	(	BigFile *	bf,
		void *	baseptr,
		const char *	name,
		const char *	dtype,
		size_t	dims[],
		ptrdiff_t	elsize,
		int	NumFiles,
		int	NumWriters,
		MPI_Comm	comm
	)

Definition at line 18 of file save.c.

{
    BigBlock block;
    BigArray array;
    BigBlockPtr ptr;
    ptrdiff_t strides[2];
    int64_t TotNumPart;
    MPI_Allreduce(&dims[0], &TotNumPart, 1, MPI_INT64, MPI_SUM, comm);
 
    strides[1] = dtype_itemsize(dtype);
    strides[0] = elsize;
 
    big_array_init(&array, baseptr, dtype, 2, dims, strides);
 
    if(0 != big_file_mpi_create_block(bf, &block, name, dtype, dims[1], NumFiles, TotNumPart, comm)) {
        endrun(0, "%s:%s\n", big_file_get_error_message(), name);
    }
 
    if(0 != big_block_seek(&block, &ptr, 0)) {
        endrun(0, "Failed to seek:%s\n", big_file_get_error_message());
    }
 
    if(0 != big_block_mpi_write(&block, &ptr, &array, NumWriters, comm)) {
        endrun(0, "Failed to write :%s\n", big_file_get_error_message());
    }
 
    if(0 != big_block_mpi_close(&block, comm)) {
        endrun(0, "%s:%s\n", big_file_get_error_message(), name);
    }
}

References endrun(), MPI_INT64, name, and TotNumPart.

Referenced by save_transfer(), and saveblock().

Here is the call graph for this function:

Here is the caller graph for this function:

compute_mass()

void compute_mass	(	double *	mass,
		int64_t	TotNumPartCDM,
		int64_t	TotNumPartGas,
		int64_t	TotNuPart,
		double	nufrac,
		const double	BoxSize,
		Cosmology *	CP,
		const struct genic_config	GenicConfig
	)

Definition at line 90 of file save.c.

{
    double UnitTime_in_s = GenicConfig.units.UnitLength_in_cm / GenicConfig.units.UnitVelocity_in_cm_per_s;
    double Grav = GRAVITY / pow(GenicConfig.units.UnitLength_in_cm, 3) * GenicConfig.units.UnitMass_in_g * pow(UnitTime_in_s, 2);
    const double OmegatoMass = 3 * CP->Hubble * CP->Hubble / (8 * M_PI * Grav) * pow(BoxSize, 3);
    double OmegaCDM = CP->Omega0;
    mass[0] = mass[2] = mass[3] = mass[4] = mass[5] = 0;
    if (TotNumPartGas > 0) {
        mass[0] = (CP->OmegaBaryon) * 3 * CP->Hubble * CP->Hubble / (8 * M_PI * Grav) * pow(BoxSize, 3) / TotNumPartGas;
        OmegaCDM -= CP->OmegaBaryon;
    }
    if(CP->MNu[0] + CP->MNu[1] + CP->MNu[2] > 0) {
        double OmegaNu = get_omega_nu(&CP->ONu, 1);
        OmegaCDM -= OmegaNu;
        if(TotNuPart > 0) {
            mass[2] = nufrac * OmegaNu * OmegatoMass / TotNuPart;
        }
    }
    mass[1] = OmegaCDM * OmegatoMass / TotNumPartCDM;
}

References CP, get_omega_nu(), GRAVITY, Cosmology::Hubble, Cosmology::MNu, Cosmology::Omega0, Cosmology::OmegaBaryon, Cosmology::ONu, UnitSystem::UnitLength_in_cm, UnitSystem::UnitMass_in_g, genic_config::units, and UnitSystem::UnitVelocity_in_cm_per_s.

Referenced by main(), and saveheader().

Here is the call graph for this function:

Here is the caller graph for this function:

displacement_fields()

void displacement_fields	(	PetaPM *	pm,
		enum TransferType	Type,
		struct ic_part_data *	dispICP,
		const int	NumPart,
		Cosmology *	CP,
		const struct genic_config	GenicConfig
	)

Definition at line 150 of file zeldovich.c.

                                                                                                                                                                       {
 
    /*MUST set this before doing force.*/
    ptype = Type;
    curICP = dispICP;
    PetaPMParticleStruct pstruct = {
        curICP,
        sizeof(curICP[0]),
        ((char*) &curICP[0].Pos[0]) - (char*) curICP,
        ((char*) &curICP[0].Mass) - (char*) curICP,
        NULL,
        NULL,
        NumPart,
    };
 
    int i;
    #pragma omp parallel for
    for(i = 0; i < NumPart; i++)
    {
        memset(&curICP[i].Disp[0], 0, sizeof(curICP[i].Disp));
        memset(&curICP[i].Vel[0], 0, sizeof(curICP[i].Vel));
        curICP[i].Density = 0;
    }
 
 
    /* This reads out the displacements into P.Disp and the velocities into P.Vel.
     * Disp is used to avoid changing the particle positions mid-way through.
     * Note that for the velocities we do NOT just use the velocity transfer functions.
     * The reason is because of the gauge: velocity transfer is in synchronous gauge for CLASS,
     * newtonian gauge for CAMB. But we want N-body gauge, which we get by taking the time derivative
     * of the synchronous gauge density perturbations. See arxiv:1505.04756*/
    PetaPMFunctions functions[] = {
        {"Density", density_transfer, readout_density},
        {"DispX", disp_x_transfer, readout_disp_x},
        {"DispY", disp_y_transfer, readout_disp_y},
        {"DispZ", disp_z_transfer, readout_disp_z},
        {"VelX", vel_x_transfer, readout_vel_x},
        {"VelY", vel_y_transfer, readout_vel_y},
        {"VelZ", vel_z_transfer, readout_vel_z},
        {NULL, NULL, NULL },
    };
 
    /*Set up the velocity pre-factors*/
    const double hubble_a = hubble_function(CP, GenicConfig.TimeIC);
 
    double vel_prefac = GenicConfig.TimeIC * hubble_a;
 
    if(GenicConfig.UsePeculiarVelocity) {
        /* already for peculiar velocity */
        message(0, "Producing Peculiar Velocity in the output.\n");
    } else {
        vel_prefac /= sqrt(GenicConfig.TimeIC); /* converts to Gadget velocity */
    }
 
    if(!GenicConfig.PowerP.ScaleDepVelocity) {
        vel_prefac *= F_Omega(CP, GenicConfig.TimeIC);
        /* If different transfer functions are disabled, we can copy displacements to velocities
         * and we don't need the extra transfers.*/
        functions[4].name = NULL;
    }
 
    int Nregions;
    struct ic_prep_data icprep = {dispICP, NumPart};
    PetaPMRegion * regions = petapm_force_init(pm,
           makeregion,
           &pstruct,
           &Nregions,
           &icprep);
 
    /*This allocates the memory*/
    pfft_complex * rho_k = petapm_alloc_rhok(pm);
 
    gaussian_fill(pm->Nmesh, petapm_get_fourier_region(pm),
          rho_k, GenicConfig.UnitaryAmplitude, GenicConfig.InvertPhase, GenicConfig.Seed);
 
    petapm_force_c2r(pm, rho_k, regions, Nregions, functions);
 
    myfree(rho_k);
    myfree(regions);
    petapm_force_finish(pm);
 
    double maxdisp = 0, maxvel = 0;
 
    #pragma omp parallel for reduction(max:maxdisp, maxvel)
    for(i = 0; i < NumPart; i++)
    {
        int k;
        double absv = 0;
        for(k = 0; k < 3; k++)
        {
            double dis = curICP[i].Disp[k];
            if(dis > maxdisp)
                maxdisp = dis;
            /*Copy displacements to positions.*/
            curICP[i].Pos[k] += curICP[i].Disp[k];
            /*Copy displacements to velocities if not done already*/
            if(!GenicConfig.PowerP.ScaleDepVelocity)
                curICP[i].Vel[k] = curICP[i].Disp[k];
            curICP[i].Vel[k] *= vel_prefac;
            absv += curICP[i].Vel[k] * curICP[i].Vel[k];
            curICP[i].Pos[k] = periodic_wrap(curICP[i].Pos[k], pm->BoxSize);
        }
        if(absv > maxvel)
            maxvel = absv;
    }
    MPI_Allreduce(MPI_IN_PLACE, &maxdisp, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
    message(0, "Type = %d max disp = %g in units of cell sep %g \n", ptype, maxdisp, maxdisp / (pm->BoxSize / pm->Nmesh) );
 
    MPI_Allreduce(MPI_IN_PLACE, &maxvel, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
    message(0, "Max vel=%g km/s, vel_prefac= %g  hubble_a=%g fom=%g \n", sqrt(maxvel), vel_prefac, hubble_a, F_Omega(CP, GenicConfig.TimeIC));
 
    walltime_measure("/Disp/Finalize");
    MPIU_Barrier(MPI_COMM_WORLD);
}

References PetaPM::BoxSize, CP, curICP, ic_part_data::Density, density_transfer(), ic_part_data::Disp, disp_x_transfer(), disp_y_transfer(), disp_z_transfer(), F_Omega(), functions, gaussian_fill(), hubble_function(), genic_config::InvertPhase, makeregion(), ic_part_data::Mass, message(), MPIU_Barrier, myfree, PetaPMFunctions::name, PetaPM::Nmesh, ic_prep_data::NumPart, periodic_wrap(), petapm_alloc_rhok(), petapm_force_c2r(), petapm_force_finish(), petapm_force_init(), petapm_get_fourier_region(), ic_part_data::Pos, genic_config::PowerP, ptype, readout_density(), readout_disp_x(), readout_disp_y(), readout_disp_z(), readout_vel_x(), readout_vel_y(), readout_vel_z(), power_params::ScaleDepVelocity, genic_config::Seed, genic_config::TimeIC, genic_config::UnitaryAmplitude, genic_config::UsePeculiarVelocity, ic_part_data::Vel, vel_x_transfer(), vel_y_transfer(), vel_z_transfer(), and walltime_measure.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

glass_evolve()

void glass_evolve	(	PetaPM *	pm,
		int	nsteps,
		const char *	pkoutname,
		struct ic_part_data *	ICP,
		const int	NumPart,
		const double	UnitLength_in_cm,
		const char *	OutputDir
	)

Definition at line 73 of file glass.c.

{
    int i;
    int step = 0;
    double t_x = 0;
    double t_v = 0;
    double t_f = 0;
 
    /*Allocate memory for a power spectrum*/
    powerspectrum_alloc(pm->ps, pm->Nmesh, omp_get_max_threads(), 0, pm->BoxSize*UnitLength_in_cm);
 
    glass_force(pm, t_x, ICP, NumPart);
 
    /* Our pick of the units ensures there is an oscillation period of 2 * M_PI.
     *
     * (don't ask me how this worked -- I think I failed this problem in physics undergrad. )
     * We use 4 steps per oscillation, and end at
     *
     * 12 + 1 = 13, the first time phase is M_PI / 2, a close encounter to the minimum.
     *
     * */
    for(step = 0; step < nsteps; step++) {
        /* leap-frog, K D D F K */
        double dt = M_PI / 2; /* step size */
        double hdt = 0.5 * dt; /* half a step */
        int d;
        /*
         * Use inverted gravity with a damping term proportional to the velocity.
         *
         * The magic setup was studied in
         *      https://github.com/rainwoodman/fastpm-python/blob/1be020b/Example/Glass.ipynb
         * */
 
        /* Kick */
        for(i = 0; i < NumPart; i ++) {
            for(d = 0; d < 3; d ++) {
                /* mind the damping term */
                ICP[i].Vel[d] += (ICP[i].Disp[d] - ICP[i].Vel[d]) * hdt;
            }
        }
        t_x += hdt;
 
        /* Drift */
        for(i = 0; i < NumPart; i ++) {
            for(d = 0; d < 3; d ++) {
                ICP[i].Pos[d] += ICP[i].Vel[d] * dt;
           }
        }
        t_v += dt;
 
        glass_force(pm, t_x, ICP, NumPart);
        t_f = t_x;
 
        /* Kick */
        for(i = 0; i < NumPart; i ++) {
            for(d = 0; d < 3; d ++) {
                /* mind the damping term */
                ICP[i].Vel[d] += (ICP[i].Disp[d] - ICP[i].Vel[d]) * hdt;
            }
        }
 
        t_x += hdt;
        message(0, "Generating glass, step = %d, t_f= %g, t_v = %g, t_x = %g\n", step, t_f / (2 * M_PI), t_v / (2 *M_PI), t_x / (2 * M_PI));
        glass_stats(ICP, NumPart);
 
        /*Now save the power spectrum*/
        powerspectrum_save(pm->ps, OutputDir, pkoutname, t_f, 1.0);
    }
 
    /*We are done with the power spectrum, free it*/
    powerspectrum_free(pm->ps);
}

References PetaPM::BoxSize, ic_part_data::Disp, glass_force(), glass_stats(), message(), PetaPM::Nmesh, ic_part_data::Pos, powerspectrum_alloc(), powerspectrum_free(), powerspectrum_save(), PetaPM::ps, UnitLength_in_cm, and ic_part_data::Vel.

Referenced by main(), and setup_glass().

Here is the call graph for this function:

Here is the caller graph for this function:

idgen_create_id_from_index()

uint64_t idgen_create_id_from_index	(	IDGenerator *	idgen,
		int	index
	)

Definition at line 69 of file zeldovich.c.

{
    int i = index / (idgen->size[2] * idgen->size[1]) + idgen->offset[0];
    int j = (index % (idgen->size[1] * idgen->size[2])) / idgen->size[2] + idgen->offset[1];
    int k = (index % idgen->size[2]) + idgen->offset[2];
    uint64_t id = ((uint64_t) i) * idgen->Ngrid * idgen->Ngrid + ((uint64_t)j) * idgen->Ngrid + k + 1;
    return id;
}

References IDGenerator::Ngrid, IDGenerator::offset, and IDGenerator::size.

Referenced by main(), and write_particle_data().

Here is the caller graph for this function:

idgen_create_pos_from_index()

void idgen_create_pos_from_index	(	IDGenerator *	idgen,
		int	index,
		double	pos[3]
	)

Definition at line 79 of file zeldovich.c.

{
    int x = index / (idgen->size[2] * idgen->size[1]) + idgen->offset[0];
    int y = (index % (idgen->size[1] * idgen->size[2])) / idgen->size[2] + idgen->offset[1];
    int z = (index % idgen->size[2]) + idgen->offset[2];
 
    pos[0] = x * idgen->BoxSize / idgen->Ngrid;
    pos[1] = y * idgen->BoxSize / idgen->Ngrid;
    pos[2] = z * idgen->BoxSize / idgen->Ngrid;
}

References IDGenerator::BoxSize, IDGenerator::Ngrid, IDGenerator::offset, and IDGenerator::size.

Referenced by setup_glass(), and setup_grid().

Here is the caller graph for this function:

idgen_init()

void idgen_init	(	IDGenerator *	idgen,
		PetaPM *	pm,
		int	Ngrid,
		double	BoxSize
	)

Definition at line 48 of file zeldovich.c.

{
 
    int * ThisTask2d = petapm_get_thistask2d(pm);
    int * NTask2d = petapm_get_ntask2d(pm);
    idgen->NumPart = 1;
    int k;
    for(k = 0; k < 2; k ++) {
        idgen->offset[k] = (ThisTask2d[k]) * Ngrid / NTask2d[k];
        idgen->size[k] = (ThisTask2d[k] + 1) * Ngrid / NTask2d[k];
        idgen->size[k] -= idgen->offset[k];
        idgen->NumPart *= idgen->size[k];
    }
    idgen->offset[2] = 0;
    idgen->size[2] = Ngrid;
    idgen->NumPart *= idgen->size[2];
    idgen->Ngrid = Ngrid;
    idgen->BoxSize = BoxSize;
}

References IDGenerator::BoxSize, IDGenerator::Ngrid, IDGenerator::NumPart, IDGenerator::offset, petapm_get_ntask2d(), petapm_get_thistask2d(), and IDGenerator::size.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

read_parameterfile()

void read_parameterfile	(	char *	fname,
		struct genic_config *	GenicConfig,
		int *	ShowBacktrace,
		double *	MaxMemSizePerNode,
		Cosmology *	CP
	)

Definition at line 75 of file params.c.

{
 
    /* read parameter file on all processes for simplicty */
    ParameterSet * ps = create_parameters();
    int ThisTask;
 
    if(0 != param_parse_file(ps, fname)) {
        endrun(0, "Parsing %s failed. \n", fname);
    }
    if(0 != param_validate(ps)) {
        endrun(0, "Validation of %s failed.\n", fname);
    }
 
    message(0, "----------- Running with Parameters ----------\n");
 
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0)
        param_dump(ps, stdout);
 
    message(0, "----------------------------------------------\n");
 
    /*Cosmology*/
    CP->Omega0 = param_get_double(ps, "Omega0");
    CP->OmegaLambda = param_get_double(ps, "OmegaLambda");
    CP->OmegaBaryon = param_get_double(ps, "OmegaBaryon");
    CP->HubbleParam = param_get_double(ps, "HubbleParam");
    CP->Omega_fld = param_get_double(ps, "Omega_fld");
    CP->w0_fld = param_get_double(ps,"w0_fld");
    CP->wa_fld = param_get_double(ps,"wa_fld");
    CP->Omega_ur = param_get_double(ps, "Omega_ur");
    if(CP->OmegaLambda > 0 && CP->Omega_fld > 0)
        endrun(0, "Cannot have OmegaLambda and Omega_fld (evolving dark energy) at the same time!\n");
    CP->CMBTemperature = param_get_double(ps, "CMBTemperature");
    CP->RadiationOn = param_get_double(ps, "RadiationOn");
    CP->MNu[0] = param_get_double(ps, "MNue");
    CP->MNu[1] = param_get_double(ps, "MNum");
    CP->MNu[2] = param_get_double(ps, "MNut");
    GenicConfig->WDM_therm_mass = param_get_double(ps, "MWDM_therm");
    *MaxMemSizePerNode = param_get_double(ps, "MaxMemSizePerNode");
    if(*MaxMemSizePerNode <= 1) {
        (*MaxMemSizePerNode) *= get_physmem_bytes() / (1024 * 1024);
    }
 
    GenicConfig->ProduceGas = param_get_int(ps, "ProduceGas");
    GenicConfig->InvertPhase = param_get_int(ps, "InvertPhase");
    /*Unit system*/
    GenicConfig->units.UnitVelocity_in_cm_per_s = param_get_double(ps, "UnitVelocity_in_cm_per_s");
    GenicConfig->units.UnitLength_in_cm = param_get_double(ps, "UnitLength_in_cm");
    GenicConfig->units.UnitMass_in_g = param_get_double(ps, "UnitMass_in_g");
 
 
    *ShowBacktrace = param_get_int(ps, "ShowBacktrace");
 
    double Redshift = param_get_double(ps, "Redshift");
 
    /*Parameters of the power spectrum*/
    GenicConfig->PowerP.InputPowerRedshift = param_get_double(ps, "InputPowerRedshift");
    if(GenicConfig->PowerP.InputPowerRedshift < 0)
        GenicConfig->PowerP.InputPowerRedshift = Redshift;
    GenicConfig->PowerP.Sigma8 = param_get_double(ps, "Sigma8");
    /*Always specify Sigm8 at z=0*/
    if(GenicConfig->PowerP.Sigma8 > 0)
        GenicConfig->PowerP.InputPowerRedshift = 0;
    GenicConfig->PowerP.FileWithInputSpectrum = param_get_string(ps, "FileWithInputSpectrum");
    GenicConfig->PowerP.FileWithTransferFunction = param_get_string(ps, "FileWithTransferFunction");
    GenicConfig->PowerP.DifferentTransferFunctions = param_get_int(ps, "DifferentTransferFunctions");
    GenicConfig->PowerP.ScaleDepVelocity = param_get_int(ps, "ScaleDepVelocity");
    /* By default ScaleDepVelocity follows DifferentTransferFunctions.*/
    if(GenicConfig->PowerP.ScaleDepVelocity < 0) {
        GenicConfig->PowerP.ScaleDepVelocity = GenicConfig->PowerP.DifferentTransferFunctions;
    }
    GenicConfig->PowerP.WhichSpectrum = param_get_int(ps, "WhichSpectrum");
    GenicConfig->PowerP.PrimordialIndex = param_get_double(ps, "PrimordialIndex");
    GenicConfig->PowerP.PrimordialRunning = param_get_double(ps, "PrimordialRunning");
 
    /*Simulation parameters*/
    GenicConfig->UsePeculiarVelocity = param_get_int(ps, "UsePeculiarVelocity");
    GenicConfig->SavePrePos = param_get_int(ps, "SavePrePos");
    GenicConfig->PrePosGridCenter = param_get_int(ps, "PrePosGridCenter");
    GenicConfig->BoxSize = param_get_double(ps, "BoxSize");
    GenicConfig->Nmesh = param_get_int(ps, "Nmesh");
    GenicConfig->Ngrid = param_get_int(ps, "Ngrid");
    GenicConfig->NgridGas = param_get_int(ps, "NgridGas");
    if(GenicConfig->NgridGas < 0)
        GenicConfig->NgridGas = GenicConfig->Ngrid;
    if(!GenicConfig->ProduceGas)
        GenicConfig->NgridGas = 0;
    /*Enable 'hybrid' neutrinos*/
    GenicConfig->NGridNu = param_get_int(ps, "NgridNu");
    /* Convert physical km/s at z=0 in an unperturbed universe to
     * internal gadget (comoving) velocity units at starting redshift.*/
    GenicConfig->Max_nuvel = param_get_double(ps, "Max_nuvel") * pow(1+Redshift, 1.5) * (GenicConfig->units.UnitVelocity_in_cm_per_s/1e5);
    GenicConfig->Seed = param_get_int(ps, "Seed");
    GenicConfig->UnitaryAmplitude = param_get_int(ps, "UnitaryAmplitude");
    param_get_string2(ps, "OutputDir", GenicConfig->OutputDir, sizeof(GenicConfig->OutputDir));
    param_get_string2(ps, "FileBase", GenicConfig->InitCondFile, sizeof(GenicConfig->InitCondFile));
    GenicConfig->MakeGlassGas = param_get_int(ps, "MakeGlassGas");
    /* We want to use a baryon glass by default if we have different transfer functions,
     * since that is the way we reproduce the linear growth. Otherwise use a grid by default.*/
    if(GenicConfig->MakeGlassGas < 0) {
        if(GenicConfig->PowerP.DifferentTransferFunctions)
            GenicConfig->MakeGlassGas = 1;
        else
            GenicConfig->MakeGlassGas = 0;
    }
    GenicConfig->MakeGlassCDM = param_get_int(ps, "MakeGlassCDM");
 
    int64_t NumPartPerFile = param_get_int(ps, "NumPartPerFile");
 
    int64_t Ngrid = GenicConfig->Ngrid;
    if(Ngrid < GenicConfig->NgridGas)
        Ngrid = GenicConfig->NgridGas;
    GenicConfig->NumFiles = ( Ngrid*Ngrid*Ngrid + NumPartPerFile - 1) / NumPartPerFile;
    GenicConfig->NumWriters = param_get_int(ps, "NumWriters");
    if(GenicConfig->PowerP.DifferentTransferFunctions && GenicConfig->PowerP.InputPowerRedshift != Redshift
        && (GenicConfig->ProduceGas || CP->MNu[0] + CP->MNu[1] + CP->MNu[2]))
        message(0, "WARNING: Using different transfer functions but also rescaling power to account for linear growth. NOT what you want!\n");
    if((CP->MNu[0] + CP->MNu[1] + CP->MNu[2] > 0) || GenicConfig->PowerP.DifferentTransferFunctions || GenicConfig->PowerP.ScaleDepVelocity)
        if(0 == strlen(GenicConfig->PowerP.FileWithTransferFunction))
            endrun(0,"For massive neutrinos, different transfer functions, or scale dependent growth functions you must specify a transfer function file\n");
    if(!CP->RadiationOn && (CP->MNu[0] + CP->MNu[1] + CP->MNu[2] > 0))
        endrun(0,"You want massive neutrinos but no background radiation: this will give an inconsistent cosmology.\n");
 
    if(GenicConfig->Nmesh == 0) {
        GenicConfig->Nmesh = 2*Ngrid;
    }
    /*Set some units*/
    GenicConfig->TimeIC = 1 / (1 + Redshift);
    double UnitTime_in_s = GenicConfig->units.UnitLength_in_cm / GenicConfig->units.UnitVelocity_in_cm_per_s;
 
    CP->Hubble = HUBBLE * UnitTime_in_s;
}

References genic_config::BoxSize, Cosmology::CMBTemperature, CP, create_parameters(), power_params::DifferentTransferFunctions, endrun(), power_params::FileWithInputSpectrum, power_params::FileWithTransferFunction, get_physmem_bytes(), Cosmology::Hubble, HUBBLE, Cosmology::HubbleParam, genic_config::InitCondFile, power_params::InputPowerRedshift, genic_config::InvertPhase, genic_config::MakeGlassCDM, genic_config::MakeGlassGas, genic_config::Max_nuvel, message(), Cosmology::MNu, genic_config::Ngrid, genic_config::NgridGas, genic_config::NGridNu, genic_config::Nmesh, genic_config::NumFiles, genic_config::NumWriters, Cosmology::Omega0, Cosmology::Omega_fld, Cosmology::Omega_ur, Cosmology::OmegaBaryon, Cosmology::OmegaLambda, genic_config::OutputDir, param_dump(), param_get_double(), param_get_int(), param_get_string(), param_get_string2(), param_parse_file(), param_validate(), genic_config::PowerP, genic_config::PrePosGridCenter, power_params::PrimordialIndex, power_params::PrimordialRunning, genic_config::ProduceGas, Cosmology::RadiationOn, genic_config::SavePrePos, power_params::ScaleDepVelocity, genic_config::Seed, ShowBacktrace, power_params::Sigma8, ThisTask, genic_config::TimeIC, genic_config::UnitaryAmplitude, UnitSystem::UnitLength_in_cm, UnitSystem::UnitMass_in_g, genic_config::units, UnitSystem::UnitVelocity_in_cm_per_s, genic_config::UsePeculiarVelocity, Cosmology::w0_fld, Cosmology::wa_fld, genic_config::WDM_therm_mass, and power_params::WhichSpectrum.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

saveheader()

void saveheader	(	BigFile *	bf,
		int64_t	TotNumPartCDM,
		int64_t	TotNumPartGas,
		int64_t	TotNuPart,
		double	nufrac,
		const double	BoxSize,
		Cosmology *	CP,
		const struct genic_config	GenicConfig
	)

Definition at line 112 of file save.c.

                                                                                                                                                                                           {
    BigBlock bheader;
    if(0 != big_file_mpi_create_block(bf, &bheader, "Header", NULL, 0, 0, 0, MPI_COMM_WORLD)) {
        endrun(0, "failed to create block %s:%s\n", "Header",
                big_file_get_error_message());
    }
 
    int64_t totnumpart[6] = {TotNumPartGas, TotNumPartCDM, TotNuPart, 0, 0, 0};
    double mass[6] = {0};
    compute_mass(mass, TotNumPartCDM, TotNumPartGas, TotNuPart, nufrac, BoxSize, CP, GenicConfig);
 
    double redshift = 1.0 / GenicConfig.TimeIC - 1.;
 
    int rt =(0 != big_block_set_attr(&bheader, "TotNumPart", totnumpart, "i8", 6)) ||
            (0 != big_block_set_attr(&bheader, "MassTable", mass, "f8", 6)) ||
            (big_block_set_attr(&bheader, "Time", &GenicConfig.TimeIC, "f8", 1)) ||
            (big_block_set_attr(&bheader, "Redshift", &redshift, "f8", 1)) ||
            (big_block_set_attr(&bheader, "BoxSize", &BoxSize, "f8", 1)) ||
            (big_block_set_attr(&bheader, "UsePeculiarVelocity", &GenicConfig.UsePeculiarVelocity, "i4", 1)) ||
            (big_block_set_attr(&bheader, "Omega0", &CP->Omega0, "f8", 1)) ||
            (big_block_set_attr(&bheader, "FractionNuInParticles", &nufrac, "f8", 1)) ||
            (big_block_set_attr(&bheader, "OmegaBaryon", &CP->OmegaBaryon, "f8", 1)) ||
            (big_block_set_attr(&bheader, "OmegaLambda", &CP->OmegaLambda, "f8", 1)) ||
            (big_block_set_attr(&bheader, "UnitLength_in_cm", &GenicConfig.units.UnitLength_in_cm, "f8", 1)) ||
            (big_block_set_attr(&bheader, "UnitMass_in_g", &GenicConfig.units.UnitMass_in_g, "f8", 1)) ||
            (big_block_set_attr(&bheader, "UnitVelocity_in_cm_per_s", &GenicConfig.units.UnitVelocity_in_cm_per_s, "f8", 1)) ||
            (big_block_set_attr(&bheader, "HubbleParam", &CP->HubbleParam, "f8", 1)) ||
            (big_block_set_attr(&bheader, "InvertPhase", &GenicConfig.InvertPhase, "i4", 1)) ||
            (big_block_set_attr(&bheader, "Seed", &GenicConfig.Seed, "i8", 1)) ||
            (big_block_set_attr(&bheader, "UnitaryAmplitude", &GenicConfig.UnitaryAmplitude, "i4", 1));
    if(rt) {
        endrun(0, "failed to create attr %s",
                big_file_get_error_message());
    }
 
    if(0 != big_block_mpi_close(&bheader, MPI_COMM_WORLD)) {
        endrun(0, "failed to close block %s:%s", "Header",
                    big_file_get_error_message());
    }
}

References compute_mass(), CP, endrun(), Cosmology::HubbleParam, genic_config::InvertPhase, Cosmology::Omega0, Cosmology::OmegaBaryon, Cosmology::OmegaLambda, genic_config::Seed, genic_config::TimeIC, genic_config::UnitaryAmplitude, UnitSystem::UnitLength_in_cm, UnitSystem::UnitMass_in_g, genic_config::units, UnitSystem::UnitVelocity_in_cm_per_s, and genic_config::UsePeculiarVelocity.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

setup_glass()

int setup_glass	(	IDGenerator *	idgen,
		PetaPM *	pm,
		double	shift,
		int	seed,
		double	mass,
		struct ic_part_data *	ICP,
		const double	UnitLength_in_cm,
		const char *	OutputDir
	)

Definition at line 41 of file glass.c.

{
    gsl_rng * rng = gsl_rng_alloc(gsl_rng_ranlxd1);
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    gsl_rng_set(rng, seed + ThisTask);
    memset(ICP, 0, idgen->NumPart*sizeof(struct ic_part_data));
 
    int i;
    /* Note: this loop should nto be omp because
     * of the call to gsl_rng_uniform*/
    for(i = 0; i < idgen->NumPart; i ++) {
        int k;
        idgen_create_pos_from_index(idgen, i, &ICP[i].Pos[0]);
        /* a spread of 3 will kill most of the grid anisotropy structure;
         * and still being local */
        for(k = 0; k < 3; k++) {
            double rand = idgen->BoxSize / idgen->Ngrid * 3 * (gsl_rng_uniform(rng) - 0.5);
            ICP[i].Pos[k] += shift + rand;
        }
        ICP[i].Mass = mass;
    }
 
    gsl_rng_free(rng);
 
    char * fn = fastpm_strdup_printf("powerspectrum-glass-%08X", seed);
    glass_evolve(pm, 14, fn, ICP, idgen->NumPart, UnitLength_in_cm, OutputDir);
    myfree(fn);
 
    return idgen->NumPart;
}

References IDGenerator::BoxSize, fastpm_strdup_printf(), glass_evolve(), idgen_create_pos_from_index(), ic_part_data::Mass, myfree, IDGenerator::Ngrid, IDGenerator::NumPart, ic_part_data::Pos, ThisTask, and UnitLength_in_cm.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

setup_grid()

int setup_grid	(	IDGenerator *	idgen,
		double	shift,
		double	mass,
		struct ic_part_data *	ICP
	)

Definition at line 91 of file zeldovich.c.

{
    memset(ICP, 0, idgen->NumPart*sizeof(struct ic_part_data));
 
    int i;
    #pragma omp parallel for
    for(i = 0; i < idgen->NumPart; i ++) {
        idgen_create_pos_from_index(idgen, i, &ICP[i].Pos[0]);
        ICP[i].Pos[0] += shift;
        ICP[i].Pos[1] += shift;
        ICP[i].Pos[2] +=  shift;
        ICP[i].Mass = mass;
    }
    return idgen->NumPart;
}

References idgen_create_pos_from_index(), ic_part_data::Mass, IDGenerator::NumPart, and ic_part_data::Pos.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

write_particle_data()

void write_particle_data	(	IDGenerator *	idgen,
		const int	Type,
		BigFile *	bf,
		const uint64_t	FirstID,
		const int	SavePrePos,
		int	NumFiles,
		int	NumWriters,
		struct ic_part_data *	curICP
	)

Definition at line 61 of file save.c.

{
    /* Write particles */
    if(SavePrePos)
        saveblock(bf, &curICP[0].PrePos, Type, "PrePosition", "f8", 3, idgen->NumPart, sizeof(curICP[0]), NumFiles, NumWriters);
    saveblock(bf, &curICP[0].Density, Type, "ICDensity", "f4", 1, idgen->NumPart, sizeof(curICP[0]), NumFiles, NumWriters);
    saveblock(bf, &curICP[0].Pos, Type, "Position", "f8", 3, idgen->NumPart, sizeof(curICP[0]), NumFiles, NumWriters);
    saveblock(bf, &curICP[0].Vel, Type, "Velocity", "f4", 3, idgen->NumPart, sizeof(curICP[0]), NumFiles, NumWriters);
    /*Generate and write IDs*/
    uint64_t * ids = (uint64_t *) mymalloc("IDs", idgen->NumPart * sizeof(uint64_t));
    memset(ids, 0, idgen->NumPart * sizeof(uint64_t));
    int i;
    #pragma omp parallel for
    for(i = 0; i < idgen->NumPart; i++)
    {
        ids[i] = idgen_create_id_from_index(idgen, i) + FirstID;
    }
    saveblock(bf, ids, Type, "ID", "u8", 1, idgen->NumPart, sizeof(uint64_t), NumFiles, NumWriters);
    myfree(ids);
    walltime_measure("/Write");
}

References curICP, idgen_create_id_from_index(), myfree, mymalloc, IDGenerator::NumPart, saveblock(), and walltime_measure.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function: