Method test_level_grad_phi_prev

Source/gravity/Gravity.cpp:1071–1128

}

void
Gravity::test_level_grad_phi_prev(int level)
{
    BL_PROFILE("Gravity::test_level_grad_phi_prev()");

    // Fill the RHS for the solve
    const MultiFab& S_old = LevelData[level]->get_old_data(State_Type);
    MultiFab Rhs(grids[level],dmap[level],1,0);
    MultiFab::Copy(Rhs,S_old, URHO,0,1,0);

    const Geometry& geom = parent->Geom(level);

    // This is a correction for fully periodic domains only
    if ( geom.isAllPeriodic() )
    {
       if (gravity::verbose > 1 && ParallelDescriptor::IOProcessor() && mass_offset != 0.0) {
          std::cout << " ... subtracting average density from RHS at level ... "
                    << level << " " << mass_offset << std::endl;
       }
       Rhs.plus(-mass_offset,0,1,0);
    }

    Rhs.mult(Ggravity);

    if (gravity::verbose > 1) {
       Real rhsnorm = Rhs.norm0();
       amrex::Print() << "... test_level_grad_phi_prev at level " << level << std::endl;
       amrex::Print() << "       norm of RHS             " << rhsnorm << std::endl;
    }

    auto dx     = parent->Geom(level).CellSizeArray();
    auto problo = parent->Geom(level).ProbLoArray();
    const int coord_type = geom.Coord();

#ifdef _OPENMP
#pragma omp parallel
#endif
    for (MFIter mfi(Rhs, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();

        test_residual(bx,
                      Rhs.array(mfi),
                      (*grad_phi_prev[level][0]).array(mfi),
#if AMREX_SPACEDIM >= 2
                      (*grad_phi_prev[level][1]).array(mfi),
#endif
#if AMREX_SPACEDIM == 3
                      (*grad_phi_prev[level][2]).array(mfi),
#endif
                      dx, problo, coord_type);
    }

    if (gravity::verbose > 1) {
       Real resnorm = Rhs.norm0();
       amrex::Print() << "       norm of residual        " << resnorm << std::endl;
    }
}