Commit 0861cdc1 authored by Matthieu Schaller's avatar Matthieu Schaller Committed by Matthieu Schaller
Browse files

Extract the mesh gravity forces to a separate set of functions. Get the forces...

Extract the mesh gravity forces to a separate set of functions. Get the forces directly from the mesh to the gparts. Do not use the tree here.
parent eb044806
......@@ -41,6 +41,7 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
mesh_side_length: 32
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.85 # Opening angle (Multipole acceptance criterion)
comoving_softening: 0.0026994 # Comoving softening length (in internal units).
......
......@@ -130,6 +130,7 @@ int main(int argc, char *argv[]) {
struct cooling_function_data cooling_func;
struct cosmology cosmo;
struct external_potential potential;
struct pm_mesh mesh;
struct gpart *gparts = NULL;
struct gravity_props gravity_properties;
struct hydro_props hydro_properties;
......@@ -688,13 +689,6 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
#ifndef HAVE_FFTW
/* Need the FFTW library if periodic and self gravity. */
if (with_self_gravity && periodic)
error(
"No FFTW library found. Cannot compute periodic long-range forces.");
#endif
#ifdef SWIFT_DEBUG_CHECKS
/* Check once and for all that we don't have unwanted links */
if (!with_stars && !dry_run) {
......@@ -733,6 +727,17 @@ int main(int argc, char *argv[]) {
/* Verify that the fields to dump actually exist */
if (myrank == 0) io_check_output_fields(params, N_total);
/* Initialise the long-range gravity mesh */
if (with_self_gravity && periodic) {
#ifdef HAVE_FFTW
pm_mesh_init(&mesh, &gravity_properties, dim[0]);
#else
/* Need the FFTW library if periodic and self gravity. */
error(
"No FFTW library found. Cannot compute periodic long-range forces.");
#endif
}
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
......@@ -827,7 +832,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &gravity_properties, &potential,
&hydro_properties, &gravity_properties, &mesh, &potential,
&cooling_func, &chemistry, &sourceterms);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
......@@ -1096,6 +1101,7 @@ int main(int argc, char *argv[]) {
/* Clean everything */
if (with_verbose_timers) timers_close_file();
if (with_cosmology) cosmology_clean(&cosmo);
if (with_self_gravity) pm_mesh_clean(&mesh);
engine_clean(&e);
free(params);
......
......@@ -47,7 +47,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.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 gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
mesh_gravity.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......@@ -59,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c \
collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c
chemistry.c cosmology.c restart.c mesh_gravity.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
......@@ -261,6 +261,9 @@ struct cell {
/*! Task computing long range non-periodic gravity interactions */
struct task *grav_long_range;
/*! Task propagating the mesh forces to the particles */
struct task *grav_mesh;
/*! Task propagating the multipole to the particles */
struct task *grav_down;
......
......@@ -333,8 +333,15 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
c->grav_down = scheduler_addtask(s, task_type_grav_down,
task_subtype_none, 0, 0, c, NULL);
/* Gravity mesh force propagation */
if (periodic)
c->grav_mesh = scheduler_addtask(s, task_type_grav_mesh,
task_subtype_none, 0, 0, c, NULL);
if (periodic) scheduler_addunlock(s, c->init_grav, c->grav_ghost_in);
if (periodic) scheduler_addunlock(s, c->grav_ghost_out, c->grav_down);
if (periodic) scheduler_addunlock(s, c->drift_gpart, c->grav_mesh);
if (periodic) scheduler_addunlock(s, c->grav_mesh, c->super->end_force);
scheduler_addunlock(s, c->init_grav, c->grav_long_range);
scheduler_addunlock(s, c->grav_long_range, c->grav_down);
scheduler_addunlock(s, c->grav_down, c->super->end_force);
......@@ -3810,6 +3817,10 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
/* Re-build the space. */
space_rebuild(e->s, e->verbose);
/* Re-compute the mesh forces */
if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
pm_mesh_compute_potential(e->mesh, e);
/* Re-compute the maximal RMS displacement constraint */
if (e->policy & engine_policy_cosmology)
engine_recompute_displacement_constraint(e);
......@@ -5427,7 +5438,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, const struct hydro_props *hydro,
struct gravity_props *gravity,
struct gravity_props *gravity, struct pm_mesh *mesh,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_global_data *chemistry,
......@@ -5490,6 +5501,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->cosmology = cosmo;
e->hydro_properties = hydro;
e->gravity_properties = gravity;
e->mesh = mesh;
e->external_potential = potential;
e->cooling_func = cooling_func;
e->chemistry = chemistry;
......
......@@ -39,6 +39,7 @@
#include "collectgroup.h"
#include "cooling_struct.h"
#include "gravity_properties.h"
#include "mesh_gravity.h"
#include "parser.h"
#include "partition.h"
#include "potential.h"
......@@ -302,6 +303,9 @@ struct engine {
/* Properties of the self-gravity scheme */
struct gravity_props *gravity_properties;
/* The mesh used for long-range gravity forces */
struct pm_mesh *mesh;
/* Properties of external gravitational potential */
const struct external_potential *external_potential;
......@@ -360,7 +364,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, const struct hydro_props *hydro,
struct gravity_props *gravity,
struct gravity_props *gravity, struct pm_mesh *mesh,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_global_data *chemistry,
......
......@@ -49,6 +49,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
error("Invalid tree rebuild frequency. Must be in [0., 1.]");
/* Tree-PM parameters */
p->mesh_size = parser_get_param_int(params, "Gravity:mesh_side_length");
p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
gravity_props_default_a_smooth);
p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max",
......@@ -120,6 +121,7 @@ void gravity_props_print(const struct gravity_props *p) {
p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent,
p->epsilon_max_physical);
message("Self-gravity mesh side-length: N=%d", p->mesh_size);
message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max);
......
......@@ -39,6 +39,9 @@ struct gravity_props {
/*! Frequency of tree-rebuild in units of #gpart updates. */
float rebuild_frequency;
/*! Periodic long-range mesh side-length */
int mesh_size;
/*! Mesh smoothing scale in units of top-level cell size */
float a_smooth;
......
......@@ -177,6 +177,41 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_dograv_external);
}
/**
* @brief Calculate gravity accelerations from the periodic mesh
*
* @param r runner task
* @param c cell
* @param timer 1 if the time is to be recorded.
*/
void runner_do_grav_mesh(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gparts = c->gparts;
const int gcount = c->gcount;
const struct engine *e = r->e;
#ifdef SWIFT_DEBUG_CHECKS
if (!e->s->periodic) error("Calling mesh forces in non-periodic mode.");
#endif
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_gravity(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_grav_mesh(r, c->progeny[k], 0);
} else {
/* Get the forces from the gravity mesh */
pm_mesh_interpolate_forces(e->mesh, e, gparts, gcount);
}
if (timer) TIMER_TOC(timer_dograv_mesh);
}
/**
* @brief Calculate change in thermal state of particles induced
* by radiative cooling and heating.
......@@ -2169,8 +2204,11 @@ void *runner_main(void *data) {
case task_type_grav_down:
runner_do_grav_down(r, t->ci, 1);
break;
case task_type_grav_mesh:
runner_do_grav_mesh(r, t->ci, 1);
break;
case task_type_grav_top_level:
runner_do_grav_fft(r, 1);
// runner_do_grav_fft(r, 1);
break;
case task_type_grav_long_range:
runner_do_grav_long_range(r, t->ci, 1);
......
......@@ -47,6 +47,7 @@
#include "lock.h"
#include "logger.h"
#include "map.h"
#include "mesh_gravity.h"
#include "multipole.h"
#include "parallel_io.h"
#include "parser.h"
......
......@@ -56,7 +56,8 @@ const char *taskID_names[task_type_count] = {
"kick2", "timestep", "send",
"recv", "grav_top_level", "grav_long_range",
"grav_ghost_in", "grav_ghost_out", "grav_mm",
"grav_down", "cooling", "sourceterms"};
"grav_down", "grav_mesh", "cooling",
"sourceterms"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
......
......@@ -64,6 +64,7 @@ enum task_types {
task_type_grav_ghost_out,
task_type_grav_mm,
task_type_grav_down,
task_type_grav_mesh,
task_type_cooling,
task_type_sourceterms,
task_type_count
......
......@@ -60,6 +60,7 @@ const char* timers_names[timer_count] = {
"dopair_grav_pp",
"dograv_external",
"dograv_down",
"dograv_mesh",
"dograv_top_level",
"dograv_long_range",
"dosource",
......
......@@ -61,6 +61,7 @@ enum {
timer_dopair_grav_pp,
timer_dograv_external,
timer_dograv_down,
timer_dograv_mesh,
timer_dograv_top_level,
timer_dograv_long_range,
timer_dosource,
......
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