multipole.h 3.73 KB
Newer Older
Pedro Gonnet's avatar
Pedro Gonnet committed
1
2
3
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
5
6
7
8
 * 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.
9
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
10
11
12
13
 * 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.
14
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
15
16
 * 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/>.
17
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
18
 ******************************************************************************/
19
20
#ifndef SWIFT_MULTIPOLE_H
#define SWIFT_MULTIPOLE_H
Pedro Gonnet's avatar
Pedro Gonnet committed
21

22
23
24
/* Includes. */
#include "part.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
25
/* Some constants. */
26
#define multipole_order 1
Pedro Gonnet's avatar
Pedro Gonnet committed
27
28
29
30

/* Multipole struct. */
struct multipole {

31
32
  /* Multipole location. */
  double x[3];
Pedro Gonnet's avatar
Pedro Gonnet committed
33

34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
  /* Acceleration on this multipole. */
  float a[3];

  /* Multipole coefficients. */
  float coeffs[multipole_order * multipole_order];
};

/* Multipole function prototypes. */
static void multipole_iact_mm(struct multipole *ma, struct multipole *mb,
                              double *shift);
void multipole_merge(struct multipole *ma, struct multipole *mb);
void multipole_addpart(struct multipole *m, struct gpart *p);
void multipole_addparts(struct multipole *m, struct gpart *p, int N);
void multipole_init(struct multipole *m, struct gpart *parts, int N);
void multipole_reset(struct multipole *m);
Pedro Gonnet's avatar
Pedro Gonnet committed
49
50
51

#include <math.h>
#include "kernel.h"
52

Pedro Gonnet's avatar
Pedro Gonnet committed
53
54
55
56
57
58
59
60
/**
 * @brief Compute the pairwise interaction between two multipoles.
 *
 * @param ma The first #multipole.
 * @param mb The second #multipole.
 * @param shift The periodicity correction.
 */

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
88
89
90
91
92
93
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
    struct multipole *ma, struct multipole *mb, double *shift) {

  float dx[3], ir, r, r2 = 0.0f, acc;
  int k;

  /* Compute the multipole distance. */
  for (k = 0; k < 3; k++) {
    dx[k] = ma->x[k] - mb->x[k] - shift[k];
    r2 += dx[k] * dx[k];
  }

  /* Compute the normalized distance vector. */
  ir = 1.0f / sqrtf(r2);
  r = r2 * ir;

  /* Evaluate the gravity kernel. */
  kernel_grav_eval(r, &acc);

  /* Scale the acceleration. */
  acc *= const_G * ir * ir * ir;

/* Compute the forces on both multipoles. */
#if multipole_order == 1
  float mma = ma->coeffs[0], mmb = mb->coeffs[0];
  for (k = 0; k < 3; k++) {
    ma->a[k] -= dx[k] * acc * mmb;
    mb->a[k] += dx[k] * acc * mma;
  }
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Pedro Gonnet's avatar
Pedro Gonnet committed
94
95
96
97
98
99
100
101
102

/**
 * @brief Compute the interaction of a multipole on a particle.
 *
 * @param m The #multipole.
 * @param p The #gpart.
 * @param shift The periodicity correction.
 */

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
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
    struct multipole *m, struct gpart *p, double *shift) {

  float dx[3], ir, r, r2 = 0.0f, acc;
  int k;

  /* Compute the multipole distance. */
  for (k = 0; k < 3; k++) {
    dx[k] = m->x[k] - p->x[k] - shift[k];
    r2 += dx[k] * dx[k];
  }

  /* Compute the normalized distance vector. */
  ir = 1.0f / sqrtf(r2);
  r = r2 * ir;

  /* Evaluate the gravity kernel. */
  kernel_grav_eval(r, &acc);

  /* Scale the acceleration. */
  acc *= const_G * ir * ir * ir * m->coeffs[0];

/* Compute the forces on both multipoles. */
#if multipole_order == 1
  for (k = 0; k < 3; k++) p->a[k] += dx[k] * acc;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Pedro Gonnet's avatar
Pedro Gonnet committed
132

133
#endif /* SWIFT_MULTIPOLE_H */