multipole.h 97 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
210
211
212
213
214
} SWIFT_STRUCT_ALIGN;

/**
 * @brief Reset the data of a #multipole.
 *
 * @param m The #multipole.
 */
215
INLINE static void gravity_reset(struct gravity_tensors *m) {
216
217

  /* Just bzero the struct. */
218
219
220
  bzero(m, sizeof(struct gravity_tensors));
}

221
222
223
/**
 * @brief Drifts a #multipole forward in time.
 *
224
225
226
 * This uses a first-order approximation in time. We only move the CoM
 * using the bulk velocity measured at the last rebuild.
 *
227
228
229
 * @param m The #multipole.
 * @param dt The drift time-step.
 */
230
INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
231

232
  /* Motion of the centre of mass */
233
234
235
  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;
236

237
  /* Move the whole thing according to bulk motion */
238
239
240
  m->CoM[0] += dx;
  m->CoM[1] += dy;
  m->CoM[2] += dz;
241

242
#ifdef SWIFT_DEBUG_CHECKS
243
  if (m->m_pole.vel[0] > m->m_pole.max_delta_vel[0] * 1.1)
244
    error("Invalid maximal velocity");
245
  if (m->m_pole.vel[0] < m->m_pole.min_delta_vel[0] * 1.1)
246
    error("Invalid minimal velocity");
247
  if (m->m_pole.vel[1] > m->m_pole.max_delta_vel[1] * 1.1)
248
    error("Invalid maximal velocity");
249
  if (m->m_pole.vel[1] < m->m_pole.min_delta_vel[1] * 1.1)
250
    error("Invalid minimal velocity");
251
  if (m->m_pole.vel[2] > m->m_pole.max_delta_vel[2] * 1.1)
252
    error("Invalid maximal velocity");
253
  if (m->m_pole.vel[2] < m->m_pole.min_delta_vel[2] * 1.1)
254
255
256
257
258
259
    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
260
              m->m_pole.vel[0] - m->m_pole.min_delta_vel[0]);
261
  dv[1] = max(m->m_pole.max_delta_vel[1] - m->m_pole.vel[1],
Matthieu Schaller's avatar
Matthieu Schaller committed
262
              m->m_pole.vel[1] - m->m_pole.min_delta_vel[1]);
263
  dv[2] = max(m->m_pole.max_delta_vel[2] - m->m_pole.vel[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
264
              m->m_pole.vel[2] - m->m_pole.min_delta_vel[2]);
265

Matthieu Schaller's avatar
Matthieu Schaller committed
266
267
  const float max_delta_vel =
      sqrt(dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]);
268
269
  const float x_diff = max_delta_vel * dt;

270
  /* Conservative change in maximal radius containing all gpart */
271
  m->r_max += x_diff;
272
273
}

274
275
276
277
/**
 * @brief Zeroes all the fields of a field tensor
 *
 * @param l The field tensor.
278
 * @param ti_current The current (integer) time (for debugging only).
279
 */
280
281
INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
                                              integertime_t ti_current) {
282

283
  bzero(l, sizeof(struct grav_tensor));
284
285
286
287

#ifdef SWIFT_DEBUG_CHECKS
  l->ti_init = ti_current;
#endif
288
289
}

290
/**
291
 * @brief Adds a field tensor to another one (i.e. does la += lb).
292
293
294
295
 *
 * @param la The gravity tensors to add to.
 * @param lb The gravity tensors to add.
 */
296
297
INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
                                             const struct grav_tensor *lb) {
298
#ifdef SWIFT_DEBUG_CHECKS
299
300
  if (lb->num_interacted == 0) error("Adding tensors that did not interact");
  la->num_interacted += lb->num_interacted;
301
#endif
302

303
304
  la->interacted = 1;

305
  /* Add 0th order terms */
306
  la->F_000 += lb->F_000;
307
308
309

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
310
311
312
  la->F_100 += lb->F_100;
  la->F_010 += lb->F_010;
  la->F_001 += lb->F_001;
313
314
315
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  /* Add 2nd order terms */
316
317
318
319
320
321
  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;
322
323
324
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  /* Add 3rd order terms */
325
326
327
328
329
330
331
332
333
334
  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;
335
#endif
336
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
  /* 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
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
  /* 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"
380
#endif
381
382
383
384
385
386
387
388
389
390
391
392
}

/**
 * @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");
393
  printf("Interacted: %d\n", l->interacted);
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
  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
417
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
418
419
420
421
422
423
424
425
426
427
428
429
430
  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");
431
432
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
433
#endif
434
435
}

436
437
438
439
440
/**
 * @brief Zeroes all the fields of a multipole.
 *
 * @param m The multipole
 */
441
442
443
444
445
INLINE static void gravity_multipole_init(struct multipole *m) {

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

446
447
448
449
450
451
452
/**
 * @brief Prints the content of a #multipole to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param m The #multipole to print.
 */
453
INLINE static void gravity_multipole_print(const struct multipole *m) {
454

455
456
  printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
  printf("-------------------------\n");
457
458
  printf("M_000= %12.5e\n", m->M_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
459
  printf("-------------------------\n");
460
  printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", m->M_100, m->M_010,
461
462
463
         m->M_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
464
  printf("-------------------------\n");
465
  printf("M_200= %12.5e M_020= %12.5e M_002= %12.5e\n", m->M_200, m->M_020,
466
         m->M_002);
467
  printf("M_110= %12.5e M_101= %12.5e M_011= %12.5e\n", m->M_110, m->M_101,
468
469
470
         m->M_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
471
  printf("-------------------------\n");
472
  printf("M_300= %12.5e M_030= %12.5e M_003= %12.5e\n", m->M_300, m->M_030,
473
         m->M_003);
474
  printf("M_210= %12.5e M_201= %12.5e M_120= %12.5e\n", m->M_210, m->M_201,
475
         m->M_120);
476
  printf("M_021= %12.5e M_102= %12.5e M_012= %12.5e\n", m->M_021, m->M_102,
477
         m->M_012);
478
  printf("M_111= %12.5e\n", m->M_111);
479
#endif
480
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
481
482
483
484
485
486
487
488
489
490
491
492
493
  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");
494
495
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
496
#endif
497
498
499
500
501
502
503
504
}

/**
 * @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.
 */
505
506
INLINE static void gravity_multipole_add(struct multipole *ma,
                                         const struct multipole *mb) {
507

508
509
  /* Add 0th order term */
  ma->M_000 += mb->M_000;
510
511
512
513
514
515
516

#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
517
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
518
  /* Add 2nd order terms */
519
520
521
522
523
524
525
526
  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
527
  /* Add 3rd order terms */
528
529
530
531
532
533
534
535
536
537
538
  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
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
#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
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
  /* 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"
583
#endif
584
585
586
587

#ifdef SWIFT_DEBUG_CHECKS
  ma->num_gpart += mb->num_gpart;
#endif
588
589
590
591
592
}

/**
 * @brief Verifies whether two #multipole's are equal or not.
 *
593
594
 * @param ga The first #multipole.
 * @param gb The second #multipole.
595
 * @param tolerance The maximal allowed relative difference for the fields.
596
 * @return 1 if the multipoles are equal, 0 otherwise
597
 */
598
599
INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
                                          const struct gravity_tensors *gb,
600
                                          double tolerance) {
601
602

  /* Check CoM */
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
  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;
  }
618
619
620
621

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

623
624
625
626
  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)*/
627
  if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
628
      (ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
629
      fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
630
631
632
633
          tolerance) {
    message("v[0] different");
    return 0;
  }
634
  if (fabsf(ma->vel[1] + mb->vel[1]) > 1e-10 &&
635
      (ma->vel[1] * ma->vel[1]) > 0.0001 * v2 &&
636
      fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
637
638
639
640
          tolerance) {
    message("v[1] different");
    return 0;
  }
641
  if (fabsf(ma->vel[2] + mb->vel[2]) > 1e-10 &&
642
      (ma->vel[2] * ma->vel[2]) > 0.0001 * v2 &&
643
      fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
644
645
646
647
          tolerance) {
    message("v[2] different");
    return 0;
  }
648

649
650
651
652
653
654
  /* 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) {
655
656
657
    message("M_000 term different");
    return 0;
  }
658
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
659
660
661
662
663
664
665
666
667
  /* 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");
668
669
    return 0;
  }
670
  if (fabsf(ma->M_010 + mb->M_010) > 0.01f * order1_norm &&
671
672
673
674
      fabsf(ma->M_010 - mb->M_010) / fabsf(ma->M_010 + mb->M_010) > tolerance) {
    message("M_010 term different");
    return 0;
  }
675
676
677
  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");
678
679
    return 0;
  }
680
#endif
681
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
  /* 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");
698
699
    return 0;
  }
700
  if (fabsf(ma->M_020 + mb->M_020) > 0.01f * order2_norm &&
701
702
703
704
      fabsf(ma->M_020 - mb->M_020) / fabsf(ma->M_020 + mb->M_020) > tolerance) {
    message("M_020 term different");
    return 0;
  }
705
706
707
  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");
708
709
    return 0;
  }
710
  if (fabsf(ma->M_110 + mb->M_110) > 0.01f * order2_norm &&
711
712
713
714
      fabsf(ma->M_110 - mb->M_110) / fabsf(ma->M_110 + mb->M_110) > tolerance) {
    message("M_110 term different");
    return 0;
  }
715
716
717
  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");
718
719
    return 0;
  }
720
#endif
721
  tolerance *= 10.;
722
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
723
724
725
726
727
728
729
730
731
732
733
734
  /* 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 &&
735
736
737
738
      fabsf(ma->M_003 - mb->M_003) / fabsf(ma->M_003 + mb->M_003) > tolerance) {
    message("M_003 term different");
    return 0;
  }
739
740
741
  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");
742
743
    return 0;
  }
744
  if (fabsf(ma->M_021 + mb->M_021) > 0.01f * order3_norm &&
745
746
747
748
      fabsf(ma->M_021 - mb->M_021) / fabsf(ma->M_021 + mb->M_021) > tolerance) {
    message("M_021 term different");
    return 0;
  }
749
750
751
  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");
752
753
    return 0;
  }
754
755
756
  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");
757
758
    return 0;
  }
759
  if (fabsf(ma->M_111 + mb->M_111) > 0.01f * order3_norm &&
760
761
762
763
      fabsf(ma->M_111 - mb->M_111) / fabsf(ma->M_111 + mb->M_111) > tolerance) {
    message("M_111 term different");
    return 0;
  }
764
765
766
  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");
767
768
    return 0;
  }
769
770
771
  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");
772
773
    return 0;
  }
774
775
776
  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");
777
778
    return 0;
  }
779
780
781
  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");
782
783
784
    return 0;
  }
#endif
785
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
786
787
788
789
790
791
792
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
  /* 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;
  }
875
#endif
876
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
877
878
879
880
881
882
883
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;
  }
  if (fabsf(ma->M_500 + mb->M_500) > 0.01f * order5_norm &&
      fabsf(ma->M_500 - mb->M_500) / fabsf(ma->M_500 + mb->M_500) > tolerance) {
    message("M_500 term different");
    return 0;
  }
#endif
For faster browsing, not all history is shown. View entire blame