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
32#include <dune/common/fmatrix.hh>
33#include <dune/common/version.hh>
34#include <dune/istl/superlu.hh>
35
36#include <opm/material/common/Unused.hpp>
37
38#include <opm/models/linear/istlsparsematrixbackend.hh>
39
41
43
44namespace Opm::Properties::TTag {
45
46struct SuperLULinearSolver {};
47
48} // namespace Opm::Properties::TTag
49
50namespace Opm::Linear {
51
52template <class Scalar, class TypeTag, class Matrix, class Vector>
53class SuperLUSolve_;
54
59template <class TypeTag>
60class SuperLUBackend
61{
62 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
63 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
64 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
65 using Matrix = typename SparseMatrixAdapter::block_type;
66
67 static_assert(std::is_same<SparseMatrixAdapter, IstlSparseMatrixAdapter<MatrixBlock>::value,
68 "The SuperLU linear solver backend requires the IstlSparseMatrixAdapter");
69
70public:
71 SuperLUBackend(Simulator&)
72 {}
73
74 static void registerParameters()
75 {
76 Parameters::Register<Parameters::LinearSolverVerbosity>
77 ("The verbosity level of the linear solver");
78 }
79
86 void eraseMatrix()
87 { }
88
89 void prepare(const SparseMatrixAdapter& M, const Vector& b)
90 { }
91
92 void setResidual(const Vector& b)
93 { b_ = &b; }
94
95 void getResidual(Vector& b) const
96 { b = *b_; }
97
98 void setMatrix(const SparseMatrixAdapter& M)
99 { M_ = &M; }
100
101 bool solve(Vector& x)
102 { return SuperLUSolve_<Scalar, TypeTag, Matrix, Vector>::solve_(*M_, x, *b_); }
103
104private:
105 const Matrix* M_;
106 Vector* b_;
107};
108
109template <class Scalar, class TypeTag, class Matrix, class Vector>
110class SuperLUSolve_
111{
112public:
113 static bool solve_(const Matrix& A, Vector& x, const Vector& b)
114 {
115 Vector bTmp(b);
116
117 int verbosity = Parameters::Get<Parameters::LinearSolverVerbosity>();
118 Dune::InverseOperatorResult result;
119 Dune::SuperLU<Matrix> solver(A, verbosity > 0);
120 solver.apply(x, bTmp, result);
121
122 if (result.converged) {
123 // make sure that the result only contains finite values.
124 Scalar tmp = 0;
125 for (unsigned i = 0; i < x.size(); ++i) {
126 const auto& xi = x[i];
127 for (unsigned j = 0; j < Vector::block_type::dimension; ++j)
128 tmp += xi[j];
129 }
130 result.converged = std::isfinite(tmp);
131 }
132
133 return result.converged;
134 }
135};
136
137// the following is required to make the SuperLU adapter of dune-istl happy with
138// quadruple precision math on Dune 2.4. this is because the most which SuperLU can
139// handle is double precision (i.e., the linear systems of equations are always solved
140// with at most double precision if chosing SuperLU as the linear solver...)
141#if HAVE_QUAD
142template <class TypeTag, class Matrix, class Vector>
143class SuperLUSolve_<__float128, TypeTag, Matrix, Vector>
144{
145public:
146 static bool solve_(const Matrix& A,
147 Vector& x,
148 const Vector& b)
149 {
150 static const int numEq = getPropValue<TypeTag, Properties::NumEq>();
151 using DoubleEqVector = Dune::FieldVector<double, numEq>;
152 using DoubleEqMatrix = Dune::FieldMatrix<double, numEq, numEq>;
153 using DoubleVector = Dune::BlockVector<DoubleEqVector>;
154 using DoubleMatrix = Dune::BCRSMatrix<DoubleEqMatrix>;
155
156 // copy the inputs into the double precision data structures
157 DoubleVector bDouble(b);
158 DoubleVector xDouble(x);
159 DoubleMatrix ADouble(A);
160
161 bool res =
162 SuperLUSolve_<double, TypeTag, Matrix, Vector>::solve_(ADouble,
163 xDouble,
164 bDouble);
165
166 // copy the result back into the quadruple precision vector.
167 x = xDouble;
168
169 return res;
170 }
171};
172#endif
173
174} // namespace Linear
175} // namespace Opm
176
177namespace Opm::Properties {
178
179template<class TypeTag>
180struct LinearSolverVerbosity<TypeTag, TTag::SuperLULinearSolver>
181{ static constexpr int value = 0; };
182
183template<class TypeTag>
184struct LinearSolverBackend<TypeTag, TTag::SuperLULinearSolver>
185{ using type = Opm::Linear::SuperLUBackend<TypeTag>; };
186
187} // namespace Opm::Properties
188
189#endif // HAVE_SUPERLU
190
191#endif
Declares the properties required by the black oil model.
Definition: bicgstabsolver.hh:42
The generic type tag for problems using the immiscible multi-phase model.
Definition: blackoilmodel.hh:74
Definition: blackoilmodel.hh:72
This file provides the infrastructure to retrieve run-time parameters.
UndefinedProperty type
Definition: linalgproperties.hh:38