multipole.c 4.16 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
21
22
23
24
25
26
27
28
29
30
31
32
 ******************************************************************************/

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

/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <limits.h>

/* MPI headers. */
#ifdef WITH_MPI
33
#include <mpi.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#endif

/* Local headers. */
#include "error.h"
#include "const.h"
#include "cycle.h"
#include "atomic.h"
#include "lock.h"
#include "space.h"
#include "multipole.h"
#include "cell.h"

/**
 * @brief Merge two multipoles.
 *
 * @param ma The #multipole which will contain the merged result.
 * @param mb The other #multipole.
 */

53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
void multipole_merge(struct multipole *ma, struct multipole *mb) {

#if multipole_order == 1

  /* Correct the position. */
  float mma = ma->coeffs[0], mmb = mb->coeffs[0];
  float w = 1.0f / (mma + mmb);
  for (int k = 0; k < 3; k++) ma->x[k] = (ma->x[k] * mma + mb->x[k] * mmb) * w;

  /* Add the particle to the moments. */
  ma->coeffs[0] = mma + mmb;

#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Pedro Gonnet's avatar
Pedro Gonnet committed
69
70
71
72
73
74
75
76

/**
 * @brief Add a particle to the given multipole.
 *
 * @param m The #multipole.
 * @param p The #gpart.
 */

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
void multipole_addpart(struct multipole *m, struct gpart *p) {

#if multipole_order == 1

  /* Correct the position. */
  float mm = m->coeffs[0], mp = p->mass;
  float w = 1.0f / (mm + mp);
  for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + p->x[k] * mp) * w;

  /* Add the particle to the moments. */
  m->coeffs[0] = mm + mp;

#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Pedro Gonnet's avatar
Pedro Gonnet committed
93
94
95
96
97
98
99
100
101

/**
 * @brief Add a group of particles to the given multipole.
 *
 * @param m The #multipole.
 * @param p The #gpart array.
 * @param N Number of parts to add.
 */

102
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
void multipole_addparts(struct multipole *m, struct gpart *p, int N) {

#if multipole_order == 1

  /* Get the combined mass and positions. */
  double xp[3] = {0.0, 0.0, 0.0};
  float mp = 0.0f, w;
  for (int k = 0; k < N; k++) {
    w = p[k].mass;
    mp += w;
    xp[0] += p[k].x[0] * w;
    xp[1] += p[k].x[1] * w;
    xp[2] += p[k].x[2] * w;
  }

  /* Correct the position. */
  float mm = m->coeffs[0];
  w = 1.0f / (mm + mp);
  for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + xp[k]) * w;

  /* Add the particle to the moments. */
  m->coeffs[0] = mm + mp;

#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Pedro Gonnet's avatar
Pedro Gonnet committed
129
130
131
132
133

/**
 * @brief Init a multipole from a set of particles.
 *
 * @param m The #multipole.
134
 * @param parts The #gpart.
Pedro Gonnet's avatar
Pedro Gonnet committed
135
136
137
 * @param N The number of particles.
 */

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
void multipole_init(struct multipole *m, struct gpart *parts, int N) {

#if multipole_order == 1

  float mass = 0.0f, w;
  double x[3] = {0.0, 0.0, 0.0};
  int k;

  /* Collect the particle data. */
  for (k = 0; k < N; k++) {
    w = parts[k].mass;
    mass += w;
    x[0] += parts[k].x[0] * w;
    x[1] += parts[k].x[1] * w;
    x[2] += parts[k].x[2] * w;
  }

  /* Store the data on the multipole. */
  m->coeffs[0] = mass;
  m->x[0] = x[0] / mass;
  m->x[1] = x[1] / mass;
  m->x[2] = x[2] / mass;

#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Pedro Gonnet's avatar
Pedro Gonnet committed
165
166
167
168
169
170
171

/**
 * @brief Reset the data of a #multipole.
 *
 * @param m The #multipole.
 */

172
173
174
175
176
void multipole_reset(struct multipole *m) {

  /* Just bzero the struct. */
  bzero(m, sizeof(struct multipole));
}