fvbasediscretization.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_DISCRETIZATION_HH
29#define EWOMS_FV_BASE_DISCRETIZATION_HH
30
31#include <dune/common/fmatrix.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/version.hh>
34#include <dune/istl/bvector.hh>
35
36#include <opm/material/common/MathToolbox.hpp>
37#include <opm/material/common/Valgrind.hpp>
38#include <opm/material/densead/Math.hpp>
39
55
57
60
65
68
69#include <algorithm>
70#include <array>
71#include <cstddef>
72#include <list>
73#include <memory>
74#include <stdexcept>
75#include <sstream>
76#include <string>
77#include <type_traits>
78#include <utility>
79#include <vector>
80
81namespace Opm {
82
83template<class TypeTag>
84class FvBaseDiscretizationNoAdapt;
85
86template<class TypeTag>
87class FvBaseDiscretization;
88
89} // namespace Opm
90
91namespace Opm::Properties {
92
94template<class TypeTag>
95struct Simulator<TypeTag, TTag::FvBaseDiscretization>
97
99template<class TypeTag>
100struct VertexMapper<TypeTag, TTag::FvBaseDiscretization>
101{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
102
104template<class TypeTag>
106{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
107
109template<class TypeTag>
111{
114
115public:
117};
118
119template<class TypeTag>
122
123template<class TypeTag>
126
127template<class TypeTag>
130
132template<class TypeTag>
135
139template<class TypeTag>
140struct EqVector<TypeTag, TTag::FvBaseDiscretization>
141{
142 using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
143 getPropValue<TypeTag, Properties::NumEq>()>;
144};
145
151template<class TypeTag>
152struct RateVector<TypeTag, TTag::FvBaseDiscretization>
154
158template<class TypeTag>
161
165template<class TypeTag>
166struct Constraints<TypeTag, TTag::FvBaseDiscretization>
168
172template<class TypeTag>
174{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
175
179template<class TypeTag>
181{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
182
186template<class TypeTag>
189
193template<class TypeTag>
195{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; };
196
202template<class TypeTag>
205
209template<class TypeTag>
212
213template<class TypeTag>
216
217template<class TypeTag>
220
224template<class TypeTag>
227
228template<class TypeTag>
230{ static constexpr bool value = true; };
231
235template<class TypeTag>
236struct Linearizer<TypeTag, TTag::FvBaseDiscretization>
238
240template<class TypeTag>
242{ static constexpr auto value = Dune::VTK::ascii; };
243
244// disable constraints by default
245template<class TypeTag>
247{ static constexpr bool value = false; };
248
250template<class TypeTag>
252{ static constexpr int value = 2; };
253
256template<class TypeTag>
258{ static constexpr bool value = false; };
259
260// use volumetric residuals is default
261template<class TypeTag>
263{ static constexpr bool value = true; };
264
267template<class TypeTag>
269{ static constexpr bool value = true; };
270
271template <class TypeTag, class MyTypeTag>
273
274#if !HAVE_DUNE_FEM
275template<class TypeTag>
278
279template<class TypeTag>
281{
284};
285#endif
286
287} // namespace Opm::Properties
288
289namespace Opm {
290
296template<class TypeTag>
298{
299 using Implementation = GetPropType<TypeTag, Properties::Model>;
326
329
330 enum {
331 numEq = getPropValue<TypeTag, Properties::NumEq>(),
332 historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
333 };
334
335 using IntensiveQuantitiesVector = std::vector<IntensiveQuantities,
336 aligned_allocator<IntensiveQuantities,
337 alignof(IntensiveQuantities)>>;
338
339 using Element = typename GridView::template Codim<0>::Entity;
340 using ElementIterator = typename GridView::template Codim<0>::Iterator;
341
342 using Toolbox = MathToolbox<Evaluation>;
343 using VectorBlock = Dune::FieldVector<Evaluation, numEq>;
344 using EvalEqVector = Dune::FieldVector<Evaluation, numEq>;
345
346 using LocalEvalBlockVector = typename LocalResidual::LocalEvalBlockVector;
347
348public:
350 {
351 protected:
352 SolutionVector blockVector_;
353 public:
354 BlockVectorWrapper(const std::string&, const std::size_t size)
355 : blockVector_(size)
356 {}
357
359
361 {
362 BlockVectorWrapper result("dummy", 3);
363 result.blockVector_[0] = 1.0;
364 result.blockVector_[1] = 2.0;
365 result.blockVector_[2] = 3.0;
366
367 return result;
368 }
369
370 SolutionVector& blockVector()
371 { return blockVector_; }
372
373 const SolutionVector& blockVector() const
374 { return blockVector_; }
375
376 bool operator==(const BlockVectorWrapper& wrapper) const
377 {
378 return std::equal(this->blockVector_.begin(), this->blockVector_.end(),
379 wrapper.blockVector_.begin(), wrapper.blockVector_.end());
380 }
381
382 template<class Serializer>
383 void serializeOp(Serializer& serializer)
384 {
385 serializer(blockVector_);
386 }
387 };
388
389private:
392
393public:
394 explicit FvBaseDiscretization(Simulator& simulator)
395 : simulator_(simulator)
396 , gridView_(simulator.gridView())
397 , elementMapper_(gridView_, Dune::mcmgElementLayout())
398 , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
399 , newtonMethod_(simulator)
400 , localLinearizer_(ThreadManager::maxThreads())
401 , linearizer_(std::make_unique<Linearizer>())
402 , enableGridAdaptation_(Parameters::Get<Parameters::EnableGridAdaptation>() )
403 , enableIntensiveQuantityCache_(Parameters::Get<Parameters::EnableIntensiveQuantityCache>())
404 , enableStorageCache_(Parameters::Get<Parameters::EnableStorageCache>())
405 , enableThermodynamicHints_(Parameters::Get<Parameters::EnableThermodynamicHints>())
407 {
408 const bool isEcfv = std::is_same_v<Discretization, EcfvDiscretization<TypeTag>>;
409 if (enableGridAdaptation_ && !isEcfv) {
410 throw std::invalid_argument("Grid adaptation currently only works for the "
411 "element-centered finite volume discretization (is: " +
412 Dune::className<Discretization>() + ")");
413 }
414
415 PrimaryVariables::init();
416 // Setting up the intensive quantities cache and storage cache is done in finishInit()
417 // and applyInitialSolution() to ensure that the history size is correct.
418 asImp_().registerOutputModules_();
419 }
420
421 // copying a discretization object is not a good idea
423
427 static void registerParameters()
428 {
429 Linearizer::registerParameters();
430 LocalLinearizer::registerParameters();
431 LocalResidual::registerParameters();
432 GradientCalculator::registerParameters();
433 IntensiveQuantities::registerParameters();
434 ExtensiveQuantities::registerParameters();
436 Linearizer::registerParameters();
437 PrimaryVariables::registerParameters();
438 // register runtime parameters of the output modules
440
441 Parameters::Register<Parameters::EnableGridAdaptation>
442 ("Enable adaptive grid refinement/coarsening");
443 Parameters::Register<Parameters::EnableVtkOutput>
444 ("Global switch for turning on writing VTK files");
445 Parameters::Register<Parameters::EnableThermodynamicHints>
446 ("Enable thermodynamic hints");
447 Parameters::Register<Parameters::EnableIntensiveQuantityCache>
448 ("Turn on caching of intensive quantities");
449 Parameters::Register<Parameters::EnableStorageCache>
450 ("Store previous storage terms and avoid re-calculating them.");
451 Parameters::Register<Parameters::OutputDir>
452 ("The directory to which result files are written");
453 }
454
459 {
460 // initialize the volume of the finite volumes to zero
461 const std::size_t numDof = asImp_().numGridDof();
462 dofTotalVolume_.resize(numDof);
463 std::fill(dofTotalVolume_.begin(), dofTotalVolume_.end(), 0.0);
464
465 ElementContext elemCtx(simulator_);
466 gridTotalVolume_ = 0.0;
467
468 // iterate through the grid and evaluate the initial condition
469 for (const auto& elem : elements(gridView_)) {
470 // ignore everything which is not in the interior if the
471 // current process' piece of the grid
472 if (elem.partitionType() != Dune::InteriorEntity) {
473 continue;
474 }
475
476 // deal with the current element
477 elemCtx.updateStencil(elem);
478 const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
479
480 // loop over all element vertices, i.e. sub control volumes
481 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++) {
482 // map the local degree of freedom index to the global one
483 const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
484
485 const Scalar dofVolume = stencil.subControlVolume(dofIdx).volume();
486 dofTotalVolume_[globalIdx] += dofVolume;
487 gridTotalVolume_ += dofVolume;
488 }
489 }
490
491 // determine which DOFs should be considered to lie fully in the interior of the
492 // local process grid partition: those which do not have a non-zero volume
493 // before taking the peer processes into account...
494 isLocalDof_.resize(numDof);
495 for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
496 isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0);
497 }
498
499 // add the volumes of the DOFs on the process boundaries
500 const auto sumHandle =
501 GridCommHandleFactory::template sumHandle<Scalar>(dofTotalVolume_,
502 asImp_().dofMapper());
503 gridView_.communicate(*sumHandle,
504 Dune::InteriorBorder_All_Interface,
505 Dune::ForwardCommunication);
506
507 // sum up the volumes of the grid partitions
509
510 linearizer_->init(simulator_);
511 for (unsigned threadId = 0; threadId < ThreadManager::maxThreads(); ++threadId) {
512 localLinearizer_[threadId].init(simulator_);
513 }
514
516
517 newtonMethod_.finishInit();
518 }
519
524 { return enableGridAdaptation_; }
525
531 {
532 // first set the whole domain to zero
533 SolutionVector& uCur = asImp_().solution(/*timeIdx=*/0);
534 uCur = Scalar(0.0);
535
536 ElementContext elemCtx(simulator_);
537
538 // iterate through the grid and evaluate the initial condition
539 for (const auto& elem : elements(gridView_)) {
540 // ignore everything which is not in the interior if the
541 // current process' piece of the grid
542 if (elem.partitionType() != Dune::InteriorEntity) {
543 continue;
544 }
545
546 // deal with the current element
547 elemCtx.updateStencil(elem);
548
549 // loop over all element vertices, i.e. sub control volumes
550 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
551 // map the local degree of freedom index to the global one
552 const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
553
554 // let the problem do the dirty work of nailing down
555 // the initial solution.
556 simulator_.problem().initial(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
557 asImp_().supplementInitialSolution_(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
558 uCur[globalIdx].checkDefined();
559 }
560 }
561
562 // synchronize the ghost DOFs (if necessary)
563 asImp_().syncOverlap();
564
565 // also set the solutions of the "previous" time steps to the initial solution.
566 for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx) {
567 solution(timeIdx) = solution(/*timeIdx=*/0);
568 }
569
570 // Initialize intensive quantities cache now that all problem-specific parameters are available.
571 // This ensures intensiveQuantityHistorySize is correct based on recycleFirstIterationStorage().
572 // TODO: Where this is done should perhaps be changed once finishInit() is refactored.
574
575 simulator_.problem().initialSolutionApplied();
576
577#ifndef NDEBUG
578 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
579 const auto& sol = solution(timeIdx);
580 for (unsigned dofIdx = 0; dofIdx < sol.size(); ++dofIdx) {
581 sol[dofIdx].checkDefined();
582 }
583 }
584#endif // NDEBUG
585 }
586
591 void prefetch(const Element&) const
592 {
593 // do nothing by default
594 }
595
599 NewtonMethod& newtonMethod()
600 { return newtonMethod_; }
601
605 const NewtonMethod& newtonMethod() const
606 { return newtonMethod_; }
607
623 const IntensiveQuantities* thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
624 {
626 return 0;
627 }
628
629 // the intensive quantities cache doubles as thermodynamic hint
630 return cachedIntensiveQuantities(globalIdx, timeIdx);
631 }
632
644 const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
645 {
648 !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
649 return nullptr;
650 }
651
652 // With the storage cache enabled, usually only the
653 // intensive quantities for the most recent time step are
654 // cached. However, this may be false for some Problem
655 // variants, so we should check if the cache exists for
656 // the timeIdx in question.
657 if (timeIdx > 0 && enableStorageCache_ && intensiveQuantityCache_[timeIdx].empty()) {
658 return nullptr;
659 }
660
661 return &intensiveQuantityCache_[timeIdx][globalIdx];
662 }
663
672 void updateCachedIntensiveQuantities(const IntensiveQuantities& intQuants,
673 unsigned globalIdx,
674 unsigned timeIdx) const
675 {
677 return;
678 }
679
680 intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
681 intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
682 }
683
693 unsigned timeIdx,
694 bool newValue) const
695 {
697 return;
698 }
699
700 intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue ? 1 : 0;
701 }
702
708
714 void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
715 {
717 return;
718 }
719
721 std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
722 intensiveQuantityCacheUpToDate_[timeIdx].end(),
723 /*value=*/0);
724 }
725 }
726
727 void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
728 {
730
731 // loop over all elements...
732 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
733#ifdef _OPENMP
734#pragma omp parallel
735#endif
736 {
737 ElementContext elemCtx(simulator_);
738 ElementIterator elemIt = threadedElemIt.beginParallel();
739 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
740 const Element& elem = *elemIt;
741 elemCtx.updatePrimaryStencil(elem);
742 elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
743 }
744 }
745 }
746
747 template <class GridViewType>
748 void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
749 {
750 // loop over all elements...
751 ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(gridView);
752#ifdef _OPENMP
753#pragma omp parallel
754#endif
755 {
756
757 ElementContext elemCtx(simulator_);
758 auto elemIt = threadedElemIt.beginParallel();
759 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
760 if (elemIt->partitionType() != Dune::InteriorEntity) {
761 continue;
762 }
763 const Element& elem = *elemIt;
764 elemCtx.updatePrimaryStencil(elem);
765 // Mark cache for this element as invalid.
766 const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
767 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
768 const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
769 setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false);
770 }
771 // Update for this element.
772 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
773 }
774 }
775 }
776
785 void shiftIntensiveQuantityCache(unsigned numSlots = 1)
786 {
787 if (!storeIntensiveQuantities() || numSlots <= 0) {
788 return;
789 }
790
791 if (enableStorageCache() && simulator_.problem().recycleFirstIterationStorage()) {
792 // If the storage term is cached, the intensive quantities of the previous
793 // time steps do not need to be accessed, and we can thus spare ourselves to
794 // copy the objects for the intensive quantities.
795 // However, if the storage term at the start of the timestep cannot be deduced
796 // from the primary variables, we must calculate it from the old intensive
797 // quantities, and need to shift them.
798 return;
799 }
800
801 const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
802 for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize - numSlots; ++timeIdx) {
803 intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
805 }
806
807 // the cache for the most recent time indices do not need to be invalidated
808 // because the solution for them did not change (TODO: that assumes that there is
809 // no post-processing of the solution after a time step! fix it?)
810 }
811
819 { return enableStorageCache_; }
820
829
843 const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
844 {
845 if (!enableStorageCache_ ||
846 timeIdx >= historySize ||
847 !storageCacheUpToDate_[timeIdx][globalIdx]) {
848 throw std::logic_error("Cached storage is not available or up to date for the requested "
849 "global index and time index. Make sure storage cache is enabled "
850 "and the entry is valid before calling this method.");
851 }
852
853 return storageCache_[timeIdx][globalIdx];
854 }
855
867 void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
868 {
869 if (!enableStorageCache_ || timeIdx >= historySize) {
870 return;
871 }
872
873 storageCache_[timeIdx][globalIdx] = value;
874 storageCacheUpToDate_[timeIdx][globalIdx] = 1;
875 }
876
883 bool storageCacheIsUpToDate(unsigned globalIdx, unsigned timeIdx) const
884 {
885 if (!enableStorageCache_ || timeIdx >= historySize) {
886 return false;
887 }
888 return storageCacheUpToDate_[timeIdx][globalIdx] != 0;
889 }
890
897 void invalidateStorageCacheEntry(unsigned globalIdx, unsigned timeIdx) const
898 {
899 if (enableStorageCache_ && timeIdx < historySize) {
900 storageCacheUpToDate_[timeIdx][globalIdx] = 0;
901 }
902 }
903
909 void invalidateStorageCache(unsigned timeIdx) const
910 {
911 if (enableStorageCache_ && timeIdx < historySize) {
912 std::fill(storageCacheUpToDate_[timeIdx].begin(),
913 storageCacheUpToDate_[timeIdx].end(),
914 /*value=*/0);
915 }
916 }
917
929 void shiftStorageCache(unsigned numSlots = 1) const
930 {
931 // If we cannot recycle first iteration storage, it does not make sense to shift the storage cache.
932 if (enableStorageCache_ && !simulator_.problem().recycleFirstIterationStorage()) {
933 for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++timeIdx) {
934 storageCache_[timeIdx + numSlots] = storageCache_[timeIdx];
935 storageCacheUpToDate_[timeIdx + numSlots] = storageCacheUpToDate_[timeIdx];
936 }
937
938 // should we invalidate the cache for the most recent time indices? (see shiftIntensiveQuantityCache)
939 }
940 }
941
949 Scalar globalResidual(GlobalEqVector& dest,
950 const SolutionVector& u) const
951 {
952 mutableSolution(/*timeIdx=*/0) = u;
953 const Scalar res = asImp_().globalResidual(dest);
954 mutableSolution(/*timeIdx=*/0) = asImp_().solution(/*timeIdx=*/0);
955 return res;
956 }
957
964 Scalar globalResidual(GlobalEqVector& dest) const
965 {
966 dest = 0;
967
968 std::mutex mutex;
969 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
970#ifdef _OPENMP
971#pragma omp parallel
972#endif
973 {
974 // Attention: the variables below are thread specific and thus cannot be
975 // moved in front of the #pragma!
976 const unsigned threadId = ThreadManager::threadId();
977 ElementContext elemCtx(simulator_);
978 ElementIterator elemIt = threadedElemIt.beginParallel();
979 LocalEvalBlockVector residual, storageTerm;
980
981 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
982 const Element& elem = *elemIt;
983 if (elem.partitionType() != Dune::InteriorEntity) {
984 continue;
985 }
986
987 elemCtx.updateAll(elem);
988 residual.resize(elemCtx.numDof(/*timeIdx=*/0));
989 storageTerm.resize(elemCtx.numPrimaryDof(/*timeIdx=*/0));
990 asImp_().localResidual(threadId).eval(residual, elemCtx);
991
992 const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
993 mutex.lock();
994 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
995 const unsigned globalI = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
996 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
997 dest[globalI][eqIdx] += Toolbox::value(residual[dofIdx][eqIdx]);
998 }
999 }
1000 mutex.unlock();
1001 }
1002 }
1003
1004 // add up the residuals on the process borders
1005 const auto sumHandle =
1006 GridCommHandleFactory::template sumHandle<EqVector>(dest, asImp_().dofMapper());
1007 gridView_.communicate(*sumHandle,
1008 Dune::InteriorBorder_InteriorBorder_Interface,
1009 Dune::ForwardCommunication);
1010
1011 // calculate the square norm of the residual. this is not
1012 // entirely correct, since the residual for the finite volumes
1013 // which are on the boundary are counted once for every
1014 // process. As often in life: shit happens (, we don't care)...
1015 return std::sqrt(asImp_().gridView().comm().sum(dest.two_norm2()));
1016 }
1017
1024 void globalStorage(EqVector& storage, unsigned timeIdx = 0) const
1025 {
1026 storage = 0;
1027
1028 std::mutex mutex;
1029 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
1030#ifdef _OPENMP
1031#pragma omp parallel
1032#endif
1033 {
1034 // Attention: the variables below are thread specific and thus cannot be
1035 // moved in front of the #pragma!
1036 const unsigned threadId = ThreadManager::threadId();
1037 ElementContext elemCtx(simulator_);
1038 ElementIterator elemIt = threadedElemIt.beginParallel();
1039 LocalEvalBlockVector elemStorage;
1040
1041 // in this method, we need to disable the storage cache because we want to
1042 // evaluate the storage term for other time indices than the most recent one
1043 elemCtx.setEnableStorageCache(false);
1044
1045 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
1046 const Element& elem = *elemIt;
1047 if (elem.partitionType() != Dune::InteriorEntity) {
1048 continue; // ignore ghost and overlap elements
1049 }
1050
1051 elemCtx.updateStencil(elem);
1052 elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
1053
1054 const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
1055 elemStorage.resize(numPrimaryDof);
1056
1057 localResidual(threadId).evalStorage(elemStorage, elemCtx, timeIdx);
1058
1059 mutex.lock();
1060 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
1061 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1062 storage[eqIdx] += Toolbox::value(elemStorage[dofIdx][eqIdx]);
1063 }
1064 }
1065 mutex.unlock();
1066 }
1067 }
1068
1069 storage = gridView_.comm().sum(storage);
1070 }
1071
1079 void checkConservativeness([[maybe_unused]] Scalar tolerance = -1,
1080 [[maybe_unused]] bool verbose = false) const
1081 {
1082#ifndef NDEBUG
1083 Scalar totalBoundaryArea(0.0);
1084 Scalar totalVolume(0.0);
1085 EvalEqVector totalRate(0.0);
1086
1087 // take the newton tolerance times the total volume of the grid if we're not
1088 // given an explicit tolerance...
1089 if (tolerance <= 0) {
1090 tolerance =
1091 simulator_.model().newtonMethod().tolerance() *
1092 simulator_.model().gridTotalVolume() *
1093 1000;
1094 }
1095
1096 // we assume the implicit Euler time discretization for now...
1097 assert(historySize == 2);
1098
1099 EqVector storageBeginTimeStep(0.0);
1100 globalStorage(storageBeginTimeStep, /*timeIdx=*/1);
1101
1102 EqVector storageEndTimeStep(0.0);
1103 globalStorage(storageEndTimeStep, /*timeIdx=*/0);
1104
1105 // calculate the rate at the boundary and the source rate
1106 ElementContext elemCtx(simulator_);
1107 elemCtx.setEnableStorageCache(false);
1108 for (const auto& elem : elements(simulator_.gridView())) {
1109 if (elem.partitionType() != Dune::InteriorEntity) {
1110 continue; // ignore ghost and overlap elements
1111 }
1112
1113 elemCtx.updateAll(elem);
1114
1115 // handle the boundary terms
1116 if (elemCtx.onBoundary()) {
1117 BoundaryContext boundaryCtx(elemCtx);
1118
1119 for (unsigned faceIdx = 0; faceIdx < boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); ++faceIdx) {
1120 BoundaryRateVector values;
1121 simulator_.problem().boundary(values,
1122 boundaryCtx,
1123 faceIdx,
1124 /*timeIdx=*/0);
1125 Valgrind::CheckDefined(values);
1126
1127 const unsigned dofIdx = boundaryCtx.interiorScvIndex(faceIdx, /*timeIdx=*/0);
1128 const auto& insideIntQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1129
1130 const Scalar bfArea =
1131 boundaryCtx.boundarySegmentArea(faceIdx, /*timeIdx=*/0) *
1132 insideIntQuants.extrusionFactor();
1133
1134 for (unsigned i = 0; i < values.size(); ++i) {
1135 values[i] *= bfArea;
1136 }
1137
1138 totalBoundaryArea += bfArea;
1139 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1140 totalRate[eqIdx] += values[eqIdx];
1141 }
1142 }
1143 }
1144
1145 // deal with the source terms
1146 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
1147 RateVector values;
1148 simulator_.problem().source(values,
1149 elemCtx,
1150 dofIdx,
1151 /*timeIdx=*/0);
1152 Valgrind::CheckDefined(values);
1153
1154 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1155 Scalar dofVolume =
1156 elemCtx.dofVolume(dofIdx, /*timeIdx=*/0) *
1157 intQuants.extrusionFactor();
1158 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1159 totalRate[eqIdx] += -dofVolume*Toolbox::value(values[eqIdx]);
1160 }
1161 totalVolume += dofVolume;
1162 }
1163 }
1164
1165 // summarize everything over all processes
1166 const auto& comm = simulator_.gridView().comm();
1167 totalRate = comm.sum(totalRate);
1168 totalBoundaryArea = comm.sum(totalBoundaryArea);
1169 totalVolume = comm.sum(totalVolume);
1170
1171 if (comm.rank() == 0) {
1172 EqVector storageRate = storageBeginTimeStep;
1173 storageRate -= storageEndTimeStep;
1174 storageRate /= simulator_.timeStepSize();
1175 if (verbose) {
1176 std::cout << "storage at beginning of time step: " << storageBeginTimeStep << "\n";
1177 std::cout << "storage at end of time step: " << storageEndTimeStep << "\n";
1178 std::cout << "rate based on storage terms: " << storageRate << "\n";
1179 std::cout << "rate based on source and boundary terms: " << totalRate << "\n";
1180 std::cout << "difference in rates: ";
1181 for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
1182 std::cout << (storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) << " ";
1183 }
1184 std::cout << "\n";
1185 }
1186 for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
1187 Scalar eps =
1188 (std::abs(storageRate[eqIdx]) + Toolbox::value(totalRate[eqIdx])) * tolerance;
1189 eps = std::max(tolerance, eps);
1190 assert(std::abs(storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) <= eps);
1191 }
1192 }
1193#endif // NDEBUG
1194 }
1195
1201 Scalar dofTotalVolume(unsigned globalIdx) const
1202 { return dofTotalVolume_[globalIdx]; }
1203
1209 bool isLocalDof(unsigned globalIdx) const
1210 { return isLocalDof_[globalIdx]; }
1211
1216 Scalar gridTotalVolume() const
1217 { return gridTotalVolume_; }
1218
1224 const SolutionVector& solution(unsigned timeIdx) const
1225 { return solution_[timeIdx]->blockVector(); }
1226
1230 SolutionVector& solution(unsigned timeIdx)
1231 { return solution_[timeIdx]->blockVector(); }
1232
1233 protected:
1237 SolutionVector& mutableSolution(unsigned timeIdx) const
1238 { return solution_[timeIdx]->blockVector(); }
1239
1240 public:
1245 const Linearizer& linearizer() const
1246 { return *linearizer_; }
1247
1252 Linearizer& linearizer()
1253 { return *linearizer_; }
1254
1263 const LocalLinearizer& localLinearizer(unsigned openMpThreadId) const
1264 { return localLinearizer_[openMpThreadId]; }
1265
1269 LocalLinearizer& localLinearizer(unsigned openMpThreadId)
1270 { return localLinearizer_[openMpThreadId]; }
1271
1275 const LocalResidual& localResidual(unsigned openMpThreadId) const
1276 { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1277
1281 LocalResidual& localResidual(unsigned openMpThreadId)
1282 { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1283
1291 Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
1292 {
1293 const Scalar absPv = std::abs(asImp_().solution(/*timeIdx=*/1)[globalDofIdx][pvIdx]);
1294 return 1.0 / std::max(absPv, 1.0);
1295 }
1296
1303 Scalar eqWeight(unsigned, unsigned) const
1304 { return 1.0; }
1305
1315 Scalar relativeDofError(unsigned vertexIdx,
1316 const PrimaryVariables& pv1,
1317 const PrimaryVariables& pv2) const
1318 {
1319 Scalar result = 0.0;
1320 for (unsigned j = 0; j < numEq; ++j) {
1321 const Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
1322 const Scalar eqErr = std::abs((pv1[j] - pv2[j])*weight);
1323 //Scalar eqErr = std::abs(pv1[j] - pv2[j]);
1324 //eqErr *= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
1325
1326 result = std::max(result, eqErr);
1327 }
1328 return result;
1329 }
1330
1336 bool update()
1337 {
1338 const TimerGuard prePostProcessGuard(prePostProcessTimer_);
1339
1340#ifndef NDEBUG
1341 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1342 // Make sure that the primary variables are defined. Note that because of padding
1343 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1344 // for definedness...
1345 for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1346 asImp_().solution(timeIdx)[i].checkDefined();
1347 }
1348 }
1349#endif // NDEBUG
1350
1351 // make sure all timers are prestine
1354 solveTimer_.halt();
1356
1358 asImp_().updateBegin();
1360
1361 bool converged = false;
1362
1363 try {
1364 converged = newtonMethod_.apply();
1365 }
1366 catch(...) {
1367 prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1368 linearizeTimer_ += newtonMethod_.linearizeTimer();
1369 solveTimer_ += newtonMethod_.solveTimer();
1370 updateTimer_ += newtonMethod_.updateTimer();
1371
1372 throw;
1373 }
1374
1375#ifndef NDEBUG
1376 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1377 // Make sure that the primary variables are defined. Note that because of padding
1378 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1379 // for definedness...
1380 for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1381 asImp_().solution(timeIdx)[i].checkDefined();
1382 }
1383 }
1384#endif // NDEBUG
1385
1386 prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1387 linearizeTimer_ += newtonMethod_.linearizeTimer();
1388 solveTimer_ += newtonMethod_.solveTimer();
1389 updateTimer_ += newtonMethod_.updateTimer();
1390
1392 if (converged) {
1393 asImp_().updateSuccessful();
1394 }
1395 else {
1396 asImp_().updateFailed();
1397 }
1399
1400#ifndef NDEBUG
1401 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1402 // Make sure that the primary variables are defined. Note that because of padding
1403 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1404 // for definedness...
1405 for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1406 asImp_().solution(timeIdx)[i].checkDefined();
1407 }
1408 }
1409#endif // NDEBUG
1410
1411 return converged;
1412 }
1413
1422 {}
1423
1430 {}
1431
1437 {}
1438
1443 {
1444 throw std::invalid_argument("Grid adaptation need to be implemented for "
1445 "specific settings of grid and function spaces");
1446 }
1447
1454 {
1455 // Reset the current solution to the one of the
1456 // previous time step so that we can start the next
1457 // update at a physically meaningful solution.
1458 solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
1460
1461#ifndef NDEBUG
1462 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1463 // Make sure that the primary variables are defined. Note that because of padding
1464 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1465 // for definedness...
1466 for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1467 asImp_().solution(timeIdx)[i].checkDefined();
1468 }
1469 }
1470#endif // NDEBUG
1471 }
1472
1481 {
1482 // at this point we can adapt the grid
1483 if (this->enableGridAdaptation_) {
1484 asImp_().adaptGrid();
1485 }
1486
1487 // make the current solution the previous one.
1488 solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
1489
1490 // shift the storage cache by one position in the history
1491 asImp_().shiftStorageCache(/*numSlots=*/1);
1492
1493 // shift the intensive quantities cache by one position in the
1494 // history
1495 asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
1496 }
1497
1505 template <class Restarter>
1506 void serialize(Restarter&)
1507 {
1508 throw std::runtime_error("Not implemented: The discretization chosen for this problem "
1509 "does not support restart files. (serialize() method unimplemented)");
1510 }
1511
1519 template <class Restarter>
1520 void deserialize(Restarter&)
1521 {
1522 throw std::runtime_error("Not implemented: The discretization chosen for this problem "
1523 "does not support restart files. (deserialize() method unimplemented)");
1524 }
1525
1534 template <class DofEntity>
1535 void serializeEntity(std::ostream& outstream,
1536 const DofEntity& dof)
1537 {
1538 const unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1539
1540 // write phase state
1541 if (!outstream.good()) {
1542 throw std::runtime_error("Could not serialize degree of freedom " +
1543 std::to_string(dofIdx));
1544 }
1545
1546 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1547 outstream << solution(/*timeIdx=*/0)[dofIdx][eqIdx] << " ";
1548 }
1549 }
1550
1559 template <class DofEntity>
1560 void deserializeEntity(std::istream& instream,
1561 const DofEntity& dof)
1562 {
1563 const unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1564
1565 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1566 if (!instream.good()) {
1567 throw std::runtime_error("Could not deserialize degree of freedom " +
1568 std::to_string(dofIdx));
1569 }
1570 instream >> solution(/*timeIdx=*/0)[dofIdx][eqIdx];
1571 }
1572 }
1573
1577 std::size_t numGridDof() const
1578 { throw std::logic_error("The discretization class must implement the numGridDof() method!"); }
1579
1583 std::size_t numAuxiliaryDof() const
1584 {
1585 return std::accumulate(auxEqModules_.begin(), auxEqModules_.end(),
1586 std::size_t{0},
1587 [](const auto acc, const auto& mod)
1588 { return acc + mod->numDofs(); });
1589 }
1590
1594 std::size_t numTotalDof() const
1595 { return asImp_().numGridDof() + numAuxiliaryDof(); }
1596
1601 const DofMapper& dofMapper() const
1602 { throw std::logic_error("The discretization class must implement the dofMapper() method!"); }
1603
1607 const VertexMapper& vertexMapper() const
1608 { return vertexMapper_; }
1609
1613 const ElementMapper& elementMapper() const
1614 { return elementMapper_; }
1615
1621 {
1622 linearizer_ = std::make_unique<Linearizer>();
1623 linearizer_->init(simulator_);
1624 }
1625
1629 static std::string discretizationName()
1630 { return ""; }
1631
1637 std::string primaryVarName(unsigned pvIdx) const
1638 {
1639 std::ostringstream oss;
1640 oss << "primary variable_" << pvIdx;
1641 return oss.str();
1642 }
1643
1649 std::string eqName(unsigned eqIdx) const
1650 {
1651 std::ostringstream oss;
1652 oss << "equation_" << eqIdx;
1653 return oss.str();
1654 }
1655
1662 void updatePVWeights(const ElementContext&) const
1663 {}
1664
1668 void addOutputModule(std::unique_ptr<BaseOutputModule<TypeTag>> newModule)
1669 { outputModules_.push_back(std::move(newModule)); }
1670
1679 template <class VtkMultiWriter>
1681 const SolutionVector& u,
1682 const GlobalEqVector& deltaU) const
1683 {
1684 using ScalarBuffer = std::vector<double>;
1685
1686 GlobalEqVector globalResid(u.size());
1687 asImp_().globalResidual(globalResid, u);
1688
1689 // create the required scalar fields
1690 const std::size_t numDof = asImp_().numGridDof();
1691
1692 // global defect of the two auxiliary equations
1693 std::array<ScalarBuffer*, numEq> def;
1694 std::array<ScalarBuffer*, numEq> delta;
1695 std::array<ScalarBuffer*, numEq> priVars;
1696 std::array<ScalarBuffer*, numEq> priVarWeight;
1697 ScalarBuffer* relError = writer.allocateManagedScalarBuffer(numDof);
1698 ScalarBuffer* normalizedRelError = writer.allocateManagedScalarBuffer(numDof);
1699 for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1700 priVars[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
1701 priVarWeight[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
1702 delta[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
1703 def[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
1704 }
1705
1706 Scalar minRelErr = 1e30;
1707 Scalar maxRelErr = -1e30;
1708 for (unsigned globalIdx = 0; globalIdx < numDof; ++ globalIdx) {
1709 for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1710 (*priVars[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
1711 (*priVarWeight[pvIdx])[globalIdx] = asImp_().primaryVarWeight(globalIdx, pvIdx);
1712 (*delta[pvIdx])[globalIdx] = - deltaU[globalIdx][pvIdx];
1713 (*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
1714 }
1715
1716 PrimaryVariables uOld(u[globalIdx]);
1717 PrimaryVariables uNew(uOld);
1718 uNew -= deltaU[globalIdx];
1719
1720 const Scalar err = asImp_().relativeDofError(globalIdx, uOld, uNew);
1721 (*relError)[globalIdx] = err;
1722 (*normalizedRelError)[globalIdx] = err;
1723 minRelErr = std::min(err, minRelErr);
1724 maxRelErr = std::max(err, maxRelErr);
1725 }
1726
1727 // do the normalization of the relative error
1728 const Scalar alpha = std::max(Scalar{1e-20},
1729 std::max(std::abs(maxRelErr),
1730 std::abs(minRelErr)));
1731 for (unsigned globalIdx = 0; globalIdx < numDof; ++globalIdx) {
1732 (*normalizedRelError)[globalIdx] /= alpha;
1733 }
1734
1735 DiscBaseOutputModule::attachScalarDofData_(writer, *relError, "relative error");
1736 DiscBaseOutputModule::attachScalarDofData_(writer, *normalizedRelError, "normalized relative error");
1737
1738 for (unsigned i = 0; i < numEq; ++i) {
1739 std::ostringstream oss;
1740 oss.str(""); oss << "priVar_" << asImp_().primaryVarName(i);
1741 DiscBaseOutputModule::attachScalarDofData_(writer,
1742 *priVars[i],
1743 oss.str());
1744
1745 oss.str(""); oss << "delta_" << asImp_().primaryVarName(i);
1746 DiscBaseOutputModule::attachScalarDofData_(writer,
1747 *delta[i],
1748 oss.str());
1749
1750 oss.str(""); oss << "weight_" << asImp_().primaryVarName(i);
1751 DiscBaseOutputModule::attachScalarDofData_(writer,
1752 *priVarWeight[i],
1753 oss.str());
1754
1755 oss.str(""); oss << "defect_" << asImp_().eqName(i);
1756 DiscBaseOutputModule::attachScalarDofData_(writer,
1757 *def[i],
1758 oss.str());
1759 }
1760
1761 asImp_().prepareOutputFields();
1762 asImp_().appendOutputFields(writer);
1763 }
1764
1770 {
1771 const bool needFullContextUpdate =
1772 std::any_of(outputModules_.begin(), outputModules_.end(),
1773 [](const auto& mod) { return mod->needExtensiveQuantities(); });
1774 std::for_each(outputModules_.begin(), outputModules_.end(),
1775 [](auto& mod) { mod->allocBuffers(); });
1776
1777 // iterate over grid
1778 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
1779#ifdef _OPENMP
1780#pragma omp parallel
1781#endif
1782 {
1783 ElementContext elemCtx(simulator_);
1784 ElementIterator elemIt = threadedElemIt.beginParallel();
1785 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
1786 const auto& elem = *elemIt;
1787 if (elem.partitionType() != Dune::InteriorEntity) {
1788 // ignore non-interior entities
1789 continue;
1790 }
1791
1792 if (needFullContextUpdate) {
1793 elemCtx.updateAll(elem);
1794 }
1795 else {
1796 elemCtx.updatePrimaryStencil(elem);
1797 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1798 }
1799
1800 std::for_each(outputModules_.begin(), outputModules_.end(),
1801 [&elemCtx](auto& mod) { mod->processElement(elemCtx); });
1802 }
1803 }
1804 }
1805
1811 {
1812 std::for_each(outputModules_.begin(), outputModules_.end(),
1813 [&writer](auto& mod) { mod->commitBuffers(writer); });
1814 }
1815
1819 const GridView& gridView() const
1820 { return gridView_; }
1821
1834 {
1835 auxMod->setDofOffset(numTotalDof());
1836 auxEqModules_.push_back(auxMod);
1837
1838 // resize the solutions
1839 if (enableGridAdaptation_ && !std::is_same_v<DiscreteFunction, BlockVectorWrapper>) {
1840 throw std::invalid_argument("Problems which require auxiliary modules cannot be used in"
1841 " conjunction with dune-fem");
1842 }
1843
1844 const std::size_t numDof = numTotalDof();
1845 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1846 solution(timeIdx).resize(numDof);
1847 }
1848
1849 auxMod->applyInitial();
1850 }
1851
1858 {
1859 auxEqModules_.clear();
1860 linearizer_->eraseMatrix();
1861 newtonMethod_.eraseMatrix();
1862 }
1863
1867 std::size_t numAuxiliaryModules() const
1868 { return auxEqModules_.size(); }
1869
1874 { return auxEqModules_[auxEqModIdx]; }
1875
1879 const BaseAuxiliaryModule<TypeTag>* auxiliaryModule(unsigned auxEqModIdx) const
1880 { return auxEqModules_[auxEqModIdx]; }
1881
1887
1889 { return prePostProcessTimer_; }
1890
1891 const Timer& linearizeTimer() const
1892 { return linearizeTimer_; }
1893
1894 const Timer& solveTimer() const
1895 { return solveTimer_; }
1896
1897 const Timer& updateTimer() const
1898 { return updateTimer_; }
1899
1900 template<class Serializer>
1901 void serializeOp(Serializer& serializer)
1902 {
1904 using Helper = typename BaseDiscretization::template SerializeHelper<Serializer>;
1905 Helper::serializeOp(serializer, solution_);
1906 }
1907
1908 bool operator==(const FvBaseDiscretization& rhs) const
1909 {
1910 return std::equal(this->solution_.begin(), this->solution_.end(),
1911 rhs.solution_.begin(), rhs.solution_.end(),
1912 [](const auto& x, const auto& y)
1913 {
1914 return *x == *y;
1915 });
1916 }
1917
1918protected:
1920 {
1921 // allocate the storage cache
1922 if (enableStorageCache()) {
1923 const std::size_t numDof = asImp_().numGridDof();
1924 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1925 storageCache_[timeIdx].resize(numDof);
1926 storageCacheUpToDate_[timeIdx].resize(numDof, /*value=*/0);
1927 }
1928 }
1929
1930 // allocate the intensive quantities cache
1932 const std::size_t numDof = asImp_().numGridDof();
1933 cachedIntensiveQuantityHistorySize_ = simulator_.problem().intensiveQuantityHistorySize();
1934 const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
1935
1936 // resize the vectors based on runtime history size
1937 intensiveQuantityCache_.resize(intensiveHistorySize);
1938 intensiveQuantityCacheUpToDate_.resize(intensiveHistorySize);
1939
1940 for(unsigned timeIdx = 0; timeIdx < intensiveHistorySize; ++timeIdx) {
1941 intensiveQuantityCache_[timeIdx].resize(numDof);
1942 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
1944 }
1945 }
1946 }
1947
1948 template <class Context>
1949 void supplementInitialSolution_(PrimaryVariables&,
1950 const Context&,
1951 unsigned,
1952 unsigned)
1953 {}
1954
1963 {
1964 // add the output modules available on all model
1965 this->outputModules_.push_back(std::make_unique<VtkPrimaryVarsModule<TypeTag>>(simulator_));
1966 }
1967
1971 LocalResidual& localResidual_()
1972 { return localLinearizer_.localResidual(); }
1973
1977 bool verbose_() const
1978 { return gridView_.comm().rank() == 0; }
1979
1980 Implementation& asImp_()
1981 { return *static_cast<Implementation*>(this); }
1982
1983 const Implementation& asImp_() const
1984 { return *static_cast<const Implementation*>(this); }
1985
1986 // the problem we want to solve. defines the constitutive
1987 // relations, matxerial laws, etc.
1988 Simulator& simulator_;
1989
1990 // the representation of the spatial domain of the problem
1991 GridView gridView_;
1992
1993 // the mappers for element and vertex entities to global indices
1994 ElementMapper elementMapper_;
1995 VertexMapper vertexMapper_;
1996
1997 // a vector with all auxiliary equations to be considered
1998 std::vector<BaseAuxiliaryModule<TypeTag>*> auxEqModules_;
1999
2000 NewtonMethod newtonMethod_;
2001
2006
2007 // calculates the local jacobian matrix for a given element
2008 std::vector<LocalLinearizer> localLinearizer_;
2009 // Linearizes the problem at the current time step using the
2010 // local jacobian
2011 std::unique_ptr<Linearizer> linearizer_;
2012
2013 // cur is the current iterative solution, prev the converged
2014 // solution of the previous time step
2015 mutable std::vector<IntensiveQuantitiesVector> intensiveQuantityCache_;
2016
2017 // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
2018 mutable std::vector<std::vector<unsigned char>> intensiveQuantityCacheUpToDate_;
2019
2020 mutable std::array<std::unique_ptr<DiscreteFunction>, historySize> solution_;
2021
2022 std::list<std::unique_ptr<BaseOutputModule<TypeTag>>> outputModules_;
2023
2025 std::vector<Scalar> dofTotalVolume_;
2026 std::vector<bool> isLocalDof_;
2027
2028 mutable std::array<GlobalEqVector, historySize> storageCache_;
2029
2030 // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
2031 mutable std::array<std::vector<unsigned char>, historySize> storageCacheUpToDate_;
2032
2037
2039};
2040
2046template<class TypeTag>
2048{
2052
2053 static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
2054
2055public:
2056 template<class Serializer>
2058 {
2059 template<class SolutionType>
2060 static void serializeOp(Serializer& serializer,
2061 SolutionType& solution)
2062 {
2063 for (auto& sol : solution) {
2064 serializer(*sol);
2065 }
2066 }
2067 };
2068
2069 explicit FvBaseDiscretizationNoAdapt(Simulator& simulator)
2070 : ParentType(simulator)
2071 {
2072 if (this->enableGridAdaptation_) {
2073 throw std::invalid_argument("Grid adaptation need to use"
2074 " BaseDiscretization = FvBaseDiscretizationFemAdapt"
2075 " which currently requires the presence of the"
2076 " dune-fem module");
2077 }
2078 const std::size_t numDof = this->asImp_().numGridDof();
2079 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
2080 this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", numDof);
2081 }
2082 }
2083};
2084
2085} // namespace Opm
2086
2087#endif // EWOMS_FV_BASE_DISCRETIZATION_HH
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1....
Base class for specifying auxiliary equations.
Definition: baseauxiliarymodule.hh:56
virtual void applyInitial()=0
Set the initial condition of the auxiliary module in the solution vector.
void setDofOffset(int value)
Set the offset in the global system of equations for the first degree of freedom of this auxiliary mo...
Definition: baseauxiliarymodule.hh:78
The base class for writer modules.
Definition: baseoutputmodule.hh:68
The base class for all output writers.
Definition: baseoutputwriter.hh:46
Represents all quantities which available on boundary segments.
Definition: fvbaseboundarycontext.hh:46
Represents all quantities which available for calculating constraints.
Definition: fvbaseconstraintscontext.hh:44
Class to specify constraints for a finite volume spatial discretization.
Definition: fvbaseconstraints.hh:48
Definition: fvbasediscretization.hh:350
const SolutionVector & blockVector() const
Definition: fvbasediscretization.hh:373
SolutionVector blockVector_
Definition: fvbasediscretization.hh:352
void serializeOp(Serializer &serializer)
Definition: fvbasediscretization.hh:383
static BlockVectorWrapper serializationTestObject()
Definition: fvbasediscretization.hh:360
SolutionVector & blockVector()
Definition: fvbasediscretization.hh:370
bool operator==(const BlockVectorWrapper &wrapper) const
Definition: fvbasediscretization.hh:376
BlockVectorWrapper(const std::string &, const std::size_t size)
Definition: fvbasediscretization.hh:354
The base class for the finite volume discretization schemes without adaptation.
Definition: fvbasediscretization.hh:2048
FvBaseDiscretizationNoAdapt(Simulator &simulator)
Definition: fvbasediscretization.hh:2069
The base class for the finite volume discretization schemes.
Definition: fvbasediscretization.hh:298
Timer linearizeTimer_
Definition: fvbasediscretization.hh:2003
std::vector< IntensiveQuantitiesVector > intensiveQuantityCache_
Definition: fvbasediscretization.hh:2015
LocalLinearizer & localLinearizer(unsigned openMpThreadId)
Definition: fvbasediscretization.hh:1269
void prepareOutputFields() const
Prepare the quantities relevant for the current solution to be appended to the output writers.
Definition: fvbasediscretization.hh:1769
void shiftIntensiveQuantityCache(unsigned numSlots=1)
Move the intensive quantities for a given time index to the back.
Definition: fvbasediscretization.hh:785
void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
Definition: fvbasediscretization.hh:727
void adaptGrid()
Called by the update() method when the grid should be refined.
Definition: fvbasediscretization.hh:1442
std::vector< BaseAuxiliaryModule< TypeTag > * > auxEqModules_
Definition: fvbasediscretization.hh:1998
const Implementation & asImp_() const
Definition: fvbasediscretization.hh:1983
void addAuxiliaryModule(BaseAuxiliaryModule< TypeTag > *auxMod)
Add a module for an auxiliary equation.
Definition: fvbasediscretization.hh:1833
void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx, unsigned timeIdx, bool newValue) const
Invalidate the cache for a given intensive quantities object.
Definition: fvbasediscretization.hh:692
void finishInit()
Apply the initial conditions to the model.
Definition: fvbasediscretization.hh:458
void prefetch(const Element &) const
Allows to improve the performance by prefetching all data which is associated with a given element.
Definition: fvbasediscretization.hh:591
bool enableStorageCache_
Definition: fvbasediscretization.hh:2035
void resizeAndResetIntensiveQuantitiesCache_()
Definition: fvbasediscretization.hh:1919
std::vector< Scalar > dofTotalVolume_
Definition: fvbasediscretization.hh:2025
void updateSuccessful()
Called by the update() method if it was successful.
Definition: fvbasediscretization.hh:1436
static std::string discretizationName()
Returns a string of discretization's human-readable name.
Definition: fvbasediscretization.hh:1629
unsigned cachedIntensiveQuantityHistorySize() const
Get the cached intensive quantity history size.
Definition: fvbasediscretization.hh:706
BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx)
Returns a given module for auxiliary equations.
Definition: fvbasediscretization.hh:1873
bool isLocalDof(unsigned globalIdx) const
Returns if the overlap of the volume ofa degree of freedom is non-zero.
Definition: fvbasediscretization.hh:1209
std::size_t numAuxiliaryModules() const
Returns the number of modules for auxiliary equations.
Definition: fvbasediscretization.hh:1867
bool operator==(const FvBaseDiscretization &rhs) const
Definition: fvbasediscretization.hh:1908
void serializeOp(Serializer &serializer)
Definition: fvbasediscretization.hh:1901
std::vector< bool > isLocalDof_
Definition: fvbasediscretization.hh:2026
LocalResidual & localResidual_()
Reference to the local residal object.
Definition: fvbasediscretization.hh:1971
void registerOutputModules_()
Register all output modules which make sense for the model.
Definition: fvbasediscretization.hh:1962
bool enableGridAdaptation() const
Returns whether the grid ought to be adapted to the solution during the simulation.
Definition: fvbasediscretization.hh:523
const NewtonMethod & newtonMethod() const
Returns the newton method object.
Definition: fvbasediscretization.hh:605
SolutionVector & mutableSolution(unsigned timeIdx) const
Definition: fvbasediscretization.hh:1237
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition: fvbasediscretization.hh:1480
NewtonMethod & newtonMethod()
Returns the newton method object.
Definition: fvbasediscretization.hh:599
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition: fvbasediscretization.hh:1607
const EqVector & cachedStorage(unsigned globalIdx, unsigned timeIdx) const
Retrieve an entry of the cache for the storage term.
Definition: fvbasediscretization.hh:843
void updatePVWeights(const ElementContext &) const
Update the weights of all primary variables within an element given the complete set of intensive qua...
Definition: fvbasediscretization.hh:1662
const Timer & solveTimer() const
Definition: fvbasediscretization.hh:1894
void updateFailed()
Called by the update() method if it was unsuccessful. This is primary a hook which the actual model c...
Definition: fvbasediscretization.hh:1453
void updateBegin()
Called by the update() method before it tries to apply the newton method. This is primary a hook whic...
Definition: fvbasediscretization.hh:1429
FvBaseDiscretization(const FvBaseDiscretization &)=delete
Scalar globalResidual(GlobalEqVector &dest) const
Compute the global residual for the current solution vector.
Definition: fvbasediscretization.hh:964
LocalResidual & localResidual(unsigned openMpThreadId)
Definition: fvbasediscretization.hh:1281
void serializeEntity(std::ostream &outstream, const DofEntity &dof)
Write the current solution for a degree of freedom to a restart file.
Definition: fvbasediscretization.hh:1535
void setEnableStorageCache(bool enableStorageCache)
Set the value of enable storage cache.
Definition: fvbasediscretization.hh:827
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: fvbasediscretization.hh:1649
Scalar eqWeight(unsigned, unsigned) const
Returns the relative weight of an equation.
Definition: fvbasediscretization.hh:1303
std::array< std::vector< unsigned char >, historySize > storageCacheUpToDate_
Definition: fvbasediscretization.hh:2031
const Timer & prePostProcessTimer() const
Definition: fvbasediscretization.hh:1888
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: fvbasediscretization.hh:1291
void deserialize(Restarter &)
Deserializes the state of the model.
Definition: fvbasediscretization.hh:1520
void checkConservativeness(Scalar tolerance=-1, bool verbose=false) const
Ensure that the difference between the storage terms of the last and of the current time step is cons...
Definition: fvbasediscretization.hh:1079
Scalar gridTotalVolume() const
Returns the volume of the whole grid which represents the spatial domain.
Definition: fvbasediscretization.hh:1216
Timer updateTimer_
Definition: fvbasediscretization.hh:2005
FvBaseDiscretization(Simulator &simulator)
Definition: fvbasediscretization.hh:394
const IntensiveQuantities * thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
Return the thermodynamic hint for a entity on the grid at given time.
Definition: fvbasediscretization.hh:623
bool storageCacheIsUpToDate(unsigned globalIdx, unsigned timeIdx) const
Returns true if the storage cache entry for a given DOF and time index is up to date.
Definition: fvbasediscretization.hh:883
const Timer & updateTimer() const
Definition: fvbasediscretization.hh:1897
const Linearizer & linearizer() const
Returns the operator linearizer for the global jacobian of the problem.
Definition: fvbasediscretization.hh:1245
void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector &value) const
Set an entry of the cache for the storage term.
Definition: fvbasediscretization.hh:867
void addConvergenceVtkFields(VtkMultiWriter &writer, const SolutionVector &u, const GlobalEqVector &deltaU) const
Add the vector fields for analysing the convergence of the newton method to the a VTK writer.
Definition: fvbasediscretization.hh:1680
Scalar relativeDofError(unsigned vertexIdx, const PrimaryVariables &pv1, const PrimaryVariables &pv2) const
Returns the relative error between two vectors of primary variables.
Definition: fvbasediscretization.hh:1315
bool enableGridAdaptation_
Definition: fvbasediscretization.hh:2033
Scalar globalResidual(GlobalEqVector &dest, const SolutionVector &u) const
Compute the global residual for an arbitrary solution vector.
Definition: fvbasediscretization.hh:949
SolutionVector & solution(unsigned timeIdx)
Definition: fvbasediscretization.hh:1230
std::array< GlobalEqVector, historySize > storageCache_
Definition: fvbasediscretization.hh:2028
Timer prePostProcessTimer_
Definition: fvbasediscretization.hh:2002
void deserializeEntity(std::istream &instream, const DofEntity &dof)
Reads the current solution variables for a degree of freedom from a restart file.
Definition: fvbasediscretization.hh:1560
Timer solveTimer_
Definition: fvbasediscretization.hh:2004
void supplementInitialSolution_(PrimaryVariables &, const Context &, unsigned, unsigned)
Definition: fvbasediscretization.hh:1949
const LocalLinearizer & localLinearizer(unsigned openMpThreadId) const
Returns the local jacobian which calculates the local stiffness matrix for an arbitrary element.
Definition: fvbasediscretization.hh:1263
std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_
Definition: fvbasediscretization.hh:2020
unsigned cachedIntensiveQuantityHistorySize_
Definition: fvbasediscretization.hh:2038
const Timer & linearizeTimer() const
Definition: fvbasediscretization.hh:1891
void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType &gridView) const
Definition: fvbasediscretization.hh:748
bool update()
Try to progress the model to the next timestep.
Definition: fvbasediscretization.hh:1336
std::size_t numAuxiliaryDof() const
Returns the number of degrees of freedom (DOFs) of the auxiliary equations.
Definition: fvbasediscretization.hh:1583
void clearAuxiliaryModules()
Causes the list of auxiliary equations to be cleared.
Definition: fvbasediscretization.hh:1857
bool enableIntensiveQuantityCache_
Definition: fvbasediscretization.hh:2034
std::unique_ptr< Linearizer > linearizer_
Definition: fvbasediscretization.hh:2011
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition: fvbasediscretization.hh:1613
void syncOverlap()
Syncronize the values of the primary variables on the degrees of freedom that overlap with the neighb...
Definition: fvbasediscretization.hh:1421
void appendOutputFields(BaseOutputWriter &writer) const
Append the quantities relevant for the current solution to an output writer.
Definition: fvbasediscretization.hh:1810
std::list< std::unique_ptr< BaseOutputModule< TypeTag > > > outputModules_
Definition: fvbasediscretization.hh:2022
ElementMapper elementMapper_
Definition: fvbasediscretization.hh:1994
NewtonMethod newtonMethod_
Definition: fvbasediscretization.hh:2000
std::size_t numTotalDof() const
Returns the total number of degrees of freedom (i.e., grid plux auxiliary DOFs)
Definition: fvbasediscretization.hh:1594
void applyInitialSolution()
Applies the initial solution for all degrees of freedom to which the model applies.
Definition: fvbasediscretization.hh:530
Implementation & asImp_()
Definition: fvbasediscretization.hh:1980
void updateCachedIntensiveQuantities(const IntensiveQuantities &intQuants, unsigned globalIdx, unsigned timeIdx) const
Update the intensive quantity cache for a entity on the grid at given time.
Definition: fvbasediscretization.hh:672
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: fvbasediscretization.hh:1637
void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
Invalidate the whole intensive quantity cache for time index.
Definition: fvbasediscretization.hh:714
Linearizer & linearizer()
Returns the object which linearizes the global system of equations at the current solution.
Definition: fvbasediscretization.hh:1252
const BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx) const
Returns a given module for auxiliary equations.
Definition: fvbasediscretization.hh:1879
std::vector< std::vector< unsigned char > > intensiveQuantityCacheUpToDate_
Definition: fvbasediscretization.hh:2018
void globalStorage(EqVector &storage, unsigned timeIdx=0) const
Compute the integral over the domain of the storage terms of all conservation quantities.
Definition: fvbasediscretization.hh:1024
GridView gridView_
Definition: fvbasediscretization.hh:1991
Scalar dofTotalVolume(unsigned globalIdx) const
Returns the volume of a given control volume.
Definition: fvbasediscretization.hh:1201
static void registerParameters()
Register all run-time parameters for the model.
Definition: fvbasediscretization.hh:427
const SolutionVector & solution(unsigned timeIdx) const
Reference to the solution at a given history index as a block vector.
Definition: fvbasediscretization.hh:1224
bool verbose_() const
Returns whether messages should be printed.
Definition: fvbasediscretization.hh:1977
void shiftStorageCache(unsigned numSlots=1) const
Shift storage cache by a given number of time step slots.
Definition: fvbasediscretization.hh:929
Simulator & simulator_
Definition: fvbasediscretization.hh:1988
const LocalResidual & localResidual(unsigned openMpThreadId) const
Returns the object to calculate the local residual function.
Definition: fvbasediscretization.hh:1275
std::size_t numGridDof() const
Returns the number of degrees of freedom (DOFs) for the computational grid.
Definition: fvbasediscretization.hh:1577
Scalar gridTotalVolume_
Definition: fvbasediscretization.hh:2024
const DofMapper & dofMapper() const
Mapper to convert the Dune entities of the discretization's degrees of freedoms are to indices.
Definition: fvbasediscretization.hh:1601
bool storeIntensiveQuantities() const
Returns true if the cache for intensive quantities is enabled.
Definition: fvbasediscretization.hh:1885
std::vector< LocalLinearizer > localLinearizer_
Definition: fvbasediscretization.hh:2008
void invalidateStorageCache(unsigned timeIdx) const
Invalidate the whole storage cache for a given time index.
Definition: fvbasediscretization.hh:909
void serialize(Restarter &)
Serializes the current state of the model.
Definition: fvbasediscretization.hh:1506
bool enableThermodynamicHints_
Definition: fvbasediscretization.hh:2036
void invalidateStorageCacheEntry(unsigned globalIdx, unsigned timeIdx) const
Invalidate the storage cache for a given DOF and time index.
Definition: fvbasediscretization.hh:897
VertexMapper vertexMapper_
Definition: fvbasediscretization.hh:1995
void resetLinearizer()
Resets the Jacobian matrix linearizer, so that the boundary types can be altered.
Definition: fvbasediscretization.hh:1620
bool enableStorageCache() const
Returns true iff the storage term is cached.
Definition: fvbasediscretization.hh:818
const GridView & gridView() const
Reference to the grid view of the spatial domain.
Definition: fvbasediscretization.hh:1819
void addOutputModule(std::unique_ptr< BaseOutputModule< TypeTag > > newModule)
Add an module for writing visualization output after a timestep.
Definition: fvbasediscretization.hh:1668
const IntensiveQuantities * cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
Return the cached intensive quantities for a entity on the grid at given time.
Definition: fvbasediscretization.hh:644
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
Definition: fvbaseelementcontext.hh:55
Provide the properties at a face which make sense independently of the conserved quantities.
Definition: fvbaseextensivequantities.hh:48
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition: fvbasegradientcalculator.hh:52
Base class for the model specific class which provides access to all intensive (i....
Definition: fvbaseintensivequantities.hh:45
The common code for the linearizers of non-linear systems of equations.
Definition: fvbaselinearizer.hh:78
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition: fvbaselocalresidual.hh:63
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:53
This is a grid manager which does not create any border list.
Definition: nullborderlistmanager.hh:44
static void registerParameters()
Register all run-time parameters for the Newton method.
Definition: newtonmethod.hh:136
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:84
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.
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
Definition: threadedentityiterator.hh:42
bool isFinished(const EntityIterator &it) const
Definition: threadedentityiterator.hh:67
EntityIterator increment()
Definition: threadedentityiterator.hh:80
EntityIterator beginParallel()
Definition: threadedentityiterator.hh:54
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition: timerguard.hh:42
Provides an encapsulation to measure the system time.
Definition: timer.hpp:46
void start()
Start counting the time resources used by the simulation.
void halt()
Stop the measurement reset all timing values.
double stop()
Stop counting the time resources.
Simplifies writing multi-file VTK datasets.
Definition: vtkmultiwriter.hh:65
ScalarBuffer * allocateManagedScalarBuffer(std::size_t numEntities)
Allocate a managed buffer for a scalar field.
Definition: vtkmultiwriter.hh:192
VTK output module for the fluid composition.
Definition: vtkprimaryvarsmodule.hpp:48
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkprimaryvarsmodule.hpp:74
Definition: alignedallocator.hh:97
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides data handles for parallel communication which operate on DOFs.
Declares the parameters for the black oil model.
Definition: fvbaseprimaryvariables.hh:141
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:187
Definition: blackoilmodel.hh:79
Definition: blackoilboundaryratevector.hh:39
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:233
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Definition: fvbasediscretization.hh:2058
static void serializeOp(Serializer &serializer, SolutionType &solution)
Definition: fvbasediscretization.hh:2060
Definition: fvbasediscretization.hh:272
GetPropType< TypeTag, Properties::GridView > GridView
Definition: fvbasediscretization.hh:113
GetPropType< TypeTag, Properties::DofMapper > DofMapper
Definition: fvbasediscretization.hh:112
The class which marks the border indices associated with the degrees of freedom on a process boundary...
Definition: basicproperties.hh:125
The secondary variables of a boundary segment.
Definition: fvbaseproperties.hh:143
GetPropType< TypeTag, Properties::RateVector > type
Definition: fvbasediscretization.hh:160
Type of object for specifying boundary conditions.
Definition: fvbaseproperties.hh:119
The secondary variables of a constraint degree of freedom.
Definition: fvbaseproperties.hh:146
The class which represents a constraint degree of freedom.
Definition: fvbaseproperties.hh:122
The part of the extensive quantities which is specific to the spatial discretization.
Definition: fvbaseproperties.hh:160
The discretization specific part of the intensive quantities.
Definition: fvbaseproperties.hh:136
The discretization specific part of the local residual.
Definition: fvbaseproperties.hh:91
typename BaseDiscretization::BlockVectorWrapper type
Definition: fvbasediscretization.hh:283
Definition: fvbaseproperties.hh:77
The secondary variables of all degrees of freedom in an element's stencil.
Definition: fvbaseproperties.hh:140
Dune::BlockVector< GetPropType< TypeTag, Properties::EqVector > > type
Definition: fvbasediscretization.hh:174
A vector of holding a quantity for each equation for each DOF of an element.
Definition: fvbaseproperties.hh:112
Dune::MultipleCodimMultipleGeomTypeMapper< GetPropType< TypeTag, Properties::GridView > > type
Definition: fvbasediscretization.hh:106
The mapper to find the global index of an element.
Definition: fvbaseproperties.hh:213
Specify whether the some degrees of fredom can be constraint.
Definition: fvbaseproperties.hh:199
Specify if experimental features should be enabled or not.
Definition: fvbaseproperties.hh:241
Dune::FieldVector< GetPropType< TypeTag, Properties::Scalar >, getPropValue< TypeTag, Properties::NumEq >()> type
Definition: fvbasediscretization.hh:143
A vector of holding a quantity for each equation (usually at a given spatial location)
Definition: fvbaseproperties.hh:109
Specify whether the storage terms use extensive quantities or not.
Definition: fvbaseproperties.hh:233
Dune::BlockVector< GetPropType< TypeTag, Properties::EqVector > > type
Definition: fvbasediscretization.hh:181
Vector containing a quantity of for equation for each DOF of the whole grid.
Definition: linalgproperties.hh:54
Calculates gradients of arbitrary quantities at flux integration points.
Definition: fvbaseproperties.hh:152
The secondary variables within a sub-control volume.
Definition: fvbaseproperties.hh:133
The class which linearizes the non-linear system of equations.
Definition: newtonmethodproperties.hh:36
A vector of primary variables within a sub-control volume.
Definition: fvbaseproperties.hh:130
GetPropType< TypeTag, Properties::EqVector > type
Definition: fvbasediscretization.hh:153
Vector containing volumetric or areal rates of quantities.
Definition: fvbaseproperties.hh:116
Manages the simulation time.
Definition: basicproperties.hh:116
Dune::BlockVector< GetPropType< TypeTag, Properties::PrimaryVariables > > type
Definition: fvbasediscretization.hh:195
Vector containing all primary variables of the grid.
Definition: fvbaseproperties.hh:126
The OpenMP threads manager.
Definition: fvbaseproperties.hh:174
The history size required by the time discretization.
Definition: fvbaseproperties.hh:225
a tag to mark properties as undefined
Definition: propertysystem.hh:38
Definition: fvbaseproperties.hh:181
Specify whether to use volumetric residuals or not.
Definition: fvbaseproperties.hh:237
Dune::MultipleCodimMultipleGeomTypeMapper< GetPropType< TypeTag, Properties::GridView > > type
Definition: fvbasediscretization.hh:101
The mapper to find the global index of a vertex.
Definition: fvbaseproperties.hh:207
Specify the format the VTK output is written to disk.
Definition: fvbaseproperties.hh:195