Add feature to automatically set reference energy for Davis Products EOS

[jhp-lanl]

The Davis Products EOS is most often used in combination with a reactants EOS to model reactive flows. For a detonation wave to be self-propagating, the Rayleigh line extending from the initial conditions of the reactants needs to be tangent to the release isentrope of the products at the CJ point. Most often, we want this CJ isentrope to the reference isentrope for the Davis EOS, so it makes sense to be able to automatically set the reference energy of the products EOS automatically in order to achieve this tangency condition.

This will require an iterative procedure at setup that will could look something like this:

1. Find bounding points on the reference isentrope where the slope of the Rayleigh line is larger than and smaller than the slope of the isentrope.
2. Use a bounded root finding method to iteratively find the point where the slope of the Rayleigh line equals the slope of the reference isentrope
3. Calculate the energy of the CJ point given the Rayleigh line and then subtract the energy of the reference isentrope at the CJ point to calculate the appropriate energy offset.

This procedure should be done with an arbitrary reactants reference energy.

[jhp-lanl]

This procedure requires two additional parameters: the reference pressure and reference density/volume for the reacting material. When these parameters are given, these would obviate the need for specifying E0.

Open question: should these parameters replace E0 (i.e. overload the constructor so that there is one more input) or should all three (E0, P0, and rho0) be optional parameters?

The former places the burden on the host code to determine what the user is trying to specify while the latter pushes that logic into singularity. I'm leaning towards overloading since I think that would make it more clear in the host code what singularity is doing under the hood.

[Yurlungur]

IMO different host codes may want to do things differently, so I am in favor of the former---the overload approach.

[jhp-lanl]

@Yurlungur brought up

👍 I'm in favor of the helper being static too.

A static member function would in addition to P0 and rho0 need to take vc, k, n, and a (i.e. all but two parameters for the EOS, excepting E0).

I suppose if this is all just part of the constructor though it shouldn't be an issue. Originally I was thinking maybe it shouldn't be static, but now that I think about it, I think I agree with you.

[jhp-lanl]

Now that the Davis Product EOS no longer takes E0, it seems prudent for this to NOT be static. That way the energy shift modifier can be used with something like

DavisProducts eos = DavisProducts(...);
eos = ShiftedEOS<DavisProducts>{eos, eos.CalcCJEnergyShift()};

The next question in my mind is how the Fortran interface should work. Within the current form it's a bit messy to change the modifier parameters within the functions, but it could look something like

int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b,
                          const double k, const double n, const double vc,
                          const double pc, const double Cv, int const *const enabled_in,
                          double *const vals) {
  assert(matindex >= 0);
  DavisProducts eos_base = DavisProducts{a, b, k, n, vc, pc, Cv};
  const int *const enabled[3] = {enabled_in[0], 1, enabled_in[2]};
  const double eshiftCJ = eos_base.CalcCJEnergyShift();
  if (enabled_in[1]) {
    vals[1] += eshiftCJ;
  } else {
    vals[1] = eshiftCJ;
  }
  EOS eosi = SGAPPLYMODSIMPLE(eos_base);
  if (enabled[3] == 1) {
    singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2],
                                           vals[3], vals[4], vals[5]);
  }
  EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv));
  eos[matindex] = eos_.GetOnDevice();
  return 0;
}

It would also be nice to move the enabled/values arrays to structs with corresponding fortran types at some point in the future so they're self-documenting.