Line data Source code
1 :
2 : #include <AMReX_MLNodeLinOp.H>
3 : #include <AMReX_MLCellLinOp.H>
4 : #include <AMReX_MLNodeLap_K.H>
5 : #include <AMReX_MultiFabUtil.H>
6 : #include "Util/Color.H"
7 : #include "Set/Set.H"
8 : #include "Operator.H"
9 :
10 : using namespace amrex;
11 : namespace Operator {
12 :
13 : // constexpr amrex::IntVect AMREX_D_DECL(Operator<Grid::Node>::dx,Operator<Grid::Node>::dy,Operator<Grid::Node>::dz);
14 : constexpr amrex::IntVect AMREX_D_DECL(Operator<Grid::Cell>::dx, Operator<Grid::Cell>::dy, Operator<Grid::Cell>::dz);
15 :
16 734 : void Operator<Grid::Node>::Diagonal(bool recompute)
17 : {
18 : BL_PROFILE(Color::FG::Yellow + "Operator::Diagonal()" + Color::Reset);
19 734 : if (!recompute && m_diagonal_computed) return;
20 734 : m_diagonal_computed = true;
21 :
22 1522 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
23 : {
24 2395 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
25 : {
26 1607 : Diagonal(amrlev, mglev, *m_diag[amrlev][mglev]);
27 : }
28 : }
29 : }
30 :
31 0 : void Operator<Grid::Node>::Diagonal(int amrlev, int mglev, amrex::MultiFab& diag)
32 : {
33 : BL_PROFILE("Operator::Diagonal()");
34 : //Util::Message(INFO);
35 :
36 0 : int ncomp = diag.nComp();
37 0 : int nghost = 0;
38 :
39 0 : int sep = 2;
40 0 : int num = AMREX_D_TERM(sep, *sep, *sep);
41 0 : int cntr = 0;
42 :
43 0 : amrex::MultiFab x(m_diag[amrlev][mglev]->boxArray(), m_diag[amrlev][mglev]->DistributionMap(), ncomp, nghost);
44 0 : amrex::MultiFab Ax(m_diag[amrlev][mglev]->boxArray(), m_diag[amrlev][mglev]->DistributionMap(), ncomp, nghost);
45 :
46 0 : for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
47 : {
48 0 : const Box& bx = mfi.validbox();
49 0 : amrex::FArrayBox& diagfab = diag[mfi];
50 0 : amrex::FArrayBox& xfab = x[mfi];
51 0 : amrex::FArrayBox& Axfab = Ax[mfi];
52 :
53 0 : diagfab.setVal(0.0);
54 :
55 0 : for (int i = 0; i < num; i++)
56 : {
57 0 : for (int n = 0; n < ncomp; n++)
58 : {
59 0 : xfab.setVal(0.0);
60 0 : Axfab.setVal(0.0);
61 :
62 : //BL_PROFILE_VAR("Operator::Part1", part1);
63 0 : AMREX_D_TERM(for (int m1 = bx.loVect()[0]; m1 <= bx.hiVect()[0]; m1++),
64 : for (int m2 = bx.loVect()[1]; m2 <= bx.hiVect()[1]; m2++),
65 : for (int m3 = bx.loVect()[2]; m3 <= bx.hiVect()[2]; m3++))
66 : {
67 0 : amrex::IntVect m(AMREX_D_DECL(m1, m2, m3));
68 :
69 0 : if (m1 % sep == i / sep && m2 % sep == i % sep) xfab(m, n) = 1.0;
70 0 : else xfab(m, n) = 0.0;
71 : }
72 : //BL_PROFILE_VAR_STOP(part1);
73 :
74 : BL_PROFILE_VAR("Operator::Part2", part2);
75 0 : Util::Message(INFO, "Calling fapply...", cntr++);
76 0 : Fapply(amrlev, mglev, Ax, x);
77 : BL_PROFILE_VAR_STOP(part2);
78 :
79 : //BL_PROFILE_VAR("Operator::Part3", part3);
80 0 : Axfab.mult(xfab, n, n, 1);
81 0 : diagfab.plus(Axfab, n, n, 1);
82 : //BL_PROFILE_VAR_STOP(part3);
83 : }
84 : }
85 0 : }
86 0 : }
87 :
88 30468 : void Operator<Grid::Node>::Fsmooth(int amrlev, int mglev, amrex::MultiFab& x, const amrex::MultiFab& b) const
89 : {
90 : BL_PROFILE("Operator::Fsmooth()");
91 :
92 30468 : amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
93 30468 : domain.convert(amrex::IntVect::TheNodeVector());
94 :
95 30468 : int ncomp = b.nComp();
96 30468 : int nghost = 2; //b.nGrow();
97 :
98 :
99 30468 : amrex::MultiFab Ax(x.boxArray(), x.DistributionMap(), ncomp, nghost);
100 30468 : amrex::MultiFab Dx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
101 30468 : amrex::MultiFab Rx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
102 :
103 30468 : if (!m_diagonal_computed) Util::Abort(INFO, "Operator::Diagonal() must be called before using Fsmooth");
104 :
105 : // This is a JACOBI iteration, not Gauss-Seidel.
106 : // So we need to do twice the number of iterations to get the same behavior as GS.
107 91404 : for (int ctr = 0; ctr < 2; ctr++)
108 : {
109 60936 : Fapply(amrlev, mglev, Ax, x); // find Ax
110 :
111 60936 : amrex::MultiFab::Copy(Dx, x, 0, 0, ncomp, nghost); // Dx = x
112 60936 : amrex::MultiFab::Multiply(Dx, *m_diag[amrlev][mglev], 0, 0, ncomp, nghost); // Dx *= diag (Dx = x*diag)
113 :
114 60936 : amrex::MultiFab::Copy(Rx, Ax, 0, 0, ncomp, nghost); // Rx = Ax
115 60936 : amrex::MultiFab::Subtract(Rx, Dx, 0, 0, ncomp, nghost); // Rx -= Dx (Rx = Ax - Dx)
116 :
117 132048 : for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
118 : {
119 71112 : Box bx = mfi.grownnodaltilebox();
120 :
121 71112 : amrex::Array4<amrex::Real> const& xfab = x.array(mfi);
122 71112 : amrex::Array4<const amrex::Real> const& bfab = b.array(mfi);
123 71112 : amrex::Array4<const amrex::Real> const& Rxfab = Rx.array(mfi);
124 71112 : amrex::Array4<const amrex::Real> const& diagfab = (*m_diag[amrlev][mglev]).array(mfi);
125 :
126 :
127 236280 : for (int n = 0; n < ncomp; n++)
128 : {
129 165168 : amrex::ParallelFor(bx, [&] AMREX_GPU_DEVICE(int i, int j, int k)
130 : {
131 :
132 : // Skip ghost cells outside problem domain
133 78202256 : if (!domain.contains(i,j,k))
134 : {
135 : //continue;
136 : }
137 31295088 : else if ( !bx.strictly_contains(i,j,k))
138 : {
139 4539008 : xfab(i, j, k, n) = 0.0;
140 : //continue;
141 : }
142 : else
143 : {
144 174153504 : xfab(i,j,k,n) = (1. - m_omega) * xfab(i,j,k, n) + m_omega * (bfab(i,j,k, n) - Rxfab(i,j,k, n)) / diagfab(i,j,k,n);
145 : }
146 78202256 : });
147 : }
148 60936 : }
149 60936 : amrex::Geometry geom = m_geom[amrlev][mglev];
150 60936 : x.setMultiGhost(true);
151 60936 : x.FillBoundary(geom.periodicity());
152 60936 : nodalSync(amrlev, mglev, x);
153 : }
154 30468 : }
155 :
156 170147 : void Operator<Grid::Node>::normalize(int amrlev, int mglev, MultiFab& a_x) const
157 : {
158 170147 : if (!m_diagonal_computed)
159 0 : Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
160 :
161 170147 : a_x.divide(*m_diag[amrlev][mglev],0,getNComp(),2);
162 :
163 170147 : a_x.setMultiGhost(true);
164 170147 : a_x.FillBoundaryAndSync(Geom(amrlev,mglev).periodicity());
165 170147 : }
166 :
167 0 : Operator<Grid::Node>::Operator(const Vector<Geometry>& a_geom,
168 : const Vector<BoxArray>& a_grids,
169 : const Vector<DistributionMapping>& a_dmap,
170 : const LPInfo& a_info,
171 0 : const Vector<FabFactory<FArrayBox> const*>& a_factory)
172 : {
173 : BL_PROFILE("Operator::Operator()");
174 0 : Util::Message(INFO);
175 :
176 0 : if (!(a_grids[0].ixType() == amrex::IndexType::TheNodeType()))
177 0 : Util::Abort(INFO, "Operator must be defined using CELL CENTERED boxarrays.");
178 :
179 0 : define(a_geom, a_grids, a_dmap, a_info, a_factory);
180 0 : }
181 :
182 734 : Operator<Grid::Node>::~Operator()
183 734 : {}
184 :
185 734 : void Operator<Grid::Node>::define(const Vector<Geometry>& a_geom,
186 : const Vector<BoxArray>& a_grids,
187 : const Vector<DistributionMapping>& a_dmap,
188 : const LPInfo& a_info,
189 : const Vector<FabFactory<FArrayBox> const*>& a_factory)
190 : {
191 : BL_PROFILE("Operator::~Operator()");
192 :
193 : // Make sure we're not trying to parallelize in vain.
194 734 : if (amrex::ParallelDescriptor::NProcs() > a_grids[0].size())
195 : {
196 0 : Util::Warning(INFO, "There are more processors than there are boxes in the amrlev=0 boxarray!!\n",
197 0 : "(NProcs = ", amrex::ParallelDescriptor::NProcs(), ", a_grids[0].size() = ", a_grids[0].size(), ")\n",
198 : "You should decrease max_grid_size or you will not get proper scaling!");
199 : }
200 :
201 : // This makes sure grids are node-centered;
202 734 : Vector<BoxArray> cc_grids = a_grids;
203 1532 : for (auto& ba : cc_grids) {
204 798 : ba.enclosedCells();
205 : }
206 :
207 734 : MLNodeLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
208 :
209 734 : int nghost = 2;
210 : // Resize the multifab containing the operator diagonal
211 734 : m_diag.resize(m_num_amr_levels);
212 1522 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
213 : {
214 788 : m_diag[amrlev].resize(m_num_mg_levels[amrlev]);
215 :
216 2395 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
217 : {
218 4821 : m_diag[amrlev][mglev].reset(new MultiFab(amrex::convert(m_grids[amrlev][mglev], amrex::IntVect::TheNodeVector()),
219 4821 : m_dmap[amrlev][mglev], getNComp(), nghost));
220 : }
221 : }
222 :
223 : // We need to instantiate the m_lobc objects.
224 : // WE DO NOT USE THEM - our BCs are implemented differently.
225 : // But they need to be the right size or the code will segfault.
226 734 : m_lobc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
227 734 : m_hibc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
228 734 : }
229 :
230 54 : void Operator<Grid::Node>::fixUpResidualMask(int amrlev, iMultiFab& resmsk)
231 : {
232 : BL_PROFILE("Operator::fixUpResidualMask()");
233 :
234 54 : if (!m_masks_built) buildMasks();
235 :
236 54 : const iMultiFab& cfmask = *m_nd_fine_mask[amrlev];
237 :
238 : #ifdef _OPENMP
239 : #pragma omp parallel if (Gpu::notInLaunchRegion())
240 : #endif
241 161 : for (MFIter mfi(resmsk, TilingIfNotGPU()); mfi.isValid(); ++mfi)
242 : {
243 107 : const Box& bx = mfi.tilebox();
244 107 : Array4<int> const& rmsk = resmsk.array(mfi);
245 107 : Array4<int const> const& fmsk = cfmask.const_array(mfi);
246 70552 : AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
247 : {
248 : if (fmsk(i,j,k) == amrex::nodelap_detail::crse_fine_node) rmsk(i,j,k) = 1;
249 : });
250 54 : }
251 54 : }
252 :
253 734 : void Operator<Grid::Node>::prepareForSolve()
254 : {
255 : BL_PROFILE("Operator::prepareForSolve()");
256 734 : MLNodeLinOp::prepareForSolve();
257 734 : buildMasks();
258 734 : averageDownCoeffs();
259 734 : Diagonal(true);
260 734 : }
261 :
262 6951 : void Operator<Grid::Node>::restriction(int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const
263 : {
264 : BL_PROFILE("Operator::restriction()");
265 :
266 6951 : applyBC(amrlev, cmglev - 1, fine, BCMode::Homogeneous, StateMode::Solution);
267 :
268 6951 : amrex::Box cdomain = m_geom[amrlev][cmglev].growPeriodicDomain(1);
269 6951 : cdomain = cdomain.convert(amrex::IntVect::TheNodeVector());
270 :
271 6951 : bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
272 6951 : MultiFab cfine;
273 6951 : if (need_parallel_copy) {
274 0 : const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
275 0 : cfine.define(ba, fine.DistributionMap(), fine.nComp(), fine.nGrow());
276 0 : }
277 :
278 6951 : MultiFab* pcrse = (need_parallel_copy) ? &cfine : &crse;
279 :
280 13902 : for (MFIter mfi(*pcrse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
281 : {
282 6951 : const Box& bx = mfi.grownnodaltilebox(-1,1) & cdomain;
283 :
284 6951 : amrex::Array4<const amrex::Real> const& fdata = fine.array(mfi);
285 6951 : amrex::Array4<amrex::Real> const& cdata = pcrse->array(mfi);
286 :
287 6951 : const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
288 :
289 :
290 23715 : for (int n = 0; n < crse.nComp(); n++)
291 : {
292 : // I,J,K == coarse coordinates
293 : // i,j,k == fine coordinates
294 16764 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
295 488320 : int i = 2 * I, j = 2 * J, k = 2 * K;
296 :
297 488320 : if ((I == lo.x || I == hi.x) &&
298 227624 : (J == lo.y || J == hi.y) &&
299 135744 : (K == lo.z || K == hi.z)) // Corner
300 : {
301 304200 : cdata(I, J, K, n) = fdata(i, j, k, n);
302 : }
303 386920 : else if ((J == lo.y || J == hi.y) &&
304 127824 : (K == lo.z || K == hi.z)) // X edge
305 : {
306 381540 : cdata(I, J, K, n) = 0.25 * fdata(i - 1, j, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i + 1, j, k, n);
307 : }
308 310612 : else if ((K == lo.z || K == hi.z) &&
309 233338 : (I == lo.x || I == hi.x)) // Y edge
310 : {
311 373540 : cdata(I, J, K, n) = 0.25 * fdata(i, j - 1, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j + 1, k, n);
312 : }
313 235904 : else if ((I == lo.x || I == hi.x) &&
314 51516 : (J == lo.y || J == hi.y)) // Z edge
315 : {
316 171720 : cdata(I, J, K, n) = 0.25 * fdata(i, j, k - 1, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j, k + 1, n);
317 : }
318 201560 : else if (I == lo.x || I == hi.x) // X face
319 : {
320 17172 : cdata(I, J, K, n) =
321 51516 : (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
322 51516 : + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
323 68688 : + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
324 : }
325 184388 : else if (J == lo.y || J == hi.y) // Y face
326 : {
327 17172 : cdata(I, J, K, n) =
328 51516 : (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
329 51516 : + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
330 68688 : + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
331 : }
332 167216 : else if (K == lo.z || K == hi.z) // Z face
333 : {
334 158630 : cdata(I, J, K, n) =
335 475890 : (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
336 475890 : + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
337 634520 : + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
338 : }
339 : else // Interior
340 8586 : cdata(I, J, K, n) =
341 34344 : (fdata(i - 1, j - 1, k - 1, n) + fdata(i - 1, j - 1, k + 1, n) + fdata(i - 1, j + 1, k - 1, n) + fdata(i - 1, j + 1, k + 1, n) +
342 34344 : fdata(i + 1, j - 1, k - 1, n) + fdata(i + 1, j - 1, k + 1, n) + fdata(i + 1, j + 1, k - 1, n) + fdata(i + 1, j + 1, k + 1, n)) / 64.0
343 8586 : +
344 34344 : (fdata(i, j - 1, k - 1, n) + fdata(i, j - 1, k + 1, n) + fdata(i, j + 1, k - 1, n) + fdata(i, j + 1, k + 1, n) +
345 34344 : fdata(i - 1, j, k - 1, n) + fdata(i + 1, j, k - 1, n) + fdata(i - 1, j, k + 1, n) + fdata(i + 1, j, k + 1, n) +
346 34344 : fdata(i - 1, j - 1, k, n) + fdata(i - 1, j + 1, k, n) + fdata(i + 1, j - 1, k, n) + fdata(i + 1, j + 1, k, n)) / 32.0
347 8586 : +
348 25758 : (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
349 25758 : fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
350 8586 : +
351 17172 : fdata(i, j, k, n) / 8.0;
352 488320 : });
353 : }
354 6951 : }
355 :
356 6951 : if (need_parallel_copy) {
357 0 : crse.ParallelCopy(cfine);
358 : }
359 :
360 6951 : crse.setMultiGhost(true);
361 6951 : crse.FillBoundary(Geom(amrlev,cmglev).periodicity());
362 6951 : nodalSync(amrlev, cmglev, crse);
363 6951 : }
364 :
365 6951 : void Operator<Grid::Node>::interpolation(int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const
366 : {
367 : BL_PROFILE("Operator::interpolation()");
368 6951 : amrex::Box fdomain = m_geom[amrlev][fmglev].growPeriodicDomain(2);
369 6951 : fdomain.convert(amrex::IntVect::TheNodeVector());
370 :
371 6951 : bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
372 6951 : MultiFab cfine;
373 6951 : const MultiFab* cmf = &crse;
374 6951 : if (need_parallel_copy) {
375 0 : const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
376 0 : cfine.define(ba, fine.DistributionMap(), crse.nComp(), crse.nGrow());
377 0 : cfine.ParallelCopy(crse);
378 0 : cmf = &cfine;
379 0 : }
380 :
381 13902 : for (MFIter mfi(fine, false); mfi.isValid(); ++mfi)
382 : {
383 6951 : Box fine_bx = mfi.validbox() & fdomain;
384 :
385 6951 : const Box& course_bx = amrex::coarsen(fine_bx, 2);
386 6951 : const Box& tmpbx = amrex::refine(course_bx, 2);
387 6951 : FArrayBox tmpfab;
388 6951 : tmpfab.resize(tmpbx, fine.nComp());
389 6951 : tmpfab.setVal(0.0);
390 6951 : const amrex::FArrayBox& crsefab = (*cmf)[mfi];
391 :
392 13902 : amrex::Array4<const amrex::Real> const& cdata = crsefab.const_array();
393 6951 : amrex::Array4<amrex::Real> const& fdata = tmpfab.array();
394 :
395 23715 : for (int n = 0; n < crse.nComp(); n++)
396 : {
397 : // I,J,K == coarse coordinates
398 : // i,j,k == fine coordinates
399 16764 : amrex::ParallelFor(fine_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
400 :
401 1930484 : int I = i / 2, J = j / 2, K = k / 2;
402 :
403 1930484 : if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) // Coincident
404 1441152 : fdata(i, j, k, n) = cdata(I, J, K, n);
405 1450100 : else if (j % 2 == 0 && k % 2 == 0) // X Edge
406 1471408 : fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n));
407 1082248 : else if (k % 2 == 0 && i % 2 == 0) // Y Edge
408 1463088 : fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n));
409 716476 : else if (i % 2 == 0 && j % 2 == 0) // Z Edge
410 618192 : fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n));
411 561928 : else if (i % 2 == 0) // X Face
412 309096 : fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n) +
413 309096 : cdata(I, J, K + 1, n) + cdata(I, J + 1, K + 1, n));
414 458896 : else if (j % 2 == 0) // Y Face
415 309096 : fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n) +
416 309096 : cdata(I + 1, J, K, n) + cdata(I + 1, J, K + 1, n));
417 355864 : else if (k % 2 == 0) // Z Face
418 861528 : fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n) +
419 861528 : cdata(I, J + 1, K, n) + cdata(I + 1, J + 1, K, n));
420 : else // Center
421 137376 : fdata(i, j, k, n) = 0.125 * (cdata(I, J, K, n) +
422 206064 : cdata(I + 1, J, K, n) + cdata(I, J + 1, K, n) + cdata(I, J, K + 1, n) +
423 206064 : cdata(I, J + 1, K + 1, n) + cdata(I + 1, J, K + 1, n) + cdata(I + 1, J + 1, K, n) +
424 137376 : cdata(I + 1, J + 1, K + 1, n));
425 :
426 1930484 : });
427 : }
428 6951 : fine[mfi].plus(tmpfab, fine_bx, fine_bx, 0, 0, fine.nComp());
429 13902 : }
430 :
431 6951 : fine.setMultiGhost(true);
432 6951 : fine.FillBoundary(Geom(amrlev,fmglev).periodicity());
433 6951 : nodalSync(amrlev, fmglev, fine);
434 6951 : }
435 :
436 54 : void Operator<Grid::Node>::averageDownSolutionRHS(int camrlev, MultiFab& crse_sol, MultiFab& /*crse_rhs*/,
437 : const MultiFab& fine_sol, const MultiFab& /*fine_rhs*/)
438 : {
439 : BL_PROFILE("Operator::averageDownSolutionRHS()");
440 54 : const auto& amrrr = AMRRefRatio(camrlev);
441 54 : amrex::average_down(fine_sol, crse_sol, 0, crse_sol.nComp(), amrrr);
442 :
443 54 : if (isSingular(0))
444 : {
445 0 : Util::Abort(INFO, "Singular operators not supported!");
446 : }
447 :
448 54 : }
449 :
450 0 : void Operator<Grid::Node>::realFillBoundary(MultiFab& phi, const Geometry& geom)
451 : {
452 : Util::RealFillBoundary(phi, geom);
453 0 : }
454 :
455 216995 : void Operator<Grid::Node>::applyBC(int amrlev, int mglev, MultiFab& phi, BCMode/* bc_mode*/,
456 : amrex::MLLinOp::StateMode /**/, bool skip_fillboundary) const
457 : {
458 : BL_PROFILE("Operator::applyBC()");
459 :
460 216995 : const Geometry& geom = m_geom[amrlev][mglev];
461 :
462 216995 : if (!skip_fillboundary) {
463 : //phi.FillBoundary(geom.periodicity());
464 : //phi.setMultiGhost(true);
465 216995 : phi.FillBoundaryAndSync(geom.periodicity());
466 : }
467 216995 : }
468 :
469 : const amrex::FArrayBox&
470 0 : Operator<Grid::Node>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
471 : {
472 : BL_PROFILE("Operator::GetFab()");
473 0 : Util::Message(INFO);
474 0 : return m_a_coeffs[num][amrlev][mglev][mfi];
475 : }
476 :
477 0 : void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
478 : {
479 : BL_PROFILE("Operator::RegisterNewFab()");
480 0 : Util::Message(INFO);
481 : /// \todo assertions here
482 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
483 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
484 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
485 : {
486 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
487 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
488 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
489 0 : m_dmap[amrlev][mglev],
490 0 : input[amrlev].nComp(),
491 0 : input[amrlev].nGrow());
492 :
493 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
494 0 : input[amrlev], 0, 0,
495 0 : input[amrlev].nComp(),
496 0 : input[amrlev].nGrow());
497 : }
498 0 : m_num_a_fabs++;
499 0 : }
500 :
501 :
502 0 : void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
503 : {
504 : BL_PROFILE("Operator::RegisterNewFab()");
505 0 : Util::Message(INFO);
506 : /// \todo assertions here
507 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
508 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
509 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
510 : {
511 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
512 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
513 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
514 0 : m_dmap[amrlev][mglev],
515 0 : input[amrlev]->nComp(),
516 0 : input[amrlev]->nGrow());
517 :
518 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
519 0 : *input[amrlev], 0, 0,
520 0 : input[amrlev]->nComp(),
521 0 : input[amrlev]->nGrow());
522 : }
523 0 : m_num_a_fabs++;
524 0 : }
525 :
526 416 : void Operator<Grid::Node>::reflux(int crse_amrlev,
527 : MultiFab& res, const MultiFab& /*crse_sol*/, const MultiFab& /*crse_rhs*/,
528 : MultiFab& fine_res, MultiFab& /*fine_sol*/, const MultiFab& /*fine_rhs*/) const
529 : {
530 : BL_PROFILE("Operator::Elastic::reflux()");
531 :
532 416 : int ncomp = AMREX_SPACEDIM;
533 :
534 416 : amrex::Box cdomain(m_geom[crse_amrlev][0].growPeriodicDomain(2));
535 416 : cdomain.convert(amrex::IntVect::TheNodeVector());
536 :
537 416 : const Geometry& cgeom = m_geom[crse_amrlev][0];
538 :
539 416 : const BoxArray& fba = fine_res.boxArray();
540 416 : const DistributionMapping& fdm = fine_res.DistributionMap();
541 :
542 416 : MultiFab fine_res_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
543 416 : fine_res_for_coarse.ParallelCopy(res, 0, 0, ncomp, 0, 0, cgeom.periodicity());
544 :
545 416 : applyBC(crse_amrlev + 1, 0, fine_res, BCMode::Inhomogeneous, StateMode::Solution);
546 :
547 : /// \todo Replace with Enum
548 : // const int coarse_coarse_node = 0;
549 416 : const int coarse_fine_node = 1;
550 416 : const int fine_fine_node = 2;
551 :
552 416 : amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
553 416 : nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
554 :
555 2400 : for (MFIter mfi(fine_res_for_coarse, false); mfi.isValid(); ++mfi)
556 : {
557 1984 : const Box& bx = mfi.grownnodaltilebox(-1,1) & cdomain;
558 :
559 1984 : amrex::Array4<const int> const& nmask = nodemask.array(mfi);
560 : //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
561 :
562 1984 : amrex::Array4<amrex::Real> const& cdata = fine_res_for_coarse.array(mfi);
563 1984 : amrex::Array4<const amrex::Real> const& fdata = fine_res.array(mfi);
564 :
565 1984 : const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
566 :
567 5952 : for (int n = 0; n < fine_res.nComp(); n++)
568 : {
569 : // I,J,K == coarse coordinates
570 : // i,j,k == fine coordinates
571 3968 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
572 696530 : int i = I * 2, j = J * 2, k = K * 2;
573 :
574 1652600 : if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
575 : {
576 528572 : if ((I == lo.x || I == hi.x) &&
577 11224 : (J == lo.y || J == hi.y) &&
578 1174 : (K == lo.z || K == hi.z)) // Corner
579 3522 : cdata(I, J, K, n) = fdata(i, j, k, n);
580 527398 : else if ((J == lo.y || J == hi.y) &&
581 6698 : (K == lo.z || K == hi.z)) // X edge
582 33490 : cdata(I, J, K, n) = 0.25 * fdata(i - 1, j, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i + 1, j, k, n);
583 520700 : else if ((K == lo.z || K == hi.z) &&
584 520700 : (I == lo.x || I == hi.x)) // Y edge
585 50250 : cdata(I, J, K, n) = 0.25 * fdata(i, j - 1, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j + 1, k, n);
586 510650 : else if ((I == lo.x || I == hi.x) &&
587 0 : (J == lo.y || J == hi.y)) // Z edge
588 0 : cdata(I, J, K, n) = 0.25 * fdata(i, j, k - 1, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j, k + 1, n);
589 510650 : else if (I == lo.x || I == hi.x) // X face
590 0 : cdata(I, J, K, n) =
591 0 : (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
592 0 : + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
593 0 : + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
594 510650 : else if (J == lo.y || J == hi.y) // Y face
595 0 : cdata(I, J, K, n) =
596 0 : (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
597 0 : + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
598 0 : + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
599 510650 : else if (K == lo.z || K == hi.z) // Z face
600 510650 : cdata(I, J, K, n) =
601 1531950 : (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
602 1531950 : + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
603 2042600 : + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
604 : else // Interior
605 0 : cdata(I, J, K, n) =
606 0 : (fdata(i - 1, j - 1, k - 1, n) + fdata(i - 1, j - 1, k + 1, n) + fdata(i - 1, j + 1, k - 1, n) + fdata(i - 1, j + 1, k + 1, n) +
607 0 : fdata(i + 1, j - 1, k - 1, n) + fdata(i + 1, j - 1, k + 1, n) + fdata(i + 1, j + 1, k - 1, n) + fdata(i + 1, j + 1, k + 1, n)) / 64.0
608 0 : +
609 0 : (fdata(i, j - 1, k - 1, n) + fdata(i, j - 1, k + 1, n) + fdata(i, j + 1, k - 1, n) + fdata(i, j + 1, k + 1, n) +
610 0 : fdata(i - 1, j, k - 1, n) + fdata(i + 1, j, k - 1, n) + fdata(i - 1, j, k + 1, n) + fdata(i + 1, j, k + 1, n) +
611 0 : fdata(i - 1, j - 1, k, n) + fdata(i - 1, j + 1, k, n) + fdata(i + 1, j - 1, k, n) + fdata(i + 1, j + 1, k, n)) / 32.0
612 0 : +
613 0 : (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
614 0 : fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
615 0 : +
616 0 : fdata(i, j, k, n) / 8.0;
617 : }
618 :
619 696530 : });
620 : }
621 416 : }
622 :
623 : // Copy the fine residual restricted onto the coarse grid
624 : // into the final residual.
625 416 : res.ParallelCopy(fine_res_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
626 :
627 : // Sync up ghost nodes
628 416 : res.setMultiGhost(true);
629 416 : res.FillBoundaryAndSync(Geom(crse_amrlev).periodicity());
630 416 : nodalSync(crse_amrlev, 0, res);
631 :
632 832 : return;
633 416 : }
634 :
635 : void
636 7426 : Operator<Grid::Node>::solutionResidual(int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b,
637 : const MultiFab* /*crse_bcdata*/)
638 : {
639 7426 : const int mglev = 0;
640 7426 : const int ncomp = b.nComp();
641 7426 : apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution);
642 7426 : MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, 2);
643 7426 : resid.setMultiGhost(true);
644 7426 : resid.FillBoundaryAndSync(Geom(amrlev).periodicity());
645 7426 : }
646 :
647 : void
648 20406 : Operator<Grid::Node>::correctionResidual(int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b,
649 : BCMode /*bc_mode*/, const MultiFab* /*crse_bcdata*/)
650 : {
651 20406 : resid.setVal(0.0);
652 20406 : apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction);
653 20406 : int ncomp = b.nComp();
654 20406 : MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, resid.nGrow());
655 20406 : resid.setMultiGhost(true);
656 20406 : resid.FillBoundaryAndSync(Geom(amrlev).periodicity());
657 20406 : }
658 :
659 :
660 :
661 :
662 0 : Operator<Grid::Cell>::Operator()
663 : {
664 0 : m_ixtype = amrex::IntVect::TheCellVector();
665 0 : }
666 :
667 : void
668 0 : Operator<Grid::Cell>::define(amrex::Vector<amrex::Geometry> a_geom,
669 : const amrex::Vector<amrex::BoxArray>& a_grids,
670 : const amrex::Vector<amrex::DistributionMapping>& a_dmap,
671 : BC::BC<Set::Scalar>& a_bc,
672 : const amrex::LPInfo& a_info,
673 : const amrex::Vector<amrex::FabFactory<amrex::FArrayBox> const*>& a_factory)
674 : {
675 0 : m_bc = &a_bc;
676 :
677 0 : std::array<int, AMREX_SPACEDIM> is_periodic = m_bc->IsPeriodic();
678 :
679 0 : MLCellLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
680 :
681 0 : Util::Warning(INFO, "This section of code has not been tested.");
682 0 : for (int n = 0; n < getNComp(); n++)
683 : {
684 0 : m_lobc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
685 : is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
686 : is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
687 0 : m_hibc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
688 : is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
689 : is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
690 : }
691 :
692 0 : for (int ilev = 0; ilev < a_geom.size(); ++ilev)
693 0 : setLevelBC(ilev, nullptr);
694 :
695 0 : }
696 :
697 :
698 : void
699 0 : Operator<Grid::Cell>::prepareForSolve()
700 : {
701 0 : MLCellLinOp::prepareForSolve();
702 0 : }
703 :
704 0 : Operator<Grid::Cell>::BndryCondLoc::BndryCondLoc(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm)
705 0 : : bcond(ba, dm),
706 0 : bcloc(ba, dm)
707 : {
708 0 : }
709 :
710 : void
711 0 : Operator<Grid::Cell>::BndryCondLoc::setLOBndryConds(const amrex::Geometry& /*geom*/, const amrex::Real* /*dx*/,
712 : const amrex::Array<BCType, AMREX_SPACEDIM>& /*lobc*/,
713 : const amrex::Array<BCType, AMREX_SPACEDIM>& /*hibc*/,
714 : int /*ratio*/, const amrex::RealVect& /*a_loc*/)
715 : {
716 0 : Util::Warning(INFO, "This code has not been properlyt tested");
717 0 : }
718 :
719 :
720 : void
721 0 : Operator<Grid::Cell>::averageDownCoeffs()
722 : {
723 0 : for (int i = 0; i < m_num_a_fabs; i++)
724 : {
725 0 : for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
726 : {
727 0 : auto& fine_a_coeffs = m_a_coeffs[i][amrlev];
728 0 : averageDownCoeffsSameAmrLevel(fine_a_coeffs);
729 : }
730 0 : averageDownCoeffsSameAmrLevel(m_a_coeffs[i][0]);
731 : }
732 0 : }
733 :
734 : void
735 0 : Operator<Grid::Cell>::averageDownCoeffsSameAmrLevel(amrex::Vector<amrex::MultiFab>& a)
736 : {
737 0 : int nmglevs = a.size();
738 0 : for (int mglev = 1; mglev < nmglevs; ++mglev)
739 : {
740 0 : amrex::average_down(a[mglev - 1], a[mglev], 0, a[0].nComp(), mg_coarsen_ratio);
741 : }
742 0 : }
743 :
744 :
745 :
746 : const amrex::FArrayBox&
747 0 : Operator<Grid::Cell>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
748 : {
749 0 : return m_a_coeffs[num][amrlev][mglev][mfi];
750 : }
751 :
752 :
753 : void
754 0 : Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
755 : {
756 : /// \todo assertions here
757 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
758 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
759 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
760 : {
761 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
762 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
763 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
764 0 : m_dmap[amrlev][mglev],
765 0 : input[amrlev].nComp(),
766 0 : input[amrlev].nGrow());
767 :
768 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
769 0 : input[amrlev], 0, 0,
770 0 : input[amrlev].nComp(),
771 0 : input[amrlev].nGrow());
772 : }
773 0 : m_num_a_fabs++;
774 0 : }
775 :
776 :
777 : void
778 0 : Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
779 : {
780 : /// \todo assertions here
781 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
782 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
783 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
784 : {
785 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
786 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
787 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
788 0 : m_dmap[amrlev][mglev],
789 0 : input[amrlev]->nComp(),
790 0 : input[amrlev]->nGrow());
791 :
792 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
793 0 : *input[amrlev], 0, 0,
794 0 : input[amrlev]->nComp(),
795 0 : input[amrlev]->nGrow());
796 : }
797 0 : m_num_a_fabs++;
798 0 : }
799 :
800 :
801 : }
|
