gravity_iact.h 6.18 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
5
 *
6
7
8
9
 * 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.
10
 *
11
12
13
14
 * 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.
15
 *
16
17
 * 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/>.
18
 *
19
 ******************************************************************************/
20
21
#ifndef SWIFT_DEFAULT_GRAVITY_IACT_H
#define SWIFT_DEFAULT_GRAVITY_IACT_H
22

23
/* Includes. */
24
#include "kernel_gravity.h"
25
#include "kernel_long_gravity.h"
26
#include "multipole.h"
27
28

/**
29
30
31
32
33
34
35
36
37
38
 * @brief Computes the intensity of the force at a point generated by a
 * point-mass.
 *
 * The returned quantity needs to be multiplied by the distance vector to obtain
 * the force vector.
 *
 * @param r2 Square of the distance to the point-mass.
 * @param h2 Square of the softening length.
 * @param h_inv Inverse of the softening length.
 * @param h_inv3 Cube of the inverse of the softening length.
39
 * @param mass Mass of the point-mass.
40
 * @param f_ij (return) The force intensity.
41
 * @param pot_ij (return) The potential.
42
 */
43
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
44
45
    float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
    float *pot_ij) {
46

47
48
  /* Should we soften ? */
  if (r2 >= h2) {
49

50
51
52
    /* Get the inverse distance */
    const float r_inv = 1.f / sqrtf(r2);

53
    /* Get Newtonian gravity */
54
    *f_ij = mass * r_inv * r_inv * r_inv;
55
    *pot_ij = -mass * r_inv;
56

57
  } else {
58

59
    const float r = sqrtf(r2);
60
    const float ui = r * h_inv;
61

62
63
64
    float W_f_ij, W_pot_ij;
    kernel_grav_force_eval(ui, &W_f_ij);
    kernel_grav_pot_eval(ui, &W_pot_ij);
65
66

    /* Get softened gravity */
67
68
    *f_ij = mass * h_inv3 * W_f_ij;
    *pot_ij = mass * h_inv * W_pot_ij;
69
  }
70
71
}

72
/**
73
74
75
76
77
78
79
80
81
82
 * @brief Computes the intensity of the force at a point generated by a
 * point-mass truncated for long-distance periodicity.
 *
 * The returned quantity needs to be multiplied by the distance vector to obtain
 * the force vector.
 *
 * @param r2 Square of the distance to the point-mass.
 * @param h2 Square of the softening length.
 * @param h_inv Inverse of the softening length.
 * @param h_inv3 Cube of the inverse of the softening length.
83
 * @param mass Mass of the point-mass.
84
85
 * @param rlr_inv Inverse of the mesh smoothing scale.
 * @param f_ij (return) The force intensity.
86
 * @param pot_ij (return) The potential.
87
 */
88
89
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
    float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
90
    float *f_ij, float *pot_ij) {
91

92
  float r;
93

94
95
  /* Should we soften ? */
  if (r2 >= h2) {
96

97
98
99
100
    /* Get the inverse distance */
    const float r_inv = 1.f / sqrtf(r2);
    r = r2 * r_inv;

101
    /* Get Newtonian gravity */
102
    *f_ij = mass * r_inv * r_inv * r_inv;
103
    *pot_ij = -mass * r_inv;
104
105
106

  } else {

107
108
109
    /* Get the distance */
    r = sqrtf(r2);

110
    const float ui = r * h_inv;
111
    float W_f_ij, W_pot_ij;
112

113
114
    kernel_grav_force_eval(ui, &W_f_ij);
    kernel_grav_pot_eval(ui, &W_pot_ij);
115
116

    /* Get softened gravity */
117
118
    *f_ij = mass * h_inv3 * W_f_ij;
    *pot_ij = mass * h_inv * W_pot_ij;
119
120
  }

121
  /* Get long-range correction */
122
  const float u_lr = r * rlr_inv;
123
124
125
126
127
  float corr_f_lr, corr_pot_lr;
  kernel_long_grav_force_eval(u_lr, &corr_f_lr);
  kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
  *f_ij *= corr_f_lr;
  *pot_ij *= corr_pot_lr;
128
}
129

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
/**
 * @brief Computes the force at a point generated by a multipole.
 *
 * This uses the quadrupole terms only and defaults to the monopole if
 * the code is compiled with low-order gravity only.
 *
 * @param r_x x-component of the distance vector to the multipole.
 * @param r_y y-component of the distance vector to the multipole.
 * @param r_z z-component of the distance vector to the multipole.
 * @param r2 Square of the distance vector to the multipole.
 * @param h The softening length.
 * @param h_inv Inverse of the softening length.
 * @param m The multipole.
 * @param f_x (return) The x-component of the acceleration.
 * @param f_y (return) The y-component of the acceleration.
 * @param f_z (return) The z-component of the acceleration.
146
 * @param pot (return) The potential.
147
148
149
 */
__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
    float r_x, float r_y, float r_z, float r2, float h, float h_inv,
150
151
    const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {

152
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
153
154
155
  float f_ij;
  runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
                           &f_ij, pot);
156
157
158
  *f_x = -f_ij * r_x;
  *f_y = -f_ij * r_y;
  *f_z = -f_ij * r_z;
159
#else
160
161
162
163

  /* Get the inverse distance */
  const float r_inv = 1.f / sqrtf(r2);

164
165
  struct potential_derivatives_M2P d;
  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &d);
166
167

  /* 1st order terms (monopole) */
168
169
170
171
  *f_x = m->M_000 * d.D_100;
  *f_y = m->M_000 * d.D_010;
  *f_z = m->M_000 * d.D_001;
  *pot = -m->M_000 * d.D_000;
172
173

  /* 3rd order terms (quadrupole) */
174
175
176
177
178
179
180
181
182
  *f_x += m->M_200 * d.D_300 + m->M_020 * d.D_120 + m->M_002 * d.D_102;
  *f_y += m->M_200 * d.D_210 + m->M_020 * d.D_030 + m->M_002 * d.D_012;
  *f_z += m->M_200 * d.D_201 + m->M_020 * d.D_021 + m->M_002 * d.D_003;
  *pot -= m->M_200 * d.D_100 + m->M_020 * d.D_020 + m->M_002 * d.D_002;

  *f_x += m->M_110 * d.D_210 + m->M_101 * d.D_201 + m->M_011 * d.D_111;
  *f_y += m->M_110 * d.D_120 + m->M_101 * d.D_111 + m->M_011 * d.D_021;
  *f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
  *pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
183
184
185
186

#endif
}

187
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */