Force AMG with matrix-add-well-contributions=false, for comparison pu…

…rposes
OPM · Mar 15, 2024 · dddb44e · dddb44e
1 parent 3818c75
commit dddb44e
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 13 deletions.
diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp
@@ -45,7 +45,10 @@ class FlexibleSolver : public Dune::InverseOperator<typename Operator::domain_ty
 {
 public:
     using VectorType = typename Operator::domain_type; // Assuming symmetry: domain == range
+    using MatrixType = typename Operator::matrix_type;
 
+    /// Base class type of the operator passed to the solver.
+    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
     /// Base class type of the contained preconditioner.
     using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
 
@@ -100,6 +103,7 @@ class FlexibleSolver : public Dune::InverseOperator<typename Operator::domain_ty
               std::size_t pressureIndex);
 
     Operator* linearoperator_for_solver_;
+    std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
     std::shared_ptr<AbstractPrecondType> preconditioner_;
     std::shared_ptr<AbstractScalarProductType> scalarproduct_;
     std::shared_ptr<AbstractSolverType> linsolver_;

diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -121,11 +121,25 @@ namespace Dune
         // Parallel case.
         linearoperator_for_solver_ = &op;
         auto child = prm.get_child_optional("preconditioner");
-        preconditioner_ = Opm::PreconditionerFactory<Operator, Comm>::create(op,
-                                                                             child ? *child : Opm::PropertyTree(),
-                                                                             weightsCalculator,
-                                                                             comm,
-                                                                             pressureIndex);
+        const std::string& type = child ? child->get<std::string>("type","ParOverILU0"):std::string("ParOverILU0");
+        if (type != "amg") {
+            preconditioner_ = Opm::PreconditionerFactory<Operator, Comm>::create(op,
+                                                                                 child ? *child : Opm::PropertyTree(),
+                                                                                 weightsCalculator,
+                                                                                 comm,
+                                                                                 pressureIndex);
+        } else {
+            using RangeType = VectorType; // Assuming symmetry: domain == range
+            using ParOperatorType = Opm::GhostLastMatrixAdapter<MatrixType, VectorType, RangeType, Comm>;
+            auto op_prec = std::make_shared<ParOperatorType>(op.getmat(), comm);
+            preconditioner_ = Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*op_prec,
+                                                                                        child ? *child : Opm::PropertyTree(),
+                                                                                        weightsCalculator,
+                                                                                        comm,
+                                                                                        pressureIndex);
+            linearoperator_for_precond_ = op_prec;
+        }
+
         scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
     }
 
@@ -141,10 +155,22 @@ namespace Dune
         // Sequential case.
         linearoperator_for_solver_ = &op;
         auto child = prm.get_child_optional("preconditioner");
-        preconditioner_ = Opm::PreconditionerFactory<Operator,Dune::Amg::SequentialInformation>::create(op,
-                                                                       child ? *child : Opm::PropertyTree(),
-                                                                       weightsCalculator,
-                                                                       pressureIndex);
+        const std::string& type = child ? child->get<std::string>("type","ParOverILU0"):std::string("ParOverILU0");
+        if (type != "amg") {
+            preconditioner_ = Opm::PreconditionerFactory<Operator,Dune::Amg::SequentialInformation>::create(op,
+                                                                           child ? *child : Opm::PropertyTree(),
+                                                                           weightsCalculator,
+                                                                           pressureIndex);
+        } else {
+            using RangeType = VectorType; // Assuming symmetry: domain == range
+            using SOT = Dune::MatrixAdapter<MatrixType, VectorType, RangeType>;
+            auto op_prec = std::make_shared<SOT>(op.getmat());
+            preconditioner_ = Opm::PreconditionerFactory<SOT,Dune::Amg::SequentialInformation>::create(*op_prec,
+                                                                      child ? *child : Opm::PropertyTree(),
+                                                                      weightsCalculator,
+                                                                      pressureIndex);
+            linearoperator_for_precond_ = op_prec;
+        }
         scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
     }
 
@@ -191,7 +217,6 @@ namespace Dune
                                                                                           verbosity);
 #if HAVE_SUITESPARSE_UMFPACK
         } else if (solver_type == "umfpack") {
-            using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
             linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity, false);
             direct_solver_ = true;
 #endif
@@ -223,7 +248,6 @@ namespace Dune
     recreateDirectSolver()
     {
 #if HAVE_SUITESPARSE_UMFPACK
-        using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
         linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), 0, false);
 #else
         OPM_THROW(std::logic_error, "Direct solver specified, but the FlexibleSolver class was not compiled with SuiteSparse support.");
@@ -269,9 +293,11 @@ using SeqOpW = Opm::WellModelMatrixAdapter<OBM<N>, BV<N>, BV<N>, false>;
 // Parallel communicator and operators.
 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
 template <int N>
-using ParOpM = Dune::OverlappingSchwarzOperator<OBM<N>, BV<N>, BV<N>, Comm>;
+using ParOpM = Opm::GhostLastMatrixAdapter<OBM<N>, BV<N>, BV<N>, Comm>;
 template <int N>
 using ParOpW = Opm::WellModelGhostLastMatrixAdapter<OBM<N>, BV<N>, BV<N>, true>;
+template <int N>
+using ParOpD = Dune::OverlappingSchwarzOperator<OBM<N>, BV<N>, BV<N>, Comm>;
 
 // Note: we must instantiate the constructor that is a template.
 // This is only needed in the parallel case, since otherwise the Comm type is
@@ -288,7 +314,8 @@ template Dune::FlexibleSolver<Operator>::FlexibleSolver(Operator& op,
 INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
 INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>); \
 INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<N>); \
-INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<N>);
+INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<N>); \
+INSTANTIATE_FLEXIBLESOLVER_OP(ParOpD<N>);
 
 #else // HAVE_MPI