diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index dcbf566..bc5b679 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -33,7 +33,7 @@ jobs: -DPELERAD_DIM=3 \ -DPELERAD_ENABLE_MPI=OFF \ -DPELERAD_ENABLE_OPENMP=OFF \ - -DPELERAD_ENABLE_EB=OFF \ + -DPELERAD_ENABLE_EB=ON \ -DPELERAD_ENABLE_LINEARSOLVERS=ON \ -DPELERAD_ENABLE_HYPRE=OFF \ -DPELERAD_ENABLE_TESTS=ON \ @@ -79,7 +79,7 @@ jobs: -DPELERAD_DIM=3 \ -DPELERAD_ENABLE_MPI=ON \ -DPELERAD_ENABLE_OPENMP=OFF \ - -DPELERAD_ENABLE_EB=OFF \ + -DPELERAD_ENABLE_EB=ON \ -DPELERAD_ENABLE_LINEARSOLVERS=ON \ -DPELERAD_ENABLE_HYPRE=OFF \ -DPELERAD_ENABLE_TESTS=ON \ diff --git a/inputs/inputs.tstPOneSingleEB b/inputs/inputs.tstPOneSingleEB new file mode 100644 index 0000000..f60cce1 --- /dev/null +++ b/inputs/inputs.tstPOneSingleEB @@ -0,0 +1,23 @@ +#amr +max_level = 0 +n_cell = 64 +max_grid_size = 32 +plot_file_name = "pltSingleEB" + +#mlmg +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 10 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.e-5 +abstol = 1.e-5 +bottom_reltol = 1.e-5 +bottom_abstol = 1.e-5 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +lo_bc = Periodic Periodic Periodic +hi_bc = Periodic Periodic Periodic diff --git a/inputs/inputs.tstPOneSolve b/inputs/inputs.tstPOneSolve deleted file mode 100644 index 24c0f4e..0000000 --- a/inputs/inputs.tstPOneSolve +++ /dev/null @@ -1,29 +0,0 @@ -# Problem -prob.bc_type = Dirichlet -#prob.bc_type = Neumann -#prob.bc_type = Periodic - -composite_solve = 1 - -#amr -max_level = 1 -n_cell = 128 -ref_ratio = 2 -max_grid_size = 64 -plot_file_name = "plt" - -#mlmg -use_hypre = 0 -verbose = 2 -bottom_verbose = 2 -fine_level_solve_only = 0 -max_iter = 100 -max_fmg_iter = 20 -max_bottom_iter = 1000 -max_coarsening_level = 30 -reltol = 1.e-6 -abstol = 0.0 -bottom_reltol = 1.e-4 -linop_maxorder = 2 -agglomeration = 1 -consolidation = 1 diff --git a/inputs/inputs.tstPOneSolveEB b/inputs/inputs.tstPOneSolveEB deleted file mode 100644 index 65e4ff0..0000000 --- a/inputs/inputs.tstPOneSolveEB +++ /dev/null @@ -1,26 +0,0 @@ -amrex.fpe_trap_invalid = 1 - -#mlmg -verbose = 2 -bottom_verbose = 2 -max_iter = 100 -max_fmg_iter = 0 -max_bottom_iter = 1000 -reltol = 1.e-12 -bottom_reltol = 1.e-4 -linop_maxorder = 3 -max_coarsening_level = 30 -agg_grid_size = -1 -con_grid_size = -1 - -#amr -prob_type = 1 -max_level = 0 -n_cell = 128 -ref_ratio = 2 -max_grid_size = 64 -plot_file_name = "plt" - -eb2.max_grid_size = 32 -eb2.geom_type = flower - diff --git a/scripts/build_frontier_hip.sh b/scripts/build_frontier_hip.sh index 3740a4a..6918f85 100755 --- a/scripts/build_frontier_hip.sh +++ b/scripts/build_frontier_hip.sh @@ -6,7 +6,7 @@ cmake \ -DCMAKE_CXX_COMPILER=$(which CC)\ -DPELERAD_DIM=3 \ -DPELERAD_ENABLE_MPI=OFF \ - -DPELERAD_ENABLE_EB=OFF \ + -DPELERAD_ENABLE_EB=ON \ -DPELERAD_ENABLE_LINEARSOLVERS=ON \ -DPELERAD_ENABLE_HYPRE=OFF \ -DPELERAD_ENABLE_PARTICLES=OFF \ diff --git a/src/POneSingleEB.hpp b/src/POneSingleEB.hpp new file mode 100644 index 0000000..8f8f184 --- /dev/null +++ b/src/POneSingleEB.hpp @@ -0,0 +1,168 @@ +#ifndef PONESINGLE_HPP +#define PONESINGLE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace PeleRad +{ + +class POneSingleEB +{ +private: + MLMGParam mlmgpp_; + +public: + amrex::Geometry const& geom_; + amrex::BoxArray const& grids_; + amrex::DistributionMapping const& dmap_; + + std::unique_ptr const& factory_; + + amrex::MultiFab& solution_; + amrex::MultiFab const& rhs_; + amrex::MultiFab const& acoef_; + amrex::Array& bcoef_; + + amrex::MultiFab const& robin_a_; + amrex::MultiFab const& robin_b_; + amrex::MultiFab const& robin_f_; + + amrex::Real const ascalar = 1.0; + amrex::Real const bscalar = 1.0 / 3.0; + + // constructor + POneSingleEB(MLMGParam const& mlmgpp, amrex::Geometry const& geom, + amrex::BoxArray const& grids, amrex::DistributionMapping const& dmap, + std::unique_ptr const& factory, + amrex::MultiFab& solution, amrex::MultiFab const& rhs, + amrex::MultiFab const& acoef, + amrex::Array& bcoef, + amrex::MultiFab const& robin_a, amrex::MultiFab const& robin_b, + amrex::MultiFab const& robin_f) + : mlmgpp_(mlmgpp), + geom_(geom), + grids_(grids), + dmap_(dmap), + factory_(factory), + solution_(solution), + rhs_(rhs), + acoef_(acoef), + bcoef_(bcoef), + robin_a_(robin_a), + robin_b_(robin_b), + robin_f_(robin_f) {}; + + void solve() + { + amrex::MultiFab const& vfrc = factory_->getVolFrac(); + + auto const max_coarsening_level = mlmgpp_.max_coarsening_level_; + auto const max_iter = mlmgpp_.max_iter_; + auto const max_fmg_iter = mlmgpp_.max_fmg_iter_; + auto const verbose = mlmgpp_.verbose_; + auto const bottom_verbose = mlmgpp_.bottom_verbose_; + auto const tol_rel = mlmgpp_.reltol_; + auto const tol_abs = mlmgpp_.abstol_; + auto const use_hypre = mlmgpp_.use_hypre_; + + auto const& lobc = mlmgpp_.lobc_; + auto const& hibc = mlmgpp_.hibc_; + + auto const& geom = geom_; + auto const& grids = grids_; + auto const& dmap = dmap_; + auto const& factory = factory_; + + auto& solution = solution_; + auto const& rhs = rhs_; + auto const& acoef = acoef_; + auto& bcoef = bcoef_; + auto const& robin_a = robin_a_; + auto const& robin_b = robin_b_; + auto const& robin_f = robin_f_; + + amrex::MLEBABecLap mlabec({ geom }, { grids }, { dmap }, + amrex::LPInfo().setMaxCoarseningLevel(max_coarsening_level), + { factory.get() }); + + mlabec.setDomainBC(lobc, hibc); + + // mlabec.setLevelBC(0, &solution, &robin_a, &robin_b, &robin_f); + mlabec.setLevelBC(0, nullptr); + + mlabec.setScalars(ascalar, bscalar); + + mlabec.setACoeffs(0, acoef); + + // amrex::Array face_bcoef; + /* + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + bcoef[idim].define( + amrex::convert(grids, + amrex::IntVect::TheDimensionVector(idim)), dmap, 1, 0, + amrex::MFInfo(), *factory); + + const amrex::BoxArray& ba = amrex::convert( + bcoef[idim].boxArray(), + amrex::IntVect::TheDimensionVector(idim)); + face_bcoef[idim].define(ba, + bcoef.DistributionMap(), 1, 0); + + } + */ + // amrex::average_cellcenter_to_face( + // GetArrOfPtrs(face_bcoef), bcoef, geom); + + // mlabec.setBCoeffs(0, amrex::GetArrOfConstPtrs(face_bcoef)); + mlabec.setBCoeffs(0, amrex::GetArrOfConstPtrs(bcoef)); + + // this is the BCoef associated iwth an EB face + amrex::MultiFab beta(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + beta.setVal(1.0); + + mlabec.setEBDirichlet(0, solution, beta); + + amrex::MLMG mlmg(mlabec); + mlmg.setMaxIter(max_iter); + mlmg.setMaxFmgIter(max_fmg_iter); + mlmg.setVerbose(verbose); + mlmg.setBottomVerbose(bottom_verbose); + + if (use_hypre) mlmg.setBottomSolver(amrex::MLMG::BottomSolver::hypre); + + mlmg.solve({ &solution }, { &rhs }, tol_rel, tol_abs); + } + + void calcRadSource(amrex::MultiFab& rad_src) + { + for (amrex::MFIter mfi(rad_src); mfi.isValid(); ++mfi) + { + amrex::Box const& bx = mfi.validbox(); + + auto const& rhsfab = rhs_.array(mfi); + auto const& solfab = solution_.array(mfi); + auto const& acfab = acoef_.array(mfi); + + auto radfab = rad_src.array(mfi); + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + radfab(i, j, k, 4) + = acfab(i, j, k) * solfab(i, j, k) - rhsfab(i, j, k); + }); + } + } +}; + +} + +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index beb306a..f2c073f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,6 +31,14 @@ target_link_libraries(PeleRad_POneSingle.exe PRIVATE PeleRad Boost::unit_test_fr target_compile_definitions(PeleRad_POneSingle.exe PRIVATE BOOST_TEST_DYN_LINK) add_test(NAME PeleRad_POneSingle_Test COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ./PeleRad_POneSingle.exe "${CMAKE_SOURCE_DIR}/inputs/inputs.tstPOneSingle") +if(PELERAD_ENABLE_EB) +add_executable(PeleRad_POneSingleEB.exe tstPOneSingleEB.cpp ut_main.cpp) +target_include_directories(PeleRad_POneSingleEB.exe PUBLIC ${AMREX_HOME_DIR}/Src/Base ${AMREX_HOME_DIR}/Src/LinearSolvers/MLMG ${AMREX_HOME_DIR}/Src/EB) +target_link_libraries(PeleRad_POneSingleEB.exe PRIVATE PeleRad Boost::unit_test_framework amrex) +target_compile_definitions(PeleRad_POneSingleEB.exe PRIVATE BOOST_TEST_DYN_LINK) +add_test(NAME PeleRad_POneSingleEB_Test COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ./PeleRad_POneSingleEB.exe "${CMAKE_SOURCE_DIR}/inputs/inputs.tstPOneSingleEB") +endif() + add_executable(PeleRad_POneMulti.exe tstPOneMulti.cpp ut_main.cpp) target_include_directories(PeleRad_POneMulti.exe PUBLIC ${AMREX_HOME_DIR}/Src/Base ${AMREX_HOME_DIR}/Src/LinearSolvers/MLMG) target_link_libraries(PeleRad_POneMulti.exe PRIVATE PeleRad Boost::unit_test_framework amrex) @@ -64,6 +72,14 @@ target_link_libraries(PeleRad_POneSingle.exe PRIVATE PeleRad Boost::unit_test_fr target_compile_definitions(PeleRad_POneSingle.exe PRIVATE BOOST_TEST_DYN_LINK) add_test(NAME PeleRad_POneSingle_Test COMMAND ./PeleRad_POneSingle.exe "${CMAKE_SOURCE_DIR}/inputs/inputs.tstPOneSingle" ) +if(PELERAD_ENABLE_EB) +add_executable(PeleRad_POneSingleEB.exe tstPOneSingleEB.cpp ut_main.cpp) +target_include_directories(PeleRad_POneSingleEB.exe PUBLIC ${AMREX_HOME_DIR}/Src/Base ${AMREX_HOME_DIR}/Src/LinearSolvers/MLMG ${AMREX_HOME_DIR}/Src/EB) +target_link_libraries(PeleRad_POneSingleEB.exe PRIVATE PeleRad Boost::unit_test_framework amrex) +target_compile_definitions(PeleRad_POneSingleEB.exe PRIVATE BOOST_TEST_DYN_LINK) +add_test(NAME PeleRad_POneSingleEB_Test COMMAND ./PeleRad_POneSingleEB.exe "${CMAKE_SOURCE_DIR}/inputs/inputs.tstPOneSingleEB" ) +endif() + add_executable(PeleRad_POneMulti.exe tstPOneMulti.cpp ut_main.cpp) target_include_directories(PeleRad_POneMulti.exe PUBLIC ${AMREX_HOME_DIR}/Src/Base ${AMREX_HOME_DIR}/Src/LinearSolvers/MLMG) target_link_libraries(PeleRad_POneMulti.exe PRIVATE PeleRad Boost::unit_test_framework amrex) @@ -90,6 +106,7 @@ setup_target_for_cuda_compilation(PeleRad_AmrexParallelFor.exe) setup_target_for_cuda_compilation(PeleRad_BoostReadDatabase.exe) setup_target_for_cuda_compilation(PeleRad_AmrexGetRadProp.exe) setup_target_for_cuda_compilation(PeleRad_POneSingle.exe) +setup_target_for_cuda_compilation(PeleRad_POneSingleEB.exe) setup_target_for_cuda_compilation(PeleRad_POneMulti.exe) setup_target_for_cuda_compilation(PeleRad_POneSingleAF.exe) setup_target_for_cuda_compilation(PeleRad_POneMultiAF.exe) diff --git a/tests/tstPOneSingleEB.cpp b/tests/tstPOneSingleEB.cpp new file mode 100644 index 0000000..10df21f --- /dev/null +++ b/tests/tstPOneSingleEB.cpp @@ -0,0 +1,279 @@ +#include + +#define BOOST_TEST_MODULE ponerobinsingle + +#include +#include +#include +#include + +#include +#include +#include +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void actual_init_coefs_eb(int i, int j, + int k, amrex::Array4 const& phi, + amrex::Array4 const& phi_exact, + amrex::Array4 const& rhs, + amrex::Array4 const& acoef, + amrex::Array4 const& bx, amrex::Array4 const& by, + amrex::Array4 const& bz, + amrex::Array4 const& flag, + amrex::GpuArray const& dx, + amrex::Box const& vbx) +{ + amrex::Real const L = std::sqrt(2.0); + amrex::Real const n = 3.0; + amrex::Real npioverL = n * M_PI / L; + + amrex::Real pioverfour = M_PI / 4.0; + amrex::Real cospioverfour = std::cos(pioverfour); + amrex::Real sinpioverfour = std::sin(pioverfour); + + if (vbx.contains(i, j, k)) + { + bx(i, j, k) = 1.0; + by(i, j, k) = 1.0; + bz(i, j, k) = 1.0; + + phi(i, j, k) = 0.0; + + if (flag(i, j, k).isCovered()) + { + rhs(i, j, k) = 0.0; + phi_exact(i, j, k) = 0.0; + } + else + { + amrex::Real x = dx[0] * (i + 0.5) - 0.5; + amrex::Real y = dx[1] * (j + 0.5) - 0.5; + amrex::Real z = dx[2] * (k + 0.5) - 0.5; + + amrex::Real xp = x * cospioverfour + y * sinpioverfour; + amrex::Real yp = -x * sinpioverfour + y * cospioverfour; + + double sincossin = std::sin(npioverL * xp) * std::cos(npioverL * yp) + * std::sin(npioverL * z); + + rhs(i, j, k) + = (1.0 + npioverL * npioverL) * bx(i, j, k) * sincossin; + acoef(i, j, k) = 1.0; + + phi_exact(i, j, k) = sincossin; + } + } +} + +void initProbABecLaplacian(amrex::Geometry& geom, + std::unique_ptr& factory, + amrex::MultiFab& solution, amrex::MultiFab& rhs, + amrex::MultiFab& exact_solution, amrex::MultiFab& acoef, + amrex::Array& bcoef) +{ + + amrex::FabArray const& flags + = factory->getMultiEBCellFlagFab(); + // amrex::MultiCutFab const& bcent = factory->getBndryCent(); + // amrex::MultiCutFab const& cent = factory->getCentroid(); + /* + auto const prob_lo = geom.ProbLoArray(); + auto const prob_hi = geom.ProbHiArray(); + auto const dlo = amrex::lbound(geom.Domain()); + auto const dhi = amrex::ubound(geom.Domain()); + */ + auto const dx = geom.CellSizeArray(); + amrex::Box const& domainbox = geom.Domain(); + + for (amrex::MFIter mfi(rhs); mfi.isValid(); ++mfi) + { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& nbx = amrex::surroundingNodes(bx); + // amrex::Box const& gbx = amrex::grow(bx, 1); + amrex::Array4 const& phi_arr = solution.array(mfi); + amrex::Array4 const& phi_ex_arr + = exact_solution.array(mfi); + amrex::Array4 const& rhs_arr = rhs.array(mfi); + + amrex::Array4 const& bx_arr = bcoef[0].array(mfi); + amrex::Array4 const& by_arr = bcoef[1].array(mfi); + amrex::Array4 const& bz_arr = bcoef[2].array(mfi); + + auto fabtyp = flags[mfi].getType(bx); + if (fabtyp == amrex::FabType::covered) + { + std::cout << " amrex::FabType::covered == fabtyp \n"; + } + else if (fabtyp == amrex::FabType::regular) + { + std::cout << " amrex::FabType::regular == fabtyp \n"; + } + else + { + // std::cout << " amrex::FabType else \n"; + amrex::Array4 const& acoef_arr = acoef.array(mfi); + amrex::Array4 const& flag_arr + = flags.const_array(mfi); + // amrex::Array4 const& cent_arr + // = cent.const_array(mfi); + // amrex::Array4 const& bcent_arr + // = bcent.const_array(mfi); + + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + actual_init_coefs_eb(i, j, k, phi_arr, phi_ex_arr, rhs_arr, + acoef_arr, bx_arr, by_arr, bz_arr, flag_arr, dx, bx); + }); + } + } +} + +double check_norm(amrex::MultiFab const& phi, amrex::MultiFab const& exact) +{ + amrex::MultiFab mf(phi.boxArray(), phi.DistributionMap(), 1, 0); + amrex::MultiFab::Copy(mf, phi, 0, 0, 1, 0); + + amrex::MultiFab::Subtract(mf, exact, 0, 0, 1, 0); + + double L0norm = mf.norm0(); + double L1norm = mf.norm1(); + std::cout << " L0 norm:" << L0norm << ", L1 norm:" << L1norm << std::endl; + + return L1norm; +} + +BOOST_AUTO_TEST_CASE(p1_eb) +{ + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = true; + int const n_cell = amrpp.n_cell_; + + amrex::Geometry geom; + amrex::BoxArray grids; + amrex::DistributionMapping dmap; + + amrex::MultiFab solution; + amrex::MultiFab exact_solution; + amrex::MultiFab rhs; + amrex::MultiFab acoef; + amrex::Array bcoef; + amrex::MultiFab bcoef_eb; + + amrex::MultiFab robin_a; + amrex::MultiFab robin_b; + amrex::MultiFab robin_f; + + std::unique_ptr factory; + + int const max_grid_size = amrpp.max_grid_size_; + + amrex::RealBox rb( + { AMREX_D_DECL(-1.0, -1.0, -1.0) }, { AMREX_D_DECL(1.0, 1.0, 1.0) }); + amrex::Array is_periodic { AMREX_D_DECL(1, 1, 1) }; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + amrex::Box domain0(amrex::IntVect { AMREX_D_DECL(0, 0, 0) }, + amrex::IntVect { AMREX_D_DECL(n_cell - 1, n_cell - 1, n_cell - 1) }); + // geom.define(domain0); + geom.define(domain0, rb, amrex::CoordSys::cartesian, is_periodic); + + grids.define(domain0); + grids.maxSize(max_grid_size); + + // rotated box + int const max_coarsening_level = mlmgpp.max_coarsening_level_; + double const la = std::sqrt(2.0) / 2.0; + amrex::EB2::BoxIF box({ AMREX_D_DECL(-la, -la, -la * 0.75) }, + { AMREX_D_DECL(la, la, la * 1.25) }, true); + auto gshop = amrex::EB2::makeShop(amrex::EB2::translate( + amrex::EB2::rotate( + amrex::EB2::translate(box, { AMREX_D_DECL(-0.0, -0.0, -0.0) }), + std::atan(1.0) * 1.0, 2), + { AMREX_D_DECL(0.0, 0.0, 0.0) })); + amrex::EB2::Build(gshop, geom, 0, max_coarsening_level); + + amrex::IntVect ng = amrex::IntVect { 1 }; + + dmap.define(grids); + + amrex::EB2::IndexSpace const& eb_is = amrex::EB2::IndexSpace::top(); + amrex::EB2::Level const& eb_level = eb_is.getLevel(geom); + factory = std::make_unique(eb_level, geom, grids, + dmap, amrex::Vector { 2, 2, 2 }, amrex::EBSupport::full); + + solution.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + exact_solution.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + rhs.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + acoef.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + bcoef[idim].define( + amrex::convert(grids, amrex::IntVect::TheDimensionVector(idim)), + dmap, 1, 0, amrex::MFInfo(), *factory); + } + + solution.setVal(0.0); + rhs.setVal(0.0); + acoef.setVal(1.0); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + bcoef[idim].setVal(1.0); + } + + initProbABecLaplacian( + geom, factory, solution, rhs, exact_solution, acoef, bcoef); + + // std::cout << "initialize data ... \n"; + // initMeshandData(amrpp, mlmgpp, geom, grids, dmap, factory, solution, + // rhs, + // exact_solution, acoef, bcoef); + + std::cout << "construct the PDE ... \n"; + PeleRad::POneSingleEB rte(mlmgpp, geom, grids, dmap, factory, solution, rhs, + acoef, bcoef, robin_a, robin_b, robin_f); + + // std::cout << "solve the PDE ... \n"; + rte.solve(); + + auto eps = check_norm(solution, exact_solution); + eps /= static_cast(n_cell * n_cell * n_cell); + std::cout << "n_cell=" << n_cell << ", normalized L1 norm:" << eps + << std::endl; + + // plot results + if (write) + { + std::cout << "write the results ... \n"; + amrex::MultiFab const& vfrc = factory->getVolFrac(); + amrex::MultiFab plotmf(grids, dmap, 8, 0); + amrex::MultiFab::Copy(plotmf, solution, 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf, exact_solution, 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf, rhs, 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf, acoef, 0, 3, 1, 0); + amrex::MultiFab::Copy(plotmf, bcoef[0], 0, 4, 1, 0); + amrex::MultiFab::Copy(plotmf, bcoef[1], 0, 5, 1, 0); + amrex::MultiFab::Copy(plotmf, bcoef[2], 0, 6, 1, 0); + amrex::MultiFab::Copy(plotmf, vfrc, 0, 7, 1, 0); + + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteSingleLevelPlotfile(plot_file_name, plotmf, + { "phi", "exact", "rhs", "acoef", "bcoefx", "bcoefy", "bcoefz", + "vfrac" }, + geom, 0.0, 0); + + // for amrvis + /* + amrex::writeFabs(solution, "solution"); + amrex::writeFabs(bcoef, "bcoef"); + amrex::writeFabs(robin_a, "robin_a"); + amrex::writeFabs(robin_b, "robin_b"); + amrex::writeFabs(robin_f, "robin_f"); + */ + } + + BOOST_TEST(eps < 1e-2); +}