gravity_properties.c 13.3 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/>.
 *
 ******************************************************************************/

/* This object's header. */
#include "gravity_properties.h"

/* Standard headers */
#include <float.h>
#include <math.h>
26
#include <string.h>
27
28
29
30
31
32
33
34

/* Local headers. */
#include "adiabatic_index.h"
#include "common_io.h"
#include "dimension.h"
#include "error.h"
#include "gravity.h"
#include "kernel_gravity.h"
35
#include "kernel_long_gravity.h"
36
#include "restart.h"
37
38

#define gravity_props_default_a_smooth 1.25f
39
40
#define gravity_props_default_r_cut_max 4.5f
#define gravity_props_default_r_cut_min 0.1f
41
#define gravity_props_default_rebuild_frequency 0.01f
42

Matthieu Schaller's avatar
Matthieu Schaller committed
43
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
44
                        const struct phys_const *phys_const,
45
                        const struct cosmology *cosmo, const int with_cosmology,
46
                        const int with_external_potential,
47
                        const int has_baryons, const int has_DM,
48
                        const int is_zoom_simulation, const int periodic) {
49

50
51
52
53
54
55
56
57
  /* Tree updates */
  p->rebuild_frequency =
      parser_get_opt_param_float(params, "Gravity:rebuild_frequency",
                                 gravity_props_default_rebuild_frequency);

  if (p->rebuild_frequency < 0.f || p->rebuild_frequency > 1.f)
    error("Invalid tree rebuild frequency. Must be in [0., 1.]");

58
  /* Tree-PM parameters */
59
60
61
62
63
64
65
66
  if (periodic) {
    p->mesh_size = parser_get_param_int(params, "Gravity:mesh_side_length");
    p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
                                             gravity_props_default_a_smooth);
    p->r_cut_max_ratio = parser_get_opt_param_float(
        params, "Gravity:r_cut_max", gravity_props_default_r_cut_max);
    p->r_cut_min_ratio = parser_get_opt_param_float(
        params, "Gravity:r_cut_min", gravity_props_default_r_cut_min);
67

68
69
70
    /* Some basic checks of what we read */
    if (p->mesh_size % 2 != 0)
      error("The mesh side-length must be an even number.");
71

72
73
    if (p->a_smooth <= 0.)
      error("The mesh smoothing scale 'a_smooth' must be > 0.");
74

75
76
77
78
79
80
81
82
83
    if (2. * p->a_smooth * p->r_cut_max_ratio > p->mesh_size)
      error("Mesh too small given r_cut_max. Should be at least %d cells wide.",
            (int)(2. * p->a_smooth * p->r_cut_max_ratio) + 1);
  } else {
    p->mesh_size = 0;
    p->a_smooth = 0.f;
    p->r_cut_min_ratio = 0.f;
    p->r_cut_max_ratio = 0.f;
  }
84

85
86
87
  /* Time integration */
  p->eta = parser_get_param_float(params, "Gravity:eta");

88
89
  /* Read the choice of multipole acceptance criterion */
  char buffer[32] = {0};
90
  parser_get_param_string(params, "Gravity:MAC", buffer);
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

  if (strcmp(buffer, "adaptive") == 0) {
    p->use_adaptive_tolerance = 1;
  } else if (strcmp(buffer, "geometric") == 0) {
    p->use_adaptive_tolerance = 0;
  } else {
    error(
        "Invalid choice of multipole acceptance criterion: '%s'. Should be "
        "'adaptive' or 'geometric'",
        buffer);
  }

  /* We always start with the geometric MAC */
  p->use_advanced_mac = 0;

  /* Geometric opening angle */
107
  p->theta_crit = parser_get_param_double(params, "Gravity:theta_cr");
108
  if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
109

110
111
112
  /* Adaptive opening angle tolerance */
  p->adaptive_tolerance = parser_get_param_float(params, "Gravity:epsilon_fmm");

113
114
115
116
  /* Are we allowing tree use below softening? */
  p->use_tree_below_softening =
      parser_get_opt_param_int(params, "Gravity:use_tree_below_softening", 1);

117
118
119
  /* Mesh dithering */
  if (periodic && !with_external_potential) {
    p->with_dithering =
120
        parser_get_opt_param_int(params, "Gravity:dithering", 0);
121
    if (p->with_dithering) {
122
123
      p->dithering_ratio =
          parser_get_opt_param_double(params, "Gravity:dithering_ratio", 1.0);
124
    }
125
126
  }

127
  /* Softening parameters */
128
  if (with_cosmology) {
129

130
131
132
133
    if (has_DM) {
      /* Maximal physical softening taken straight from the parameter file */
      p->epsilon_DM_max_physical =
          parser_get_param_double(params, "Gravity:max_physical_DM_softening");
134

135
136
137
138
139
140
141
142
143
144
145
146
147
148
      /* Co-moving softenings taken straight from the parameter file */
      p->epsilon_DM_comoving =
          parser_get_param_double(params, "Gravity:comoving_DM_softening");
    }

    if (has_baryons) {
      /* Maximal physical softening taken straight from the parameter file */
      p->epsilon_baryon_max_physical = parser_get_param_double(
          params, "Gravity:max_physical_baryon_softening");

      /* Co-moving softenings taken straight from the parameter file */
      p->epsilon_baryon_comoving =
          parser_get_param_double(params, "Gravity:comoving_baryon_softening");
    }
149

150
    if (is_zoom_simulation) {
151

152
153
154
155
156
157
158
      /* Compute the comoving softening length for background particles as
       * a fraction of the mean inter-particle density of the background DM
       * particles Since they have variable masses the mass factor will be
       * multiplied in later on. Note that we already multiply in the conversion
       * from Plummer -> real softening length */
      const double ratio_background =
          parser_get_param_double(params, "Gravity:softening_ratio_background");
159

160
161
162
163
164
165
166
      const double mean_matter_density =
          cosmo->Omega_m * cosmo->critical_density_0;

      p->epsilon_background_fac = kernel_gravity_softening_plummer_equivalent *
                                  ratio_background *
                                  cbrt(1. / mean_matter_density);
    }
167

168
  } else {
169

170
171
172
173
174
175
176
177
    if (has_DM) {
      p->epsilon_DM_max_physical =
          parser_get_param_double(params, "Gravity:max_physical_DM_softening");
    }
    if (has_baryons) {
      p->epsilon_baryon_max_physical = parser_get_param_double(
          params, "Gravity:max_physical_baryon_softening");
    }
178

179
180
181
182
183
184
    /* Some gravity models use the DM softening as the one and only softening
       length that exists. So, if we don't have DM (e.g. hydro test or planetary
       physics), we must have a non-zero epsilon. */
    if (!has_DM && has_baryons)
      p->epsilon_DM_max_physical = p->epsilon_baryon_max_physical;

185
186
    p->epsilon_DM_comoving = p->epsilon_DM_max_physical;
    p->epsilon_baryon_comoving = p->epsilon_baryon_max_physical;
187
  }
188

189
190
191
  /* Copy over the gravitational constant */
  p->G_Newton = phys_const->const_newton_G;

192
  /* Set the softening to the current time */
193
  gravity_props_update(p, cosmo);
194
195
}

196
void gravity_props_update_MAC_choice(struct gravity_props *p) {
197

198
199
  /* Now that we have run initial accelerations,
   * switch to the better MAC */
200
  if (p->use_adaptive_tolerance) p->use_advanced_mac = 1;
201
202
203
204
}

void gravity_props_update(struct gravity_props *p,
                          const struct cosmology *cosmo) {
205

206
207
208
209
  /* Current softening length for the high-res. DM particles. */
  double DM_softening, baryon_softening;
  if (p->epsilon_DM_comoving * cosmo->a > p->epsilon_DM_max_physical)
    DM_softening = p->epsilon_DM_max_physical / cosmo->a;
210
  else
211
212
213
214
215
216
217
    DM_softening = p->epsilon_DM_comoving;

  /* Current softening length for the high-res. baryon particles. */
  if (p->epsilon_baryon_comoving * cosmo->a > p->epsilon_baryon_max_physical)
    baryon_softening = p->epsilon_baryon_max_physical / cosmo->a;
  else
    baryon_softening = p->epsilon_baryon_comoving;
218
219

  /* Plummer equivalent -> internal */
220
221
  DM_softening *= kernel_gravity_softening_plummer_equivalent;
  baryon_softening *= kernel_gravity_softening_plummer_equivalent;
222
223

  /* Store things */
224
225
  p->epsilon_DM_cur = DM_softening;
  p->epsilon_baryon_cur = baryon_softening;
226
227
228
229
}

void gravity_props_print(const struct gravity_props *p) {

230
231
  message("Self-gravity scheme: %s", GRAVITY_IMPLEMENTATION);

232
233
  message("Self-gravity scheme: FMM-MM with m-poles of order %d",
          SELF_GRAVITY_MULTIPOLE_ORDER);
234
235
236

  message("Self-gravity time integration: eta=%.4f", p->eta);

237
  message("Self-gravity opening angle:  theta=%.4f", p->theta_crit);
238

239
240
  message("Self-gravity softening functional form: %s",
          kernel_gravity_softening_name);
241

242
  message(
243
244
245
246
247
248
249
250
251
252
253
254
255
256
      "Self-gravity DM comoving softening: epsilon=%.6f (Plummer equivalent: "
      "%.6f)",
      p->epsilon_DM_comoving * kernel_gravity_softening_plummer_equivalent,
      p->epsilon_DM_comoving);

  message(
      "Self-gravity DM maximal physical softening:    epsilon=%.6f (Plummer "
      "equivalent: %.6f)",
      p->epsilon_DM_max_physical * kernel_gravity_softening_plummer_equivalent,
      p->epsilon_DM_max_physical);

  message(
      "Self-gravity baryon comoving softening: epsilon=%.6f (Plummer "
      "equivalent: "
257
      "%.6f)",
258
259
      p->epsilon_baryon_comoving * kernel_gravity_softening_plummer_equivalent,
      p->epsilon_baryon_comoving);
260
261

  message(
262
263
      "Self-gravity baryon maximal physical softening:    epsilon=%.6f "
      "(Plummer "
264
      "equivalent: %.6f)",
265
266
267
      p->epsilon_baryon_max_physical *
          kernel_gravity_softening_plummer_equivalent,
      p->epsilon_baryon_max_physical);
268

269
  message("Self-gravity mesh side-length: N=%d", p->mesh_size);
270
  message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
271

272
273
274
  message("Self-gravity tree cut-off ratio: r_cut_max=%f", p->r_cut_max_ratio);
  message("Self-gravity truncation cut-off ratio: r_cut_min=%f",
          p->r_cut_min_ratio);
275

276
277
  message("Self-gravity mesh dithering ratio: %f", p->dithering_ratio);

278
  message("Self-gravity mesh truncation function: %s",
279
          kernel_long_gravity_truncation_name);
280
281

  message("Self-gravity tree update frequency: f=%f", p->rebuild_frequency);
282
283
284
285
286
287
}

#if defined(HAVE_HDF5)
void gravity_props_print_snapshot(hid_t h_grpgrav,
                                  const struct gravity_props *p) {

288
289
290
291
292
293
294
295
296
297
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
  io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
  io_write_attribute_s(h_grpgrav, "Softening style",
                       kernel_gravity_softening_name);

  io_write_attribute_f(
      h_grpgrav, "Comoving DM softening length [internal units]",
      p->epsilon_DM_comoving * kernel_gravity_softening_plummer_equivalent);
  io_write_attribute_f(
      h_grpgrav,
      "Comoving DM softening length (Plummer equivalent)  [internal units]",
      p->epsilon_DM_comoving);

  io_write_attribute_f(
      h_grpgrav, "Maximal physical DM softening length  [internal units]",
      p->epsilon_DM_max_physical * kernel_gravity_softening_plummer_equivalent);
  io_write_attribute_f(h_grpgrav,
                       "Maximal physical DM softening length (Plummer "
                       "equivalent) [internal units]",
                       p->epsilon_DM_max_physical);

  io_write_attribute_f(
      h_grpgrav, "Comoving baryon softening length [internal units]",
      p->epsilon_baryon_comoving * kernel_gravity_softening_plummer_equivalent);
  io_write_attribute_f(
      h_grpgrav,
      "Comoving baryon softening length (Plummer equivalent)  [internal units]",
      p->epsilon_baryon_comoving);

  io_write_attribute_f(
      h_grpgrav, "Maximal physical baryon softening length  [internal units]",
      p->epsilon_baryon_max_physical *
          kernel_gravity_softening_plummer_equivalent);
  io_write_attribute_f(h_grpgrav,
                       "Maximal physical baryon softening length (Plummer "
                       "equivalent) [internal units]",
                       p->epsilon_baryon_max_physical);

  io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
  io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
  io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
  io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
  io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
  io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio", p->r_cut_min_ratio);
331
  io_write_attribute_f(h_grpgrav, "Mesh dithering ratio", p->dithering_ratio);
332
333
334
335
  io_write_attribute_f(h_grpgrav, "Tree update frequency",
                       p->rebuild_frequency);
  io_write_attribute_s(h_grpgrav, "Mesh truncation function",
                       kernel_long_gravity_truncation_name);
336
337
}
#endif
338
339
340
341
342
343
344
345

/**
 * @brief Write a gravity_props struct to the given FILE as a stream of bytes.
 *
 * @param p the struct
 * @param stream the file stream
 */
void gravity_props_struct_dump(const struct gravity_props *p, FILE *stream) {
Peter W. Draper's avatar
Peter W. Draper committed
346
  restart_write_blocks((void *)p, sizeof(struct gravity_props), 1, stream,
347
                       "gravity", "gravity props");
348
349
350
351
352
353
354
355
356
}

/**
 * @brief Restore a gravity_props struct from the given FILE as a stream of
 * bytes.
 *
 * @param p the struct
 * @param stream the file stream
 */
357
void gravity_props_struct_restore(struct gravity_props *p, FILE *stream) {
Peter W. Draper's avatar
Peter W. Draper committed
358
359
  restart_read_blocks((void *)p, sizeof(struct gravity_props), 1, stream, NULL,
                      "gravity props");
360
}