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
25 void
solve()26 MyTest::solve ()
27 {
28 for (int ilev = 0; ilev <= max_level; ++ilev) {
29 const MultiFab& vfrc = factory[ilev]->getVolFrac();
30 MultiFab v(vfrc.boxArray(), vfrc.DistributionMap(), 1, 0,
31 MFInfo(), *factory[ilev]);
32 MultiFab::Copy(v, vfrc, 0, 0, 1, 0);
33 amrex::EB_set_covered(v, 1.0);
34 amrex::Print() << "vfrc min = " << v.min(0) << std::endl;
35 }
36
37 std::array<LinOpBCType,AMREX_SPACEDIM> mlmg_lobc;
38 std::array<LinOpBCType,AMREX_SPACEDIM> mlmg_hibc;
39 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
40 if (geom[0].isPeriodic(idim)) {
41 mlmg_lobc[idim] = LinOpBCType::Periodic;
42 mlmg_hibc[idim] = LinOpBCType::Periodic;
43 } else {
44 mlmg_lobc[idim] = LinOpBCType::Dirichlet;
45 mlmg_hibc[idim] = LinOpBCType::Dirichlet;
46 }
47 }
48
49 LPInfo info;
50 info.setMaxCoarseningLevel(max_coarsening_level);
51
52 MLEBABecLap mleb (geom, grids, dmap, info, amrex::GetVecOfConstPtrs(factory));
53 mleb.setMaxOrder(linop_maxorder);
54
55 mleb.setDomainBC(mlmg_lobc, mlmg_hibc);
56
57 for (int ilev = 0; ilev <= max_level; ++ilev) {
58 mleb.setLevelBC(ilev, &phi[ilev]);
59 }
60
61 mleb.setScalars(scalars[0], scalars[1]);
62
63 for (int ilev = 0; ilev <= max_level; ++ilev) {
64 mleb.setACoeffs(ilev, acoef[ilev]);
65 mleb.setBCoeffs(ilev, amrex::GetArrOfConstPtrs(bcoef[ilev]));
66 }
67
68 if (eb_is_dirichlet) {
69 for (int ilev = 0; ilev <= max_level; ++ilev) {
70 mleb.setEBDirichlet(ilev, phi[ilev], bcoef_eb[ilev]);
71 }
72 }
73
74 MLMG mlmg(mleb);
75 mlmg.setMaxIter(max_iter);
76 mlmg.setMaxFmgIter(max_fmg_iter);
77 mlmg.setBottomMaxIter(max_bottom_iter);
78 mlmg.setBottomTolerance(bottom_reltol);
79 mlmg.setVerbose(verbose);
80 mlmg.setBottomVerbose(bottom_verbose);
81 if (use_hypre) mlmg.setBottomSolver(MLMG::BottomSolver::hypre);
82 if (use_petsc) mlmg.setBottomSolver(MLMG::BottomSolver::petsc);
83 const Real tol_rel = reltol;
84 const Real tol_abs = 0.0;
85 mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs);
86 }
87
88 void
writePlotfile()89 MyTest::writePlotfile ()
90 {
91 Vector<MultiFab> plotmf(max_level+1);
92 for (int ilev = 0; ilev <= max_level; ++ilev) {
93 const MultiFab& vfrc = factory[ilev]->getVolFrac();
94 plotmf[ilev].define(grids[ilev],dmap[ilev],2,0);
95 MultiFab::Copy(plotmf[ilev], phi[ilev], 0, 0, 1, 0);
96 MultiFab::Copy(plotmf[ilev], vfrc, 0, 1, 1, 0);
97 }
98 WriteMultiLevelPlotfile(plot_file_name, max_level+1,
99 amrex::GetVecOfConstPtrs(plotmf),
100 {"phi","vfrac"},
101 geom, 0.0, Vector<int>(max_level+1,0),
102 Vector<IntVect>(max_level,IntVect{2}));
103
104 }
105
106 void
readParameters()107 MyTest::readParameters ()
108 {
109 ParmParse pp;
110 pp.query("max_level", max_level);
111 pp.query("n_cell", n_cell);
112 pp.query("max_grid_size", max_grid_size);
113 pp.query("is_periodic", is_periodic);
114 pp.query("eb_is_dirichlet", eb_is_dirichlet);
115
116 pp.query("plot_file", plot_file_name);
117
118 scalars.resize(2);
119 if (is_periodic) {
120 scalars[0] = 0.0;
121 scalars[1] = 1.0;
122 } else {
123 scalars[0] = 1.0;
124 scalars[1] = 1.0;
125 }
126 pp.queryarr("scalars", scalars);
127
128 pp.query("verbose", verbose);
129 pp.query("bottom_verbose", bottom_verbose);
130 pp.query("max_iter", max_iter);
131 pp.query("max_fmg_iter", max_fmg_iter);
132 pp.query("max_bottom_iter", max_bottom_iter);
133 pp.query("bottom_reltol", bottom_reltol);
134 pp.query("reltol", reltol);
135 pp.query("linop_maxorder", linop_maxorder);
136 pp.query("max_coarsening_level", max_coarsening_level);
137 #ifdef AMREX_USE_HYPRE
138 pp.query("use_hypre", use_hypre);
139 #endif
140 #ifdef AMREX_USE_PETSC
141 pp.query("use_petsc",use_petsc);
142 #endif
143 }
144
145 void
initGrids()146 MyTest::initGrids ()
147 {
148 int nlevels = max_level + 1;
149 geom.resize(nlevels);
150 grids.resize(nlevels);
151
152 RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(1.,1.,1.)});
153 std::array<int,AMREX_SPACEDIM> isperiodic{AMREX_D_DECL(is_periodic,is_periodic,is_periodic)};
154 Geometry::Setup(&rb, 0, isperiodic.data());
155 Box domain0(IntVect{AMREX_D_DECL(0,0,0)}, IntVect{AMREX_D_DECL(n_cell-1,n_cell-1,n_cell-1)});
156 Box domain = domain0;
157 for (int ilev = 0; ilev < nlevels; ++ilev)
158 {
159 geom[ilev].define(domain);
160 domain.refine(ref_ratio);
161 }
162
163 domain = domain0;
164 for (int ilev = 0; ilev < nlevels; ++ilev)
165 {
166 grids[ilev].define(domain);
167 grids[ilev].maxSize(max_grid_size);
168 domain.grow(-n_cell/4); // fine level cover the middle of the coarse domain
169 domain.refine(ref_ratio);
170 }
171 }
172
173 void
initData()174 MyTest::initData ()
175 {
176 int nlevels = max_level + 1;
177 dmap.resize(nlevels);
178 factory.resize(nlevels);
179 phi.resize(nlevels);
180 rhs.resize(nlevels);
181 acoef.resize(nlevels);
182 bcoef.resize(nlevels);
183 bcoef_eb.resize(nlevels);
184
185 for (int ilev = 0; ilev < nlevels; ++ilev)
186 {
187 dmap[ilev].define(grids[ilev]);
188 const EB2::IndexSpace& eb_is = EB2::IndexSpace::top();
189 const EB2::Level& eb_level = eb_is.getLevel(geom[ilev]);
190 factory[ilev].reset(new EBFArrayBoxFactory(eb_level, geom[ilev], grids[ilev], dmap[ilev],
191 {2,2,2}, EBSupport::full));
192
193 phi[ilev].define(grids[ilev], dmap[ilev], 1, 1, MFInfo(), *factory[ilev]);
194 rhs[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
195 acoef[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
196 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
197 bcoef[ilev][idim].define(amrex::convert(grids[ilev],IntVect::TheDimensionVector(idim)),
198 dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
199 }
200 if (eb_is_dirichlet) {
201 bcoef_eb[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
202 bcoef_eb[ilev].setVal(1.0);
203 }
204
205 phi[ilev].setVal(0.0);
206 rhs[ilev].setVal(0.0);
207 acoef[ilev].setVal(1.0);
208 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
209 bcoef[ilev][idim].setVal(1.0);
210 }
211
212 const auto dx = geom[ilev].CellSizeArray();
213
214 if (is_periodic)
215 {
216 const Real pi = 4.0*std::atan(1.0);
217
218 for (MFIter mfi(rhs[ilev]); mfi.isValid(); ++mfi)
219 {
220 const Box& bx = mfi.fabbox();
221 Array4<Real> const& fab = rhs[ilev].array(mfi);
222 amrex::ParallelFor(bx,
223 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
224 {
225 Real rx = (i+0.5)*dx[0];
226 Real ry = (j+0.5)*dx[1];
227 fab(i,j,k) = std::sin(rx*2.*pi + 43.5)*std::sin(ry*2.*pi + 89.);
228 });
229 }
230 }
231 else if (eb_is_dirichlet)
232 {
233 phi[ilev].setVal(10.0);
234 phi[ilev].setVal(0.0, 0, 1, 0); // set interior
235 }
236 else
237 {
238 // Initialize Dirichlet boundary
239 for (MFIter mfi(phi[ilev]); mfi.isValid(); ++mfi)
240 {
241 const Box& bx = mfi.fabbox();
242 Array4<Real> const& fab = phi[ilev].array(mfi);
243 amrex::ParallelFor(bx,
244 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
245 {
246 Real rx = (i+0.5)*dx[0];
247 Real ry = (j+0.5)*dx[1];
248 fab(i,j,k) = std::sqrt(0.5)*(rx + ry);
249 });
250 }
251 phi[ilev].setVal(0.0, 0, 1, 0); // set interior to zero
252 }
253
254 }
255 }
256
257