Skip to content
GitLab
Commit 52c52f4d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added the gravity downward pass.

parent 391c962b
Hide whitespace changes
Inline Side-by-side
src/cell.c
View file @ 52c52f4d
......@@ -1124,15 +1124,16 @@ void cell_check_multipole(struct cell *c, void *data) {
if (c->gcount > 0) {
/* Brute-force calculation */
multipole_P2M(&ma, c->gparts, c->gcount);
gravity_P2M(&ma, c->gparts, c->gcount);
/* Now compare the multipole expansion */
if (!multipole_equal(&ma.m_pole, &c->multipole->m_pole, tolerance)) {
if (!gravity_multipole_equal(&ma.m_pole, &c->multipole->m_pole,
tolerance)) {
message("Multipoles are not equal at depth=%d!", c->depth);
message("Correct answer:");
multipole_print(&ma.m_pole);
gravity_multipole_print(&ma.m_pole);
message("Recursive multipole:");
multipole_print(&c->multipole->m_pole);
gravity_multipole_print(&c->multipole->m_pole);
error("Aborting");
}
}
......@@ -1486,7 +1487,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
} else if (ti_current > ti_old_multipole) {
/* Drift the multipole */
multipole_drift(c->multipole, dt);
gravity_multipole_drift(c->multipole, dt);
}
/* Update the time of the last drift */
......@@ -1502,7 +1503,9 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_multipole(struct cell *c, const struct engine *e) {
error("To be implemented");
}
/**
* @brief Recursively checks that all particles in a cell have a time-step
......
src/multipole.h
View file @ 52c52f4d
......@@ -25,6 +25,7 @@
#include <string.h>
/* Includes. */
//#include "active.h"
#include "align.h"
#include "const.h"
#include "error.h"
......@@ -82,6 +83,13 @@ struct gravity_tensors {
/*! Field tensor for the potential */
struct pot_tensor pot;
#ifdef SWIFT_DEBUG_CHECKS
/* Total mass this gpart interacted with */
double mass_interacted;
#endif
};
} SWIFT_STRUCT_ALIGN;
......@@ -90,18 +98,22 @@ struct gravity_tensors {
*
* @param m The #multipole.
*/
INLINE static void multipole_reset(struct gravity_tensors *m) {
INLINE static void gravity_reset(struct gravity_tensors *m) {
/* Just bzero the struct. */
bzero(m, sizeof(struct gravity_tensors));
}
INLINE static void multipole_init(struct gravity_tensors *m) {
INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
bzero(&m->a_x, sizeof(struct acc_tensor));
bzero(&m->a_y, sizeof(struct acc_tensor));
bzero(&m->a_z, sizeof(struct acc_tensor));
bzero(&m->pot, sizeof(struct pot_tensor));
#ifdef SWIFT_DEBUG_CHECKS
m->mass_interacted = 0.;
#endif
}
/**
......@@ -111,7 +123,7 @@ INLINE static void multipole_init(struct gravity_tensors *m) {
*
* @param m The #multipole to print.
*/
INLINE static void multipole_print(const struct multipole *m) {
INLINE static void gravity_multipole_print(const struct multipole *m) {
// printf("CoM= [%12.5e %12.5e %12.5e\n", m->CoM[0], m->CoM[1], m->CoM[2]);
printf("Mass= %12.5e\n", m->mass);
......@@ -124,8 +136,8 @@ INLINE static void multipole_print(const struct multipole *m) {
* @param ma The multipole to add to.
* @param mb The multipole to add.
*/
INLINE static void multipole_add(struct multipole *ma,
const struct multipole *mb) {
INLINE static void gravity_multipole_add(struct multipole *ma,
const struct multipole *mb) {
const float mass = ma->mass + mb->mass;
const float imass = 1.f / mass;
......@@ -147,9 +159,9 @@ INLINE static void multipole_add(struct multipole *ma,
* @param tolerance The maximal allowed difference for the fields.
* @return 1 if the multipoles are equal 0 otherwise.
*/
INLINE static int multipole_equal(const struct multipole *ma,
const struct multipole *mb,
double tolerance) {
INLINE static int gravity_multipole_equal(const struct multipole *ma,
const struct multipole *mb,
double tolerance) {
/* Check CoM */
/* if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) >
......@@ -190,7 +202,8 @@ INLINE static int multipole_equal(const struct multipole *ma,
* @param m The #multipole.
* @param dt The drift time-step.
*/
INLINE static void multipole_drift(struct gravity_tensors *m, double dt) {
INLINE static void gravity_multipole_drift(struct gravity_tensors *m,
double dt) {
/* Move the whole thing according to bulk motion */
m->CoM[0] += m->m_pole.vel[0];
......@@ -206,8 +219,8 @@ INLINE static void multipole_drift(struct gravity_tensors *m, double dt) {
* @param gparts_j The #gpart that source the gravity field.
* @param gcount_j The number of sources.
*/
INLINE static void multipole_P2P(struct gpart *gparts_i, int gcount_i,
const struct gpart *gparts_j, int gcount_j) {}
INLINE static void gravity_P2P(struct gpart *gparts_i, int gcount_i,
const struct gpart *gparts_j, int gcount_j) {}
/**
* @brief Constructs the #multipole of a bunch of particles around their
......@@ -219,8 +232,8 @@ INLINE static void multipole_P2P(struct gpart *gparts_i, int gcount_i,
* @param gparts The #gpart.
* @param gcount The number of particles.
*/
INLINE static void multipole_P2M(struct gravity_tensors *m,
const struct gpart *gparts, int gcount) {
INLINE static void gravity_P2M(struct gravity_tensors *m,
const struct gpart *gparts, int gcount) {
#if const_gravity_multipole_order >= 2
#error "Implementation of P2M kernel missing for this order."
......@@ -267,10 +280,10 @@ INLINE static void multipole_P2M(struct gravity_tensors *m,
* @param pos_b The current postion of the multipole to shift.
* @param periodic Is the calculation periodic ?
*/
INLINE static void multipole_M2M(struct multipole *m_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3],
int periodic) {
INLINE static void gravity_M2M(struct multipole *m_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3],
int periodic) {
m_a->mass = m_b->mass;
......@@ -289,10 +302,10 @@ INLINE static void multipole_M2M(struct multipole *m_a,
* @param pos_b The position of the multipole.
* @param periodic Is the calculation periodic ?
*/
INLINE static void multipole_M2L(struct gravity_tensors *l_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3],
int periodic) {
INLINE static void gravity_M2L(struct gravity_tensors *l_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3],
int periodic) {
double dx, dy, dz;
if (periodic) {
......@@ -309,9 +322,48 @@ INLINE static void multipole_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;
#endif
}
/**
* @brief Creates a copy of #acc_tensor shifted to a new location.
*
* Corresponds to equation (28e).
*
* @param l_a The #acc_tensor copy (content will be overwritten).
* @param l_b The #acc_tensor to shift.
* @param pos_a The position to which m_b will be shifted.
* @param pos_b The current postion of the multipole to shift.
* @param periodic Is the calculation periodic ?
*/
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) {}
/**
* @brief Applies the #acc_tensor to a set of #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 = &gparts[i];
// if(gpart_is_active(gp, e)){
gp->mass_interacted += l->mass_interacted;
//}
}
}
#if 0
......
src/runner.c
View file @ 52c52f4d
......@@ -496,7 +496,8 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (!cell_is_active(c, e)) return;
/* Reset the gravity acceleration tensors */
if (e->policy & engine_policy_self_gravity) multipole_init(c->multipole);
if (e->policy & engine_policy_self_gravity)
gravity_field_tensor_init(c->multipole);
/* Recurse? */
if (c->split) {
......@@ -1808,7 +1809,7 @@ void *runner_main(void *data) {
// runner_do_grav_mm(r, t->ci, 1);
break;
case task_type_grav_down:
// runner_do_grav_down(r, t->ci);
runner_do_grav_down(r, t->ci);
break;
case task_type_grav_top_level:
// runner_do_grav_top_level(r);
......
src/runner_doiact_grav.h
View file @ 52c52f4d
......@@ -53,6 +53,26 @@ void runner_do_grav_up(struct runner *r, struct cell *c) {
/* } */
}
void runner_do_grav_down(struct runner *r, struct cell *c) {
if (c->split) {
for (int k = 0; k < 8; ++k) {
struct cell *cp = c->progeny[k];
struct gravity_tensors temp;
if (cp != NULL) {
gravity_L2L(&temp, c->multipole, cp->multipole->CoM, c->multipole->CoM,
1);
}
}
} else {
gravity_L2P(c->multipole, c->gparts, c->gcount);
}
}
/**
* @brief Computes the interaction of the field tensor in a cell with the
* multipole of another cell.
......@@ -68,27 +88,24 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
const struct engine *e = r->e;
const int periodic = e->s->periodic;
const struct multipole *multi_j = &cj->multipole->m_pole;
//const float a_smooth = e->gravity_properties->a_smooth;
//const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
// const float a_smooth = e->gravity_properties->a_smooth;
// const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
if (multi_j->mass == 0.0)
error("Multipole does not seem to have been set.");
if (multi_j->mass == 0.0) error("Multipole does not seem to have been set.");
#endif
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
multipole_M2L(ci->multipole, multi_j, ci->multipole->CoM, cj->multipole->CoM,
periodic);
gravity_M2L(ci->multipole, multi_j, ci->multipole->CoM, cj->multipole->CoM,
periodic);
TIMER_TOC(timer_dopair_grav_mm);
}
/**
* @brief Computes the interaction of all the particles in a cell with the
* multipole of another cell.
......
src/space.c
View file @ 52c52f4d
......@@ -2091,7 +2091,7 @@ void space_split_recursive(struct space *s, struct cell *c,
if (s->gravity) {
/* Reset everything */
multipole_reset(c->multipole);
gravity_reset(c->multipole);
/* Compute CoM of all progenies */
double CoM[3] = {0., 0., 0.};
......@@ -2116,9 +2116,9 @@ void space_split_recursive(struct space *s, struct cell *c,
if (c->progeny[k] != NULL) {
const struct cell *cp = c->progeny[k];
const struct multipole *m = &cp->multipole->m_pole;
multipole_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
s->periodic);
multipole_add(&c->multipole->m_pole, &temp);
gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
s->periodic);
gravity_multipole_add(&c->multipole->m_pole, &temp);
}
}
}
......@@ -2174,7 +2174,7 @@ void space_split_recursive(struct space *s, struct cell *c,
}
/* Construct the multipole and the centre of mass*/
if (s->gravity) multipole_P2M(c->multipole, c->gparts, c->gcount);
if (s->gravity) gravity_P2M(c->multipole, c->gparts, c->gcount);
}
/* Set the values for this cell. */
......
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