CD_KrylovBiCGStabImplem.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_KRYLOVBICGSTABIMPLEM_H
#define CD_KRYLOVBICGSTABIMPLEM_H
 
// Chombo includes
#include <CH_Timer.H>
#include <parstream.H>
 
// Our includes
#include <CD_KrylovBiCGStab.H>
#include <CD_NamespaceHeader.H>
 
template <class T>
KrylovBiCGStab<T>::~KrylovBiCGStab() noexcept
{
  this->undefine();
}
 
template <class T>
void
KrylovBiCGStab<T>::define(AMRMultigridKrylovOp<T>* a_op, const Vector<LevelData<T>*>& a_template) noexcept
{
  CH_TIME("KrylovBiCGStab::define");
 
  this->undefine();
 
  m_op = a_op;
 
  // Allocate the work vectors once; they are reused across all solves until the next define() (regrid).
  m_op->create(m_r, a_template);
  m_op->create(m_rTilde, a_template);
  m_op->create(m_e, a_template);
  m_op->create(m_p, a_template);
  m_op->create(m_pTilde, a_template);
  m_op->create(m_sTilde, a_template);
  m_op->create(m_t, a_template);
  m_op->create(m_v, a_template);
 
  m_isDefined = true;
}
 
template <class T>
void
KrylovBiCGStab<T>::undefine() noexcept
{
  if (m_isDefined) {
    m_op->clear(m_r);
    m_op->clear(m_rTilde);
    m_op->clear(m_e);
    m_op->clear(m_p);
    m_op->clear(m_pTilde);
    m_op->clear(m_sTilde);
    m_op->clear(m_t);
    m_op->clear(m_v);
 
    m_isDefined = false;
  }
}
 
template <class T>
bool
KrylovBiCGStab<T>::isDefined() const noexcept
{
  return m_isDefined;
}
 
template <class T>
bool
KrylovBiCGStab<T>::solve(Vector<LevelData<T>*>& a_phi, const Vector<LevelData<T>*>& a_rhs) noexcept
{
  CH_TIME("KrylovBiCGStab::solve");
 
  CH_assert(m_isDefined);
 
  AMRMultigridKrylovOp<T>& op = *m_op;
 
  m_exitStatus = -1;
 
  // Initialize the (homogeneous) residual r = rhs - A*phi and the BiCGStab state. The correction is accumulated in
  // m_e and added to a_phi at the very end, so warm starts (a_phi != 0) are honoured.
  auto computeResidual = [&]() -> void {
    op.applyOp(m_r, a_phi, true);        // r = A*phi
    op.axby(m_r, a_rhs, m_r, 1.0, -1.0); // r = rhs - A*phi
  };
 
  computeResidual();
  op.assign(m_rTilde, m_r);
  op.setToZero(m_e);
  op.setToZero(m_p);
  op.setToZero(m_v);
 
  const Real rnorm0 = op.norm(m_r, 2);
 
  if (rnorm0 == 0.0) {
    m_exitStatus   = 1;
    m_residualNorm = 0.0;
    m_iterations   = 0;
    return true;
  }
 
  // Tolerance scale. m_residualScale is the host's zero-residual ||rhs - L(0)||; falling back to rnorm0 when the
  // host did not supply one reproduces a plain reduction relative to this solve's own initial residual.
  const Real scale     = (m_residualScale > 0.0) ? m_residualScale : rnorm0;
  const Real epsScaled = m_eps * scale;
 
  Real rnorm    = rnorm0;
  Real rho      = 1.0;
  Real alpha    = 1.0;
  Real omega    = 1.0;
  int  iter     = 0;
  int  restarts = 0;
 
  // Fold the work-so-far into a_phi and restart the iteration from the current solution. Used on breakdown.
  auto restart = [&]() -> void {
    op.incr(a_phi, m_e, 1.0);
    op.setToZero(m_e);
 
    computeResidual();
 
    op.assign(m_rTilde, m_r);
    op.setToZero(m_p);
    op.setToZero(m_v);
 
    rho   = 1.0;
    alpha = 1.0;
    omega = 1.0;
    rnorm = op.norm(m_r, 2);
 
    restarts++;
  };
 
  // Treat a breakdown: restart if we still can, otherwise stop with status 2.
  auto breakdown = [&]() -> bool {
    if (restarts < m_numRestarts) {
      restart();
 
      return false; // continue
    }
 
    m_exitStatus = 2;
 
    return true; // stop
  };
 
  while (iter < m_maxIter && rnorm > epsScaled) {
    const Real rnormStart = rnorm; // residual at the start of this iteration (for the convergence-rate report)
 
    const Real rhoNew = op.dotProduct(m_rTilde, m_r);
 
    if (rhoNew == 0.0) {
      if (breakdown()) {
        break;
      }
 
      continue;
    }
 
    // p = r + beta*(p - omega*v). With the initial/post-restart state (p = v = 0, rho = alpha = omega = 1) this
    // correctly reduces to p = r, so no special first-iteration branch is needed.
    const Real beta = (rhoNew / rho) * (alpha / omega);
    op.incr(m_p, m_v, -omega);
    op.scale(m_p, beta);
    op.incr(m_p, m_r, 1.0);
 
    op.preCond(m_pTilde, m_p);       // p_tilde = K^{-1} p
    op.applyOp(m_v, m_pTilde, true); // v = A p_tilde
 
    const Real rtv = op.dotProduct(m_rTilde, m_v);
    if (rtv == 0.0) {
      if (breakdown()) {
        break;
      }
 
      continue;
    }
    alpha = rhoNew / rtv;
 
    op.incr(m_e, m_pTilde, alpha); // e += alpha*p_tilde
    op.incr(m_r, m_v, -alpha);     // r = s = r - alpha*v
 
    rnorm = op.norm(m_r, 2);
    if (rnorm <= epsScaled) { // converged on the half-step s
      m_exitStatus = 1;
      iter++; // this (half) iteration was performed -- count it before breaking
 
      break;
    }
 
    op.preCond(m_sTilde, m_r);       // s_tilde = K^{-1} s
    op.applyOp(m_t, m_sTilde, true); // t = A s_tilde
 
    // omega = (t,s)/(t,t) -- both inner products in a single fused MPI reduction.
    Real                        dots[2];
    const Vector<LevelData<T>*> rhsPair[2] = {m_r, m_t};
 
    op.mDotProduct(m_t, 2, rhsPair, dots);
 
    const Real tDotT = dots[1];
 
    if (tDotT == 0.0) {
      if (breakdown()) {
        break;
      }
 
      continue;
    }
 
    omega = dots[0] / tDotT;
 
    op.incr(m_e, m_sTilde, omega); // e += omega*s_tilde
    op.incr(m_r, m_t, -omega);     // r = s - omega*t
 
    rnorm = op.norm(m_r, 2);
 
    rho = rhoNew;
 
    if (m_verbosity >= 3) {
      // Rate is the iteration-over-iteration reduction factor r_{k-1}/r_k (>1 = converging), matching AMRMultiGrid.
      pout() << "      KrylovBiCGStab:: iteration = " << iter << ", rel. residual = " << (rnorm / scale)
             << ", rate = " << ((rnorm > 0.0) ? rnormStart / rnorm : 0.0) << endl;
    }
 
    if (rnorm <= epsScaled) {
      m_exitStatus = 1;
      iter++; // this iteration was performed -- count it before breaking
 
      break;
    }
 
    if (omega == 0.0) {
      if (breakdown()) {
        break;
      }
 
      continue;
    }
 
    iter++;
  }
 
  if (m_exitStatus == -1) {
    m_exitStatus = (rnorm <= epsScaled) ? 1 : 3; // converged on the loop guard, else max iterations
  }
 
  m_residualNorm = rnorm;
  m_iterations   = iter;
 
  op.incr(a_phi, m_e, 1.0); // commit the accumulated correction
 
  return (m_exitStatus == 1);
}
 
#include <CD_NamespaceFooter.H>
 
#endif