density.c File Reference

SPH density computation and smoothing length determination. More...

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <omp.h>
#include "physconst.h"
#include "walltime.h"
#include "cooling.h"
#include "density.h"
#include "treewalk.h"
#include "timefac.h"
#include "slotsmanager.h"
#include "timestep.h"
#include "utils.h"
#include "gravity.h"
#include "winds.h"

Include dependency graph for density.c:

Go to the source code of this file.

Classes
struct	TreeWalkNgbIterDensity
struct	TreeWalkQueryDensity
struct	TreeWalkResultDensity
struct	DensityPriv

Macros
#define	DENSITY_GET_PRIV(tw)   ((struct DensityPriv*) ((tw)->priv))

Functions
void	set_densitypar (struct density_params dp)
void	set_density_params (ParameterSet *ps)
double	GetNumNgb (enum DensityKernelType KernelType)
enum DensityKernelType	GetDensityKernelType (void)
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)
static void	density_ngbiter (TreeWalkQueryDensity *I, TreeWalkResultDensity *O, TreeWalkNgbIterDensity *iter, LocalTreeWalk *lv)
static int	density_haswork (int n, TreeWalk *tw)
static void	density_postprocess (int i, TreeWalk *tw)
static void	density_check_neighbours (int i, TreeWalk *tw)
static void	density_reduce (int place, TreeWalkResultDensity *remote, enum TreeWalkReduceMode mode, TreeWalk *tw)
static void	density_copy (int place, TreeWalkQueryDensity *I, TreeWalk *tw)
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)
struct sph_pred_data	slots_allocate_sph_pred_data (int nsph)
void	slots_free_sph_pred_data (struct sph_pred_data *sph_scratch)

Variables
static struct density_params	DensityParams

Detailed Description

SPH density computation and smoothing length determination.

This file contains the "first SPH loop", where the SPH densities and some auxiliary quantities are computed. There is also functionality that corrects the smoothing length if needed.

Definition in file density.c.

Macro Definition Documentation

DENSITY_GET_PRIV

#define DENSITY_GET_PRIV

(

tw

)

((struct DensityPriv*) ((tw)->priv))

Definition at line 166 of file density.c.

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:

density_check_neighbours()

void density_check_neighbours ( int  i, TreeWalk *  tw  )

static

Definition at line 554 of file density.c.

{
    /* now check whether we had enough neighbours */
    int tid = omp_get_thread_num();
    double desnumngb = DENSITY_GET_PRIV(tw)->DesNumNgb;
 
    if(DENSITY_GET_PRIV(tw)->BlackHoleOn && P[i].Type == 5)
        desnumngb = desnumngb * DensityParams.BlackHoleNgbFactor;
 
    MyFloat * Left = DENSITY_GET_PRIV(tw)->Left;
    MyFloat * Right = DENSITY_GET_PRIV(tw)->Right;
    MyFloat * NumNgb = DENSITY_GET_PRIV(tw)->NumNgb;
 
    if(NumNgb[i] < (desnumngb - DensityParams.MaxNumNgbDeviation) ||
            (NumNgb[i] > (desnumngb + DensityParams.MaxNumNgbDeviation)))
    {
        /* This condition is here to prevent the density code looping forever if it encounters
         * multiple particles at the same position. If this happens you likely have worse
         * problems anyway, so warn also. */
        if((Right[i] - Left[i]) < 1.0e-3 * Left[i])
        {
            /* If this happens probably the exchange is screwed up and all your particles have moved to (0,0,0)*/
            message(1, "Very tight Hsml bounds for i=%d ID=%lu Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g pos=(%g|%g|%g)\n",
             i, P[i].ID, P[i].Hsml, Left[i], Right[i], NumNgb[i], Right[i] - Left[i], P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
            P[i].Hsml = Right[i];
            return;
        }
 
        /* If we need more neighbours, move the lower bound up. If we need fewer, move the upper bound down.*/
        if(NumNgb[i] < desnumngb) {
                Left[i] = P[i].Hsml;
        } else {
                Right[i] = P[i].Hsml;
        }
 
        /* Next step is geometric mean of previous. */
        if((Right[i] < tw->tree->BoxSize && Left[i] > 0) || (P[i].Hsml * 1.26 > 0.99 * tw->tree->BoxSize))
            P[i].Hsml = pow(0.5 * (pow(Left[i], 3) + pow(Right[i], 3)), 1.0 / 3);
        else
        {
            if(!(Right[i] < tw->tree->BoxSize) && Left[i] == 0)
                endrun(8188, "Cannot occur. Check for memory corruption: i=%d L = %g R = %g N=%g. Type %d, Pos %g %g %g hsml %g Box %g\n", i, Left[i], Right[i], NumNgb[i], P[i].Type, P[i].Pos[0], P[i].Pos[1], P[i].Pos[2], P[i].Hsml, tw->tree->BoxSize);
 
            MyFloat DensFac = DENSITY_GET_PRIV(tw)->DhsmlDensityFactor[i];
            double fac = 1.26;
            if(NumNgb[i] > 0)
                fac = 1 - (NumNgb[i] - desnumngb) / (NUMDIMS * NumNgb[i]) * DensFac;
 
            /* Find the initial bracket using the kernel gradients*/
            if(Right[i] > 0.99 * tw->tree->BoxSize && Left[i] > 0)
                if(DensFac <= 0 || fabs(NumNgb[i] - desnumngb) >= 0.5 * desnumngb || fac > 1.26)
                    fac = 1.26;
 
            if(Right[i] < 0.99*tw->tree->BoxSize && Left[i] == 0)
                if(DensFac <=0 || fac < 1./3)
                    fac = 1./3;
 
            P[i].Hsml *= fac;
        }
 
        if(DENSITY_GET_PRIV(tw)->BlackHoleOn && P[i].Type == 5)
            if(Left[i] > DensityParams.BlackHoleMaxAccretionRadius)
            {
                P[i].Hsml = DensityParams.BlackHoleMaxAccretionRadius;
                return;
            }
 
        if(Right[i] < DENSITY_GET_PRIV(tw)->MinGasHsml) {
            P[i].Hsml = DENSITY_GET_PRIV(tw)->MinGasHsml;
            return;
        }
        /* More work needed: add this particle to the redo queue*/
        tw->NPRedo[tid][tw->NPLeft[tid]] = i;
        tw->NPLeft[tid] ++;
        if(tw->NPLeft[tid] > tw->Redo_thread_alloc)
            endrun(5, "Particle %d on thread %d exceeded allocated size of redo queue %ld\n", tw->NPLeft[tid], tid, tw->Redo_thread_alloc);
    }
    else {
        /* We might have got here by serendipity, without bounding.*/
        if(DENSITY_GET_PRIV(tw)->BlackHoleOn && P[i].Type == 5)
            if(P[i].Hsml > DensityParams.BlackHoleMaxAccretionRadius)
                P[i].Hsml = DensityParams.BlackHoleMaxAccretionRadius;
        if(P[i].Hsml < DENSITY_GET_PRIV(tw)->MinGasHsml)
            P[i].Hsml = DENSITY_GET_PRIV(tw)->MinGasHsml;
    }
    if(tw->maxnumngb[tid] < NumNgb[i])
        tw->maxnumngb[tid] = NumNgb[i];
    if(tw->minnumngb[tid] > NumNgb[i])
        tw->minnumngb[tid] = NumNgb[i];
 
    if(tw->Niteration >= MAXITER - 10)
    {
         message(1, "i=%d ID=%lu Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n   pos=(%g|%g|%g)\n",
             i, P[i].ID, P[i].Hsml, Left[i], Right[i],
             NumNgb[i], Right[i] - Left[i], P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
    }
}

References density_params::BlackHoleMaxAccretionRadius, density_params::BlackHoleNgbFactor, ForceTree::BoxSize, DENSITY_GET_PRIV, DensityParams, endrun(), MAXITER, TreeWalk::maxnumngb, density_params::MaxNumNgbDeviation, message(), TreeWalk::minnumngb, TreeWalk::Niteration, TreeWalk::NPLeft, TreeWalk::NPRedo, NUMDIMS, P, TreeWalk::Redo_thread_alloc, and TreeWalk::tree.

Referenced by density_postprocess().

Here is the call graph for this function:

Here is the caller graph for this function:

density_copy()

static void density_copy ( int  place, TreeWalkQueryDensity *  I, TreeWalk *  tw  )

static

Definition at line 316 of file density.c.

{
    I->Hsml = P[place].Hsml;
 
    I->Type = P[place].Type;
 
    if(P[place].Type != 0)
    {
        I->Vel[0] = P[place].Vel[0];
        I->Vel[1] = P[place].Vel[1];
        I->Vel[2] = P[place].Vel[2];
    }
    else
    {
        MyFloat * velpred = DENSITY_GET_PRIV(tw)->SPH_predicted->VelPred;
        I->Vel[0] = velpred[3 * P[place].PI];
        I->Vel[1] = velpred[3 * P[place].PI + 1];
        I->Vel[2] = velpred[3 * P[place].PI + 2];
    }
 
}

References DENSITY_GET_PRIV, TreeWalkQueryDensity::Hsml, P, TreeWalkQueryDensity::Type, and TreeWalkQueryDensity::Vel.

Referenced by density().

Here is the caller graph for this function:

density_haswork()

static int density_haswork ( int  n, TreeWalk *  tw  )

static

Definition at line 501 of file density.c.

{
    /* Don't want a density for swallowed black hole particles*/
    if(P[n].Swallowed)
        return 0;
    if(P[n].Type == 0 || P[n].Type == 5)
        return 1;
    return 0;
}

References P.

Referenced by density().

Here is the caller graph for this function:

density_ngbiter()

static void density_ngbiter ( TreeWalkQueryDensity *  I, TreeWalkResultDensity *  O, TreeWalkNgbIterDensity *  iter, LocalTreeWalk *  lv  )

static

Definition at line 388 of file density.c.

{
    if(iter->base.other == -1) {
        const double h = I->Hsml;
        density_kernel_init(&iter->kernel, h, DensityParams.DensityKernelType);
        iter->kernel_volume = density_kernel_volume(&iter->kernel);
 
        iter->base.Hsml = h;
        iter->base.mask = 1; /* gas only */
        iter->base.symmetric = NGB_TREEFIND_ASYMMETRIC;
        return;
    }
    const int other = iter->base.other;
    const double r = iter->base.r;
    const double r2 = iter->base.r2;
    const double * dist = iter->base.dist;
 
    if(P[other].Mass == 0) {
        endrun(12, "Density found zero mass particle %d type %d id %ld pos %g %g %g\n",
               other, P[other].Type, P[other].ID, P[other].Pos[0], P[other].Pos[1], P[other].Pos[2]);
    }
 
    if(r2 < iter->kernel.HH)
    {
        /* For the BH we wish to exclude wind particles from the density,
         * because they are excluded from the accretion treewalk.*/
        if(I->Type == 5 && winds_is_particle_decoupled(other))
            return;
 
        const double u = r * iter->kernel.Hinv;
        const double wk = density_kernel_wk(&iter->kernel, u);
        O->Ngb += wk * iter->kernel_volume;
 
        const double dwk = density_kernel_dwk(&iter->kernel, u);
 
        const double mass_j = P[other].Mass;
 
        O->Rho += (mass_j * wk);
 
        /* Hinv is here because O->DhsmlDensity is drho / dH.
         * nothing to worry here */
        double density_dW = density_kernel_dW(&iter->kernel, u, wk, dwk);
        O->DhsmlDensity += mass_j * density_dW;
 
        /* For the BH and stars only density and dhsmldensity is used.*/
        if(I->Type != 0)
            return;
 
        struct sph_pred_data * SphP_scratch = DENSITY_GET_PRIV(lv->tw)->SPH_predicted;
 
        double EntVarPred;
        MyFloat VelPred[3];
        #pragma omp atomic read
        EntVarPred = SphP_scratch->EntVarPred[P[other].PI];
        /* Lazily compute the predicted quantities. We can do this
         * with minimal locking since nothing happens should we compute them twice.
         * Zero can be the special value since there should never be zero entropy.*/
        if(EntVarPred == 0) {
            struct DensityPriv * priv = DENSITY_GET_PRIV(lv->tw);
            int bin = P[other].TimeBin;
            double dloga = dloga_from_dti(priv->times->Ti_Current - priv->times->Ti_kick[bin], priv->times->Ti_Current);
            EntVarPred = SPH_EntVarPred(P[other].PI, priv->MinEgySpec, priv->a3inv, dloga);
            SPH_VelPred(other, VelPred, priv->FgravkickB, priv->gravkicks[bin], priv->hydrokicks[bin]);
            /* Note this goes first to avoid threading issues: EntVarPred will only be set after this is done.
             * The worst that can happen is that some data points get copied twice.*/
            int i;
            for(i = 0; i < 3; i++) {
                #pragma omp atomic write
                SphP_scratch->VelPred[3 * P[other].PI + i] = VelPred[i];
            }
            #pragma omp atomic write
            SphP_scratch->EntVarPred[P[other].PI] = EntVarPred;
        }
        else {
            int i;
            for(i = 0; i < 3; i++) {
                #pragma omp atomic read
                VelPred[i] = SphP_scratch->VelPred[3 * P[other].PI + i];
            }
        }
        if(DENSITY_GET_PRIV(lv->tw)->DoEgyDensity) {
            O->EgyRho += mass_j * EntVarPred * wk;
            O->DhsmlEgyDensity += mass_j * EntVarPred * density_dW;
        }
 
        if(r > 0)
        {
            double fac = mass_j * dwk / r;
            double dv[3];
            double rot[3];
            int d;
            for(d = 0; d < 3; d ++) {
                dv[d] = I->Vel[d] - VelPred[d];
            }
            O->Div += -fac * dotproduct(dist, dv);
 
            crossproduct(dv, dist, rot);
            for(d = 0; d < 3; d ++) {
                O->Rot[d] += fac * rot[d];
            }
            if(DENSITY_GET_PRIV(lv->tw)->GradRho) {
                for (d = 0; d < 3; d ++)
                    O->GradRho[d] += fac * dist[d];
            }
        }
    }
}

References DensityPriv::a3inv, TreeWalkNgbIterDensity::base, crossproduct(), DENSITY_GET_PRIV, density_kernel_dW(), density_kernel_dwk(), density_kernel_init(), density_kernel_volume(), density_kernel_wk(), density_params::DensityKernelType, DensityParams, TreeWalkResultDensity::DhsmlDensity, TreeWalkResultDensity::DhsmlEgyDensity, TreeWalkNgbIterBase::dist, TreeWalkResultDensity::Div, dloga, dloga_from_dti(), dotproduct(), dwk, TreeWalkResultDensity::EgyRho, endrun(), sph_pred_data::EntVarPred, DensityPriv::FgravkickB, TreeWalkResultDensity::GradRho, DensityPriv::gravkicks, DensityKernel::Hinv, TreeWalkQueryDensity::Hsml, TreeWalkNgbIterBase::Hsml, DensityPriv::hydrokicks, TreeWalkNgbIterDensity::kernel, kernel(), TreeWalkNgbIterDensity::kernel_volume, TreeWalkNgbIterBase::mask, DensityPriv::MinEgySpec, TreeWalkResultDensity::Ngb, NGB_TREEFIND_ASYMMETRIC, TreeWalkNgbIterBase::other, P, TreeWalkNgbIterBase::r, TreeWalkNgbIterBase::r2, TreeWalkResultDensity::Rho, TreeWalkResultDensity::Rot, SPH_EntVarPred(), SPH_VelPred(), TreeWalkNgbIterBase::symmetric, DriftKickTimes::Ti_Current, DriftKickTimes::Ti_kick, DensityPriv::times, LocalTreeWalk::tw, TreeWalkQueryDensity::Type, TreeWalkQueryDensity::Vel, sph_pred_data::VelPred, winds_is_particle_decoupled(), and wk.

Referenced by density().

Here is the call graph for this function:

Here is the caller graph for this function:

density_postprocess()

static void density_postprocess ( int  i, TreeWalk *  tw  )

static

Definition at line 512 of file density.c.

{
    MyFloat * DhsmlDens = &(DENSITY_GET_PRIV(tw)->DhsmlDensityFactor[i]);
    double density = -1;
    if(P[i].Type == 0)
        density = SPHP(i).Density;
    else if(P[i].Type == 5)
        density = BHP(i).Density;
    if(density <= 0 && DENSITY_GET_PRIV(tw)->NumNgb[i] > 0) {
        endrun(12, "Particle %d type %d has bad density: %g\n", i, P[i].Type, density);
    }
    *DhsmlDens *= P[i].Hsml / (NUMDIMS * density);
    *DhsmlDens = 1 / (1 + *DhsmlDens);
 
    /* Uses DhsmlDensityFactor and changes Hsml, hence the location.*/
    if(DENSITY_GET_PRIV(tw)->update_hsml)
        density_check_neighbours(i, tw);
 
    if(P[i].Type == 0)
    {
        int PI = P[i].PI;
        /*Compute the EgyWeight factors, which are only useful for density independent SPH */
        if(DENSITY_GET_PRIV(tw)->DoEgyDensity) {
            struct sph_pred_data * SphP_scratch = DENSITY_GET_PRIV(tw)->SPH_predicted;
            const double EntPred = SphP_scratch->EntVarPred[P[i].PI];
            if(EntPred <= 0 || SPHP(i).EgyWtDensity <=0)
                endrun(12, "Particle %d has bad predicted entropy: %g or EgyWtDensity: %g\n", i, EntPred, SPHP(i).EgyWtDensity);
            SPHP(i).DhsmlEgyDensityFactor *= P[i].Hsml/ (NUMDIMS * SPHP(i).EgyWtDensity);
            SPHP(i).DhsmlEgyDensityFactor *= - (*DhsmlDens);
            SPHP(i).EgyWtDensity /= EntPred;
        }
        else
            SPHP(i).DhsmlEgyDensityFactor = *DhsmlDens;
 
        MyFloat * Rot = DENSITY_GET_PRIV(tw)->Rot[PI];
        SPHP(i).CurlVel = sqrt(Rot[0] * Rot[0] + Rot[1] * Rot[1] + Rot[2] * Rot[2]) / SPHP(i).Density;
 
        SPHP(i).DivVel /= SPHP(i).Density;
        P[i].DtHsml = (1.0 / NUMDIMS) * SPHP(i).DivVel * P[i].Hsml;
    }
}

References BHP, density(), density_check_neighbours(), DENSITY_GET_PRIV, endrun(), sph_pred_data::EntVarPred, NUMDIMS, DensityPriv::NumNgb, P, and SPHP.

Referenced by density().

Here is the call graph for this function:

Here is the caller graph for this function:

density_reduce()

static void density_reduce ( int  place, TreeWalkResultDensity *  remote, enum TreeWalkReduceMode  mode, TreeWalk *  tw  )

static

Definition at line 339 of file density.c.

{
    TREEWALK_REDUCE(DENSITY_GET_PRIV(tw)->NumNgb[place], remote->Ngb);
    TREEWALK_REDUCE(DENSITY_GET_PRIV(tw)->DhsmlDensityFactor[place], remote->DhsmlDensity);
 
    if(P[place].Type == 0)
    {
        TREEWALK_REDUCE(SPHP(place).Density, remote->Rho);
 
        TREEWALK_REDUCE(SPHP(place).DivVel, remote->Div);
        int pi = P[place].PI;
        TREEWALK_REDUCE(DENSITY_GET_PRIV(tw)->Rot[pi][0], remote->Rot[0]);
        TREEWALK_REDUCE(DENSITY_GET_PRIV(tw)->Rot[pi][1], remote->Rot[1]);
        TREEWALK_REDUCE(DENSITY_GET_PRIV(tw)->Rot[pi][2], remote->Rot[2]);
 
        MyFloat * gradrho = DENSITY_GET_PRIV(tw)->GradRho;
 
        if(gradrho) {
            TREEWALK_REDUCE(gradrho[3*pi], remote->GradRho[0]);
            TREEWALK_REDUCE(gradrho[3*pi+1], remote->GradRho[1]);
            TREEWALK_REDUCE(gradrho[3*pi+2], remote->GradRho[2]);
        }
 
        /*Only used for density independent SPH*/
        if(DENSITY_GET_PRIV(tw)->DoEgyDensity) {
            TREEWALK_REDUCE(SPHP(place).EgyWtDensity, remote->EgyRho);
            TREEWALK_REDUCE(SPHP(place).DhsmlEgyDensityFactor, remote->DhsmlEgyDensity);
        }
    }
    else if(P[place].Type == 5)
    {
        TREEWALK_REDUCE(BHP(place).Density, remote->Rho);
    }
}

References BHP, DENSITY_GET_PRIV, TreeWalkResultDensity::DhsmlDensity, DensityPriv::DhsmlDensityFactor, TreeWalkResultDensity::DhsmlEgyDensity, TreeWalkResultDensity::Div, DensityPriv::DoEgyDensity, TreeWalkResultDensity::EgyRho, TreeWalkResultDensity::GradRho, TreeWalkResultDensity::Ngb, DensityPriv::NumNgb, P, TreeWalkResultDensity::Rho, DensityPriv::Rot, TreeWalkResultDensity::Rot, SPHP, and TREEWALK_REDUCE.

Referenced by density().

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_scratch

)

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:

Variable Documentation

DensityParams

struct density_params DensityParams

static

Definition at line 1 of file density.c.

Referenced by density(), density_check_neighbours(), density_ngbiter(), GetDensityKernelType(), set_density_params(), and set_densitypar().