Commit 4ad965ea authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Re-instated an estimate of the evolution of multipole radii for the multipole drift.

parent e5738a31
......@@ -2707,7 +2707,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
/* Drift the multipole */
if (ti_current > ti_old_multipole)
gravity_drift(c->multipole, dt_drift, 0.f /*c->dx_max_gpart*/);
gravity_drift(c->multipole, dt_drift);
/* Are we not in a leaf ? */
if (c->split) {
......@@ -2749,7 +2749,7 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
dt_drift = (ti_current - ti_old_multipole) * e->time_base;
if (ti_current > ti_old_multipole)
gravity_drift(c->multipole, dt_drift, 0.f /*c->dx_max_gpart*/);
gravity_drift(c->multipole, dt_drift);
/* Update the time of the last drift */
c->ti_old_multipole = ti_current;
......
......@@ -108,10 +108,16 @@ struct grav_tensor {
struct multipole {
/* Bulk velocity */
/*! Bulk velocity */
float vel[3];
/* 0th order terms */
/*! Maximal velocity along each axis of all #gpart */
float max_delta_vel[3];
/*! Minimal velocity along each axis of all #gpart */
float min_delta_vel[3];
/* 0th order term */
float M_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
......@@ -215,14 +221,15 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
/**
* @brief Drifts a #multipole forward in time.
*
* This uses a first-order approximation in time. We only move the CoM
* using the bulk velocity measured at the last rebuild.
*
* @param m The #multipole.
* @param dt The drift time-step.
* @param x_diff The maximal distance moved by any particle since the last
* rebuild.
*/
INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
float x_diff) {
INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
/* Motion of the centre of mass */
const double dx = m->m_pole.vel[0] * dt;
const double dy = m->m_pole.vel[1] * dt;
const double dz = m->m_pole.vel[2] * dt;
......@@ -232,8 +239,35 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
m->CoM[1] += dy;
m->CoM[2] += dz;
#ifdef SWIFT_DEBUG_CHECKS
if(m->m_pole.vel[0] > m->m_pole.max_delta_vel[0])
error("Invalid maximal velocity");
if(m->m_pole.vel[0] < m->m_pole.min_delta_vel[0])
error("Invalid minimal velocity");
if(m->m_pole.vel[1] > m->m_pole.max_delta_vel[1])
error("Invalid maximal velocity");
if(m->m_pole.vel[1] < m->m_pole.min_delta_vel[1])
error("Invalid minimal velocity");
if(m->m_pole.vel[2] > m->m_pole.max_delta_vel[2])
error("Invalid maximal velocity");
if(m->m_pole.vel[2] < m->m_pole.min_delta_vel[2])
error("Invalid minimal velocity");
#endif
/* Maximal distance covered by any particle */
float dv[3];
dv[0] = max(m->m_pole.max_delta_vel[0] - m->m_pole.vel[0],
m->m_pole.vel[0] - m->m_pole.min_delta_vel[0]);
dv[1] = max(m->m_pole.max_delta_vel[1] - m->m_pole.vel[1],
m->m_pole.vel[1] - m->m_pole.min_delta_vel[1]);
dv[2] = max(m->m_pole.max_delta_vel[2] - m->m_pole.vel[2],
m->m_pole.vel[2] - m->m_pole.min_delta_vel[2]);
const float max_delta_vel = sqrt(dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]);
const float x_diff = max_delta_vel * dt;
/* Conservative change in maximal radius containing all gpart */
m->r_max = m->r_max_rebuild + x_diff;
m->r_max += x_diff;
}
/**
......@@ -470,16 +504,8 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
INLINE static void gravity_multipole_add(struct multipole *ma,
const struct multipole *mb) {
const float M_000 = ma->M_000 + mb->M_000;
const float inv_M_000 = 1.f / M_000;
/* Add the bulk velocities */
ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;
/* Add 0th order terms */
ma->M_000 = M_000;
/* Add 0th order term */
ma->M_000 += mb->M_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* Add 1st order terms */
......@@ -555,11 +581,6 @@ INLINE static void gravity_multipole_add(struct multipole *ma,
#error "Missing implementation for order >5"
#endif
// MATTHIEU
ma->M_100 = 0.f;
ma->M_010 = 0.f;
ma->M_001 = 0.f;
#ifdef SWIFT_DEBUG_CHECKS
ma->num_gpart += mb->num_gpart;
#endif
......@@ -1026,6 +1047,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
/* Prepare some local counters */
double r_max2 = 0.;
float max_delta_vel[3] = {0., 0., 0.};
float min_delta_vel[3] = {0., 0., 0.};
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
double M_100 = 0., M_010 = 0., M_001 = 0.;
#endif
......@@ -1068,6 +1092,16 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
/* Maximal distance CoM<->gpart */
r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
/* Store the vector of the maximal vel difference */
max_delta_vel[0] = max(gparts[k].v_full[0], max_delta_vel[0]);
max_delta_vel[1] = max(gparts[k].v_full[1], max_delta_vel[1]);
max_delta_vel[2] = max(gparts[k].v_full[2], max_delta_vel[2]);
/* Store the vector of the minimal vel difference */
min_delta_vel[0] = min(gparts[k].v_full[0], min_delta_vel[0]);
min_delta_vel[1] = min(gparts[k].v_full[1], min_delta_vel[1]);
min_delta_vel[2] = min(gparts[k].v_full[2], min_delta_vel[2]);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
const double m = gparts[k].mass;
......@@ -1150,7 +1184,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
}
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
M_100 = M_010 = M_001 = 0.f; /* Matthieu */
/* We know the first-order multipole (dipole) is 0. */
M_100 = M_010 = M_001 = 0.f;
#endif
/* Store the data on the multipole. */
......@@ -1162,6 +1198,12 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
multi->m_pole.vel[0] = vel[0];
multi->m_pole.vel[1] = vel[1];
multi->m_pole.vel[2] = vel[2];
multi->m_pole.max_delta_vel[0] = max_delta_vel[0];
multi->m_pole.max_delta_vel[1] = max_delta_vel[1];
multi->m_pole.max_delta_vel[2] = max_delta_vel[2];
multi->m_pole.min_delta_vel[0] = min_delta_vel[0];
multi->m_pole.min_delta_vel[1] = min_delta_vel[1];
multi->m_pole.min_delta_vel[2] = min_delta_vel[2];
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
......@@ -1260,10 +1302,6 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
INLINE static void gravity_M2M(struct multipole *m_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3]) {
/* Shift bulk velocity */
m_a->vel[0] = m_b->vel[0];
m_a->vel[1] = m_b->vel[1];
m_a->vel[2] = m_b->vel[2];
/* Shift 0th order term */
m_a->M_000 = m_b->M_000;
......
......@@ -1886,22 +1886,53 @@ void space_split_recursive(struct space *s, struct cell *c,
/* Reset everything */
gravity_reset(c->multipole);
/* Compute CoM of all progenies */
/* Compute CoM and bulk velocity from all progenies */
double CoM[3] = {0., 0., 0.};
double vel[3] = {0., 0., 0.};
float max_delta_vel[3] = {0.f, 0.f, 0.f};
float min_delta_vel[3] = {0.f, 0.f, 0.f};
double mass = 0.;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct gravity_tensors *m = c->progeny[k]->multipole;
mass += m->m_pole.M_000;
CoM[0] += m->CoM[0] * m->m_pole.M_000;
CoM[1] += m->CoM[1] * m->m_pole.M_000;
CoM[2] += m->CoM[2] * m->m_pole.M_000;
mass += m->m_pole.M_000;
vel[0] += m->m_pole.vel[0] * m->m_pole.M_000;
vel[1] += m->m_pole.vel[1] * m->m_pole.M_000;
vel[2] += m->m_pole.vel[2] * m->m_pole.M_000;
max_delta_vel[0] = max(m->m_pole.max_delta_vel[0], max_delta_vel[0]);
max_delta_vel[1] = max(m->m_pole.max_delta_vel[1], max_delta_vel[1]);
max_delta_vel[2] = max(m->m_pole.max_delta_vel[2], max_delta_vel[2]);
min_delta_vel[0] = min(m->m_pole.min_delta_vel[0], min_delta_vel[0]);
min_delta_vel[1] = min(m->m_pole.min_delta_vel[1], min_delta_vel[1]);
min_delta_vel[2] = min(m->m_pole.min_delta_vel[2], min_delta_vel[2]);
}
}
c->multipole->CoM[0] = CoM[0] / mass;
c->multipole->CoM[1] = CoM[1] / mass;
c->multipole->CoM[2] = CoM[2] / mass;
/* Final operation on the CoM and bulk velocity */
const double inv_mass = 1. / mass;
c->multipole->CoM[0] = CoM[0] * inv_mass;
c->multipole->CoM[1] = CoM[1] * inv_mass;
c->multipole->CoM[2] = CoM[2] * inv_mass;
c->multipole->m_pole.vel[0] = vel[0] * inv_mass;
c->multipole->m_pole.vel[1] = vel[1] * inv_mass;
c->multipole->m_pole.vel[2] = vel[2] * inv_mass;
/* Min max velocity along each axis */
c->multipole->m_pole.max_delta_vel[0] = max_delta_vel[0];
c->multipole->m_pole.max_delta_vel[1] = max_delta_vel[1];
c->multipole->m_pole.max_delta_vel[2] = max_delta_vel[2];
c->multipole->m_pole.min_delta_vel[0] = min_delta_vel[0];
c->multipole->m_pole.min_delta_vel[1] = min_delta_vel[1];
c->multipole->m_pole.min_delta_vel[2] = min_delta_vel[2];
/* Now shift progeny multipoles and add them up */
struct multipole temp;
......@@ -1923,6 +1954,7 @@ void space_split_recursive(struct space *s, struct cell *c,
r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
}
}
/* Alternative upper limit of max CoM<->gpart distance */
const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
? c->multipole->CoM[0] - c->loc[0]
......@@ -1943,6 +1975,11 @@ void space_split_recursive(struct space *s, struct cell *c,
c->multipole->CoM_rebuild[1] = c->multipole->CoM[1];
c->multipole->CoM_rebuild[2] = c->multipole->CoM[2];
/* We know the first-order multipole (dipole) is 0. */
c->multipole->m_pole.M_100 = 0.f;
c->multipole->m_pole.M_010 = 0.f;
c->multipole->m_pole.M_001 = 0.f;
} /* Deal with gravity */
} /* Split or let it be? */
......
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