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

/* Includes. */
#include "cell.h"
25
#include "gravity.h"
26
#include "inline.h"
27
28
#include "part.h"

Matthieu Schaller's avatar
Matthieu Schaller committed
29
30
31
32
33
34
35
36
/**
 * @brief Recursively propagate the multipoles down the tree by applying the
 * L2L and L2P kernels.
 *
 * @param r The #runner.
 * @param c The #cell we are working on.
 * @param timer Are we timing this ?
 */
37
38
void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {

39
  /* Some constants */
40
  const struct engine *e = r->e;
41
42

  /* Cell properties */
43
44
  struct gpart *gparts = c->gparts;
  const int gcount = c->gcount;
45

46
  TIMER_TIC;
47

48
49
50
51
#ifdef SWIFT_DEBUG_CHECKS
  if (c->ti_old_multipole != e->ti_current) error("c->multipole not drifted.");
#endif

52
  if (c->split) { /* Node case */
53

54
    /* Add the field-tensor to all the 8 progenitors */
55
56
57
    for (int k = 0; k < 8; ++k) {
      struct cell *cp = c->progeny[k];

58
59
      /* Do we have a progenitor with any active g-particles ? */
      if (cp != NULL && cell_is_active(cp, e)) {
60

61
62
63
64
#ifdef SWIFT_DEBUG_CHECKS
        if (cp->ti_old_multipole != e->ti_current)
          error("cp->multipole not drifted.");
#endif
65
        struct grav_tensor shifted_tensor;
66

67
        /* Shift the field tensor */
68
69
        gravity_L2L(&shifted_tensor, &c->multipole->pot, cp->multipole->CoM,
                    c->multipole->CoM);
70

71
        /* Add it to this level's tensor */
72
        gravity_field_tensors_add(&cp->multipole->pot, &shifted_tensor);
73

74
        /* Recurse */
75
        runner_do_grav_down(r, cp, 0);
76
77
78
      }
    }

79
  } else { /* Leaf case */
80

81
82
    /* Apply accelerations to the particles */
    for (int i = 0; i < gcount; ++i) {
83
84

      /* Get a handle on the gpart */
85
      struct gpart *gp = &gparts[i];
86
87

      /* Update if active */
88
89
90
91
92
93
94
95
      if (gpart_is_active(gp, e)) {

#ifdef SWIFT_DEBUG_CHECKS
        /* Check that particles have been drifted to the current time */
        if (gp->ti_drift != e->ti_current)
          error("gpart not drifted to current time");
#endif

96
        /* Apply the kernel */
97
        gravity_L2P(&c->multipole->pot, c->multipole->CoM, gp);
98
      }
99
    }
100
  }
101
102

  if (timer) TIMER_TOC(timer_dograv_down);
103
104
}

105
106
107
108
109
110
111
112
/**
 * @brief Computes the interaction of the field tensor in a cell with the
 * multipole of another cell.
 *
 * @param r The #runner.
 * @param ci The #cell with field tensor to interact.
 * @param cj The #cell with the multipole.
 */
113
114
void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
                           struct cell *restrict cj) {
115

116
  /* Some constants */
117
  const struct engine *e = r->e;
118
119
120
  const struct space *s = e->s;
  const int periodic = s->periodic;
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
121
  const struct gravity_props *props = e->gravity_properties;
122
123
  // const float a_smooth = e->gravity_properties->a_smooth;
  // const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
124
125
126

  TIMER_TIC;

127
128
129
  /* Anything to do here? */
  if (!cell_is_active(ci, e)) return;

130
131
132
  /* Short-cut to the multipole */
  const struct multipole *multi_j = &cj->multipole->m_pole;

133
#ifdef SWIFT_DEBUG_CHECKS
134
135
  if (ci == cj) error("Interacting a cell with itself using M2L");

136
  if (multi_j->M_000 == 0.f) error("Multipole does not seem to have been set.");
137

138
139
140
  if (ci->ti_old_multipole != e->ti_current)
    error("ci->multipole not drifted.");
#endif
141

142
143
144
145
  /* Do we need to drift the multipole ? */
  if (cj->ti_old_multipole != e->ti_current) cell_drift_multipole(cj, e);

  /* Let's interact at this level */
146
  gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
147
              cj->multipole->CoM, props, periodic, dim);
148
149
150
151

  TIMER_TOC(timer_dopair_grav_mm);
}

152
153
154
155
156
157
158
159
/**
 * @brief Computes the interaction of all the particles in a cell with the
 * multipole of another cell.
 *
 * @param r The #runner.
 * @param ci The #cell with particles to interct.
 * @param cj The #cell with the multipole.
 */
160
161
162
void runner_dopair_grav_pm(const struct runner *r,
                           const struct cell *restrict ci,
                           const struct cell *restrict cj) {
163

164
  error("Function should not be called");
165
166
167
168
169
170
171
172
173
}

/**
 * @brief Computes the interaction of all the particles in a cell with all the
 * particles of another cell.
 *
 * @param r The #runner.
 * @param ci The first #cell.
 * @param cj The other #cell.
174
175
 *
 * @todo Use a local cache for the particles.
176
 */
177
void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
178

179
  /* Some constants */
180
  const struct engine *e = r->e;
181
182
183
184
185
186
187
  const struct space *s = e->s;
  const int periodic = s->periodic;
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
  const float a_smooth = e->gravity_properties->a_smooth;
  const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);

  /* Cell properties */
188
189
190
191
192
193
194
  const int gcount_i = ci->gcount;
  const int gcount_j = cj->gcount;
  struct gpart *restrict gparts_i = ci->gparts;
  struct gpart *restrict gparts_j = cj->gparts;

  TIMER_TIC;

195
  /* Anything to do here? */
196
  if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
197

198
  /* Let's start by drifting things */
199
200
  if (!cell_are_gpart_drifted(ci, e)) cell_drift_gpart(ci, e);
  if (!cell_are_gpart_drifted(cj, e)) cell_drift_gpart(cj, e);
201

202
203
#if ICHECK > 0
  for (int pid = 0; pid < gcount_i; pid++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
204

205
206
    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gp = &gparts_i[pid];
Matthieu Schaller's avatar
Matthieu Schaller committed
207

Matthieu Schaller's avatar
Matthieu Schaller committed
208
209
210
    if (gp->id_or_neg_offset == ICHECK)
      message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
              gp->id_or_neg_offset, cj->loc[0], cj->loc[1], cj->loc[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
211
              cj->width[0], cj->gcount);
212
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
213

214
  for (int pid = 0; pid < gcount_j; pid++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
215

216
217
    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gp = &gparts_j[pid];
Matthieu Schaller's avatar
Matthieu Schaller committed
218

Matthieu Schaller's avatar
Matthieu Schaller committed
219
220
221
    if (gp->id_or_neg_offset == ICHECK)
      message("id=%lld loc=[ %f %f %f ] size= %f count=%d",
              gp->id_or_neg_offset, ci->loc[0], ci->loc[1], ci->loc[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
222
              ci->width[0], ci->gcount);
223
224
  }
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
225

226
227
228
229
230
231
232
233
  /* Get the relative distance between the pairs, wrapping. */
  double shift[3] = {0.0, 0.0, 0.0};
  if (periodic) {
    shift[0] = nearest(cj->loc[0] - ci->loc[0], dim[0]);
    shift[1] = nearest(cj->loc[1] - ci->loc[1], dim[1]);
    shift[2] = nearest(cj->loc[2] - ci->loc[2], dim[2]);
  }

234
235
  /* MATTHIEU: Should we use local DP accumulators ? */

236
  /* Loop over all particles in ci... */
237
238
  if (cell_is_active(ci, e)) {
    for (int pid = 0; pid < gcount_i; pid++) {
239

240
241
      /* Get a hold of the ith part in ci. */
      struct gpart *restrict gpi = &gparts_i[pid];
242

243
      if (!gpart_is_active(gpi, e)) continue;
244

245
246
247
248
      /* Apply boundary condition */
      const double pix[3] = {gpi->x[0] - shift[0], gpi->x[1] - shift[1],
                             gpi->x[2] - shift[2]};

249
250
      /* Loop over every particle in the other cell. */
      for (int pjd = 0; pjd < gcount_j; pjd++) {
251

252
253
        /* Get a hold of the jth part in cj. */
        const struct gpart *restrict gpj = &gparts_j[pjd];
254

255
        /* Compute the pairwise distance. */
256
257
258
        const float dx[3] = {pix[0] - gpj->x[0],   // x
                             pix[1] - gpj->x[1],   // y
                             pix[2] - gpj->x[2]};  // z
259
        const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
260

261
262
263
264
265
266
267
268
269
270
#ifdef SWIFT_DEBUG_CHECKS
        /* Check that particles have been drifted to the current time */
        if (gpi->ti_drift != e->ti_current)
          error("gpi not drifted to current time");
        if (gpj->ti_drift != e->ti_current)
          error("gpj not drifted to current time");
#endif

        /* Interact ! */
        runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
271

272
#ifdef SWIFT_DEBUG_CHECKS
273
        gpi->num_interacted++;
274
#endif
275
      }
276
277
    }
  }
278

279
  /* Loop over all particles in cj... */
280
281
  if (cell_is_active(cj, e)) {
    for (int pjd = 0; pjd < gcount_j; pjd++) {
282

283
284
      /* Get a hold of the ith part in ci. */
      struct gpart *restrict gpj = &gparts_j[pjd];
285

286
      if (!gpart_is_active(gpj, e)) continue;
287

288
289
290
291
      /* Apply boundary condition */
      const double pjx[3] = {gpj->x[0] + shift[0], gpj->x[1] + shift[1],
                             gpj->x[2] + shift[2]};

292
293
      /* Loop over every particle in the other cell. */
      for (int pid = 0; pid < gcount_i; pid++) {
294

295
296
        /* Get a hold of the ith part in ci. */
        const struct gpart *restrict gpi = &gparts_i[pid];
297

298
        /* Compute the pairwise distance. */
299
300
301
        const float dx[3] = {pjx[0] - gpi->x[0],   // x
                             pjx[1] - gpi->x[1],   // y
                             pjx[2] - gpi->x[2]};  // z
302
        const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
303

304
305
306
307
308
309
310
311
312
313
#ifdef SWIFT_DEBUG_CHECKS
        /* Check that particles have been drifted to the current time */
        if (gpi->ti_drift != e->ti_current)
          error("gpi not drifted to current time");
        if (gpj->ti_drift != e->ti_current)
          error("gpj not drifted to current time");
#endif

        /* Interact ! */
        runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
314
315

#ifdef SWIFT_DEBUG_CHECKS
316
        gpj->num_interacted++;
317
#endif
318
      }
319
320
321
    }
  }

322
  TIMER_TOC(timer_dopair_grav_pp);
323
324
}

325
326
327
328
/**
 * @brief Computes the interaction of all the particles in a cell directly
 *
 * @param r The #runner.
Matthieu Schaller's avatar
Matthieu Schaller committed
329
 * @param c The #cell.
330
331
332
 *
 * @todo Use a local cache for the particles.
 */
333
void runner_doself_grav_pp(struct runner *r, struct cell *c) {
334

335
  /* Some constants */
336
  const struct engine *e = r->e;
337
338
  const float a_smooth = e->gravity_properties->a_smooth;
  const float rlr_inv = 1. / (a_smooth * c->super->width[0]);
339

340
341
342
343
  /* Cell properties */
  const int gcount = c->gcount;
  struct gpart *restrict gparts = c->gparts;

344
345
  TIMER_TIC;

Matthieu Schaller's avatar
Matthieu Schaller committed
346
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
347
  if (c->gcount == 0) error("Doing self gravity on an empty cell !");
348
349
#endif

350
  /* Anything to do here? */
351
  if (!cell_is_active(c, e)) return;
352

353
  /* Do we need to start by drifting things ? */
354
  if (!cell_are_gpart_drifted(c, e)) cell_drift_gpart(c, e);
355

356
357
#if ICHECK > 0
  for (int pid = 0; pid < gcount; pid++) {
358

359
360
    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gp = &gparts[pid];
361

Matthieu Schaller's avatar
Matthieu Schaller committed
362
363
    if (gp->id_or_neg_offset == ICHECK)
      message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
364
365
              gp->id_or_neg_offset, c->loc[0], c->loc[1], c->loc[2],
              c->width[0], c->gcount);
366
367
  }
#endif
368

369
370
  /* MATTHIEU: Should we use local DP accumulators ? */

371
372
373
374
375
376
377
378
379
380
381
382
383
  /* Loop over all particles in ci... */
  for (int pid = 0; pid < gcount; pid++) {

    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gpi = &gparts[pid];

    /* Loop over every particle in the other cell. */
    for (int pjd = pid + 1; pjd < gcount; pjd++) {

      /* Get a hold of the jth part in ci. */
      struct gpart *restrict gpj = &gparts[pjd];

      /* Compute the pairwise distance. */
384
385
386
      float dx[3] = {gpi->x[0] - gpj->x[0],   // x
                     gpi->x[1] - gpj->x[1],   // y
                     gpi->x[2] - gpj->x[2]};  // z
387
388
      const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];

389
390
391
392
393
394
395
396
#ifdef SWIFT_DEBUG_CHECKS
      /* Check that particles have been drifted to the current time */
      if (gpi->ti_drift != e->ti_current)
        error("gpi not drifted to current time");
      if (gpj->ti_drift != e->ti_current)
        error("gpj not drifted to current time");
#endif

397
      /* Interact ! */
398
      if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
399

400
        runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
401

402
#ifdef SWIFT_DEBUG_CHECKS
403
404
        gpi->num_interacted++;
        gpj->num_interacted++;
405
406
#endif

407
408
      } else {

409
        if (gpart_is_active(gpi, e)) {
410

411
          runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
412

413
#ifdef SWIFT_DEBUG_CHECKS
414
          gpi->num_interacted++;
415
416
#endif

417
        } else if (gpart_is_active(gpj, e)) {
418
419
420
421

          dx[0] = -dx[0];
          dx[1] = -dx[1];
          dx[2] = -dx[2];
422
          runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
423
424

#ifdef SWIFT_DEBUG_CHECKS
425
          gpj->num_interacted++;
426
#endif
427
428
        }
      }
429
430
431
    }
  }

432
  TIMER_TOC(timer_doself_grav_pp);
433
434
}

Matthieu Schaller's avatar
Matthieu Schaller committed
435
436
437
438
439
440
441
/**
 * @brief Computes the interaction of all the particles in a cell with all the
 * particles of another cell.
 *
 * @param r The #runner.
 * @param ci The first #cell.
 * @param cj The other #cell.
442
 * @param gettimer Are we timing this ?
Matthieu Schaller's avatar
Matthieu Schaller committed
443
444
445
 *
 * @todo Use a local cache for the particles.
 */
446
447
void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
                        int gettimer) {
Matthieu Schaller's avatar
Matthieu Schaller committed
448

449
450
  /* Some constants */
  const struct engine *e = r->e;
451
452
453
  const struct space *s = e->s;
  const int periodic = s->periodic;
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
454
455
456
  const struct gravity_props *props = e->gravity_properties;
  const double theta_crit_inv = props->theta_crit_inv;

Matthieu Schaller's avatar
Matthieu Schaller committed
457
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
458

Matthieu Schaller's avatar
Matthieu Schaller committed
459
460
461
462
  const int gcount_i = ci->gcount;
  const int gcount_j = cj->gcount;

  /* Early abort? */
Matthieu Schaller's avatar
Matthieu Schaller committed
463
464
  if (gcount_i == 0 || gcount_j == 0)
    error("Doing pair gravity on an empty cell !");
Matthieu Schaller's avatar
Matthieu Schaller committed
465
466

  /* Sanity check */
467
  if (ci == cj) error("Pair interaction between a cell and itself.");
468
469
470
471
472

  if (cell_is_active(ci, e) && ci->ti_old_multipole != e->ti_current)
    error("ci->multipole not drifted.");
  if (cell_is_active(cj, e) && cj->ti_old_multipole != e->ti_current)
    error("cj->multipole not drifted.");
473
474
#endif

475
#if ICHECK > 0
Matthieu Schaller's avatar
Matthieu Schaller committed
476
  for (int pid = 0; pid < ci->gcount; pid++) {
477

478
479
    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gp = &ci->gparts[pid];
480

Matthieu Schaller's avatar
Matthieu Schaller committed
481
482
483
    if (gp->id_or_neg_offset == ICHECK)
      message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
              gp->id_or_neg_offset, cj->loc[0], cj->loc[1], cj->loc[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
484
              cj->width[0], cj->gcount);
485
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
486

Matthieu Schaller's avatar
Matthieu Schaller committed
487
  for (int pid = 0; pid < cj->gcount; pid++) {
488

489
490
    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gp = &cj->gparts[pid];
491

Matthieu Schaller's avatar
Matthieu Schaller committed
492
493
494
    if (gp->id_or_neg_offset == ICHECK)
      message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
              gp->id_or_neg_offset, ci->loc[0], ci->loc[1], ci->loc[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
495
              ci->width[0], ci->gcount);
496
497
  }
#endif
498

499
500
  TIMER_TIC;

501
  /* Can we use M-M interactions ? */
502
  if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
503
504
                               periodic, dim)) {

505
    /* MATTHIEU: make a symmetric M-M interaction function ! */
506
507
508
509
510
511
512
513
514
    runner_dopair_grav_mm(r, ci, cj);
    runner_dopair_grav_mm(r, cj, ci);
  }
  /* We have two leaves. Go P-P. */
  else if (!ci->split && !cj->split) {
    runner_dopair_grav_pp(r, ci, cj);
  }
  /* Alright, we'll have to split and recurse. */
  else {
Matthieu Schaller's avatar
Matthieu Schaller committed
515

516
517
518
519
520
521
522
523
524
525
526
527
528
529
    const double ri_max = ci->multipole->r_max;
    const double rj_max = cj->multipole->r_max;

    /* Split the larger of the two cells and start over again */
    if (ri_max > rj_max) {

      /* Can we actually split that interaction ? */
      if (ci->split) {

        /* Loop over ci's children */
        for (int k = 0; k < 8; k++) {
          if (ci->progeny[k] != NULL)
            runner_dopair_grav(r, ci->progeny[k], cj, 0);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
530

531
      } else if (cj->split) {
532
        /* MATTHIEU: This could maybe be replaced by P-M interactions ?  */
533
534

        /* Loop over cj's children */
Matthieu Schaller's avatar
Matthieu Schaller committed
535
        for (int k = 0; k < 8; k++) {
536
537
538
          if (cj->progeny[k] != NULL)
            runner_dopair_grav(r, ci, cj->progeny[k], 0);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
539

540
      } else {
541
        error("Fundamental error in the logic");
542
543
      }
    } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
544

545
546
547
548
549
550
551
552
      /* Can we actually split that interaction ? */
      if (cj->split) {

        /* Loop over cj's children */
        for (int k = 0; k < 8; k++) {
          if (cj->progeny[k] != NULL)
            runner_dopair_grav(r, ci, cj->progeny[k], 0);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
553

554
      } else if (ci->split) {
555
        /* MATTHIEU: This could maybe be replaced by P-M interactions ?  */
Matthieu Schaller's avatar
Matthieu Schaller committed
556

557
558
559
560
        /* Loop over ci's children */
        for (int k = 0; k < 8; k++) {
          if (ci->progeny[k] != NULL)
            runner_dopair_grav(r, ci->progeny[k], cj, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
561
        }
562
563

      } else {
564
        error("Fundamental error in the logic");
Matthieu Schaller's avatar
Matthieu Schaller committed
565
      }
566
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
567
  }
568
569

  if (gettimer) TIMER_TOC(timer_dosub_pair_grav);
Matthieu Schaller's avatar
Matthieu Schaller committed
570
571
}

572
573
574
575
576
577
578
579
580
/**
 * @brief Computes the interaction of all the particles in a cell
 *
 * @param r The #runner.
 * @param c The first #cell.
 * @param gettimer Are we timing this ?
 *
 * @todo Use a local cache for the particles.
 */
581
void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
Matthieu Schaller's avatar
Matthieu Schaller committed
582

Matthieu Schaller's avatar
Matthieu Schaller committed
583
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
584
  /* Early abort? */
Matthieu Schaller's avatar
Matthieu Schaller committed
585
  if (c->gcount == 0) error("Doing self gravity on an empty cell !");
Matthieu Schaller's avatar
Matthieu Schaller committed
586
587
#endif

588
589
  TIMER_TIC;

Matthieu Schaller's avatar
Matthieu Schaller committed
590
591
592
593
594
595
596
  /* If the cell is split, interact each progeny with itself, and with
     each of its siblings. */
  if (c->split) {

    for (int j = 0; j < 8; j++) {
      if (c->progeny[j] != NULL) {

597
        runner_doself_grav(r, c->progeny[j], 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
598
599
600
601

        for (int k = j + 1; k < 8; k++) {
          if (c->progeny[k] != NULL) {

602
            runner_dopair_grav(r, c->progeny[j], c->progeny[k], 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
603
604
605
606
607
608
609
610
611
612
613
          }
        }
      }
    }
  }

  /* If the cell is not split, then just go for it... */
  else {

    runner_doself_grav_pp(r, c);
  }
614
615

  if (gettimer) TIMER_TOC(timer_dosub_self_grav);
Matthieu Schaller's avatar
Matthieu Schaller committed
616
617
}

618
619
void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
                       int timer) {
620
621
622
623

  /* Is this a single cell? */
  if (cj == NULL) {

624
    runner_doself_grav(r, ci, 1);
625
626
627

  } else {

628
    runner_dopair_grav(r, ci, cj, 1);
629
630
631
  }
}

632
633
634
635
636
637
638
639
/**
 * @brief Performs all M-M interactions between a given top-level cell and all
 * the other top-levels that are far enough.
 *
 * @param r The thread #runner.
 * @param ci The #cell of interest.
 * @param timer Are we timing this ?
 */
640
void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
641

642
#if ICHECK > 0
643
  for (int pid = 0; pid < ci->gcount; pid++) {
644

645
646
    /* Get a hold of the ith part in ci. */
    struct gpart *restrict gp = &ci->gparts[pid];
647

Matthieu Schaller's avatar
Matthieu Schaller committed
648
649
650
    if (gp->id_or_neg_offset == ICHECK)
      message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
              gp->id_or_neg_offset, ci->loc[0], ci->loc[1], ci->loc[2],
Matthieu Schaller's avatar
Matthieu Schaller committed
651
              ci->width[0], ci->gcount);
652
  }
653
#endif
654

655
656
  /* Some constants */
  const struct engine *e = r->e;
657
658
659
  const struct space *s = e->s;
  const int periodic = s->periodic;
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
660
661
662
  const struct gravity_props *props = e->gravity_properties;
  const double theta_crit_inv = props->theta_crit_inv;

663
664
  TIMER_TIC;

665
  /* Recover the list of top-level cells */
666
  struct cell *cells = e->s->cells_top;
667
  const int nr_cells = e->s->nr_cells;
668

669
  /* Anything to do here? */
670
  if (!cell_is_active(ci, e)) return;
671

672
673
674
  /* Check multipole has been drifted */
  if (ci->ti_old_multipole != e->ti_current)
    error("Interacting un-drifted multipole");
675

676
677
  /* Loop over all the top-level cells and go for a M-M interaction if
   * well-separated */
678
679
  for (int i = 0; i < nr_cells; ++i) {

680
    /* Handle on the top-level cell */
681
682
    struct cell *cj = &cells[i];

683
684
    /* Avoid stupid cases */
    if (ci == cj || cj->gcount == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
685

686
    /* Check the multipole acceptance criterion */
687
    if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
688
                                 periodic, dim)) {
689

690
691
692
      /* Go for a (non-symmetric) M-M calculation */
      runner_dopair_grav_mm(r, ci, cj);
    }
693
    /* Is the criterion violated now but was OK at the last rebuild ? */
694
695
    else if (gravity_multipole_accept_rebuild(ci->multipole, cj->multipole,
                                              theta_crit_inv, periodic, dim)) {
696
697

      /* Alright, we have to take charge of that pair in a different way. */
Matthieu Schaller's avatar
Matthieu Schaller committed
698
      // MATTHIEU: We should actually open the tree-node here and recurse.
699
700
      runner_dopair_grav_mm(r, ci, cj);
    }
701
  }
702
703

  if (timer) TIMER_TOC(timer_dograv_long_range);
704
}
705

706
#endif /* SWIFT_RUNNER_DOIACT_GRAV_H */