multipole.h 97.9 KB
Newer Older
Pedro Gonnet's avatar
Pedro Gonnet committed
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *               2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
5
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
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
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
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
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
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
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
19
 ******************************************************************************/
20
21
#ifndef SWIFT_MULTIPOLE_H
#define SWIFT_MULTIPOLE_H
Pedro Gonnet's avatar
Pedro Gonnet committed
22

23
24
25
/* Config parameters. */
#include "../config.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
/* Some standard headers. */
#include <math.h>
28
#include <string.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
29

30
/* Includes. */
31
#include "align.h"
32
#include "const.h"
33
#include "error.h"
34
#include "gravity.h"
35
#include "gravity_derivatives.h"
36
#include "gravity_properties.h"
37
#include "gravity_softened_derivatives.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
38
#include "inline.h"
39
#include "kernel_gravity.h"
40
#include "part.h"
41
#include "periodic.h"
42
#include "vector_power.h"
43

44
45
#define multipole_align 128

46
struct grav_tensor {
47

48
  /* 0th order terms */
49
  float F_000;
50

51
52
53
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
54
  float F_100, F_010, F_001;
55
56
57
58
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
59
60
  float F_200, F_020, F_002;
  float F_110, F_101, F_011;
61
62
63
64
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
65
66
67
68
69
  float F_300, F_030, F_003;
  float F_210, F_201;
  float F_120, F_021;
  float F_102, F_012;
  float F_111;
70
71
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
72
73

  /* 4th order terms */
74
75
76
77
78
79
  float F_400, F_040, F_004;
  float F_310, F_301;
  float F_130, F_031;
  float F_103, F_013;
  float F_220, F_202, F_022;
  float F_211, F_121, F_112;
80
81
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
82
83

  /* 5th order terms */
84
85
86
87
88
89
90
  float F_005, F_014, F_023;
  float F_032, F_041, F_050;
  float F_104, F_113, F_122;
  float F_131, F_140, F_203;
  float F_212, F_221, F_230;
  float F_302, F_311, F_320;
  float F_401, F_410, F_500;
91
92
93
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
94
95
96
#endif
#ifdef SWIFT_DEBUG_CHECKS

97
98
  /* Total number of gpart this field tensor interacted with */
  long long num_interacted;
99

100
101
102
  /* Last time this tensor was zeroed */
  integertime_t ti_init;

103
#endif
104
105
106

  /* Has this tensor received any contribution? */
  char interacted;
107
108
};

109
struct multipole {
110

111
  /*! Bulk velocity */
112
  float vel[3];
113

114
115
116
117
118
119
120
  /*! Maximal velocity along each axis of all #gpart */
  float max_delta_vel[3];

  /*! Minimal velocity along each axis of all #gpart */
  float min_delta_vel[3];

  /* 0th order term */
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
  float M_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
  float M_100, M_010, M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
  float M_200, M_020, M_002;
  float M_110, M_101, M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
  float M_300, M_030, M_003;
  float M_210, M_201;
  float M_120, M_021;
  float M_102, M_012;
  float M_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
144
145
146
147
148
149
150
151
152
153

  /* 4th order terms */
  float M_400, M_040, M_004;
  float M_310, M_301;
  float M_130, M_031;
  float M_103, M_013;
  float M_220, M_202, M_022;
  float M_211, M_121, M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
154
155

  /* 5th order terms */
156
157
158
159
160
161
162
  float M_005, M_014, M_023;
  float M_032, M_041, M_050;
  float M_104, M_113, M_122;
  float M_131, M_140, M_203;
  float M_212, M_221, M_230;
  float M_302, M_311, M_320;
  float M_401, M_410, M_500;
163
164
165
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
166
#endif
167
168
169
170
171
172
173

#ifdef SWIFT_DEBUG_CHECKS

  /* Total number of gpart in this multipole */
  long long num_gpart;

#endif
174
};
175

176
177
178
179
/**
 * @brief The multipole expansion of a mass distribution.
 */
struct gravity_tensors {
180

181
  union {
182

183
184
    /*! Linking pointer for "memory management". */
    struct gravity_tensors *next;
185

186
187
    /*! The actual content */
    struct {
188

189
190
191
      /*! Field tensor for the potential */
      struct grav_tensor pot;

192
193
194
      /*! Multipole mass */
      struct multipole m_pole;

195
196
      /*! Centre of mass of the matter dsitribution */
      double CoM[3];
197

198
199
200
      /*! Centre of mass of the matter dsitribution at the last rebuild */
      double CoM_rebuild[3];

201
202
203
      /*! Upper limit of the CoM<->gpart distance */
      double r_max;

204
205
      /*! Upper limit of the CoM<->gpart distance at the last rebuild */
      double r_max_rebuild;
206
    };
207
  };
208
209
} SWIFT_STRUCT_ALIGN;

210
211
212
213
214
215
216
#ifdef WITH_MPI
/* MPI datatypes for transfers */
extern MPI_Datatype multipole_mpi_type;
extern MPI_Op multipole_mpi_reduce_op;
void multipole_create_mpi_types(void);
#endif

217
218
219
220
221
/**
 * @brief Reset the data of a #multipole.
 *
 * @param m The #multipole.
 */
222
INLINE static void gravity_reset(struct gravity_tensors *m) {
223
224

  /* Just bzero the struct. */
225
226
227
  bzero(m, sizeof(struct gravity_tensors));
}

228
229
230
/**
 * @brief Drifts a #multipole forward in time.
 *
231
232
233
 * This uses a first-order approximation in time. We only move the CoM
 * using the bulk velocity measured at the last rebuild.
 *
234
235
236
 * @param m The #multipole.
 * @param dt The drift time-step.
 */
237
INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
238

239
  /* Motion of the centre of mass */
240
241
242
  const double dx = m->m_pole.vel[0] * dt;
  const double dy = m->m_pole.vel[1] * dt;
  const double dz = m->m_pole.vel[2] * dt;
243

244
  /* Move the whole thing according to bulk motion */
245
246
247
  m->CoM[0] += dx;
  m->CoM[1] += dy;
  m->CoM[2] += dz;
248

249
#ifdef SWIFT_DEBUG_CHECKS
250
  if (m->m_pole.vel[0] > m->m_pole.max_delta_vel[0] * 1.1)
251
    error("Invalid maximal velocity");
252
  if (m->m_pole.vel[0] < m->m_pole.min_delta_vel[0] * 1.1)
253
    error("Invalid minimal velocity");
254
  if (m->m_pole.vel[1] > m->m_pole.max_delta_vel[1] * 1.1)
255
    error("Invalid maximal velocity");
256
  if (m->m_pole.vel[1] < m->m_pole.min_delta_vel[1] * 1.1)
257
    error("Invalid minimal velocity");
258
  if (m->m_pole.vel[2] > m->m_pole.max_delta_vel[2] * 1.1)
259
    error("Invalid maximal velocity");
260
  if (m->m_pole.vel[2] < m->m_pole.min_delta_vel[2] * 1.1)
261
262
263
264
265
266
    error("Invalid minimal velocity");
#endif

  /* Maximal distance covered by any particle */
  float dv[3];
  dv[0] = max(m->m_pole.max_delta_vel[0] - m->m_pole.vel[0],
Matthieu Schaller's avatar
Matthieu Schaller committed
267
              m->m_pole.vel[0] - m->m_pole.min_delta_vel[0]);
268
  dv[1] = max(m->m_pole.max_delta_vel[1] - m->m_pole.vel[1],
Matthieu Schaller's avatar
Matthieu Schaller committed
269
              m->m_pole.vel[1] - m->m_pole.min_delta_vel[1]);
270
  dv[2] = max(m->m_pole.max_delta_vel[2] - m->m_pole.vel[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
271
              m->m_pole.vel[2] - m->m_pole.min_delta_vel[2]);
272

Matthieu Schaller's avatar
Matthieu Schaller committed
273
274
  const float max_delta_vel =
      sqrt(dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]);
275
276
  const float x_diff = max_delta_vel * dt;

277
  /* Conservative change in maximal radius containing all gpart */
278
  m->r_max += x_diff;
279
280
}

281
282
283
284
/**
 * @brief Zeroes all the fields of a field tensor
 *
 * @param l The field tensor.
285
 * @param ti_current The current (integer) time (for debugging only).
286
 */
287
288
INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
                                              integertime_t ti_current) {
289

290
  bzero(l, sizeof(struct grav_tensor));
291
292
293
294

#ifdef SWIFT_DEBUG_CHECKS
  l->ti_init = ti_current;
#endif
295
296
}

297
/**
298
 * @brief Adds a field tensor to another one (i.e. does la += lb).
299
300
301
302
 *
 * @param la The gravity tensors to add to.
 * @param lb The gravity tensors to add.
 */
303
304
INLINE static void gravity_field_tensors_add(
    struct grav_tensor *restrict la, const struct grav_tensor *restrict lb) {
305
#ifdef SWIFT_DEBUG_CHECKS
306
307
  if (lb->num_interacted == 0) error("Adding tensors that did not interact");
  la->num_interacted += lb->num_interacted;
308
#endif
309

310
311
  la->interacted = 1;

312
  /* Add 0th order terms */
313
  la->F_000 += lb->F_000;
314
315
316

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
317
318
319
  la->F_100 += lb->F_100;
  la->F_010 += lb->F_010;
  la->F_001 += lb->F_001;
320
321
322
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  /* Add 2nd order terms */
323
324
325
326
327
328
  la->F_200 += lb->F_200;
  la->F_020 += lb->F_020;
  la->F_002 += lb->F_002;
  la->F_110 += lb->F_110;
  la->F_101 += lb->F_101;
  la->F_011 += lb->F_011;
329
330
331
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  /* Add 3rd order terms */
332
333
334
335
336
337
338
339
340
341
  la->F_300 += lb->F_300;
  la->F_030 += lb->F_030;
  la->F_003 += lb->F_003;
  la->F_210 += lb->F_210;
  la->F_201 += lb->F_201;
  la->F_120 += lb->F_120;
  la->F_021 += lb->F_021;
  la->F_102 += lb->F_102;
  la->F_012 += lb->F_012;
  la->F_111 += lb->F_111;
342
#endif
343
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
  /* Add 4th order terms */
  la->F_400 += lb->F_400;
  la->F_040 += lb->F_040;
  la->F_004 += lb->F_004;
  la->F_310 += lb->F_310;
  la->F_301 += lb->F_301;
  la->F_130 += lb->F_130;
  la->F_031 += lb->F_031;
  la->F_103 += lb->F_103;
  la->F_013 += lb->F_013;
  la->F_220 += lb->F_220;
  la->F_202 += lb->F_202;
  la->F_022 += lb->F_022;
  la->F_211 += lb->F_211;
  la->F_121 += lb->F_121;
  la->F_112 += lb->F_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
  /* 5th order terms */
  la->F_005 += lb->F_005;
  la->F_014 += lb->F_014;
  la->F_023 += lb->F_023;
  la->F_032 += lb->F_032;
  la->F_041 += lb->F_041;
  la->F_050 += lb->F_050;
  la->F_104 += lb->F_104;
  la->F_113 += lb->F_113;
  la->F_122 += lb->F_122;
  la->F_131 += lb->F_131;
  la->F_140 += lb->F_140;
  la->F_203 += lb->F_203;
  la->F_212 += lb->F_212;
  la->F_221 += lb->F_221;
  la->F_230 += lb->F_230;
  la->F_302 += lb->F_302;
  la->F_311 += lb->F_311;
  la->F_320 += lb->F_320;
  la->F_401 += lb->F_401;
  la->F_410 += lb->F_410;
  la->F_500 += lb->F_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
387
#endif
388
389
390
391
392
393
394
395
396
397
398
399
}

/**
 * @brief Prints the content of a #grav_tensor to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param l The #grav_tensor to print.
 */
INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {

  printf("-------------------------\n");
400
  printf("Interacted: %d\n", l->interacted);
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
  printf("F_000= %12.5e\n", l->F_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  printf("-------------------------\n");
  printf("F_100= %12.5e F_010= %12.5e F_001= %12.5e\n", l->F_100, l->F_010,
         l->F_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  printf("-------------------------\n");
  printf("F_200= %12.5e F_020= %12.5e F_002= %12.5e\n", l->F_200, l->F_020,
         l->F_002);
  printf("F_110= %12.5e F_101= %12.5e F_011= %12.5e\n", l->F_110, l->F_101,
         l->F_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  printf("-------------------------\n");
  printf("F_300= %12.5e F_030= %12.5e F_003= %12.5e\n", l->F_300, l->F_030,
         l->F_003);
  printf("F_210= %12.5e F_201= %12.5e F_120= %12.5e\n", l->F_210, l->F_201,
         l->F_120);
  printf("F_021= %12.5e F_102= %12.5e F_012= %12.5e\n", l->F_021, l->F_102,
         l->F_012);
  printf("F_111= %12.5e\n", l->F_111);
#endif
424
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
425
426
427
428
429
430
431
432
433
434
435
436
437
  printf("-------------------------\n");
  printf("F_400= %12.5e F_040= %12.5e F_004= %12.5e\n", l->F_400, l->F_040,
         l->F_004);
  printf("F_310= %12.5e F_301= %12.5e F_130= %12.5e\n", l->F_310, l->F_301,
         l->F_130);
  printf("F_031= %12.5e F_103= %12.5e F_013= %12.5e\n", l->F_031, l->F_103,
         l->F_013);
  printf("F_220= %12.5e F_202= %12.5e F_022= %12.5e\n", l->F_220, l->F_202,
         l->F_022);
  printf("F_211= %12.5e F_121= %12.5e F_112= %12.5e\n", l->F_211, l->F_121,
         l->F_112);
#endif
  printf("-------------------------\n");
438
439
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
440
#endif
441
442
}

443
444
445
446
447
/**
 * @brief Zeroes all the fields of a multipole.
 *
 * @param m The multipole
 */
448
449
450
451
452
INLINE static void gravity_multipole_init(struct multipole *m) {

  bzero(m, sizeof(struct multipole));
}

453
454
455
456
457
458
459
/**
 * @brief Prints the content of a #multipole to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param m The #multipole to print.
 */
460
INLINE static void gravity_multipole_print(const struct multipole *m) {
461

462
463
  printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
  printf("-------------------------\n");
464
465
  printf("M_000= %12.5e\n", m->M_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
466
  printf("-------------------------\n");
467
  printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", m->M_100, m->M_010,
468
469
470
         m->M_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
471
  printf("-------------------------\n");
472
  printf("M_200= %12.5e M_020= %12.5e M_002= %12.5e\n", m->M_200, m->M_020,
473
         m->M_002);
474
  printf("M_110= %12.5e M_101= %12.5e M_011= %12.5e\n", m->M_110, m->M_101,
475
476
477
         m->M_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
478
  printf("-------------------------\n");
479
  printf("M_300= %12.5e M_030= %12.5e M_003= %12.5e\n", m->M_300, m->M_030,
480
         m->M_003);
481
  printf("M_210= %12.5e M_201= %12.5e M_120= %12.5e\n", m->M_210, m->M_201,
482
         m->M_120);
483
  printf("M_021= %12.5e M_102= %12.5e M_012= %12.5e\n", m->M_021, m->M_102,
484
         m->M_012);
485
  printf("M_111= %12.5e\n", m->M_111);
486
#endif
487
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
488
489
490
491
492
493
494
495
496
497
498
499
500
  printf("-------------------------\n");
  printf("M_400= %12.5e M_040= %12.5e M_004= %12.5e\n", m->M_400, m->M_040,
         m->M_004);
  printf("M_310= %12.5e M_301= %12.5e M_130= %12.5e\n", m->M_310, m->M_301,
         m->M_130);
  printf("M_031= %12.5e M_103= %12.5e M_013= %12.5e\n", m->M_031, m->M_103,
         m->M_013);
  printf("M_220= %12.5e M_202= %12.5e M_022= %12.5e\n", m->M_220, m->M_202,
         m->M_022);
  printf("M_211= %12.5e M_121= %12.5e M_112= %12.5e\n", m->M_211, m->M_121,
         m->M_112);
#endif
  printf("-------------------------\n");
501
502
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
503
#endif
504
505
506
507
508
509
510
511
}

/**
 * @brief Adds a #multipole to another one (i.e. does ma += mb).
 *
 * @param ma The multipole to add to.
 * @param mb The multipole to add.
 */
512
513
INLINE static void gravity_multipole_add(struct multipole *restrict ma,
                                         const struct multipole *restrict mb) {
514

515
516
  /* Add 0th order term */
  ma->M_000 += mb->M_000;
517
518
519
520
521
522
523

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
  ma->M_100 += mb->M_100;
  ma->M_010 += mb->M_010;
  ma->M_001 += mb->M_001;
#endif
524
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
525
  /* Add 2nd order terms */
526
527
528
529
530
531
532
533
  ma->M_200 += mb->M_200;
  ma->M_020 += mb->M_020;
  ma->M_002 += mb->M_002;
  ma->M_110 += mb->M_110;
  ma->M_101 += mb->M_101;
  ma->M_011 += mb->M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
534
  /* Add 3rd order terms */
535
536
537
538
539
540
541
542
543
544
545
  ma->M_300 += mb->M_300;
  ma->M_030 += mb->M_030;
  ma->M_003 += mb->M_003;
  ma->M_210 += mb->M_210;
  ma->M_201 += mb->M_201;
  ma->M_120 += mb->M_120;
  ma->M_021 += mb->M_021;
  ma->M_102 += mb->M_102;
  ma->M_012 += mb->M_012;
  ma->M_111 += mb->M_111;
#endif
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Add 4th order terms */
  ma->M_400 += mb->M_400;
  ma->M_040 += mb->M_040;
  ma->M_004 += mb->M_004;
  ma->M_310 += mb->M_310;
  ma->M_301 += mb->M_301;
  ma->M_130 += mb->M_130;
  ma->M_031 += mb->M_031;
  ma->M_103 += mb->M_103;
  ma->M_013 += mb->M_013;
  ma->M_220 += mb->M_220;
  ma->M_202 += mb->M_202;
  ma->M_022 += mb->M_022;
  ma->M_211 += mb->M_211;
  ma->M_121 += mb->M_121;
  ma->M_112 += mb->M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
  /* 5th order terms */
  ma->M_005 += mb->M_005;
  ma->M_014 += mb->M_014;
  ma->M_023 += mb->M_023;
  ma->M_032 += mb->M_032;
  ma->M_041 += mb->M_041;
  ma->M_050 += mb->M_050;
  ma->M_104 += mb->M_104;
  ma->M_113 += mb->M_113;
  ma->M_122 += mb->M_122;
  ma->M_131 += mb->M_131;
  ma->M_140 += mb->M_140;
  ma->M_203 += mb->M_203;
  ma->M_212 += mb->M_212;
  ma->M_221 += mb->M_221;
  ma->M_230 += mb->M_230;
  ma->M_302 += mb->M_302;
  ma->M_311 += mb->M_311;
  ma->M_320 += mb->M_320;
  ma->M_401 += mb->M_401;
  ma->M_410 += mb->M_410;
  ma->M_500 += mb->M_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
590
#endif
591
592
593
594

#ifdef SWIFT_DEBUG_CHECKS
  ma->num_gpart += mb->num_gpart;
#endif
595
596
597
598
599
}

/**
 * @brief Verifies whether two #multipole's are equal or not.
 *
600
601
 * @param ga The first #multipole.
 * @param gb The second #multipole.
602
 * @param tolerance The maximal allowed relative difference for the fields.
603
 * @return 1 if the multipoles are equal, 0 otherwise
604
 */
605
606
INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
                                          const struct gravity_tensors *gb,
607
                                          double tolerance) {
608
609

  /* Check CoM */
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
  if (fabs(ga->CoM[0] - gb->CoM[0]) / fabs(ga->CoM[0] + gb->CoM[0]) >
      tolerance) {
    message("CoM[0] different");
    return 0;
  }
  if (fabs(ga->CoM[1] - gb->CoM[1]) / fabs(ga->CoM[1] + gb->CoM[1]) >
      tolerance) {
    message("CoM[1] different");
    return 0;
  }
  if (fabs(ga->CoM[2] - gb->CoM[2]) / fabs(ga->CoM[2] + gb->CoM[2]) >
      tolerance) {
    message("CoM[2] different");
    return 0;
  }
625
626
627
628

  /* Helper pointers */
  const struct multipole *ma = &ga->m_pole;
  const struct multipole *mb = &gb->m_pole;
629

630
631
632
633
  const double v2 = ma->vel[0] * ma->vel[0] + ma->vel[1] * ma->vel[1] +
                    ma->vel[2] * ma->vel[2];

  /* Check bulk velocity (if non-zero and component > 1% of norm)*/
634
  if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
635
      (ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
636
      fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
637
638
639
640
          tolerance) {
    message("v[0] different");
    return 0;
  }
641
  if (fabsf(ma->vel[1] + mb->vel[1]) > 1e-10 &&
642
      (ma->vel[1] * ma->vel[1]) > 0.0001 * v2 &&
643
      fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
644
645
646
647
          tolerance) {
    message("v[1] different");
    return 0;
  }
648
  if (fabsf(ma->vel[2] + mb->vel[2]) > 1e-10 &&
649
      (ma->vel[2] * ma->vel[2]) > 0.0001 * v2 &&
650
      fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
651
652
653
654
          tolerance) {
    message("v[2] different");
    return 0;
  }
655

656
657
658
659
660
661
  /* Manhattan Norm of 0th order terms */
  const float order0_norm = fabsf(ma->M_000) + fabsf(mb->M_000);

  /* Compare 0th order terms above 1% of norm */
  if (fabsf(ma->M_000 + mb->M_000) > 0.01f * order0_norm &&
      fabsf(ma->M_000 - mb->M_000) / fabsf(ma->M_000 + mb->M_000) > tolerance) {
662
663
664
    message("M_000 term different");
    return 0;
  }
665
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
666
667
668
669
670
671
672
673
674
  /* Manhattan Norm of 1st order terms */
  const float order1_norm = fabsf(ma->M_001) + fabsf(mb->M_001) +
                            fabsf(ma->M_010) + fabsf(mb->M_010) +
                            fabsf(ma->M_100) + fabsf(mb->M_100);

  /* Compare 1st order terms above 1% of norm */
  if (fabsf(ma->M_001 + mb->M_001) > 0.01f * order1_norm &&
      fabsf(ma->M_001 - mb->M_001) / fabsf(ma->M_001 + mb->M_001) > tolerance) {
    message("M_001 term different");
675
676
    return 0;
  }
677
  if (fabsf(ma->M_010 + mb->M_010) > 0.01f * order1_norm &&
678
679
680
681
      fabsf(ma->M_010 - mb->M_010) / fabsf(ma->M_010 + mb->M_010) > tolerance) {
    message("M_010 term different");
    return 0;
  }
682
683
684
  if (fabsf(ma->M_100 + mb->M_100) > 0.01f * order1_norm &&
      fabsf(ma->M_100 - mb->M_100) / fabsf(ma->M_100 + mb->M_100) > tolerance) {
    message("M_100 term different");
685
686
    return 0;
  }
687
#endif
688
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
  /* Manhattan Norm of 2nd order terms */
  const float order2_norm =
      fabsf(ma->M_002) + fabsf(mb->M_002) + fabsf(ma->M_011) +
      fabsf(mb->M_011) + fabsf(ma->M_020) + fabsf(mb->M_020) +
      fabsf(ma->M_101) + fabsf(mb->M_101) + fabsf(ma->M_110) +
      fabsf(mb->M_110) + fabsf(ma->M_200) + fabsf(mb->M_200);

  /* Compare 2nd order terms above 1% of norm */
  if (fabsf(ma->M_002 + mb->M_002) > 0.01f * order2_norm &&
      fabsf(ma->M_002 - mb->M_002) / fabsf(ma->M_002 + mb->M_002) > tolerance) {
    message("M_002 term different");
    return 0;
  }
  if (fabsf(ma->M_011 + mb->M_011) > 0.01f * order2_norm &&
      fabsf(ma->M_011 - mb->M_011) / fabsf(ma->M_011 + mb->M_011) > tolerance) {
    message("M_011 term different");
705
706
    return 0;
  }
707
  if (fabsf(ma->M_020 + mb->M_020) > 0.01f * order2_norm &&
708
709
710
711
      fabsf(ma->M_020 - mb->M_020) / fabsf(ma->M_020 + mb->M_020) > tolerance) {
    message("M_020 term different");
    return 0;
  }
712
713
714
  if (fabsf(ma->M_101 + mb->M_101) > 0.01f * order2_norm &&
      fabsf(ma->M_101 - mb->M_101) / fabsf(ma->M_101 + mb->M_101) > tolerance) {
    message("M_101 term different");
715
716
    return 0;
  }
717
  if (fabsf(ma->M_110 + mb->M_110) > 0.01f * order2_norm &&
718
719
720
721
      fabsf(ma->M_110 - mb->M_110) / fabsf(ma->M_110 + mb->M_110) > tolerance) {
    message("M_110 term different");
    return 0;
  }
722
723
724
  if (fabsf(ma->M_200 + mb->M_200) > 0.01f * order2_norm &&
      fabsf(ma->M_200 - mb->M_200) / fabsf(ma->M_200 + mb->M_200) > tolerance) {
    message("M_200 term different");
725
726
    return 0;
  }
727
#endif
728
  tolerance *= 10.;
729
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
730
731
732
733
734
735
736
737
738
739
740
741
  /* Manhattan Norm of 3rd order terms */
  const float order3_norm =
      fabsf(ma->M_003) + fabsf(mb->M_003) + fabsf(ma->M_012) +
      fabsf(mb->M_012) + fabsf(ma->M_021) + fabsf(mb->M_021) +
      fabsf(ma->M_030) + fabsf(mb->M_030) + fabsf(ma->M_102) +
      fabsf(mb->M_102) + fabsf(ma->M_111) + fabsf(mb->M_111) +
      fabsf(ma->M_120) + fabsf(mb->M_120) + fabsf(ma->M_201) +
      fabsf(mb->M_201) + fabsf(ma->M_210) + fabsf(mb->M_210) +
      fabsf(ma->M_300) + fabsf(mb->M_300);

  /* Compare 3rd order terms above 1% of norm */
  if (fabsf(ma->M_003 + mb->M_003) > 0.01f * order3_norm &&
742
743
744
745
      fabsf(ma->M_003 - mb->M_003) / fabsf(ma->M_003 + mb->M_003) > tolerance) {
    message("M_003 term different");
    return 0;
  }
746
747
748
  if (fabsf(ma->M_012 + mb->M_012) > 0.01f * order3_norm &&
      fabsf(ma->M_012 - mb->M_012) / fabsf(ma->M_012 + mb->M_012) > tolerance) {
    message("M_012 term different");
749
750
    return 0;
  }
751
  if (fabsf(ma->M_021 + mb->M_021) > 0.01f * order3_norm &&
752
753
754
755
      fabsf(ma->M_021 - mb->M_021) / fabsf(ma->M_021 + mb->M_021) > tolerance) {
    message("M_021 term different");
    return 0;
  }
756
757
758
  if (fabsf(ma->M_030 + mb->M_030) > 0.01f * order3_norm &&
      fabsf(ma->M_030 - mb->M_030) / fabsf(ma->M_030 + mb->M_030) > tolerance) {
    message("M_030 term different");
759
760
    return 0;
  }
761
762
763
  if (fabsf(ma->M_102 + mb->M_102) > 0.01f * order3_norm &&
      fabsf(ma->M_102 - mb->M_102) / fabsf(ma->M_102 + mb->M_102) > tolerance) {
    message("M_102 term different");
764
765
    return 0;
  }
766
  if (fabsf(ma->M_111 + mb->M_111) > 0.01f * order3_norm &&
767
768
769
770
      fabsf(ma->M_111 - mb->M_111) / fabsf(ma->M_111 + mb->M_111) > tolerance) {
    message("M_111 term different");
    return 0;
  }
771
772
773
  if (fabsf(ma->M_120 + mb->M_120) > 0.01f * order3_norm &&
      fabsf(ma->M_120 - mb->M_120) / fabsf(ma->M_120 + mb->M_120) > tolerance) {
    message("M_120 term different");
774
775
    return 0;
  }
776
777
778
  if (fabsf(ma->M_201 + mb->M_201) > 0.01f * order3_norm &&
      fabsf(ma->M_201 - mb->M_201) / fabsf(ma->M_201 + mb->M_201) > tolerance) {
    message("M_201 term different");
779
780
    return 0;
  }
781
782
783
  if (fabsf(ma->M_210 + mb->M_210) > 0.01f * order3_norm &&
      fabsf(ma->M_210 - mb->M_210) / fabsf(ma->M_210 + mb->M_210) > tolerance) {
    message("M_210 term different");
784
785
    return 0;
  }
786
787
788
  if (fabsf(ma->M_300 + mb->M_300) > 0.01f * order3_norm &&
      fabsf(ma->M_300 - mb->M_300) / fabsf(ma->M_300 + mb->M_300) > tolerance) {
    message("M_300 term different");
789
790
791
    return 0;
  }
#endif
792
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
  /* Manhattan Norm of 4th order terms */
  const float order4_norm =
      fabsf(ma->M_004) + fabsf(mb->M_004) + fabsf(ma->M_013) +
      fabsf(mb->M_013) + fabsf(ma->M_022) + fabsf(mb->M_022) +
      fabsf(ma->M_031) + fabsf(mb->M_031) + fabsf(ma->M_040) +
      fabsf(mb->M_040) + fabsf(ma->M_103) + fabsf(mb->M_103) +
      fabsf(ma->M_112) + fabsf(mb->M_112) + fabsf(ma->M_121) +
      fabsf(mb->M_121) + fabsf(ma->M_130) + fabsf(mb->M_130) +
      fabsf(ma->M_202) + fabsf(mb->M_202) + fabsf(ma->M_211) +
      fabsf(mb->M_211) + fabsf(ma->M_220) + fabsf(mb->M_220) +
      fabsf(ma->M_301) + fabsf(mb->M_301) + fabsf(ma->M_310) +
      fabsf(mb->M_310) + fabsf(ma->M_400) + fabsf(mb->M_400);

  /* Compare 4th order terms above 1% of norm */
  if (fabsf(ma->M_004 + mb->M_004) > 0.01f * order4_norm &&
      fabsf(ma->M_004 - mb->M_004) / fabsf(ma->M_004 + mb->M_004) > tolerance) {
    message("M_004 term different");
    return 0;
  }
  if (fabsf(ma->M_013 + mb->M_013) > 0.01f * order4_norm &&
      fabsf(ma->M_013 - mb->M_013) / fabsf(ma->M_013 + mb->M_013) > tolerance) {
    message("M_013 term different");
    return 0;
  }
  if (fabsf(ma->M_022 + mb->M_022) > 0.01f * order4_norm &&
      fabsf(ma->M_022 - mb->M_022) / fabsf(ma->M_022 + mb->M_022) > tolerance) {
    message("M_022 term different");
    return 0;
  }
  if (fabsf(ma->M_031 + mb->M_031) > 0.01f * order4_norm &&
      fabsf(ma->M_031 - mb->M_031) / fabsf(ma->M_031 + mb->M_031) > tolerance) {
    message("M_031 term different");
    return 0;
  }
  if (fabsf(ma->M_040 + mb->M_040) > 0.01f * order4_norm &&
      fabsf(ma->M_040 - mb->M_040) / fabsf(ma->M_040 + mb->M_040) > tolerance) {
    message("M_040 term different");
    return 0;
  }
  if (fabsf(ma->M_103 + mb->M_103) > 0.01f * order4_norm &&
      fabsf(ma->M_103 - mb->M_103) / fabsf(ma->M_103 + mb->M_103) > tolerance) {
    message("M_103 term different");
    return 0;
  }
  if (fabsf(ma->M_112 + mb->M_112) > 0.01f * order4_norm &&
      fabsf(ma->M_112 - mb->M_112) / fabsf(ma->M_112 + mb->M_112) > tolerance) {
    message("M_112 term different");
    return 0;
  }
  if (fabsf(ma->M_121 + mb->M_121) > 0.01f * order4_norm &&
      fabsf(ma->M_121 - mb->M_121) / fabsf(ma->M_121 + mb->M_121) > tolerance) {
    message("M_121 term different");
    return 0;
  }
  if (fabsf(ma->M_130 + mb->M_130) > 0.01f * order4_norm &&
      fabsf(ma->M_130 - mb->M_130) / fabsf(ma->M_130 + mb->M_130) > tolerance) {
    message("M_130 term different");
    return 0;
  }
  if (fabsf(ma->M_202 + mb->M_202) > 0.01f * order4_norm &&
      fabsf(ma->M_202 - mb->M_202) / fabsf(ma->M_202 + mb->M_202) > tolerance) {
    message("M_202 term different");
    return 0;
  }
  if (fabsf(ma->M_211 + mb->M_211) > 0.01f * order4_norm &&
      fabsf(ma->M_211 - mb->M_211) / fabsf(ma->M_211 + mb->M_211) > tolerance) {
    message("M_211 term different");
    return 0;
  }
  if (fabsf(ma->M_220 + mb->M_220) > 0.01f * order4_norm &&
      fabsf(ma->M_220 - mb->M_220) / fabsf(ma->M_220 + mb->M_220) > tolerance) {
    message("M_220 term different");
    return 0;
  }
  if (fabsf(ma->M_301 + mb->M_301) > 0.01f * order4_norm &&
      fabsf(ma->M_301 - mb->M_301) / fabsf(ma->M_301 + mb->M_301) > tolerance) {
    message("M_301 term different");
    return 0;
  }
  if (fabsf(ma->M_310 + mb->M_310) > 0.01f * order4_norm &&
      fabsf(ma->M_310 - mb->M_310) / fabsf(ma->M_310 + mb->M_310) > tolerance) {
    message("M_310 term different");
    return 0;
  }
  if (fabsf(ma->M_400 + mb->M_400) > 0.01f * order4_norm &&
      fabsf(ma->M_400 - mb->M_400) / fabsf(ma->M_400 + mb->M_400) > tolerance) {
    message("M_400 term different");
    return 0;
  }
882
#endif
883
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
  /* Manhattan Norm of 5th order terms */
  const float order5_norm =
      fabsf(ma->M_005) + fabsf(mb->M_005) + fabsf(ma->M_014) +
      fabsf(mb->M_014) + fabsf(ma->M_023) + fabsf(mb->M_023) +
      fabsf(ma->M_032) + fabsf(mb->M_032) + fabsf(ma->M_041) +
      fabsf(mb->M_041) + fabsf(ma->M_050) + fabsf(mb->M_050) +
      fabsf(ma->M_104) + fabsf(mb->M_104) + fabsf(ma->M_113) +
      fabsf(mb->M_113) + fabsf(ma->M_122) + fabsf(mb->M_122) +
      fabsf(ma->M_131) + fabsf(mb->M_131) + fabsf(ma->M_140) +
      fabsf(mb->M_140) + fabsf(ma->M_203) + fabsf(mb->M_203) +
      fabsf(ma->M_212) + fabsf(mb->M_212) + fabsf(ma->M_221) +
      fabsf(mb->M_221) + fabsf(ma->M_230) + fabsf(mb->M_230) +
      fabsf(ma->M_302) + fabsf(mb->M_302) + fabsf(ma->M_311) +
      fabsf(mb->M_311) + fabsf(ma->M_320) + fabsf(mb->M_320) +
      fabsf(ma->M_401) + fabsf(mb->M_401) + fabsf(ma->M_410) +
      fabsf(mb->M_410) + fabsf(ma->M_500) + fabsf(mb->M_500);

  /* Compare 5th order terms above 1% of norm */
  if (fabsf(ma->M_005 + mb->M_005) > 0.01f * order5_norm &&
      fabsf(ma->M_005 - mb->M_005) / fabsf(ma->M_005 + mb->M_005) > tolerance) {
    message("M_005 term different");
    return 0;
  }
  if (fabsf(ma->M_014 + mb->M_014) > 0.01f * order5_norm &&
      fabsf(ma->M_014 - mb->M_014) / fabsf(ma->M_014 + mb->M_014) > tolerance) {
    message("M_014 term different");
    return 0;
  }
  if (fabsf(ma->M_023 + mb->M_023) > 0.01f * order5_norm &&
      fabsf(ma->M_023 - mb->M_023) / fabsf(ma->M_023 + mb->M_023) > tolerance) {
    message("M_023 term different");
    return 0;
  }
  if (fabsf(ma->M_032 + mb->M_032) > 0.01f * order5_norm &&
      fabsf(ma->M_032 - mb->M_032) / fabsf(ma->M_032 + mb->M_032) > tolerance) {
    message("M_032 term different");
    return 0;
  }
  if (fabsf(ma->M_041 + mb->M_041) > 0.01f * order5_norm &&
      fabsf(ma->M_041 - mb->M_041) / fabsf(ma->M_041 + mb->M_041) > tolerance) {
    message("M_041 term different");
    return 0;
  }
  if (fabsf(ma->M_050 + mb->M_050) > 0.01f * order5_norm &&
      fabsf(ma->M_050 - mb->M_050) / fabsf(ma->M_050 + mb->M_050) > tolerance) {
    message("M_050 term different");
    return 0;
  }
  if (fabsf(ma->M_104 + mb->M_104) > 0.01f * order5_norm &&
      fabsf(ma->M_104 - mb->M_104) / fabsf(ma->M_104 + mb->M_104) > tolerance) {
    message("M_104 term different");
    return 0;
  }
  if (fabsf(ma->M_113 + mb->M_113) > 0.01f * order5_norm &&
      fabsf(ma->M_113 - mb->M_113) / fabsf(ma->M_113 + mb->M_113) > tolerance) {
    message("M_113 term different");
    return 0;
  }
  if (fabsf(ma->M_122 + mb->M_122) > 0.01f * order5_norm &&
      fabsf(ma->M_122 - mb->M_122) / fabsf(ma->M_122 + mb->M_122) > tolerance) {
    message("M_122 term different");
    return 0;
  }
  if (fabsf(ma->M_131 + mb->M_131) > 0.01f * order5_norm &&
      fabsf(ma->M_131 - mb->M_131) / fabsf(ma->M_131 + mb->M_131) > tolerance) {
    message("M_131 term different");
    return 0;
  }
  if (fabsf(ma->M_140 + mb->M_140) > 0.01f * order5_norm &&
      fabsf(ma->M_140 - mb->M_140) / fabsf(ma->M_140 + mb->M_140) > tolerance) {
    message("M_140 term different");
    return 0;
  }
  if (fabsf(ma->M_203 + mb->M_203) > 0.01f * order5_norm &&
      fabsf(ma->M_203 - mb->M_203) / fabsf(ma->M_203 + mb->M_203) > tolerance) {
    message("M_203 term different");
    return 0;
  }
  if (fabsf(ma->M_212 + mb->M_212) > 0.01f * order5_norm &&
      fabsf(ma->M_212 - mb->M_212) / fabsf(ma->M_212 + mb->M_212) > tolerance) {
    message("M_212 term different");
    return 0;
  }
  if (fabsf(ma->M_221 + mb->M_221) > 0.01f * order5_norm &&
      fabsf(ma->M_221 - mb->M_221) / fabsf(ma->M_221 + mb->M_221) > tolerance) {
    message("M_221 term different");
    return 0;
  }
  if (fabsf(ma->M_230 + mb->M_230) > 0.01f * order5_norm &&
      fabsf(ma->M_230 - mb->M_230) / fabsf(ma->M_230 + mb->M_230) > tolerance) {
    message("M_230 term different");
    return 0;
  }
  if (fabsf(ma->M_302 + mb->M_302) > 0.01f * order5_norm &&
      fabsf(ma->M_302 - mb->M_302) / fabsf(ma->M_302 + mb->M_302) > tolerance) {
    message("M_302 term different");
    return 0;
  }
  if (fabsf(ma->M_311 + mb->M_311) > 0.01f * order5_norm &&
      fabsf(ma->M_311 - mb->M_311) / fabsf(ma->M_311 + mb->M_311) > tolerance) {
    message("M_311 term different");
    return 0;
  }
  if (fabsf(ma->M_320 + mb->M_320) > 0.01f * order5_norm &&
      fabsf(ma->M_320 - mb->M_320) / fabsf(ma->M_320 + mb->M_320) > tolerance) {
    message("M_320 term different");
    return 0;
  }
  if (fabsf(ma->M_401 + mb->M_401) > 0.01f * order5_norm &&
      fabsf(ma->M_401 - mb->M_401) / fabsf(ma->M_401 + mb->M_401) > tolerance) {
    message("M_401 term different");
    return 0;
  }
  if (fabsf(ma->M_410 + mb->M_410) > 0.01f * order5_norm &&
      fabsf(ma->M_410 - mb->M_410) / fabsf(ma->M_410 + mb->M_410) > tolerance) {
    message("M_410 term different");
    return 0;
For faster browsing, not all history is shown. View entire blame