boolean scale_by_resolution "scale error estimate by geometric mean of grid spacing"

BOOLEAN scale_by_resolution "Scale error estimate by geometric mean of grid spacing"

CCTK_REAL epsilon "Epsilon to avoid infinity when dividing by zero"
{
  0.0:* :: ""
} 1.0e-10

#include <algorithm>

  const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
  const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
  const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
  const CCTK_REAL scalefactor = scale_by_resolution ? cbrt(dx * dy * dz) : 1;

  const CCTK_REAL scalefactor = scale_by_resolution ?
    pow(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(2), 1./3.) :
    1.;

    auto r = sqrt(pow2(p.x) + pow2(p.y) + pow2(p.z));
    regrid_error_(p.I) = scalefactor / fmax(r, CCTK_REAL(1.0e-10));

    CCTK_REAL maxerr = 0;
    // Sample the neighbourhood of this point
    for (int k = -1; k <= +1; ++k) {
      for (int j = -1; j <= +1; ++j) {
        for (int i = -1; i <= +1; ++i) {
          const auto r = sqrt(pow2(p.x + i * dx / 2) + pow2(p.y + j * dy / 2) +
                              pow2(p.z + k * dz / 2));
          const auto err = scalefactor / fmax(r, epsilon);
          maxerr = fmax(maxerr, err);
        }
      }
    }
    regrid_error_(p.I) = maxerr;