Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLMG.H
Go to the documentation of this file.
1#ifndef AMREX_ML_MG_H_
2#define AMREX_ML_MG_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_MLLinOp.H>
6#include <AMReX_MLCGSolver.H>
7
8namespace amrex {
9
10template <typename MF>
11class MLMGT
12{
13public:
14
15 class error
16 : public std::runtime_error
17 {
18 public :
19 using std::runtime_error::runtime_error;
20 };
21
22 template <typename T> friend class MLCGSolverT;
23 template <typename M> friend class GMRESMLMGT;
24
25 using MFType = MF;
26 using FAB = typename MLLinOpT<MF>::FAB;
27 using RT = typename MLLinOpT<MF>::RT;
28
29 using BCMode = typename MLLinOpT<MF>::BCMode;
31
33 enum class CFStrategy : int {none,ghostnodes};
34
35 MLMGT (MLLinOpT<MF>& a_lp);
37
38 MLMGT (MLMGT<MF> const&) = delete;
39 MLMGT (MLMGT<MF> &&) = delete;
40 MLMGT<MF>& operator= (MLMGT<MF> const&) = delete;
42
43 // Optional argument checkpoint_file is for debugging only.
44 template <typename AMF>
45 RT solve (const Vector<AMF*>& a_sol, const Vector<AMF const*>& a_rhs,
46 RT a_tol_rel, RT a_tol_abs, const char* checkpoint_file = nullptr);
47
48 template <typename AMF>
49 RT solve (std::initializer_list<AMF*> a_sol,
50 std::initializer_list<AMF const*> a_rhs,
51 RT a_tol_rel, RT a_tol_abs, const char* checkpoint_file = nullptr);
52
53 RT precond (Vector<MF*> const& a_sol, Vector<MF const*> const& a_rhs,
54 RT a_tol_rel, RT a_tol_abs);
55
56 template <typename AMF>
57 void getGradSolution (const Vector<Array<AMF*,AMREX_SPACEDIM> >& a_grad_sol,
58 Location a_loc = Location::FaceCenter);
59
60 template <typename AMF>
61 void getGradSolution (std::initializer_list<Array<AMF*,AMREX_SPACEDIM>> a_grad_sol,
62 Location a_loc = Location::FaceCenter);
63
67 template <typename AMF>
68 void getFluxes (const Vector<Array<AMF*,AMREX_SPACEDIM> >& a_flux,
69 Location a_loc = Location::FaceCenter);
70
71 template <typename AMF>
72 void getFluxes (std::initializer_list<Array<AMF*,AMREX_SPACEDIM>> a_flux,
73 Location a_loc = Location::FaceCenter);
74
75 template <typename AMF>
76 void getFluxes (const Vector<Array<AMF*,AMREX_SPACEDIM> >& a_flux,
77 const Vector<AMF*> & a_sol,
78 Location a_loc = Location::FaceCenter);
79
80 template <typename AMF>
81 void getFluxes (std::initializer_list<Array<AMF*,AMREX_SPACEDIM>> a_flux,
82 std::initializer_list<AMF*> a_sol,
83 Location a_loc = Location::FaceCenter);
84
85 template <typename AMF>
86 void getFluxes (const Vector<AMF*> & a_flux,
87 Location a_loc = Location::CellCenter);
88
89 template <typename AMF>
90 void getFluxes (std::initializer_list<AMF*> a_flux,
91 Location a_loc = Location::CellCenter);
92
93 template <typename AMF>
94 void getFluxes (const Vector<AMF*> & a_flux,
95 const Vector<AMF*> & a_sol,
96 Location a_loc = Location::CellCenter);
97
98 template <typename AMF>
99 void getFluxes (std::initializer_list<AMF*> a_flux,
100 std::initializer_list<AMF*> a_sol,
101 Location a_loc = Location::CellCenter);
102
103 void compResidual (const Vector<MF*>& a_res, const Vector<MF*>& a_sol,
104 const Vector<MF const*>& a_rhs);
105
106#ifdef AMREX_USE_EB
107 // Flux into the EB wall
108 void getEBFluxes (const Vector<MF*>& a_eb_flux);
109 void getEBFluxes (const Vector<MF*>& a_eb_flux, const Vector<MF*> & a_sol);
110#endif
111
117 void apply (const Vector<MF*>& out, const Vector<MF*>& in);
118
120 void applyPrecond (const Vector<MF*>& out, const Vector<MF*>& in);
121
122 [[nodiscard]] int getVerbose () const { return verbose; }
123 [[nodiscard]] int getBottomVerbose () const { return bottom_verbose; }
124
125 void incPrintIdentation ();
126 void decPrintIdentation ();
127
128 void setThrowException (bool t) noexcept { throw_exception = t; }
129 void setVerbose (int v) noexcept { verbose = v; }
130 void setMaxIter (int n) noexcept { max_iters = n; }
131 void setMaxFmgIter (int n) noexcept { max_fmg_iters = n; }
132 void setFixedIter (int nit) noexcept { do_fixed_number_of_iters = nit; }
133 void setPrecondIter (int nit) noexcept { max_precond_iters = nit; }
134
135 void setPreSmooth (int n) noexcept { nu1 = n; }
136 void setPostSmooth (int n) noexcept { nu2 = n; }
137 void setFinalSmooth (int n) noexcept { nuf = n; }
138 void setBottomSmooth (int n) noexcept { nub = n; }
139
140 void setBottomSolver (BottomSolver s) noexcept { bottom_solver = s; }
141 [[nodiscard]] BottomSolver getBottomSolver () const noexcept { return bottom_solver; }
142 void setCFStrategy (CFStrategy a_cf_strategy) noexcept {cf_strategy = a_cf_strategy;}
143 void setBottomVerbose (int v) noexcept { bottom_verbose = v; }
144 void setBottomMaxIter (int n) noexcept { bottom_maxiter = n; }
145 void setBottomTolerance (RT t) noexcept { bottom_reltol = t; }
146 void setBottomToleranceAbs (RT t) noexcept { bottom_abstol = t;}
147 [[nodiscard]] RT getBottomToleranceAbs () const noexcept{ return bottom_abstol; }
148
149 void setAlwaysUseBNorm (int flag) noexcept { always_use_bnorm = flag; }
150
151 void setFinalFillBC (int flag) noexcept { final_fill_bc = flag; }
152
153 [[nodiscard]] int numAMRLevels () const noexcept { return namrlevs; }
154
155 void setNSolve (int flag) noexcept { do_nsolve = flag; }
156 void setNSolveGridSize (int s) noexcept { nsolve_grid_size = s; }
157
158#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
159 void setHypreInterface (Hypre::Interface f) noexcept {
160 // must use ij interface for EB
161#ifndef AMREX_USE_EB
162 hypre_interface = f;
163#else
165#endif
166 }
167
169 void setHypreOptionsNamespace(const std::string& prefix) noexcept
170 {
171 hypre_options_namespace = prefix;
172 }
173
174 void setHypreOldDefault (bool l) noexcept {hypre_old_default = l;}
175 void setHypreRelaxType (int n) noexcept {hypre_relax_type = n;}
176 void setHypreRelaxOrder (int n) noexcept {hypre_relax_order = n;}
177 void setHypreNumSweeps (int n) noexcept {hypre_num_sweeps = n;}
178 void setHypreStrongThreshold (Real t) noexcept {hypre_strong_threshold = t;}
179#endif
180
181 void prepareForFluxes (Vector<MF const*> const& a_sol);
182
183 template <typename AMF>
184 void prepareForSolve (Vector<AMF*> const& a_sol, Vector<AMF const*> const& a_rhs);
185
186 void prepareForNSolve ();
187
188 void prepareLinOp ();
189
190 void preparePrecond ();
191
192 void oneIter (int iter);
193
194 void miniCycle (int amrlev);
195
196 void mgVcycle (int amrlev, int mglev);
197 void mgFcycle ();
198
199 void bottomSolve ();
200 void NSolve (MLMGT<MF>& a_solver, MF& a_sol, MF& a_rhs);
201 void actualBottomSolve ();
202
203 void computeMLResidual (int amrlevmax);
204 void computeResidual (int alev);
205 void computeResWithCrseSolFineCor (int calev, int falev);
206 void computeResWithCrseCorFineCor (int falev);
207 void interpCorrection (int alev);
208 void interpCorrection (int alev, int mglev);
209 void addInterpCorrection (int alev, int mglev);
210
211 void computeResOfCorrection (int amrlev, int mglev);
212
213 RT ResNormInf (int alev, bool local = false);
214 RT MLResNormInf (int alevmax, bool local = false);
215 RT MLRhsNormInf (bool local = false);
216
217 void makeSolvable ();
218 void makeSolvable (int amrlev, int mglev, MF& mf);
219
220#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
221 template <class TMF=MF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,int> = 0>
222 void bottomSolveWithHypre (MF& x, const MF& b);
223#endif
224
225#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
226 template <class TMF=MF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,int> = 0>
227 void bottomSolveWithPETSc (MF& x, const MF& b);
228#endif
229
230 int bottomSolveWithCG (MF& x, const MF& b, typename MLCGSolverT<MF>::Type type);
231
232 [[nodiscard]] RT getInitRHS () const noexcept { return m_rhsnorm0; }
233 // Initial composite residual
234 [[nodiscard]] RT getInitResidual () const noexcept { return m_init_resnorm0; }
235 // Final composite residual
236 [[nodiscard]] RT getFinalResidual () const noexcept { return m_final_resnorm0; }
237 // Residuals on the *finest* AMR level after each iteration
238 [[nodiscard]] Vector<RT> const& getResidualHistory () const noexcept { return m_iter_fine_resnorm0; }
239 [[nodiscard]] int getNumIters () const noexcept { return m_iter_fine_resnorm0.size(); }
240 [[nodiscard]] Vector<int> const& getNumCGIters () const noexcept { return m_niters_cg; }
241
243
244private:
245
246 bool precond_mode = false;
247 bool throw_exception = false;
248 int verbose = 1;
249
250 int max_iters = 200;
253
254 int nu1 = 2;
255 int nu2 = 2;
256 int nuf = 8;
257 int nub = 0;
258
260
264 int bottom_maxiter = 200;
265 RT bottom_reltol = std::is_same<RT,double>() ? RT(1.e-4) : RT(1.e-3);
267
269
271
273 int ncomp;
276
277 bool linop_prepared = false;
278 Long solve_called = 0;
279
281 int do_nsolve = false;
283 std::unique_ptr<MLLinOpT<MF>> ns_linop;
284 std::unique_ptr<MLMGT<MF>> ns_mlmg;
285 std::unique_ptr<MF> ns_sol;
286 std::unique_ptr<MF> ns_rhs;
287
288 std::string print_ident;
289
291#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
292 // Hypre::Interface hypre_interface = Hypre::Interface::structed;
293 // Hypre::Interface hypre_interface = Hypre::Interface::semi_structed;
294 Hypre::Interface hypre_interface = Hypre::Interface::ij;
295
296 std::unique_ptr<Hypre> hypre_solver;
297 std::unique_ptr<MLMGBndryT<MF>> hypre_bndry;
298 std::unique_ptr<HypreNodeLap> hypre_node_solver;
299
300 std::string hypre_options_namespace = "hypre";
301 bool hypre_old_default = true; // Falgout coarsening with modified classical interpolation
302 int hypre_relax_type = 6; // G-S/Jacobi hybrid relaxation
303 int hypre_relax_order = 1; // uses C/F relaxation
304 int hypre_num_sweeps = 2; // Sweeps on each level
305 Real hypre_strong_threshold = 0.25; // Hypre default is 0.25
306#endif
307
309#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
310 std::unique_ptr<PETScABecLap> petsc_solver;
311 std::unique_ptr<MLMGBndryT<MF>> petsc_bndry;
312#endif
313
321
323
333
336
337 RT m_rhsnorm0 = RT(-1.0);
341 Vector<RT> m_iter_fine_resnorm0; // Residual for each iteration at the finest level
342
343 void checkPoint (const Vector<MultiFab*>& a_sol,
344 const Vector<MultiFab const*>& a_rhs,
345 RT a_tol_rel, RT a_tol_abs, const char* a_file_name) const;
346
347};
348
349template <typename MF>
351 : linop(a_lp), ncomp(a_lp.getNComp()), namrlevs(a_lp.NAMRLevels()),
352 finest_amr_lev(a_lp.NAMRLevels()-1)
353{}
354
355template <typename MF> MLMGT<MF>::~MLMGT () = default;
356
357template <typename MF>
358template <typename AMF>
359auto
360MLMGT<MF>::solve (std::initializer_list<AMF*> a_sol,
361 std::initializer_list<AMF const*> a_rhs,
362 RT a_tol_rel, RT a_tol_abs, const char* checkpoint_file) -> RT
363{
364 return solve(Vector<AMF*>(std::move(a_sol)),
365 Vector<AMF const*>(std::move(a_rhs)),
366 a_tol_rel, a_tol_abs, checkpoint_file);
367}
368
369template <typename MF>
370template <typename AMF>
371auto
373 RT a_tol_rel, RT a_tol_abs, const char* checkpoint_file) -> RT
374{
375 BL_PROFILE("MLMG::solve()");
376
377 if constexpr (std::is_same<AMF,MultiFab>()) {
378 if (checkpoint_file != nullptr) {
379 checkPoint(a_sol, a_rhs, a_tol_rel, a_tol_abs, checkpoint_file);
380 }
381 }
382
383 if (bottom_solver == BottomSolver::Default) {
384 bottom_solver = linop.getDefaultBottomSolver();
385 }
386
387#if (defined(AMREX_USE_HYPRE) || defined(AMREX_USE_PETSC)) && (AMREX_SPACEDIM > 1)
388 if (bottom_solver == BottomSolver::hypre || bottom_solver == BottomSolver::petsc) {
389 int mo = linop.getMaxOrder();
390 if (a_sol[0]->hasEBFabFactory()) {
391 linop.setMaxOrder(2);
392 } else {
393 linop.setMaxOrder(std::min(3,mo)); // maxorder = 4 not supported
394 }
395 }
396#endif
397
398 bool is_nsolve = linop.m_parent;
399
400 auto solve_start_time = amrex::second();
401
402 RT& composite_norminf = m_final_resnorm0;
403
404 m_niters_cg.clear();
405 m_iter_fine_resnorm0.clear();
406
407 prepareForSolve(a_sol, a_rhs);
408
409 computeMLResidual(finest_amr_lev);
410
411 bool local = true;
412 RT resnorm0 = MLResNormInf(finest_amr_lev, local);
413 RT rhsnorm0 = MLRhsNormInf(local);
414 if (!is_nsolve) {
415 ParallelAllReduce::Max<RT>({resnorm0, rhsnorm0}, ParallelContext::CommunicatorSub());
416
417 if (verbose >= 1)
418 {
419 amrex::Print() << print_ident << "MLMG: Initial rhs = " << rhsnorm0 << "\n"
420 << print_ident << "MLMG: Initial residual (resid0) = " << resnorm0 << "\n";
421 }
422 }
423
424 m_init_resnorm0 = resnorm0;
425 m_rhsnorm0 = rhsnorm0;
426
427 RT max_norm;
428 std::string norm_name;
429 if (always_use_bnorm || rhsnorm0 >= resnorm0) {
430 norm_name = "bnorm";
431 max_norm = rhsnorm0;
432 } else {
433 norm_name = "resid0";
434 max_norm = resnorm0;
435 }
436 const RT res_target = std::max(a_tol_abs, std::max(a_tol_rel,RT(1.e-16))*max_norm);
437
438 if (!is_nsolve && resnorm0 <= res_target) {
439 composite_norminf = resnorm0;
440 if (verbose >= 1) {
441 amrex::Print() << print_ident << "MLMG: No iterations needed\n";
442 }
443 } else {
444 auto iter_start_time = amrex::second();
445 bool converged = false;
446
447 const int niters = do_fixed_number_of_iters ? do_fixed_number_of_iters : max_iters;
448 for (int iter = 0; iter < niters; ++iter)
449 {
450 oneIter(iter);
451
452 converged = false;
453
454 // Test convergence on the fine amr level
455 computeResidual(finest_amr_lev);
456
457 if (is_nsolve) { continue; }
458
459 RT fine_norminf = ResNormInf(finest_amr_lev);
460 m_iter_fine_resnorm0.push_back(fine_norminf);
461 composite_norminf = fine_norminf;
462 if (verbose >= 2) {
463 amrex::Print() << print_ident << "MLMG: Iteration " << std::setw(3) << iter+1 << " Fine resid/"
464 << norm_name << " = " << fine_norminf/max_norm << "\n";
465 }
466 bool fine_converged = (fine_norminf <= res_target);
467
468 if (namrlevs == 1 && fine_converged) {
469 converged = true;
470 } else if (fine_converged) {
471 // finest level is converged, but we still need to test the coarse levels
472 computeMLResidual(finest_amr_lev-1);
473 RT crse_norminf = MLResNormInf(finest_amr_lev-1);
474 if (verbose >= 2) {
475 amrex::Print() << print_ident << "MLMG: Iteration " << std::setw(3) << iter+1
476 << " Crse resid/" << norm_name << " = "
477 << crse_norminf/max_norm << "\n";
478 }
479 converged = (crse_norminf <= res_target);
480 composite_norminf = std::max(fine_norminf, crse_norminf);
481 } else {
482 converged = false;
483 }
484
485 if (converged) {
486 if (verbose >= 1) {
487 amrex::Print() << print_ident << "MLMG: Final Iter. " << iter+1
488 << " resid, resid/" << norm_name << " = "
489 << composite_norminf << ", "
490 << composite_norminf/max_norm << "\n";
491 }
492 break;
493 } else {
494 if (composite_norminf > RT(1.e20)*max_norm)
495 {
496 if (verbose > 0) {
497 amrex::Print() << print_ident << "MLMG: Failing to converge after " << iter+1 << " iterations."
498 << " resid, resid/" << norm_name << " = "
499 << composite_norminf << ", "
500 << composite_norminf/max_norm << "\n";
501 }
502
503 if ( throw_exception ) {
504 throw error("MLMG blew up.");
505 } else {
506 amrex::Abort("MLMG failing so lets stop here");
507 }
508 }
509 }
510 }
511
512 if (!converged && do_fixed_number_of_iters == 0) {
513 if (verbose > 0) {
514 amrex::Print() << print_ident << "MLMG: Failed to converge after " << max_iters << " iterations."
515 << " resid, resid/" << norm_name << " = "
516 << composite_norminf << ", "
517 << composite_norminf/max_norm << "\n";
518 }
519
520 if ( throw_exception ) {
521 throw error("MLMG failed to converge.");
522 } else {
523 amrex::Abort("MLMG failed.");
524 }
525 }
526 timer[iter_time] = amrex::second() - iter_start_time;
527 }
528
529 linop.postSolve(GetVecOfPtrs(sol));
530
531 IntVect ng_back = final_fill_bc ? IntVect(1) : IntVect(0);
532 if (linop.hasHiddenDimension()) {
533 ng_back[linop.hiddenDirection()] = 0;
534 }
535 for (int alev = 0; alev < namrlevs; ++alev)
536 {
537 if (!sol_is_alias[alev]) {
538 LocalCopy(*a_sol[alev], sol[alev], 0, 0, ncomp, ng_back);
539 }
540 }
541
542 timer[solve_time] = amrex::second() - solve_start_time;
543 if (verbose >= 1) {
544 ParallelReduce::Max<double>(timer.data(), timer.size(), 0,
547 {
548 amrex::AllPrint() << print_ident << "MLMG: Timers: Solve = " << timer[solve_time]
549 << " Iter = " << timer[iter_time]
550 << " Bottom = " << timer[bottom_time] << "\n";
551 }
552 }
553
554 ++solve_called;
555
556 return composite_norminf;
557}
558
559template <typename MF>
560auto
562 RT a_tol_rel, RT a_tol_abs) -> RT
563{
564 precond_mode = true;
565 std::swap(max_precond_iters, do_fixed_number_of_iters);
566 linop.beginPrecondBC();
567
568 auto r = solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
569
570 linop.endPrecondBC();
571 std::swap(max_precond_iters, do_fixed_number_of_iters);
572 precond_mode = false;
573
574 return r;
575}
576
577template <typename MF>
578void
580{
581 for (int alev = finest_amr_lev; alev >= 0; --alev) {
582 const MF* crse_bcdata = (alev > 0) ? a_sol[alev-1] : nullptr;
583 linop.prepareForFluxes(alev, crse_bcdata);
584 }
585}
586
587template <typename MF>
588template <typename AMF>
589void
591{
592 BL_PROFILE("MLMG::getGradSolution()");
593 for (int alev = 0; alev <= finest_amr_lev; ++alev) {
594 if constexpr (std::is_same<AMF,MF>()) {
595 linop.compGrad(alev, a_grad_sol[alev], sol[alev], a_loc);
596 } else {
598 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
599 auto const& amf = *(a_grad_sol[alev][idim]);
600 grad_sol[idim].define(boxArray(amf), DistributionMap(amf), ncomp, 0);
601 }
602 linop.compGrad(alev, GetArrOfPtrs(grad_sol), sol[alev], a_loc);
603 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
604 LocalCopy(*a_grad_sol[alev][idim], grad_sol[idim], 0, 0, ncomp, IntVect(0));
605 }
606 }
607 }
608}
609
610template <typename MF>
611template <typename AMF>
612void
613MLMGT<MF>::getGradSolution (std::initializer_list<Array<AMF*,AMREX_SPACEDIM>> a_grad_sol, Location a_loc)
614{
615 getGradSolution(Vector<Array<AMF*,AMREX_SPACEDIM>>(std::move(a_grad_sol)), a_loc);
616}
617
618template <typename MF>
619template <typename AMF>
620void
622 Location a_loc)
623{
624 if (!linop.isCellCentered()) {
625 amrex::Abort("Calling wrong getFluxes for nodal solver");
626 }
627
628 AMREX_ASSERT(sol.size() == a_flux.size());
629
630 if constexpr (std::is_same<AMF,MF>()) {
631 getFluxes(a_flux, GetVecOfPtrs(sol), a_loc);
632 } else {
633 Vector<Array<MF,AMREX_SPACEDIM>> fluxes(namrlevs);
634 for (int ilev = 0; ilev < namrlevs; ++ilev) {
635 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
636 auto const& amf = *(a_flux[ilev][idim]);
637 fluxes[ilev][idim].define(boxArray(amf), DistributionMap(amf), ncomp, 0);
638 }
639 }
640 getFluxes(GetVecOfArrOfPtrs(fluxes), GetVecOfPtrs(sol), a_loc);
641 for (int ilev = 0; ilev < namrlevs; ++ilev) {
642 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
643 LocalCopy(*a_flux[ilev][idim], fluxes[ilev][idim], 0, 0, ncomp, IntVect(0));
644 }
645 }
646 }
647}
648
649template <typename MF>
650template <typename AMF>
651void
653 Location a_loc)
654{
655 getFluxes(Vector<Array<AMF*,AMREX_SPACEDIM>>(std::move(a_flux)), a_loc);
656}
657
658template <typename MF>
659template <typename AMF>
660void
662 const Vector<AMF*>& a_sol, Location a_loc)
663{
664 BL_PROFILE("MLMG::getFluxes()");
665
666 if (!linop.isCellCentered()) {
667 amrex::Abort("Calling wrong getFluxes for nodal solver");
668 }
669
670 if constexpr (std::is_same<AMF,MF>()) {
671 linop.getFluxes(a_flux, a_sol, a_loc);
672 } else {
673 Vector<Array<MF,AMREX_SPACEDIM>> fluxes(namrlevs);
674 for (int ilev = 0; ilev < namrlevs; ++ilev) {
675 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
676 auto const& amf = *(a_flux[ilev][idim]);
677 fluxes[ilev][idim].define(boxArray(amf), DistributionMap(amf), ncomp, 0);
678 }
679 LocalCopy(sol[ilev], *a_sol[ilev], 0, 0, ncomp, nGrowVect(sol[ilev]));
680 }
681 linop.getFluxes(GetVecOfArrOfPtrs(fluxes), GetVecOfPtrs(sol), a_loc);
682 for (int ilev = 0; ilev < namrlevs; ++ilev) {
683 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
684 LocalCopy(*a_flux[ilev][idim], fluxes[ilev][idim], 0, 0, ncomp, IntVect(0));
685 }
686 }
687 }
688}
689
690template <typename MF>
691template <typename AMF>
692void
694 std::initializer_list<AMF*> a_sol, Location a_loc)
695{
696 getFluxes(Vector<Array<AMF*,AMREX_SPACEDIM>>(std::move(a_flux)),
697 Vector<AMF*>(std::move(a_sol)), a_loc);
698}
699
700template <typename MF>
701template <typename AMF>
702void
704{
705 AMREX_ASSERT(sol.size() == a_flux.size());
706 if constexpr (std::is_same<AMF,MF>()) {
707 getFluxes(a_flux, GetVecOfPtrs(sol), a_loc);
708 } else {
709 Vector<MF> fluxes(namrlevs);
710 for (int ilev = 0; ilev < namrlevs; ++ilev) {
711 auto const& amf = *a_flux[ilev];
712 fluxes[ilev].define(boxArray(amf), DistributionMap(amf), ncomp, 0);
713 }
714 getFluxes(GetVecOfPtrs(fluxes), GetVecOfPtrs(sol), a_loc);
715 for (int ilev = 0; ilev < namrlevs; ++ilev) {
716 LocalCopy(*a_flux[ilev], fluxes[ilev], 0, 0, ncomp, IntVect(0));
717 }
718 }
719}
720
721template <typename MF>
722template <typename AMF>
723void
724MLMGT<MF>::getFluxes (std::initializer_list<AMF*> a_flux, Location a_loc)
725{
726 getFluxes(Vector<AMF*>(std::move(a_flux)), a_loc);
727}
728
729template <typename MF>
730template <typename AMF>
731void
733 const Vector<AMF*>& a_sol, Location /*a_loc*/)
734{
735 AMREX_ASSERT(nComp(*a_flux[0]) >= AMREX_SPACEDIM);
736
737 if constexpr (! std::is_same<AMF,MF>()) {
738 for (int alev = 0; alev < namrlevs; ++alev) {
739 LocalCopy(sol[alev], *a_sol[alev], 0, 0, ncomp, nGrowVect(sol[alev]));
740 }
741 }
742
743 if (linop.isCellCentered())
744 {
745 Vector<Array<MF,AMREX_SPACEDIM> > ffluxes(namrlevs);
746 for (int alev = 0; alev < namrlevs; ++alev) {
747 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
748 const int mglev = 0;
749 int nghost = 0;
750 if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(alev); }
751 ffluxes[alev][idim].define(amrex::convert(linop.m_grids[alev][mglev],
753 linop.m_dmap[alev][mglev], ncomp, nghost, MFInfo(),
754 *linop.m_factory[alev][mglev]);
755 }
756 }
757 if constexpr (std::is_same<AMF,MF>()) {
758 getFluxes(amrex::GetVecOfArrOfPtrs(ffluxes), a_sol, Location::FaceCenter);
759 } else {
760 getFluxes(amrex::GetVecOfArrOfPtrs(ffluxes), GetVecOfPtrs(sol), Location::FaceCenter);
761 }
762 for (int alev = 0; alev < namrlevs; ++alev) {
763#ifdef AMREX_USE_EB
764 EB_average_face_to_cellcenter(*a_flux[alev], 0, amrex::GetArrOfConstPtrs(ffluxes[alev]));
765#else
766 average_face_to_cellcenter(*a_flux[alev], 0, amrex::GetArrOfConstPtrs(ffluxes[alev]));
767#endif
768 }
769
770 } else {
771 if constexpr (std::is_same<AMF,MF>()) {
772 linop.getFluxes(a_flux, a_sol);
773 } else {
774 Vector<MF> fluxes(namrlevs);
775 for (int ilev = 0; ilev < namrlevs; ++ilev) {
776 auto const& amf = *a_flux[ilev];
777 fluxes[ilev].define(boxArray(amf), DistributionMap(amf), ncomp, 0);
778 }
779 linop.getFluxes(GetVecOfPtrs(fluxes), GetVecOfPtrs(sol));
780 for (int ilev = 0; ilev < namrlevs; ++ilev) {
781 LocalCopy(*a_flux[ilev], fluxes[ilev], 0, 0, ncomp, IntVect(0));
782 }
783 }
784 }
785}
786
787template <typename MF>
788template <typename AMF>
789void
790MLMGT<MF>::getFluxes (std::initializer_list<AMF*> a_flux,
791 std::initializer_list<AMF*> a_sol, Location a_loc)
792{
793 getFluxes(Vector<AMF*>(std::move(a_flux)),
794 Vector<AMF*>(std::move(a_sol)), a_loc);
795}
796
797#ifdef AMREX_USE_EB
798template <typename MF>
799void
800MLMGT<MF>::getEBFluxes (const Vector<MF*>& a_eb_flux)
801{
802 if (!linop.isCellCentered()) {
803 amrex::Abort("getEBFluxes is for cell-centered only");
804 }
805
806 AMREX_ASSERT(sol.size() == a_eb_flux.size());
807 getEBFluxes(a_eb_flux, GetVecOfPtrs(sol));
808}
809
810template <typename MF>
811void
812MLMGT<MF>::getEBFluxes (const Vector<MF*>& a_eb_flux, const Vector<MF*>& a_sol)
813{
814 BL_PROFILE("MLMG::getEBFluxes()");
815
816 if (!linop.isCellCentered()) {
817 amrex::Abort("getEBFluxes is for cell-centered only");
818 }
819
820 linop.getEBFluxes(a_eb_flux, a_sol);
821}
822#endif
823
824template <typename MF>
825void
827 const Vector<MF const*>& a_rhs)
828{
829 BL_PROFILE("MLMG::compResidual()");
830
831 IntVect ng_sol(1);
832 if (linop.hasHiddenDimension()) { ng_sol[linop.hiddenDirection()] = 0; }
833
834 sol.resize(namrlevs);
835 sol_is_alias.resize(namrlevs,true);
836 for (int alev = 0; alev < namrlevs; ++alev)
837 {
838 if (cf_strategy == CFStrategy::ghostnodes || nGrowVect(*a_sol[alev]) == ng_sol)
839 {
840 sol[alev] = linop.makeAlias(*a_sol[alev]);
841 sol_is_alias[alev] = true;
842 }
843 else
844 {
845 if (sol_is_alias[alev])
846 {
847 sol[alev] = linop.make(alev, 0, ng_sol);
848 }
849 LocalCopy(sol[alev], *a_sol[alev], 0, 0, ncomp, IntVect(0));
850 }
851 }
852
853 prepareLinOp();
854
855 const auto& amrrr = linop.AMRRefRatio();
856
857 for (int alev = finest_amr_lev; alev >= 0; --alev) {
858 const MF* crse_bcdata = (alev > 0) ? &(sol[alev-1]) : nullptr;
859 const MF* prhs = a_rhs[alev];
860#if (AMREX_SPACEDIM != 3)
861 int nghost = (cf_strategy == CFStrategy::ghostnodes) ? linop.getNGrow(alev) : 0;
862 MF rhstmp(boxArray(*prhs), DistributionMap(*prhs), ncomp, nghost,
863 MFInfo(), *linop.Factory(alev));
864 LocalCopy(rhstmp, *prhs, 0, 0, ncomp, IntVect(nghost));
865 linop.applyMetricTerm(alev, 0, rhstmp);
866 linop.unimposeNeumannBC(alev, rhstmp);
867 linop.applyInhomogNeumannTerm(alev, rhstmp);
868 prhs = &rhstmp;
869#endif
870 linop.solutionResidual(alev, *a_res[alev], sol[alev], *prhs, crse_bcdata);
871 if (alev < finest_amr_lev) {
872 linop.reflux(alev, *a_res[alev], sol[alev], *prhs,
873 *a_res[alev+1], sol[alev+1], *a_rhs[alev+1]);
874 if (linop.isCellCentered()) {
875#ifdef AMREX_USE_EB
876 EB_average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]);
877#else
878 average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]);
879#endif
880 }
881 }
882 }
883
884
885#if (AMREX_SPACEDIM != 3)
886 for (int alev = 0; alev <= finest_amr_lev; ++alev) {
887 linop.unapplyMetricTerm(alev, 0, *a_res[alev]);
888 }
889#endif
890}
891
892template <typename MF>
893void
894MLMGT<MF>::apply (const Vector<MF*>& out, const Vector<MF*>& a_in)
895{
896 BL_PROFILE("MLMG::apply()");
897
898 Vector<MF*> in(namrlevs);
899 Vector<MF> in_raii(namrlevs);
900 Vector<MF> rh(namrlevs);
901 int nghost = 0;
902 IntVect ng_sol(1);
903 if (linop.hasHiddenDimension()) { ng_sol[linop.hiddenDirection()] = 0; }
904
905 for (int alev = 0; alev < namrlevs; ++alev)
906 {
907 if (cf_strategy == CFStrategy::ghostnodes)
908 {
909 nghost = linop.getNGrow(alev);
910 in[alev] = a_in[alev];
911 }
912 else if (nGrowVect(*a_in[alev]) == ng_sol)
913 {
914 in[alev] = a_in[alev];
915 }
916 else
917 {
918 IntVect ng = ng_sol;
919 if (cf_strategy == CFStrategy::ghostnodes) { ng = IntVect(nghost); }
920 in_raii[alev] = linop.make(alev, 0, ng);
921 LocalCopy(in_raii[alev], *a_in[alev], 0, 0, ncomp, IntVect(nghost));
922 in[alev] = &(in_raii[alev]);
923 }
924 rh[alev] = linop.make(alev, 0, IntVect(nghost));
925 setVal(rh[alev], RT(0.0));
926 }
927
928 prepareLinOp();
929
930 for (int alev = 0; alev < namrlevs; ++alev) {
931 linop.applyInhomogNeumannTerm(alev, rh[alev]);
932 }
933
934 const auto& amrrr = linop.AMRRefRatio();
935
936 for (int alev = finest_amr_lev; alev >= 0; --alev) {
937 const MF* crse_bcdata = (alev > 0) ? in[alev-1] : nullptr;
938 linop.solutionResidual(alev, *out[alev], *in[alev], rh[alev], crse_bcdata);
939 if (alev < finest_amr_lev) {
940 linop.reflux(alev, *out[alev], *in[alev], rh[alev],
941 *out[alev+1], *in[alev+1], rh[alev+1]);
942 if (linop.isCellCentered()) {
943 if constexpr (IsMultiFabLike_v<MF>) {
944#ifdef AMREX_USE_EB
945 EB_average_down(*out[alev+1], *out[alev], 0, nComp(*out[alev]), amrrr[alev]);
946#else
947 average_down(*out[alev+1], *out[alev], 0, nComp(*out[alev]), amrrr[alev]);
948#endif
949 } else {
950 amrex::Abort("MLMG: TODO average_down for non-MultiFab");
951 }
952 }
953 }
954 }
955
956#if (AMREX_SPACEDIM != 3)
957 for (int alev = 0; alev <= finest_amr_lev; ++alev) {
958 linop.unapplyMetricTerm(alev, 0, *out[alev]);
959 }
960#endif
961
962 for (int alev = 0; alev <= finest_amr_lev; ++alev) {
963 if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(alev); }
964 Scale(*out[alev], RT(-1), 0, nComp(*out[alev]), nghost);
965 }
966}
967
968template <typename MF>
969void
971{
972 precond_mode = true;
973 linop.beginPrecondBC();
974 apply(out, in);
975 linop.endPrecondBC();
976 precond_mode = false;
977}
978
979template <typename MF>
980template <typename AMF>
981void
983{
984 BL_PROFILE("MLMG::prepareForSolve()");
985
986 AMREX_ASSERT(namrlevs <= a_sol.size());
987 AMREX_ASSERT(namrlevs <= a_rhs.size());
988
989 timer.assign(ntimers, 0.0);
990
991 IntVect ng_rhs(0);
992 IntVect ng_sol(1);
993 if (linop.hasHiddenDimension()) { ng_sol[linop.hiddenDirection()] = 0; }
994
995 if (!linop_prepared) {
996 linop.prepareForSolve();
997 linop_prepared = true;
998 } else if (linop.needsUpdate()) {
999 linop.update();
1000
1001#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
1002 hypre_solver.reset();
1003 hypre_bndry.reset();
1004 hypre_node_solver.reset();
1005#endif
1006
1007#ifdef AMREX_USE_PETSC
1008 petsc_solver.reset();
1009 petsc_bndry.reset();
1010#endif
1011 }
1012
1013 sol.resize(namrlevs);
1014 sol_is_alias.resize(namrlevs,false);
1015 for (int alev = 0; alev < namrlevs; ++alev)
1016 {
1017 if (cf_strategy == CFStrategy::ghostnodes)
1018 {
1019 if constexpr (std::is_same<AMF,MF>()) {
1020 sol[alev] = linop.makeAlias(*a_sol[alev]);
1021 sol_is_alias[alev] = true;
1022 } else {
1023 amrex::Abort("Type conversion not supported for CFStrategy::ghostnodes");
1024 }
1025 }
1026 else
1027 {
1028 if (nGrowVect(*a_sol[alev]) == ng_sol) {
1029 if constexpr (std::is_same<AMF,MF>()) {
1030 sol[alev] = linop.makeAlias(*a_sol[alev]);
1031 sol_is_alias[alev] = true;
1032 }
1033 }
1034 if (!sol_is_alias[alev]) {
1035 if (!solve_called) {
1036 sol[alev] = linop.make(alev, 0, ng_sol);
1037 }
1038 LocalCopy(sol[alev], *a_sol[alev], 0, 0, ncomp, IntVect(0));
1039 setBndry(sol[alev], RT(0.0), 0, ncomp);
1040 }
1041 }
1042 }
1043
1044 rhs.resize(namrlevs);
1045 for (int alev = 0; alev < namrlevs; ++alev)
1046 {
1047 if (cf_strategy == CFStrategy::ghostnodes) { ng_rhs = IntVect(linop.getNGrow(alev)); }
1048 if (!solve_called) {
1049 rhs[alev] = linop.make(alev, 0, ng_rhs);
1050 }
1051 LocalCopy(rhs[alev], *a_rhs[alev], 0, 0, ncomp, ng_rhs);
1052 linop.applyMetricTerm(alev, 0, rhs[alev]);
1053 linop.unimposeNeumannBC(alev, rhs[alev]);
1054 linop.applyInhomogNeumannTerm(alev, rhs[alev]);
1055 linop.applyOverset(alev, rhs[alev]);
1056 if ( ! precond_mode) {
1057 bool r = linop.scaleRHS(alev, &(rhs[alev]));
1059 }
1060
1061#ifdef AMREX_USE_EB
1062 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(linop.Factory(alev));
1063 if (factory && !factory->isAllRegular()) {
1064 if constexpr (std::is_same<MF,MultiFab>()) {
1065 EB_set_covered(rhs[alev], 0, ncomp, 0, RT(0.0));
1066 EB_set_covered(sol[alev], 0, ncomp, 0, RT(0.0));
1067 } else {
1068 amrex::Abort("TODO: MLMG with EB only works with MultiFab");
1069 }
1070 }
1071#endif
1072 }
1073
1074 for (int falev = finest_amr_lev; falev > 0; --falev)
1075 {
1076 linop.averageDownSolutionRHS(falev-1, sol[falev-1], rhs[falev-1], sol[falev], rhs[falev]);
1077 }
1078
1079 // enforce solvability if appropriate
1080 if (linop.isSingular(0) && linop.getEnforceSingularSolvable())
1081 {
1082 makeSolvable();
1083 }
1084
1085 IntVect ng = linop.getNGrowVectRestriction();
1086 if (cf_strategy == CFStrategy::ghostnodes) { ng = ng_rhs; }
1087 if (!solve_called) {
1088 linop.make(res, ng);
1089 linop.make(rescor, ng);
1090 }
1091 for (int alev = 0; alev <= finest_amr_lev; ++alev)
1092 {
1093 const int nmglevs = linop.NMGLevels(alev);
1094 for (int mglev = 0; mglev < nmglevs; ++mglev)
1095 {
1096 setVal(res [alev][mglev], RT(0.0));
1097 setVal(rescor[alev][mglev], RT(0.0));
1098 }
1099 }
1100
1101 if (cf_strategy != CFStrategy::ghostnodes) { ng = ng_sol; }
1102 cor.resize(namrlevs);
1103 for (int alev = 0; alev <= finest_amr_lev; ++alev)
1104 {
1105 const int nmglevs = linop.NMGLevels(alev);
1106 cor[alev].resize(nmglevs);
1107 for (int mglev = 0; mglev < nmglevs; ++mglev)
1108 {
1109 if (!solve_called) {
1110 IntVect _ng = ng;
1111 if (cf_strategy == CFStrategy::ghostnodes) { _ng=IntVect(linop.getNGrow(alev,mglev)); }
1112 cor[alev][mglev] = linop.make(alev, mglev, _ng);
1113 }
1114 setVal(cor[alev][mglev], RT(0.0));
1115 }
1116 }
1117
1118 cor_hold.resize(std::max(namrlevs-1,1));
1119 {
1120 const int alev = 0;
1121 const int nmglevs = linop.NMGLevels(alev);
1122 cor_hold[alev].resize(nmglevs);
1123 for (int mglev = 0; mglev < nmglevs-1; ++mglev)
1124 {
1125 if (!solve_called) {
1126 IntVect _ng = ng;
1127 if (cf_strategy == CFStrategy::ghostnodes) { _ng=IntVect(linop.getNGrow(alev,mglev)); }
1128 cor_hold[alev][mglev] = linop.make(alev, mglev, _ng);
1129 }
1130 setVal(cor_hold[alev][mglev], RT(0.0));
1131 }
1132 }
1133 for (int alev = 1; alev < finest_amr_lev; ++alev)
1134 {
1135 cor_hold[alev].resize(1);
1136 if (!solve_called) {
1137 IntVect _ng = ng;
1138 if (cf_strategy == CFStrategy::ghostnodes) { _ng=IntVect(linop.getNGrow(alev)); }
1139 cor_hold[alev][0] = linop.make(alev, 0, _ng);
1140 }
1141 setVal(cor_hold[alev][0], RT(0.0));
1142 }
1143
1144 if (linop.m_parent // no embedded N-Solve
1145 || !linop.supportNSolve())
1146 {
1147 do_nsolve = false;
1148 }
1149
1150 if (do_nsolve && ns_linop == nullptr)
1151 {
1152 prepareForNSolve();
1153 }
1154
1155 if (verbose >= 2) {
1156 amrex::Print() << print_ident << "MLMG: # of AMR levels: " << namrlevs << "\n"
1157 << print_ident << " # of MG levels on the coarsest AMR level: " << linop.NMGLevels(0)
1158 << "\n";
1159 if (ns_linop) {
1160 amrex::Print() << print_ident << " # of MG levels in N-Solve: " << ns_linop->NMGLevels(0) << "\n"
1161 << print_ident << " # of grids in N-Solve: " << ns_linop->m_grids[0][0].size() << "\n";
1162 }
1163 }
1164}
1165
1166template <typename MF>
1167void
1169{
1170 if (!linop_prepared) {
1171 linop.prepareForSolve();
1172 linop_prepared = true;
1173 } else if (linop.needsUpdate()) {
1174 linop.update();
1175 }
1176}
1177
1178template <typename MF>
1179void
1181{
1182 prepareLinOp();
1183 linop.preparePrecond();
1184}
1185
1186template <typename MF>
1187void
1189{
1190 if constexpr (IsMultiFabLike_v<MF>) {
1191 ns_linop = linop.makeNLinOp(nsolve_grid_size);
1192
1193 int nghost = 0;
1194 if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(); }
1195
1196 const BoxArray& ba = (*ns_linop).m_grids[0][0];
1197 const DistributionMapping& dm =(*ns_linop).m_dmap[0][0];
1198
1199 int ng = 1;
1200 if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; }
1201 ns_sol = std::make_unique<MF>(ba, dm, ncomp, ng, MFInfo(), *(ns_linop->Factory(0,0)));
1202 ng = 0;
1203 if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; }
1204 ns_rhs = std::make_unique<MF>(ba, dm, ncomp, ng, MFInfo(), *(ns_linop->Factory(0,0)));
1205 setVal(*ns_sol, RT(0.0));
1206 setVal(*ns_rhs, RT(0.0));
1207
1208 ns_linop->setLevelBC(0, ns_sol.get());
1209
1210 ns_mlmg = std::make_unique<MLMGT<MF>>(*ns_linop);
1211 ns_mlmg->setVerbose(0);
1212 ns_mlmg->setFixedIter(1);
1213 ns_mlmg->setMaxFmgIter(20);
1214 ns_mlmg->setBottomSolver(BottomSolver::smoother);
1215 }
1216}
1217
1218// in : Residual (res) on the finest AMR level
1219// out : sol on all AMR levels
1220template <typename MF>
1221void MLMGT<MF>::oneIter (int iter)
1222{
1223 BL_PROFILE("MLMG::oneIter()");
1224
1225 for (int alev = finest_amr_lev; alev > 0; --alev)
1226 {
1227 miniCycle(alev);
1228
1229 IntVect nghost(0);
1230 if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(alev)); }
1231 LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost);
1232
1233 // compute residual for the coarse AMR level
1234 computeResWithCrseSolFineCor(alev-1,alev);
1235
1236 if (alev != finest_amr_lev) {
1237 std::swap(cor_hold[alev][0], cor[alev][0]); // save it for the up cycle
1238 }
1239 }
1240
1241 // coarsest amr level
1242 {
1243 // enforce solvability if appropriate
1244 if (linop.isSingular(0) && linop.getEnforceSingularSolvable())
1245 {
1246 makeSolvable(0,0,res[0][0]);
1247 }
1248
1249 if (iter < max_fmg_iters) {
1250 mgFcycle();
1251 } else {
1252 mgVcycle(0, 0);
1253 }
1254
1255 IntVect nghost(0);
1256 if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(0)); }
1257 LocalAdd(sol[0], cor[0][0], 0, 0, ncomp, nghost);
1258 }
1259
1260 for (int alev = 1; alev <= finest_amr_lev; ++alev)
1261 {
1262 // (Fine AMR correction) = I(Coarse AMR correction)
1263 interpCorrection(alev);
1264
1265 IntVect nghost(0);
1266 if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(alev)); }
1267 LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost);
1268
1269 if (alev != finest_amr_lev) {
1270 LocalAdd(cor_hold[alev][0], cor[alev][0], 0, 0, ncomp, nghost);
1271 }
1272
1273 // Update fine AMR level correction
1274 computeResWithCrseCorFineCor(alev);
1275
1276 miniCycle(alev);
1277
1278 LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost);
1279
1280 if (alev != finest_amr_lev) {
1281 LocalAdd(cor[alev][0], cor_hold[alev][0], 0, 0, ncomp, nghost);
1282 }
1283 }
1284
1285 linop.averageDownAndSync(sol);
1286}
1287
1288template <typename MF>
1289void
1291{
1292 BL_PROFILE("MLMG::miniCycle()");
1293 const int mglev = 0;
1294 mgVcycle(amrlev, mglev);
1295}
1296
1297// in : Residual (res)
1298// out : Correction (cor) from bottom to this function's local top
1299template <typename MF>
1300void
1301MLMGT<MF>::mgVcycle (int amrlev, int mglev_top)
1302{
1303 BL_PROFILE("MLMG::mgVcycle()");
1304
1305 const int mglev_bottom = linop.NMGLevels(amrlev) - 1;
1306
1307 for (int mglev = mglev_top; mglev < mglev_bottom; ++mglev)
1308 {
1309 BL_PROFILE_VAR("MLMG::mgVcycle_down::"+std::to_string(mglev), blp_mgv_down_lev);
1310
1311 if (verbose >= 4)
1312 {
1313 RT norm = norminf(res[amrlev][mglev],0,ncomp,IntVect(0));
1314 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev
1315 << " DN: Norm before smooth " << norm << "\n";
1316 }
1317
1318 setVal(cor[amrlev][mglev], RT(0.0));
1319 bool skip_fillboundary = true;
1320 for (int i = 0; i < nu1; ++i) {
1321 linop.smooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev], skip_fillboundary);
1322 skip_fillboundary = false;
1323 }
1324
1325 // rescor = res - L(cor)
1326 computeResOfCorrection(amrlev, mglev);
1327
1328 if (verbose >= 4)
1329 {
1330 RT norm = norminf(rescor[amrlev][mglev],0,ncomp,IntVect(0));
1331 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev
1332 << " DN: Norm after smooth " << norm << "\n";
1333 }
1334
1335 // res_crse = R(rescor_fine); this provides res/b to the level below
1336 linop.restriction(amrlev, mglev+1, res[amrlev][mglev+1], rescor[amrlev][mglev]);
1337 }
1338
1339 BL_PROFILE_VAR("MLMG::mgVcycle_bottom", blp_bottom);
1340 if (amrlev == 0)
1341 {
1342 if (verbose >= 4)
1343 {
1344 RT norm = norminf(res[amrlev][mglev_bottom],0,ncomp,IntVect(0));
1345 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev_bottom
1346 << " DN: Norm before bottom " << norm << "\n";
1347 }
1348 bottomSolve();
1349 if (verbose >= 4)
1350 {
1351 computeResOfCorrection(amrlev, mglev_bottom);
1352 RT norm = norminf(rescor[amrlev][mglev_bottom],0,ncomp,IntVect(0));
1353 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev_bottom
1354 << " UP: Norm after bottom " << norm << "\n";
1355 }
1356 }
1357 else
1358 {
1359 if (verbose >= 4)
1360 {
1361 RT norm = norminf(res[amrlev][mglev_bottom],0,ncomp,IntVect(0));
1362 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev_bottom
1363 << " Norm before smooth " << norm << "\n";
1364 }
1365 setVal(cor[amrlev][mglev_bottom], RT(0.0));
1366 bool skip_fillboundary = true;
1367 for (int i = 0; i < nu1; ++i) {
1368 linop.smooth(amrlev, mglev_bottom, cor[amrlev][mglev_bottom],
1369 res[amrlev][mglev_bottom], skip_fillboundary);
1370 skip_fillboundary = false;
1371 }
1372 if (verbose >= 4)
1373 {
1374 computeResOfCorrection(amrlev, mglev_bottom);
1375 RT norm = norminf(rescor[amrlev][mglev_bottom],0,ncomp,IntVect(0));
1376 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev_bottom
1377 << " Norm after smooth " << norm << "\n";
1378 }
1379 }
1380 BL_PROFILE_VAR_STOP(blp_bottom);
1381
1382 for (int mglev = mglev_bottom-1; mglev >= mglev_top; --mglev)
1383 {
1384 BL_PROFILE_VAR("MLMG::mgVcycle_up::"+std::to_string(mglev), blp_mgv_up_lev);
1385 // cor_fine += I(cor_crse)
1386 addInterpCorrection(amrlev, mglev);
1387 if (verbose >= 4)
1388 {
1389 computeResOfCorrection(amrlev, mglev);
1390 RT norm = norminf(rescor[amrlev][mglev],0,ncomp,IntVect(0));
1391 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev
1392 << " UP: Norm before smooth " << norm << "\n";
1393 }
1394 for (int i = 0; i < nu2; ++i) {
1395 linop.smooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev]);
1396 }
1397
1398 if (cf_strategy == CFStrategy::ghostnodes) { computeResOfCorrection(amrlev, mglev); }
1399
1400 if (verbose >= 4)
1401 {
1402 computeResOfCorrection(amrlev, mglev);
1403 RT norm = norminf(rescor[amrlev][mglev],0,ncomp,IntVect(0));
1404 amrex::Print() << print_ident << "AT LEVEL " << amrlev << " " << mglev
1405 << " UP: Norm after smooth " << norm << "\n";
1406 }
1407 }
1408}
1409
1410// FMG cycle on the coarsest AMR level.
1411// in: Residual on the top MG level (i.e., 0)
1412// out: Correction (cor) on all MG levels
1413template <typename MF>
1414void
1416{
1417 BL_PROFILE("MLMG::mgFcycle()");
1418
1419#ifdef AMREX_USE_EB
1420 AMREX_ASSERT(linop.isCellCentered());
1421#endif
1422
1423 const int amrlev = 0;
1424 const int mg_bottom_lev = linop.NMGLevels(amrlev) - 1;
1425 IntVect nghost(0);
1426 if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(amrlev)); }
1427
1428 for (int mglev = 1; mglev <= mg_bottom_lev; ++mglev)
1429 {
1430 linop.avgDownResMG(mglev, res[amrlev][mglev], res[amrlev][mglev-1]);
1431 }
1432
1433 bottomSolve();
1434
1435 for (int mglev = mg_bottom_lev-1; mglev >= 0; --mglev)
1436 {
1437 // cor_fine = I(cor_crse)
1438 interpCorrection(amrlev, mglev);
1439
1440 // rescor = res - L(cor)
1441 computeResOfCorrection(amrlev, mglev);
1442 // res = rescor; this provides b to the vcycle below
1443 LocalCopy(res[amrlev][mglev], rescor[amrlev][mglev], 0, 0, ncomp, nghost);
1444
1445 // save cor; do v-cycle; add the saved to cor
1446 std::swap(cor[amrlev][mglev], cor_hold[amrlev][mglev]);
1447 mgVcycle(amrlev, mglev);
1448 LocalAdd(cor[amrlev][mglev], cor_hold[amrlev][mglev], 0, 0, ncomp, nghost);
1449 }
1450}
1451
1452// At the true bottom of the coarsest AMR level.
1453// in : Residual (res) as b
1454// out : Correction (cor) as x
1455template <typename MF>
1456void
1458{
1459 if (do_nsolve)
1460 {
1461 NSolve(*ns_mlmg, *ns_sol, *ns_rhs);
1462 }
1463 else
1464 {
1465 actualBottomSolve();
1466 }
1467}
1468
1469template <typename MF>
1470void
1471MLMGT<MF>::NSolve (MLMGT<MF>& a_solver, MF& a_sol, MF& a_rhs)
1472{
1473 BL_PROFILE("MLMG::NSolve()");
1474
1475 setVal(a_sol, RT(0.0));
1476
1477 MF const& res_bottom = res[0].back();
1478 if (BoxArray::SameRefs(boxArray(a_rhs),boxArray(res_bottom)) &&
1480 {
1481 LocalCopy(a_rhs, res_bottom, 0, 0, ncomp, IntVect(0));
1482 } else {
1483 setVal(a_rhs, RT(0.0));
1484 ParallelCopy(a_rhs, res_bottom, 0, 0, ncomp);
1485 }
1486
1487 a_solver.solve(Vector<MF*>{&a_sol}, Vector<MF const*>{&a_rhs},
1488 RT(-1.0), RT(-1.0));
1489
1490 linop.copyNSolveSolution(cor[0].back(), a_sol);
1491}
1492
1493template <typename MF>
1494void
1496{
1497 BL_PROFILE("MLMG::actualBottomSolve()");
1498
1499 if (!linop.isBottomActive()) { return; }
1500
1501 auto bottom_start_time = amrex::second();
1502
1503 ParallelContext::push(linop.BottomCommunicator());
1504
1505 const int amrlev = 0;
1506 const int mglev = linop.NMGLevels(amrlev) - 1;
1507 auto& x = cor[amrlev][mglev];
1508 auto& b = res[amrlev][mglev];
1509
1510 setVal(x, RT(0.0));
1511
1512 if (bottom_solver == BottomSolver::smoother)
1513 {
1514 bool skip_fillboundary = true;
1515 for (int i = 0; i < nuf; ++i) {
1516 linop.smooth(amrlev, mglev, x, b, skip_fillboundary);
1517 skip_fillboundary = false;
1518 }
1519 }
1520 else
1521 {
1522 MF* bottom_b = &b;
1523 MF raii_b;
1524 if (linop.isBottomSingular() && linop.getEnforceSingularSolvable())
1525 {
1526 const IntVect ng = nGrowVect(b);
1527 raii_b = linop.make(amrlev, mglev, ng);
1528 LocalCopy(raii_b, b, 0, 0, ncomp, ng);
1529 bottom_b = &raii_b;
1530
1531 makeSolvable(amrlev,mglev,*bottom_b);
1532 }
1533
1534 if (bottom_solver == BottomSolver::hypre)
1535 {
1536#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
1537 if constexpr (std::is_same<MF,MultiFab>()) {
1538 bottomSolveWithHypre(x, *bottom_b);
1539 } else
1540#endif
1541 {
1542 amrex::Abort("Using Hypre as bottom solver not supported in this case");
1543 }
1544 }
1545 else if (bottom_solver == BottomSolver::petsc)
1546 {
1547#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
1548 if constexpr (std::is_same<MF,MultiFab>()) {
1549 bottomSolveWithPETSc(x, *bottom_b);
1550 } else
1551#endif
1552 {
1553 amrex::Abort("Using PETSc as bottom solver not supported in this case");
1554 }
1555 }
1556 else
1557 {
1558 typename MLCGSolverT<MF>::Type cg_type;
1559 if (bottom_solver == BottomSolver::cg ||
1560 bottom_solver == BottomSolver::cgbicg) {
1561 cg_type = MLCGSolverT<MF>::Type::CG;
1562 } else {
1564 }
1565
1566 int ret = bottomSolveWithCG(x, *bottom_b, cg_type);
1567
1568 if (ret != 0 && (bottom_solver == BottomSolver::cgbicg ||
1569 bottom_solver == BottomSolver::bicgcg))
1570 {
1571 if (bottom_solver == BottomSolver::cgbicg) {
1572 cg_type = MLCGSolverT<MF>::Type::BiCGStab; // switch to bicg
1573 } else {
1574 cg_type = MLCGSolverT<MF>::Type::CG; // switch to cg
1575 }
1576 setVal(cor[amrlev][mglev], RT(0.0));
1577 ret = bottomSolveWithCG(x, *bottom_b, cg_type);
1578 if (ret == 0) { // switch permanently
1579 if (cg_type == MLCGSolverT<MF>::Type::CG) {
1580 bottom_solver = BottomSolver::cg;
1581 } else {
1582 bottom_solver = BottomSolver::bicgstab;
1583 }
1584 }
1585 }
1586
1587 // If the bottom solve failed then set the correction to zero
1588 if (ret != 0 && ret != 9) {
1589 setVal(cor[amrlev][mglev], RT(0.0));
1590 }
1591 const int n = (ret==0) ? nub : nuf;
1592 for (int i = 0; i < n; ++i) {
1593 linop.smooth(amrlev, mglev, x, b);
1594 }
1595 }
1596 }
1597
1599
1600 if (! timer.empty()) {
1601 timer[bottom_time] += amrex::second() - bottom_start_time;
1602 }
1603}
1604
1605template <typename MF>
1606int
1608{
1609 MLCGSolverT<MF> cg_solver(linop);
1610 cg_solver.setSolver(type);
1611 cg_solver.setVerbose(bottom_verbose);
1612 cg_solver.setPrintIdentation(print_ident);
1613 cg_solver.setMaxIter(bottom_maxiter);
1614 cg_solver.setInitSolnZeroed(true);
1615 if (cf_strategy == CFStrategy::ghostnodes) { cg_solver.setNGhost(linop.getNGrow()); }
1616
1617 int ret = cg_solver.solve(x, b, bottom_reltol, bottom_abstol);
1618 if (ret != 0 && verbose > 1) {
1619 amrex::Print() << print_ident << "MLMG: Bottom solve failed.\n";
1620 }
1621 m_niters_cg.push_back(cg_solver.getNumIters());
1622 return ret;
1623}
1624
1625// Compute multi-level Residual (res) up to amrlevmax.
1626template <typename MF>
1627void
1629{
1630 BL_PROFILE("MLMG::computeMLResidual()");
1631
1632 const int mglev = 0;
1633 for (int alev = amrlevmax; alev >= 0; --alev) {
1634 const MF* crse_bcdata = (alev > 0) ? &(sol[alev-1]) : nullptr;
1635 linop.solutionResidual(alev, res[alev][mglev], sol[alev], rhs[alev], crse_bcdata);
1636 if (alev < finest_amr_lev) {
1637 linop.reflux(alev, res[alev][mglev], sol[alev], rhs[alev],
1638 res[alev+1][mglev], sol[alev+1], rhs[alev+1]);
1639 }
1640 }
1641}
1642
1643// Compute single AMR level residual without masking.
1644template <typename MF>
1645void
1647{
1648 BL_PROFILE("MLMG::computeResidual()");
1649 const MF* crse_bcdata = (alev > 0) ? &(sol[alev-1]) : nullptr;
1650 linop.solutionResidual(alev, res[alev][0], sol[alev], rhs[alev], crse_bcdata);
1651}
1652
1653// Compute coarse AMR level composite residual with coarse solution and fine correction
1654template <typename MF>
1655void
1657{
1658 BL_PROFILE("MLMG::computeResWithCrseSolFineCor()");
1659
1660 IntVect nghost(0);
1661 if (cf_strategy == CFStrategy::ghostnodes) {
1662 nghost = IntVect(std::min(linop.getNGrow(falev),linop.getNGrow(calev)));
1663 }
1664
1665 MF& crse_sol = sol[calev];
1666 const MF& crse_rhs = rhs[calev];
1667 MF& crse_res = res[calev][0];
1668
1669 MF& fine_sol = sol[falev];
1670 const MF& fine_rhs = rhs[falev];
1671 MF& fine_cor = cor[falev][0];
1672 MF& fine_res = res[falev][0];
1673 MF& fine_rescor = rescor[falev][0];
1674
1675 const MF* crse_bcdata = (calev > 0) ? &(sol[calev-1]) : nullptr;
1676 linop.solutionResidual(calev, crse_res, crse_sol, crse_rhs, crse_bcdata);
1677
1678 linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res, BCMode::Homogeneous);
1679 LocalCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost);
1680
1681 linop.reflux(calev, crse_res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs);
1682
1683 linop.avgDownResAmr(calev, crse_res, fine_res);
1684}
1685
1686// Compute fine AMR level residual fine_res = fine_res - L(fine_cor) with coarse providing BC.
1687template <typename MF>
1688void
1690{
1691 BL_PROFILE("MLMG::computeResWithCrseCorFineCor()");
1692
1693 IntVect nghost(0);
1694 if (cf_strategy == CFStrategy::ghostnodes) {
1695 nghost = IntVect(linop.getNGrow(falev));
1696 }
1697
1698 const MF& crse_cor = cor[falev-1][0];
1699
1700 MF& fine_cor = cor [falev][0];
1701 MF& fine_res = res [falev][0];
1702 MF& fine_rescor = rescor[falev][0];
1703
1704 // fine_rescor = fine_res - L(fine_cor)
1705 linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res,
1706 BCMode::Inhomogeneous, &crse_cor);
1707 LocalCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost);
1708}
1709
1710// Interpolate correction from coarse to fine AMR level.
1711template <typename MF>
1712void
1714{
1715 BL_PROFILE("MLMG::interpCorrection_1");
1716
1717 IntVect nghost(0);
1718 if (cf_strategy == CFStrategy::ghostnodes) {
1719 nghost = IntVect(linop.getNGrow(alev));
1720 }
1721
1722 MF const& crse_cor = cor[alev-1][0];
1723 MF & fine_cor = cor[alev ][0];
1724
1725 const Geometry& crse_geom = linop.Geom(alev-1,0);
1726
1727 int ng_src = 0;
1728 int ng_dst = linop.isCellCentered() ? 1 : 0;
1729 if (cf_strategy == CFStrategy::ghostnodes)
1730 {
1731 ng_src = linop.getNGrow(alev-1);
1732 ng_dst = linop.getNGrow(alev-1);
1733 }
1734
1735 MF cfine = linop.makeCoarseAmr(alev, IntVect(ng_dst));
1736 setVal(cfine, RT(0.0));
1737 ParallelCopy(cfine, crse_cor, 0, 0, ncomp, IntVect(ng_src), IntVect(ng_dst),
1738 crse_geom.periodicity());
1739
1740 linop.interpolationAmr(alev, fine_cor, cfine, nghost); // NOLINT(readability-suspicious-call-argument)
1741}
1742
1743// Interpolate correction between MG levels
1744// inout: Correction (cor) on coarse MG lev. (out due to FillBoundary)
1745// out : Correction (cor) on fine MG lev.
1746template <typename MF>
1747void
1748MLMGT<MF>::interpCorrection (int alev, int mglev)
1749{
1750 BL_PROFILE("MLMG::interpCorrection_2");
1751
1752 MF& crse_cor = cor[alev][mglev+1];
1753 MF& fine_cor = cor[alev][mglev ];
1754 linop.interpAssign(alev, mglev, fine_cor, crse_cor);
1755}
1756
1757// (Fine MG level correction) += I(Coarse MG level correction)
1758template <typename MF>
1759void
1761{
1762 BL_PROFILE("MLMG::addInterpCorrection()");
1763
1764 const MF& crse_cor = cor[alev][mglev+1];
1765 MF& fine_cor = cor[alev][mglev ];
1766
1767 MF cfine;
1768 const MF* cmf;
1769
1770 if (linop.isMFIterSafe(alev, mglev, mglev+1))
1771 {
1772 cmf = &crse_cor;
1773 }
1774 else
1775 {
1776 cfine = linop.makeCoarseMG(alev, mglev, IntVect(0));
1777 ParallelCopy(cfine, crse_cor, 0, 0, ncomp);
1778 cmf = &cfine;
1779 }
1780
1781 linop.interpolation(alev, mglev, fine_cor, *cmf);
1782}
1783
1784// Compute rescor = res - L(cor)
1785// in : res
1786// inout: cor (out due to FillBoundary in linop.correctionResidual)
1787// out : rescor
1788template <typename MF>
1789void
1791{
1792 BL_PROFILE("MLMG:computeResOfCorrection()");
1793 MF & x = cor[amrlev][mglev];
1794 const MF& b = res[amrlev][mglev];
1795 MF & r = rescor[amrlev][mglev];
1796 linop.correctionResidual(amrlev, mglev, r, x, b, BCMode::Homogeneous);
1797}
1798
1799// Compute single-level masked inf-norm of Residual (res).
1800template <typename MF>
1801auto
1802MLMGT<MF>::ResNormInf (int alev, bool local) -> RT
1803{
1804 BL_PROFILE("MLMG::ResNormInf()");
1805 return linop.normInf(alev, res[alev][0], local);
1806}
1807
1808// Computes multi-level masked inf-norm of Residual (res).
1809template <typename MF>
1810auto
1811MLMGT<MF>::MLResNormInf (int alevmax, bool local) -> RT
1812{
1813 BL_PROFILE("MLMG::MLResNormInf()");
1814 RT r = RT(0.0);
1815 for (int alev = 0; alev <= alevmax; ++alev)
1816 {
1817 r = std::max(r, ResNormInf(alev,true));
1818 }
1820 return r;
1821}
1822
1823// Compute multi-level masked inf-norm of RHS (rhs).
1824template <typename MF>
1825auto
1827{
1828 BL_PROFILE("MLMG::MLRhsNormInf()");
1829 RT r = RT(0.0);
1830 for (int alev = 0; alev <= finest_amr_lev; ++alev) {
1831 auto t = linop.normInf(alev, rhs[alev], true);
1832 r = std::max(r, t);
1833 }
1835 return r;
1836}
1837
1838template <typename MF>
1839void
1841{
1842 auto const& offset = linop.getSolvabilityOffset(0, 0, rhs[0]);
1843 if (verbose >= 4) {
1844 for (int c = 0; c < ncomp; ++c) {
1845 amrex::Print() << print_ident << "MLMG: Subtracting " << offset[c] << " from rhs component "
1846 << c << "\n";
1847 }
1848 }
1849 for (int alev = 0; alev < namrlevs; ++alev) {
1850 linop.fixSolvabilityByOffset(alev, 0, rhs[alev], offset);
1851 }
1852}
1853
1854template <typename MF>
1855void
1856MLMGT<MF>::makeSolvable (int amrlev, int mglev, MF& mf)
1857{
1858 auto const& offset = linop.getSolvabilityOffset(amrlev, mglev, mf);
1859 if (verbose >= 4) {
1860 for (int c = 0; c < ncomp; ++c) {
1861 amrex::Print() << print_ident << "MLMG: Subtracting " << offset[c]
1862 << " from mf component c = " << c
1863 << " on level (" << amrlev << ", " << mglev << ")\n";
1864 }
1865 }
1866 linop.fixSolvabilityByOffset(amrlev, mglev, mf, offset);
1867}
1868
1869#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
1870template <typename MF>
1871template <class TMF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,int>>
1872void
1873MLMGT<MF>::bottomSolveWithHypre (MF& x, const MF& b)
1874{
1875 const int amrlev = 0;
1876 const int mglev = linop.NMGLevels(amrlev) - 1;
1877
1878 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ncomp == 1, "bottomSolveWithHypre doesn't work with ncomp > 1");
1879
1880 if (linop.isCellCentered())
1881 {
1882 if (hypre_solver == nullptr) // We should reuse the setup
1883 {
1884 hypre_solver = linop.makeHypre(hypre_interface);
1885
1886 hypre_solver->setVerbose(bottom_verbose);
1887 if (hypre_interface == amrex::Hypre::Interface::ij) {
1888 hypre_solver->setHypreOptionsNamespace(hypre_options_namespace);
1889 } else {
1890 hypre_solver->setHypreOldDefault(hypre_old_default);
1891 hypre_solver->setHypreRelaxType(hypre_relax_type);
1892 hypre_solver->setHypreRelaxOrder(hypre_relax_order);
1893 hypre_solver->setHypreNumSweeps(hypre_num_sweeps);
1894 hypre_solver->setHypreStrongThreshold(hypre_strong_threshold);
1895 }
1896
1897 const BoxArray& ba = linop.m_grids[amrlev].back();
1898 const DistributionMapping& dm = linop.m_dmap[amrlev].back();
1899 const Geometry& geom = linop.m_geom[amrlev].back();
1900
1901 hypre_bndry = std::make_unique<MLMGBndryT<MF>>(ba, dm, ncomp, geom);
1902 hypre_bndry->setHomogValues();
1903 const Real* dx = linop.m_geom[0][0].CellSize();
1904 IntVect crse_ratio = linop.m_coarse_data_crse_ratio.allGT(0) ? linop.m_coarse_data_crse_ratio : IntVect(1);
1905 RealVect bclocation(AMREX_D_DECL(0.5*dx[0]*crse_ratio[0],
1906 0.5*dx[1]*crse_ratio[1],
1907 0.5*dx[2]*crse_ratio[2]));
1908 hypre_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc, IntVect(-1), bclocation,
1909 linop.m_coarse_fine_bc_type);
1910 }
1911
1912 // IJ interface understands absolute tolerance API of hypre
1913 amrex::Real hypre_abstol =
1914 (hypre_interface == amrex::Hypre::Interface::ij)
1915 ? bottom_abstol : Real(-1.0);
1916 hypre_solver->solve(
1917 x, b, bottom_reltol, hypre_abstol, bottom_maxiter, *hypre_bndry,
1918 linop.getMaxOrder());
1919 }
1920 else
1921 {
1922 if (hypre_node_solver == nullptr)
1923 {
1924 hypre_node_solver =
1925 linop.makeHypreNodeLap(bottom_verbose, hypre_options_namespace);
1926 }
1927 hypre_node_solver->solve(x, b, bottom_reltol, bottom_abstol, bottom_maxiter);
1928 }
1929
1930 // For singular problems there may be a large constant added to all values of the solution
1931 // For precision reasons we enforce that the average of the correction from hypre is 0
1932 if (linop.isSingular(amrlev) && linop.getEnforceSingularSolvable())
1933 {
1934 makeSolvable(amrlev, mglev, x);
1935 }
1936}
1937#endif
1938
1939#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
1940template <typename MF>
1941template <class TMF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,int>>
1942void
1943MLMGT<MF>::bottomSolveWithPETSc (MF& x, const MF& b)
1944{
1945 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ncomp == 1, "bottomSolveWithPETSc doesn't work with ncomp > 1");
1946
1947 if(petsc_solver == nullptr)
1948 {
1949 petsc_solver = linop.makePETSc();
1950 petsc_solver->setVerbose(bottom_verbose);
1951
1952 const BoxArray& ba = linop.m_grids[0].back();
1953 const DistributionMapping& dm = linop.m_dmap[0].back();
1954 const Geometry& geom = linop.m_geom[0].back();
1955
1956 petsc_bndry = std::make_unique<MLMGBndryT<MF>>(ba, dm, ncomp, geom);
1957 petsc_bndry->setHomogValues();
1958 const Real* dx = linop.m_geom[0][0].CellSize();
1959 auto crse_ratio = linop.m_coarse_data_crse_ratio.allGT(0) ? linop.m_coarse_data_crse_ratio : IntVect(1);
1960 RealVect bclocation(AMREX_D_DECL(0.5*dx[0]*crse_ratio[0],
1961 0.5*dx[1]*crse_ratio[1],
1962 0.5*dx[2]*crse_ratio[2]));
1963 petsc_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc, IntVect(-1), bclocation,
1964 linop.m_coarse_fine_bc_type);
1965 }
1966 petsc_solver->solve(x, b, bottom_reltol, Real(-1.), bottom_maxiter, *petsc_bndry,
1967 linop.getMaxOrder());
1968}
1969#endif
1970
1971template <typename MF>
1972void
1974 const Vector<MultiFab const*>& a_rhs,
1975 RT a_tol_rel, RT a_tol_abs, const char* a_file_name) const
1976{
1977 std::string file_name(a_file_name);
1978 UtilCreateCleanDirectory(file_name, false);
1979
1981 {
1982 std::string HeaderFileName(std::string(a_file_name)+"/Header");
1983 std::ofstream HeaderFile;
1984 HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out |
1985 std::ofstream::trunc |
1986 std::ofstream::binary);
1987 if( ! HeaderFile.good()) {
1988 FileOpenFailed(HeaderFileName);
1989 }
1990
1991 HeaderFile.precision(17);
1992
1993 HeaderFile << linop.name() << "\n"
1994 << "a_tol_rel = " << a_tol_rel << "\n"
1995 << "a_tol_abs = " << a_tol_abs << "\n"
1996 << "verbose = " << verbose << "\n"
1997 << "max_iters = " << max_iters << "\n"
1998 << "nu1 = " << nu1 << "\n"
1999 << "nu2 = " << nu2 << "\n"
2000 << "nuf = " << nuf << "\n"
2001 << "nub = " << nub << "\n"
2002 << "max_fmg_iters = " << max_fmg_iters << "\n"
2003 << "bottom_solver = " << static_cast<int>(bottom_solver) << "\n"
2004 << "bottom_verbose = " << bottom_verbose << "\n"
2005 << "bottom_maxiter = " << bottom_maxiter << "\n"
2006 << "bottom_reltol = " << bottom_reltol << "\n"
2007 << "always_use_bnorm = " << always_use_bnorm << "\n"
2008 << "namrlevs = " << namrlevs << "\n"
2009 << "finest_amr_lev = " << finest_amr_lev << "\n"
2010 << "linop_prepared = " << linop_prepared << "\n"
2011 << "solve_called = " << solve_called << "\n";
2012
2013 for (int ilev = 0; ilev <= finest_amr_lev; ++ilev) {
2014 UtilCreateCleanDirectory(file_name+"/Level_"+std::to_string(ilev), false);
2015 }
2016 }
2017
2019
2020 for (int ilev = 0; ilev <= finest_amr_lev; ++ilev) {
2021 VisMF::Write(*a_sol[ilev], file_name+"/Level_"+std::to_string(ilev)+"/sol");
2022 VisMF::Write(*a_rhs[ilev], file_name+"/Level_"+std::to_string(ilev)+"/rhs");
2023 }
2024
2025 linop.checkPoint(file_name+"/linop");
2026}
2027
2028template <typename MF>
2029void
2031{
2032 print_ident.resize(print_ident.size()+4, ' ');
2033}
2034
2035template <typename MF>
2036void
2038{
2039 if (print_ident.size() > 4) {
2040 print_ident.resize(print_ident.size()-4, ' ');
2041 } else {
2042 print_ident.clear();
2043 }
2044}
2045
2046extern template class MLMGT<MultiFab>;
2047
2049
2050}
2051
2052#endif
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define BL_PROFILE_VAR_STOP(vname)
Definition AMReX_BLProfiler.H:563
#define BL_PROFILE_VAR(fname, vname)
Definition AMReX_BLProfiler.H:560
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
Print on all processors of the default communicator.
Definition AMReX_Print.H:117
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:550
static bool SameRefs(const BoxArray &lhs, const BoxArray &rhs)
whether two BoxArrays share the same data
Definition AMReX_BoxArray.H:820
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:41
static bool SameRefs(const DistributionMapping &lhs, const DistributionMapping &rhs)
Definition AMReX_DistributionMapping.H:187
Definition AMReX_EBFabFactory.H:24
Solve using GMRES with multigrid as preconditioner.
Definition AMReX_GMRES_MLMG.H:21
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:73
Periodicity periodicity() const noexcept
Definition AMReX_Geometry.H:355
Interface
Definition AMReX_Hypre.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool allGT(const IntVectND< dim > &rhs) const noexcept
Returns true if this is greater than argument for all components. NOTE: This is NOT a strict weak ord...
Definition AMReX_IntVect.H:418
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr IntVectND< dim > TheDimensionVector(int d) noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:691
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE IntVectND< new_dim > resize(int fill_extra=0) const noexcept
Returns a new IntVectND of size new_dim by either shrinking or expanding this IntVectND.
Definition AMReX_IntVect.H:768
Definition AMReX_MLCGSolver.H:12
void setSolver(Type _typ) noexcept
Definition AMReX_MLCGSolver.H:28
void setVerbose(int _verbose)
Definition AMReX_MLCGSolver.H:39
int getNumIters() const noexcept
Definition AMReX_MLCGSolver.H:64
void setInitSolnZeroed(bool _sol_zeroed)
Definition AMReX_MLCGSolver.H:53
void setPrintIdentation(std::string s)
Definition AMReX_MLCGSolver.H:45
int solve(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition AMReX_MLCGSolver.H:89
Type
Definition AMReX_MLCGSolver.H:18
void setNGhost(int _nghost)
Definition AMReX_MLCGSolver.H:56
void setMaxIter(int _maxiter)
Definition AMReX_MLCGSolver.H:42
Definition AMReX_MLLinOp.H:98
typename FabDataType< MF >::fab_type FAB
Definition AMReX_MLLinOp.H:108
typename FabDataType< MF >::value_type RT
Definition AMReX_MLLinOp.H:109
Definition AMReX_MLMG.H:17
Definition AMReX_MLMG.H:12
void prepareForFluxes(Vector< MF const * > const &a_sol)
Definition AMReX_MLMG.H:579
void setBottomVerbose(int v) noexcept
Definition AMReX_MLMG.H:143
void setMaxFmgIter(int n) noexcept
Definition AMReX_MLMG.H:131
RT MLResNormInf(int alevmax, bool local=false)
Definition AMReX_MLMG.H:1811
RT MLRhsNormInf(bool local=false)
Definition AMReX_MLMG.H:1826
MLMGT(MLMGT< MF > &&)=delete
void actualBottomSolve()
Definition AMReX_MLMG.H:1495
Vector< Vector< MF > > rescor
Definition AMReX_MLMG.H:331
MF MFType
Definition AMReX_MLMG.H:25
BottomSolver getBottomSolver() const noexcept
Definition AMReX_MLMG.H:141
bool linop_prepared
Definition AMReX_MLMG.H:277
int bottom_verbose
Definition AMReX_MLMG.H:263
void checkPoint(const Vector< MultiFab * > &a_sol, const Vector< MultiFab const * > &a_rhs, RT a_tol_rel, RT a_tol_abs, const char *a_file_name) const
Definition AMReX_MLMG.H:1973
void setPreSmooth(int n) noexcept
Definition AMReX_MLMG.H:135
void setBottomToleranceAbs(RT t) noexcept
Definition AMReX_MLMG.H:146
Vector< Vector< MF > > cor
= rhs - L(sol)
Definition AMReX_MLMG.H:329
RT getFinalResidual() const noexcept
Definition AMReX_MLMG.H:236
int do_nsolve
N Solve.
Definition AMReX_MLMG.H:281
void interpCorrection(int alev)
Definition AMReX_MLMG.H:1713
Long solve_called
Definition AMReX_MLMG.H:278
void setBottomSmooth(int n) noexcept
Definition AMReX_MLMG.H:138
int final_fill_bc
Definition AMReX_MLMG.H:270
bool precond_mode
Definition AMReX_MLMG.H:246
void setNSolve(int flag) noexcept
Definition AMReX_MLMG.H:155
int getBottomVerbose() const
Definition AMReX_MLMG.H:123
std::unique_ptr< MF > ns_sol
Definition AMReX_MLMG.H:285
std::string print_ident
Definition AMReX_MLMG.H:288
Vector< Vector< MF > > cor_hold
Definition AMReX_MLMG.H:330
void computeResOfCorrection(int amrlev, int mglev)
Definition AMReX_MLMG.H:1790
void applyPrecond(const Vector< MF * > &out, const Vector< MF * > &in)
out = L(in) as a preconditioner
Definition AMReX_MLMG.H:970
void setCFStrategy(CFStrategy a_cf_strategy) noexcept
Definition AMReX_MLMG.H:142
void computeResWithCrseCorFineCor(int falev)
Definition AMReX_MLMG.H:1689
void NSolve(MLMGT< MF > &a_solver, MF &a_sol, MF &a_rhs)
Definition AMReX_MLMG.H:1471
typename MLLinOpT< MF >::Location Location
Definition AMReX_MLMG.H:30
RT m_final_resnorm0
Definition AMReX_MLMG.H:339
void apply(const Vector< MF * > &out, const Vector< MF * > &in)
out = L(in). Note that, if no actual solve is needed, one could turn off multigrid coarsening by cons...
Definition AMReX_MLMG.H:894
Vector< MF > rhs
Definition AMReX_MLMG.H:319
void setNSolveGridSize(int s) noexcept
Definition AMReX_MLMG.H:156
int nuf
when smoother is used as bottom solver
Definition AMReX_MLMG.H:256
void setVerbose(int v) noexcept
Definition AMReX_MLMG.H:129
int nu1
pre
Definition AMReX_MLMG.H:254
void computeMLResidual(int amrlevmax)
Definition AMReX_MLMG.H:1628
Vector< RT > m_iter_fine_resnorm0
Definition AMReX_MLMG.H:341
RT getInitResidual() const noexcept
Definition AMReX_MLMG.H:234
int getNumIters() const noexcept
Definition AMReX_MLMG.H:239
void setPostSmooth(int n) noexcept
Definition AMReX_MLMG.H:136
BottomSolver bottom_solver
Definition AMReX_MLMG.H:261
void mgVcycle(int amrlev, int mglev)
Definition AMReX_MLMG.H:1301
MLLinOpT< MF > & linop
Definition AMReX_MLMG.H:272
int ncomp
Definition AMReX_MLMG.H:273
void prepareForNSolve()
Definition AMReX_MLMG.H:1188
RT precond(Vector< MF * > const &a_sol, Vector< MF const * > const &a_rhs, RT a_tol_rel, RT a_tol_abs)
Definition AMReX_MLMG.H:561
int always_use_bnorm
Definition AMReX_MLMG.H:268
int max_fmg_iters
Definition AMReX_MLMG.H:259
void makeSolvable()
Definition AMReX_MLMG.H:1840
void setBottomSolver(BottomSolver s) noexcept
Definition AMReX_MLMG.H:140
int namrlevs
Definition AMReX_MLMG.H:274
void preparePrecond()
Definition AMReX_MLMG.H:1180
int nub
additional smoothing after bottom cg solver
Definition AMReX_MLMG.H:257
void incPrintIdentation()
Definition AMReX_MLMG.H:2030
typename MLLinOpT< MF >::RT RT
Definition AMReX_MLMG.H:27
int max_iters
Definition AMReX_MLMG.H:250
int do_fixed_number_of_iters
Definition AMReX_MLMG.H:251
void setThrowException(bool t) noexcept
Definition AMReX_MLMG.H:128
void decPrintIdentation()
Definition AMReX_MLMG.H:2037
void setFixedIter(int nit) noexcept
Definition AMReX_MLMG.H:132
int finest_amr_lev
Definition AMReX_MLMG.H:275
void getGradSolution(const Vector< Array< AMF *, AMREX_SPACEDIM > > &a_grad_sol, Location a_loc=Location::FaceCenter)
Definition AMReX_MLMG.H:590
Vector< RT > const & getResidualHistory() const noexcept
Definition AMReX_MLMG.H:238
void prepareLinOp()
Definition AMReX_MLMG.H:1168
RT m_rhsnorm0
Definition AMReX_MLMG.H:337
void setPrecondIter(int nit) noexcept
Definition AMReX_MLMG.H:133
CFStrategy
Definition AMReX_MLMG.H:33
void prepareForSolve(Vector< AMF * > const &a_sol, Vector< AMF const * > const &a_rhs)
Definition AMReX_MLMG.H:982
int bottomSolveWithCG(MF &x, const MF &b, typename MLCGSolverT< MF >::Type type)
Definition AMReX_MLMG.H:1607
std::unique_ptr< MLMGT< MF > > ns_mlmg
Definition AMReX_MLMG.H:284
void setAlwaysUseBNorm(int flag) noexcept
Definition AMReX_MLMG.H:149
void compResidual(const Vector< MF * > &a_res, const Vector< MF * > &a_sol, const Vector< MF const * > &a_rhs)
Definition AMReX_MLMG.H:826
void miniCycle(int amrlev)
Definition AMReX_MLMG.H:1290
RT solve(std::initializer_list< AMF * > a_sol, std::initializer_list< AMF const * > a_rhs, RT a_tol_rel, RT a_tol_abs, const char *checkpoint_file=nullptr)
int bottom_maxiter
Definition AMReX_MLMG.H:264
bool throw_exception
Definition AMReX_MLMG.H:247
RT bottom_reltol
Definition AMReX_MLMG.H:265
Vector< Vector< MF > > res
First Vector: Amr levels. 0 is the coarest level Second Vector: MG levels. 0 is the finest level.
Definition AMReX_MLMG.H:328
void setFinalFillBC(int flag) noexcept
Definition AMReX_MLMG.H:151
typename MLLinOpT< MF >::BCMode BCMode
Definition AMReX_MLMG.H:29
void computeResWithCrseSolFineCor(int calev, int falev)
Definition AMReX_MLMG.H:1656
MLMGT< MF > & operator=(MLMGT< MF > const &)=delete
Vector< MF > sol
Hypre.
Definition AMReX_MLMG.H:318
int nsolve_grid_size
Definition AMReX_MLMG.H:282
MLMGT(MLLinOpT< MF > &a_lp)
Definition AMReX_MLMG.H:350
int verbose
Definition AMReX_MLMG.H:248
void computeResidual(int alev)
Definition AMReX_MLMG.H:1646
MLLinOpT< MF > & getLinOp()
Definition AMReX_MLMG.H:242
std::unique_ptr< MF > ns_rhs
Definition AMReX_MLMG.H:286
timer_types
Definition AMReX_MLMG.H:334
@ ntimers
Definition AMReX_MLMG.H:334
@ bottom_time
Definition AMReX_MLMG.H:334
@ iter_time
Definition AMReX_MLMG.H:334
@ solve_time
Definition AMReX_MLMG.H:334
CFStrategy cf_strategy
Definition AMReX_MLMG.H:262
typename MLLinOpT< MF >::FAB FAB
Definition AMReX_MLMG.H:26
Vector< int > sol_is_alias
Definition AMReX_MLMG.H:322
RT getBottomToleranceAbs() const noexcept
Definition AMReX_MLMG.H:147
int numAMRLevels() const noexcept
Definition AMReX_MLMG.H:153
RT bottom_abstol
Definition AMReX_MLMG.H:266
Vector< int > m_niters_cg
Definition AMReX_MLMG.H:340
MLMGT(MLMGT< MF > const &)=delete
void mgFcycle()
Definition AMReX_MLMG.H:1415
Vector< double > timer
Definition AMReX_MLMG.H:335
RT getInitRHS() const noexcept
Definition AMReX_MLMG.H:232
RT ResNormInf(int alev, bool local=false)
Definition AMReX_MLMG.H:1802
Vector< int > const & getNumCGIters() const noexcept
Definition AMReX_MLMG.H:240
int nu2
post
Definition AMReX_MLMG.H:255
int max_precond_iters
Definition AMReX_MLMG.H:252
void bottomSolve()
Definition AMReX_MLMG.H:1457
void getFluxes(const Vector< Array< AMF *, AMREX_SPACEDIM > > &a_flux, Location a_loc=Location::FaceCenter)
For (alpha * a - beta * (del dot b grad)) phi = rhs, flux means -b grad phi
Definition AMReX_MLMG.H:621
void setBottomTolerance(RT t) noexcept
Definition AMReX_MLMG.H:145
void setFinalSmooth(int n) noexcept
Definition AMReX_MLMG.H:137
std::unique_ptr< MLLinOpT< MF > > ns_linop
Definition AMReX_MLMG.H:283
void addInterpCorrection(int alev, int mglev)
Definition AMReX_MLMG.H:1760
RT m_init_resnorm0
Definition AMReX_MLMG.H:338
int getVerbose() const
Definition AMReX_MLMG.H:122
RT solve(const Vector< AMF * > &a_sol, const Vector< AMF const * > &a_rhs, RT a_tol_rel, RT a_tol_abs, const char *checkpoint_file=nullptr)
void setBottomMaxIter(int n) noexcept
Definition AMReX_MLMG.H:144
void oneIter(int iter)
Definition AMReX_MLMG.H:1221
void setMaxIter(int n) noexcept
Definition AMReX_MLMG.H:130
This class provides the user with a few print options.
Definition AMReX_Print.H:35
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Long size() const noexcept
Definition AMReX_Vector.H:50
static Long Write(const FabArray< FArrayBox > &mf, const std::string &name, VisMF::How how=NFiles, bool set_ghost=false)
Write a FabArray<FArrayBox> to disk in a "smart" way. Returns the total number of bytes written on th...
Definition AMReX_VisMF.cpp:975
void Max(KeyValuePair< K, V > &vi, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:126
void push(MPI_Comm c)
Definition AMReX_ParallelContext.H:102
void BarrierSub() noexcept
Definition AMReX_ParallelContext.H:88
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
int MyProcSub() noexcept
my sub-rank in current frame
Definition AMReX_ParallelContext.H:76
void pop()
Note that it's the user's responsibility to free the MPI_Comm.
Definition AMReX_ParallelContext.H:108
bool IOProcessorSub() noexcept
Am IO processor for current frame?
Definition AMReX_ParallelContext.H:80
Definition AMReX_Amr.cpp:49
MF::value_type norminf(MF const &mf, int scomp, int ncomp, IntVect const &nghost, bool local=false)
Definition AMReX_FabArrayUtility.H:1883
int nComp(FabArrayBase const &fa)
void FileOpenFailed(const std::string &file)
Output a message and abort when couldn't open the file.
Definition AMReX_Utility.cpp:131
DistributionMapping const & DistributionMap(FabArrayBase const &fa)
void average_down(const MultiFab &S_fine, MultiFab &S_crse, const Geometry &fgeom, const Geometry &cgeom, int scomp, int ncomp, int rr)
Definition AMReX_MultiFabUtil.cpp:314
IntVect nGrowVect(FabArrayBase const &fa)
void EB_average_face_to_cellcenter(MultiFab &ccmf, int dcomp, const Array< MultiFab const *, AMREX_SPACEDIM > &fmf)
Definition AMReX_EBMultiFabUtil.cpp:806
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Returns a BoxND with different type.
Definition AMReX_Box.H:1435
void LocalCopy(DMF &dst, SMF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst = src
Definition AMReX_FabArrayUtility.H:1831
void EB_set_covered(MultiFab &mf, Real val)
Definition AMReX_EBMultiFabUtil.cpp:21
void LocalAdd(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst += src
Definition AMReX_FabArrayUtility.H:1839
double second() noexcept
Definition AMReX_Utility.cpp:922
void UtilCreateCleanDirectory(const std::string &path, bool callbarrier=true)
Create a new directory, renaming the old one if it exists.
Definition AMReX_Utility.cpp:161
void EB_average_down(const MultiFab &S_fine, MultiFab &S_crse, const MultiFab &vol_fine, const MultiFab &vfrac_fine, int scomp, int ncomp, const IntVect &ratio)
Definition AMReX_EBMultiFabUtil.cpp:336
Vector< std::array< T *, AMREX_SPACEDIM > > GetVecOfArrOfPtrs(const Vector< std::array< std::unique_ptr< T >, AMREX_SPACEDIM > > &a)
Definition AMReX_Vector.H:138
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
BottomSolver
Definition AMReX_MLLinOp.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:111
void average_face_to_cellcenter(MultiFab &cc, int dcomp, const Vector< const MultiFab * > &fc, int ngrow)
Average face-based MultiFab onto cell-centered MultiFab.
Definition AMReX_MultiFabUtil.cpp:141
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T norm(const GpuComplex< T > &a_z) noexcept
Return the norm (magnitude squared) of a complex number.
Definition AMReX_GpuComplex.H:344
void setBndry(MF &dst, typename MF::value_type val, int scomp, int ncomp)
dst = val in ghost cells.
Definition AMReX_FabArrayUtility.H:1815
Vector< T * > GetVecOfPtrs(Vector< T > &a)
Definition AMReX_Vector.H:61
std::array< T *, AMREX_SPACEDIM > GetArrOfPtrs(std::array< T, AMREX_SPACEDIM > &a) noexcept
Definition AMReX_Array.H:852
void Scale(MF &dst, typename MF::value_type val, int scomp, int ncomp, int nghost)
dst *= val
Definition AMReX_FabArrayUtility.H:1822
std::array< T const *, AMREX_SPACEDIM > GetArrOfConstPtrs(const std::array< T, AMREX_SPACEDIM > &a) noexcept
Definition AMReX_Array.H:864
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:225
void ParallelCopy(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &ng_src=IntVect(0), IntVect const &ng_dst=IntVect(0), Periodicity const &period=Periodicity::NonPeriodic())
dst = src w/ MPI communication
Definition AMReX_FabArrayUtility.H:1873
int verbose
Definition AMReX_DistributionMapping.cpp:36
void setVal(MF &dst, typename MF::value_type val)
dst = val
Definition AMReX_FabArrayUtility.H:1808
BoxArray const & boxArray(FabArrayBase const &fa)
std::array< T, N > Array
Definition AMReX_Array.H:24
BCMode
Definition AMReX_MLLinOp.H:85
Location
Definition AMReX_MLLinOp.H:87
FabArray memory allocation information.
Definition AMReX_FabArray.H:66