My Project
Loading...
Searching...
No Matches
fvbaselinearizer.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*/
28#ifndef EWOMS_FV_BASE_LINEARIZER_HH
29#define EWOMS_FV_BASE_LINEARIZER_HH
30
31#include "fvbaseproperties.hh"
32#include "linearizationtype.hh"
33
34#include <opm/common/Exceptions.hpp>
35#include <opm/common/TimingMacros.hpp>
36#include <opm/grid/utility/SparseTable.hpp>
37
42
43#include <dune/common/version.hh>
44#include <dune/common/fvector.hh>
45#include <dune/common/fmatrix.hh>
46
47#include <type_traits>
48#include <iostream>
49#include <vector>
50#include <thread>
51#include <set>
52#include <exception> // current_exception, rethrow_exception
53#include <mutex>
54
55namespace Opm {
56// forward declarations
57template<class TypeTag>
58class EcfvDiscretization;
59
69template<class TypeTag>
71{
83
91
93
94 using Toolbox = MathToolbox<Evaluation>;
95
96 using Element = typename GridView::template Codim<0>::Entity;
97 using ElementIterator = typename GridView::template Codim<0>::Iterator;
98
99 using Vector = GlobalEqVector;
100
101 using IstlMatrix = typename SparseMatrixAdapter::IstlMatrix;
102
105
106 using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
107 using VectorBlock = Dune::FieldVector<Scalar, numEq>;
108
110
111 // copying the linearizer is not a good idea
114
115public:
117 : jacobian_()
118 {
119 simulatorPtr_ = 0;
120 }
121
123 {
124 auto it = elementCtx_.begin();
125 const auto& endIt = elementCtx_.end();
126 for (; it != endIt; ++it)
127 delete *it;
128 }
129
133 static void registerParameters()
134 { }
135
145 void init(Simulator& simulator)
146 {
147 simulatorPtr_ = &simulator;
148 eraseMatrix();
149 auto it = elementCtx_.begin();
150 const auto& endIt = elementCtx_.end();
151 for (; it != endIt; ++it){
152 delete *it;
153 }
154 elementCtx_.resize(0);
155 fullDomain_ = std::make_unique<FullDomain>(simulator.gridView());
156 }
157
166 {
167 jacobian_.reset();
168 }
169
180 {
183 }
184
196 {
197 linearizeDomain(*fullDomain_);
198 }
199
200 template <class SubDomainType>
202 {
204 // we defer the initialization of the Jacobian matrix until here because the
205 // auxiliary modules usually assume the problem, model and grid to be fully
206 // initialized...
207 if (!jacobian_)
208 initFirstIteration_();
209
210 // Called here because it is no longer called from linearize_().
211 if (static_cast<std::size_t>(domain.view.size(0)) == model_().numTotalDof()) {
212 // We are on the full domain.
213 resetSystem_();
214 } else {
215 resetSystem_(domain);
216 }
217
218 int succeeded;
219 try {
220 linearize_(domain);
221 succeeded = 1;
222 }
223 catch (const std::exception& e)
224 {
225 std::cout << "rank " << simulator_().gridView().comm().rank()
226 << " caught an exception while linearizing:" << e.what()
227 << "\n" << std::flush;
228 succeeded = 0;
229 }
230 catch (...)
231 {
232 std::cout << "rank " << simulator_().gridView().comm().rank()
233 << " caught an exception while linearizing"
234 << "\n" << std::flush;
235 succeeded = 0;
236 }
237 succeeded = simulator_().gridView().comm().min(succeeded);
238
239 if (!succeeded)
240 throw NumericalProblem("A process did not succeed in linearizing the system");
241 }
242
243 void finalize()
244 { jacobian_->finalize(); }
245
251 {
253 // flush possible local caches into matrix structure
254 jacobian_->commit();
255
256 auto& model = model_();
257 const auto& comm = simulator_().gridView().comm();
258 for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
259 bool succeeded = true;
260 try {
261 model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
262 }
263 catch (const std::exception& e) {
264 succeeded = false;
265
266 std::cout << "rank " << simulator_().gridView().comm().rank()
267 << " caught an exception while linearizing:" << e.what()
268 << "\n" << std::flush;
269 }
270
271 succeeded = comm.min(succeeded);
272
273 if (!succeeded)
274 throw NumericalProblem("linearization of an auxiliary equation failed");
275 }
276 }
277
281 const SparseMatrixAdapter& jacobian() const
282 { return *jacobian_; }
283
284 SparseMatrixAdapter& jacobian()
285 { return *jacobian_; }
286
290 const GlobalEqVector& residual() const
291 { return residual_; }
292
293 GlobalEqVector& residual()
294 { return residual_; }
295
296 void setLinearizationType(LinearizationType linearizationType){
297 linearizationType_ = linearizationType;
298 };
299
300 const LinearizationType& getLinearizationType() const{
301 return linearizationType_;
302 };
303
304 void updateDiscretizationParameters()
305 {
306 // This linearizer stores no such parameters.
307 }
308
309 void updateBoundaryConditionData()
310 {
311 // This linearizer stores no such data.
312 }
313
314 void updateFlowsInfo() {
315 // This linearizer stores no such data.
316 }
317
323 const std::map<unsigned, Constraints>& constraintsMap() const
324 { return constraintsMap_; }
325
331 const auto& getFlowsInfo() const
332 {return flowsInfo_;}
333
339 const auto& getFloresInfo() const
340 {return floresInfo_;}
341
342 template <class SubDomainType>
343 void resetSystem_(const SubDomainType& domain)
344 {
345 if (!jacobian_) {
346 initFirstIteration_();
347 }
348
349 // loop over selected elements
350 using GridViewType = decltype(domain.view);
351 ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(domain.view);
352#ifdef _OPENMP
353#pragma omp parallel
354#endif
355 {
356 unsigned threadId = ThreadManager::threadId();
357 auto elemIt = threadedElemIt.beginParallel();
358 MatrixBlock zeroBlock;
359 zeroBlock = 0.0;
360 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
361 const Element& elem = *elemIt;
362 ElementContext& elemCtx = *elementCtx_[threadId];
363 elemCtx.updatePrimaryStencil(elem);
364 // Set to zero the relevant residual and jacobian parts.
365 for (unsigned primaryDofIdx = 0;
366 primaryDofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0);
368 {
369 unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, /*timeIdx=*/0);
370 residual_[globI] = 0.0;
371 jacobian_->clearRow(globI, 0.0);
372 }
373 }
374 }
375 }
376
377private:
378 Simulator& simulator_()
379 { return *simulatorPtr_; }
380 const Simulator& simulator_() const
381 { return *simulatorPtr_; }
382
383 Problem& problem_()
384 { return simulator_().problem(); }
385 const Problem& problem_() const
386 { return simulator_().problem(); }
387
388 Model& model_()
389 { return simulator_().model(); }
390 const Model& model_() const
391 { return simulator_().model(); }
392
393 const GridView& gridView_() const
394 { return problem_().gridView(); }
395
396 const ElementMapper& elementMapper_() const
397 { return model_().elementMapper(); }
398
399 const DofMapper& dofMapper_() const
400 { return model_().dofMapper(); }
401
402 void initFirstIteration_()
403 {
404 // initialize the BCRS matrix for the Jacobian of the residual function
405 createMatrix_();
406
407 // initialize the Jacobian matrix and the vector for the residual function
408 residual_.resize(model_().numTotalDof());
409 resetSystem_();
410
411 // create the per-thread context objects
412 elementCtx_.resize(ThreadManager::maxThreads());
413 for (unsigned threadId = 0; threadId != ThreadManager::maxThreads(); ++ threadId)
414 elementCtx_[threadId] = new ElementContext(simulator_());
415 }
416
417 // Construct the BCRS matrix for the Jacobian of the residual function
418 void createMatrix_()
419 {
420 const auto& model = model_();
421 Stencil stencil(gridView_(), model_().dofMapper());
422
423 // for the main model, find out the global indices of the neighboring degrees of
424 // freedom of each primary degree of freedom
425 sparsityPattern_.clear();
426 sparsityPattern_.resize(model.numTotalDof());
427
428 for (const auto& elem : elements(gridView_())) {
429 stencil.update(elem);
430
431 for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
432 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
433
434 for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
435 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
436 sparsityPattern_[myIdx].insert(neighborIdx);
437 }
438 }
439 }
440
441 // add the additional neighbors and degrees of freedom caused by the auxiliary
442 // equations
443 size_t numAuxMod = model.numAuxiliaryModules();
444 for (unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx)
445 model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern_);
446
447 // allocate raw matrix
448 jacobian_.reset(new SparseMatrixAdapter(simulator_()));
449
450 // create matrix structure based on sparsity pattern
451 jacobian_->reserve(sparsityPattern_);
452 }
453
454 // reset the global linear system of equations.
455 void resetSystem_()
456 {
457 residual_ = 0.0;
458 // zero all matrix entries
459 jacobian_->clear();
460 }
461
462 // query the problem for all constraint degrees of freedom. note that this method is
463 // quite involved and is thus relatively slow.
464 void updateConstraintsMap_()
465 {
466 if (!enableConstraints_())
467 // constraints are not explictly enabled, so we don't need to consider them!
468 return;
469
470 constraintsMap_.clear();
471
472 // loop over all elements...
473 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_());
474#ifdef _OPENMP
475#pragma omp parallel
476#endif
477 {
478 unsigned threadId = ThreadManager::threadId();
479 ElementIterator elemIt = threadedElemIt.beginParallel();
480 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
481 // create an element context (the solution-based quantities are not
482 // available here!)
483 const Element& elem = *elemIt;
484 ElementContext& elemCtx = *elementCtx_[threadId];
485 elemCtx.updateStencil(elem);
486
487 // check if the problem wants to constrain any degree of the current
488 // element's freedom. if yes, add the constraint to the map.
489 for (unsigned primaryDofIdx = 0;
490 primaryDofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0);
491 ++ primaryDofIdx)
492 {
493 Constraints constraints;
494 elemCtx.problem().constraints(constraints,
495 elemCtx,
497 /*timeIdx=*/0);
498 if (constraints.isActive()) {
499 unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, /*timeIdx=*/0);
500 constraintsMap_[globI] = constraints;
501 continue;
502 }
503 }
504 }
505 }
506 }
507
508 // linearize the whole or part of the system
509 template <class SubDomainType>
510 void linearize_(const SubDomainType& domain)
511 {
512 OPM_TIMEBLOCK(linearize_);
513
514 // We do not call resetSystem_() here, since that will set
515 // the full system to zero, not just our part.
516 // Instead, that must be called before starting the linearization.
517
518 // before the first iteration of each time step, we need to update the
519 // constraints. (i.e., we assume that constraints can be time dependent, but they
520 // can't depend on the solution.)
521 if (model_().newtonMethod().numIterations() == 0)
522 updateConstraintsMap_();
523
524 applyConstraintsToSolution_();
525
526 // to avoid a race condition if two threads handle an exception at the same time,
527 // we use an explicit lock to control access to the exception storage object
528 // amongst thread-local handlers
529 std::mutex exceptionLock;
530
531 // storage to any exception that needs to be bridged out of the
532 // parallel block below. initialized to null to indicate no exception
533 std::exception_ptr exceptionPtr = nullptr;
534
535 // relinearize the elements...
536 using GridViewType = decltype(domain.view);
537 ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(domain.view);
538#ifdef _OPENMP
539#pragma omp parallel
540#endif
541 {
542 auto elemIt = threadedElemIt.beginParallel();
543 auto nextElemIt = elemIt;
544 try {
545 for (; !threadedElemIt.isFinished(elemIt); elemIt = nextElemIt) {
546 // give the model and the problem a chance to prefetch the data required
547 // to linearize the next element, but only if we need to consider it
548 nextElemIt = threadedElemIt.increment();
549 if (!threadedElemIt.isFinished(nextElemIt)) {
550 const auto& nextElem = *nextElemIt;
552 || nextElem.partitionType() == Dune::InteriorEntity)
553 {
554 model_().prefetch(nextElem);
555 problem_().prefetch(nextElem);
556 }
557 }
558
559 const auto& elem = *elemIt;
560 if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
561 continue;
562
563 linearizeElement_(elem);
564 }
565 }
566 // If an exception occurs in the parallel block, it won't escape the
567 // block; terminate() is called instead of a handler outside! hence, we
568 // tuck any exceptions that occur away in the pointer. If an exception
569 // occurs in more than one thread at the same time, we must pick one of
570 // them to be rethrown as we cannot have two active exceptions at the
571 // same time. This solution essentially picks one at random. This will
572 // only be a problem if two different kinds of exceptions are thrown, for
573 // instance if one thread experiences a (recoverable) numerical issue
574 // while another is out of memory.
575 catch(...) {
576 std::lock_guard<std::mutex> take(exceptionLock);
577 exceptionPtr = std::current_exception();
578 threadedElemIt.setFinished();
579 }
580 } // parallel block
581
582 // after reduction from the parallel block, exceptionPtr will point to
583 // a valid exception if one occurred in one of the threads; rethrow
584 // it here to let the outer handler take care of it properly
585 if(exceptionPtr) {
586 std::rethrow_exception(exceptionPtr);
587 }
588
589 applyConstraintsToLinearization_();
590 }
591
592
593 // linearize an element in the interior of the process' grid partition
594 template <class ElementType>
595 void linearizeElement_(const ElementType& elem)
596 {
597 unsigned threadId = ThreadManager::threadId();
598
599 ElementContext *elementCtx = elementCtx_[threadId];
600 auto& localLinearizer = model_().localLinearizer(threadId);
601
602 // the actual work of linearization is done by the local linearizer class
603 localLinearizer.linearize(*elementCtx, elem);
604
605 // update the right hand side and the Jacobian matrix
607 globalMatrixMutex_.lock();
608
609 size_t numPrimaryDof = elementCtx->numPrimaryDof(/*timeIdx=*/0);
610 for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
611 unsigned globI = elementCtx->globalSpaceIndex(/*spaceIdx=*/primaryDofIdx, /*timeIdx=*/0);
612
613 // update the right hand side
614 residual_[globI] += localLinearizer.residual(primaryDofIdx);
615
616 // update the global Jacobian matrix
617 for (unsigned dofIdx = 0; dofIdx < elementCtx->numDof(/*timeIdx=*/0); ++ dofIdx) {
618 unsigned globJ = elementCtx->globalSpaceIndex(/*spaceIdx=*/dofIdx, /*timeIdx=*/0);
619
620 jacobian_->addToBlock(globJ, globI, localLinearizer.jacobian(dofIdx, primaryDofIdx));
621 }
622 }
623
625 globalMatrixMutex_.unlock();
626 }
627
628 // apply the constraints to the solution. (i.e., the solution of constraint degrees
629 // of freedom is set to the value of the constraint.)
630 void applyConstraintsToSolution_()
631 {
632 if (!enableConstraints_())
633 return;
634
635 // TODO: assuming a history size of 2 only works for Euler time discretizations!
636 auto& sol = model_().solution(/*timeIdx=*/0);
637 auto& oldSol = model_().solution(/*timeIdx=*/1);
638
639 auto it = constraintsMap_.begin();
640 const auto& endIt = constraintsMap_.end();
641 for (; it != endIt; ++it) {
642 sol[it->first] = it->second;
643 oldSol[it->first] = it->second;
644 }
645 }
646
647 // apply the constraints to the linearization. (i.e., for constrain degrees of
648 // freedom the Jacobian matrix maps to identity and the residual is zero)
649 void applyConstraintsToLinearization_()
650 {
651 if (!enableConstraints_())
652 return;
653
654 auto it = constraintsMap_.begin();
655 const auto& endIt = constraintsMap_.end();
656 for (; it != endIt; ++it) {
657 unsigned constraintDofIdx = it->first;
658
659 // reset the column of the Jacobian matrix
660 // put an identity matrix on the main diagonal of the Jacobian
661 jacobian_->clearRow(constraintDofIdx, Scalar(1.0));
662
663 // make the right-hand side of constraint DOFs zero
664 residual_[constraintDofIdx] = 0.0;
665 }
666 }
667
668 static bool enableConstraints_()
670
671 Simulator *simulatorPtr_;
672 std::vector<ElementContext*> elementCtx_;
673
674 // The constraint equations (only non-empty if the
675 // EnableConstraints property is true)
676 std::map<unsigned, Constraints> constraintsMap_;
677
678
679 struct FlowInfo
680 {
681 int faceId;
682 VectorBlock flow;
683 unsigned int nncId;
684 };
685 SparseTable<FlowInfo> flowsInfo_;
686 SparseTable<FlowInfo> floresInfo_;
687
688 // the jacobian matrix
689 std::unique_ptr<SparseMatrixAdapter> jacobian_;
690
691 // the right-hand side
692 GlobalEqVector residual_;
693
694 LinearizationType linearizationType_;
695
696 std::mutex globalMatrixMutex_;
697
698 std::vector<std::set<unsigned int>> sparsityPattern_;
699
700 struct FullDomain
701 {
702 explicit FullDomain(const GridView& v) : view (v) {}
703 GridView view;
704 std::vector<bool> interior; // Should remain empty.
705 };
706 // Simple domain object used for full-domain linearization, it allows
707 // us to have the same interface for sub-domain and full-domain work.
708 // Pointer since it must defer construction, due to GridView member.
709 std::unique_ptr<FullDomain> fullDomain_;
710};
711
712} // namespace Opm
713
714#endif
Base class for specifying auxiliary equations.
The common code for the linearizers of non-linear systems of equations.
Definition fvbaselinearizer.hh:71
const auto & getFlowsInfo() const
Return constant reference to the flowsInfo.
Definition fvbaselinearizer.hh:331
const GlobalEqVector & residual() const
Return constant reference to global residual vector.
Definition fvbaselinearizer.hh:290
static void registerParameters()
Register all run-time parameters for the Jacobian linearizer.
Definition fvbaselinearizer.hh:133
void linearize()
Linearize the full system of non-linear equations.
Definition fvbaselinearizer.hh:179
void init(Simulator &simulator)
Initialize the linearizer.
Definition fvbaselinearizer.hh:145
void linearizeDomain()
Linearize the part of the non-linear system of equations that is associated with the spatial domain.
Definition fvbaselinearizer.hh:195
const SparseMatrixAdapter & jacobian() const
Return constant reference to global Jacobian matrix backend.
Definition fvbaselinearizer.hh:281
const auto & getFloresInfo() const
Return constant reference to the floresInfo.
Definition fvbaselinearizer.hh:339
void linearizeAuxiliaryEquations()
Linearize the part of the non-linear system of equations that is associated with the spatial domain.
Definition fvbaselinearizer.hh:250
const std::map< unsigned, Constraints > & constraintsMap() const
Returns the map of constraint degrees of freedom.
Definition fvbaselinearizer.hh:323
void eraseMatrix()
Causes the Jacobian matrix to be recreated from scratch before the next iteration.
Definition fvbaselinearizer.hh:165
Definition matrixblock.hh:227
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition simulator.hh:291
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition simulator.hh:272
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:278
Simplifies multi-threaded capabilities.
Definition threadmanager.hpp:36
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition threadmanager.hpp:66
static unsigned threadId()
Return the index of the current OpenMP thread.
Definition threadmanager.cpp:84
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides data handles for parallel communication which operate on DOFs.
The common code for the linearizers of non-linear systems of equations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
ElementType
The types of reference elements available.
Definition vcfvstencil.hh:52
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
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
Simplifies multi-threaded capabilities.