gravity.h File Reference

#include "forcetree.h"
#include "petapm.h"
#include "powerspectrum.h"
#include "timestep.h"

Include dependency graph for gravity.h:

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

Go to the source code of this file.

Classes
struct	gravshort_tree_params

Enumerations
enum	ShortRangeForceWindowType { SHORTRANGE_FORCE_WINDOW_TYPE_EXACT = 0 , SHORTRANGE_FORCE_WINDOW_TYPE_ERFC = 1 }

Functions
void	gravshort_fill_ntab (const enum ShortRangeForceWindowType ShortRangeForceWindowType, const double Asmth)
void	gravshort_set_softenings (double MeanDMSeparation)
double	FORCE_SOFTENING (int i, int type)
void	gravpm_init_periodic (PetaPM *pm, double BoxSize, double Asmth, int Nmesh, double G)
int	grav_apply_short_range_window (double r, double *fac, double *pot, const double cellsize)
void	set_gravshort_tree_params (ParameterSet *ps)
void	set_gravshort_treepar (struct gravshort_tree_params tree_params)
struct gravshort_tree_params	get_gravshort_treepar (void)
void	gravpm_force (PetaPM *pm, ForceTree *tree, Cosmology *CP, double Time, double UnitLength_in_cm, const char *PowerOutputDir, double TimeIC, int FastParticleType)
void	grav_short_pair (const ActiveParticles *act, PetaPM *pm, ForceTree *tree, double Rcut, double rho0, int NeutrinoTracer, int FastParticleType)
void	grav_short_tree (const ActiveParticles *act, PetaPM *pm, ForceTree *tree, double rho0, int NeutrinoTracer, int FastParticleType)
void	measure_power_spectrum (PetaPM *pm, int64_t k2, int kpos[3], pfft_complex *value)
void	powerspectrum_add_mode (Power *PowerSpectrum, const int64_t k2, const int kpos[3], pfft_complex *const value, const double invwindow, double Nmesh)

Enumeration Type Documentation

ShortRangeForceWindowType

enum ShortRangeForceWindowType

Enumerator
SHORTRANGE_FORCE_WINDOW_TYPE_EXACT
SHORTRANGE_FORCE_WINDOW_TYPE_ERFC

Definition at line 23 of file gravity.h.

                               {
    SHORTRANGE_FORCE_WINDOW_TYPE_EXACT = 0,
    SHORTRANGE_FORCE_WINDOW_TYPE_ERFC = 1,
};

Function Documentation

FORCE_SOFTENING()

double FORCE_SOFTENING	(	int	i,
		int	type
	)

Definition at line 36 of file gravshort-tree.c.

{
    if (TreeParams.AdaptiveSoftening == 1 && type == 0) {
        return P[i].Hsml;
    }
    /* Force is Newtonian beyond this.*/
    return 2.8 * GravitySoftening;
}

References gravshort_tree_params::AdaptiveSoftening, GravitySoftening, P, and TreeParams.

Referenced by blackhole_accretion_ngbiter(), density(), force_treeev_shortrange(), get_timestep_dloga(), grav_force(), grav_short_copy(), grav_short_pair_ngbiter(), and grav_short_postprocess().

Here is the caller graph for this function:

get_gravshort_treepar()

struct gravshort_tree_params get_gravshort_treepar

(

void

)

Definition at line 55 of file gravshort-tree.c.

{
    return TreeParams;
}

References TreeParams.

Referenced by run(), and run_gravity_test().

Here is the caller graph for this function:

grav_apply_short_range_window()

int grav_apply_short_range_window	(	double	r,
		double *	fac,
		double *	pot,
		const double	cellsize
	)

Definition at line 55 of file gravity.c.

{
    const double dx = shortrange_force_kernels[1][0];
    double i = (r / cellsize / dx);
    size_t tabindex = floor(i);
    if(tabindex >= NTAB - 1)
        return 1;
    /* use a linear interpolation; */
    *fac *= (tabindex + 1 - i) * shortrange_table[tabindex] + (i - tabindex) * shortrange_table[tabindex + 1];
    *pot *= (tabindex + 1 - i) * shortrange_table_potential[tabindex] + (i - tabindex) * shortrange_table_potential[tabindex];
    return 0;
}

References NTAB, shortrange_force_kernels, shortrange_table, and shortrange_table_potential.

Referenced by apply_accn_to_output(), and grav_short_pair_ngbiter().

Here is the caller graph for this function:

grav_short_pair()

void grav_short_pair	(	const ActiveParticles *	act,
		PetaPM *	pm,
		ForceTree *	tree,
		double	Rcut,
		double	rho0,
		int	NeutrinoTracer,
		int	FastParticleType
	)

Definition at line 24 of file gravshort-pair.c.

{
    TreeWalk tw[1] = {{0}};
 
    struct GravShortPriv priv;
    priv.cellsize = tree->BoxSize / pm->Nmesh;
    priv.Rcut = Rcut * pm->Asmth * priv.cellsize;
    priv.FastParticleType = FastParticleType;
    priv.NeutrinoTracer = NeutrinoTracer;
    priv.G = pm->G;
    priv.cbrtrho0 = pow(rho0, 1.0 / 3);
 
    message(0, "Starting pair-wise short range gravity...\n");
 
    tw->ev_label = "GRAV_SHORT";
    tw->visit = (TreeWalkVisitFunction) treewalk_visit_ngbiter;
    tw->ngbiter_type_elsize = sizeof(TreeWalkNgbIterGravShort);
    tw->ngbiter = (TreeWalkNgbIterFunction) grav_short_pair_ngbiter;
 
    tw->haswork = NULL;
    tw->fill = (TreeWalkFillQueryFunction) grav_short_copy;
    tw->reduce = (TreeWalkReduceResultFunction) grav_short_reduce;
    tw->postprocess = (TreeWalkProcessFunction) grav_short_postprocess;
    tw->query_type_elsize = sizeof(TreeWalkQueryGravShort);
    tw->result_type_elsize = sizeof(TreeWalkResultGravShort);
    tw->tree = tree;
    tw->priv = &priv;
 
    walltime_measure("/Misc");
 
    treewalk_run(tw, act->ActiveParticle, act->NumActiveParticle);
 
    walltime_measure("/Tree/Pairwise");
}

References ActiveParticles::ActiveParticle, PetaPM::Asmth, ForceTree::BoxSize, GravShortPriv::cbrtrho0, GravShortPriv::cellsize, TreeWalk::ev_label, GravShortPriv::FastParticleType, TreeWalk::fill, GravShortPriv::G, PetaPM::G, grav_short_copy(), grav_short_pair_ngbiter(), grav_short_postprocess(), grav_short_reduce(), TreeWalk::haswork, message(), GravShortPriv::NeutrinoTracer, TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, PetaPM::Nmesh, ActiveParticles::NumActiveParticle, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, GravShortPriv::Rcut, TreeWalk::reduce, TreeWalk::result_type_elsize, TreeWalk::tree, treewalk_run(), treewalk_visit_ngbiter(), TreeWalk::visit, and walltime_measure.

Referenced by run(), and run_gravity_test().

Here is the call graph for this function:

Here is the caller graph for this function:

grav_short_tree()

void grav_short_tree	(	const ActiveParticles *	act,
		PetaPM *	pm,
		ForceTree *	tree,
		double	rho0,
		int	NeutrinoTracer,
		int	FastParticleType
	)

This function computes the gravitational forces for all active particles. If needed, a new tree is constructed, otherwise the dynamically updated tree is used. Particles are only exported to other processors when really needed, thereby allowing a good use of the communication buffer. NeutrinoTracer = All.HybridNeutrinosOn && (atime <= All.HybridNuPartTime); rho0 = CP.Omega0 * 3 * CP.Hubble * CP.Hubble / (8 * M_PI * G)

Definition at line 105 of file gravshort-tree.c.

{
    double timeall = 0;
    double timetree, timewait, timecomm;
 
    TreeWalk tw[1] = {{0}};
    struct GravShortPriv priv;
    priv.cellsize = tree->BoxSize / pm->Nmesh;
    priv.Rcut = TreeParams.Rcut * pm->Asmth * priv.cellsize;;
    priv.ErrTolForceAcc = TreeParams.ErrTolForceAcc;
    priv.TreeUseBH = TreeParams.TreeUseBH;
    priv.BHOpeningAngle = TreeParams.BHOpeningAngle;
    priv.FastParticleType = FastParticleType;
    priv.NeutrinoTracer = NeutrinoTracer;
    priv.G = pm->G;
    priv.cbrtrho0 = pow(rho0, 1.0 / 3);
 
    if(!tree->moments_computed_flag)
        endrun(2, "Gravtree called before tree moments computed!\n");
 
    tw->ev_label = "GRAVTREE";
    tw->visit = (TreeWalkVisitFunction) force_treeev_shortrange;
    /* gravity applies to all particles. Including Tracer particles to enhance numerical stability. */
    tw->haswork = NULL;
    tw->reduce = (TreeWalkReduceResultFunction) grav_short_reduce;
    tw->postprocess = (TreeWalkProcessFunction) grav_short_postprocess;
 
    tw->query_type_elsize = sizeof(TreeWalkQueryGravShort);
    tw->result_type_elsize = sizeof(TreeWalkResultGravShort);
    tw->fill = (TreeWalkFillQueryFunction) grav_short_copy;
    tw->tree = tree;
    tw->priv = &priv;
 
    walltime_measure("/Misc");
 
    /* allocate buffers to arrange communication */
    MPIU_Barrier(MPI_COMM_WORLD);
    message(0, "Begin tree force.  (presently allocated=%g MB)\n", mymalloc_usedbytes() / (1024.0 * 1024.0));
 
    walltime_measure("/Misc");
 
    treewalk_run(tw, act->ActiveParticle, act->NumActiveParticle);
 
    /* Now the force computation is finished */
    /*  gather some diagnostic information */
 
    timetree = tw->timecomp1 + tw->timecomp2 + tw->timecomp3;
    timewait = tw->timewait1 + tw->timewait2;
    timecomm= tw->timecommsumm1 + tw->timecommsumm2;
 
    walltime_add("/Tree/Walk1", tw->timecomp1);
    walltime_add("/Tree/Walk2", tw->timecomp2);
    walltime_add("/Tree/PostProcess", tw->timecomp3);
    walltime_add("/Tree/Send", tw->timecommsumm1);
    walltime_add("/Tree/Recv", tw->timecommsumm2);
    walltime_add("/Tree/Wait1", tw->timewait1);
    walltime_add("/Tree/Wait2", tw->timewait2);
 
    timeall = walltime_measure(WALLTIME_IGNORE);
 
    walltime_add("/Tree/Misc", timeall - (timetree + timewait + timecomm));
 
    /* TreeUseBH > 1 means use the BH criterion on the initial timestep only,
     * avoiding the fully open O(N^2) case.*/
    if(TreeParams.TreeUseBH > 1)
        TreeParams.TreeUseBH = 0;
}

References ActiveParticles::ActiveParticle, PetaPM::Asmth, gravshort_tree_params::BHOpeningAngle, GravShortPriv::BHOpeningAngle, ForceTree::BoxSize, GravShortPriv::cbrtrho0, GravShortPriv::cellsize, endrun(), gravshort_tree_params::ErrTolForceAcc, GravShortPriv::ErrTolForceAcc, TreeWalk::ev_label, GravShortPriv::FastParticleType, TreeWalk::fill, force_treeev_shortrange(), GravShortPriv::G, PetaPM::G, grav_short_copy(), grav_short_postprocess(), grav_short_reduce(), TreeWalk::haswork, message(), ForceTree::moments_computed_flag, MPIU_Barrier, mymalloc_usedbytes, GravShortPriv::NeutrinoTracer, PetaPM::Nmesh, ActiveParticles::NumActiveParticle, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, gravshort_tree_params::Rcut, GravShortPriv::Rcut, TreeWalk::reduce, TreeWalk::result_type_elsize, TreeWalk::timecommsumm1, TreeWalk::timecommsumm2, TreeWalk::timecomp1, TreeWalk::timecomp2, TreeWalk::timecomp3, TreeWalk::timewait1, TreeWalk::timewait2, TreeWalk::tree, TreeParams, gravshort_tree_params::TreeUseBH, GravShortPriv::TreeUseBH, treewalk_run(), TreeWalk::visit, walltime_add, WALLTIME_IGNORE, and walltime_measure.

Referenced by do_force_test(), run(), and run_gravity_test().

Here is the call graph for this function:

Here is the caller graph for this function:

gravpm_force()

void gravpm_force	(	PetaPM *	pm,
		ForceTree *	tree,
		Cosmology *	CP,
		double	Time,
		double	UnitLength_in_cm,
		const char *	PowerOutputDir,
		double	TimeIC,
		int	FastParticleType
	)

Definition at line 62 of file gravpm.c.

                                                                                                                                                                    {
    PetaPMParticleStruct pstruct = {
        P,
        sizeof(P[0]),
        (char*) &P[0].Pos[0]  - (char*) P,
        (char*) &P[0].Mass  - (char*) P,
        /* Regions allocated inside _prepare*/
        NULL,
        /* By default all particles are active. For hybrid neutrinos set below.*/
        NULL,
        PartManager->NumPart,
    };
 
    PetaPMGlobalFunctions global_functions = {NULL, NULL, potential_transfer};
 
    /*Initialise the kspace neutrino code if it is enabled.
     * Mpc units are used to match power spectrum code.*/
    if(CP->MassiveNuLinRespOn) {
        global_functions.global_readout = measure_power_spectrum;
        global_functions.global_analysis = compute_neutrino_power;
    }
 
    if(CP->HybridNeutrinosOn && particle_nu_fraction(&(CP->ONu.hybnu), Time, 0) == 0.)
        pstruct.active = &hybrid_nu_gravpm_is_active;
 
    int i;
    #pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++)
    {
        P[i].GravPM[0] = P[i].GravPM[1] = P[i].GravPM[2] = 0;
    }
 
    /* Set up parameters*/
    GravPM.Time = Time;
    GravPM.TimeIC = TimeIC;
    GravPM.CP = CP;
    GravPM.FastParticleType = FastParticleType;
    GravPM.UnitLength_in_cm = UnitLength_in_cm;
    /*
     * we apply potential transfer immediately after the R2C transform,
     * Therefore the force transfer functions are based on the potential,
     * not the density.
     * */
    petapm_force(pm, _prepare, &global_functions, functions, &pstruct, tree);
    powerspectrum_sum(pm->ps);
    /*Now save the power spectrum*/
    powerspectrum_save(pm->ps, PowerOutputDir, "powerspectrum", Time, GrowthFactor(CP, Time, 1.0));
    /* Save the neutrino power if it is allocated*/
    if(pm->ps->logknu)
        powerspectrum_nu_save(pm->ps, PowerOutputDir, "powerspectrum-nu", Time);
    /*We are done with the power spectrum, free it*/
    powerspectrum_free(pm->ps);
    walltime_measure("/LongRange");
}

References _prepare(), PetaPMParticleStruct::active, compute_neutrino_power(), gravpm_params::CP, CP, gravpm_params::FastParticleType, functions, PetaPMGlobalFunctions::global_analysis, global_functions, PetaPMGlobalFunctions::global_readout, GravPM, GrowthFactor(), _omega_nu::hybnu, hybrid_nu_gravpm_is_active(), Cosmology::HybridNeutrinosOn, _powerspectrum::logknu, Cosmology::MassiveNuLinRespOn, measure_power_spectrum(), part_manager_type::NumPart, Cosmology::ONu, P, particle_nu_fraction(), PartManager, petapm_force(), potential_transfer(), powerspectrum_free(), powerspectrum_nu_save(), powerspectrum_save(), powerspectrum_sum(), PetaPM::ps, gravpm_params::Time, gravpm_params::TimeIC, gravpm_params::UnitLength_in_cm, UnitLength_in_cm, and walltime_measure.

Referenced by do_force_test(), run(), run_gravity_test(), and runpower().

Here is the call graph for this function:

Here is the caller graph for this function:

gravpm_init_periodic()

void gravpm_init_periodic	(	PetaPM *	pm,
		double	BoxSize,
		double	Asmth,
		int	Nmesh,
		double	G
	)

Definition at line 53 of file gravpm.c.

                                                                                     {
    petapm_init(pm, BoxSize, Asmth, Nmesh, G, MPI_COMM_WORLD);
}

References G, and petapm_init().

Referenced by do_force_test(), run(), run_gravity_test(), runfof(), and runpower().

Here is the call graph for this function:

Here is the caller graph for this function:

gravshort_fill_ntab()

void gravshort_fill_ntab	(	const enum ShortRangeForceWindowType	ShortRangeForceWindowType,
		const double	Asmth
	)

Definition at line 23 of file gravity.c.

{
    if (ShortRangeForceWindowType == SHORTRANGE_FORCE_WINDOW_TYPE_EXACT) {
        if(Asmth != 1.5) {
            endrun(0, "The short range force window is calibrated for Asmth = 1.5, but running with %g\n", Asmth);
        }
    }
 
    size_t i;
    for(i = 0; i < NTAB; i++)
    {
        /* force_kernels is in units of mesh points; */
        double u = shortrange_force_kernels[i][0] * 0.5 / Asmth;
        switch (ShortRangeForceWindowType) {
            case SHORTRANGE_FORCE_WINDOW_TYPE_EXACT:
                /* Notice that the table is only calibrated for smth of 1.25*/
                shortrange_table[i] = shortrange_force_kernels[i][2]; /* ~ erfc(u) + 2.0 * u / sqrt(M_PI) * exp(-u * u); */
                /* The potential of the calibrated kernel is a bit off, so we still use erfc here; we do not use potential anyways.*/
                shortrange_table_potential[i] = shortrange_force_kernels[i][1];
            break;
            case SHORTRANGE_FORCE_WINDOW_TYPE_ERFC:
                shortrange_table[i] = erfc(u) + 2.0 * u / sqrt(M_PI) * exp(-u * u);
                shortrange_table_potential[i] = erfc(u);
            break;
        }
        /* we don't have a table for that and don't use it anyways. */
        shortrange_table_tidal[i] = 4.0 * u * u * u / sqrt(M_PI) * exp(-u * u);
    }
}

References endrun(), NTAB, shortrange_force_kernels, SHORTRANGE_FORCE_WINDOW_TYPE_ERFC, SHORTRANGE_FORCE_WINDOW_TYPE_EXACT, shortrange_table, shortrange_table_potential, and shortrange_table_tidal.

Referenced by begrun(), and do_force_test().

Here is the call graph for this function:

Here is the caller graph for this function:

gravshort_set_softenings()

void gravshort_set_softenings

(

double

MeanSeparation

)

Sets the (comoving) softening length, converting from units of the mean DM separation to comoving internal units.

Sets the (comoving) softening length, converting from units of the mean separation to comoving internal units.

Definition at line 47 of file gravshort-tree.c.

{
    GravitySoftening = TreeParams.FractionalGravitySoftening * MeanSeparation;
    /* 0: Gas is collisional */
    message(0, "GravitySoftening = %g\n", GravitySoftening);
}

References gravshort_tree_params::FractionalGravitySoftening, GravitySoftening, message(), and TreeParams.

Referenced by do_force_test(), init(), and setup_density().

Here is the call graph for this function:

Here is the caller graph for this function:

measure_power_spectrum()

void measure_power_spectrum	(	PetaPM *	pm,
		int64_t	k2,
		int	kpos[3],
		pfft_complex *	value
	)

Definition at line 361 of file gravpm.c.

                                                                                  {
    double f = 1.0;
    /* the CIC deconvolution kernel is
     *
     * sinc_unnormed(k_x L / 2 Nmesh) ** 2
     *
     * k_x = kpos * 2pi / L
     *
     * */
    int k;
    for(k = 0; k < 3; k ++) {
        double tmp = (kpos[k] * M_PI) / pm->Nmesh;
        tmp = sinc_unnormed(tmp);
        f *= 1. / (tmp * tmp);
    }
    powerspectrum_add_mode(pm->ps, k2, kpos, value, f, pm->Nmesh);
}

References NODE::f, PetaPM::Nmesh, powerspectrum_add_mode(), PetaPM::ps, and sinc_unnormed().

Referenced by gravpm_force().

Here is the call graph for this function:

Here is the caller graph for this function:

powerspectrum_add_mode()

void powerspectrum_add_mode	(	Power *	PowerSpectrum,
		const int64_t	k2,
		const int	kpos[3],
		pfft_complex *const	value,
		const double	invwindow,
		double	Nmesh
	)

Definition at line 327 of file gravpm.c.

{
    if(k2 == 0) {
        /* Save zero mode corresponding to the mean as the normalisation factor.*/
        PowerSpectrum->Norm = (value[0][0] * value[0][0] + value[0][1] * value[0][1]);
        return;
    }
    /* Measure power spectrum: we don't want the zero mode.
     * Some modes with k_z = 0 or N/2 have weight 1, the rest have weight 2.
     * This is because of the symmetry of the real fft. */
    if(k2 > 0) {
        /*How many bins per unit (log) interval in k?*/
        const double binsperunit=(PowerSpectrum->size-1)/log(sqrt(3) * Nmesh/2.0);
        int kint=floor(binsperunit*log(k2)/2.);
        int w;
        const double keff = sqrt(kpos[0]*kpos[0]+kpos[1]*kpos[1]+kpos[2]*kpos[2]);
        const double m = (value[0][0] * value[0][0] + value[0][1] * value[0][1]);
        /*Make sure we do not overflow (although this should never happen)*/
        if(kint >= PowerSpectrum->size)
            return;
        if(kpos[2] == 0 || kpos[2] == Nmesh/2) w = 1;
        else w = 2;
        /*Make sure we use thread-local memory to avoid racing.*/
        const int index = kint + omp_get_thread_num() * PowerSpectrum->size;
        /*Multiply P(k) by inverse window function*/
        PowerSpectrum->Power[index] += w * m * invwindow * invwindow;
        PowerSpectrum->Nmodes[index] += w;
        PowerSpectrum->kk[index] += w * keff;
    }
 
}

References _powerspectrum::kk, _powerspectrum::Nmodes, _powerspectrum::Norm, _powerspectrum::Power, and _powerspectrum::size.

Referenced by measure_power_spectrum(), and potential_transfer().

Here is the caller graph for this function:

set_gravshort_tree_params()

void set_gravshort_tree_params

(

ParameterSet *

ps

)

Definition at line 67 of file gravshort-tree.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        TreeParams.BHOpeningAngle = param_get_double(ps, "BHOpeningAngle");
        TreeParams.ErrTolForceAcc = param_get_double(ps, "ErrTolForceAcc");
        TreeParams.BHOpeningAngle = param_get_double(ps, "BHOpeningAngle");
        TreeParams.TreeUseBH= param_get_int(ps, "TreeUseBH");
        TreeParams.Rcut = param_get_double(ps, "TreeRcut");
        TreeParams.FractionalGravitySoftening = param_get_double(ps, "GravitySoftening");
        TreeParams.AdaptiveSoftening = !param_get_int(ps, "GravitySofteningGas");
 
 
    }
    MPI_Bcast(&TreeParams, sizeof(struct gravshort_tree_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References gravshort_tree_params::AdaptiveSoftening, gravshort_tree_params::BHOpeningAngle, gravshort_tree_params::ErrTolForceAcc, gravshort_tree_params::FractionalGravitySoftening, param_get_double(), param_get_int(), gravshort_tree_params::Rcut, ThisTask, TreeParams, and gravshort_tree_params::TreeUseBH.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

set_gravshort_treepar()

void set_gravshort_treepar

(

struct gravshort_tree_params

tree_params

)

Definition at line 55 of file gravshort-tree.c.

{
    TreeParams = tree_params;
}

Referenced by do_force_test(), run_gravity_test(), and setup_density().

Here is the caller graph for this function: