Evolution of grid arrays with Mol no longer works

[Former user]

In the releases prior to ET_2020_05 I was able to evolve grid arrays using RK4 with a kludge algorithm coupled to MoL. (The issue with MoL and grid arrays is that Mol attempts to updates the arrays for all components of the grid owned by that process. The kludge was to check to see that if the "correct" component is being updated, and if not, re-write the old values of the grid array). I used this to evolve geodesics and the puncture locations. From ET_2020_05, I get the error

terminate called after throwing an instance of 'std::out_of_range' what(): vector::_M_range_check

1. CarpetLib::signal_handler(int) [./cactus_test(_ZN9CarpetLib14signal_handlerEi+0xf4) [0x6b1c74]]
2. /lib64/libc.so.6(+0x35270) [0x7f4fbc077270]
3. /lib64/libc.so.6(gsignal+0x37) [0x7f4fbc0771f7]
4. /lib64/libc.so.6(abort+0x148) [0x7f4fbc0788e8]
5. __gnu_cxx::__verbose_terminate_handler() [/lib64/libstdc++.so.6(_ZN9__gnu_cxx27__verbose_terminate_handlerEv+0x165) [0x7f4fbca9bac5]]
6. /lib64/libstdc++.so.6(+0x5ea36) [0x7f4fbca99a36]
7. /lib64/libstdc++.so.6(+0x5ea63) [0x7f4fbca99a63]
8. /lib64/libstdc++.so.6(+0x5ec83) [0x7f4fbca99c83]
9. std::__throw_out_of_range(char const*) [/lib64/libstdc++.so.6(_ZSt20__throw_out_of_rangePKc+0x77) [0x7f4fbcaeea97]]
10. Carpet::PoisonGroup(_cGH const*, int, Carpet::checktimes) [./cactus_test(_ZN6Carpet11PoisonGroupEPK4_cGHiNS_10checktimesE+0x551) [0x5241f1]]
11. Carpet::Poison(_cGH const*, Carpet::checktimes, int) [./cactus_test(_ZN6Carpet6PoisonEPK4_cGHNS_10checktimesEi+0xcd) [0x523bbd]]
12. ./cactus_test() [0x51881d]
13. Carpet::Evolve(tFleshConfig*) [./cactus_test(_ZN6Carpet6EvolveEP12tFleshConfig+0x191) [0x5162e1]]
14. ./cactus_test(main+0xa5) [0xd560e5]
15. /lib64/libc.so.6(__libc_start_main+0xf5) [0x7f4fbc063c05]
16. ./cactus_test() [0x4b6ce9]

[Yosef Zlochower]

  <div class="preview-container wiki-content"><!-- loaded via ajax --></div>
  <div class="mask"></div>
</div>

</div> </form>

[Yosef Zlochower]

  <div class="preview-container wiki-content"><!-- loaded via ajax --></div>
  <div class="mask"></div>
</div>

</div> </form>

[Yosef Zlochower]

The attached thorn simply registers and array with for evolution with MoL. The thorn is not intended to actually evolve a real array. The issue also only shows up if the number of refinement levels is larger than 1

‌

[Yosef Zlochower]

The following patch seems to “fix” the problem.

diff --git a/Carpet/src/Poison.cc b/Carpet/src/Poison.cc
index 12fc3d3..39a989a 100644
--- a/Carpet/src/Poison.cc
+++ b/Carpet/src/Poison.cc
@@ -131,12 +131,16 @@ void PoisonGroup(cGH const *const cctkGH, int const group,
     END_LOCAL_MAP_LOOP;

     // mark as invalid
-    for (int var = 0; var < nvar; ++var) {
-      int const map0 = 0;
-      ggf *const ff = arrdata.AT(group).AT(map0).data.AT(var);
-      assert(ff);
-      for (int tl = min_tl; tl <= max_tl; ++tl) {
-        ff->set_valid(mglevel, reflevel, tl, CCTK_VALID_NOWHERE);
+    // Exclude CCTK_ARRAY from this
+    if (CCTK_GroupTypeI(group) != CCTK_ARRAY)
+    {
+      for (int var = 0; var < nvar; ++var) {
+        int const map0 = 0;
+        ggf *const ff = arrdata.AT(group).AT(map0).data.AT(var);
+        assert(ff);
+        for (int tl = min_tl; tl <= max_tl; ++tl) {
+          ff->set_valid(mglevel, reflevel, tl, CCTK_VALID_NOWHERE);
+        }
       }
     }

‌

[Roland Haas]

I cannot reproduce this with the provided thorn and parfile. The parfile seems to be for a different setup. Eg the thorn name differs, the thorn in the parfile has parameters while the thorn in the tarball does not, the parfile contains @FOO@ placeholders, the parfile has no mesh refinement or MoL active but the ticket description mentions both being required.

I have a guess what is going on, but without a reproducer will not be able to make headway fixing this.

[Roland Haas]

Writing my own parfile:

ActiveThorns = "Time MoL"
ActiveThorns = "Coordbase CartGrid3d Boundary StaticConformal"
ActiveThorns = "SymBase InitBase"
ActiveThorns = "IOUtil"
ActiveThorns = "Carpet CarpetLib CarpetRegrid2"

driver::ghost_size = 1
CoordBase::boundary_size_x_lower = 1
CoordBase::boundary_size_y_lower = 1
CoordBase::boundary_size_z_lower = 1
CoordBase::boundary_size_x_upper = 1
CoordBase::boundary_size_y_upper = 1
CoordBase::boundary_size_z_upper = 1

CartGrid3D::type      = "coordbase"
CoordBase::domainsize = "minmax"
CoordBase::xmin = -32.0
CoordBase::ymin = -32.0
CoordBase::zmin = -32.0
CoordBase::xmax =  32.0
CoordBase::ymax =  32.0
CoordBase::zmax =  32.0
CoordBase::dx   =   8.0
CoordBase::dy   =   8.0
CoordBase::dz   =   8.0

cactus::cctk_itlast = 1
Cactus::terminate           = "iteration"

Carpet::domain_from_coordbase = "yes"
Carpet::use_buffer_zones         = "yes"


Carpet::init_fill_timelevels     = "yes"

CarpetLib::poison_new_memory  = "yes"
CarpetLib::poison_value       = 114
Carpet::poison_new_timelevels = "yes"

Carpet::max_refinement_levels    = 2

CarpetRegrid2::regrid_every = 1
CarpetRegrid2::num_centres  = 1
CarpetRegrid2::num_levels_1 = 2
CarpetRegrid2::radius_1[1]  =16.0

time::dtfac = 0.25
MoL::ODE_Method             = "rk4"
MoL::MoL_Intermediate_Steps = 4
MoL::MoL_Num_Scratch_Levels = 1

ActiveThorns = EvolveArray

I can get a error similar to what is being described in the ticket.
My guess is that this branch rhaas/poisonarray (https://bitbucket.org/eschnett/carpet/branch/rhaas/poisonarray) will fix things since it applies the same logic as is used otherwhere in PreSync.cc to Poison.cc when it come to handling non-grid functions.

[Roland Haas]

Somewhat unrelated: instead of the current workaround for arrays of

The kludge was to check to see that if the "correct" component is being updated, and if not, re-write the old values of the grid array).

wouldn’t it be possible ot make this work by scheduling the RHS function as GLOBAL ? That way it would only execute once per iteration and not once per iteration per refinement level. MoL would still update the state vector for the array but b/c it initializes the RHS to zero before the CalcRHS schedule group the fact that the RHS routine does not run would leave the RHS unchanged (ie 0) and thus effectively no update would be done.

[Yosef Zlochower]

To maintain accuracy for geodesics, I needed to evolve the geodesics at the same time as the underlying grid that they happen to be on.

[Roland Haas]

I see, that would require extra logic even if MoL was being ore accommodating about evolving arrays since you will have to query Carpet about the finest refinement level a given location is contained in anyway. GLOBAL routines are called every timestep (same frequency as eg routines will execute on the finest timestep). Anyway, I was just curious.

[Roland Haas]

@Yosef Zlochower confirmed (private email) that the branch rhaas/poisonarray seems to work.

Applied as git hash 9d0eaa4c "Carpet: handle marking grid-arrays invalid correctly" of carpet