Commit b3d0dc61 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added fields to the multipole structure to hold the acceleration tensors

parent 1eacbd41
......@@ -323,13 +323,14 @@ int main(int argc, char *argv[]) {
/* How large are the parts? */
if (myrank == 0) {
message("sizeof(part) is %4zi bytes.", sizeof(struct part));
message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
message("sizeof(task) is %4zi bytes.", sizeof(struct task));
message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
message("sizeof(part) is %4zi bytes.", sizeof(struct part));
message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
message("sizeof(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor));
message("sizeof(task) is %4zi bytes.", sizeof(struct task));
message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
}
/* Read the parameter file */
......
......@@ -1113,7 +1113,7 @@ int cell_are_neighbours(const struct cell *restrict ci,
void cell_check_multipole(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
struct multipole ma;
struct gravity_tensors ma;
const double tolerance = 1e-5; /* Relative */
/* First recurse */
......@@ -1127,12 +1127,12 @@ void cell_check_multipole(struct cell *c, void *data) {
multipole_P2M(&ma, c->gparts, c->gcount);
/* Now compare the multipole expansion */
if (!multipole_equal(&ma, c->multipole, tolerance)) {
if (!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);
multipole_print(&ma.m_pole);
message("Recursive multipole:");
multipole_print(c->multipole);
multipole_print(&c->multipole->m_pole);
error("Aborting");
}
}
......
......@@ -104,7 +104,7 @@ struct cell {
double h_max;
/*! This cell's multipole. */
struct multipole *multipole;
struct gravity_tensors *multipole;
/*! Linking pointer for "memory management". */
struct cell *next;
......
......@@ -36,31 +36,6 @@
#define multipole_align 128
/**
* @brief The multipole expansion of a mass distribution.
*/
struct multipole {
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;
......@@ -71,30 +46,43 @@ struct pot_tensor {
double F_000;
};
struct field_tensors {
struct multipole {
union {
float mass;
/*! Linking pointer for "memory management". */
struct field_tensors *next;
/*! Bulk velocity */
float vel[3];
};
/*! The actual content */
struct {
/**
* @brief The multipole expansion of a mass distribution.
*/
struct gravity_tensors {
/*! Field tensor for acceleration along x */
struct acc_tensor a_x;
/*! Linking pointer for "memory management". */
struct gravity_tensors *next;
/*! Field tensor for acceleration along y */
struct acc_tensor a_y;
/*! Centre of mass of the matter dsitribution */
double CoM[3];
/*! Field tensor for acceleration along z */
struct acc_tensor a_z;
/*! The actual content */
struct {
/*! Field tensor for the potential */
struct pot_tensor pot;
};
};
/*! Multipole mass */
struct multipole m_pole;
/*! 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;
/**
......@@ -102,10 +90,18 @@ struct field_tensors {
*
* @param m The #multipole.
*/
INLINE static void multipole_init(struct multipole *m) {
INLINE static void multipole_reset(struct gravity_tensors *m) {
/* Just bzero the struct. */
bzero(m, sizeof(struct multipole));
bzero(m, sizeof(struct gravity_tensors));
}
INLINE static void multipole_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));
}
/**
......@@ -117,7 +113,7 @@ INLINE static void multipole_init(struct multipole *m) {
*/
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("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]);
}
......@@ -156,12 +152,15 @@ INLINE static int multipole_equal(const struct multipole *ma,
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;
/* 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 &&
......@@ -191,12 +190,12 @@ 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 multipole *m, double dt) {
INLINE static void multipole_drift(struct gravity_tensors *m, double dt) {
/* Move the whole thing according to bulk motion */
m->CoM[0] += m->vel[0];
m->CoM[1] += m->vel[1];
m->CoM[2] += m->vel[2];
m->CoM[0] += m->m_pole.vel[0];
m->CoM[1] += m->m_pole.vel[1];
m->CoM[2] += m->m_pole.vel[2];
}
/**
......@@ -220,7 +219,7 @@ 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 multipole *m,
INLINE static void multipole_P2M(struct gravity_tensors *m,
const struct gpart *gparts, int gcount) {
#if const_gravity_multipole_order >= 2
......@@ -248,13 +247,13 @@ INLINE static void multipole_P2M(struct multipole *m,
const double imass = 1.0 / mass;
/* Store the data on the multipole. */
m->mass = mass;
m->m_pole.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;
m->vel[1] = vel[1] * imass;
m->vel[2] = vel[2] * imass;
m->m_pole.vel[0] = vel[0] * imass;
m->m_pole.vel[1] = vel[1] * imass;
m->m_pole.vel[2] = vel[2] * imass;
}
/**
......@@ -290,38 +289,39 @@ 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 field_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) { */
/* dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.); */
/* dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.); */
/* dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.); */
/* } else { */
/* dx = pos_a[0] - pos_b[0]; */
/* dy = pos_a[1] - pos_b[1]; */
/* dz = pos_a[2] - pos_b[2]; */
/* } */
/* const double r2 = dx * dx + dy * dy + dz * dz; */
/* const double r_inv = 1. / sqrt(r2); */
/* /\* 1st order multipole term *\/ */
/* l_a->x.F_000 = D_100(dx, dy, dz, r_inv); */
/* l_a->y.F_000 = D_010(dx, dy, dz, r_inv); */
/* l_a->z.F_000 = D_001(dx, dy, dz, r_inv); */
}
/* INLINE static void multipole_M2L(struct field_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) { *\/ */
/* /\* dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.); *\/ */
/* /\* dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.); *\/ */
/* /\* dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.); *\/ */
/* /\* } else { *\/ */
/* /\* dx = pos_a[0] - pos_b[0]; *\/ */
/* /\* dy = pos_a[1] - pos_b[1]; *\/ */
/* /\* dz = pos_a[2] - pos_b[2]; *\/ */
/* /\* } *\/ */
/* /\* const double r2 = dx * dx + dy * dy + dz * dz; *\/ */
/* /\* const double r_inv = 1. / sqrt(r2); *\/ */
/* /\* /\\* 1st order multipole term *\\/ *\/ */
/* /\* l_a->x.F_000 = D_100(dx, dy, dz, r_inv); *\/ */
/* /\* l_a->y.F_000 = D_010(dx, dy, dz, r_inv); *\/ */
/* /\* l_a->z.F_000 = D_001(dx, dy, dz, r_inv); *\/ */
/* } */
#if 0
/* Multipole function prototypes. */
void multipole_add(struct multipole *m_sum, const struct multipole *m_term);
void multipole_init(struct multipole *m, const struct gpart *gparts,
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 multipole *m);
void multipole_reset(struct gravity_tensors *m);
/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
/* double *shift); */
......@@ -336,7 +336,7 @@ void multipole_reset(struct multipole *m);
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
struct multipole *ma, struct multipole *mb, double *shift) {
struct gravity_tensors *ma, struct gravity_tensors *mb, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
......@@ -378,7 +378,7 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mm(
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
struct multipole *m, struct gpart *p, double *shift) {
struct gravity_tensors *m, struct gpart *p, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
......
......@@ -267,8 +267,9 @@ void part_create_mpi_types() {
MPI_Type_commit(&spart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for sparts.");
}
if (MPI_Type_contiguous(sizeof(struct multipole) / sizeof(unsigned char),
MPI_BYTE, &multipole_mpi_type) != MPI_SUCCESS ||
if (MPI_Type_contiguous(
sizeof(struct gravity_tensors) / sizeof(unsigned char), MPI_BYTE,
&multipole_mpi_type) != MPI_SUCCESS ||
MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for multipole.");
}
......
......@@ -495,6 +495,9 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
/* Reset the gravity acceleration tensors */
if (c->multipole) multipole_init(c->multipole);
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
......
......@@ -68,7 +68,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const struct engine *e = r->e;
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
const struct multipole *multi = cj->multipole;
const struct gravity_tensors *multi = cj->multipole;
const float a_smooth = e->gravity_properties->a_smooth;
const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
......@@ -77,7 +77,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
#ifdef SWIFT_DEBUG_CHECKS
if (gcount == 0) error("Empty cell!");
if (multi->mass == 0.0) error("Multipole does not seem to have been set.");
if (multi->m_pole.mass == 0.0)
error("Multipole does not seem to have been set.");
#endif
/* Anything to do here? */
......@@ -111,10 +112,10 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact !*/
runner_iact_grav_pm(rlr_inv, r2, dx, gp, multi);
runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi->m_pole);
#ifdef SWIFT_DEBUG_CHECKS
gp->mass_interacted += multi->mass;
gp->mass_interacted += multi->m_pole.mass;
#endif
}
......
......@@ -186,8 +186,8 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
struct cell **cell_rec_begin,
struct cell **cell_rec_end,
struct multipole **multipole_rec_begin,
struct multipole **multipole_rec_end) {
struct gravity_tensors **multipole_rec_begin,
struct gravity_tensors **multipole_rec_end) {
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
......@@ -221,7 +221,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
for (int k = 0; k < num_elements; k++) {
struct cell *c = &cells[k];
struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
struct multipole *multipole_rec_begin = NULL, *multipole_rec_end = NULL;
struct gravity_tensors *multipole_rec_begin = NULL,
*multipole_rec_end = NULL;
space_rebuild_recycle_rec(s, c, &cell_rec_begin, &cell_rec_end,
&multipole_rec_begin, &multipole_rec_end);
if (cell_rec_begin != NULL)
......@@ -412,9 +413,9 @@ void space_regrid(struct space *s, int verbose) {
/* Allocate the multipoles for the top-level cells. */
if (s->gravity) {
if (posix_memalign((void *)&s->multipoles_top, multipole_align,
s->nr_cells * sizeof(struct multipole)) != 0)
s->nr_cells * sizeof(struct gravity_tensors)) != 0)
error("Failed to allocate top-level multipoles.");
bzero(s->multipoles_top, s->nr_cells * sizeof(struct multipole));
bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
}
/* Set the cells' locks */
......@@ -2090,7 +2091,7 @@ void space_split_recursive(struct space *s, struct cell *c,
if (s->gravity) {
/* Reset everything */
multipole_init(c->multipole);
multipole_reset(c->multipole);
/* Compute CoM of all progenies */
double CoM[3] = {0., 0., 0.};
......@@ -2098,11 +2099,11 @@ void space_split_recursive(struct space *s, struct cell *c,
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct multipole *m = c->progeny[k]->multipole;
CoM[0] += m->CoM[0] * m->mass;
CoM[1] += m->CoM[1] * m->mass;
CoM[2] += m->CoM[2] * m->mass;
mass += m->mass;
const struct gravity_tensors *m = c->progeny[k]->multipole;
CoM[0] += m->CoM[0] * m->m_pole.mass;
CoM[1] += m->CoM[1] * m->m_pole.mass;
CoM[2] += m->CoM[2] * m->m_pole.mass;
mass += m->m_pole.mass;
}
}
c->multipole->CoM[0] = CoM[0] / mass;
......@@ -2113,9 +2114,11 @@ void space_split_recursive(struct space *s, struct cell *c,
struct multipole temp;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct multipole *m = c->progeny[k]->multipole;
multipole_M2M(&temp, m, c->multipole->CoM, m->CoM, s->periodic);
multipole_add(c->multipole, &temp);
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);
}
}
}
......@@ -2284,8 +2287,8 @@ void space_recycle(struct space *s, struct cell *c) {
*/
void space_recycle_list(struct space *s, struct cell *cell_list_begin,
struct cell *cell_list_end,
struct multipole *multipole_list_begin,
struct multipole *multipole_list_end) {
struct gravity_tensors *multipole_list_begin,
struct gravity_tensors *multipole_list_end) {
int count = 0;
......@@ -2354,8 +2357,9 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* Is the multipole buffer empty? */
if (s->gravity && s->multipoles_sub == NULL) {
if (posix_memalign((void *)&s->multipoles_sub, multipole_align,
space_cellallocchunk * sizeof(struct multipole)) != 0)
if (posix_memalign(
(void *)&s->multipoles_sub, multipole_align,
space_cellallocchunk * sizeof(struct gravity_tensors)) != 0)
error("Failed to allocate more multipoles.");
/* Constructed a linked list */
......@@ -2381,7 +2385,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* Init some things in the cell we just got. */
for (int j = 0; j < nr_cells; j++) {
struct multipole *temp = cells[j]->multipole;
struct gravity_tensors *temp = cells[j]->multipole;
bzero(cells[j], sizeof(struct cell));
cells[j]->multipole = temp;
cells[j]->nodeID = -1;
......
......@@ -106,10 +106,10 @@ struct space {
struct cell *cells_sub;
/*! The multipoles associated with the top-level (level 0) cells */
struct multipole *multipoles_top;
struct gravity_tensors *multipoles_top;
/*! Buffer of unused multipoles for the sub-cells. */
struct multipole *multipoles_sub;
struct gravity_tensors *multipoles_sub;
/*! The total number of parts in the space. */
size_t nr_parts, size_parts;
......@@ -197,8 +197,8 @@ void space_rebuild(struct space *s, int verbose);
void space_recycle(struct space *s, struct cell *c);
void space_recycle_list(struct space *s, struct cell *cell_list_begin,
struct cell *cell_list_end,
struct multipole *multipole_list_begin,
struct multipole *multipole_list_end);
struct gravity_tensors *multipole_list_begin,
struct gravity_tensors *multipole_list_end);
void space_split(struct space *s, struct cell *cells, int nr_cells,
int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
......
Markdown is supported
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