Actual source code: MyTest.cpp

      1: #include "MyTest.H"

      3: #include <AMReX_MLEBABecLap.H>
      4: #include <AMReX_ParmParse.H>
      5: #include <AMReX_MultiFabUtil.H>
      6: #include <AMReX_EBMultiFabUtil.H>
      7: #include <AMReX_PlotFileUtil.H>
      8: #include <AMReX_EB2.H>

     10: #include <cmath>

     12: using namespace amrex;

     14: MyTest::MyTest()
     15: {
     16:   readParameters();

     18:   initGrids();

     20:   initializeEB();

     22:   initData();
     23: }

     25: void MyTest::solve()
     26: {
     27:   for (int ilev = 0; ilev <= max_level; ++ilev) {
     28:     const MultiFab &vfrc = factory[ilev]->getVolFrac();
     29:     MultiFab        v(vfrc.boxArray(), vfrc.DistributionMap(), 1, 0, MFInfo(), *factory[ilev]);
     30:     MultiFab::Copy(v, vfrc, 0, 0, 1, 0);
     31:     amrex::EB_set_covered(v, 1.0);
     32:     amrex::Print() << "vfrc min = " << v.min(0) << std::endl;
     33:   }

     35:   std::array<LinOpBCType, AMREX_SPACEDIM> mlmg_lobc;
     36:   std::array<LinOpBCType, AMREX_SPACEDIM> mlmg_hibc;
     37:   for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
     38:     if (geom[0].isPeriodic(idim)) {
     39:       mlmg_lobc[idim] = LinOpBCType::Periodic;
     40:       mlmg_hibc[idim] = LinOpBCType::Periodic;
     41:     } else {
     42:       mlmg_lobc[idim] = LinOpBCType::Dirichlet;
     43:       mlmg_hibc[idim] = LinOpBCType::Dirichlet;
     44:     }
     45:   }

     47:   LPInfo info;
     48:   info.setMaxCoarseningLevel(max_coarsening_level);

     50:   MLEBABecLap mleb(geom, grids, dmap, info, amrex::GetVecOfConstPtrs(factory));
     51:   mleb.setMaxOrder(linop_maxorder);

     53:   mleb.setDomainBC(mlmg_lobc, mlmg_hibc);

     55:   for (int ilev = 0; ilev <= max_level; ++ilev) mleb.setLevelBC(ilev, &phi[ilev]);

     57:   mleb.setScalars(scalars[0], scalars[1]);

     59:   for (int ilev = 0; ilev <= max_level; ++ilev) {
     60:     mleb.setACoeffs(ilev, acoef[ilev]);
     61:     mleb.setBCoeffs(ilev, amrex::GetArrOfConstPtrs(bcoef[ilev]));
     62:   }

     64:   if (eb_is_dirichlet) {
     65:     for (int ilev = 0; ilev <= max_level; ++ilev) mleb.setEBDirichlet(ilev, phi[ilev], bcoef_eb[ilev]);
     66:   }

     68:   MLMG mlmg(mleb);
     69:   mlmg.setMaxIter(max_iter);
     70:   mlmg.setMaxFmgIter(max_fmg_iter);
     71:   mlmg.setBottomMaxIter(max_bottom_iter);
     72:   mlmg.setBottomTolerance(bottom_reltol);
     73:   mlmg.setVerbose(verbose);
     74:   mlmg.setBottomVerbose(bottom_verbose);
     75:   if (use_hypre) mlmg.setBottomSolver(MLMG::BottomSolver::hypre);
     76:   if (use_petsc) mlmg.setBottomSolver(MLMG::BottomSolver::petsc);
     77:   const Real tol_rel = reltol;
     78:   const Real tol_abs = 0.0;
     79:   mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs);
     80: }

     82: void MyTest::writePlotfile()
     83: {
     84:   Vector<MultiFab> plotmf(max_level + 1);
     85:   for (int ilev = 0; ilev <= max_level; ++ilev) {
     86:     const MultiFab &vfrc = factory[ilev]->getVolFrac();
     87:     plotmf[ilev].define(grids[ilev], dmap[ilev], 2, 0);
     88:     MultiFab::Copy(plotmf[ilev], phi[ilev], 0, 0, 1, 0);
     89:     MultiFab::Copy(plotmf[ilev], vfrc, 0, 1, 1, 0);
     90:   }
     91:   WriteMultiLevelPlotfile(plot_file_name, max_level + 1, amrex::GetVecOfConstPtrs(plotmf), {"phi", "vfrac"}, geom, 0.0, Vector<int>(max_level + 1, 0), Vector<IntVect>(max_level, IntVect{2}));
     92: }

     94: void MyTest::readParameters()
     95: {
     96:   ParmParse pp;
     97:   pp.query("max_level", max_level);
     98:   pp.query("n_cell", n_cell);
     99:   pp.query("max_grid_size", max_grid_size);
    100:   pp.query("is_periodic", is_periodic);
    101:   pp.query("eb_is_dirichlet", eb_is_dirichlet);

    103:   pp.query("plot_file", plot_file_name);

    105:   scalars.resize(2);
    106:   if (is_periodic) {
    107:     scalars[0] = 0.0;
    108:     scalars[1] = 1.0;
    109:   } else {
    110:     scalars[0] = 1.0;
    111:     scalars[1] = 1.0;
    112:   }
    113:   pp.queryarr("scalars", scalars);

    115:   pp.query("verbose", verbose);
    116:   pp.query("bottom_verbose", bottom_verbose);
    117:   pp.query("max_iter", max_iter);
    118:   pp.query("max_fmg_iter", max_fmg_iter);
    119:   pp.query("max_bottom_iter", max_bottom_iter);
    120:   pp.query("bottom_reltol", bottom_reltol);
    121:   pp.query("reltol", reltol);
    122:   pp.query("linop_maxorder", linop_maxorder);
    123:   pp.query("max_coarsening_level", max_coarsening_level);
    124: #ifdef AMREX_USE_HYPRE
    125:   pp.query("use_hypre", use_hypre);
    126: #endif
    127: #ifdef AMREX_USE_PETSC
    128:   pp.query("use_petsc", use_petsc);
    129: #endif
    130: }

    132: void MyTest::initGrids()
    133: {
    134:   int nlevels = max_level + 1;
    135:   geom.resize(nlevels);
    136:   grids.resize(nlevels);

    138:   RealBox                         rb({AMREX_D_DECL(0., 0., 0.)}, {AMREX_D_DECL(1., 1., 1.)});
    139:   std::array<int, AMREX_SPACEDIM> isperiodic{AMREX_D_DECL(is_periodic, is_periodic, is_periodic)};
    140:   Geometry::Setup(&rb, 0, isperiodic.data());
    141:   Box domain0(IntVect{AMREX_D_DECL(0, 0, 0)}, IntVect{AMREX_D_DECL(n_cell - 1, n_cell - 1, n_cell - 1)});
    142:   Box domain = domain0;
    143:   for (int ilev = 0; ilev < nlevels; ++ilev) {
    144:     geom[ilev].define(domain);
    145:     domain.refine(ref_ratio);
    146:   }

    148:   domain = domain0;
    149:   for (int ilev = 0; ilev < nlevels; ++ilev) {
    150:     grids[ilev].define(domain);
    151:     grids[ilev].maxSize(max_grid_size);
    152:     domain.grow(-n_cell / 4); // fine level cover the middle of the coarse domain
    153:     domain.refine(ref_ratio);
    154:   }
    155: }

    157: void MyTest::initData()
    158: {
    159:   int nlevels = max_level + 1;
    160:   dmap.resize(nlevels);
    161:   factory.resize(nlevels);
    162:   phi.resize(nlevels);
    163:   rhs.resize(nlevels);
    164:   acoef.resize(nlevels);
    165:   bcoef.resize(nlevels);
    166:   bcoef_eb.resize(nlevels);

    168:   for (int ilev = 0; ilev < nlevels; ++ilev) {
    169:     dmap[ilev].define(grids[ilev]);
    170:     const EB2::IndexSpace &eb_is    = EB2::IndexSpace::top();
    171:     const EB2::Level      &eb_level = eb_is.getLevel(geom[ilev]);
    172:     factory[ilev].reset(new EBFArrayBoxFactory(eb_level, geom[ilev], grids[ilev], dmap[ilev], {2, 2, 2}, EBSupport::full));

    174:     phi[ilev].define(grids[ilev], dmap[ilev], 1, 1, MFInfo(), *factory[ilev]);
    175:     rhs[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
    176:     acoef[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
    177:     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) bcoef[ilev][idim].define(amrex::convert(grids[ilev], IntVect::TheDimensionVector(idim)), dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
    178:     if (eb_is_dirichlet) {
    179:       bcoef_eb[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
    180:       bcoef_eb[ilev].setVal(1.0);
    181:     }

    183:     phi[ilev].setVal(0.0);
    184:     rhs[ilev].setVal(0.0);
    185:     acoef[ilev].setVal(1.0);
    186:     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) bcoef[ilev][idim].setVal(1.0);

    188:     const auto dx = geom[ilev].CellSizeArray();

    190:     if (is_periodic) {
    191:       const Real pi = 4.0 * std::atan(1.0);

    193:       for (MFIter mfi(rhs[ilev]); mfi.isValid(); ++mfi) {
    194:         const Box          &bx  = mfi.fabbox();
    195:         Array4<Real> const &fab = rhs[ilev].array(mfi);
    196:         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    197:           Real rx      = (i + 0.5) * dx[0];
    198:           Real ry      = (j + 0.5) * dx[1];
    199:           fab(i, j, k) = std::sin(rx * 2. * pi + 43.5) * std::sin(ry * 2. * pi + 89.);
    200:         });
    201:       }
    202:     } else if (eb_is_dirichlet) {
    203:       phi[ilev].setVal(10.0);
    204:       phi[ilev].setVal(0.0, 0, 1, 0); // set interior
    205:     } else {
    206:       // Initialize Dirichlet boundary
    207:       for (MFIter mfi(phi[ilev]); mfi.isValid(); ++mfi) {
    208:         const Box          &bx  = mfi.fabbox();
    209:         Array4<Real> const &fab = phi[ilev].array(mfi);
    210:         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    211:           Real rx      = (i + 0.5) * dx[0];
    212:           Real ry      = (j + 0.5) * dx[1];
    213:           fab(i, j, k) = std::sqrt(0.5) * (rx + ry);
    214:         });
    215:       }
    216:       phi[ilev].setVal(0.0, 0, 1, 0); // set interior to zero
    217:     }
    218:   }
    219: }