[DR] Implement full adaptivity for mesh refinement functionality

[JonathanGorard]

Thanks to the recent PR #365, all of the core infrastructure for handling static/fixed mesh refinement in the G0 moment app now exists in the /amr directory (e.g. methods for creating the relevant block topologies for representing the full mesh hierarchy; methods for applying inter-block boundary conditions, including the necessary projection and restriction operators to accommodate coarse-fine boundaries; methods for computing stable coarse and fine grid time-steps, and for performing hierarchical time-stepping in which the fine grids are evolved recursively to "catch up" with coarse grids; etc.). In general, we have followed the data structure and time-stepping algorithm design presented by Berger and Oliger for hyperbolic PDEs, with appropriate modifications to the finite-volume update algorithm in accordance with the approach of Berger and Colella for shock hydrodynamics (we also use their form of the coarse-to-fine projection operators and the fine-to-coarse restriction operators). Moreover, thanks to the inherent flexibility of the gkyl_block_topo data structure, we are effectively able to simulate both block-based and patch-based mesh refinement techniques with a single set of algorithms.

With this static functionality now in place, the only critical piece of infrastructure currently preventing us from progressing to doing fully adaptive mesh refinement (in which refined blocks/patches are recursively added or deleted in accordance with a predetermined refinement or tagging function) is the ability to perform dynamic regridding, in which fluid variables are dynamically projected downwards from coarse-to-fine, or restricted upwards from fine-to-coarse, mid-simulation. This regridding step requires one to perform copies of fluid variables both to-and-from buffers at run-time, where these copies are performed across entire blocks (i.e. not just in skin/ghost regions, as done currently), and must be performed every time the AMR hierarchy is modified - likely with a predefined frequency (e.g. every 100 time-steps). This, in turn, requires a small and hopefully reasonably conservative modification to the moment_update_one_step() method in apps/mom_update_one_step.c, hence this request for a design review. Since only one minor modification to the App code is being proposed, and therefore all of the relevant changes to the App code may be fit into the text of the Issue itself, I have not yet produced a prototype version of this modification; however, I will shortly begin working in a resurrected version of the amr-prototyping branch.

My proposed design is to introduce the following recursively-defined tree data structure (to be defined within zero/gkyl_block_topo.h, so as to prevent unnecessary proliferation of additional header files), for the purpose representing the AMR hierarchy in a manner which can easily be modified, with minimal reshuffling of memory, in order to add/remove refined regions during the regridding step as required:

struct gkyl_block_hierarchy {
  int bid;
  struct gkyl_block_hierarchy *parent;
  struct gkyl_block_hierarchy **child;
}

with bid being the ID of the block/patch, *parent being a pointer to the (coarser) parent level in the block hierarchy (to be set to null if we are already at the coarsest level), and **child being an array of pointers to the (finer) child blocks/patches in the block hierarchy (to be set to null if we are already at the finest level). We may also wish to accommodate the case in which a single block/patch has multiple parent blocks/patches (e.g. a finer grid is contained within two coarser grid patches), in which case we should consider having **parent instead be an array of pointers, just like **child. This depends upon how flexible we want our AMR topologies to be going forward

Then, in apps/mom_update_one_step.c, we would add an additional REGRIDDING time-stepper state to the state enum in moment_update_one_step(), thus yielding:

  enum {
    UPDATE_DONE = 0,
    PRE_UPDATE,
    POST_UPDATE,
    FIRST_COUPLING_UPDATE,
    FIELD_UPDATE,
    SPECIES_UPDATE,
    SECOND_COUPLING_UPDATE,
    REGRIDDING,
    UPDATE_REDO,
  } state = PRE_UPDATE;

along with appropriate modification to the switch statement:

      case POST_UPDATE:
        state = REGRIDDING;

        // copy solution in prep for next time-step
        for (int i=0; i<ns; ++i) {
          // check for nans before copying
          if (check_for_nans(app->species[i].f[ndim], app->local))
            have_nans_occured = true;
          else // only copy in case no nans, so old solution can be written out
            gkyl_array_copy(app->species[i].f[0], app->species[i].f[ndim]);
        }

        if (app->has_field)
          gkyl_array_copy(app->field.f[0], app->field.f[ndim]);

        break;

      case REGRIDDING:
        state = UPDATE_DONE;

        for (int i=0; i<ns; ++i) {
          // regridding logic for each species to go here
        }

        break;

We would also add a corresponding function pointer to the gkyl_moment_app struct in apps/gkyl_moment_priv.h to represent the regridding function - for any non-AMR simulation, this function pointer will default to null, and so the REGRIDDING time-stepper state will fall straight through to the UPDATE_DONE state.

Finally the input specification for an AMR simulation will also take in a function pointer representing the "tagging" or refinement criterion function (e.g. to specify whether you want to refine based upon the first derivative of density, second derivative of pressure, etc.), along with a tolerance level wherein, if the absolute value of the tagging function exceeds that tolerance value within some region of the domain, then it will flag that region of the domain for refinement (and will continue to do so recursively until the absolute value of the tagging function falls below that predefined tolerance value). However, since this change does not require modification to any code in either the /apps or /zero directories of G0, we will not worry too much about its precise design for now (not least because we fully expect the design of the current code in the /amr directory to change considerably over the coming months as this AMR capability is progressively expanded and refactored). The eventual aim will be to move the /amr code fully into the /apps directory, either to exist alongside the existing moment app time-stepper, or to supersede it entirely (thereby treating unigrid moment simulations as a special case of AMR moment simulations in which the grid hierarchy is trivial). But this is a lower-priority, longer-term, and significantly more complex design discussion, which we do not aim to have here.

[ammarhakim]

@JonathanGorard if the recent PR addressed this Issue, please close it as completed.