Skip to content
GitLab
Commit 2a597d23 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

0th order calculation now correct.

parent 0b3df5fe
Hide whitespace changes
Inline Side-by-side
src/engine.c
View file @ 2a597d23
......@@ -2980,7 +2980,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to all-reduce total mass in the system.");
#endif
message("Total mass in the system: %e", e->s->total_mass);
if (e->nodeID == 0) message("Total mass in the system: %e", e->s->total_mass);
#endif
/* Now time to get ready for the first time-step */
......@@ -3079,7 +3079,7 @@ void engine_step(struct engine *e) {
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Check the accuracy of the gravity calculation */
gravity_exact_force_check(e->s, e, 1e-1);
gravity_exact_force_check(e->s, e, 1e1);
#endif
/* Collect the values of rebuild from all nodes. */
......
src/gravity.c
View file @ 2a597d23
......@@ -149,9 +149,12 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
/* Check that we are not above tolerance */
if (fabsf(err_rel[0]) > rel_tol || fabsf(err_rel[1]) > rel_tol ||
fabsf(err_rel[2]) > rel_tol)
error("Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e]",
gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
error(
"Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e], "
"gp->mass_interacted=%e",
gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
gpi->mass_interacted);
/* Construct some statistics */
for (int k = 0; k < 3; ++k) {
......
src/gravity/Default/gravity.h
View file @ 2a597d23
......@@ -113,7 +113,7 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
gp->time_bin = 0;
gp->epsilon = 0.; // MATTHIEU
gp->epsilon = 0.1; // MATTHIEU
gravity_init_gpart(gp);
}
......
src/kernel_long_gravity.h
View file @ 2a597d23
......@@ -37,14 +37,15 @@
__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
float u, float *const W) {
const float arg1 = u * 0.5f;
const float arg2 = u * one_over_sqrt_pi;
const float arg3 = -arg1 * arg1;
/* const float arg1 = u * 0.5f; */
/* const float arg2 = u * one_over_sqrt_pi; */
/* const float arg3 = -arg1 * arg1; */
const float term1 = erfcf(arg1);
const float term2 = arg2 * expf(arg3);
/* const float term1 = erfcf(arg1); */
/* const float term2 = arg2 * expf(arg3); */
*W = term1 + term2;
/* *W = term1 + term2; */
*W = 1.f;
}
#endif // SWIFT_KERNEL_LONG_GRAVITY_H
src/multipole.h
View file @ 2a597d23
......@@ -104,7 +104,7 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
bzero(m, sizeof(struct gravity_tensors));
}
INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
INLINE static void gravity_field_tensors_init(struct gravity_tensors *m) {
bzero(&m->a_x, sizeof(struct acc_tensor));
bzero(&m->a_y, sizeof(struct acc_tensor));
......@@ -116,6 +116,19 @@ INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
#endif
}
/**
* @brief Adds field tensrs to other ones (i.e. does la += lb).
*
* @param la The gravity tensors to add to.
* @param lb The gravity tensors to add.
*/
INLINE static void gravity_field_tensors_add(struct gravity_tensors *la,
const struct gravity_tensors *lb) {
la->a_x.F_000 += lb->a_x.F_000;
la->a_y.F_000 += lb->a_y.F_000;
la->a_z.F_000 += lb->a_z.F_000;
}
/**
* @brief Prints the content of a #multipole to stdout.
*
......@@ -322,9 +335,9 @@ INLINE static void gravity_M2L(struct gravity_tensors *l_a,
const double r_inv = 1. / sqrt(r2);
/* 1st order multipole term */
l_a->a_x.F_000 = D_100(dx, dy, dz, r_inv) * m_b->mass;
l_a->a_y.F_000 = D_010(dx, dy, dz, r_inv) * m_b->mass;
l_a->a_z.F_000 = D_001(dx, dy, dz, r_inv) * m_b->mass;
l_a->a_x.F_000 += D_100(dx, dy, dz, r_inv) * m_b->mass;
l_a->a_y.F_000 += D_010(dx, dy, dz, r_inv) * m_b->mass;
l_a->a_z.F_000 += D_001(dx, dy, dz, r_inv) * m_b->mass;
#ifdef SWIFT_DEBUG_CHECKS
l_a->mass_interacted += m_b->mass;
......@@ -345,123 +358,28 @@ INLINE static void gravity_M2L(struct gravity_tensors *l_a,
INLINE static void gravity_L2L(struct gravity_tensors *l_a,
const struct gravity_tensors *l_b,
const double pos_a[3], const double pos_b[3],
int periodic) {}
int periodic) {
error("Not implemented yet");
}
/**
* @brief Applies the #acc_tensor to a set of #gpart.
* @brief Applies the #acc_tensor to a #gpart.
*
* Corresponds to equation (28a).
*/
INLINE static void gravity_L2P(const struct gravity_tensors *l,
struct gpart *gparts, int gcount) {
for (int i = 0; i < gcount; ++i) {
struct gpart *gp) {
#ifdef SWIFT_DEBUG_CHECKS
struct gpart *gp = &gparts[i];
// if(gpart_is_active(gp, e)){
gp->mass_interacted += l->mass_interacted;
gp->mass_interacted += l->mass_interacted;
#endif
//}
}
}
#if 0
/* Multipole function prototypes. */
void multipole_add(struct gravity_tensors *m_sum, const struct gravity_tensors *m_term);
void multipole_init(struct gravity_tensors *m, const struct gpart *gparts,
int gcount);
void multipole_reset(struct gravity_tensors *m);
/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
/* double *shift); */
/* void multipole_addpart(struct multipole *m, struct gpart *p); */
/* void multipole_addparts(struct multipole *m, struct gpart *p, int N); */
// message("a=[%e %e %e]", l->a_x.F_000, l->a_y.F_000, l->a_z.F_000);
/**
* @brief Compute the pairwise interaction between two multipoles.
*
* @param ma The first #multipole.
* @param mb The second #multipole.
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
struct gravity_tensors *ma, struct gravity_tensors *mb, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
/* /\* Compute the multipole distance. *\/ */
/* for (k = 0; k < 3; k++) { */
/* dx[k] = ma->x[k] - mb->x[k] - shift[k]; */
/* r2 += dx[k] * dx[k]; */
/* } */
/* /\* Compute the normalized distance vector. *\/ */
/* ir = 1.0f / sqrtf(r2); */
/* r = r2 * ir; */
/* /\* Evaluate the gravity kernel. *\/ */
/* kernel_grav_eval(r, &acc); */
/* /\* Scale the acceleration. *\/ */
/* acc *= const_G * ir * ir * ir; */
/* /\* Compute the forces on both multipoles. *\/ */
/* #if const_gravity_multipole_order == 1 */
/* float mma = ma->coeffs[0], mmb = mb->coeffs[0]; */
/* for (k = 0; k < 3; k++) { */
/* ma->a[k] -= dx[k] * acc * mmb; */
/* mb->a[k] += dx[k] * acc * mma; */
/* } */
/* #else */
/* #error( "Multipoles of order %i not yet implemented." ,
* const_gravity_multipole_order )
*/
/* #endif */
/* 0th order interaction */
gp->a_grav[0] += l->a_x.F_000;
gp->a_grav[1] += l->a_y.F_000;
gp->a_grav[2] += l->a_z.F_000;
}
/**
* @brief Compute the interaction of a multipole on a particle.
*
* @param m The #multipole.
* @param p The #gpart.
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
struct gravity_tensors *m, struct gpart *p, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
/* /\* Compute the multipole distance. *\/ */
/* for (k = 0; k < 3; k++) { */
/* dx[k] = m->x[k] - p->x[k] - shift[k]; */
/* r2 += dx[k] * dx[k]; */
/* } */
/* /\* Compute the normalized distance vector. *\/ */
/* ir = 1.0f / sqrtf(r2); */
/* r = r2 * ir; */
/* /\* Evaluate the gravity kernel. *\/ */
/* kernel_grav_eval(r, &acc); */
/* /\* Scale the acceleration. *\/ */
/* acc *= const_G * ir * ir * ir * m->coeffs[0]; */
/* /\* Compute the forces on both multipoles. *\/ */
/* #if const_gravity_multipole_order == 1 */
/* for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc; */
/* #else */
/* #error( "Multipoles of order %i not yet implemented." ,
* const_gravity_multipole_order )
*/
/* #endif */
}
#endif
#endif /* SWIFT_MULTIPOLE_H */
src/runner.c
View file @ 2a597d23
......@@ -497,7 +497,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Reset the gravity acceleration tensors */
if (e->policy & engine_policy_self_gravity)
gravity_field_tensor_init(c->multipole);
gravity_field_tensors_init(c->multipole);
/* Recurse? */
if (c->split) {
......
src/runner_doiact_grav.h
View file @ 2a597d23
......@@ -64,12 +64,22 @@ void runner_do_grav_down(struct runner *r, struct cell *c) {
if (cp != NULL) {
gravity_L2L(&temp, c->multipole, cp->multipole->CoM, c->multipole->CoM,
1);
gravity_field_tensors_add(cp->multipole, &temp);
}
}
} else {
gravity_L2P(c->multipole, c->gparts, c->gcount);
const struct engine *e = r->e;
struct gpart *gparts = c->gparts;
const int gcount = c->gcount;
/* Apply accelerations to the particles */
for (int i = 0; i < gcount; ++i) {
struct gpart *gp = &gparts[i];
if (gpart_is_active(gp, e)) gravity_L2P(c->multipole, gp);
}
}
}
......@@ -237,49 +247,55 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
/* Get a hold of the ith part in ci. */
struct gpart *restrict gpi = &gparts_i[pid];
if (!gpart_is_active(gpi, e)) continue;
/* Loop over every particle in the other cell. */
for (int pjd = 0; pjd < gcount_j; pjd++) {
/* Get a hold of the jth part in cj. */
struct gpart *restrict gpj = &gparts_j[pjd];
const struct gpart *restrict gpj = &gparts_j[pjd];
/* Compute the pairwise distance. */
float dx[3] = {gpi->x[0] - gpj->x[0], // x
gpi->x[1] - gpj->x[1], // y
gpi->x[2] - gpj->x[2]}; // z
const float dx[3] = {gpi->x[0] - gpj->x[0], // x
gpi->x[1] - gpj->x[1], // y
gpi->x[2] - gpj->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact ! */
if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
gpi->mass_interacted += gpj->mass;
gpj->mass_interacted += gpi->mass;
gpi->mass_interacted += gpj->mass;
#endif
}
}
} else {
/* Loop over all particles in cj... */
for (int pjd = 0; pjd < gcount_j; pjd++) {
if (gpart_is_active(gpi, e)) {
/* Get a hold of the ith part in ci. */
struct gpart *restrict gpj = &gparts_j[pjd];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
if (!gpart_is_active(gpj, e)) continue;
#ifdef SWIFT_DEBUG_CHECKS
gpi->mass_interacted += gpj->mass;
#endif
} else if (gpart_is_active(gpj, e)) {
/* Loop over every particle in the other cell. */
for (int pid = 0; pid < gcount_i; pid++) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
/* Get a hold of the ith part in ci. */
const struct gpart *restrict gpi = &gparts_i[pid];
/* Compute the pairwise distance. */
const float dx[3] = {gpj->x[0] - gpi->x[0], // x
gpj->x[1] - gpi->x[1], // y
gpj->x[2] - gpi->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact ! */
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
gpj->mass_interacted += gpi->mass;
gpj->mass_interacted += gpi->mass;
#endif
}
}
}
}
......@@ -468,9 +484,9 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
} else {
/* Ok, here we can go for particle-multipole interactions */
runner_dopair_grav_pm(r, ci->progeny[j], cj->progeny[k]);
runner_dopair_grav_pm(r, cj->progeny[k], ci->progeny[j]);
/* Ok, here we can go for multipole-multipole interactions */
runner_dopair_grav_mm(r, ci->progeny[j], cj->progeny[k]);
runner_dopair_grav_mm(r, cj->progeny[k], ci->progeny[j]);
}
}
}
......
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