ISTL_Solver.hpp 3.33 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#pragma once

#include <dune/istl/solvers.hh>
#include <dune/istl/umfpack.hh>

namespace AMDiS
{ 
  // creators for solvers, depending on matrix-type
  // ---------------------------------------------------------------------------

  template <class Matrix, class VectorX, class VectorB = VectorX>
  class ISTLSolverCreator;
  
  template <class... MatrixRows, class VectorX, class VectorB>
  class ISTLSolverCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
  {
    using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
    
  public:
    ISTLSolverCreator(std::string prefix)
      : prefix(prefix)
      , info(2)
    {
      Parameters::get(prefix + "->info", info);
    }
    
    template <class Matrix, class Precon>
    std::unique_ptr<BaseSolver> create(std::string name, Matrix& A, Precon& P)
    {
      size_t maxIter = 500;
      Parameters::get(prefix + "->max iteration", maxIter);
      
      double rTol = 1.e-6;
      Parameters::get(prefix + "->reduction", rTol);
      
      if (name == "cg")
        return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
      else if (name == "bicgstab")
        return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
      else if (name == "minres")
        return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
      else if (name == "gmres") 
      {
        int restart = 30;
        Parameters::get(prefix + "->restart", restart);
        return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info); 
      }
      else
        AMDIS_ERROR_EXIT("Unknown solver name!");
    }
    
  private:
    std::string prefix;
    int         info;
  };
  
  
  template <class Block, class Alloc, class VectorX, class VectorB>
  class ISTLSolverCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
  {
    using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
    using Matrix     = Dune::BCRSMatrix<Block,Alloc>;
    
  public:
    ISTLSolverCreator(std::string prefix)
      : prefix(prefix)
      , info(2)
    {
      Parameters::get(prefix + "->info", info);
    }
    
    template <class LinOperator, class Precon>
    std::unique_ptr<BaseSolver> create(std::string name, LinOperator& A, Precon& P)
    {
      size_t maxIter = 500;
      Parameters::get(prefix + "->max iteration", maxIter);
      
      double rTol = 1.e-6;
      Parameters::get(prefix + "->relative tolerance", rTol);
      
      if (name == "cg")
        return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
      else if (name == "bicgstab")
        return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
      else if (name == "minres")
        return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
      else if (name == "gmres") 
      {
        int restart = 30;
        Parameters::get(prefix + "->restart", restart);
        return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info); 
      }
#if HAVE_SUITESPARSE_UMFPACK
      else if (name == "umfpack" || name == "direct")
        return std::make_unique<Dune::UMFPack<Matrix>>(A, info);
#endif
      else
        AMDIS_ERROR_EXIT("Unknown solver name!");
    }
    
  private:
    std::string prefix;
    int         info;
  };

} // end namespace AMDiS