diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml index 54e56794c85c08ee3cf7f4a25e6e530ea434c77d..1ea1518825debe56cb8462c4a1b398c03c257bfe 100644 --- a/examples/EAGLE_100/eagle_100.yml +++ b/examples/EAGLE_100/eagle_100.yml @@ -28,8 +28,6 @@ Gravity: eta: 0.025 # Constant dimensionless multiplier for time integration. epsilon: 0.0001 # Softening length (in internal units). theta: 0.7 # Opening angle (Multipole acceptance criterion) - a_smooth: 1.25 # PM smoothing scale in top-level cell-size units - r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size units # Parameters for the hydrodynamics scheme SPH: diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml index 4bb04e757f8e43e9653e2dde066b8ca560956ce4..fd9569c50f4ad0d87eccd54e34636180c0337596 100644 --- a/examples/EAGLE_12/eagle_12.yml +++ b/examples/EAGLE_12/eagle_12.yml @@ -28,8 +28,6 @@ Gravity: eta: 0.025 # Constant dimensionless multiplier for time integration. epsilon: 0.0001 # Softening length (in internal units). theta: 0.7 # Opening angle (Multipole acceptance criterion) - a_smooth: 1.25 # PM smoothing scale in top-level cell-size units - r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size un # Parameters for the hydrodynamics scheme SPH: diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml index 8920fb89a9966d1de3ebbb0e3b043f4ae85df3f2..5dee9dad0b5d7f694c61fa4c983ead0f1cd6e5e2 100644 --- a/examples/EAGLE_25/eagle_25.yml +++ b/examples/EAGLE_25/eagle_25.yml @@ -28,8 +28,6 @@ Gravity: eta: 0.025 # Constant dimensionless multiplier for time integration. epsilon: 0.0001 # Softening length (in internal units). theta: 0.7 # Opening angle (Multipole acceptance criterion) - a_smooth: 1.25 # PM smoothing scale in top-level cell-size units - r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size units # Parameters for the hydrodynamics scheme SPH: diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml index b84b1eb7c362f85d8cd6a08ff2a15f72d1337396..898c28935abd02ec115ce107bdcfa4006c41dc48 100644 --- a/examples/EAGLE_50/eagle_50.yml +++ b/examples/EAGLE_50/eagle_50.yml @@ -23,6 +23,12 @@ Snapshots: Statistics: delta_time: 1e-2 # Time between statistics output +# Parameters for the self-gravity scheme +Gravity: + eta: 0.025 # Constant dimensionless multiplier for time integration. + epsilon: 0.0001 # Softening length (in internal units). + theta: 0.7 # Opening angle (Multipole acceptance criterion) + # Parameters for the hydrodynamics scheme SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). diff --git a/examples/UniformDMBox/uniformBox.yml b/examples/UniformDMBox/uniformBox.yml index dd8ea53a20b14e9fe7206e1ac692d16f6524976d..e59d677b308ca70f212f74c7e4d8b79f015c77a9 100644 --- a/examples/UniformDMBox/uniformBox.yml +++ b/examples/UniformDMBox/uniformBox.yml @@ -28,8 +28,7 @@ Gravity: eta: 0.025 # Constant dimensionless multiplier for time integration. theta: 0.7 # Opening angle (Multipole acceptance criterion) epsilon: 0.00001 # Softening length (in internal units). - a_smooth: 1.25 # PM smoothing scale in top-level cell-size units - r_cut: 4.5 # Distance beyong wwhich FMM forces are zero in top-level cell-size un + # Parameters governing the conserved quantities statistics Statistics: delta_time: 1e-2 # Time between statistics output diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 7ca15576bc9da20a858db6213bb9f79d5a880bf7..9c3cee7630edf1be1e161a3e70547f06e6108ebd 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -51,11 +51,12 @@ SPH: # Parameters for the self-gravity scheme Gravity: - eta: 0.025 # Constant dimensionless multiplier for time integration. - theta: 0.7 # Opening angle (Multipole acceptance criterion) - epsilon: 0.1 # Softening length (in internal units). - a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value). - r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value). + eta: 0.025 # Constant dimensionless multiplier for time integration. + theta: 0.7 # Opening angle (Multipole acceptance criterion) + epsilon: 0.1 # Softening length (in internal units). + a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value). + r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value). + r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value). # Parameters related to the initial conditions InitialConditions: diff --git a/src/gravity_properties.c b/src/gravity_properties.c index b1098888b96cdef2205ed513e60a3799c63e8b9f..18cf044434f7840a5a76f483540bb924a2365e26 100644 --- a/src/gravity_properties.c +++ b/src/gravity_properties.c @@ -33,7 +33,8 @@ #include "kernel_gravity.h" #define gravity_props_default_a_smooth 1.25f -#define gravity_props_default_r_cut 4.5f +#define gravity_props_default_r_cut_max 4.5f +#define gravity_props_default_r_cut_min 0.1f void gravity_props_init(struct gravity_props *p, const struct swift_params *params) { @@ -41,8 +42,10 @@ void gravity_props_init(struct gravity_props *p, /* Tree-PM parameters */ p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth", gravity_props_default_a_smooth); - p->r_cut = parser_get_opt_param_float(params, "Gravity:r_cut", - gravity_props_default_r_cut); + p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max", + gravity_props_default_r_cut_max); + p->r_cut_min = parser_get_opt_param_float(params, "Gravity:r_cut_min", + gravity_props_default_r_cut_min); /* Time integration */ p->eta = parser_get_param_float(params, "Gravity:eta"); @@ -69,9 +72,10 @@ void gravity_props_print(const struct gravity_props *p) { message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)", p->epsilon, p->epsilon / 3.); - message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth); + message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth); - message("Self-gravity MM cut-off: r_cut=%f", p->r_cut); + message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max); + message("Self-gravity truncation cut-off: r_cut_min=%f", p->r_cut_min); } #if defined(HAVE_HDF5) @@ -84,7 +88,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav, p->epsilon / 3.); io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit); io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER); - io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth); - io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut); + io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth); + io_write_attribute_f(h_grpgrav, "Mesh r_cut_max", p->r_cut_max); + io_write_attribute_f(h_grpgrav, "Mesh r_cut_min", p->r_cut_min); } #endif diff --git a/src/gravity_properties.h b/src/gravity_properties.h index be26f0d1d23b8cec71fa3cbbeedac9f61f337b2c..2a5e4cb1e07ea591e2e3821704ec55abe7980360 100644 --- a/src/gravity_properties.h +++ b/src/gravity_properties.h @@ -34,9 +34,16 @@ */ struct gravity_props { - /* Tree-PM parameters */ + /*! Mesh smoothing scale in units of top-level cell size */ float a_smooth; - float r_cut; + + /*! Distance below which the truncated mesh force is Newtonian in units of + * a_smooth */ + float r_cut_min; + + /*! Distance above which the truncated mesh force is negligible in units of + * a_smooth */ + float r_cut_max; /*! Time integration dimensionless multiplier */ float eta; diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 1de491fa01aa65026d377a6abf14aa1e7564a0d4..acf6478a09d6a4d9670ada8431fd288628296bf4 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -664,7 +664,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) { const double cell_width = s->width[0]; const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; const double theta_crit_inv = props->theta_crit_inv; - const double max_distance = props->a_smooth * props->r_cut * cell_width; + const double max_distance = props->a_smooth * props->r_cut_max * cell_width; const double max_distance2 = max_distance * max_distance; struct gravity_tensors *const mi = ci->multipole; const double CoM[3] = {mi->CoM[0], mi->CoM[1], mi->CoM[2]}; @@ -701,9 +701,11 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) { const double dz = nearest(CoM[2] - mj->CoM[2], dim[2]); const double r2 = dx * dx + dy * dy + dz * dz; + /* Are we beyond the distance where the truncated forces are 0 ?*/ if (r2 > max_distance2) { #ifdef SWIFT_DEBUG_CHECKS + /* Need to account for the interactions we missed */ mi->pot.num_interacted += mj->m_pole.num_gpart; #endif continue; diff --git a/src/tools.c b/src/tools.c index d701afc5fd1c81960ff014c7f45af2b8debd0123..0f423d5aad620edad3a430e7ea962218430da86d 100644 --- a/src/tools.c +++ b/src/tools.c @@ -755,7 +755,7 @@ void gravity_n2(struct gpart *gparts, const int gcount, const struct gravity_props *gravity_properties, float rlr) { const float rlr_inv = 1. / rlr; - const float r_cut = gravity_properties->r_cut; + const float r_cut = gravity_properties->r_cut_max; const float max_d = r_cut * rlr; const float max_d2 = max_d * max_d;