`CheckRhoSieFromPT` is a bit weird

[bgclayton-lanl]

The subroutine in question is:

template <typename EOS, typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION bool
CheckRhoSieFromPT(EOS eos, Real rho, Real T,
                  Indexer_t &&lambda = static_cast<Real *>(nullptr)) {
  const Real P = eos.PressureFromDensityTemperature(rho, T, lambda);
  const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T, lambda);
  Real rtest = 12; // set these to something
  Real etest = 1;
  eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest);
  Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda);
  Real residual = P_test - P;
  Real frac_residual =
      std::min(std::abs(singularity::robust::ratio(residual, P)), std::abs(residual));
  bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8))
                      // This is as good as it will get sometimes.
                      || (std::abs(frac_residual) <= 1e-8);
  if (!results_good) {
    printf("RhoSie of PT failure!\n"
           "\trho_true = %.14e\n"
           "\tsie_true = %.14e\n"
           "\tP        = %.14e\n"
           "\tT        = %.14e\n"
           "\trho      = %.14e\n"
           "\tsie      = %.14e\n"
           "\tP_test   = %.14e\n"
           "\tresidual = %.14e\n"
           "\tfracres  = %.14e\n",
           rho, sie, P, T, rtest, etest, P_test, residual, frac_residual);
  }
  return results_good;
}

As I understand it, is seems like this subroutine want to do the following. Initialize a state $(P_0, T_0)$, then to compute $\rho_0 := \rho(P_0,T_0)$ and $e_0 := e(P_0, T_0)$, and then check that the inversion is consistent $P(\rho_0, e_0) = P_0$ and $T(\rho_0, e_0) = T_0$. Basically, just check that inversion is consistent and there are no mistakes in the formula. I guess you could also start with $(\rho_0, e_0)$ and then do the reverse process.

I'm assuming the reason this might be happening this way is because certain EOS member function calls aren't available? What do we think?

[Yurlungur]

Generically, most equations of state are expessed with density and either energy or temperature as independent variables. Currently none of our EOS's depend on pressure as an independent variable. So $(r\rho, e)(P, T)$ is an inversion, possibly twice. As you say, this method is to test that that inversion works properly. You should be able to go forward or backwards and get the same state (within Newton solver tolerances).

If an equation of state had pressure as an independent variable, I would suggest an inverted process, so that the thing that is checked is always the inversion. Does that make sense?