My Project
Loading...
Searching...
No Matches
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
45#include <opm/simulators/linalg/matrixblock.hh>
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:
71
72public:
74};
75
76} // namespace Opm::Properties
77
78namespace Opm {
79namespace Linear {
107template <class TypeTag>
108class ParallelBaseBackend
109{
110protected:
112
120
124
126 using SequentialPreconditioner = typename PreconditionerWrapper::SequentialPreconditioner;
127
130 using ParallelOperator = Opm::Linear::OverlappingOperator<OverlappingMatrix,
131 OverlappingVector,
132 OverlappingVector>;
133
134 enum { dimWorld = GridView::dimensionworld };
135
136public:
137 explicit ParallelBaseBackend(const Simulator& simulator)
138 : simulator_(simulator)
139 , gridSequenceNumber_( -1 )
140 , lastIterations_( -1 )
141 {
142 overlappingMatrix_ = nullptr;
143 overlappingb_ = nullptr;
144 overlappingx_ = nullptr;
145 }
146
147 ~ParallelBaseBackend()
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
173 void eraseMatrix()
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
206 overlappingb_ = new OverlappingVector(overlappingMatrix_->overlap());
207 overlappingx_ = new OverlappingVector(*overlappingb_);
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
241 void setMatrix(const SparseMatrixAdapter& M)
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());
262 ParallelOperator parOperator(*overlappingMatrix_);
263
264 // retrieve the linear solver
265 auto solver = asImp_().prepareSolver_(parOperator,
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:
293 Implementation& asImp_()
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
306 overlappingMatrix_ = 0;
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
316 precWrapper_.prepare(*overlappingMatrix_);
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;
323 }
324
325 // make sure that the preconditioner is also ready on all peer
326 // ranks.
327 preconditionerIsReady = simulator_.gridView().comm().min(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
335 void cleanupPreconditioner_()
336 {
337 precWrapper_.cleanup();
338 }
339
340 void writeOverlapToVTK_()
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);
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
377 const Simulator& simulator_;
378 int gridSequenceNumber_;
379 size_t lastIterations_;
380
381 OverlappingMatrix *overlappingMatrix_;
382 OverlappingVector *overlappingb_;
383 OverlappingVector *overlappingx_;
384
385 PreconditionerWrapper precWrapper_;
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{
437 using type = Opm::Linear::OverlappingOperator<OverlappingMatrix, OverlappingVector,
438 OverlappingVector>;
439};
440
441template<class TypeTag>
442struct PreconditionerWrapper<TypeTag, TTag::ParallelBaseLinearSolver>
444
445} // namespace Opm::Properties
446
447#endif
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
Definition istlpreconditionerwrappers.hh:153
Definition matrixblock.hh:227
Manages the initializing and running of time dependent problems.
Definition simulator.hh:97
A simple class which makes sure that a cleanup function is called once the object is destroyed.
Provides wrapper classes for the (non-AMG) preconditioners provided by dune-istl.
A sparse matrix interface backend for BCRSMatrix from dune-istl.
Declares the parameters for the black oil model.
Declares the properties required by the black oil model.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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
An overlap aware block-compressed row storage (BCRS) matrix.
An overlap aware block vector.
An overlap aware linear operator usable by ISTL.
An overlap aware preconditioner for any ISTL linear solver.
An overlap aware ISTL scalar product.
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
The floating point type used internally by the linear solver.
Definition linalgproperties.hh:46
Definition linalgproperties.hh:60
Definition linalgproperties.hh:63
Definition linalgproperties.hh:66
Definition linalgproperties.hh:69
Definition linalgproperties.hh:72
the preconditioner used by the linear solver
Definition linalgproperties.hh:42
The class that allows to manipulate sparse matrices.
Definition linalgproperties.hh:50
Definition parallelbasebackend.hh:59