1 #include "MyTest.H"
2
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>
9
10 #include <cmath>
11
12 using namespace amrex;
13
MyTest()14 MyTest::MyTest()
15 {
16 readParameters();
17
18 initGrids();
19
20 initializeEB();
21
22 initData();
23 }
24
solve()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 }
34
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 }
46
47 LPInfo info;
48 info.setMaxCoarseningLevel(max_coarsening_level);
49
50 MLEBABecLap mleb(geom, grids, dmap, info, amrex::GetVecOfConstPtrs(factory));
51 mleb.setMaxOrder(linop_maxorder);
52
53 mleb.setDomainBC(mlmg_lobc, mlmg_hibc);
54
55 for (int ilev = 0; ilev <= max_level; ++ilev) mleb.setLevelBC(ilev, &phi[ilev]);
56
57 mleb.setScalars(scalars[0], scalars[1]);
58
59 for (int ilev = 0; ilev <= max_level; ++ilev) {
60 mleb.setACoeffs(ilev, acoef[ilev]);
61 mleb.setBCoeffs(ilev, amrex::GetArrOfConstPtrs(bcoef[ilev]));
62 }
63
64 if (eb_is_dirichlet) {
65 for (int ilev = 0; ilev <= max_level; ++ilev) mleb.setEBDirichlet(ilev, phi[ilev], bcoef_eb[ilev]);
66 }
67
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 }
81
writePlotfile()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 }
93
readParameters()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);
102
103 pp.query("plot_file", plot_file_name);
104
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);
114
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 }
131
initGrids()132 void MyTest::initGrids()
133 {
134 int nlevels = max_level + 1;
135 geom.resize(nlevels);
136 grids.resize(nlevels);
137
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 }
147
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 }
156
initData()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);
167
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));
173
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 }
182
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);
187
188 const auto dx = geom[ilev].CellSizeArray();
189
190 if (is_periodic) {
191 const Real pi = 4.0 * std::atan(1.0);
192
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 }
220