gravity.h 9.45 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2016 Tom Theuns (tom.theuns@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/>.
 *
 ******************************************************************************/
20
#ifndef SWIFT_MULTI_SOFTENING_GRAVITY_H
Josh Borrow's avatar
Josh Borrow committed
21
#define SWIFT_MULTI_SOFTENING_GRAVITY_H
22
23
24
25
26

#include <float.h>

/* Local includes. */
#include "cosmology.h"
27
#include "error.h"
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include "gravity_properties.h"
#include "kernel_gravity.h"
#include "minmax.h"

/**
 * @brief Returns the mass of a particle
 *
 * @param gp The particle of interest
 */
__attribute__((always_inline)) INLINE static float gravity_get_mass(
    const struct gpart* restrict gp) {

  return gp->mass;
}

/**
 * @brief Returns the current co-moving softening of a particle
 *
 * @param gp The particle of interest
 * @param grav_props The global gravity properties.
 */
__attribute__((always_inline)) INLINE static float gravity_get_softening(
    const struct gpart* gp, const struct gravity_props* restrict grav_props) {

52
  return gp->epsilon;
53
54
55
56
57
58
59
60
61
}

/**
 * @brief Add a contribution to this particle's potential.
 *
 * @param gp The particle.
 * @param pot The contribution to add.
 */
__attribute__((always_inline)) INLINE static void
62
63
gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {

64
  gp->potential += pot;
65
}
66
67
68
69
70
71
72
73
74

/**
 * @brief Returns the comoving potential of a particle.
 *
 * @param gp The particle of interest
 */
__attribute__((always_inline)) INLINE static float
gravity_get_comoving_potential(const struct gpart* restrict gp) {

75
  return gp->potential;
76
77
78
79
80
81
82
83
84
85
86
87
}

/**
 * @brief Returns the physical potential of a particle
 *
 * @param gp The particle of interest.
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
gravity_get_physical_potential(const struct gpart* restrict gp,
                               const struct cosmology* cosmo) {

88
  return gp->potential * cosmo->a_inv;
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
}

/**
 * @brief Computes the gravity time-step of a given particle due to self-gravity
 *
 * We use Gadget-2's type 0 time-step criterion.
 *
 * @param gp Pointer to the g-particle data.
 * @param a_hydro The accelerations coming from the hydro scheme (can be 0).
 * @param grav_props Constants used in the gravity scheme.
 * @param cosmo The current cosmological model.
 */
__attribute__((always_inline)) INLINE static float
gravity_compute_timestep_self(const struct gpart* const gp,
                              const float a_hydro[3],
                              const struct gravity_props* restrict grav_props,
                              const struct cosmology* cosmo) {

  /* Get physical acceleration (gravity contribution) */
  float a_phys_x = gp->a_grav[0] * cosmo->a_factor_grav_accel;
  float a_phys_y = gp->a_grav[1] * cosmo->a_factor_grav_accel;
  float a_phys_z = gp->a_grav[2] * cosmo->a_factor_grav_accel;

  /* Get physical acceleration (hydro contribution) */
  a_phys_x += a_hydro[0] * cosmo->a_factor_hydro_accel;
  a_phys_y += a_hydro[1] * cosmo->a_factor_hydro_accel;
  a_phys_z += a_hydro[2] * cosmo->a_factor_hydro_accel;

  const float ac2 =
      a_phys_x * a_phys_x + a_phys_y * a_phys_y + a_phys_z * a_phys_z;

  const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;

  const float epsilon = gravity_get_softening(gp, grav_props);

  const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
                         cosmo->a * grav_props->eta * epsilon * ac_inv);

  return dt;
}

/**
 * @brief Prepares a g-particle for the gravity calculation
 *
 * Zeroes all the relevant arrays in preparation for the sums taking place in
 * the variaous tasks
 *
 * @param gp The particle to act upon
 */
__attribute__((always_inline)) INLINE static void gravity_init_gpart(
    struct gpart* gp) {

  /* Zero the acceleration */
  gp->a_grav[0] = 0.f;
  gp->a_grav[1] = 0.f;
  gp->a_grav[2] = 0.f;
145
  gp->potential = 0.f;
146
147
148

#ifdef SWIFT_GRAVITY_FORCE_CHECKS
  gp->potential_PM = 0.f;
149
150

  /* Track accelerations of each component. */
151
  for (int i = 0; i < 3; i++) {
152
153
154
155
156
157
158
159
160
161
    gp->a_grav_PM[i] = 0.f;
    gp->a_grav_p2p[i] = 0.f;
    gp->a_grav_m2p[i] = 0.f;
    gp->a_grav_m2l[i] = 0.f;
  }

  /* Interaction counters. */
  gp->num_interacted_m2p = 0;
  gp->num_interacted_m2l = 0;
  gp->num_interacted_p2p = 0;
162
  gp->num_interacted_pm = 0;
163
164
165
#endif

#ifdef SWIFT_DEBUG_CHECKS
166
  gp->num_interacted = 0;
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
  gp->initialised = 1;
#endif
}

/**
 * @brief Finishes the gravity calculation.
 *
 * Multiplies the forces and accelerations by the appropiate constants.
 * Applies cosmological correction for periodic BCs.
 *
 * No need to apply the potential normalisation correction for periodic
 * BCs here since we do not compute the potential.
 *
 * @param gp The particle to act upon
 * @param const_G Newton's constant in internal units.
 * @param potential_normalisation Term to be added to all the particles.
 * @param periodic Are we using periodic BCs?
184
 * @param with_self_gravity Are we running with self-gravity?
185
186
 */
__attribute__((always_inline)) INLINE static void gravity_end_force(
187
188
    struct gpart* gp, const float const_G, const float potential_normalisation,
    const int periodic, const int with_self_gravity) {
189

190
191
192
  /* Apply the periodic correction to the peculiar potential */
  if (periodic) gp->potential += potential_normalisation;

193
194
195
196
197
198
199
200
  /* Record the norm of the acceleration for the adaptive opening criteria.
   * Will always be an (active) timestep behind. */
  gp->old_a_grav_norm = gp->a_grav[0] * gp->a_grav[0] +
                        gp->a_grav[1] * gp->a_grav[1] +
                        gp->a_grav[2] * gp->a_grav[2];

  gp->old_a_grav_norm = sqrtf(gp->old_a_grav_norm);

201
#ifdef SWIFT_DEBUG_CHECKS
202
203
  if (with_self_gravity && gp->old_a_grav_norm == 0.f)
    error("Old acceleration is 0!");
204
205
#endif

206
207
208
209
  /* Let's get physical... */
  gp->a_grav[0] *= const_G;
  gp->a_grav[1] *= const_G;
  gp->a_grav[2] *= const_G;
210
  gp->potential *= const_G;
211
212
213

#ifdef SWIFT_GRAVITY_FORCE_CHECKS
  gp->potential_PM *= const_G;
214
215
216
217
218
219
  for (int i = 0; i < 3; i++) {
    gp->a_grav_PM[i] *= const_G;
    gp->a_grav_p2p[i] *= const_G;
    gp->a_grav_m2p[i] *= const_G;
    gp->a_grav_m2l[i] *= const_G;
  }
220
221
222
223
224
225
226
#endif

#ifdef SWIFT_DEBUG_CHECKS
  gp->initialised = 0; /* Ready for next step */
#endif
}

227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
/**
 * @brief Update the #gpart after a drift step.
 *
 * This is typically used to update the softening lengths.
 *
 * @param gp The particle to act upon
 * @param grav_props The global properties of the gravity calculation.
 */
__attribute__((always_inline)) INLINE static void gravity_predict_extra(
    struct gpart* gp, const struct gravity_props* grav_props) {

  switch (gp->type) {
    case swift_type_dark_matter:
      gp->epsilon = grav_props->epsilon_DM_cur;
      break;
    case swift_type_stars:
      gp->epsilon = grav_props->epsilon_baryon_cur;
      break;
    case swift_type_gas:
      gp->epsilon = grav_props->epsilon_baryon_cur;
      break;
    case swift_type_black_hole:
      gp->epsilon = grav_props->epsilon_baryon_cur;
      break;
    case swift_type_dark_matter_background:
      gp->epsilon = grav_props->epsilon_background_fac * cbrtf(gp->mass);
      break;
    default:
#ifdef SWIFT_DEBUG_CHECKS
      error("Invalid gpart type!");
#endif
      break;
  }
}

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
/**
 * @brief Kick the additional variables
 *
 * @param gp The particle to act upon
 * @param dt The time-step for this kick
 */
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
    struct gpart* gp, float dt) {}

/**
 * @brief Sets the values to be predicted in the drifts to their values at a
 * kick time
 *
 * @param gp The particle.
 */
__attribute__((always_inline)) INLINE static void
gravity_reset_predicted_values(struct gpart* gp) {}

/**
 * @brief Initialises the g-particles for the first time
 *
 * This function is called only once just after the ICs have been
 * read in to do some conversions.
 *
 * @param gp The particle to act upon
 * @param grav_props The global properties of the gravity calculation.
 */
__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
    struct gpart* gp, const struct gravity_props* grav_props) {

  gp->time_bin = 0;
293
  gp->old_a_grav_norm = 0.f;
294

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
  switch (gp->type) {
    case swift_type_dark_matter:
      gp->epsilon = grav_props->epsilon_DM_cur;
      break;
    case swift_type_stars:
      gp->epsilon = grav_props->epsilon_baryon_cur;
      break;
    case swift_type_gas:
      gp->epsilon = grav_props->epsilon_baryon_cur;
      break;
    case swift_type_black_hole:
      gp->epsilon = grav_props->epsilon_baryon_cur;
      break;
    case swift_type_dark_matter_background:
      gp->epsilon = grav_props->epsilon_background_fac * cbrtf(gp->mass);
      break;
    default:
#ifdef SWIFT_DEBUG_CHECKS
      error("Invalid gpart type!");
#endif
      break;
  }

318
319
320
  gravity_init_gpart(gp);
}

321
#endif /* SWIFT_MULTI_SOFTENING_GRAVITY_H */