// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey
 *
 * This program is free software, distributed under the terms of the
 * Apache License, Version 2.0
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Conjugate gradient solver, for pressure and viscosity
 *
 ******************************************************************************/

#include "conjugategrad.h"
#include "commonkernels.h"

using namespace std;
namespace Manta {

const int CG_DEBUGLEVEL = 3;

//*****************************************************************************
//  Precondition helpers

//! Preconditioning a la Wavelet Turbulence (needs 4 add. grids)
void InitPreconditionIncompCholesky(const FlagGrid &flags,
                                    Grid<Real> &A0,
                                    Grid<Real> &Ai,
                                    Grid<Real> &Aj,
                                    Grid<Real> &Ak,
                                    Grid<Real> &orgA0,
                                    Grid<Real> &orgAi,
                                    Grid<Real> &orgAj,
                                    Grid<Real> &orgAk)
{
  // compute IC according to Golub and Van Loan
  A0.copyFrom(orgA0);
  Ai.copyFrom(orgAi);
  Aj.copyFrom(orgAj);
  Ak.copyFrom(orgAk);

  FOR_IJK(A0)
  {
    if (flags.isFluid(i, j, k)) {
      const IndexInt idx = A0.index(i, j, k);
      A0[idx] = sqrt(A0[idx]);

      // correct left and top stencil in other entries
      // for i = k+1:n
      //  if (A(i,k) != 0)
      //    A(i,k) = A(i,k) / A(k,k)
      Real invDiagonal = 1.0f / A0[idx];
      Ai[idx] *= invDiagonal;
      Aj[idx] *= invDiagonal;
      Ak[idx] *= invDiagonal;

      // correct the right and bottom stencil in other entries
      // for j = k+1:n
      //   for i = j:n
      //      if (A(i,j) != 0)
      //        A(i,j) = A(i,j) - A(i,k) * A(j,k)
      A0(i + 1, j, k) -= square(Ai[idx]);
      A0(i, j + 1, k) -= square(Aj[idx]);
      A0(i, j, k + 1) -= square(Ak[idx]);
    }
  }

  // invert A0 for faster computation later
  InvertCheckFluid(flags, A0);
};

//! Preconditioning using modified IC ala Bridson (needs 1 add. grid)
void InitPreconditionModifiedIncompCholesky2(const FlagGrid &flags,
                                             Grid<Real> &Aprecond,
                                             Grid<Real> &A0,
                                             Grid<Real> &Ai,
                                             Grid<Real> &Aj,
                                             Grid<Real> &Ak)
{
  // compute IC according to Golub and Van Loan
  Aprecond.clear();

  FOR_IJK(flags)
  {
    if (!flags.isFluid(i, j, k))
      continue;

    const Real tau = 0.97;
    const Real sigma = 0.25;

    // compute modified incomplete cholesky
    Real e = 0.;
    e = A0(i, j, k) - square(Ai(i - 1, j, k) * Aprecond(i - 1, j, k)) -
        square(Aj(i, j - 1, k) * Aprecond(i, j - 1, k)) -
        square(Ak(i, j, k - 1) * Aprecond(i, j, k - 1));
    e -= tau *
         (Ai(i - 1, j, k) * (Aj(i - 1, j, k) + Ak(i - 1, j, k)) * square(Aprecond(i - 1, j, k)) +
          Aj(i, j - 1, k) * (Ai(i, j - 1, k) + Ak(i, j - 1, k)) * square(Aprecond(i, j - 1, k)) +
          Ak(i, j, k - 1) * (Ai(i, j, k - 1) + Aj(i, j, k - 1)) * square(Aprecond(i, j, k - 1)) +
          0.);

    // stability cutoff
    if (e < sigma * A0(i, j, k))
      e = A0(i, j, k);

    Aprecond(i, j, k) = 1. / sqrt(e);
  }
};

//! Preconditioning using multigrid ala Dick et al.
void InitPreconditionMultigrid(
    GridMg *MG, Grid<Real> &A0, Grid<Real> &Ai, Grid<Real> &Aj, Grid<Real> &Ak, Real mAccuracy)
{
  // build multigrid hierarchy if necessary
  if (!MG->isASet())
    MG->setA(&A0, &Ai, &Aj, &Ak);
  MG->setCoarsestLevelAccuracy(mAccuracy * 1E-4);
  MG->setSmoothing(1, 1);
};

//! Apply WT-style ICP
void ApplyPreconditionIncompCholesky(Grid<Real> &dst,
                                     Grid<Real> &Var1,
                                     const FlagGrid &flags,
                                     Grid<Real> &A0,
                                     Grid<Real> &Ai,
                                     Grid<Real> &Aj,
                                     Grid<Real> &Ak,
                                     Grid<Real> &orgA0,
                                     Grid<Real> &orgAi,
                                     Grid<Real> &orgAj,
                                     Grid<Real> &orgAk)
{

  // forward substitution
  FOR_IJK(dst)
  {
    if (!flags.isFluid(i, j, k))
      continue;
    dst(i, j, k) = A0(i, j, k) *
                   (Var1(i, j, k) - dst(i - 1, j, k) * Ai(i - 1, j, k) -
                    dst(i, j - 1, k) * Aj(i, j - 1, k) - dst(i, j, k - 1) * Ak(i, j, k - 1));
  }

  // backward substitution
  FOR_IJK_REVERSE(dst)
  {
    const IndexInt idx = A0.index(i, j, k);
    if (!flags.isFluid(idx))
      continue;
    dst[idx] = A0[idx] * (dst[idx] - dst(i + 1, j, k) * Ai[idx] - dst(i, j + 1, k) * Aj[idx] -
                          dst(i, j, k + 1) * Ak[idx]);
  }
}

//! Apply Bridson-style mICP
void ApplyPreconditionModifiedIncompCholesky2(Grid<Real> &dst,
                                              Grid<Real> &Var1,
                                              const FlagGrid &flags,
                                              Grid<Real> &Aprecond,
                                              Grid<Real> &A0,
                                              Grid<Real> &Ai,
                                              Grid<Real> &Aj,
                                              Grid<Real> &Ak)
{
  // forward substitution
  FOR_IJK(dst)
  {
    if (!flags.isFluid(i, j, k))
      continue;
    const Real p = Aprecond(i, j, k);
    dst(i, j, k) = p *
                   (Var1(i, j, k) - dst(i - 1, j, k) * Ai(i - 1, j, k) * Aprecond(i - 1, j, k) -
                    dst(i, j - 1, k) * Aj(i, j - 1, k) * Aprecond(i, j - 1, k) -
                    dst(i, j, k - 1) * Ak(i, j, k - 1) * Aprecond(i, j, k - 1));
  }

  // backward substitution
  FOR_IJK_REVERSE(dst)
  {
    const IndexInt idx = A0.index(i, j, k);
    if (!flags.isFluid(idx))
      continue;
    const Real p = Aprecond[idx];
    dst[idx] = p * (dst[idx] - dst(i + 1, j, k) * Ai[idx] * p - dst(i, j + 1, k) * Aj[idx] * p -
                    dst(i, j, k + 1) * Ak[idx] * p);
  }
}

//! Perform one Multigrid VCycle
void ApplyPreconditionMultigrid(GridMg *pMG, Grid<Real> &dst, Grid<Real> &Var1)
{
  // one VCycle on "A*dst = Var1" with initial guess dst=0
  pMG->setRhs(Var1);
  pMG->doVCycle(dst);
}

//*****************************************************************************
// Kernels

//! Kernel: Compute the dot product between two Real grids
/*! Uses double precision internally */

struct GridDotProduct : public KernelBase {
  GridDotProduct(const Grid<Real> &a, const Grid<Real> &b)
      : KernelBase(&a, 0), a(a), b(b), result(0.0)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const Grid<Real> &a, const Grid<Real> &b, double &result)
  {
    result += (a[idx] * b[idx]);
  }
  inline operator double()
  {
    return result;
  }
  inline double &getRet()
  {
    return result;
  }
  inline const Grid<Real> &getArg0()
  {
    return a;
  }
  typedef Grid<Real> type0;
  inline const Grid<Real> &getArg1()
  {
    return b;
  }
  typedef Grid<Real> type1;
  void runMessage()
  {
    debMsg("Executing kernel GridDotProduct ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, a, b, result);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  GridDotProduct(GridDotProduct &o, tbb::split) : KernelBase(o), a(o.a), b(o.b), result(0.0)
  {
  }
  void join(const GridDotProduct &o)
  {
    result += o.result;
  }
  const Grid<Real> &a;
  const Grid<Real> &b;
  double result;
};
;

//! Kernel: compute residual (init) and add to sigma

struct InitSigma : public KernelBase {
  InitSigma(const FlagGrid &flags, Grid<Real> &dst, Grid<Real> &rhs, Grid<Real> &temp)
      : KernelBase(&flags, 0), flags(flags), dst(dst), rhs(rhs), temp(temp), sigma(0)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const FlagGrid &flags,
                 Grid<Real> &dst,
                 Grid<Real> &rhs,
                 Grid<Real> &temp,
                 double &sigma)
  {
    const double res = rhs[idx] - temp[idx];
    dst[idx] = (Real)res;

    // only compute residual in fluid region
    if (flags.isFluid(idx))
      sigma += res * res;
  }
  inline operator double()
  {
    return sigma;
  }
  inline double &getRet()
  {
    return sigma;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<Real> &getArg1()
  {
    return dst;
  }
  typedef Grid<Real> type1;
  inline Grid<Real> &getArg2()
  {
    return rhs;
  }
  typedef Grid<Real> type2;
  inline Grid<Real> &getArg3()
  {
    return temp;
  }
  typedef Grid<Real> type3;
  void runMessage()
  {
    debMsg("Executing kernel InitSigma ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, flags, dst, rhs, temp, sigma);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  InitSigma(InitSigma &o, tbb::split)
      : KernelBase(o), flags(o.flags), dst(o.dst), rhs(o.rhs), temp(o.temp), sigma(0)
  {
  }
  void join(const InitSigma &o)
  {
    sigma += o.sigma;
  }
  const FlagGrid &flags;
  Grid<Real> &dst;
  Grid<Real> &rhs;
  Grid<Real> &temp;
  double sigma;
};
;

//! Kernel: update search vector

struct UpdateSearchVec : public KernelBase {
  UpdateSearchVec(Grid<Real> &dst, Grid<Real> &src, Real factor)
      : KernelBase(&dst, 0), dst(dst), src(src), factor(factor)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, Grid<Real> &dst, Grid<Real> &src, Real factor) const
  {
    dst[idx] = src[idx] + factor * dst[idx];
  }
  inline Grid<Real> &getArg0()
  {
    return dst;
  }
  typedef Grid<Real> type0;
  inline Grid<Real> &getArg1()
  {
    return src;
  }
  typedef Grid<Real> type1;
  inline Real &getArg2()
  {
    return factor;
  }
  typedef Real type2;
  void runMessage()
  {
    debMsg("Executing kernel UpdateSearchVec ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, dst, src, factor);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  Grid<Real> &dst;
  Grid<Real> &src;
  Real factor;
};

//*****************************************************************************
//  CG class

template<class APPLYMAT>
GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
                         Grid<Real> &rhs,
                         Grid<Real> &residual,
                         Grid<Real> &search,
                         const FlagGrid &flags,
                         Grid<Real> &tmp,
                         Grid<Real> *pA0,
                         Grid<Real> *pAi,
                         Grid<Real> *pAj,
                         Grid<Real> *pAk)
    : GridCgInterface(),
      mInited(false),
      mIterations(0),
      mDst(dst),
      mRhs(rhs),
      mResidual(residual),
      mSearch(search),
      mFlags(flags),
      mTmp(tmp),
      mpA0(pA0),
      mpAi(pAi),
      mpAj(pAj),
      mpAk(pAk),
      mPcMethod(PC_None),
      mpPCA0(nullptr),
      mpPCAi(nullptr),
      mpPCAj(nullptr),
      mpPCAk(nullptr),
      mMG(nullptr),
      mSigma(0.),
      mAccuracy(VECTOR_EPSILON),
      mResNorm(1e20)
{
}

template<class APPLYMAT> void GridCg<APPLYMAT>::doInit()
{
  mInited = true;
  mIterations = 0;

  mDst.clear();
  mResidual.copyFrom(mRhs);  // p=0, residual = b

  if (mPcMethod == PC_ICP) {
    assertMsg(mDst.is3D(), "ICP only supports 3D grids so far");
    InitPreconditionIncompCholesky(
        mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
    ApplyPreconditionIncompCholesky(
        mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
  }
  else if (mPcMethod == PC_mICP) {
    assertMsg(mDst.is3D(), "mICP only supports 3D grids so far");
    InitPreconditionModifiedIncompCholesky2(mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
    ApplyPreconditionModifiedIncompCholesky2(
        mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
  }
  else if (mPcMethod == PC_MGP) {
    InitPreconditionMultigrid(mMG, *mpA0, *mpAi, *mpAj, *mpAk, mAccuracy);
    ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
  }
  else {
    mTmp.copyFrom(mResidual);
  }

  mSearch.copyFrom(mTmp);

  mSigma = GridDotProduct(mTmp, mResidual);
}

template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate()
{
  if (!mInited)
    doInit();

  mIterations++;

  // create matrix application operator passed as template argument,
  // this could reinterpret the mpA pointers (not so clean right now)
  // tmp = applyMat(search)

  APPLYMAT(mFlags, mTmp, mSearch, *mpA0, *mpAi, *mpAj, *mpAk);

  // alpha = sigma/dot(tmp, search)
  Real dp = GridDotProduct(mTmp, mSearch);
  Real alpha = 0.;
  if (fabs(dp) > 0.)
    alpha = mSigma / (Real)dp;

  gridScaledAdd<Real, Real>(mDst, mSearch, alpha);     // dst += search * alpha
  gridScaledAdd<Real, Real>(mResidual, mTmp, -alpha);  // residual += tmp * -alpha

  if (mPcMethod == PC_ICP)
    ApplyPreconditionIncompCholesky(
        mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
  else if (mPcMethod == PC_mICP)
    ApplyPreconditionModifiedIncompCholesky2(
        mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
  else if (mPcMethod == PC_MGP)
    ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
  else
    mTmp.copyFrom(mResidual);

  // use the l2 norm of the residual for convergence check? (usually max norm is recommended
  // instead)
  if (this->mUseL2Norm) {
    mResNorm = GridSumSqr(mResidual).sum;
  }
  else {
    mResNorm = mResidual.getMaxAbs();
  }

  // abort here to safe some work...
  if (mResNorm < mAccuracy) {
    mSigma = mResNorm;  // this will be returned later on to the caller...
    return false;
  }

  Real sigmaNew = GridDotProduct(mTmp, mResidual);
  Real beta = sigmaNew / mSigma;

  // search =  tmp + beta * search
  UpdateSearchVec(mSearch, mTmp, beta);

  debMsg("GridCg::iterate i=" << mIterations << " sigmaNew=" << sigmaNew << " sigmaLast=" << mSigma
                              << " alpha=" << alpha << " beta=" << beta << " ",
         CG_DEBUGLEVEL);
  mSigma = sigmaNew;

  if (!(mResNorm < 1e35)) {
    if (mPcMethod == PC_MGP) {
      // diverging solves can be caused by the static multigrid mode, we cannot detect this here,
      // though only the pressure solve call "knows" whether the MG is static or dynamics...
      debMsg(
          "GridCg::iterate: Warning - this diverging solve can be caused by the 'static' mode of "
          "the MG preconditioner. If the static mode is active, try switching to dynamic.",
          1);
    }
    errMsg("GridCg::iterate: The CG solver diverged, residual norm > 1e30, stopping.");
  }

  // debMsg("PB-CG-Norms::p"<<sqrt( GridOpNormNosqrt(mpDst, mpFlags).getValue() ) <<"
  // search"<<sqrt( GridOpNormNosqrt(mpSearch, mpFlags).getValue(), CG_DEBUGLEVEL ) <<" res"<<sqrt(
  // GridOpNormNosqrt(mpResidual, mpFlags).getValue() ) <<" tmp"<<sqrt( GridOpNormNosqrt(mpTmp,
  // mpFlags).getValue() ), CG_DEBUGLEVEL ); // debug
  return true;
}

template<class APPLYMAT> void GridCg<APPLYMAT>::solve(int maxIter)
{
  for (int iter = 0; iter < maxIter; iter++) {
    if (!iterate())
      iter = maxIter;
  }
  return;
}

static bool gPrint2dWarning = true;
template<class APPLYMAT>
void GridCg<APPLYMAT>::setICPreconditioner(
    PreconditionType method, Grid<Real> *A0, Grid<Real> *Ai, Grid<Real> *Aj, Grid<Real> *Ak)
{
  assertMsg(method == PC_ICP || method == PC_mICP,
            "GridCg<APPLYMAT>::setICPreconditioner: Invalid method specified.");

  mPcMethod = method;
  if ((!A0->is3D())) {
    if (gPrint2dWarning) {
      debMsg("ICP/mICP pre-conditioning only supported in 3D for now, disabling it.", 1);
      gPrint2dWarning = false;
    }
    mPcMethod = PC_None;
  }
  mpPCA0 = A0;
  mpPCAi = Ai;
  mpPCAj = Aj;
  mpPCAk = Ak;
}

template<class APPLYMAT>
void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg *MG)
{
  assertMsg(method == PC_MGP, "GridCg<APPLYMAT>::setMGPreconditioner: Invalid method specified.");

  mPcMethod = method;

  mMG = MG;
}

// explicit instantiation
template class GridCg<ApplyMatrix>;
template class GridCg<ApplyMatrix2D>;

//*****************************************************************************
// diffusion for real and vec grids, e.g. for viscosity

//! do a CG solve for diffusion; note: diffusion coefficient alpha given in grid space,
//  rescale in python file for discretization independence (or physical correspondence)
//  see lidDrivenCavity.py for an example

void cgSolveDiffusion(const FlagGrid &flags,
                      GridBase &grid,
                      Real alpha = 0.25,
                      Real cgMaxIterFac = 1.0,
                      Real cgAccuracy = 1e-4)
{
  // reserve temp grids
  FluidSolver *parent = flags.getParent();
  Grid<Real> rhs(parent);
  Grid<Real> residual(parent), search(parent), tmp(parent);
  Grid<Real> A0(parent), Ai(parent), Aj(parent), Ak(parent);

  // setup matrix and boundaries
  FlagGrid flagsDummy(parent);
  flagsDummy.setConst(FlagGrid::TypeFluid);
  MakeLaplaceMatrix(flagsDummy, A0, Ai, Aj, Ak);

  FOR_IJK(flags)
  {
    if (flags.isObstacle(i, j, k)) {
      Ai(i, j, k) = Aj(i, j, k) = Ak(i, j, k) = 0.0;
      A0(i, j, k) = 1.0;
    }
    else {
      Ai(i, j, k) *= alpha;
      Aj(i, j, k) *= alpha;
      Ak(i, j, k) *= alpha;
      A0(i, j, k) *= alpha;
      A0(i, j, k) += 1.;
    }
  }

  GridCgInterface *gcg;
  // note , no preconditioning for now...
  const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);

  if (grid.getType() & GridBase::TypeReal) {
    Grid<Real> &u = ((Grid<Real> &)grid);
    rhs.copyFrom(u);
    if (flags.is3D())
      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
    else
      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);

    gcg->setAccuracy(cgAccuracy);
    gcg->solve(maxIter);

    debMsg("FluidSolver::solveDiffusion iterations:" << gcg->getIterations()
                                                     << ", res:" << gcg->getSigma(),
           CG_DEBUGLEVEL);
  }
  else if ((grid.getType() & GridBase::TypeVec3) || (grid.getType() & GridBase::TypeMAC)) {
    Grid<Vec3> &vec = ((Grid<Vec3> &)grid);
    Grid<Real> u(parent);

    // core solve is same as for a regular real grid
    if (flags.is3D())
      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
    else
      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
    gcg->setAccuracy(cgAccuracy);

    // diffuse every component separately
    for (int component = 0; component < (grid.is3D() ? 3 : 2); ++component) {
      getComponent(vec, u, component);
      gcg->forceReinit();

      rhs.copyFrom(u);
      gcg->solve(maxIter);
      debMsg("FluidSolver::solveDiffusion vec3, iterations:" << gcg->getIterations()
                                                             << ", res:" << gcg->getSigma(),
             CG_DEBUGLEVEL);

      setComponent(u, vec, component);
    }
  }
  else {
    errMsg("cgSolveDiffusion: Grid Type is not supported (only Real, Vec3, MAC, or Levelset)");
  }

  delete gcg;
}
static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "cgSolveDiffusion", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      GridBase &grid = *_args.getPtr<GridBase>("grid", 1, &_lock);
      Real alpha = _args.getOpt<Real>("alpha", 2, 0.25, &_lock);
      Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 3, 1.0, &_lock);
      Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 4, 1e-4, &_lock);
      _retval = getPyNone();
      cgSolveDiffusion(flags, grid, alpha, cgMaxIterFac, cgAccuracy);
      _args.check();
    }
    pbFinalizePlugin(parent, "cgSolveDiffusion", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("cgSolveDiffusion", e.what());
    return 0;
  }
}
static const Pb::Register _RP_cgSolveDiffusion("", "cgSolveDiffusion", _W_0);
extern "C" {
void PbRegister_cgSolveDiffusion()
{
  KEEP_UNUSED(_RP_cgSolveDiffusion);
}
}

};  // namespace Manta