Kranc's automatic Boundary code does not work well with Llama

[rhaas80]

Kranc's automated boundary code for non-MoL code looks like this Auxiliary/Cactus/SourceFiles/Kranc.cc:

/*********************************************************************
 * GetBoundaryWidths
 *********************************************************************/

void GetBoundaryWidths(cGH const * restrict const cctkGH, CCTK_INT nboundaryzones[6])
{
  CCTK_INT is_internal[6];
  CCTK_INT is_staggered[6];
  CCTK_INT shiftout[6];
  int ierr = -1;

  if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) {
    int const map = MultiPatch_GetMap (cctkGH);
    /* This doesn't make sense in level mode */
    if (map < 0)
    {
      static int have_warned = 0;
      if (!have_warned)
      {
        CCTK_WARN(1, "GetBoundaryWidths: Could not determine current map (can be caused by calling in LEVEL mode)");
        have_warned = 1;
      }
      for (int i = 0; i < 6; i++)
        nboundaryzones[i] = 0;
      return;
    }
    ierr = MultiPatch_GetBoundarySpecification
      (map, 6, nboundaryzones, is_internal, is_staggered, shiftout);
    if (ierr != 0)
      CCTK_WARN(0, "Could not obtain boundary specification");
  } else if (CCTK_IsFunctionAliased ("GetBoundarySpecification")) {
    ierr = GetBoundarySpecification
      (6, nboundaryzones, is_internal, is_staggered, shiftout);
    if (ierr != 0)
      CCTK_WARN(0, "Could not obtain boundary specification");
  } else {
    CCTK_WARN(0, "Could not obtain boundary specification");
  }
}

/*********************************************************************
 * GetBoundaryWidth
 *********************************************************************/

int GetBoundaryWidth(cGH const * restrict const cctkGH)
{
  CCTK_INT nboundaryzones[6];
  GetBoundaryWidths(cctkGH, nboundaryzones);

  int bw = nboundaryzones[0];

  for (int i = 1; i < 6; i++)
    if (nboundaryzones[i] != bw)
    CCTK_WARN(0, "Number of boundary points is different on different faces");

  return bw;
}

which is called from routines scheduled in LEVEL mode eg Tools/CodeGen/CalculationBoundaries.m:

CalculationBoundariesFunction[calc_List] :=
  Module[{funcName, selectGroup, groups, groupNamesSet, imp},

  funcName = lookup[calc, Name];
  imp = lookup[calc,Implementation];

  (* If the calculation sets every grid point, we don't need to apply
     any boundary condition *)
  If[lookupDefault[calc, Where, Everywhere] === Everywhere,
     Return[{}]];

    (* The boundary interface allows you to give different boundary
       widths for each boundary face.  However, this is not compatible
       with ReflectionSymmetry, so we don't allow it here.*)
    selectGroup[g_String] :=
      Module[{},
        {(* "table = GenericFD_BoundaryWidthTable(cctkGH);\n", *)
         "ierr = KrancBdy_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GetBoundaryWidth(cctkGH), -1 /* no table */, ",
                 Quote[g], ",", Quote["flat"], ");\n",
        "if (ierr < 0)\n",
        "  CCTK_WARN(CCTK_WARN_ALERT, " <> Quote["Failed to register flat BC for "<>g<>"."] <> ");\n" }];

    groups = lookup[calc, Groups];
    groupNamesSet = qualifyGroupName[#, imp] & /@ groupsSetInCalc[calc, groups];

    DefineCCTKFunction[funcName <> "_SelectBCs", "void",
    {
      "if (cctk_iteration % " <> funcName <> "_calc_every != " <> funcName <> "_calc_offset)\n",
      "  return;\n",
      DefineVariable["ierr",   "CCTK_INT", "0"],
(*      DefineVariable["table",  "CCTK_INT", "-1"],*)
      Map[selectGroup, groupNamesSet],
      "return;\n"
    }]];

and scheduled in LEVEL mode (see Tools/CodeGen/Schedule.m):

      selbcSched = {
        Name               -> lookup[calc, Name] <> "_SelectBCs",
        SchedulePoint      -> "in " <> bcGroupName,
        SynchronizedGroups -> groupsToSync,
        Options               -> "level",
        Language           -> CodeGenC`SOURCELANGUAGE,
        Comment            -> lookup[calc, Name] <> "_SelectBCs"
      };

Since is called in LEVEL mode map is -1 and all boundary sizes are set to 0. This is indeed what happens in ML_ADMConstraints when used with a Llama run (eg the multipatch test arrangements/McLachlan/ML_BSSN_Test/test/ML_BSSN_MP_O8_bh.par).

On the MoL side things are somehow even worse. The code in will always ask for 1 boundary point (Tools/CodeGen/MoL.m):

    {"\n",
     "if (CCTK_EQUALS(" <> boundpar <> ", \"none\"  ) ||\n",
     "    CCTK_EQUALS(" <> boundpar <> ", \"static\") ||\n",
     "    CCTK_EQUALS(" <> boundpar <> ", \"flat\"  ) ||\n",
     "    CCTK_EQUALS(" <> boundpar <> ", \"zero\"  ) )\n",
     "{\n",

     "  ierr = KrancBdy_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,\n",
     "                    \"" <> fullgroupname <> "\", " <> boundpar <> ");\n",

     "  if (ierr < 0)\n",
     "     CCTK_ERROR(\"Failed to register "<>boundpar<>" BC for "<>fullgroupname<>"!\");\n",

     "}\n"}];

which really only works since the physical boundary condition is usually NewRad which is not called via Boundary but as part of the RHS and b/c the symmetry thorns (which are called via Boundary) use cctk_nghostzones rather than the boundary width (which is ok since they, or at least the Reflection thorn, check that the number of ghost zones agrees with the boundary width). This however would fail is one used eg a scalar boundary condition with an outer boundary width of more than 1 (or flat).

For the non-MoL code it seems to me as if one should collect all boundary widths from all maps and use that, eg code somewhat like this one (which is utterly untested):

  if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) {
    int const maps = MultiPatch_GetMaps (cctkGH);
    for (int i = 0; i < 6; i++)
      nboundaryzones[i] = -1;
    // collect all boundary information from all maps, and get 
    // outer boundary information only
    for (int m = 0; m < maps; m++) {
      CCTK_INT nbounds[6];
      ierr = MultiPatch_GetBoundarySpecification
        (m, 6, nboundaryzones, is_internal, is_staggered, shiftout);
      if (ierr != 0)
        CCTK_WARN(0, "Could not obtain boundary specification");
      for(int i = 0; i < 6; i++) {
        if(!is_internal[i]) {
          if(nboundaryzones[i] == -1) {
            nboundaryzones[i] = nbounds[i];
          } else {
            if(nboundaryzones[i] != nbounds[i]) {
              CCTK_WARN(0,
                        "Inconsistent number of outer boundary points: %d != %d in the %s %c direction between map %d and an earlier map",
                        nboundaryzones[i], nbounds[i], i % 2 ? "upper" : "lower", "xyz" (i/2), m);
            } else {
              // identical boundary sizes, all is fine
            }
          }
        }
      }
    }
  }