Commit d400531c authored by Matthieu Schaller's avatar Matthieu Schaller Committed by Matthieu Schaller
Browse files

The testPotential unit test now also tests the forces.

parent 3a6d5d10
......@@ -35,7 +35,8 @@
* @param c The #cell we are working on.
* @param timer Are we timing this ?
*/
static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
int timer) {
/* Some constants */
const struct engine *e = r->e;
......@@ -126,8 +127,9 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, int tim
* @param ci The #cell with field tensor to interact.
* @param cj The #cell with the multipole.
*/
static INLINE void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
static INLINE void runner_dopair_grav_mm(const struct runner *r,
struct cell *restrict ci,
struct cell *restrict cj) {
/* Some constants */
const struct engine *e = r->e;
......@@ -440,7 +442,8 @@ static INLINE void runner_dopair_grav_pm(
* @param ci The first #cell.
* @param cj The other #cell.
*/
static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
struct cell *cj) {
const struct engine *e = r->e;
......@@ -641,7 +644,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, stru
*
* @todo Use a local cache for the particles.
*/
static INLINE void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
static INLINE void runner_doself_grav_pp_full(struct runner *r,
struct cell *c) {
/* Some constants */
const struct engine *const e = r->e;
......@@ -763,7 +767,8 @@ static INLINE void runner_doself_grav_pp_full(struct runner *r, struct cell *c)
*
* @todo Use a local cache for the particles.
*/
static INLINE void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
struct cell *c) {
/* Some constants */
const struct engine *const e = r->e;
......@@ -945,8 +950,8 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
*
* @todo Use a local cache for the particles.
*/
static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
int gettimer) {
static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci,
struct cell *cj, int gettimer) {
/* Some constants */
const struct engine *e = r->e;
......@@ -1097,7 +1102,8 @@ static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci, struct
*
* @todo Use a local cache for the particles.
*/
static INLINE void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
static INLINE void runner_doself_grav(struct runner *r, struct cell *c,
int gettimer) {
/* Some constants */
const struct engine *e = r->e;
......@@ -1148,7 +1154,8 @@ static INLINE void runner_doself_grav(struct runner *r, struct cell *c, int gett
* @param ci The #cell of interest.
* @param timer Are we timing this ?
*/
static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
int timer) {
/* Some constants */
const struct engine *e = r->e;
......
......@@ -19,39 +19,69 @@
#include "../config.h"
/* Some standard headers. */
#include <math.h>
#include <fenv.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
/* Local headers. */
#include "swift.h"
#include "runner_doiact_grav.h"
#include "swift.h"
const int num_tests = 100;
const double eps = 0.01;
const double eps = 0.02;
/**
* @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
*/
void check_value(double a, double b, const char *s) {
if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
}
/* Definitions of the potential and force that match
exactly the theory document */
double S(double x) {
return exp(x) / (1. + exp(x));
}
double S(double x) { return exp(x) / (1. + exp(x)); }
double potential(double r, double H, double rlr) {
double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }
double potential(double mass, double r, double H, double rlr) {
const double u = r / H;
const double x = r / rlr;
double pot;
if (u > 1.)
pot = -1. / r;
pot = -mass / r;
else
pot = -(-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;
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;
return pot * (2. - 2.*S(2.*x));
return pot * (2. - 2. * S(2. * x));
}
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
acc = -mass * (21. * u * u * u * u * u - 90. * u * u * u * u +
140. * u * u * u - 84. * u * u + 14.) /
(H * H * H);
return r * acc * (4. * x * S_prime(2 * x) - 2. * S(2. * x) + 2.);
}
int main() {
/* Initialize CPU frequency, this also starts time. */
......@@ -70,32 +100,37 @@ int main() {
s.width[1] = 0.2;
s.width[2] = 0.2;
e.s = &s;
struct gravity_props props;
props.a_smooth = 1.25;
e.gravity_properties = &props;
struct runner r;
bzero(&r, sizeof(struct runner));
r.e = &e;
const double rlr = props.a_smooth * s.width[0];
/* Init the cache for gravity interaction */
gravity_cache_init(&r.ci_gravity_cache, num_tests * 2);
/* Let's create one cell with a massive particle and a bunch of test particles */
/* Let's create one cell with a massive particle and a bunch of test particles
*/
struct cell c;
bzero(&c, sizeof(struct cell));
c.width[0] = 1.;
c.width[1] = 1.;
c.width[2] = 1.;
c.loc[0] = 0.;
c.loc[1] = 0.;
c.loc[2] = 0.;
c.gcount = 1 + num_tests;
c.ti_old_gpart = 8;
c.ti_gravity_end_min = 8;
c.ti_gravity_end_max = 8;
posix_memalign((void**) &c.gparts, gpart_align, c.gcount * sizeof(struct gpart));
posix_memalign((void **)&c.gparts, gpart_align,
c.gcount * sizeof(struct gpart));
bzero(c.gparts, c.gcount * sizeof(struct gpart));
/* Create the massive particle */
......@@ -110,13 +145,13 @@ int main() {
#ifdef SWIFT_DEBUG_CHECKS
c.gparts[0].ti_drift = 8;
#endif
/* Create the mass-less particles */
for (int n = 1; n < num_tests + 1; ++n) {
struct gpart *gp = &c.gparts[n];
gp->x[0] = n / ((double) num_tests);
gp->x[0] = n / ((double)num_tests);
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
......@@ -128,16 +163,24 @@ int main() {
gp->ti_drift = 8;
#endif
}
/* Now compute the forces */
runner_doself_grav_pp_truncated(&r, &c);
/* Verify everything */
for (int n = 1; n < num_tests + 1; ++n) {
const struct gpart *gp = &c.gparts[n];
message("x=%f pot=%f true=%f",
gp->x[0], gp->potential, potential(gp->x[0], gp->epsilon, rlr));
double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
double acc_true =
acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
// message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
check_value(gp->potential, pot_true, "potential");
check_value(gp->a_grav[0], acc_true, "acceleration");
}
free(c.gparts);
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment