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