potential.h 6.8 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
34
35
36
37
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
 *                    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_DISC_PATCH_H
#define SWIFT_DISC_PATCH_H

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

/* Some standard headers. */
#include <float.h>
#include <math.h>

/* Local includes. */
#include "const.h"
#include "error.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"

38
39
40
41
42
/**
 * @brief External Potential Properties - Disc patch case
 *
 * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
 */
43
44
struct external_potential {

45
  /*! Surface density of the disc */
46
  double surface_density;
47
48

  /*! Disc scale-height */
49
  double scale_height;
50
51

  /*! Position of the disc along the z-axis */
52
  double z_disc;
53
54

  /*! Dynamical time of the system */
55
  double dynamical_time;
56
57

  /*! Time over which to grow the disk in units of the dynamical time */
58
  double growth_time;
59
60

  /*! Time-step condition pre-factor */
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
  double timestep_mult;
};

/**
 * @brief Computes the time-step from the acceleration due to a hydrostatic
 * disc.
 *
 * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
 *
 * @param time The current time.
 * @param potential The properties of the potential.
 * @param phys_const The physical constants in internal units.
 * @param g Pointer to the g-particle data.
 */
__attribute__((always_inline)) INLINE static float external_gravity_timestep(
    double time, const struct external_potential* restrict potential,
    const struct phys_const* restrict phys_const,
    const struct gpart* restrict g) {

  /* initilize time step to disc dynamical time */
  const float dt_dyn = potential->dynamical_time;
  float dt = dt_dyn;

  /* absolute value of height above disc */
  const float dz = fabs(g->x[2] - potential->z_disc);

  /* vertical cceleration */
88
  const float z_accel = 2.f * M_PI * phys_const->const_newton_G *
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
                        potential->surface_density *
                        tanh(dz / potential->scale_height);

  /* demand that dt * velocity <  fraction of scale height of disc */
  float dt1 = FLT_MAX;
  if (fabs(g->v_full[2]) > 0) {
    dt1 = potential->scale_height / fabs(g->v_full[2]);
    if (dt1 < dt) dt = dt1;
  }

  /* demand that dt^2 * acceleration < fraction of scale height of disc */
  float dt2 = FLT_MAX;
  if (fabs(z_accel) > 0) {
    dt2 = potential->scale_height / fabs(z_accel);
    if (dt2 < dt * dt) dt = sqrt(dt2);
  }

  /* demand that dt^3 jerk < fraction of scale height of disc */
  float dt3 = FLT_MAX;
  if (abs(g->v_full[2]) > 0) {
    const float dz_accel_over_dt =
110
        2.f * M_PI * phys_const->const_newton_G * potential->surface_density /
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
        potential->scale_height / cosh(dz / potential->scale_height) /
        cosh(dz / potential->scale_height) * fabs(g->v_full[2]);

    dt3 = potential->scale_height / fabs(dz_accel_over_dt);
    if (dt3 < dt * dt * dt) dt = pow(dt3, 1. / 3.);
  }

  return potential->timestep_mult * dt;
}

/**
 * @brief Computes the gravitational acceleration along z due to a hydrostatic
 * disc
 *
 * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
 *
 * @param time The current time in internal units.
 * @param potential The properties of the potential.
 * @param phys_const The physical constants in internal units.
 * @param g Pointer to the g-particle data.
 */
__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
    double time, const struct external_potential* restrict potential,
    const struct phys_const* restrict phys_const, struct gpart* restrict g) {

  const float dz = g->x[2] - potential->z_disc;
  const float t_dyn = potential->dynamical_time;

139
  float reduction_factor = 1.f;
140
141
142
  if (time < potential->growth_time * t_dyn)
    reduction_factor = time / (potential->growth_time * t_dyn);

143
144
  /* Accelerations. Note that they are multiplied by G later on */
  const float z_accel = reduction_factor * 2.f * M_PI *
145
146
147
                        potential->surface_density *
                        tanh(fabs(dz) / potential->scale_height);

148
149
  if (dz > 0) g->a_grav[2] -= z_accel;
  if (dz < 0) g->a_grav[2] += z_accel;
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
}

/**
 * @brief Initialises the external potential properties in the internal system
 * of units.
 *
 * @param parameter_file The parsed parameter file
 * @param phys_const Physical constants in internal units
 * @param us The current internal system of units
 * @param potential The external potential properties to initialize
 */
static INLINE void potential_init_backend(
    const struct swift_params* parameter_file,
    const struct phys_const* phys_const, const struct UnitSystem* us,
    struct external_potential* potential) {

166
167
168
169
  potential->surface_density = parser_get_param_double(
      parameter_file, "DiscPatchPotential:surface_density");
  potential->scale_height = parser_get_param_double(
      parameter_file, "DiscPatchPotential:scale_height");
170
171
  potential->z_disc =
      parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
172
173
174
175
  potential->timestep_mult = parser_get_param_double(
      parameter_file, "DiscPatchPotential:timestep_mult");
  potential->growth_time = parser_get_opt_param_double(
      parameter_file, "DiscPatchPotential:growth_time", 0.);
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
  potential->dynamical_time =
      sqrt(potential->scale_height /
           (phys_const->const_newton_G * potential->surface_density));
}

/**
 * @brief Prints the properties of the external potential to stdout.
 *
 * @param  potential The external potential properties.
 */
static INLINE void potential_print_backend(
    const struct external_potential* potential) {

  message(
      "External potential is 'Disk-patch' with properties surface_density = %e "
      "disc height= %e scale height = %e timestep multiplier = %e.",
      potential->surface_density, potential->z_disc, potential->scale_height,
      potential->timestep_mult);

  if (potential->growth_time > 0.)
196
    message("Disc will grow for %f dynamical times.", potential->growth_time);
197
198
199
}

#endif /* SWIFT_DISC_PATCH_H */