Method test_level_grad_phi_curr

Source/gravity/Gravity.cpp:1130–1188

}

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

    // Fill the RHS for the solve
    const MultiFab& S_new = LevelData[level]->get_new_data(State_Type);
    MultiFab Rhs(grids[level],dmap[level],1,0);
    MultiFab::Copy(Rhs,S_new, 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 in solve ... " << mass_offset << std::endl;
       }
       Rhs.plus(-mass_offset,0,1,0);
    }

    Rhs.mult(Ggravity);

    if (gravity::verbose > 1) {
       Real rhsnorm = Rhs.norm0();
       if (ParallelDescriptor::IOProcessor()) {
          std::cout << "... test_level_grad_phi_curr at level " << level << std::endl;
          std::cout << "       norm of RHS             " << rhsnorm << std::endl;
        }
    }

    auto dx     = geom.CellSizeArray();
    auto problo = geom.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_curr[level][0]).array(mfi),
#if AMREX_SPACEDIM >= 2
                      (*grad_phi_curr[level][1]).array(mfi),
#endif
#if AMREX_SPACEDIM == 3
                      (*grad_phi_curr[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;
    }