CD_EllipticSolverChain.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_ELLIPTICSOLVERCHAIN_H
#define CD_ELLIPTICSOLVERCHAIN_H
 
// Std includes
#include <algorithm>
#include <string>
 
// Chombo includes
#include <AMRMultiGrid.H>
#include <LevelData.H>
#include <Vector.H>
#include <parstream.H>
 
// Our includes
#include <CD_AMRMultigridKrylovOp.H>
#include <CD_KrylovBiCGStab.H>
#include <CD_KrylovGMRES.H>
#include <CD_NamespaceHeader.H>
 
namespace EllipticSolverChain {
 
enum class SolverType
{
  GMG,
 
  GMRES,
 
  BiCGStab
};
 
struct Settings
{
  Vector<SolverType> solvers = Vector<SolverType>(1, SolverType::GMG);
 
  Real eps = 1.E-10;
 
  int maxIter = 50;
 
  int restart = 30;
 
  int numVCycles = 1;
 
  int verbosity = -1;
 
  bool
  usesKrylov() const
  {
    for (int i = 0; i < solvers.size(); i++) {
      if (solvers[i] != SolverType::GMG) {
        return true;
      }
    }
 
    return false;
  }
};
 
struct FallbackPolicy
{
  int retryInterval = 1;
 
  int reprobeIters = 0;
 
  int
  startIndex(const Vector<SolverType>& a_solvers)
  {
    if (!m_latched) {
      return 0;
    }
 
    if (m_reprobeNext || m_sinceProbe + 1 >= std::max(1, retryInterval)) {
      m_sinceProbe  = 0;
      m_reprobeNext = false;
 
      return 0; // re-probe GMG (periodic or adaptive)
    }
 
    m_sinceProbe++; // another solve skipped GMG; advance toward the next periodic re-probe
 
    return this->firstNonGmg(a_solvers);
  }
 
  void
  recordOutcome(const bool         a_gmgAttempted,
                const bool         a_gmgConverged,
                const int          a_krylovIters = -1,
                const int          a_verbosity   = -1,
                const std::string& a_name        = "")
  {
    if (a_gmgAttempted) {
      const bool wasLatched = m_latched;
 
      m_latched = !a_gmgConverged;
      if (m_latched) {
        m_sinceProbe = 0;
      }
 
      if (a_verbosity >= 1 && wasLatched && !m_latched) {
        pout() << a_name << " - multigrid re-probe succeeded; releasing the Krylov fallback latch" << endl;
      }
    }
    else if (m_latched && reprobeIters > 0 && a_krylovIters >= 0 && a_krylovIters <= reprobeIters) {
      // The Krylov fallback found the problem easy, so GMG alone may now succeed -- re-probe it next solve.
      m_reprobeNext = true;
    }
  }
 
  bool
  isLatched() const
  {
    return m_latched;
  }
 
protected:
  bool m_latched = false;
 
  int m_sinceProbe = 0;
 
  bool m_reprobeNext = false;
 
  int
  firstNonGmg(const Vector<SolverType>& a_solvers) const
  {
    for (int i = 0; i < a_solvers.size(); i++) {
      if (a_solvers[i] != SolverType::GMG) {
        return i;
      }
    }
 
    return 0;
  }
};
 
inline SolverType
toSolverType(const std::string& a_str)
{
  if (a_str == "gmg") {
    return SolverType::GMG;
  }
  else if (a_str == "gmres") {
    return SolverType::GMRES;
  }
  else if (a_str == "bicgstab") {
    return SolverType::BiCGStab;
  }
  else {
    MayDay::Error("EllipticSolverChain::toSolverType - solver must be 'gmg', 'gmres', or 'bicgstab'");
  }
 
  return SolverType::GMG;
}
 
inline std::string
solverTypeName(const SolverType a_type)
{
  switch (a_type) {
  case SolverType::GMG: {
    return "gmg";
  }
  case SolverType::GMRES: {
    return "gmres";
  }
  case SolverType::BiCGStab: {
    return "bicgstab";
  }
  }
 
  return "unknown";
}
 
template <class T>
class KrylovSolver
{
public:
  using Vec = Vector<LevelData<T>*>;
 
  KrylovSolver() noexcept = default;
 
  ~KrylovSolver() noexcept;
 
  void
  define(AMRMultiGrid<LevelData<T>>*                            a_mg,
         const Vector<RefCountedPtr<LevelData<BaseFab<bool>>>>& a_validCells,
         const Vector<Real>&                                    a_dx,
         const int                                              a_lbase,
         const int                                              a_lmax,
         const int                                              a_numVCycles);
 
  bool
  solve(Vec& a_phi, const Vec& a_rhs, const bool a_zeroPhi, const SolverType a_type, const Settings& a_settings);
 
  int
  lastKrylovIterations() const
  {
    return m_lastKrylovIters;
  }
 
  void
  snapshotGuess(Vec& a_phi, const Vec& a_rhs);
 
  bool
  keepOrRestore(Vec& a_phi, const Vec& a_rhs);
 
protected:
  void
  allocateWork(const Vec& a_template);
 
  void
  freeWork();
 
  AMRMultiGrid<LevelData<T>>* m_mg = nullptr;
 
  AMRMultigridKrylovOp<T> m_op;
 
  KrylovBiCGStab<T> m_bicg;
 
  KrylovGMRES<T> m_gmres;
 
  int m_lastKrylovIters = -1;
 
  int m_lbase = 0;
 
  int m_lmax = 0;
 
  bool m_isDefined = false;
 
  bool m_workAllocated = false;
 
  Vec m_r0;
 
  Vec m_dphi;
 
  Vec m_phiInit;
 
  Real m_initResidual = 0.0;
};
 
} // namespace EllipticSolverChain
 
#include <CD_NamespaceFooter.H>
 
#include <CD_EllipticSolverChainImplem.H>
 
#endif