timestep.h 12.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 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_H
#define SWIFT_TIMESTEP_H

/* Config parameters. */
#include "../config.h"

/* Local headers. */
Tom Theuns's avatar
Tom Theuns committed
26
#include "cooling.h"
27
#include "debug.h"
28 29
#include "timeline.h"

Matthieu Schaller's avatar
Matthieu Schaller committed
30 31 32
/**
 * @brief Compute a valid integer time-step form a given time-step
 *
33
 * We consider the minimal time-bin of any neighbours and prevent particles
34
 * to differ from it by a fixed constant `time_bin_neighbour_max_delta_bin`.
35 36 37 38
 *
 * If min_ngb_bin is set to `num_time_bins`, then no limit from the neighbours
 * is imposed.
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
39
 * @param new_dt The time-step to convert.
40
 * @param old_bin The old time bin.
41
 * @param min_ngb_bin Minimal time-bin of any neighbour of this particle.
42
 * @param ti_current The current time on the integer time-line.
43
 * @param time_base_inv The inverse of the system's minimal time-step.
Matthieu Schaller's avatar
Matthieu Schaller committed
44
 */
45 46 47 48 49
__attribute__((always_inline, const)) INLINE static integertime_t
make_integer_timestep(const float new_dt, const timebin_t old_bin,
                      const timebin_t min_ngb_bin,
                      const integertime_t ti_current,
                      const double time_base_inv) {
50 51

  /* Convert to integer time */
52
  integertime_t new_dti = (integertime_t)(new_dt * time_base_inv);
53

54 55
  /* Are we allowed to use this bin given the neighbours? */
  timebin_t new_bin = get_time_bin(new_dti);
56
  new_bin = min(new_bin, min_ngb_bin + time_bin_neighbour_max_delta_bin);
57 58
  new_dti = get_integer_timestep(new_bin);

59
  /* Current time-step */
60 61
  const integertime_t current_dti = get_integer_timestep(old_bin);
  const integertime_t ti_end = get_integer_time_end(ti_current, old_bin);
62 63

  /* Limit timestep increase */
64
  if (old_bin > 0) new_dti = min(new_dti, 2 * current_dti);
65 66

  /* Put this timestep on the time line */
67
  integertime_t dti_timeline = max_nr_timesteps;
68
  while (new_dti < dti_timeline) dti_timeline /= ((integertime_t)2);
69 70 71 72
  new_dti = dti_timeline;

  /* Make sure we are allowed to increase the timestep size */
  if (new_dti > current_dti) {
73
    if ((max_nr_timesteps - ti_end) % new_dti > 0) new_dti = current_dti;
74
  }
75 76 77 78 79

#ifdef SWIFT_DEBUG_CHECKS
  if (new_dti == 0) error("Computed an integer time-step of size 0");
#endif

80 81 82
  return new_dti;
}

Matthieu Schaller's avatar
Matthieu Schaller committed
83 84 85 86 87 88
/**
 * @brief Compute the new (integer) time-step of a given #gpart
 *
 * @param gp The #gpart.
 * @param e The #engine (used to get some constants).
 */
89
__attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
90
    const struct gpart *restrict gp, const struct engine *restrict e) {
91

92
  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
93

94
  if (e->policy & engine_policy_external_gravity)
95 96
    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
                                           e->physical_constants, gp);
97

98
  const float a_hydro[3] = {0.f, 0.f, 0.f};
99
  if (e->policy & engine_policy_self_gravity)
100
    new_dt_self = gravity_compute_timestep_self(
101
        gp, a_hydro, e->gravity_properties, e->cosmology);
102 103 104 105

  /* Take the minimum of all */
  float new_dt = min(new_dt_self, new_dt_ext);

106 107 108
  /* Apply the maximal displacement constraint (FLT_MAX  if non-cosmological)*/
  new_dt = min(new_dt, e->dt_max_RMS_displacement);

109
  /* Apply cosmology correction (This is 1 if non-cosmological) */
110
  new_dt *= e->cosmology->time_step_factor;
111 112

  /* Limit timestep within the allowed range */
113
  new_dt = min(new_dt, e->dt_max);
114 115 116
  if (new_dt < e->dt_min)
    error("gpart (id=%lld) wants a time-step (%e) below dt_min (%e)",
          gp->id_or_neg_offset, new_dt, e->dt_min);
117 118

  /* Convert to integer time */
119
  const integertime_t new_dti = make_integer_timestep(
120
      new_dt, gp->time_bin, num_time_bins, e->ti_current, e->time_base_inv);
121 122 123 124

  return new_dti;
}

Matthieu Schaller's avatar
Matthieu Schaller committed
125 126 127 128 129 130 131
/**
 * @brief Compute the new (integer) time-step of a given #part
 *
 * @param p The #part.
 * @param xp The #xpart partner of p.
 * @param e The #engine (used to get some constants).
 */
132
__attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
133 134
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct engine *restrict e) {
135 136

  /* Compute the next timestep (hydro condition) */
137 138
  const float new_dt_hydro =
      hydro_compute_timestep(p, xp, e->hydro_properties, e->cosmology);
139

140 141 142
  /* Compute the next timestep (cooling condition) */
  float new_dt_cooling = FLT_MAX;
  if (e->policy & engine_policy_cooling)
143 144
    new_dt_cooling =
        cooling_timestep(e->cooling_func, e->physical_constants, e->cosmology,
145
                         e->internal_units, e->hydro_properties, p, xp);
146

147
  /* Compute the next timestep (gravity condition) */
148 149
  float new_dt_grav = FLT_MAX, new_dt_self_grav = FLT_MAX,
        new_dt_ext_grav = FLT_MAX;
150 151
  if (p->gpart != NULL) {

152
    if (e->policy & engine_policy_external_gravity)
153 154
      new_dt_ext_grav = external_gravity_timestep(
          e->time, e->external_potential, e->physical_constants, p->gpart);
155

156
    if (e->policy & engine_policy_self_gravity)
157
      new_dt_self_grav = gravity_compute_timestep_self(
158
          p->gpart, p->a_hydro, e->gravity_properties, e->cosmology);
159 160

    new_dt_grav = min(new_dt_self_grav, new_dt_ext_grav);
161 162
  }

163 164 165 166 167 168 169 170
  /* Compute the next timestep (chemistry condition, e.g. diffusion) */
  const float new_dt_chemistry =
      chemistry_timestep(e->physical_constants, e->cosmology, e->internal_units,
                         e->hydro_properties, e->chemistry, p);

  /* Final time-step is minimum of hydro, gravity and subgrid */
  float new_dt =
      min4(new_dt_hydro, new_dt_cooling, new_dt_grav, new_dt_chemistry);
171

172
  /* Limit change in smoothing length */
173
  const float dt_h_change =
174 175
      (p->force.h_dt != 0.0f)
          ? fabsf(e->hydro_properties->log_max_h_change * p->h / p->force.h_dt)
176
          : FLT_MAX;
177

178
  new_dt = min(new_dt, dt_h_change);
179

180
  /* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/
181 182
  new_dt = min(new_dt, e->dt_max_RMS_displacement);

183
  /* Apply cosmology correction (This is 1 if non-cosmological) */
184
  new_dt *= e->cosmology->time_step_factor;
185

186
  /* Limit timestep within the allowed range */
187
  new_dt = min(new_dt, e->dt_max);
188

189 190 191
  if (new_dt < e->dt_min)
    error("part (id=%lld) wants a time-step (%e) below dt_min (%e)", p->id,
          new_dt, e->dt_min);
192 193

  /* Convert to integer time */
194 195 196
  const integertime_t new_dti = make_integer_timestep(
      new_dt, p->time_bin, p->limiter_data.min_ngb_time_bin, e->ti_current,
      e->time_base_inv);
197

198
  return new_dti;
199 200
}

201 202 203 204 205 206 207 208 209
/**
 * @brief Compute the new (integer) time-step of a given #spart
 *
 * @param sp The #spart.
 * @param e The #engine (used to get some constants).
 */
__attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
    const struct spart *restrict sp, const struct engine *restrict e) {

210
  /* Stellar time-step */
211 212 213
  float new_dt_stars = stars_compute_timestep(
      sp, e->stars_properties, (e->policy & engine_policy_cosmology),
      e->cosmology, e->time);
214 215 216

  /* Gravity time-step */
  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
217 218

  if (e->policy & engine_policy_external_gravity)
219 220
    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
                                           e->physical_constants, sp->gpart);
221

222
  const float a_hydro[3] = {0.f, 0.f, 0.f};
223
  if (e->policy & engine_policy_self_gravity)
224
    new_dt_self = gravity_compute_timestep_self(
225
        sp->gpart, a_hydro, e->gravity_properties, e->cosmology);
226 227

  /* Take the minimum of all */
228
  float new_dt = min3(new_dt_stars, new_dt_self, new_dt_ext);
229

230 231 232
  /* Apply the maximal displacement constraint (FLT_MAX  if non-cosmological)*/
  new_dt = min(new_dt, e->dt_max_RMS_displacement);

233
  /* Apply cosmology correction (This is 1 if non-cosmological) */
234
  new_dt *= e->cosmology->time_step_factor;
235 236

  /* Limit timestep within the allowed range */
237
  new_dt = min(new_dt, e->dt_max);
238
  if (new_dt < e->dt_min) {
239 240
    error("spart (id=%lld) wants a time-step (%e) below dt_min (%e)", sp->id,
          new_dt, e->dt_min);
241
  }
242 243 244

  /* Convert to integer time */
  const integertime_t new_dti = make_integer_timestep(
245
      new_dt, sp->time_bin, num_time_bins, e->ti_current, e->time_base_inv);
246 247 248 249

  return new_dti;
}

250 251 252 253 254 255 256 257 258
/**
 * @brief Compute the new (integer) time-step of a given #bpart
 *
 * @param bp The #bpart.
 * @param e The #engine (used to get some constants).
 */
__attribute__((always_inline)) INLINE static integertime_t get_bpart_timestep(
    const struct bpart *restrict bp, const struct engine *restrict e) {

259 260 261
  /* Black hole internal time-step */
  float new_dt_black_holes = black_holes_compute_timestep(
      bp, e->black_holes_properties, e->physical_constants);
262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292

  /* Gravity time-step */
  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;

  if (e->policy & engine_policy_external_gravity)
    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
                                           e->physical_constants, bp->gpart);

  const float a_hydro[3] = {0.f, 0.f, 0.f};
  if (e->policy & engine_policy_self_gravity)
    new_dt_self = gravity_compute_timestep_self(
        bp->gpart, a_hydro, e->gravity_properties, e->cosmology);

  /* Take the minimum of all */
  float new_dt = min3(new_dt_black_holes, new_dt_self, new_dt_ext);

  /* Apply the maximal dibslacement constraint (FLT_MAX  if non-cosmological)*/
  new_dt = min(new_dt, e->dt_max_RMS_displacement);

  /* Apply cosmology correction (This is 1 if non-cosmological) */
  new_dt *= e->cosmology->time_step_factor;

  /* Limit timestep within the allowed range */
  new_dt = min(new_dt, e->dt_max);
  if (new_dt < e->dt_min) {
    error("bpart (id=%lld) wants a time-step (%e) below dt_min (%e)", bp->id,
          new_dt, e->dt_min);
  }

  /* Convert to integer time */
  const integertime_t new_dti = make_integer_timestep(
293
      new_dt, bp->time_bin, num_time_bins, e->ti_current, e->time_base_inv);
294 295 296 297

  return new_dti;
}

298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344
/**
 * @brief Compute the new (integer) time-step of a given #sink.
 *
 * @param sink The #sink.
 * @param e The #engine (used to get some constants).
 */
__attribute__((always_inline)) INLINE static integertime_t get_sink_timestep(
    const struct sink *restrict sink, const struct engine *restrict e) {

  /* Sink time-step */
  float new_dt_sink = sink_compute_timestep(sink);

  /* Gravity time-step */
  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;

  if (e->policy & engine_policy_external_gravity)
    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
                                           e->physical_constants, sink->gpart);

  const float a_hydro[3] = {0.f, 0.f, 0.f};
  if (e->policy & engine_policy_self_gravity)
    new_dt_self = gravity_compute_timestep_self(
        sink->gpart, a_hydro, e->gravity_properties, e->cosmology);

  /* Take the minimum of all */
  float new_dt = min3(new_dt_sink, new_dt_self, new_dt_ext);

  /* Apply the maximal dibslacement constraint (FLT_MAX  if non-cosmological)*/
  new_dt = min(new_dt, e->dt_max_RMS_displacement);

  /* Apply cosmology correction (This is 1 if non-cosmological) */
  new_dt *= e->cosmology->time_step_factor;

  /* Limit timestep within the allowed range */
  new_dt = min(new_dt, e->dt_max);
  if (new_dt < e->dt_min) {
    error("sink (id=%lld) wants a time-step (%e) below dt_min (%e)", sink->id,
          new_dt, e->dt_min);
  }

  /* Convert to integer time */
  const integertime_t new_dti = make_integer_timestep(
      new_dt, sink->time_bin, num_time_bins, e->ti_current, e->time_base_inv);

  return new_dti;
}

345
#endif /* SWIFT_TIMESTEP_H */