Meudon_BinNS fails if initial_dtshift or initial_dtlapse is set to an unexpected value

[Roland Haas]

Meudon_BinNS contains code like this:

‌

    if (CCTK_EQUALS(initial_lapse, "Meudon_Bin_NS")) { 
      if (CCTK_EQUALS (initial_dtlapse, "Meudon_Bin_NS")) {
        CCTK_INFO ("Calculating time derivatives of lapse");
        set_dt_from_domega (CCTK_PASS_CTOC, alp, dtalp, omega);
      } else if (CCTK_EQUALS (initial_dtlapse, "none") or CCTK_EQUALS(initial_dtlapse,"zero")) {
        // do nothing
      } else {
        CCTK_WARN (CCTK_WARN_ABORT, "internal error");
      }
    }

which fails if initial_lapse is set to Meudon_Bin_NS but initial_dtlapse is set to an unexpected value (and similar for shift). It also only sets dtalp if it sets alp.

Neither is correct. The correct way, documented here https://www.einsteintoolkit.org/thornguide/EinsteinBase/ADMBase/documentation.html#x1-80003.2, is to inspect each parameter independently and then, if the the parameter has a value matching one of the values that Meudon_Bin_NS added to the list of allowed keyword values, set the variable accordingly.

That is something along the lines:

  for (int i=0; i<npoints; ++i) {
    if (CCTK_EQUALS(initial_lapse, "Meudon_Bin_NS")) { 
      alp[i] = bin_ns.nnn[i];
    }
  }
  if (CCTK_EQUALS (initial_dtlapse, "Meudon_Bin_NS")) {
    // use bin_ns.nnn to compute derivs of lapse in case alp is not set by us
    set_dt_from_domega (CCTK_PASS_CTOC, bin_ns.nnn, dtalp, omega);
  }

‌

[Roland Haas]

One may also have to check if this will not fail when cctk_lsh != cctk_ash ie when padding is used. Currently the code uses:

  int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
  vector<double> xx(npoints), yy(npoints), zz(npoints);

and

  for (int i=0; i<npoints; ++i) {
    if (CCTK_EQUALS(initial_lapse, "Meudon_Bin_NS")) { 
      alp[i] = bin_ns.nnn[i];
    }
  }

which is incorrect when padding is used.

‌

[Roland Haas]

There’s also a bug in that setting dtalp requires computing spatial derivatives (via SummationByParts Diff_Gv aliased function) so that it can only be computed in the interior, but no SYNC or boundary condition is applied. The easiest fix to this (since boundary conditions applied to initial data are tricky) would likely be to compute alp on a grid larger than the Cactus grid so that dx_alp can be computed on the whole Cactus grid.