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 748 : void Operator<Grid::Node>::Diagonal(bool recompute)
17 : {
18 : BL_PROFILE(Color::FG::Yellow + "Operator::Diagonal()" + Color::Reset);
19 748 : if (!recompute && m_diagonal_computed) return;
20 748 : m_diagonal_computed = true;
21 :
22 1665 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
23 : {
24 2744 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
25 : {
26 1827 : 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 96375 : 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 96375 : amrex::Box domain(m_geom[amrlev][mglev].Domain());
93 :
94 96375 : int ncomp = b.nComp();
95 96375 : int nghost = 2; //b.nGrow();
96 :
97 :
98 96375 : amrex::MultiFab Ax(x.boxArray(), x.DistributionMap(), ncomp, nghost);
99 96375 : amrex::MultiFab Dx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
100 96375 : amrex::MultiFab Rx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
101 :
102 96375 : if (!m_diagonal_computed) Util::Abort(INFO, "Operator::Diagonal() must be called before using Fsmooth");
103 :
104 : // This is a JACOBI iteration, not Gauss-Seidel.
105 : // So we need to do twice the number of iterations to get the same behavior as GS.
106 289111 : for (int ctr = 0; ctr < 2; ctr++)
107 : {
108 192748 : Fapply(amrlev, mglev, Ax, x); // find Ax
109 :
110 192738 : amrex::MultiFab::Copy(Dx, x, 0, 0, ncomp, nghost); // Dx = x
111 192738 : amrex::MultiFab::Multiply(Dx, *m_diag[amrlev][mglev], 0, 0, ncomp, nghost); // Dx *= diag (Dx = x*diag)
112 :
113 192738 : amrex::MultiFab::Copy(Rx, Ax, 0, 0, ncomp, nghost); // Rx = Ax
114 192738 : amrex::MultiFab::Subtract(Rx, Dx, 0, 0, ncomp, nghost); // Rx -= Dx (Rx = Ax - Dx)
115 :
116 629049 : for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
117 : {
118 436313 : const Box& bx = mfi.validbox();
119 436313 : amrex::FArrayBox& xfab = x[mfi];
120 436313 : const amrex::FArrayBox& bfab = b[mfi];
121 : //const amrex::FArrayBox &Axfab = Ax[mfi];
122 436313 : const amrex::FArrayBox& Rxfab = Rx[mfi];
123 436313 : const amrex::FArrayBox& diagfab = (*m_diag[amrlev][mglev])[mfi];
124 :
125 1390335 : for (int n = 0; n < ncomp; n++)
126 : {
127 399518390 : AMREX_D_TERM(for (int m1 = bx.loVect()[0] - 2; m1 <= bx.hiVect()[0] + 2; m1++),
128 : for (int m2 = bx.loVect()[1] - 2; m2 <= bx.hiVect()[1] + 2; m2++),
129 : for (int m3 = bx.loVect()[2] - 2; m3 <= bx.hiVect()[2] + 2; m3++))
130 : {
131 :
132 365197649 : amrex::IntVect m(AMREX_D_DECL(m1, m2, m3));
133 :
134 : // Skip ghost cells outside problem domain
135 810690541 : if (AMREX_D_TERM(m[0] < domain.loVect()[0], ||
136 : m[1] < domain.loVect()[1], ||
137 197518109 : m[2] < domain.loVect()[2])) continue;
138 551035304 : if (AMREX_D_TERM(m[0] > domain.hiVect()[0] + 1, ||
139 : m[1] > domain.hiVect()[1] + 1, ||
140 63330126 : m[2] > domain.hiVect()[2] + 1)) continue;
141 :
142 835036331 : if (AMREX_D_TERM(m[0] == bx.loVect()[0] - nghost || m[0] == bx.hiVect()[0] + nghost, ||
143 : m[1] == bx.loVect()[1] - nghost || m[1] == bx.hiVect()[1] + nghost, ||
144 : m[2] == bx.loVect()[2] - nghost || m[2] == bx.hiVect()[2] + nghost))
145 : {
146 28961498 : xfab(m, n) = 0.0;
147 28961498 : continue;
148 : }
149 :
150 : //xfab(m,n) = xfab(m,n) + omega*(bfab(m,n) - Axfab(m,n))/diagfab(m,n);
151 1006077230 : xfab(m, n) = (1. - m_omega) * xfab(m, n) + m_omega * (bfab(m, n) - Rxfab(m, n)) / diagfab(m, n);
152 : }
153 : }
154 192736 : }
155 : }
156 96363 : amrex::Geometry geom = m_geom[amrlev][mglev];
157 96363 : realFillBoundary(x, geom);
158 96363 : nodalSync(amrlev, mglev, x);
159 96363 : }
160 :
161 694666 : void Operator<Grid::Node>::normalize(int amrlev, int mglev, MultiFab& a_x) const
162 : {
163 : BL_PROFILE("Operator::normalize()");
164 694666 : amrex::Box domain(m_geom[amrlev][mglev].Domain());
165 694666 : domain.convert(amrex::IntVect::TheNodeVector());
166 :
167 694666 : int ncomp = getNComp();
168 694666 : int nghost = 1; //x.nGrow();
169 :
170 694666 : if (!m_diagonal_computed)
171 0 : Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
172 :
173 1389332 : for (MFIter mfi(a_x, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
174 : {
175 :
176 694666 : Box bx = mfi.tilebox();
177 694666 : bx.grow(nghost);
178 694666 : bx = bx & domain;
179 :
180 694666 : amrex::Array4<amrex::Real> const& x = a_x.array(mfi);
181 694666 : amrex::Array4<const amrex::Real> const& diag = m_diag[amrlev][mglev]->array(mfi);
182 :
183 2618509 : for (int n = 0; n < ncomp; n++)
184 : {
185 1923843 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
186 :
187 92356362 : x(i, j, k) /= diag(i, j, k);
188 :
189 46178181 : });
190 : }
191 694666 : }
192 694666 : }
193 :
194 0 : Operator<Grid::Node>::Operator(const Vector<Geometry>& a_geom,
195 : const Vector<BoxArray>& a_grids,
196 : const Vector<DistributionMapping>& a_dmap,
197 : const LPInfo& a_info,
198 0 : const Vector<FabFactory<FArrayBox> const*>& a_factory)
199 : {
200 : BL_PROFILE("Operator::Operator()");
201 0 : Util::Message(INFO);
202 :
203 0 : if (!(a_grids[0].ixType() == amrex::IndexType::TheNodeType()))
204 0 : Util::Abort(INFO, "Operator must be defined using CELL CENTERED boxarrays.");
205 :
206 0 : define(a_geom, a_grids, a_dmap, a_info, a_factory);
207 0 : }
208 :
209 729 : Operator<Grid::Node>::~Operator()
210 729 : {}
211 :
212 749 : void Operator<Grid::Node>::define(const Vector<Geometry>& a_geom,
213 : const Vector<BoxArray>& a_grids,
214 : const Vector<DistributionMapping>& a_dmap,
215 : const LPInfo& a_info,
216 : const Vector<FabFactory<FArrayBox> const*>& a_factory)
217 : {
218 : BL_PROFILE("Operator::~Operator()");
219 :
220 : // Make sure we're not trying to parallelize in vain.
221 749 : if (amrex::ParallelDescriptor::NProcs() > a_grids[0].size())
222 : {
223 0 : Util::Warning(INFO, "There are more processors than there are boxes in the amrlev=0 boxarray!!\n",
224 0 : "(NProcs = ", amrex::ParallelDescriptor::NProcs(), ", a_grids[0].size() = ", a_grids[0].size(), ")\n",
225 : "You should decrease max_grid_size or you will not get proper scaling!");
226 : }
227 :
228 : // This makes sure grids are node-centered;
229 749 : Vector<BoxArray> cc_grids = a_grids;
230 1688 : for (auto& ba : cc_grids) {
231 939 : ba.enclosedCells();
232 : }
233 :
234 749 : MLNodeLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
235 :
236 749 : int nghost = 2;
237 : // Resize the multifab containing the operator diagonal
238 749 : m_diag.resize(m_num_amr_levels);
239 1668 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
240 : {
241 919 : m_diag[amrlev].resize(m_num_mg_levels[amrlev]);
242 :
243 2748 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
244 : {
245 5487 : m_diag[amrlev][mglev].reset(new MultiFab(amrex::convert(m_grids[amrlev][mglev], amrex::IntVect::TheNodeVector()),
246 5487 : m_dmap[amrlev][mglev], getNComp(), nghost));
247 : }
248 : }
249 :
250 : // We need to instantiate the m_lobc objects.
251 : // WE DO NOT USE THEM - our BCs are implemented differently.
252 : // But they need to be the right size or the code will segfault.
253 749 : m_lobc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
254 749 : m_hibc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
255 749 : }
256 :
257 170 : void Operator<Grid::Node>::fixUpResidualMask(int amrlev, iMultiFab& resmsk)
258 : {
259 : BL_PROFILE("Operator::fixUpResidualMask()");
260 :
261 170 : if (!m_masks_built) buildMasks();
262 :
263 170 : const iMultiFab& cfmask = *m_nd_fine_mask[amrlev];
264 :
265 : #ifdef _OPENMP
266 : #pragma omp parallel if (Gpu::notInLaunchRegion())
267 : #endif
268 2114 : for (MFIter mfi(resmsk, TilingIfNotGPU()); mfi.isValid(); ++mfi)
269 : {
270 1944 : const Box& bx = mfi.tilebox();
271 1944 : Array4<int> const& rmsk = resmsk.array(mfi);
272 1944 : Array4<int const> const& fmsk = cfmask.const_array(mfi);
273 1034873 : AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
274 : {
275 : if (fmsk(i,j,k) == amrex::nodelap_detail::crse_fine_node) rmsk(i,j,k) = 1;
276 : });
277 170 : }
278 170 : }
279 :
280 748 : void Operator<Grid::Node>::prepareForSolve()
281 : {
282 : BL_PROFILE("Operator::prepareForSolve()");
283 748 : MLNodeLinOp::prepareForSolve();
284 748 : buildMasks();
285 748 : averageDownCoeffs();
286 748 : Diagonal(true);
287 747 : }
288 :
289 20123 : void Operator<Grid::Node>::restriction(int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const
290 : {
291 : BL_PROFILE("Operator::restriction()");
292 :
293 20123 : applyBC(amrlev, cmglev - 1, fine, BCMode::Homogeneous, StateMode::Solution);
294 :
295 20123 : amrex::Box cdomain = m_geom[amrlev][cmglev].Domain();
296 20123 : cdomain.convert(amrex::IntVect::TheNodeVector());
297 :
298 20123 : bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
299 20123 : MultiFab cfine;
300 20123 : if (need_parallel_copy) {
301 345 : const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
302 345 : cfine.define(ba, fine.DistributionMap(), fine.nComp(), fine.nGrow());
303 345 : }
304 :
305 20123 : MultiFab* pcrse = (need_parallel_copy) ? &cfine : &crse;
306 :
307 41845 : for (MFIter mfi(*pcrse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
308 : {
309 21722 : const Box& bx = mfi.tilebox();
310 :
311 21722 : amrex::Array4<const amrex::Real> const& fdata = fine.array(mfi);
312 21722 : amrex::Array4<amrex::Real> const& cdata = pcrse->array(mfi);
313 :
314 21722 : const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
315 :
316 :
317 75303 : for (int n = 0; n < crse.nComp(); n++)
318 : {
319 : // I,J,K == coarse coordinates
320 : // i,j,k == fine coordinates
321 53581 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
322 1812907 : int i = 2 * I, j = 2 * J, k = 2 * K;
323 :
324 1812907 : if ((I == lo.x || I == hi.x) &&
325 751354 : (J == lo.y || J == hi.y) &&
326 446020 : (K == lo.z || K == hi.z)) // Corner
327 : {
328 968952 : cdata(I, J, K, n) = fdata(i, j, k, n);
329 : }
330 1489923 : else if ((J == lo.y || J == hi.y) &&
331 467418 : (K == lo.z || K == hi.z)) // X edge
332 : {
333 1376880 : 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);
334 : }
335 1214547 : else if ((K == lo.z || K == hi.z) &&
336 873492 : (I == lo.x || I == hi.x)) // Y edge
337 : {
338 1181640 : 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);
339 : }
340 978219 : else if ((I == lo.x || I == hi.x) &&
341 192042 : (J == lo.y || J == hi.y)) // Z edge
342 : {
343 615180 : 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);
344 : }
345 855183 : else if (I == lo.x || I == hi.x) // X face
346 : {
347 69006 : cdata(I, J, K, n) =
348 207018 : (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
349 207018 : + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
350 276024 : + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
351 : }
352 786177 : else if (J == lo.y || J == hi.y) // Y face
353 : {
354 69006 : cdata(I, J, K, n) =
355 207018 : (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
356 207018 : + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
357 276024 : + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
358 : }
359 717171 : else if (K == lo.z || K == hi.z) // Z face
360 : {
361 637164 : cdata(I, J, K, n) =
362 1911492 : (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
363 1911492 : + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
364 2548656 : + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
365 : }
366 : else // Interior
367 80007 : cdata(I, J, K, n) =
368 320028 : (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) +
369 320028 : 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
370 80007 : +
371 320028 : (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) +
372 320028 : 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) +
373 320028 : 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
374 80007 : +
375 240021 : (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
376 240021 : fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
377 80007 : +
378 160014 : fdata(i, j, k, n) / 8.0;
379 1812907 : });
380 : }
381 20123 : }
382 :
383 20123 : if (need_parallel_copy) {
384 345 : crse.ParallelCopy(cfine);
385 : }
386 :
387 20123 : amrex::Geometry geom = m_geom[amrlev][cmglev];
388 20123 : realFillBoundary(crse, geom);
389 20123 : nodalSync(amrlev, cmglev, crse);
390 20123 : }
391 :
392 20121 : void Operator<Grid::Node>::interpolation(int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const
393 : {
394 : BL_PROFILE("Operator::interpolation()");
395 : //int nghost = getNGrow();
396 40242 : amrex::Box fdomain = m_geom[amrlev][fmglev].Domain(); fdomain.convert(amrex::IntVect::TheNodeVector());
397 :
398 20121 : bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
399 20121 : MultiFab cfine;
400 20121 : const MultiFab* cmf = &crse;
401 20121 : if (need_parallel_copy) {
402 0 : const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
403 0 : cfine.define(ba, fine.DistributionMap(), crse.nComp(), crse.nGrow());
404 0 : cfine.ParallelCopy(crse);
405 0 : cmf = &cfine;
406 0 : }
407 :
408 41829 : for (MFIter mfi(fine, false); mfi.isValid(); ++mfi)
409 : {
410 21708 : const Box& fine_bx = mfi.validbox() & fdomain;
411 21708 : const Box& course_bx = amrex::coarsen(fine_bx, 2);
412 21708 : const Box& tmpbx = amrex::refine(course_bx, 2);
413 21708 : FArrayBox tmpfab;
414 21708 : tmpfab.resize(tmpbx, fine.nComp());
415 21708 : tmpfab.setVal(0.0);
416 21708 : const amrex::FArrayBox& crsefab = (*cmf)[mfi];
417 :
418 43416 : amrex::Array4<const amrex::Real> const& cdata = crsefab.const_array();
419 21708 : amrex::Array4<amrex::Real> const& fdata = tmpfab.array();
420 :
421 75248 : for (int n = 0; n < crse.nComp(); n++)
422 : {
423 : // I,J,K == coarse coordinates
424 : // i,j,k == fine coordinates
425 53540 : amrex::ParallelFor(fine_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
426 :
427 7553132 : int I = i / 2, J = j / 2, K = k / 2;
428 :
429 7553132 : if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) // Coincident
430 5435112 : fdata(i, j, k, n) = cdata(I, J, K, n);
431 5741428 : else if (j % 2 == 0 && k % 2 == 0) // X Edge
432 5686208 : fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n));
433 4319876 : else if (k % 2 == 0 && i % 2 == 0) // Y Edge
434 5608240 : fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n));
435 2917816 : else if (i % 2 == 0 && j % 2 == 0) // Z Edge
436 2480160 : fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n));
437 2297776 : else if (i % 2 == 0) // X Face
438 1297584 : fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n) +
439 1297584 : cdata(I, J, K + 1, n) + cdata(I, J + 1, K + 1, n));
440 1865248 : else if (j % 2 == 0) // Y Face
441 1297584 : fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n) +
442 1297584 : cdata(I + 1, J, K, n) + cdata(I + 1, J, K + 1, n));
443 1432720 : else if (k % 2 == 0) // Z Face
444 3379728 : fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n) +
445 3379728 : cdata(I, J + 1, K, n) + cdata(I + 1, J + 1, K, n));
446 : else // Center
447 612288 : fdata(i, j, k, n) = 0.125 * (cdata(I, J, K, n) +
448 918432 : cdata(I + 1, J, K, n) + cdata(I, J + 1, K, n) + cdata(I, J, K + 1, n) +
449 918432 : cdata(I, J + 1, K + 1, n) + cdata(I + 1, J, K + 1, n) + cdata(I + 1, J + 1, K, n) +
450 612288 : cdata(I + 1, J + 1, K + 1, n));
451 :
452 7553132 : });
453 : }
454 21708 : fine[mfi].plus(tmpfab, fine_bx, fine_bx, 0, 0, fine.nComp());
455 41829 : }
456 :
457 20121 : amrex::Geometry geom = m_geom[amrlev][fmglev];
458 20121 : realFillBoundary(fine, geom);
459 20121 : nodalSync(amrlev, fmglev, fine);
460 20121 : }
461 :
462 188 : void Operator<Grid::Node>::averageDownSolutionRHS(int camrlev, MultiFab& crse_sol, MultiFab& /*crse_rhs*/,
463 : const MultiFab& fine_sol, const MultiFab& /*fine_rhs*/)
464 : {
465 : BL_PROFILE("Operator::averageDownSolutionRHS()");
466 188 : const auto& amrrr = AMRRefRatio(camrlev);
467 188 : amrex::average_down(fine_sol, crse_sol, 0, crse_sol.nComp(), amrrr);
468 :
469 188 : if (isSingular(0))
470 : {
471 0 : Util::Abort(INFO, "Singular operators not supported!");
472 : }
473 :
474 188 : }
475 :
476 1054467 : void Operator<Grid::Node>::realFillBoundary(MultiFab& phi, const Geometry& geom)
477 : {
478 : Util::RealFillBoundary(phi, geom);
479 1054467 : }
480 :
481 837441 : void Operator<Grid::Node>::applyBC(int amrlev, int mglev, MultiFab& phi, BCMode/* bc_mode*/,
482 : amrex::MLLinOp::StateMode /**/, bool skip_fillboundary) const
483 : {
484 : BL_PROFILE("Operator::applyBC()");
485 :
486 837441 : const Geometry& geom = m_geom[amrlev][mglev];
487 :
488 837441 : if (!skip_fillboundary) {
489 :
490 837441 : realFillBoundary(phi, geom);
491 : }
492 837441 : }
493 :
494 : const amrex::FArrayBox&
495 0 : Operator<Grid::Node>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
496 : {
497 : BL_PROFILE("Operator::GetFab()");
498 0 : Util::Message(INFO);
499 0 : return m_a_coeffs[num][amrlev][mglev][mfi];
500 : }
501 :
502 0 : void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<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 :
527 0 : void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
528 : {
529 : BL_PROFILE("Operator::RegisterNewFab()");
530 0 : Util::Message(INFO);
531 : /// \todo assertions here
532 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
533 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
534 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
535 : {
536 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
537 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
538 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
539 0 : m_dmap[amrlev][mglev],
540 0 : input[amrlev]->nComp(),
541 0 : input[amrlev]->nGrow());
542 :
543 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
544 0 : *input[amrlev], 0, 0,
545 0 : input[amrlev]->nComp(),
546 0 : input[amrlev]->nGrow());
547 : }
548 0 : m_num_a_fabs++;
549 0 : }
550 :
551 2180 : void Operator<Grid::Node>::reflux(int crse_amrlev,
552 : MultiFab& res, const MultiFab& /*crse_sol*/, const MultiFab& /*crse_rhs*/,
553 : MultiFab& fine_res, MultiFab& /*fine_sol*/, const MultiFab& /*fine_rhs*/) const
554 : {
555 : BL_PROFILE("Operator::Elastic::reflux()");
556 :
557 2180 : int ncomp = AMREX_SPACEDIM;
558 :
559 2180 : amrex::Box cdomain(m_geom[crse_amrlev][0].Domain());
560 2180 : cdomain.convert(amrex::IntVect::TheNodeVector());
561 :
562 2180 : const Geometry& cgeom = m_geom[crse_amrlev][0];
563 :
564 2180 : const BoxArray& fba = fine_res.boxArray();
565 2180 : const DistributionMapping& fdm = fine_res.DistributionMap();
566 :
567 2180 : MultiFab fine_res_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
568 2180 : fine_res_for_coarse.ParallelCopy(res, 0, 0, ncomp, 0, 0, cgeom.periodicity());
569 :
570 2180 : applyBC(crse_amrlev + 1, 0, fine_res, BCMode::Inhomogeneous, StateMode::Solution);
571 :
572 : /// \todo Replace with Enum
573 : // const int coarse_coarse_node = 0;
574 2180 : const int coarse_fine_node = 1;
575 2180 : const int fine_fine_node = 2;
576 :
577 2180 : amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
578 2180 : nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
579 :
580 4360 : amrex::iMultiFab cellmask(amrex::convert(amrex::coarsen(fba, 2), amrex::IntVect::TheCellVector()), fdm, 1, 2);
581 2180 : cellmask.ParallelCopy(*m_cc_fine_mask[crse_amrlev], 0, 0, 1, 1, 1, cgeom.periodicity());
582 :
583 37547 : for (MFIter mfi(fine_res_for_coarse, false); mfi.isValid(); ++mfi)
584 : {
585 35367 : const Box& bx = mfi.validbox();
586 :
587 35367 : amrex::Array4<const int> const& nmask = nodemask.array(mfi);
588 : //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
589 :
590 35367 : amrex::Array4<amrex::Real> const& cdata = fine_res_for_coarse.array(mfi);
591 35367 : amrex::Array4<const amrex::Real> const& fdata = fine_res.array(mfi);
592 :
593 35367 : const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
594 :
595 106192 : for (int n = 0; n < fine_res.nComp(); n++)
596 : {
597 : // I,J,K == coarse coordinates
598 : // i,j,k == fine coordinates
599 70825 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
600 3514703 : int i = I * 2, j = J * 2, k = K * 2;
601 :
602 7720960 : if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
603 : {
604 3514703 : if ((I == lo.x || I == hi.x) &&
605 33210 : (J == lo.y || J == hi.y) &&
606 2752 : (K == lo.z || K == hi.z)) // Corner
607 8256 : cdata(I, J, K, n) = fdata(i, j, k, n);
608 3511951 : else if ((J == lo.y || J == hi.y) &&
609 46834 : (K == lo.z || K == hi.z)) // X edge
610 234170 : 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);
611 3465117 : else if ((K == lo.z || K == hi.z) &&
612 3109584 : (I == lo.x || I == hi.x)) // Y edge
613 152290 : 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);
614 3434659 : else if ((I == lo.x || I == hi.x) &&
615 0 : (J == lo.y || J == hi.y)) // Z edge
616 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);
617 3434659 : else if (I == lo.x || I == hi.x) // X face
618 0 : cdata(I, J, K, n) =
619 0 : (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
620 0 : + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
621 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;
622 3434659 : else if (J == lo.y || J == hi.y) // Y face
623 0 : cdata(I, J, K, n) =
624 0 : (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
625 0 : + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
626 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;
627 3434659 : else if (K == lo.z || K == hi.z) // Z face
628 3079126 : cdata(I, J, K, n) =
629 9237378 : (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
630 9237378 : + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
631 12316504 : + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
632 : else // Interior
633 355533 : cdata(I, J, K, n) =
634 1422132 : (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) +
635 1422132 : 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
636 355533 : +
637 1422132 : (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) +
638 1422132 : 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) +
639 1422132 : 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
640 355533 : +
641 1066599 : (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
642 1066599 : fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
643 355533 : +
644 711066 : fdata(i, j, k, n) / 8.0;
645 : }
646 :
647 3514703 : });
648 : }
649 2180 : }
650 :
651 : // Copy the fine residual restricted onto the coarse grid
652 : // into the final residual.
653 2180 : res.ParallelCopy(fine_res_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
654 :
655 2180 : const int mglev = 0;
656 :
657 : // Sync up ghost nodes
658 2180 : amrex::Geometry geom = m_geom[crse_amrlev][mglev];
659 2180 : realFillBoundary(res, geom);
660 2180 : nodalSync(crse_amrlev, mglev, res);
661 4360 : return;
662 2180 : }
663 :
664 : void
665 17995 : Operator<Grid::Node>::solutionResidual(int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b,
666 : const MultiFab* /*crse_bcdata*/)
667 : {
668 17995 : const int mglev = 0;
669 17995 : const int ncomp = b.nComp();
670 17995 : apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution);
671 17994 : MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, 2);
672 17994 : amrex::Geometry geom = m_geom[amrlev][mglev];
673 17994 : realFillBoundary(resid, geom);
674 17994 : }
675 :
676 : void
677 60247 : Operator<Grid::Node>::correctionResidual(int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b,
678 : BCMode /*bc_mode*/, const MultiFab* /*crse_bcdata*/)
679 : {
680 60247 : resid.setVal(0.0);
681 60247 : apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction);
682 60245 : int ncomp = b.nComp();
683 60245 : MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, resid.nGrow());
684 60245 : amrex::Geometry geom = m_geom[amrlev][mglev];
685 60245 : realFillBoundary(resid, geom);
686 60245 : }
687 :
688 :
689 :
690 :
691 0 : Operator<Grid::Cell>::Operator()
692 : {
693 0 : m_ixtype = amrex::IntVect::TheCellVector();
694 0 : }
695 :
696 : void
697 0 : Operator<Grid::Cell>::define(amrex::Vector<amrex::Geometry> a_geom,
698 : const amrex::Vector<amrex::BoxArray>& a_grids,
699 : const amrex::Vector<amrex::DistributionMapping>& a_dmap,
700 : BC::BC<Set::Scalar>& a_bc,
701 : const amrex::LPInfo& a_info,
702 : const amrex::Vector<amrex::FabFactory<amrex::FArrayBox> const*>& a_factory)
703 : {
704 0 : m_bc = &a_bc;
705 :
706 0 : std::array<int, AMREX_SPACEDIM> is_periodic = m_bc->IsPeriodic();
707 :
708 0 : MLCellLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
709 :
710 0 : Util::Warning(INFO, "This section of code has not been tested.");
711 0 : for (int n = 0; n < getNComp(); n++)
712 : {
713 0 : m_lobc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
714 : is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
715 : is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
716 0 : m_hibc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
717 : is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
718 : is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
719 : }
720 :
721 0 : for (int ilev = 0; ilev < a_geom.size(); ++ilev)
722 0 : setLevelBC(ilev, nullptr);
723 :
724 0 : }
725 :
726 :
727 : void
728 0 : Operator<Grid::Cell>::prepareForSolve()
729 : {
730 0 : MLCellLinOp::prepareForSolve();
731 0 : }
732 :
733 0 : Operator<Grid::Cell>::BndryCondLoc::BndryCondLoc(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm)
734 0 : : bcond(ba, dm),
735 0 : bcloc(ba, dm)
736 : {
737 0 : }
738 :
739 : void
740 0 : Operator<Grid::Cell>::BndryCondLoc::setLOBndryConds(const amrex::Geometry& /*geom*/, const amrex::Real* /*dx*/,
741 : const amrex::Array<BCType, AMREX_SPACEDIM>& /*lobc*/,
742 : const amrex::Array<BCType, AMREX_SPACEDIM>& /*hibc*/,
743 : int /*ratio*/, const amrex::RealVect& /*a_loc*/)
744 : {
745 0 : Util::Warning(INFO, "This code has not been properlyt tested");
746 0 : }
747 :
748 :
749 : void
750 0 : Operator<Grid::Cell>::averageDownCoeffs()
751 : {
752 0 : for (int i = 0; i < m_num_a_fabs; i++)
753 : {
754 0 : for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
755 : {
756 0 : auto& fine_a_coeffs = m_a_coeffs[i][amrlev];
757 0 : averageDownCoeffsSameAmrLevel(fine_a_coeffs);
758 : }
759 0 : averageDownCoeffsSameAmrLevel(m_a_coeffs[i][0]);
760 : }
761 0 : }
762 :
763 : void
764 0 : Operator<Grid::Cell>::averageDownCoeffsSameAmrLevel(amrex::Vector<amrex::MultiFab>& a)
765 : {
766 0 : int nmglevs = a.size();
767 0 : for (int mglev = 1; mglev < nmglevs; ++mglev)
768 : {
769 0 : amrex::average_down(a[mglev - 1], a[mglev], 0, a[0].nComp(), mg_coarsen_ratio);
770 : }
771 0 : }
772 :
773 :
774 :
775 : const amrex::FArrayBox&
776 0 : Operator<Grid::Cell>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
777 : {
778 0 : return m_a_coeffs[num][amrlev][mglev][mfi];
779 : }
780 :
781 :
782 : void
783 0 : Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
784 : {
785 : /// \todo assertions here
786 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
787 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
788 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
789 : {
790 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
791 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
792 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
793 0 : m_dmap[amrlev][mglev],
794 0 : input[amrlev].nComp(),
795 0 : input[amrlev].nGrow());
796 :
797 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
798 0 : input[amrlev], 0, 0,
799 0 : input[amrlev].nComp(),
800 0 : input[amrlev].nGrow());
801 : }
802 0 : m_num_a_fabs++;
803 0 : }
804 :
805 :
806 : void
807 0 : Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
808 : {
809 : /// \todo assertions here
810 0 : m_a_coeffs.resize(m_a_coeffs.size() + 1);
811 0 : m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
812 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
813 : {
814 0 : m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
815 0 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
816 0 : m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
817 0 : m_dmap[amrlev][mglev],
818 0 : input[amrlev]->nComp(),
819 0 : input[amrlev]->nGrow());
820 :
821 0 : amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
822 0 : *input[amrlev], 0, 0,
823 0 : input[amrlev]->nComp(),
824 0 : input[amrlev]->nGrow());
825 : }
826 0 : m_num_a_fabs++;
827 0 : }
828 :
829 :
830 : }
|
