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).
Update (Jun 16, 2020): The cctk_nghostzones idea in the previous paragraph (combined with calling the function in LOCAL) mode results in issues when certain numbers of MPI processes are used (thanks to Tanmayee Gupte for pointing this out!) I ended up using Erik’s suggestion (run in LEVEL mode) below. 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.