Interacting unsorted cells about swift HOT 15 CLOSED

Could you describe how the decoupling is implemented?

Also, how did you guaranty that the BH action takes place after the stellar properties are computed?

from swift.

In hydro_iact.h for SPHENIX we have:

if (pi->feedback_data.decoupling_delay_time <= 0.f && pj->feedback_data.decoupling_delay_time <= 0.f) {
      pi->a_hydro[0] -= mj * acc * dx[0];
      pi->a_hydro[1] -= mj * acc * dx[1];
      pi->a_hydro[2] -= mj * acc * dx[2];`
      
     pj->a_hydro[0] += mi * acc * dx[0];
      pj->a_hydro[1] += mi * acc * dx[1];
      pj->a_hydro[2] += mi * acc * dx[2];
  }

The decoupling_delay_time parameter is set during a stellar feedback/black hole feedback kick event to some value. In this function we decrease the delay time until it reaches <0:

__attribute__((always_inline)) INLINE static void feedback_update_part(
    struct part* p, struct xpart* xp, const struct engine* e, 
    const int with_cosmology) {

  /* No reason to do this is the decoupling time is zero */
  if (p->feedback_data.decoupling_delay_time > 0.f) {
    const integertime_t ti_step = get_integer_timestep(p->time_bin);
    const integertime_t ti_begin =
        get_integer_time_begin(e->ti_current - 1, p->time_bin);

    /* Get particle time-step */
    double dt_part;
    if (with_cosmology) {
      dt_part = cosmology_get_delta_time(e->cosmology, ti_begin,
                                          ti_begin + ti_step);
    } else {
      dt_part = get_timestep(p->time_bin, e->time_base);
    }

    p->feedback_data.decoupling_delay_time -= dt_part;
    if (p->feedback_data.decoupling_delay_time <= 0.f) {
      p->feedback_data.decoupling_delay_time = 0.f;
    }
  } else {
    p->feedback_data.decoupling_delay_time = 0.f;
  }
}

feedback_update_part is called right before drifting. The loop over star particles for each black hole is done in the black hole density loop (I hope?) by using the following:

#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
      /* Only do bh-stars loop if necessary for the model. */
      if (bh_stars_loop_is_active(bi, e)) {
        for (int sjd = 0; sjd < scount; sjd++) {
          struct spart *restrict sj = &sparts[sjd];

          /* Compute the pairwise distance. */
          const float sjx[3] = {(float)(sj->x[0] - c->loc[0]),
                                (float)(sj->x[1] - c->loc[1]),
                                (float)(sj->x[2] - c->loc[2])};
          const float dx[3] = {bix[0] - sjx[0], bix[1] - sjx[1], bix[2] - sjx[2]};
          const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];

          /* Star is within smoothing length of black hole */
          if (r2 < hig2) {
            runner_iact_nonsym_bh_stars_density(r2, dx, bi, sj);
          }
        }
        /* TODO: One possible speed-up is just to flag the BH id 
        * that each star is associated with in the previous loop,
        * and then just use that to loop instead of doing the distance
        * calculation everytime.
        */

        /* Now that we have the angular momentum, find the bulge mass */
        for (int sjd = 0; sjd < scount; sjd++) {
          struct spart *restrict sj = &sparts[sjd];

          /* Compute the pairwise distance. */
          const float sjx[3] = {(float)(sj->x[0] - c->loc[0]),
                                (float)(sj->x[1] - c->loc[1]),
                                (float)(sj->x[2] - c->loc[2])};
          const float dx[3] = {bix[0] - sjx[0], bix[1] - sjx[1], bix[2] - sjx[2]};
          const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];

          /* Star is within smoothing length of black hole */
          if (r2 < hig2) {
            runner_iact_nonsym_bh_stars_bulge(r2, dx, bi, sj);
          }
        }
      }
#endif

That loop is done right before the loop over all gas neighbours. That code is run in 3 different functions in the black hole functions runner:

DOSELF1_BH
DO_NONSYM_PAIR1_BH_NAIVE
DOSELF1_SUBSET_BH

from swift.

The BH bits seem sensible.

Why not also change DOPAIR1_SUBSET_BH_NAIVE?

The decoupling is not doing what is usually understood by decoupling in the literature I think. All you are doing here is not adding hydro accelerations to the particle. It will still participate in all the other loops and pieces of code. For instance in the density estimate or feedback.
I don't know whether this leads to issues however. Likely not but it means interactions are not consistent between the loops.

You could also try switching on all the internal consistency checks by configuring with --enable-debugging-checks.

from swift.

Thank you for noticing that, I will also add it into that function.

I looked deeper into the decoupling and noticed that it was NOT done in runner_iact_nonsym_force for nonsymmetric interactions. I changed it to make sure that the hydro accelerations and energy transfer were not done there. I do not get the error now, so hopefully that was the solution. I also made sure that all of the other sub-grid physics loops always ignore decoupled particles, just in case.

I will leave this open in case it pops up again.

It does seem like the code is sensitive to decoupling (as it should be!). One question came out of me digger deeper:

What about the time derivative of h?

p->force.h_dt

In the default decoupling implementation that I began with that was still being computed and updated even for decoupled particles. I ran a test where I did not update h_dt for particles interacting with a wind, and the simulation did not crash. That should be the correct thing to do, right?

Thanks.

from swift.

dh/dt is not super important as it is mainly used to get a good guess at h for the next step.

If you do ignore the decoupled particles in the density loop (as I think you should) then I would also not update dh/dt. Such that the estimate of h is (more) consistent with the density field.

Note that you should also not update du/dt in the force loops for the decoupled particles.

from swift.

Right, thanks for the info. Indeed we are not computing du/dt for decoupled particles.

I forgot to mention that I tried decoupling particles from the main hydro density loop. However, that led to serious issues with particles having zero neighbours and maximum softening lengths, and then to having a timestep of 0 and crashing swift.

It is possible for me to go through and simply ignore these errors until recoupling though, so I might test that out.

The tests I ran made it to z = 0 and nothing looks amiss, so for now I will leave it as a TODO.

from swift.

Ok. I can imagine that it will go wrong indeed. As mentioned above, there is a lot more that should be done to properly define decoupled particles in the way usually assumed.

But if this now works for you, then that's the main bit. Glad it's resolved. :)

from swift.

This could possibly be a compiler problem. I was seeing this with gcc/11.x compilers on some machines in Edinburgh, but I have not seen the problem since I moved to our Canadian cluster and have been using the latest Intel compilers. I haven't tested extensively, but will update if I get a chance. Definitely not fixed yet!

from swift.

That sounds very suspicious. If that is really what is going on then I'd rather think that something somewhere relies on not-well-defined behaviour or on some form of thread-race condition that different compilers implement in a different way. But you are then relying on bad coding to get things to work in the one compiler that seemingly does not crash the code.

from swift.

Indeed, this is just conjecture at this point. I will note I have not seen the problem on any system but the one in Edinburgh. I did run the EAGLE-XL model on the CCA cluster, and did not see the issue with a gcc/11.x compiler.

from swift.

One small note so that I don't forget. The implementation of the DM vel. disp. won't be robust over MPI. That's not a problem for now but if you intend to scale things up additional communications will be required in various places.

from swift.

Hi Doug,

We may have fixed something related to the original problem last week. I don't know whether it is actually the same situation you are facing here but there have been code logic changes around that section.
Note that we are still confused about the GCC vs. ICC difference.

May or may not help.

from swift.

Are you still encountering this issue?

from swift.

We have not seen this issue since last time we spoke.

from swift.

Let's re-open or create a new issue if that problem reappears then.

from swift.