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: }