gravity_properties.c 8.08 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
29
30
31
32
33
/*******************************************************************************
 * 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>

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

#define gravity_props_default_a_smooth 1.25f
37
38
#define gravity_props_default_r_cut_max 4.5f
#define gravity_props_default_r_cut_min 0.1f
39
#define gravity_props_default_rebuild_frequency 0.01f
40

Matthieu Schaller's avatar
Matthieu Schaller committed
41
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
42
43
                        const struct cosmology *cosmo, int with_cosmology,
                        int periodic) {
44

45
46
47
48
49
50
51
52
  /* 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.]");

53
  /* Tree-PM parameters */
54
55
56
57
58
59
60
61
  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);
62

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

67
68
    if (p->a_smooth <= 0.)
      error("The mesh smoothing scale 'a_smooth' must be > 0.");
69

70
71
72
73
74
75
76
77
78
    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;
  }
79

80
81
82
  /* Time integration */
  p->eta = parser_get_param_float(params, "Gravity:eta");

83
  /* Opening angle */
84
  p->theta_crit = parser_get_param_double(params, "Gravity:theta");
85
86
  if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
  p->theta_crit2 = p->theta_crit * p->theta_crit;
87
  p->theta_crit_inv = 1. / p->theta_crit;
88

89
  /* Softening parameters */
90
91
92
93
94
95
96
97
98
99
  if (with_cosmology) {
    p->epsilon_comoving =
        parser_get_param_double(params, "Gravity:comoving_softening");
    p->epsilon_max_physical =
        parser_get_param_double(params, "Gravity:max_physical_softening");
  } else {
    p->epsilon_max_physical =
        parser_get_param_double(params, "Gravity:max_physical_softening");
    p->epsilon_comoving = p->epsilon_max_physical;
  }
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

  /* Set the softening to the current time */
  gravity_update(p, cosmo);
}

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

  /* Current softening lengths */
  double softening;
  if (p->epsilon_comoving * cosmo->a > p->epsilon_max_physical)
    softening = p->epsilon_max_physical / cosmo->a;
  else
    softening = p->epsilon_comoving;

  /* Plummer equivalent -> internal */
  softening *= kernel_gravity_softening_plummer_equivalent;

  /* Store things */
  p->epsilon_cur = softening;

  /* Other factors */
  p->epsilon_cur2 = softening * softening;
  p->epsilon_cur_inv = 1. / softening;
  p->epsilon_cur_inv3 = 1. / (softening * softening * softening);
124
125
126
127
}

void gravity_props_print(const struct gravity_props *p) {

128
129
  message("Self-gravity scheme: %s", GRAVITY_IMPLEMENTATION);

130
131
  message("Self-gravity scheme: FMM-MM with m-poles of order %d",
          SELF_GRAVITY_MULTIPOLE_ORDER);
132
133
134

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

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

137
138
  message("Self-gravity softening functional form: %s",
          kernel_gravity_softening_name);
139

140
141
142
143
144
145
146
147
148
149
150
  message(
      "Self-gravity comoving softening:    epsilon=%.4f (Plummer equivalent: "
      "%.4f)",
      p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent,
      p->epsilon_comoving);

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

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

155
156
157
  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);
158
159

  message("Self-gravity mesh truncation function: %s",
160
          kernel_long_gravity_truncation_name);
161
162

  message("Self-gravity tree update frequency: f=%f", p->rebuild_frequency);
163
164
165
166
167
168
169
}

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

  io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
170
171
  io_write_attribute_s(h_grpgrav, "Softening style",
                       kernel_gravity_softening_name);
172
  io_write_attribute_f(
173
      h_grpgrav, "Comoving softening length [internal units]",
174
175
      p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
  io_write_attribute_f(
176
177
178
179
180
      h_grpgrav,
      "Comoving Softening length (Plummer equivalent)  [internal units]",
      p->epsilon_comoving);
  io_write_attribute_f(
      h_grpgrav, "Maximal physical softening length  [internal units]",
181
182
      p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
  io_write_attribute_f(h_grpgrav,
183
184
                       "Maximal physical softening length (Plummer equivalent) "
                       " [internal units]",
185
                       p->epsilon_max_physical);
186
  io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
187
  io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
188
  io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
189
  io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
190
191
  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);
192
193
194
  io_write_attribute_f(h_grpgrav, "Tree update frequency",
                       p->rebuild_frequency);
  io_write_attribute_s(h_grpgrav, "Mesh truncation function",
195
                       kernel_long_gravity_truncation_name);
196
197
}
#endif
198
199
200
201
202
203
204
205

/**
 * @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
206
  restart_write_blocks((void *)p, sizeof(struct gravity_props), 1, stream,
207
                       "gravity", "gravity props");
208
209
210
211
212
213
214
215
216
}

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