From 50d0486b28b45100cc120ea04b7677d094c14668 Mon Sep 17 00:00:00 2001
From: Tong Dong Qiu <tongdong.qiu@bigdataccelerate.com>
Date: Mon, 29 Nov 2021 15:55:56 +0100
Subject: [PATCH] Allow a variable number of pre/post smooth steps during AMG

---
 opm/simulators/linalg/bda/CPR.cpp | 16 +++++++++++-----
 opm/simulators/linalg/bda/CPR.hpp |  2 ++
 2 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/opm/simulators/linalg/bda/CPR.cpp b/opm/simulators/linalg/bda/CPR.cpp
index 5bf52decc..4f0524630 100644
--- a/opm/simulators/linalg/bda/CPR.cpp
+++ b/opm/simulators/linalg/bda/CPR.cpp
@@ -279,6 +279,8 @@ void CPR<block_size>::create_preconditioner_amg(BlockedMatrix *mat_) {
             using CriterionBase = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricDependency<DuneMat, Dune::Amg::FirstDiagonal>>;
             using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
             const Criterion c = Opm::PreconditionerFactory<MatrixOperator, Dune::Amg::SequentialInformation>::amgCriterion(property_tree);
+            num_pre_smooth_steps = c.getNoPreSmoothSteps();
+            num_post_smooth_steps = c.getNoPostSmoothSteps();
 
             dune_amg->build<OverlapFlags>(c);
 
@@ -451,9 +453,11 @@ void CPR<block_size>::amg_cycle_gpu(const int level, cl::Buffer &y, cl::Buffer &
     cl::Buffer& u = d_u[level]; // u was 0-initialized earlier
 
     // presmooth
-    double jacobi_damping = 0.72; // default value in amgcl: 0.72
-    OpenclKernels::residual(A->nnzValues, A->colIndices, A->rowPointers, x, y, t, Ncur, 1);
-    OpenclKernels::vmul(jacobi_damping, d_invDiags[level], t, x, Ncur);
+    double jacobi_damping = 0.65; // default value in amgcl: 0.72
+    for (unsigned i = 0; i < num_pre_smooth_steps; ++i){
+        OpenclKernels::residual(A->nnzValues, A->colIndices, A->rowPointers, x, y, t, Ncur, 1);
+        OpenclKernels::vmul(jacobi_damping, d_invDiags[level], t, x, Ncur);
+    }
 
     // move to coarser level
     OpenclKernels::residual(A->nnzValues, A->colIndices, A->rowPointers, x, y, t, Ncur, 1);
@@ -462,8 +466,10 @@ void CPR<block_size>::amg_cycle_gpu(const int level, cl::Buffer &y, cl::Buffer &
     OpenclKernels::prolongate_vector(u, x, d_PcolIndices[level], Ncur);
 
     // postsmooth
-    OpenclKernels::residual(A->nnzValues, A->colIndices, A->rowPointers, x, y, t, Ncur, 1);
-    OpenclKernels::vmul(jacobi_damping, d_invDiags[level], t, x, Ncur);
+    for (unsigned i = 0; i < num_post_smooth_steps; ++i){
+        OpenclKernels::residual(A->nnzValues, A->colIndices, A->rowPointers, x, y, t, Ncur, 1);
+        OpenclKernels::vmul(jacobi_damping, d_invDiags[level], t, x, Ncur);
+    }
 }
 
 
diff --git a/opm/simulators/linalg/bda/CPR.hpp b/opm/simulators/linalg/bda/CPR.hpp
index ac8157055..1d3160bcf 100644
--- a/opm/simulators/linalg/bda/CPR.hpp
+++ b/opm/simulators/linalg/bda/CPR.hpp
@@ -89,6 +89,8 @@ private:
     bool always_recalculate_aggregates = false; // OPM always reuses the aggregates by default
     bool recalculate_aggregates = true;         // only rerecalculate if true
     const int pressure_idx = 1;                 // hardcoded to mimic OPM
+    unsigned num_pre_smooth_steps;              // number of Jacobi smooth steps before restriction
+    unsigned num_post_smooth_steps;             // number of Jacobi smooth steps after prolongation
 
     std::unique_ptr<openclSolverBackend<1> > coarse_solver; // coarse solver is scalar
     ILUReorder opencl_ilu_reorder;                          // reordering strategy for ILU0 in coarse solver