In developing BaikalETK (a NRPy+-based BSSN thorn), I set up a schedule.ccl file that follows almost exactly the same pattern as Lean.
I noticed that when BaikalETK evaluated the constraints without initializing the variables to zero (as is done in Lean), I would get uninitialized on the boundaries **after ApplyBCs was called**. This confused me since BaikalETK, like Lean, calls for "flat" boundary conditions on the constraints. Sure enough when I performed the same run with Lean, most boundary points were filled with zeros (unlike BaikalETK, the constraints in Lean are initialized to zero). I spoke to Roland (Haas) about this, and he was also confused until he noticed that when BaikalETK/Lean calls Boundary_SelectVarForBC/Boundary_SelectGroupForBC a la:
Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, +1, -1, "BaikalETK::somegridfunction", "flat") ,
the "+1" indicates that only one ghost zone point should be filled in. You might then wonder what would happen if we just replaced that "+1" with "cctk_nghostzones". Well, since this function in Lean & BaikalETK is called in "LEVEL" mode, "cctk_nghostzones" does not exist (actually it compiles and gives a garbage value).
Roland's brilliant solution: call the function in "Local" mode instead. This applies "flat" boundary conditions on the constraints to *all* ghostzones. It's currently implemented in BaikalETK (here's the schedule.ccl: https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-BaikalETK.ipynb#scheduleccl) and seems to work great.