parallelbasebackend.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_PARALLEL_BASE_BACKEND_HH
28#define EWOMS_PARALLEL_BASE_BACKEND_HH
29
30#include <dune/common/fvector.hh>
31#include <dune/common/version.hh>
32
33#include <dune/grid/io/file/vtk/vtkwriter.hh>
34
35#include <opm/common/Exceptions.hpp>
36
40
51
52#include <iostream>
53#include <memory>
54#include <sstream>
55
56namespace Opm::Properties {
57
58namespace TTag {
60}
61
64template<class TypeTag>
65struct SparseMatrixAdapter<TypeTag, TTag::ParallelBaseLinearSolver>
66{
67private:
69 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
71
72public:
74};
75
76} // namespace Opm::Properties
77
78namespace Opm {
79namespace Linear {
107template <class TypeTag>
109{
110protected:
112
120
124
126 using SequentialPreconditioner = typename PreconditionerWrapper::SequentialPreconditioner;
127
133
134 enum { dimWorld = GridView::dimensionworld };
135
136public:
138 : simulator_(simulator)
139 , gridSequenceNumber_( -1 )
140 , lastIterations_( -1 )
141 {
142 overlappingMatrix_ = nullptr;
143 overlappingb_ = nullptr;
144 overlappingx_ = nullptr;
145 }
146
148 { cleanup_(); }
149
153 static void registerParameters()
154 {
155 Parameters::Register<Parameters::LinearSolverTolerance<Scalar>>
156 ("The maximum allowed error between of the linear solver");
157 Parameters::Register<Parameters::LinearSolverAbsTolerance<Scalar>>
158 ("The maximum accepted error of the norm of the residual");
159 Parameters::Register<Parameters::LinearSolverOverlapSize>
160 ("The size of the algebraic overlap for the linear solver");
161 Parameters::Register<Parameters::LinearSolverMaxIterations>
162 ("The maximum number of iterations of the linear solver");
163 Parameters::Register<Parameters::LinearSolverVerbosity>
164 ("The verbosity level of the linear solver");
165
166 PreconditionerWrapper::registerParameters();
167 }
168
174 { cleanup_(); }
175
182 void prepare(const SparseMatrixAdapter& M, const Vector& )
183 {
184 // if grid has changed the sequence number has changed too
185 int curSeqNum = simulator_.vanguard().gridSequenceNumber();
186 if (gridSequenceNumber_ == curSeqNum && overlappingMatrix_)
187 // the grid has not changed since the overlappingMatrix_has been created, so
188 // there's noting to do
189 return;
190
191 asImp_().cleanup_();
192 gridSequenceNumber_ = curSeqNum;
193
194 BorderListCreator borderListCreator(simulator_.gridView(),
195 simulator_.model().dofMapper());
196
197 // create the overlapping Jacobian matrix
198 unsigned overlapSize = Parameters::Get<Parameters::LinearSolverOverlapSize>();
199 overlappingMatrix_ = new OverlappingMatrix(M.istlMatrix(),
200 borderListCreator.borderList(),
201 borderListCreator.blackList(),
202 overlapSize);
203
204 // create the overlapping vectors for the residual and the
205 // solution
208
209 // writeOverlapToVTK_();
210 }
211
217 void setResidual(const Vector& b)
218 {
219 // copy the interior values of the non-overlapping residual vector to the
220 // overlapping one
221 overlappingb_->assignAddBorder(b);
222 }
223
229 void getResidual(Vector& b) const
230 {
231 // update the non-overlapping vector with the overlapping one
232 overlappingb_->assignTo(b);
233 }
234
242 {
243 overlappingMatrix_->assignFromNative(M.istlMatrix());
244 overlappingMatrix_->syncAdd();
245 }
246
252 bool solve(Vector& x)
253 {
254 (*overlappingx_) = 0.0;
255
256 auto parPreCond = asImp_().preparePreconditioner_();
257 auto precondCleanupFn = [this]() -> void
258 { this->asImp_().cleanupPreconditioner_(); };
259 auto precondCleanupGuard = Opm::make_guard(precondCleanupFn);
260 // create the parallel scalar product and the parallel operator
261 ParallelScalarProduct parScalarProduct(overlappingMatrix_->overlap());
263
264 // retrieve the linear solver
265 auto solver = asImp_().prepareSolver_(parOperator,
266 parScalarProduct,
267 *parPreCond);
268
269 auto cleanupSolverFn =
270 [this]() -> void
271 { this->asImp_().cleanupSolver_(); };
272 GenericGuard<decltype(cleanupSolverFn)> solverGuard(cleanupSolverFn);
273
274 // run the linear solver and have some fun
275 auto result = asImp_().runSolver_(solver);
276 // store number of iterations used
277 lastIterations_ = result.second;
278
279 // copy the result back to the non-overlapping vector
280 overlappingx_->assignTo(x);
281
282 // return the result of the solver
283 return result.first;
284 }
285
289 size_t iterations () const
290 { return lastIterations_; }
291
292protected:
294 { return *static_cast<Implementation *>(this); }
295
296 const Implementation& asImp_() const
297 { return *static_cast<const Implementation *>(this); }
298
299 void cleanup_()
300 {
301 // create the overlapping Jacobian matrix and vectors
302 delete overlappingMatrix_;
303 delete overlappingb_;
304 delete overlappingx_;
305
307 overlappingb_ = 0;
308 overlappingx_ = 0;
309 }
310
311 std::shared_ptr<ParallelPreconditioner> preparePreconditioner_()
312 {
313 int preconditionerIsReady = 1;
314 try {
315 // update sequential preconditioner
317 }
318 catch (const Dune::Exception& e) {
319 std::cout << "Preconditioner threw exception \"" << e.what()
320 << " on rank " << overlappingMatrix_->overlap().myRank()
321 << "\n" << std::flush;
322 preconditionerIsReady = 0;
323 }
324
325 // make sure that the preconditioner is also ready on all peer
326 // ranks.
327 preconditionerIsReady = simulator_.gridView().comm().min(preconditionerIsReady);
328 if (!preconditionerIsReady)
329 throw NumericalProblem("Creating the preconditioner failed");
330
331 // create the parallel preconditioner
332 return std::make_shared<ParallelPreconditioner>(precWrapper_.get(), overlappingMatrix_->overlap());
333 }
334
336 {
337 precWrapper_.cleanup();
338 }
339
341 {
342 for (int lookedAtRank = 0;
343 lookedAtRank < simulator_.gridView().comm().size(); ++lookedAtRank) {
344 std::cout << "writing overlap for rank " << lookedAtRank << "\n" << std::flush;
345 using VtkField = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
346 int n = simulator_.gridView().size(/*codim=*/dimWorld);
347 VtkField isInOverlap(n);
348 VtkField rankField(n);
349 isInOverlap = 0.0;
350 rankField = 0.0;
351 assert(rankField.two_norm() == 0.0);
352 assert(isInOverlap.two_norm() == 0.0);
353 auto vIt = simulator_.gridView().template begin</*codim=*/dimWorld>();
354 const auto& vEndIt = simulator_.gridView().template end</*codim=*/dimWorld>();
355 const auto& overlap = overlappingMatrix_->overlap();
356 for (; vIt != vEndIt; ++vIt) {
357 int nativeIdx = simulator_.model().vertexMapper().map(*vIt);
358 int localIdx = overlap.foreignOverlap().nativeToLocal(nativeIdx);
359 if (localIdx < 0)
360 continue;
361 rankField[nativeIdx] = simulator_.gridView().comm().rank();
362 if (overlap.peerHasIndex(lookedAtRank, localIdx))
363 isInOverlap[nativeIdx] = 1.0;
364 }
365
366 using VtkWriter = Dune::VTKWriter<GridView>;
367 VtkWriter writer(simulator_.gridView(), Dune::VTK::conforming);
368 writer.addVertexData(isInOverlap, "overlap");
369 writer.addVertexData(rankField, "rank");
370
371 std::ostringstream oss;
372 oss << "overlap_rank=" << lookedAtRank;
373 writer.write(oss.str().c_str(), Dune::VTK::ascii);
374 }
375 }
376
380
384
386};
387}} // namespace Linear, Opm
388
389namespace Opm::Properties {
390
393template<class TypeTag>
394struct LinearSolverScalar<TypeTag, TTag::ParallelBaseLinearSolver>
396
397template<class TypeTag>
398struct OverlappingMatrix<TypeTag, TTag::ParallelBaseLinearSolver>
399{
400private:
401 static constexpr int numEq = getPropValue<TypeTag, Properties::NumEq>();
404 using NonOverlappingMatrix = Dune::BCRSMatrix<MatrixBlock>;
405
406public:
408};
409
410template<class TypeTag>
411struct Overlap<TypeTag, TTag::ParallelBaseLinearSolver>
413
414template<class TypeTag>
415struct OverlappingVector<TypeTag, TTag::ParallelBaseLinearSolver>
416{
417 static constexpr int numEq = getPropValue<TypeTag, Properties::NumEq>();
419 using VectorBlock = Dune::FieldVector<LinearSolverScalar, numEq>;
422};
423
424template<class TypeTag>
425struct OverlappingScalarProduct<TypeTag, TTag::ParallelBaseLinearSolver>
426{
430};
431
432template<class TypeTag>
433struct OverlappingLinearOperator<TypeTag, TTag::ParallelBaseLinearSolver>
434{
439};
440
441template<class TypeTag>
442struct PreconditionerWrapper<TypeTag, TTag::ParallelBaseLinearSolver>
444
445} // namespace Opm::Properties
446
447#endif
A simple class which makes sure that a cleanup function is called once the object is destroyed.
Definition: genericguard.hh:43
A sparse matrix interface backend for BCRSMatrix from dune-istl.
Definition: istlsparsematrixadapter.hh:43
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:54
An overlap aware block vector.
Definition: overlappingblockvector.hh:50
An overlap aware linear operator usable by ISTL.
Definition: overlappingoperator.hh:42
An overlap aware preconditioner for any ISTL linear solver.
Definition: overlappingpreconditioner.hh:48
An overlap aware ISTL scalar product.
Definition: overlappingscalarproduct.hh:42
Provides the common code which is required by most linear solvers.
Definition: parallelbasebackend.hh:109
void getResidual(Vector &b) const
Retrieve the synchronized internal residual vector.
Definition: parallelbasebackend.hh:229
GetPropType< TypeTag, Properties::OverlappingVector > OverlappingVector
Definition: parallelbasebackend.hh:122
int gridSequenceNumber_
Definition: parallelbasebackend.hh:378
size_t iterations() const
Return number of iterations used during last solve.
Definition: parallelbasebackend.hh:289
ParallelBaseBackend(const Simulator &simulator)
Definition: parallelbasebackend.hh:137
const Implementation & asImp_() const
Definition: parallelbasebackend.hh:296
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: parallelbasebackend.hh:114
PreconditionerWrapper precWrapper_
Definition: parallelbasebackend.hh:385
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: parallelbasebackend.hh:113
GetPropType< TypeTag, Properties::GlobalEqVector > Vector
Definition: parallelbasebackend.hh:117
OverlappingVector * overlappingx_
Definition: parallelbasebackend.hh:383
typename PreconditionerWrapper::SequentialPreconditioner SequentialPreconditioner
Definition: parallelbasebackend.hh:126
void writeOverlapToVTK_()
Definition: parallelbasebackend.hh:340
std::shared_ptr< ParallelPreconditioner > preparePreconditioner_()
Definition: parallelbasebackend.hh:311
GetPropType< TypeTag, Properties::SparseMatrixAdapter > SparseMatrixAdapter
Definition: parallelbasebackend.hh:116
void eraseMatrix()
Causes the solve() method to discared the structure of the linear system of equations the next time i...
Definition: parallelbasebackend.hh:173
void setResidual(const Vector &b)
Assign values to the internal data structure for the residual vector.
Definition: parallelbasebackend.hh:217
Implementation & asImp_()
Definition: parallelbasebackend.hh:293
void prepare(const SparseMatrixAdapter &M, const Vector &)
Set up the internal data structures required for the linear solver.
Definition: parallelbasebackend.hh:182
void setMatrix(const SparseMatrixAdapter &M)
Sets the values of the residual's Jacobian matrix.
Definition: parallelbasebackend.hh:241
void cleanupPreconditioner_()
Definition: parallelbasebackend.hh:335
GetPropType< TypeTag, Properties::LinearSolverScalar > LinearSolverScalar
Definition: parallelbasebackend.hh:115
GetPropType< TypeTag, Properties::LinearSolverBackend > Implementation
Definition: parallelbasebackend.hh:111
static void registerParameters()
Register all run-time parameters for the linear solver.
Definition: parallelbasebackend.hh:153
void cleanup_()
Definition: parallelbasebackend.hh:299
GetPropType< TypeTag, Properties::PreconditionerWrapper > PreconditionerWrapper
Definition: parallelbasebackend.hh:125
size_t lastIterations_
Definition: parallelbasebackend.hh:379
@ dimWorld
Definition: parallelbasebackend.hh:134
OverlappingVector * overlappingb_
Definition: parallelbasebackend.hh:382
OverlappingMatrix * overlappingMatrix_
Definition: parallelbasebackend.hh:381
bool solve(Vector &x)
Actually solve the linear system of equations.
Definition: parallelbasebackend.hh:252
GetPropType< TypeTag, Properties::GridView > GridView
Definition: parallelbasebackend.hh:119
~ParallelBaseBackend()
Definition: parallelbasebackend.hh:147
GetPropType< TypeTag, Properties::Overlap > Overlap
Definition: parallelbasebackend.hh:121
GetPropType< TypeTag, Properties::BorderListCreator > BorderListCreator
Definition: parallelbasebackend.hh:118
GetPropType< TypeTag, Properties::OverlappingMatrix > OverlappingMatrix
Definition: parallelbasebackend.hh:123
const Simulator & simulator_
Definition: parallelbasebackend.hh:377
Definition: istlpreconditionerwrappers.hh:153
Definition: matrixblock.hh:227
Provides wrapper classes for the (non-AMG) preconditioners provided by dune-istl.
Declares the parameters for the black oil model.
Declares the properties required by the black oil model.
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
GenericGuard< Callback > make_guard(Callback &callback)
Definition: genericguard.hh:88
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
GetPropType< TypeTag, Properties::Scalar > type
Definition: parallelbasebackend.hh:395
The floating point type used internally by the linear solver.
Definition: linalgproperties.hh:46
typename GetPropType< TypeTag, Properties::OverlappingMatrix >::Overlap type
Definition: parallelbasebackend.hh:412
Definition: linalgproperties.hh:60
GetPropType< TypeTag, Properties::OverlappingVector > OverlappingVector
Definition: parallelbasebackend.hh:436
GetPropType< TypeTag, Properties::OverlappingMatrix > OverlappingMatrix
Definition: parallelbasebackend.hh:435
Definition: linalgproperties.hh:63
Definition: linalgproperties.hh:66
GetPropType< TypeTag, Properties::Overlap > Overlap
Definition: parallelbasebackend.hh:428
GetPropType< TypeTag, Properties::OverlappingVector > OverlappingVector
Definition: parallelbasebackend.hh:427
Definition: linalgproperties.hh:69
GetPropType< TypeTag, Properties::LinearSolverScalar > LinearSolverScalar
Definition: parallelbasebackend.hh:418
Dune::FieldVector< LinearSolverScalar, numEq > VectorBlock
Definition: parallelbasebackend.hh:419
GetPropType< TypeTag, Properties::Overlap > Overlap
Definition: parallelbasebackend.hh:420
Definition: linalgproperties.hh:72
the preconditioner used by the linear solver
Definition: linalgproperties.hh:42
typename Opm::Linear::IstlSparseMatrixAdapter< Block > type
Definition: parallelbasebackend.hh:73
The class that allows to manipulate sparse matrices.
Definition: linalgproperties.hh:50
Definition: parallelbasebackend.hh:59