All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
superlubackend.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef EWOMS_SUPER_LU_BACKEND_HH
28 #define EWOMS_SUPER_LU_BACKEND_HH
29 
30 #if HAVE_SUPERLU
31 
33 
34 #include <opm/common/Unused.hpp>
35 
36 #include <dune/istl/superlu.hh>
37 #include <dune/common/fmatrix.hh>
38 #include <dune/common/version.hh>
39 
40 namespace Ewoms {
41 namespace Properties {
42 // forward declaration of the required property tags
43 NEW_PROP_TAG(Scalar);
44 NEW_PROP_TAG(NumEq);
45 NEW_PROP_TAG(Simulator);
46 NEW_PROP_TAG(JacobianMatrix);
47 NEW_PROP_TAG(GlobalEqVector);
48 NEW_PROP_TAG(LinearSolverVerbosity);
49 NEW_PROP_TAG(LinearSolverBackend);
50 NEW_TYPE_TAG(SuperLULinearSolver);
51 } // namespace Properties
52 } // namespace Ewoms
53 
54 namespace Ewoms {
55 namespace Linear {
56 template <class Scalar, class TypeTag, class Matrix, class Vector>
57 class SuperLUSolve_;
58 
63 template <class TypeTag>
64 class SuperLUBackend
65 {
66  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
67  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
68  typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix;
69  typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
70 
71 public:
72  SuperLUBackend(Simulator& simulator OPM_UNUSED)
73  {}
74 
75  static void registerParameters()
76  {
77  EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity,
78  "The verbosity level of the linear solver");
79  }
80 
87  void eraseMatrix()
88  { }
89 
90  void prepareMatrix(const Matrix& M)
91  {
92  M_ = &M;
93  }
94 
95  void prepareRhs(const Matrix& M OPM_UNUSED, Vector& b)
96  {
97  b_ = &b;
98  }
99 
100  bool solve(Vector& x)
101  { return SuperLUSolve_<Scalar, TypeTag, Matrix, Vector>::solve_(*M_, x, *b_); }
102 
103 private:
104  const Matrix* M_;
105  Vector* b_;
106 };
107 
108 template <class Scalar, class TypeTag, class Matrix, class Vector>
109 class SuperLUSolve_
110 {
111 public:
112  static bool solve_(const Matrix& A, Vector& x, const Vector& b)
113  {
114  Vector bTmp(b);
115 
116  int verbosity = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity);
117  Dune::InverseOperatorResult result;
118  Dune::SuperLU<Matrix> solver(A, verbosity > 0);
119  solver.apply(x, bTmp, result);
120 
121  if (result.converged) {
122  // make sure that the result only contains finite values.
123  Scalar tmp = 0;
124  for (unsigned i = 0; i < x.size(); ++i) {
125  const auto& xi = x[i];
126  for (unsigned j = 0; j < Vector::block_type::dimension; ++j)
127  tmp += xi[j];
128  }
129  result.converged = std::isfinite(tmp);
130  }
131 
132  return result.converged;
133  }
134 };
135 
136 // the following is required to make the SuperLU adapter of dune-istl happy with
137 // quadruple precision math on Dune 2.4. this is because the most which SuperLU can
138 // handle is double precision (i.e., the linear systems of equations are always solved
139 // with at most double precision if chosing SuperLU as the linear solver...)
140 #if HAVE_QUAD
141 template <class TypeTag, class Matrix, class Vector>
142 class SuperLUSolve_<__float128, TypeTag, Matrix, Vector>
143 {
144 public:
145  static bool solve_(const Matrix& A,
146  Vector& x,
147  const Vector& b)
148  {
149  static const int numEq = GET_PROP_VALUE(TypeTag, NumEq);
150  typedef Dune::FieldVector<double, numEq> DoubleEqVector;
151  typedef Dune::FieldMatrix<double, numEq, numEq> DoubleEqMatrix;
152  typedef Dune::BlockVector<DoubleEqVector> DoubleVector;
153  typedef Dune::BCRSMatrix<DoubleEqMatrix> DoubleMatrix;
154 
155  // copy the inputs into the double precision data structures
156  DoubleVector bDouble(b);
157  DoubleVector xDouble(x);
158  DoubleMatrix ADouble(A);
159 
160  bool res =
161  SuperLUSolve_<double, TypeTag, Matrix, Vector>::solve_(ADouble,
162  xDouble,
163  bDouble);
164 
165  // copy the result back into the quadruple precision vector.
166  x = xDouble;
167 
168  return res;
169  }
170 };
171 #endif
172 
173 } // namespace Linear
174 } // namespace Ewoms
175 
176 namespace Ewoms {
177 namespace Properties {
178 SET_INT_PROP(SuperLULinearSolver, LinearSolverVerbosity, 0);
179 SET_TYPE_PROP(SuperLULinearSolver, LinearSolverBackend,
180  Ewoms::Linear::SuperLUBackend<TypeTag>);
181 } // namespace Properties
182 } // namespace Ewoms
183 
184 #endif // HAVE_SUPERLU
185 
186 #endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
This file provides the infrastructure to retrieve run-time parameters.
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377