Mistake in short range gravity task creation for non-periodic case

[yuyttenhove]

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

[MatthieuSchaller]

@stuartmcalpine will be interested.

[MatthieuSchaller]

For the remaining cases you mention at the end, could you let me know what the size of the top-level grid is?

[yuyttenhove]

There is at least a 2D case with a 6x6 top level grid.

The problem seems to be that some pairs of cells are too far for the short range gravity tasks to be constructed in engine_make_self_gravity_tasks_mapper while the gravity_M2L_accept (which uses the not-purely geometric criterion for my configuration) also fails.

I do not now where the factor 2.5 in this code bloc comes from, but simply changing it to 3 (so that more pairs are considered to create the short range gravity interaction tasks for) seems to fix the issue...

I seem to remember the issue also popped up in at least one 3D case before my holiday break, but unfortunately I didn't document it and I have not found another example (yet), but I will keep you updated.

[MatthieuSchaller]

I agree with the changes you proposed above.

I am confused by this last comment. Are you running 2D gravity?

[yuyttenhove]

I was running some 2D simulations with gravity, yes. Not because I'm interested in the results, but because I had to make some changes to make the gravity work with zero mass particles (used to represent a vacuum in the moving mesh hydro code) and just wanted a quick way to check if there were no crashes and if the right interactions happened...

[MatthieuSchaller]

Right, I have no real intention to support 2D gravity... I am even surprised that it does not collapse completely into hell. :)

[MatthieuSchaller]

If you have some 3D cases that break then I am still very very interested. :) I have been banging my head for a while now on something that looks similar but can't find the cause, so additional broken examples can maybe help.

Your fixes to the non-periodic case are good. I'll merge them.