diff --git a/src/cell.c b/src/cell.c index ab5b99c03024b1cc218c358b6d7051632496c665..753bdd55061ebd2b2cdd2067691bd8319ca5a623 100644 --- a/src/cell.c +++ b/src/cell.c @@ -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 diff --git a/src/multipole.h b/src/multipole.h index 6065a8563359a888f7ebd11db55e9b124ae97a43..f166a6078b16ec0ef297217ab1c19e8071bf8d16 100644 --- a/src/multipole.h +++ b/src/multipole.h @@ -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 diff --git a/src/runner.c b/src/runner.c index dc831855930f01dd745273e553e8f593c6238b92..521914bfdecb18f6c5e4e0a55ec2f0e759424797 100644 --- a/src/runner.c +++ b/src/runner.c @@ -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); diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index c9e95dc47e8cd7c31b8e8bf76f1733d7f0d6bdcc..9988b7d553a962ee541f20cab64cc534c991957a 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -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. diff --git a/src/space.c b/src/space.c index 6199522d94131c9f2ebc25306aebce4556aa00fa..625fe944c488c871d02b14c72d15051b3e53b3c4 100644 --- a/src/space.c +++ b/src/space.c @@ -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. */