Assertion raised deep inside MLMG solver when initializing MLPoisson

[BenWibking]

I have adapted the gravity code from Castro to work with my AMRCore-based code (quokka-astro/quokka#20) but I get a runtime assertion when calling the MLPoisson constructor

MLPoisson mlpoisson(gmv, bav, dmv, info);

that I can't figure out:

AMReX (22.01-1-ga73ce6759c91) initialized
Initializing gravity solver...
Initialized gravity solver.
Starting solve...
Gravity<T>::fill_multipole_BCs() time = 0.137941189

terminate called after throwing an instance of 'std::runtime_error'
  what():  Assertion `rhs.ixType().cellCentered()' failed, file "/home/benwibking/quokka-code/extern/amrex/Src/Base/AMReX_BoxArray.cpp", line 324

There are no memory errors or warnings reported by valgrind. I get the following backtrace:

(gdb) bt
#0  0x00007ffff7aab18b in raise () from /usr/lib/x86_64-linux-gnu/libc.so.6
#1  0x00007ffff7a8a859 in abort () from /usr/lib/x86_64-linux-gnu/libc.so.6
#2  0x00007ffff7e84911 in ?? () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#3  0x00007ffff7e9038c in ?? () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#4  0x00007ffff7e903f7 in std::terminate() () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#5  0x00007ffff7e906a9 in __cxa_throw () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#6  0x0000555555689722 in amrex::Assert_host (EX=0x5555559da072 "rhs.ixType().cellCentered()", 
    file=0x5555559da008 "/home/benwibking/quokka-code/extern/amrex/Src/Base/AMReX_BoxArray.cpp", line=324, msg=0x0)
    at /home/benwibking/quokka-code/extern/amrex/Src/Base/AMReX.cpp:275
#7  0x000055555574c02e in amrex::Assert (msg=0x0, line=324, 
    file=0x5555559da008 "/home/benwibking/quokka-code/extern/amrex/Src/Base/AMReX_BoxArray.cpp", 
    EX=0x5555559da072 "rhs.ixType().cellCentered()") at /home/benwibking/quokka-code/extern/amrex/Src/Base/AMReX.H:182
#8  amrex::BoxArray::BoxArray (this=0x7fffffff8cc0, rhs=..., trans=...)
    at /home/benwibking/quokka-code/extern/amrex/Src/Base/AMReX_BoxArray.cpp:324
#9  0x000055555596117d in amrex::BndryRegister::define (this=0x555555b98f00, _face=..., _typ=..., _in_rad=1, _out_rad=0, 
    _extent_rad=0, _ncomp=1, dmap=...)
    at /home/benwibking/quokka-code/extern/amrex/Src/Boundary/AMReX_BndryRegister.cpp:117
#10 0x0000555555960605 in amrex::BndryRegister::define (this=0x555555b98f00, grids_=..., dmap=..., in_rad=1, out_rad=0, 
    extent_rad=0, ncomp=1) at /home/benwibking/quokka-code/extern/amrex/Src/Boundary/AMReX_BndryRegister.cpp:46
#11 0x0000555555913190 in amrex::MLCellLinOp::defineAuxData (this=0x7fffffff93c0)
    at /home/benwibking/quokka-code/extern/amrex/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp:108
#12 0x0000555555912f28 in amrex::MLCellLinOp::define (this=0x7fffffff93c0, a_geom=..., a_grids=..., a_dmap=..., 
    a_info=..., a_factory=...)
    at /home/benwibking/quokka-code/extern/amrex/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp:88
#13 0x00005555559997f1 in amrex::MLCellABecLap::define (this=0x7fffffff93c0, a_geom=..., a_grids=..., a_dmap=..., 
    a_info=..., a_factory=...)
    at /home/benwibking/quokka-code/extern/amrex/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp:26
#14 0x000055555594883b in amrex::MLPoisson::define (this=0x7fffffff93c0, a_geom=..., a_grids=..., a_dmap=..., 
--Type <RET> for more, q to quit, c to continue without paging--
    a_info=..., a_factory=...) at /home/benwibking/quokka-code/extern/amrex/Src/LinearSolvers/MLMG/AMReX_MLPoisson.cpp:35
#15 0x00005555559486cb in amrex::MLPoisson::MLPoisson (this=0x7fffffff93c0, a_geom=..., a_grids=..., a_dmap=..., 
    a_info=..., a_factory=...) at /home/benwibking/quokka-code/extern/amrex/Src/LinearSolvers/MLMG/AMReX_MLPoisson.cpp:14
#16 0x000055555558ea22 in Gravity<PoissonProblem>::actual_solve_with_mlmg (this=0x7fffffff9bf0, crse_level=0, 
    fine_level=0, phi=..., rhs=..., grad_phi=..., res=..., crse_bcdata=0x0, rel_eps=0, abs_eps=0)
    at /home/benwibking/quokka-code/src/Gravity_impl.hpp:1031
#17 0x0000555555582b37 in Gravity<PoissonProblem>::solve_phi_with_mlmg (this=0x7fffffff9bf0, crse_level=0, fine_level=0, 
    phi=..., rhs=..., grad_phi=..., res=..., time=0) at /home/benwibking/quokka-code/src/Gravity_impl.hpp:967
#18 0x000055555557970e in Gravity<PoissonProblem>::solve_for_phi (this=0x7fffffff9bf0, level=0, phi=..., grad_phi=..., 
    is_new=1) at /home/benwibking/quokka-code/src/Gravity_impl.hpp:321
#19 0x0000555555567f5c in problem_main () at /home/benwibking/quokka-code/src/PoissonFreeSpace/test_poisson.cpp:119
#20 0x0000555555565512 in main (argc=2, argv=0x7fffffffa538) at /home/benwibking/quokka-code/src/main.cpp:41

Does this point to something relatively obvious? I'm afraid I don't understand what is going on behind the scenes enough to step through each frame of the backtrace to see where it starts to go wrong.

[WeiqunZhang]

It looks like the BoxArray passed to MLPoisson is not cell-centered. MLPoisson is for cell-centered only.

[BenWibking]

Yes, they have the same name by accident. But this rhs is passed to MLPoisson via the rhs argument of solve_phi_with_mlmg. Walking back through the call stack, you can see that the BoxArray passed to the constructor of MLPoisson is actually the same BoxArray as the one above.

[WeiqunZhang]

Then that's very strange. You can add more assertions on the call stack to see what happens.

[BenWibking]

Ok, I think I've figured it out. I had inadvertently called sim->grids[lev].surroundingNodes(dir), which modified sim->grids[lev] in place, rather that returning a new BoxArray. However, when I change it to call sim->boxArray(lev).surroundingNodes(dir) instead, I get a compiler error:

[build] ../src/Gravity_impl.hpp: In instantiation of ‘void Gravity::install_level(int) [with T = PoissonProblem]’:
[build] ../src/PoissonFreeSpace/test_poisson.cpp:114:33: required from here
[build] ../src/Gravity_impl.hpp:219:12: error: passing ‘const amrex::BoxArray’ as ‘this’ argument discards qualifiers [-fpermissive]
[build] 219 | auto ba = sim->boxArray(level).surroundingNodes(n);
[build] | ^~
[build] In file included from ../extern/amrex/Src/Base/AMReX_Utility.H:11,
[build] from ../extern/amrex/Src/Base/AMReX_GpuError.H:5,
[build] from ../extern/amrex/Src/Base/AMReX_GpuLaunch.H:9,
[build] from ../extern/amrex/Src/Base/AMReX_PODVector.H:6,
[build] from ../extern/amrex/Src/Base/AMReX_GpuContainers.H:6,
[build] from ../src/Gravity.H:16,
[build] from ../src/PoissonFreeSpace/test_poisson.cpp:17:
[build] ../extern/amrex/Src/Base/AMReX_BoxArray.H:686:15: note: in call to ‘amrex::BoxArray& amrex::BoxArray::surroundingNodes(int)’
[build] 686 | BoxArray& surroundingNodes (int dir);
[build] | ^~~~~~~~~~~~~~~~
[build] ninja: build stopped: subcommand failed.

However, calling this elsewhere in the code works fine. And if I call it on a separate line, it also works:

    const DistributionMapping &dm = sim->DistributionMap(level);

    grad_phi_prev[level].resize(AMREX_SPACEDIM);
    for (int n = 0; n < AMREX_SPACEDIM; ++n) {
      auto ba = sim->boxArray(level).surroundingNodes(n);

      grad_phi_prev[level][n] = std::make_unique<MultiFab>(ba, dm, 1, 1);
    }

    grad_phi_curr[level].resize(AMREX_SPACEDIM);
    for (int n = 0; n < AMREX_SPACEDIM; ++n) {
      auto ba = sim->boxArray(level);
      ba.surroundingNodes(n);

      grad_phi_curr[level][n] = std::make_unique<MultiFab>(ba, dm, 1, 1);
    }

The first loop leads to a compiler error, whereas the code compiles and runs correctly if I do it the second way. Does this make any sense?

[WeiqunZhang]

Here is the member function, boxArray.

const BoxArray& boxArray (int lev) const noexcept { return grids[lev]; }

It returns const &. So boxArray(level).surroundingNodes(n) does not compile because surroundingNodes is not a const funciton. The second approach works because auto ba = .. is BoxArray ba = .... So it makes a copy.

There is another approach,

... = std::make_unique<MultiFab>(amrex::convert(sim->boxArray(lev), IntVect::TheDimensionVector(n)), dm, 1, 1);

[BenWibking]

Ah, right. I've switched to the new approach. Thanks!