Skip to content
Snippets Groups Projects
Commit 1b526d3c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

New arguments to the feedback_iact function. New file to apply the synchronization of particles.

parent b28de4e5
Branches
Tags
1 merge request!956Fully implement the timestep limiter (non MPI)
Hide whitespace changes
Inline Side-by-side
src/feedback/EAGLE/feedback_iact.h
+ 18
− 18
View file @ 1b526d3c
......@@ -21,6 +21,7 @@
/* Local includes */
#include "random.h"
#include "timestep_sync.h"
/**
* @brief Density interaction between two particles (non-symmetric).
......@@ -36,13 +37,11 @@
* @param ti_current Current integer time value
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_feedback_density(const float r2, const float *dx,
const float hi, const float hj,
struct spart *restrict si,
const struct part *restrict pj,
const struct xpart *restrict xpj,
const struct cosmology *restrict cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_feedback_density(
const float r2, const float *dx, const float hi, const float hj,
struct spart *restrict si, const struct part *restrict pj,
const struct xpart *restrict xpj, const struct cosmology *restrict cosmo,
const struct engine *e, const integertime_t ti_current) {
/* Get the gas mass. */
const float mj = hydro_get_mass(pj);
......@@ -85,13 +84,11 @@ runner_iact_nonsym_feedback_density(const float r2, const float *dx,
* generator
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
const float hi, const float hj,
const struct spart *restrict si,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *restrict cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_feedback_apply(
const float r2, const float *dx, const float hi, const float hj,
const struct spart *restrict si, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *restrict cosmo,
const struct engine *e, const integertime_t ti_current) {
/* Get r and 1/r. */
const float r_inv = 1.0f / sqrtf(r2);
......@@ -292,10 +289,13 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
/* Impose maximal viscosity */
hydro_diffusive_feedback_reset(pj);
/* message( */
/* "We did some heating! id %llu star id %llu probability %.5e " */
/* "random_num %.5e du %.5e du/ini %.5e", */
/* pj->id, si->id, prob, rand, delta_u, delta_u / u_init); */
message(
"We did some heating! id %llu star id %llu probability %.5e "
"random_num %.5e du %.5e du/ini %.5e",
pj->id, si->id, prob, rand, delta_u, delta_u / u_init);
/* Synchronize the particle on the timeline */
timestep_sync_part(pj, xpj, e, cosmo);
}
}
}
......
src/runner_doiact_functions_stars.h
+ 16
− 16
View file @ 1b526d3c
......@@ -110,9 +110,9 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
IACT_STARS(r2, dx, hi, hj, si, pj, a, H);
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, NULL, cosmo,
ti_current);
e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo,
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo, e,
ti_current);
#endif
}
......@@ -219,9 +219,9 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, NULL, cosmo,
ti_current);
e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo,
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo, e,
ti_current);
#endif
}
......@@ -392,10 +392,10 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, ti_current);
cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
ti_current);
e, ti_current);
#endif
}
} /* loop over the parts in cj. */
......@@ -521,10 +521,10 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hj, hi, spj, pi, xpi,
cosmo, ti_current);
cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hj, hi, spj, pi, xpi, cosmo,
ti_current);
e, ti_current);
#endif
}
} /* loop over the parts in ci. */
......@@ -649,10 +649,10 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, ti_current);
cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
ti_current);
e, ti_current);
#endif
}
} /* loop over the parts in cj. */
......@@ -711,10 +711,10 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, ti_current);
cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
ti_current);
e, ti_current);
#endif
}
} /* loop over the parts in cj. */
......@@ -812,10 +812,10 @@ void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, ti_current);
cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
ti_current);
e, ti_current);
#endif
}
} /* loop over the parts in cj. */
......@@ -903,10 +903,10 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H);
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, pj->h, spi, pj, NULL,
cosmo, ti_current);
cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, pj->h, spi, pj, xpj,
cosmo, ti_current);
cosmo, e, ti_current);
#endif
}
} /* loop over the parts in cj. */
......
src/timestep_sync.h 0 → 100644
+ 128
− 0
View file @ 1b526d3c
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_TIMESTEP_SYNC_H
#define SWIFT_TIMESTEP_SYNC_H
/* Config parameters. */
#include "../config.h"
#include "kick.h"
INLINE static void timestep_sync_part(struct part *p, struct xpart *xp, const struct engine *e,
const struct cosmology *cosmo) {
const int with_cosmology = (e->policy & engine_policy_cosmology);
const integertime_t ti_current = e->ti_current;
const timebin_t max_active_bin = e->max_active_bin;
const double time_base = e->time_base;
/* This particle is already active. Nothing to do here... */
if (p->time_bin <= max_active_bin) return;
message("Synchronizing particle %lld", p->id);
/* We want to make the particle finish it's time-step now. */
/* Start by recovering the start and end point of the particle's time-step. */
const integertime_t old_ti_beg = get_integer_time_begin(ti_current, p->time_bin);
const integertime_t old_ti_end = get_integer_time_end(ti_current, p->time_bin);
/* Old time-step length on the time-line */
const integertime_t old_dti = old_ti_end - old_ti_beg;
/* The actual time-step size this particle will use */
const integertime_t new_ti_beg = old_ti_beg;
const integertime_t new_ti_end = ti_current;
const integertime_t new_dti = new_ti_end - new_ti_beg;
#ifdef SWIFT_DEBUG_CHECKS
/* Some basic safety checks */
if (old_ti_beg >= ti_current)
error(
"Incorrect value for old time-step beginning ti_current=%lld, "
"old_ti_beg=%lld",
ti_current, old_ti_beg);
if (old_ti_end <= ti_current)
error(
"Incorrect value for old time-step end ti_current=%lld, "
"old_ti_end=%lld",
ti_current, old_ti_end);
if (new_ti_end > old_ti_end) error("New end of time-step after the old one");
if (new_dti > old_dti) error("New time-step larger than old one");
#endif
double dt_kick_grav = 0., dt_kick_hydro = 0., dt_kick_therm = 0.,
dt_kick_corr = 0.;
/* Now we need to reverse the kick1... (the dt are negative here) */
if (with_cosmology) {
dt_kick_hydro = -cosmology_get_hydro_kick_factor(cosmo, old_ti_beg,
old_ti_beg + old_dti / 2);
dt_kick_grav = -cosmology_get_grav_kick_factor(cosmo, old_ti_beg,
old_ti_beg + old_dti / 2);
dt_kick_therm = -cosmology_get_therm_kick_factor(cosmo, old_ti_beg,
old_ti_beg + old_dti / 2);
dt_kick_corr = -cosmology_get_corr_kick_factor(cosmo, old_ti_beg,
old_ti_beg + old_dti / 2);
} else {
dt_kick_hydro = -(old_dti / 2) * time_base;
dt_kick_grav = -(old_dti / 2) * time_base;
dt_kick_therm = -(old_dti / 2) * time_base;
dt_kick_corr = -(old_dti / 2) * time_base;
}
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, dt_kick_corr,
e->cosmology, e->hydro_properties, e->entropy_floor,
old_ti_beg + old_dti / 2, old_ti_beg);
/* ...and apply the new one (dt is positiive) */
if (with_cosmology) {
dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, new_ti_beg,
new_ti_beg + new_dti );
dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, new_ti_beg,
new_ti_beg + new_dti );
dt_kick_therm = cosmology_get_therm_kick_factor(cosmo, new_ti_beg,
new_ti_beg + new_dti );
dt_kick_corr = cosmology_get_corr_kick_factor(cosmo, new_ti_beg,
new_ti_beg + new_dti );
} else {
dt_kick_hydro = (new_dti ) * time_base;
dt_kick_grav = (new_dti ) * time_base;
dt_kick_therm = (new_dti ) * time_base;
dt_kick_corr = (new_dti ) * time_base;
}
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, dt_kick_corr,
e->cosmology, e->hydro_properties, e->entropy_floor, new_ti_beg,
new_ti_beg + new_dti);
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_kick != ti_current) error("Particle has not been synchronized correctly.");
#endif
/* The particle is now ready to compute its new time-step size and for the next kick */
}
#endif /* SWIFT_TIMESTEP_SYNC_H */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment