testPotentialPair.c 13.6 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
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2015 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/>.
 *
 ******************************************************************************/
#include "../config.h"

/* Some standard headers. */
#include <fenv.h>
23
#include <math.h>
24 25 26 27 28 29 30
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

/* Local headers. */
#include "runner_doiact_grav.h"
31
#include "swift.h"
32 33

const int num_tests = 100;
34
const double eps = 0.1;
35 36 37 38 39 40 41

/**
 * @brief Check that a and b are consistent (up to some relative error)
 *
 * @param a First value
 * @param b Second value
 * @param s String used to identify this check in messages
42 43
 * @param rel_tol Maximal relative error
 * @param limit Minimal value to consider in the tests
44
 */
45 46 47
void check_value_backend(double a, double b, const char *s, double rel_tol,
                         double limit) {
  if (fabs(a - b) / fabs(a + b) > rel_tol && fabs(a - b) > limit)
48
    error("Values are inconsistent: SWIFT:%12.15e true:%12.15e (%s)!", a, b, s);
49
}
50

51 52 53 54
void check_value(double a, double b, const char *s) {
  check_value_backend(a, b, s, 2e-6, 1e-6);
}

55 56
/* Definitions of the potential and force that match
   exactly the theory document */
57
double S(double x) { return exp(x) / (1. + exp(x)); }
58

59 60 61
double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }

double potential(double mass, double r, double H, double rlr) {
62 63 64 65 66

  const double u = r / H;
  const double x = r / rlr;
  double pot;
  if (u > 1.)
67
    pot = -mass / r;
68
  else
69 70 71 72
    pot = -mass *
          (-3. * u * u * u * u * u * u * u + 15. * u * u * u * u * u * u -
           28. * u * u * u * u * u + 21. * u * u * u * u - 7. * u * u + 3.) /
          H;
73

74
  return pot * (2. - 2. * S(2. * x));
75
}
76 77 78 79 80 81 82 83 84

double acceleration(double mass, double r, double H, double rlr) {

  const double u = r / H;
  const double x = r / rlr;
  double acc;
  if (u > 1.)
    acc = -mass / (r * r * r);
  else
85 86 87
    acc = -mass *
          (21. * u * u * u * u * u - 90. * u * u * u * u + 140. * u * u * u -
           84. * u * u + 14.) /
88 89 90 91 92
          (H * H * H);

  return r * acc * (4. * x * S_prime(2 * x) - 2. * S(2. * x) + 2.);
}

93
int main(int argc, char *argv[]) {
94 95 96 97

  /* Initialize CPU frequency, this also starts time. */
  unsigned long long cpufreq = 0;
  clocks_set_cpufreq(cpufreq);
98 99

  /* Initialise a few things to get us going */
100 101 102 103

  /* Non-truncated forces first */
  double rlr = FLT_MAX;

104 105 106 107
  struct engine e;
  e.max_active_bin = num_time_bins;
  e.time = 0.1f;
  e.ti_current = 8;
108
  e.time_base = 1e-10;
109
  e.nodeID = 0;
110

111 112 113 114
  struct space s;
  s.periodic = 0;
  e.s = &s;

115 116 117 118 119
  struct pm_mesh mesh;
  mesh.periodic = 0;
  mesh.dim[0] = 10.;
  mesh.dim[1] = 10.;
  mesh.dim[2] = 10.;
120 121
  mesh.r_s = rlr;
  mesh.r_s_inv = 1. / rlr;
122
  mesh.r_cut_min = 0.;
123
  mesh.r_cut_max = FLT_MAX;
124
  e.mesh = &mesh;
125

126
  struct gravity_props props;
127
  props.theta_crit = 0.;
128 129
  props.epsilon_DM_cur = eps;
  props.epsilon_baryon_cur = eps;
130
  e.gravity_properties = &props;
131

132 133 134 135
  struct runner r;
  bzero(&r, sizeof(struct runner));
  r.e = &e;

136
  /* Init the cache for gravity interaction */
137 138
  gravity_cache_init(&r.ci_gravity_cache, num_tests);
  gravity_cache_init(&r.cj_gravity_cache, num_tests);
139 140 141

  /* Let's create one cell with a massive particle and a bunch of test particles
   */
142 143 144
  struct cell ci, cj;
  bzero(&ci, sizeof(struct cell));
  bzero(&cj, sizeof(struct cell));
145

146
  ci.nodeID = 0;
147 148 149 150 151 152
  ci.width[0] = 1.;
  ci.width[1] = 1.;
  ci.width[2] = 1.;
  ci.loc[0] = 0.;
  ci.loc[1] = 0.;
  ci.loc[2] = 0.;
153 154
  ci.grav.count = 1;
  ci.grav.ti_old_part = 8;
155 156 157
  ci.grav.ti_old_multipole = 8;
  ci.grav.ti_end_min = 8;
  ci.grav.ti_end_max = 8;
158

159
  cj.nodeID = 0;
160 161 162 163 164 165
  cj.width[0] = 1.;
  cj.width[1] = 1.;
  cj.width[2] = 1.;
  cj.loc[0] = 1.;
  cj.loc[1] = 0.;
  cj.loc[2] = 0.;
166 167
  cj.grav.count = num_tests;
  cj.grav.ti_old_part = 8;
168 169 170
  cj.grav.ti_old_multipole = 8;
  cj.grav.ti_end_min = 8;
  cj.grav.ti_end_max = 8;
171 172

  /* Allocate multipoles */
173
  ci.grav.multipole =
174
      (struct gravity_tensors *)malloc(sizeof(struct gravity_tensors));
175
  cj.grav.multipole =
176
      (struct gravity_tensors *)malloc(sizeof(struct gravity_tensors));
177 178
  bzero(ci.grav.multipole, sizeof(struct gravity_tensors));
  bzero(cj.grav.multipole, sizeof(struct gravity_tensors));
179 180

  /* Set the multipoles */
181 182
  ci.grav.multipole->r_max = 0.1;
  cj.grav.multipole->r_max = 0.1;
183

184
  /* Allocate the particles */
185 186
  if (posix_memalign((void **)&ci.grav.parts, gpart_align,
                     ci.grav.count * sizeof(struct gpart)) != 0)
187
    error("Error allocating gparts for cell ci");
188
  bzero(ci.grav.parts, ci.grav.count * sizeof(struct gpart));
189

190 191
  if (posix_memalign((void **)&cj.grav.parts, gpart_align,
                     cj.grav.count * sizeof(struct gpart)) != 0)
192
    error("Error allocating gparts for cell ci");
193
  bzero(cj.grav.parts, cj.grav.count * sizeof(struct gpart));
194

195
  /* Create the mass-less test particles */
196
  for (int n = 0; n < num_tests; ++n) {
197

198
    struct gpart *gp = &cj.grav.parts[n];
199

200
    gp->x[0] = 1. + (n + 1) / ((double)num_tests);
201 202 203 204 205 206
    gp->x[1] = 0.5;
    gp->x[2] = 0.5;
    gp->mass = 0.;
    gp->time_bin = 1;
    gp->type = swift_type_dark_matter;
    gp->id_or_neg_offset = n + 1;
207 208 209
#ifdef MULTI_SOFTENING_GRAVITY
    gp->epsilon = eps;
#endif
210 211
#ifdef SWIFT_DEBUG_CHECKS
    gp->ti_drift = 8;
212
    gp->initialised = 1;
213 214
#endif
  }
215

216 217 218 219 220
  /***********************************************/
  /* Let's start by testing the P-P interactions */
  /***********************************************/

  /* Create the massive particle */
221 222 223 224 225 226 227
  ci.grav.parts[0].x[0] = 0.;
  ci.grav.parts[0].x[1] = 0.5;
  ci.grav.parts[0].x[2] = 0.5;
  ci.grav.parts[0].mass = 1.;
  ci.grav.parts[0].time_bin = 1;
  ci.grav.parts[0].type = swift_type_dark_matter;
  ci.grav.parts[0].id_or_neg_offset = 1;
228 229 230
#ifdef MULTI_SOFTENING_GRAVITY
  ci.grav.parts[0].epsilon = eps;
#endif
231
#ifdef SWIFT_DEBUG_CHECKS
232
  ci.grav.parts[0].ti_drift = 8;
233
  ci.grav.parts[0].initialised = 1;
234 235
#endif

236
  /* Now compute the forces */
237
  runner_dopair_grav_pp(&r, &ci, &cj, 1, 1);
238

239
  /* Verify everything */
240
  for (int n = 0; n < num_tests; ++n) {
241 242
    const struct gpart *gp = &cj.grav.parts[n];
    const struct gpart *gp2 = &ci.grav.parts[0];
243
    const double epsilon = gravity_get_softening(gp, &props);
244

245
#if defined(POTENTIAL_GRAVITY)
246
    double pot_true =
247
        potential(ci.grav.parts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
248 249 250
    check_value(gp->potential, pot_true, "potential");
#endif

251 252
    double acc_true =
        acceleration(ci.grav.parts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
253

254
    /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0],
255 256 257
       gp->a_grav[0], acc_true, gp->potential, pot_true); */

    check_value(gp->a_grav[0], acc_true, "acceleration");
258
  }
259

260 261 262
  message("\n\t\t P-P interactions all good\n");

  /* Reset the accelerations */
263
  for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.grav.parts[n]);
264

265 266 267
  /**********************************/
  /* Test the basic PM interactions */
  /**********************************/
268

269
  /* Set an opening angle that allows P-M interactions */
270
  props.theta_crit = 1.;
271

272
  ci.grav.parts[0].mass = 0.;
273 274 275
  ci.grav.multipole->CoM[0] = 0.;
  ci.grav.multipole->CoM[1] = 0.5;
  ci.grav.multipole->CoM[2] = 0.5;
276

277 278
  bzero(&ci.grav.multipole->m_pole, sizeof(struct multipole));
  bzero(&cj.grav.multipole->m_pole, sizeof(struct multipole));
279 280
  cj.grav.multipole->m_pole.M_000 = 1.;
  cj.grav.multipole->m_pole.max_softening = eps;
281 282

  /* Now compute the forces */
283
  runner_dopair_grav_pp(&r, &ci, &cj, /*symmetric*/ 1, /*allow_mpoles=*/1);
284 285 286

  /* Verify everything */
  for (int n = 0; n < num_tests; ++n) {
287
    const struct gpart *gp = &cj.grav.parts[n];
288
    const struct gravity_tensors *mpole = ci.grav.multipole;
289
    const double epsilon = gravity_get_softening(gp, &props);
290

291
#if defined(POTENTIAL_GRAVITY)
292 293
    double pot_true =
        potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
294 295 296
    check_value(gp->potential, pot_true, "potential");
#endif

297 298
    double acc_true = acceleration(mpole->m_pole.M_000,
                                   gp->x[0] - mpole->CoM[0], epsilon, rlr);
299
    check_value(gp->a_grav[0], acc_true, "acceleration");
300

301 302
    /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] -
     * mpole->CoM[0], gp->a_grav[0], acc_true, gp->potential, pot_true); */
303 304
  }

305 306
  message("\n\t\t basic P-M interactions all good\n");

307 308
#ifndef GADGET2_LONG_RANGE_CORRECTION

309
  /* Reset the accelerations */
310
  for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.grav.parts[n]);
311

312 313 314 315 316 317 318 319 320 321 322
  /***************************************/
  /* Test the truncated PM interactions  */
  /***************************************/
  rlr = 2.;
  mesh.r_s = rlr;
  mesh.r_s_inv = 1. / rlr;
  mesh.periodic = 1;
  s.periodic = 1;
  props.epsilon_cur = FLT_MIN; /* No softening */

  /* Now compute the forces */
323
  runner_dopair_grav_pp(&r, &ci, &cj, 1, 1);
324 325 326

  /* Verify everything */
  for (int n = 0; n < num_tests; ++n) {
327
    const struct gpart *gp = &cj.grav.parts[n];
328
    const struct gravity_tensors *mpole = ci.grav.multipole;
329 330
    const double epsilon = gravity_get_softening(gp, &props);

331
#if defined(POTENTIAL_GRAVITY)
332 333
    double pot_true =
        potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
334 335 336
    check_value(gp->potential, pot_true, "potential");
#endif

337 338 339
    double acc_true = acceleration(mpole->m_pole.M_000,
                                   gp->x[0] - mpole->CoM[0], epsilon, rlr);
    check_value(gp->a_grav[0], acc_true, "acceleration");
340

Matthieu Schaller's avatar
Matthieu Schaller committed
341 342
    /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] -
     * mpole->CoM[0], */
343
    /*         gp->a_grav[0], acc_true, gp->potential, pot_true); */
344 345 346 347
  }

  message("\n\t\t truncated P-M interactions all good\n");

348 349
#endif

350 351 352
  /************************************************/
  /* Test the high-order periodic PM interactions */
  /************************************************/
353 354

  /* Reset the accelerations */
355
  for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.grav.parts[n]);
356

357
#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3
358 359

  /* Let's make ci more interesting */
360 361 362 363
  free(ci.grav.parts);
  ci.grav.count = 8;
  if (posix_memalign((void **)&ci.grav.parts, gpart_align,
                     ci.grav.count * sizeof(struct gpart)) != 0)
364
    error("Error allocating gparts for cell ci");
365
  bzero(ci.grav.parts, ci.grav.count * sizeof(struct gpart));
366 367 368 369

  /* Place particles on a simple cube of side-length 0.2 */
  for (int n = 0; n < 8; ++n) {
    if (n & 1)
370
      ci.grav.parts[n].x[0] = 0.0 - 0.1;
371
    else
372
      ci.grav.parts[n].x[0] = 0.0 + 0.1;
373 374

    if (n & 2)
375
      ci.grav.parts[n].x[1] = 0.5 - 0.1;
376
    else
377
      ci.grav.parts[n].x[1] = 0.5 + 0.1;
378 379

    if (n & 2)
380
      ci.grav.parts[n].x[2] = 0.5 - 0.1;
381
    else
382
      ci.grav.parts[n].x[2] = 0.5 + 0.1;
383

384
    ci.grav.parts[n].mass = 1. / 8.;
385

386 387 388
    ci.grav.parts[n].time_bin = 1;
    ci.grav.parts[n].type = swift_type_dark_matter;
    ci.grav.parts[n].id_or_neg_offset = 1;
389 390 391
#ifdef MULTI_SOFTENING_GRAVITY
    ci.grav.parts[0].epsilon = eps;
#endif
392
#ifdef SWIFT_DEBUG_CHECKS
393
    ci.grav.parts[n].ti_drift = 8;
394
    ci.grav.parts[n].initialised = 1;
395 396 397 398
#endif
  }

  /* Now let's make a multipole out of it. */
399
  gravity_reset(ci.grav.multipole);
400
  gravity_P2M(ci.grav.multipole, ci.grav.parts, ci.grav.count, &props);
401

402
  gravity_multipole_print(&ci.grav.multipole->m_pole);
403 404

  /* Compute the forces */
405
  runner_dopair_grav_pp(&r, &ci, &cj, 1, 1);
406 407 408

  /* Verify everything */
  for (int n = 0; n < num_tests; ++n) {
409
    const struct gpart *gp = &cj.grav.parts[n];
410

411 412 413 414
#if defined(POTENTIAL_GRAVITY)
    double pot_true = 0;
#endif
    double acc_true[3] = {0., 0., 0.};
415 416

    for (int i = 0; i < 8; ++i) {
417
      const struct gpart *gp2 = &ci.grav.parts[i];
418
      const double epsilon = gravity_get_softening(gp, &props);
419 420 421 422 423

      const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1],
                            gp2->x[2] - gp->x[2]};
      const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);

424
#if defined(POTENTIAL_GRAVITY)
425
      pot_true += potential(gp2->mass, d, epsilon, rlr);
426 427
#endif

428 429 430
      acc_true[0] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[0] / d;
      acc_true[1] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[1] / d;
      acc_true[2] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[2] / d;
431 432
    }

433 434 435 436 437
#if defined(POTENTIAL_GRAVITY)
    check_value_backend(gp->potential, pot_true, "potential", 1e-2, 1e-6);
#endif
    check_value_backend(gp->a_grav[0], acc_true[0], "acceleration", 1e-2, 1e-6);

438
    /* const struct gravity_tensors *mpole = ci.grav.multipole; */
439 440 441 442
    /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e", */
    /*         gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true[0],
     * gp->potential, */
    /*         pot_true, acc_true[1], acc_true[2]); */
443 444 445 446 447
  }

  message("\n\t\t high-order P-M interactions all good\n");

#endif
448

449 450
  free(ci.grav.multipole);
  free(cj.grav.multipole);
451 452
  free(ci.grav.parts);
  free(cj.grav.parts);
453 454 455 456 457

  /* Clean up the caches */
  gravity_cache_clean(&r.ci_gravity_cache);
  gravity_cache_clean(&r.cj_gravity_cache);

458 459
  return 0;
}