GK passive advection & twistshift fixes

[manauref]

This goes with:

• PR 105 of gkylcas.
• PR 224 of postgkyl.

Changes

We add a new passive advection collisionless equation type, GKYL_GK_COLLISIONLESS_PASSIVE. This can be invoked from the input file along with user-specified advection speeds. See rt_gk_passive_3x2v_p1.c, for example, where we use

    .collisionless = {
      .type = GKYL_GK_COLLISIONLESS_PASSIVE,
      .passive_speeds = passive_velocity_elc,
      .passive_speeds_ctx = &ctx,
    },

and passive_velocity_elc is a function that returns 3 values for the x,y,z components of the advection velocity.

This object has been quite useful for debugging some twist-shift BC issues. For example, we can now reproduce the test in figure 14 of M.Francisquez, et al. CPC 298 (2024) 109109:

and now there's a unit test with this set up.

These kinds of tests allowed us to make 3 improvements in twistshift:

1. There were a few subcell integrals that were not checking if the subcell region had zero volume, which was leading to NaNs for the grid and shift in the published passive advection case. Zero volume subcell regions can happen due to floating-point precision, but should not contribute anything.
2. There were some cases in which points defining the subcell region were attributed to different periodic copies of the domain, leading to checkerboarding. This was fixed by ensuring subcell regions are always attributed to the same periodic copy of the domain.
3. The twistshift updater did not complain when the shift wrapped around Ly in a single x-cell. But it was not designed to support this, and it was leading to silent checkerboarding. Now it complains and prompts the user to change q, Ly, or Nx.

However, because of 3 now the d3d and tcv IWL reg tests don't run (because they are too coarse in x). Maybe we should increase the resolution of those.

Other small but consequential changes:

• We fixed the CBC regression tests, they were setting the y domain incorrectly, leading to erroneous shifts at the twist-shift boundary.
• We fixed a small bug in gkyl_rect_grid_find_cell. Parts of it assumed the indexing of cells went from 1 to N, but in one place it placed the lower cell at 0.

Additional changes:

• We renamed gkyl_bc_twistshift_new to gkyl_bc_twistshift_inew, and implemented a new gkyl_bc_twistshift_new which takes the arguments in the function signature, not in an input struct.
• We also changed some inputs to ``gkyl_bc_twistshift` so that they are passed by address.

These changes make it easier to call bc_twistshift from python/postgkyl.

Results

In addition to reproducing the passive advection test in the twistshift paper, we were also able to check advection across the twistshift boundary in the CBC geometry (see rt_gk_cbc_passive_3x2v_p1.c). Here's a movie made with the gk-rz command in PR 214 of postgkyl:

rt_gk_cbc_passive_3x2v_p1-elc_M0.1.mp4

Moreover, with the changes made to the bc_twistshift API, we created a Marimo notebook that can plot the effect of the shift for a specific input file (see postgkyl PR 224). It allows us to visualize the shifted lower-y cell, e.g.

and a specified quantity, twistshifted by the bc_twistshift operator:

Lastly, the pre-existent bc_twistshift unit tests continue to pass.

[Antoinehoff]

Great addition and fixes. Can you share your results of production-like simulation with that branch? Once this is done I think we can merge.

[manauref]

Great addition and fixes. Can you share your results of production-like simulation with that branch? Once this is done I think we can merge.

I checked a production TCV and D3D setup, and for those resolutions and shifts/q-profiles, nothing seemed to change. The TS fixes here appears to have afflicted other scenarios, but not those two.

I'll valgrind check when Perlmutter comes back online.

[manauref]

• the TS BC unit test is valgrind and compute-sanitizer clean.
• the d3d and tcv IWL reg tests are compute-sanitizer clean.
• the d3d and tcv IWL reg test is valgrind clean.

[Antoinehoff]

thanks !