drift.c File Reference

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include "drift.h"
#include "blackhole.h"
#include "walltime.h"
#include "timefac.h"
#include "timestep.h"
#include "utils.h"

Include dependency graph for drift.c:

Go to the source code of this file.

Functions
void	real_drift_particle (struct particle_data *pp, struct slots_manager_type *sman, const double ddrift, const double BoxSize, const double random_shift[3])
void	drift_all_particles (inttime_t ti0, inttime_t ti1, Cosmology *CP, const double random_shift[3])

Function Documentation

drift_all_particles()

void drift_all_particles	(	inttime_t	ti0,
		inttime_t	ti1,
		Cosmology *	CP,
		const double	random_shift[3]
	)

Definition at line 91 of file drift.c.

{
    int i;
    walltime_measure("/Misc");
    if(ti1 < ti0) {
        endrun(12, "Trying to reverse time: ti0=%d ti1=%d\n", ti0, ti1);
    }
    const double ddrift = get_exact_drift_factor(CP, ti0, ti1);
 
#pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++) {
#ifdef DEBUG
        if(PartManager->Base[i].Ti_drift != ti0)
            endrun(10, "Drift time mismatch: (ids = %ld %ld) %d != %d\n",PartManager->Base[0].ID, PartManager->Base[i].ID, ti0,  PartManager->Base[i].Ti_drift);
#endif
        real_drift_particle(&PartManager->Base[i], SlotsManager, ddrift, PartManager->BoxSize, random_shift);
        PartManager->Base[i].Ti_drift = ti1;
    }
 
    walltime_measure("/Drift/All");
}

References part_manager_type::Base, part_manager_type::BoxSize, CP, endrun(), get_exact_drift_factor(), particle_data::ID, part_manager_type::NumPart, PartManager, real_drift_particle(), SlotsManager, particle_data::Ti_drift, and walltime_measure.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

real_drift_particle()

void real_drift_particle	(	struct particle_data *	pp,
		struct slots_manager_type *	sman,
		const double	ddrift,
		const double	BoxSize,
		const double	random_shift[3]
	)

Definition at line 22 of file drift.c.

{
    int j;
    if(pp->IsGarbage || pp->Swallowed) {
        /* Keep the random shift updated so the
         * physical position of swallowed particles remains unchanged.*/
        for(j = 0; j < 3; j++) {
            pp->Pos[j] += random_shift[j];
            while(pp->Pos[j] > BoxSize) pp->Pos[j] -= BoxSize;
            while(pp->Pos[j] <= 0) pp->Pos[j] += BoxSize;
        }
        /* Swallowed particles still need a peano key.*/
        if(pp->Swallowed)
            pp->Key = PEANO(pp->Pos, BoxSize);
        return;
    }
 
    /* Jumping of BH */
    if(BHGetRepositionEnabled() && pp->Type == 5) {
        int k;
        int pi = pp->PI;
        struct bh_particle_data * BH = (struct bh_particle_data *) sman->info[5].ptr;
        if (BH[pi].JumpToMinPot) {
            for(k = 0; k < 3; k++) {
                double dx = NEAREST(pp->Pos[k] - BH[pi].MinPotPos[k], BoxSize);
                if(dx > 0.1 * BoxSize) {
                    endrun(1, "Drifting blackhole very far, from %g %g %g to %g %g %g id = %ld. Likely due to the time step is too sparse.\n",
                        pp->Pos[0],
                        pp->Pos[1],
                        pp->Pos[2],
                        BH[pi].MinPotPos[0],
                        BH[pi].MinPotPos[1],
                        BH[pi].MinPotPos[2], pp->ID);
                }
                pp->Pos[k] = BH[pi].MinPotPos[k];
                pp->Vel[k] = BH[pi].MinPotVel[k];
            }
        }
        BH[pi].JumpToMinPot = 0;
    }
    else if(pp->Type == 0)
    {
        /* DtHsml is 1/3 DivVel * Hsml evaluated at the last active timestep for this particle.
         * This predicts Hsml during the current timestep in the way used in Gadget-4, more accurate
         * than the Gadget-2 prediction which could run away in deep timesteps. */
        pp->Hsml += pp->DtHsml * ddrift;
        if(pp->Hsml <= 0)
            endrun(5, "Part id %ld has bad Hsml %g with DtHsml %g vel %g %g %g\n",
                   pp->ID, pp->Hsml, pp->DtHsml, pp->Vel[0], pp->Vel[1], pp->Vel[2]);
        /* Cap the Hsml just in case: if DivVel is large for a particle with a long timestep
         * at one point Hsml could rarely run away.*/
        const double Maxhsml = BoxSize /2.;
        if(pp->Hsml > Maxhsml)
            pp->Hsml = Maxhsml;
    }
 
    for(j = 0; j < 3; j++) {
        pp->Pos[j] += pp->Vel[j] * ddrift + random_shift[j];
    }
 
    for(j = 0; j < 3; j ++) {
        while(pp->Pos[j] > BoxSize) pp->Pos[j] -= BoxSize;
        while(pp->Pos[j] <= 0) pp->Pos[j] += BoxSize;
    }
    /* avoid recomputing them during layout and force tree build.*/
    pp->Key = PEANO(pp->Pos, BoxSize);
}

References BHGetRepositionEnabled(), particle_data::DtHsml, endrun(), particle_data::Hsml, particle_data::ID, slots_manager_type::info, particle_data::IsGarbage, bh_particle_data::JumpToMinPot, particle_data::Key, bh_particle_data::MinPotPos, bh_particle_data::MinPotVel, NEAREST, PEANO(), particle_data::PI, particle_data::Pos, slot_info::ptr, particle_data::Swallowed, particle_data::Type, and particle_data::Vel.

Referenced by domain_build_exchange_list(), and drift_all_particles().

Here is the call graph for this function:

Here is the caller graph for this function: