Skip to content
GitLab
Commit 68ce4a36 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Recursive multipole construction now fully implemented

parent 8d8aae94
Hide whitespace changes
Inline Side-by-side
src/Makefile.am
View file @ 68ce4a36
......@@ -45,7 +45,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......
src/cell.c
View file @ 68ce4a36
......@@ -1073,50 +1073,27 @@ int cell_are_neighbours(const struct cell *restrict ci,
void cell_check_multipole(struct cell *c, void *data) {
struct multipole ma;
const double tolerance = 1e-5; /* Relative */
/* First recurse */
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k], NULL);
if (c->gcount > 0) {
/* Brute-force calculation */
multipole_init(&ma, c->gparts, c->gcount);
/* Compare with recursive one */
struct multipole mb = c->multipole;
if (fabsf(ma.mass - mb.mass) / fabsf(ma.mass + mb.mass) > 1e-5)
error("Multipole masses are different (%12.15e vs. %12.15e)", ma.mass,
mb.mass);
for (int k = 0; k < 3; ++k)
if (fabs(ma.CoM[k] - mb.CoM[k]) / fabs(ma.CoM[k] + mb.CoM[k]) > 1e-5)
error("Multipole CoM are different (%12.15e vs. %12.15e", ma.CoM[k],
mb.CoM[k]);
#if const_gravity_multipole_order >= 2
if (fabsf(ma.I_xx - mb.I_xx) / fabsf(ma.I_xx + mb.I_xx) > 1e-5 &&
ma.I_xx > 1e-9)
error("Multipole I_xx are different (%12.15e vs. %12.15e)", ma.I_xx,
mb.I_xx);
if (fabsf(ma.I_yy - mb.I_yy) / fabsf(ma.I_yy + mb.I_yy) > 1e-5 &&
ma.I_yy > 1e-9)
error("Multipole I_yy are different (%12.15e vs. %12.15e)", ma.I_yy,
mb.I_yy);
if (fabsf(ma.I_zz - mb.I_zz) / fabsf(ma.I_zz + mb.I_zz) > 1e-5 &&
ma.I_zz > 1e-9)
error("Multipole I_zz are different (%12.15e vs. %12.15e)", ma.I_zz,
mb.I_zz);
if (fabsf(ma.I_xy - mb.I_xy) / fabsf(ma.I_xy + mb.I_xy) > 1e-5 &&
ma.I_xy > 1e-9)
error("Multipole I_xy are different (%12.15e vs. %12.15e)", ma.I_xy,
mb.I_xy);
if (fabsf(ma.I_xz - mb.I_xz) / fabsf(ma.I_xz + mb.I_xz) > 1e-5 &&
ma.I_xz > 1e-9)
error("Multipole I_xz are different (%12.15e vs. %12.15e)", ma.I_xz,
mb.I_xz);
if (fabsf(ma.I_yz - mb.I_yz) / fabsf(ma.I_yz + mb.I_yz) > 1e-5 &&
ma.I_yz > 1e-9)
error("Multipole I_yz are different (%12.15e vs. %12.15e)", ma.I_yz,
mb.I_yz);
#endif
multipole_P2M(&ma, c->gparts, c->gcount);
/* Now compare the multipole expansion */
if (!multipole_equal(&ma, c->multipole, tolerance)) {
message("Multipoles are not equal at depth=%d!", c->depth);
message("Correct answer:");
multipole_print(&ma);
message("Recursive multipole:");
multipole_print(c->multipole);
error("Aborting");
}
}
}
......
src/cell.h
View file @ 68ce4a36
......@@ -94,9 +94,6 @@ struct pcell {
*/
struct cell {
/*! This cell's multipole. */
struct multipole multipole;
/*! The cell location on the grid. */
double loc[3];
......@@ -106,6 +103,9 @@ struct cell {
/*! Max smoothing length in this cell. */
double h_max;
/*! This cell's multipole. */
struct multipole *multipole;
/*! Linking pointer for "memory management". */
struct cell *next;
......
src/engine.c
View file @ 68ce4a36
......@@ -3448,6 +3448,7 @@ void engine_unpin() {
* @param internal_units The system of units used internally.
* @param physical_constants The #phys_const used for this run.
* @param hydro The #hydro_props used for this run.
* @param gravity The #gravity_props used for this run.
* @param potential The properties of the external potential.
* @param cooling_func The properties of the cooling function.
* @param sourceterms The properties of the source terms function.
......
src/gravity_derivatives.h 0 → 100644
View file @ 68ce4a36
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_GRAVITY_DERIVATIVE_H
#define SWIFT_GRAVITY_DERIVATIVE_H
/* Some standard headers. */
#include <math.h>
/* Local headers. */
#include "inline.h"
/**
* @brief \f$ \phi(r_x, r_y, r_z) \f$.
*/
__attribute__((always_inline)) INLINE static double D_000(double r_x,
double r_y,
double r_z,
double r_inv) {
return r_inv;
}
/**
* @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
*/
__attribute__((always_inline)) INLINE static double D_100(double r_x,
double r_y,
double r_z,
double r_inv) {
return -r_x * r_inv * r_inv * r_inv;
}
/**
* @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
*/
__attribute__((always_inline)) INLINE static double D_010(double r_x,
double r_y,
double r_z,
double r_inv) {
return -r_y * r_inv * r_inv * r_inv;
}
/**
* @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
*/
__attribute__((always_inline)) INLINE static double D_001(double r_x,
double r_y,
double r_z,
double r_inv) {
return -r_z * r_inv * r_inv * r_inv;
}
#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
src/multipole.c
View file @ 68ce4a36
......@@ -20,205 +20,3 @@
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <strings.h>
/* This object's header. */
#include "multipole.h"
/**
* @brief Reset the data of a #multipole.
*
* @param m The #multipole.
*/
void multipole_reset(struct multipole *m) {
/* Just bzero the struct. */
bzero(m, sizeof(struct multipole));
}
/**
* @brief Init a multipole from a set of particles.
*
* @param m The #multipole.
* @param gparts The #gpart.
* @param gcount The number of particles.
*/
void multipole_init(struct multipole *m, const struct gpart *gparts,
int gcount) {
#if const_gravity_multipole_order > 2
#error "Multipoles of order >2 not yet implemented."
#endif
/* Zero everything */
multipole_reset(m);
/* Temporary variables */
double mass = 0.0;
double com[3] = {0.0, 0.0, 0.0};
#if const_gravity_multipole_order >= 2
double I_xx = 0.0, I_yy = 0.0, I_zz = 0.0;
double I_xy = 0.0, I_xz = 0.0, I_yz = 0.0;
#endif
/* Collect the particle data. */
for (int k = 0; k < gcount; k++) {
const float w = gparts[k].mass;
mass += w;
com[0] += gparts[k].x[0] * w;
com[1] += gparts[k].x[1] * w;
com[2] += gparts[k].x[2] * w;
#if const_gravity_multipole_order >= 2
I_xx += gparts[k].x[0] * gparts[k].x[0] * w;
I_yy += gparts[k].x[1] * gparts[k].x[1] * w;
I_zz += gparts[k].x[2] * gparts[k].x[2] * w;
I_xy += gparts[k].x[0] * gparts[k].x[1] * w;
I_xz += gparts[k].x[0] * gparts[k].x[2] * w;
I_yz += gparts[k].x[1] * gparts[k].x[2] * w;
#endif
}
const double imass = 1.0 / mass;
/* Store the data on the multipole. */
m->mass = mass;
m->CoM[0] = com[0] * imass;
m->CoM[1] = com[1] * imass;
m->CoM[2] = com[2] * imass;
#if const_gravity_multipole_order >= 2
m->I_xx = I_xx - imass * com[0] * com[0];
m->I_yy = I_yy - imass * com[1] * com[1];
m->I_zz = I_zz - imass * com[2] * com[2];
m->I_xy = I_xy - imass * com[0] * com[1];
m->I_xz = I_xz - imass * com[0] * com[2];
m->I_yz = I_yz - imass * com[1] * com[2];
#endif
}
/**
* @brief Add the second multipole to the first one.
*
* @param ma The #multipole which will contain the sum.
* @param mb The #multipole to add.
*/
void multipole_add(struct multipole *ma, const struct multipole *mb) {
#if const_gravity_multipole_order > 2
#error "Multipoles of order >2 not yet implemented."
#endif
/* Correct the position. */
const double ma_mass = ma->mass;
const double mb_mass = mb->mass;
const double M_tot = ma_mass + mb_mass;
const double M_tot_inv = 1.0 / M_tot;
const double ma_CoM[3] = {ma->CoM[0], ma->CoM[1], ma->CoM[2]};
const double mb_CoM[3] = {mb->CoM[0], mb->CoM[1], mb->CoM[2]};
#if const_gravity_multipole_order >= 2
const double ma_I_xx = (double)ma->I_xx + ma_mass * ma_CoM[0] * ma_CoM[0];
const double ma_I_yy = (double)ma->I_yy + ma_mass * ma_CoM[1] * ma_CoM[1];
const double ma_I_zz = (double)ma->I_zz + ma_mass * ma_CoM[2] * ma_CoM[2];
const double ma_I_xy = (double)ma->I_xy + ma_mass * ma_CoM[0] * ma_CoM[1];
const double ma_I_xz = (double)ma->I_xz + ma_mass * ma_CoM[0] * ma_CoM[2];
const double ma_I_yz = (double)ma->I_yz + ma_mass * ma_CoM[1] * ma_CoM[2];
const double mb_I_xx = (double)mb->I_xx + mb_mass * mb_CoM[0] * mb_CoM[0];
const double mb_I_yy = (double)mb->I_yy + mb_mass * mb_CoM[1] * mb_CoM[1];
const double mb_I_zz = (double)mb->I_zz + mb_mass * mb_CoM[2] * mb_CoM[2];
const double mb_I_xy = (double)mb->I_xy + mb_mass * mb_CoM[0] * mb_CoM[1];
const double mb_I_xz = (double)mb->I_xz + mb_mass * mb_CoM[0] * mb_CoM[2];
const double mb_I_yz = (double)mb->I_yz + mb_mass * mb_CoM[1] * mb_CoM[2];
#endif
/* New mass */
ma->mass = M_tot;
/* New CoM */
ma->CoM[0] = (ma_CoM[0] * ma_mass + mb_CoM[0] * mb_mass) * M_tot_inv;
ma->CoM[1] = (ma_CoM[1] * ma_mass + mb_CoM[1] * mb_mass) * M_tot_inv;
ma->CoM[2] = (ma_CoM[2] * ma_mass + mb_CoM[2] * mb_mass) * M_tot_inv;
/* New quadrupole */
#if const_gravity_multipole_order >= 2
ma->I_xx = (ma_I_xx + mb_I_xx) - M_tot * ma->CoM[0] * ma->CoM[0];
ma->I_yy = (ma_I_yy + mb_I_yy) - M_tot * ma->CoM[1] * ma->CoM[1];
ma->I_zz = (ma_I_zz + mb_I_zz) - M_tot * ma->CoM[2] * ma->CoM[2];
ma->I_xy = (ma_I_xy + mb_I_xy) - M_tot * ma->CoM[0] * ma->CoM[1];
ma->I_xz = (ma_I_xz + mb_I_xz) - M_tot * ma->CoM[0] * ma->CoM[2];
ma->I_yz = (ma_I_yz + mb_I_yz) - M_tot * ma->CoM[1] * ma->CoM[2];
#endif
}
/**
* @brief Add a particle to the given multipole.
*
* @param m The #multipole.
* @param p The #gpart.
*/
void multipole_addpart(struct multipole *m, struct gpart *p) {
/* #if const_gravity_multipole_order == 1 */
/* /\* Correct the position. *\/ */
/* float mm = m->coeffs[0], mp = p->mass; */
/* float w = 1.0f / (mm + mp); */
/* for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + p->x[k] * mp) * w;
*/
/* /\* Add the particle to the moments. *\/ */
/* m->coeffs[0] = mm + mp; */
/* #else */
/* #error( "Multipoles of order %i not yet implemented." ,
* const_gravity_multipole_order )
*/
/* #endif */
}
/**
* @brief Add a group of particles to the given multipole.
*
* @param m The #multipole.
* @param p The #gpart array.
* @param N Number of parts to add.
*/
void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
/* #if const_gravity_multipole_order == 1 */
/* /\* Get the combined mass and positions. *\/ */
/* double xp[3] = {0.0, 0.0, 0.0}; */
/* float mp = 0.0f, w; */
/* for (int k = 0; k < N; k++) { */
/* w = p[k].mass; */
/* mp += w; */
/* xp[0] += p[k].x[0] * w; */
/* xp[1] += p[k].x[1] * w; */
/* xp[2] += p[k].x[2] * w; */
/* } */
/* /\* Correct the position. *\/ */
/* float mm = m->coeffs[0]; */
/* w = 1.0f / (mm + mp); */
/* for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + xp[k]) * w; */
/* /\* Add the particle to the moments. *\/ */
/* m->coeffs[0] = mm + mp; */
/* #else */
/* #error( "Multipoles of order %i not yet implemented." ,
* const_gravity_multipole_order )
*/
/* #endif */
}
src/multipole.h
View file @ 68ce4a36
......@@ -22,28 +22,286 @@
/* Some standard headers. */
#include <math.h>
#include <string.h>
/* Includes. */
#include "align.h"
#include "const.h"
#include "error.h"
#include "gravity_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
#include "minmax.h"
#include "part.h"
/* Multipole struct. */
#define multipole_align 128
/**
* @brief The multipole expansion of a mass distribution.
*/
struct multipole {
/* Multipole location. */
double CoM[3];
union {
/*! Linking pointer for "memory management". */
struct multipole *next;
/*! The actual content */
struct {
/*! Multipole mass */
float mass;
/*! Centre of mass of the matter dsitribution */
double CoM[3];
/*! Bulk velocity */
float vel[3];
};
};
} SWIFT_STRUCT_ALIGN;
struct acc_tensor {
double F_000;
};
struct pot_tensor {
double F_000;
};
struct field_tensors {
union {
/*! Linking pointer for "memory management". */
struct field_tensors *next;
/*! The actual content */
struct {
/*! Field tensor for acceleration along x */
struct acc_tensor a_x;
/*! Field tensor for acceleration along y */
struct acc_tensor a_y;
/*! Field tensor for acceleration along z */
struct acc_tensor a_z;
/*! Field tensor for the potential */
struct pot_tensor pot;
};
};
} SWIFT_STRUCT_ALIGN;
/**
* @brief Reset the data of a #multipole.
*
* @param m The #multipole.
*/
INLINE static void multipole_init(struct multipole *m) {
/* Just bzero the struct. */
bzero(m, sizeof(struct multipole));
}
/**
* @brief Prints the content of a #multipole to stdout.
*
* Note: Uses directly printf(), not a call to message().
*
* @param m The #multipole to print.
*/
INLINE static void 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);
printf("Vel= [%12.5e %12.5e %12.5e\n", m->vel[0], m->vel[1], m->vel[2]);
}
/**
* @brief Adds a #multipole to another one (i.e. does ma += mb).
*
* @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) {
const float mass = ma->mass + mb->mass;
const float imass = 1.f / mass;
/* Add the bulk velocities */
ma->vel[0] = (ma->vel[0] * ma->mass + mb->vel[0] * mb->mass) * imass;
ma->vel[1] = (ma->vel[1] * ma->mass + mb->vel[1] * mb->mass) * imass;
ma->vel[2] = (ma->vel[2] * ma->mass + mb->vel[2] * mb->mass) * imass;
/* Add the masses */
ma->mass = mass;
}
/**
* @brief Verifies whether two #multipole's are equal or not.
*
* @param ma The first #multipole.
* @param mb The second #multipole.
* @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) {
/* Check CoM */
if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) > tolerance)
return 0;
if (fabs(ma->CoM[1] - mb->CoM[1]) / fabs(ma->CoM[1] + mb->CoM[1]) > tolerance)
return 0;
if (fabs(ma->CoM[2] - mb->CoM[2]) / fabs(ma->CoM[2] + mb->CoM[2]) > tolerance)
return 0;
/* Check bulk velocity (if non-zero)*/
if (fabsf(ma->vel[0] + mb->vel[0]) > 0.f &&
fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
tolerance)
return 0;
if (fabsf(ma->vel[1] + mb->vel[1]) > 0.f &&
fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
tolerance)
return 0;
if (fabsf(ma->vel[2] + mb->vel[2]) > 0.f &&
fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
tolerance)
return 0;
/* Check mass */
if (fabsf(ma->mass - mb->mass) / fabsf(ma->mass + mb->mass) > tolerance)
return 0;
/* Multipole mass */
float mass;
/* All is good */
return 1;
}
/**
* @brief Applies the forces due to particles j onto particles i directly.
*
* @param gparts_i The #gpart to update.
* @param gcount_i The number of particles to update.
* @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) {}
/**
* @brief Constructs the #multipole of a bunch of particles around their
* centre of mass.
*
* Corresponds to equation (28c).
*
* @param m The #multipole (content will be overwritten).
* @param gparts The #gpart.
* @param gcount The number of particles.
*/
INLINE static void multipole_P2M(struct multipole *m,
const struct gpart *gparts, int gcount) {
#if const_gravity_multipole_order >= 2
/* Quadrupole terms */
float I_xx, I_yy, I_zz;
float I_xy, I_xz, I_yz;
#error "Implementation of P2M kernel missing for this order."
#endif
};
/* Temporary variables */
double mass = 0.0;
double com[3] = {0.0, 0.0, 0.0};
float vel[3] = {0.f, 0.f, 0.f};
/* Collect the particle data. */
for (int k = 0; k < gcount; k++) {
const float m = gparts[k].mass;
mass += m;
com[0] += gparts[k].x[0] * m;
com[1] += gparts[k].x[1] * m;
com[2] += gparts[k].x[2] * m;
vel[0] += gparts[k].v_full[0] * m;
vel[1] += gparts[k].v_full[1] * m;
vel[2] += gparts[k].v_full[2] * m;
}
const double imass = 1.0 / mass;
/* Store the data on the multipole. */
m->mass = mass;
m->CoM[0] = com[0] * imass;
m->CoM[1] = com[1] * imass;
m->CoM[2] = com[2] * imass;
m->vel[0] = vel[0] * imass;