[DR] Refactor the handling of geometric sources in GR (and other things)

Currently, source terms in the moment app are handled through the gkyl_moment_em_coupling struct, which populates a forcing vector based on specified couplings between (multiple) fluid species and electromagnetic fields. However, this design is clearly inadequate for various projects that we are currently pursuing and intend to pursue in the near future, since in many cases (such as general relativity), there may exist geometric sources that are inherent to a single fluid species, even in the absence of any other fluids or fields that it is coupled to. The same would also be true if we ever intend to support reactive fluid species, viscous hydrodynamics, etc. in Gkeyll.

One possible design, which was initially discussed, would be to implement a gkyl_moment_spacetime_coupling, in the same kind of vein as the existing gkyl_moment_em_coupling. However, it seems to me that a more minimalistic and elegant formulation (not to mention more general) would be to introduce a gkyl_<species>_source function that computes the forcing vector for a given fluid species, in the same way that we currently have a gkyl_<species>_flux function for computing that fluid's flux vector. For all fluid species currently implemented (with the exception of GR Euler), which are governed by homogeneous conservation laws in the absence of other fluids and fields, this forcing vector would default to being the zero vector. Then, we could introduce a pair of new function pointers to the gkyl_wv_eqn object: the first called source_func (which would default to being a "do nothing" function) which solves the forcing part of the problem using whatever combination of explicit SSP RK3, semi-implicit, fully implicit, etc. is desired; the second called split_func (which would simply default to calling the hyperbolic updater in the homogeneous case) which performs whatever kind of operator splitting is desired (Strang splitting, Godunov splitting, etc., in which the hyperbolic updater and the source updater specified via the source_func pointer may be called in any order, and with any specified combination of time-steps). This would necessitate adding only two additional lines in the gkyl_wv_<species>_inew method for each fluid species (which, as mentioned, could be populated by default so that no existing code would need to be changed):

struct gkyl_wv_eqn*
gkyl_wv_<species>_inew(const struct gkyl_wv_<species>_inp *inp)
{
  struct wv_<species> *<species> = gkyl_malloc(sizeof(struct wv_<species>));

  ...

  <species>->eqn.source_func = source_func;
  <species>->eqn.update_func = update_func;

  ...

  return &<species>->eqn;
}

For the time being, the GR Euler equations are the only ones that will have any non-default values for these function pointers (and all other fluid species' implementations may therefore remain unchanged), but as discussed with @JunoRavin on Friday, we may also wish to add non-default values for (e.g.) the 10-moment equations to facilitate some of the work that Jada is doing on the expanding box formalism. This will also be needed very soon to handle geometric sources appearing in the curved spacetime Maxwell equations from the work that Kiran is doing, and in several other places.

My hope is that these modifications will tesselate nicely with the changes being made to the handling of sources in the Vlasov and PKPM apps proposed (and implemented) as part of Design Review #386 - specifically the decision to pass an update_func function pointer as part of the app's input specification to specify which source solver to use - and pave the way to having a more-or-less completely unified design for the handling of source terms between all four apps.

ammarhakim / gkylzero

[DR] Refactor the handling of geometric sources in GR (and other things) #390