density.h File Reference

#include "forcetree.h"
#include "timestep.h"
#include "densitykernel.h"
#include "utils/paramset.h"

Include dependency graph for density.h:

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

Go to the source code of this file.

Classes
struct	density_params
struct	sph_pred_data

Functions
void	set_density_params (ParameterSet *ps)
void	set_densitypar (struct density_params dp)
void	density (const ActiveParticles *act, int update_hsml, int DoEgyDensity, int BlackHoleOn, double MinEgySpec, const DriftKickTimes times, Cosmology *CP, struct sph_pred_data *SPH_predicted, MyFloat *GradRho, const ForceTree *const tree)
double	GetNumNgb (enum DensityKernelType KernelType)
enum DensityKernelType	GetDensityKernelType (void)
struct sph_pred_data	slots_allocate_sph_pred_data (int nsph)
void	slots_free_sph_pred_data (struct sph_pred_data *sph_pred)
MyFloat	SPH_EntVarPred (int PI, double MinEgySpec, double a3inv, double dloga)
void	SPH_VelPred (int i, MyFloat *VelPred, const double FgravkickB, double gravkick, double hydrokick)

Function Documentation

density()

void density	(	const ActiveParticles *	act,
		int	update_hsml,
		int	DoEgyDensity,
		int	BlackHoleOn,
		double	MinEgySpec,
		const DriftKickTimes	times,
		Cosmology *	CP,
		struct sph_pred_data *	SPH_predicted,
		MyFloat *	GradRho,
		const ForceTree *const	tree
	)

This function computes the local density for each active SPH particle, the number of neighbours in the current smoothing radius, and the divergence and rotation of the velocity field. The pressure is updated as well. If a particle with its smoothing region is fully inside the local domain, it is not exported to the other processors. The function also detects particles that have a number of neighbours outside the allowed tolerance range. For these particles, the smoothing length is adjusted accordingly, and the density() computation is called again. Note that the smoothing length is not allowed to fall below the lower bound set by MinGasHsml (this may mean that one has to deal with substantially more than normal number of neighbours.)

Definition at line 203 of file density.c.

{
    TreeWalk tw[1] = {{0}};
    struct DensityPriv priv[1];
 
    tw->ev_label = "DENSITY";
    tw->visit = (TreeWalkVisitFunction) treewalk_visit_nolist_ngbiter;
    tw->NoNgblist = 1;
    tw->ngbiter_type_elsize = sizeof(TreeWalkNgbIterDensity);
    tw->ngbiter = (TreeWalkNgbIterFunction) density_ngbiter;
    tw->haswork = density_haswork;
    tw->fill = (TreeWalkFillQueryFunction) density_copy;
    tw->reduce = (TreeWalkReduceResultFunction) density_reduce;
    tw->postprocess = (TreeWalkProcessFunction) density_postprocess;
    tw->query_type_elsize = sizeof(TreeWalkQueryDensity);
    tw->result_type_elsize = sizeof(TreeWalkResultDensity);
    tw->priv = priv;
    tw->tree = tree;
 
    int i;
 
    double timeall = 0;
    double timecomp, timecomm, timewait;
    const double atime = exp(loga_from_ti(times.Ti_Current));
 
    walltime_measure("/Misc");
 
    DENSITY_GET_PRIV(tw)->Left = (MyFloat *) mymalloc("DENS_PRIV->Left", PartManager->NumPart * sizeof(MyFloat));
    DENSITY_GET_PRIV(tw)->Right = (MyFloat *) mymalloc("DENS_PRIV->Right", PartManager->NumPart * sizeof(MyFloat));
    DENSITY_GET_PRIV(tw)->NumNgb = (MyFloat *) mymalloc("DENS_PRIV->NumNgb", PartManager->NumPart * sizeof(MyFloat));
    DENSITY_GET_PRIV(tw)->Rot = (MyFloat (*) [3]) mymalloc("DENS_PRIV->Rot", SlotsManager->info[0].size * sizeof(priv->Rot[0]));
    /* This one stores the gradient for h finding. The factor stored in SPHP->DhsmlEgyDensityFactor depends on whether PE SPH is enabled.*/
    DENSITY_GET_PRIV(tw)->DhsmlDensityFactor = (MyFloat *) mymalloc("DENSITY_GET_PRIV(tw)->DhsmlDensity", PartManager->NumPart * sizeof(MyFloat));
 
    DENSITY_GET_PRIV(tw)->update_hsml = update_hsml;
    DENSITY_GET_PRIV(tw)->DoEgyDensity = DoEgyDensity;
 
    DENSITY_GET_PRIV(tw)->DesNumNgb = GetNumNgb(DensityParams.DensityKernelType);
    DENSITY_GET_PRIV(tw)->MinGasHsml = DensityParams.MinGasHsmlFractional * (FORCE_SOFTENING(1, 1)/2.8);
 
    DENSITY_GET_PRIV(tw)->BlackHoleOn = BlackHoleOn;
    DENSITY_GET_PRIV(tw)->SPH_predicted = SPH_predicted;
    DENSITY_GET_PRIV(tw)->GradRho = GradRho;
 
    /* Init Left and Right: this has to be done before treewalk */
    #pragma omp parallel for
    for(i = 0; i < act->NumActiveParticle; i++)  {
        int p_i = act->ActiveParticle ? act->ActiveParticle[i] : i;
        /* We only really need active particles with work
         * but I don't want to read the particle table here*/
        DENSITY_GET_PRIV(tw)->Right[p_i] = tree->BoxSize;
        DENSITY_GET_PRIV(tw)->NumNgb[p_i] = 0;
        DENSITY_GET_PRIV(tw)->Left[p_i] = 0;
    }
 
    priv->a3inv = pow(atime, -3);
    priv->MinEgySpec = MinEgySpec;
    /* Factor this out since all particles have the same drift time*/
    priv->FgravkickB = get_exact_gravkick_factor(CP, times.PM_kick, times.Ti_Current);
    memset(priv->gravkicks, 0, sizeof(priv->gravkicks[0])*(TIMEBINS+1));
    memset(priv->hydrokicks, 0, sizeof(priv->hydrokicks[0])*(TIMEBINS+1));
    /* Compute the factors to move a current kick times velocity to the drift time velocity.
     * We need to do the computation for all timebins up to the maximum because even inactive
     * particles may have interactions. */
    #pragma omp parallel for
    for(i = times.mintimebin; i <= TIMEBINS; i++)
    {
        priv->gravkicks[i] = get_exact_gravkick_factor(CP, times.Ti_kick[i], times.Ti_Current);
        priv->hydrokicks[i] = get_exact_hydrokick_factor(CP, times.Ti_kick[i], times.Ti_Current);
    }
    priv->times = &times;
 
    #pragma omp parallel for
    for(i = 0; i < act->NumActiveParticle; i++)
    {
        int p_i = act->ActiveParticle ? act->ActiveParticle[i] : i;
        if(P[p_i].Type == 0 && !P[p_i].IsGarbage) {
            int bin = P[p_i].TimeBin;
            double dloga = dloga_from_dti(priv->times->Ti_Current - priv->times->Ti_kick[bin], priv->times->Ti_Current);
            priv->SPH_predicted->EntVarPred[P[p_i].PI] = SPH_EntVarPred(P[p_i].PI, priv->MinEgySpec, priv->a3inv, dloga);
            SPH_VelPred(p_i, priv->SPH_predicted->VelPred + 3 * P[p_i].PI, priv->FgravkickB, priv->gravkicks[bin], priv->hydrokicks[bin]);
        }
    }
 
    /* allocate buffers to arrange communication */
 
    walltime_measure("/SPH/Density/Init");
 
    /* Do the treewalk with looping for hsml*/
    treewalk_do_hsml_loop(tw, act->ActiveParticle, act->NumActiveParticle, update_hsml);
 
    myfree(DENSITY_GET_PRIV(tw)->DhsmlDensityFactor);
    myfree(DENSITY_GET_PRIV(tw)->Rot);
    myfree(DENSITY_GET_PRIV(tw)->NumNgb);
    myfree(DENSITY_GET_PRIV(tw)->Right);
    myfree(DENSITY_GET_PRIV(tw)->Left);
 
 
    /* collect some timing information */
 
    timeall = walltime_measure(WALLTIME_IGNORE);
 
    timecomp = tw->timecomp3 + tw->timecomp1 + tw->timecomp2;
    timewait = tw->timewait1 + tw->timewait2;
    timecomm = tw->timecommsumm1 + tw->timecommsumm2;
 
    walltime_add("/SPH/Density/Compute", timecomp);
    walltime_add("/SPH/Density/Wait", timewait);
    walltime_add("/SPH/Density/Comm", timecomm);
    walltime_add("/SPH/Density/Misc", timeall - (timecomp + timewait + timecomm));
}

References DensityPriv::a3inv, ActiveParticles::ActiveParticle, DensityPriv::BlackHoleOn, ForceTree::BoxSize, CP, density_copy(), DENSITY_GET_PRIV, density_haswork(), density_ngbiter(), density_postprocess(), density_reduce(), density_params::DensityKernelType, DensityParams, DensityPriv::DhsmlDensityFactor, dloga, dloga_from_dti(), DensityPriv::DoEgyDensity, sph_pred_data::EntVarPred, TreeWalk::ev_label, DensityPriv::FgravkickB, TreeWalk::fill, FORCE_SOFTENING(), get_exact_gravkick_factor(), get_exact_hydrokick_factor(), GetNumNgb(), DensityPriv::GradRho, DensityPriv::gravkicks, TreeWalk::haswork, DensityPriv::hydrokicks, slots_manager_type::info, DensityPriv::Left, loga_from_ti(), DensityPriv::MinEgySpec, density_params::MinGasHsmlFractional, DriftKickTimes::mintimebin, myfree, mymalloc, TreeWalk::ngbiter, TreeWalk::ngbiter_type_elsize, TreeWalk::NoNgblist, ActiveParticles::NumActiveParticle, DensityPriv::NumNgb, part_manager_type::NumPart, P, PartManager, DriftKickTimes::PM_kick, TreeWalk::postprocess, TreeWalk::priv, TreeWalk::query_type_elsize, TreeWalk::reduce, TreeWalk::result_type_elsize, DensityPriv::Right, DensityPriv::Rot, slot_info::size, SlotsManager, SPH_EntVarPred(), DensityPriv::SPH_predicted, SPH_VelPred(), DriftKickTimes::Ti_Current, DriftKickTimes::Ti_kick, TIMEBINS, TreeWalk::timecommsumm1, TreeWalk::timecommsumm2, TreeWalk::timecomp1, TreeWalk::timecomp2, TreeWalk::timecomp3, DensityPriv::times, TreeWalk::timewait1, TreeWalk::timewait2, TreeWalk::tree, treewalk_do_hsml_loop(), treewalk_visit_nolist_ngbiter(), DensityPriv::update_hsml, sph_pred_data::VelPred, TreeWalk::visit, walltime_add, WALLTIME_IGNORE, and walltime_measure.

Referenced by density_postprocess(), do_density_test(), get_compton_cooling(), get_equilib_ne(), get_heatingcooling_rate(), get_helium_ion_phys_cgs(), get_individual_cooling(), get_ne_by_nh(), get_neutral_fraction_phys_cgs(), get_temp(), run(), runfof(), setup_density_indep_entropy(), and setup_smoothinglengths().

Here is the call graph for this function:

Here is the caller graph for this function:

GetDensityKernelType()

enum DensityKernelType GetDensityKernelType

(

void

)

Definition at line 55 of file density.c.

{
    return DensityParams.DensityKernelType;
}

References density_kernel_desnumngb(), density_kernel_init(), DensityParams, density_params::DensityResolutionEta, and kernel().

Referenced by blackhole_accretion_ngbiter(), blackhole_dynfric_ngbiter(), blackhole_feedback_ngbiter(), do_density_test(), hydro_ngbiter(), metal_return_ngbiter(), petaio_write_header(), set_density_params(), set_init_hsml(), setup_smoothinglengths(), stellar_density(), and stellar_density_ngbiter().

Here is the call graph for this function:

Here is the caller graph for this function:

GetNumNgb()

double GetNumNgb

(

enum DensityKernelType

KernelType

)

Definition at line 55 of file density.c.

{
    DensityKernel kernel;
    density_kernel_init(&kernel, 1.0, KernelType);
    return density_kernel_desnumngb(&kernel, DensityParams.DensityResolutionEta);
}

Referenced by density(), do_density_test(), set_density_params(), set_init_hsml(), setup_smoothinglengths(), and stellar_density().

Here is the caller graph for this function:

set_density_params()

void set_density_params

(

ParameterSet *

ps

)

Definition at line 32 of file density.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        DensityParams.DensityKernelType = (enum DensityKernelType) param_get_enum(ps, "DensityKernelType");
        DensityParams.MaxNumNgbDeviation = param_get_double(ps, "MaxNumNgbDeviation");
        DensityParams.DensityResolutionEta = param_get_double(ps, "DensityResolutionEta");
        DensityParams.MinGasHsmlFractional = param_get_double(ps, "MinGasHsmlFractional");
 
        DensityKernel kernel;
        density_kernel_init(&kernel, 1.0, DensityParams.DensityKernelType);
        message(1, "The Density Kernel type is %s\n", kernel.name);
        message(1, "The Density resolution is %g * mean separation, or %g neighbours\n",
                    DensityParams.DensityResolutionEta, GetNumNgb(GetDensityKernelType()));
        /*These two look like black hole parameters but they are really neighbour finding parameters*/
        DensityParams.BlackHoleNgbFactor = param_get_double(ps, "BlackHoleNgbFactor");
        DensityParams.BlackHoleMaxAccretionRadius = param_get_double(ps, "BlackHoleMaxAccretionRadius");
    }
    MPI_Bcast(&DensityParams, sizeof(struct density_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References density_params::BlackHoleMaxAccretionRadius, density_params::BlackHoleNgbFactor, density_kernel_init(), density_params::DensityKernelType, DensityParams, density_params::DensityResolutionEta, GetDensityKernelType(), GetNumNgb(), kernel(), density_params::MaxNumNgbDeviation, message(), density_params::MinGasHsmlFractional, param_get_double(), param_get_enum(), and ThisTask.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

set_densitypar()

void set_densitypar

(

struct density_params

dp

)

Definition at line 25 of file density.c.

{
    DensityParams = dp;
}

References DensityParams.

Referenced by do_density_test(), and setup_density().

Here is the caller graph for this function:

slots_allocate_sph_pred_data()

struct sph_pred_data slots_allocate_sph_pred_data

(

int

nsph

)

Definition at line 554 of file density.c.

{
    struct sph_pred_data sph_scratch;
    /*Data is allocated high so that we can free the tree around it*/
    sph_scratch.EntVarPred = (MyFloat *) mymalloc2("EntVarPred", sizeof(MyFloat) * nsph);
    memset(sph_scratch.EntVarPred, 0, sizeof(sph_scratch.EntVarPred[0]) * nsph);
    sph_scratch.VelPred = (MyFloat *) mymalloc2("VelPred", sizeof(MyFloat) * 3 * nsph);
    return sph_scratch;
}

Referenced by run(), runfof(), setup_density(), and setup_smoothinglengths().

Here is the caller graph for this function:

slots_free_sph_pred_data()

void slots_free_sph_pred_data

(

struct sph_pred_data *

sph_pred

)

Definition at line 665 of file density.c.

{
    myfree(sph_scratch->VelPred);
    sph_scratch->VelPred = NULL;
    myfree(sph_scratch->EntVarPred);
    sph_scratch->EntVarPred = NULL;
}

References sph_pred_data::EntVarPred, myfree, and sph_pred_data::VelPred.

Referenced by run(), runfof(), and setup_smoothinglengths().

Here is the caller graph for this function:

SPH_EntVarPred()

MyFloat SPH_EntVarPred	(	int	PI,
		double	MinEgySpec,
		double	a3inv,
		double	dloga
	)

Definition at line 71 of file density.c.

{
        double EntVarPred = SphP[PI].Entropy + SphP[PI].DtEntropy * dloga;
        /*Entropy limiter for the predicted entropy: makes sure entropy stays positive. */
        if(dloga > 0 && EntVarPred < 0.5*SphP[PI].Entropy)
            EntVarPred = 0.5 * SphP[PI].Entropy;
        const double enttou = pow(SphP[PI].Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
        if(EntVarPred < MinEgySpec / enttou)
            EntVarPred = MinEgySpec / enttou;
        EntVarPred = pow(EntVarPred, 1/GAMMA);
        return EntVarPred;
}

References dloga, GAMMA, GAMMA_MINUS1, and SphP.

Referenced by density(), density_ngbiter(), and hydro_ngbiter().

Here is the caller graph for this function:

SPH_VelPred()

void SPH_VelPred	(	int	i,
		MyFloat *	VelPred,
		const double	FgravkickB,
		double	gravkick,
		double	hydrokick
	)

Definition at line 89 of file density.c.

{
    int j;
    for(j = 0; j < 3; j++) {
        VelPred[j] = P[i].Vel[j] + gravkick * P[i].GravAccel[j]
            + P[i].GravPM[j] * FgravkickB + hydrokick * SPHP(i).HydroAccel[j];
    }
}

References P, and SPHP.

Referenced by density(), density_ngbiter(), and hydro_ngbiter().

Here is the caller graph for this function: