Register DILU as MPI parallel smoother

OPM · Oct 17, 2023 · e514e43 · e514e43
1 parent bc97a01
commit e514e43
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 47 deletions.
diff --git a/opm/simulators/linalg/DILU.hpp b/opm/simulators/linalg/DILU.hpp
@@ -70,22 +70,21 @@ class SeqDilu : public PreconditionerWithUpdate<X, Y>
     virtual void update() override
     {
         OPM_TIMEBLOCK(update);
-        using Block = typename M::block_type;
 
         auto endi = A_.end();
         for ( auto row_i = A_.begin(); row_i != endi; ++row_i) {
             Dinv_[row_i.index()] = A_[row_i.index()][row_i.index()];
         }
-
+        
         for ( auto row_i = A_.begin(); row_i != endi; ++row_i)
         {
-            Block Dinv_temp = Dinv_[row_i.index()];
+           auto Dinv_temp = Dinv_[row_i.index()];
             for (auto a_ij = row_i->begin(); a_ij.index() < row_i.index(); ++a_ij) {
                 auto a_ji = A_[a_ij.index()].find(row_i.index());
                 // if A[i, j] != 0 and A[j, i] != 0
                 if (a_ji != A_[a_ij.index()].end()) {
                     // Dinv_temp -= A[i, j] * d[j] * A[j, i]
-                    Dinv_temp -= (*a_ij)*Dinv_[a_ij.index()]*(*a_ji);
+                    Dinv_temp -= (*a_ij) * Dune::FieldMatrix(Dinv_[a_ij.index()]) * (*a_ji);
                 }
             }
             Dinv_temp.invert();

diff --git a/opm/simulators/linalg/ExtraSmoothers.hpp b/opm/simulators/linalg/ExtraSmoothers.hpp
@@ -10,32 +10,6 @@ namespace Dune
 
 namespace Amg
 {
-
-    template <class T>
-    class ConstructionTraits;
-
-    template <class T>
-    struct SmootherTraits;
-
-
-    template <class F>
-    struct DiluSmootherArgs : public Dune::Amg::DefaultSmootherArgs<F> {
-        //bool leftPrecond;
-        //DiluSmootherArgs();
-	  //:
-	  //iterations(1),
-	  //relaxationFactor(1.0),
-	  //leftPrecond(true)
-    //   {
-    //   }
-    };
-  template <class M, class X, class Y>
-  struct SmootherTraits< Dune::SeqDilu<M, X, Y>> {
-    //typedef DiluSmootherArgs< Dune::SeqDilu<M, X, Y> > Arguments;
-   typedef DiluSmootherArgs< double > Arguments;
-  };
-
-
     /**
      * @brief Policy for the construction of the SeqDilu smoother
      */
@@ -46,15 +20,13 @@ namespace Amg
         static inline std::shared_ptr<SeqDilu<M, X, Y>> construct(Arguments& args)
         {
             return std::make_shared<SeqDilu<M, X, Y>>(
-							  //args.getMatrix(), args.getArgs().iterations, args.getArgs().relaxationFactor,args.getArgs().leftPrecond);
 							  args.getMatrix());
         }
 
 #else
         static inline SeqDilu<M, X, Y>* construct(Arguments& args)
         {
             return new SeqDilu<M, X, Y>(
-					    //args.getMatrix(), args.getArgs().iterations, args.getArgs().relaxationFactor,args.getArgs().leftPrecond);
 					    args.getMatrix());
         }
 

diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp
@@ -94,20 +94,6 @@ struct AMGSmootherArgsHelper<Opm::ParallelOverlappingILU0<M,V,V,C>>
     }
 };
 
-  template<class M, class V>
-  struct AMGSmootherArgsHelper<Dune::SeqDilu<M,V,V>>
-{
-    static auto args(const PropertyTree& prm)
-    {
-        using Smoother = Dune::SeqDilu<M, V, V>;
-        using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
-        SmootherArgs smootherArgs;
-	//smootherArgs.iterations = prm.get<int>("iterations", 1);
-      //  smootherArgs.relaxationFactor = prm.get<double>("relaxation", 1.0);
-	//smootherArgs.leftPrecond = prm.get<bool>("left_precond", true);       
-        return smootherArgs;
-    }
-};
 
 template <class Operator, class Comm, class Matrix, class Vector>
 typename AMGHelper<Operator, Comm, Matrix, Vector>::Criterion
@@ -208,12 +194,23 @@ struct StandardPreconditioners
         // with the AMG hierarchy construction.
         if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
           F::addCreator("amg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
+            using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<V, V>>;
             const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
             if (smoother == "ILU0" || smoother == "ParOverILU0") {
               using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
               auto crit = AMGHelper<O,C,M,V>::criterion(prm);
               auto sargs = AMGSmootherArgsHelper<Smoother>::args(prm);
-              return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
+              PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
+              return prec;
+            }
+            else if (smoother == "DILU") {
+              using Smoother = Dune::SeqDilu<M, V, V>;
+              using ParSmoother = Dune::BlockPreconditioner<V, V, C, Smoother>;
+              using SmootherArgs = typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments;
+              SmootherArgs sargs;
+              auto crit = AMGHelper<O,C,M,V>::criterion(prm);
+              PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, ParSmoother, C>>(op, crit, sargs, comm);
+              return prec;
             } else {
               OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
             }