Hi SWIFT team,
I have come accross the following in getting gravity to work with the moving mesh hydro scheme.
In the function engine_make_self_gravity_tasks_mapper
there is the following code block:
/* Special case where every cell is in range of every other one */
if (delta >= cdim[0] / 2) {
if (cdim[0] % 2 == 0) {
delta_m = cdim[0] / 2;
delta_p = cdim[0] / 2 - 1;
} else {
delta_m = cdim[0] / 2;
delta_p = cdim[0] / 2;
}
}
However, it is only true that every cell is in range of every other one if delta >= cdim[0] / 2
with periodic boundary conditions. This resulted in some particles not interacting with all the particles they are supposed to interact with. The following error is produced when I run a non-periodic simulation with self gravity with a very small IC of a grid of 15x15x15 particles all with the same non-zero mass:
[00002.1] runner_others.c:runner_do_end_grav_force():836: g-particle (id=1936, type=Gas) did not interact gravitationally with all other gparts gp->num_interacted=2250, total_gparts=3375 (local num_gparts=3375 inhibited_gparts=0)
I would suggest replacing the above code block: with the following to fix this:
/* Special case where every cell is in range of every other one */
if (periodic) {
if (delta >= cdim[0] / 2) {
if (cdim[0] % 2 == 0) {
delta_m = cdim[0] / 2;
delta_p = cdim[0] / 2 - 1;
} else {
delta_m = cdim[0] / 2;
delta_p = cdim[0] / 2;
}
}
} else {
if (delta > cdim[0]) {
delta_m = cdim[0];
delta_p = cdim[0];
}
}
On a related note:
I have also noticed that the function cell_min_dist2_same_size
(which is used in engine_make_self_gravity_tasks_mapper
) contains the following code for the non-periodic case:
const double dx = min(fabs(cix_max - cjx_min), fabs(cix_min - cjx_max));
const double dy = min(fabs(ciy_max - cjy_min), fabs(ciy_min - cjy_max));
const double dz = min(fabs(ciz_max - cjz_min), fabs(ciz_min - cjz_max));
I also think this is not true in general and should be replaced with:
const double dx = min4(fabs(cix_min - cjx_min), fabs(cix_min - cjx_max), fabs(cix_max - cjx_min), fabs(cix_max - cjx_max));
const double dy = min4(fabs(ciy_min - cjy_min), fabs(ciy_min - cjy_max), fabs(ciy_max - cjy_min), fabs(ciy_max - cjy_max));
const double dz = min4(fabs(ciz_min - cjz_min), fabs(ciz_min - cjz_max), fabs(ciz_max - cjz_min), fabs(ciz_max - cjz_max));
This then also follows the same structure of the periodic case. Since the cells are assumed to be the same size in this function, fabs(cix_min - cjx_min)
will be equal to fabs(cix_max - cjx_max)
(and similarly for y and z of course), so this could maybe be simplified a bit.
The above seems to have fixed g-particle did not interact gravitationally with all other gparts
in the majority of my test cases, but in some cases there seems to be one cell whose gparts do not interact with other gparts. I have yet to find the cause of that. I can provide more details about those cases here, but they are probably more related to this issue.
Thanks in advance for your input.
Yolan Uyttenhove