timestep.h 4.84 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
26
27
28
/*******************************************************************************
 * 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. */
#include "const.h"
#include "debug.h"

Matthieu Schaller's avatar
Matthieu Schaller committed
29
30
31
32
33
34
35
36
/**
 * @brief Compute a valid integer time-step form a given time-step
 *
 * @param new_dt The time-step to convert.
 * @param ti_begin The (integer) start of the previous time-step.
 * @param ti_end The (integer) end of the previous time-step.
 * @param timeBase_inv The inverse of the system's minimal time-step.
 */
37
38
__attribute__((always_inline)) INLINE static int get_integer_timestep(
    float new_dt, int ti_begin, int ti_end, double timeBase_inv) {
39
40

  /* Convert to integer time */
41
  int new_dti = (int)(new_dt * timeBase_inv);
42
43

  /* Recover the current timestep */
44
  const int current_dti = ti_end - ti_begin;
45
46
47
48
49
50
51
52
53
54
55

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

  /* Put this timestep on the time line */
  int dti_timeline = max_nr_timesteps;
  while (new_dti < dti_timeline) dti_timeline /= 2;
  new_dti = dti_timeline;

  /* Make sure we are allowed to increase the timestep size */
  if (new_dti > current_dti) {
56
    if ((max_nr_timesteps - ti_end) % new_dti > 0) new_dti = current_dti;
57
58
59
60
61
  }

  return new_dti;
}

Matthieu Schaller's avatar
Matthieu Schaller committed
62
63
64
65
66
67
/**
 * @brief Compute the new (integer) time-step of a given #gpart
 *
 * @param gp The #gpart.
 * @param e The #engine (used to get some constants).
 */
68
__attribute__((always_inline)) INLINE static int get_gpart_timestep(
69
    const struct gpart *restrict gp, const struct engine *restrict e) {
70
71
72

  const float new_dt_external = gravity_compute_timestep_external(
      e->external_potential, e->physical_constants, gp);
73
74
75
  /* const float new_dt_self = */
  /*     gravity_compute_timestep_self(e->physical_constants, gp); */
  const float new_dt_self = FLT_MAX;  // MATTHIEU
76

77
  float new_dt = min(new_dt_external, new_dt_self);
78
79

  /* Limit timestep within the allowed range */
80
81
  new_dt = min(new_dt, e->dt_max);
  new_dt = max(new_dt, e->dt_min);
82
83
84
85
86
87
88
89

  /* Convert to integer time */
  const int new_dti =
      get_integer_timestep(new_dt, gp->ti_begin, gp->ti_end, e->timeBase_inv);

  return new_dti;
}

Matthieu Schaller's avatar
Matthieu Schaller committed
90
91
92
93
94
95
96
/**
 * @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).
 */
97
__attribute__((always_inline)) INLINE static int get_part_timestep(
98
99
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct engine *restrict e) {
100
101

  /* Compute the next timestep (hydro condition) */
102
  const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
103

104
105
106
107
108
109
  /* Compute the next timestep (cooling condition) */
  float new_dt_cooling = FLT_MAX;
  if (e->policy & engine_policy_cooling)
    new_dt_cooling = cooling_timestep(e->cooling_data, e->physical_constants,
                                      e->internalUnits, p);

110
111
112
113
  /* Compute the next timestep (gravity condition) */
  float new_dt_grav = FLT_MAX;
  if (p->gpart != NULL) {

114
115
    const float new_dt_external = gravity_compute_timestep_external(
        e->external_potential, e->physical_constants, p->gpart);
116
117
    /* const float new_dt_self = */
    /*     gravity_compute_timestep_self(e->physical_constants, p->gpart); */
118
    const float new_dt_self = FLT_MAX;  // MATTHIEU
119

120
    new_dt_grav = min(new_dt_external, new_dt_self);
121
122
123
  }

  /* Final time-step is minimum of hydro and gravity */
124
  float new_dt = min(min(new_dt_hydro, new_dt_cooling), new_dt_grav);
125
126
127

  /* Limit change in h */
  const float dt_h_change =
128
129
      (p->force.h_dt != 0.0f)
          ? fabsf(e->hydro_properties->log_max_h_change * p->h / p->force.h_dt)
130
          : FLT_MAX;
131

132
  new_dt = min(new_dt, dt_h_change);
133
134

  /* Limit timestep within the allowed range */
135
136
  new_dt = min(new_dt, e->dt_max);
  new_dt = max(new_dt, e->dt_min);
137
138

  /* Convert to integer time */
139
140
  const int new_dti =
      get_integer_timestep(new_dt, p->ti_begin, p->ti_end, e->timeBase_inv);
141
142
143
144
145

  return new_dti;
}

#endif /* SWIFT_TIMESTEP_H */