My Project
Loading...
Searching...
No Matches
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/version.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/fmatrix.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 <cstddef>
71#include <list>
72#include <stdexcept>
73#include <sstream>
74#include <string>
75#include <type_traits>
76#include <vector>
77
78namespace Opm {
79template<class TypeTag>
80class FvBaseDiscretizationNoAdapt;
81template<class TypeTag>
82class FvBaseDiscretization;
83
84} // namespace Opm
85
86namespace Opm::Properties {
87
89template<class TypeTag>
90struct Simulator<TypeTag, TTag::FvBaseDiscretization>
92
94template<class TypeTag>
95struct VertexMapper<TypeTag, TTag::FvBaseDiscretization>
96{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
97
99template<class TypeTag>
101{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
102
104template<class TypeTag>
113
114template<class TypeTag>
117
118template<class TypeTag>
121
122template<class TypeTag>
125
127template<class TypeTag>
130
134template<class TypeTag>
135struct EqVector<TypeTag, TTag::FvBaseDiscretization>
136{
137 using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
139};
140
146template<class TypeTag>
149
153template<class TypeTag>
156
160template<class TypeTag>
161struct Constraints<TypeTag, TTag::FvBaseDiscretization>
163
167template<class TypeTag>
169{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
170
174template<class TypeTag>
176{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
177
181template<class TypeTag>
184
188template<class TypeTag>
190{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; };
191
197template<class TypeTag>
200
204template<class TypeTag>
207
208template<class TypeTag>
211
212template<class TypeTag>
215
219template<class TypeTag>
221{ using type = ::Opm::ThreadManager; };
222
223template<class TypeTag>
225{ static constexpr bool value = true; };
226
230template<class TypeTag>
231struct Linearizer<TypeTag, TTag::FvBaseDiscretization>
233
235template<class TypeTag>
237{ static constexpr int value = Dune::VTK::ascii; };
238
239// disable constraints by default
240template<class TypeTag>
242{ static constexpr bool value = false; };
243
245template<class TypeTag>
247{ static constexpr int value = 2; };
248
251template<class TypeTag>
253{ static constexpr bool value = false; };
254
255// use volumetric residuals is default
256template<class TypeTag>
258{ static constexpr bool value = true; };
259
262template<class TypeTag>
264{ static constexpr bool value = true; };
265
266template <class TypeTag, class MyTypeTag>
268
269#if !HAVE_DUNE_FEM
270template<class TypeTag>
273
274template<class TypeTag>
280#endif
281
282} // namespace Opm::Properties
283
284namespace Opm {
285
291template<class TypeTag>
293{
294 using Implementation = GetPropType<TypeTag, Properties::Model>;
321
324
325 enum {
328 };
329
330 using IntensiveQuantitiesVector = std::vector<IntensiveQuantities, aligned_allocator<IntensiveQuantities, alignof(IntensiveQuantities)> >;
331
332 using Element = typename GridView::template Codim<0>::Entity;
333 using ElementIterator = typename GridView::template Codim<0>::Iterator;
334
335 using Toolbox = MathToolbox<Evaluation>;
336 using VectorBlock = Dune::FieldVector<Evaluation, numEq>;
337 using EvalEqVector = Dune::FieldVector<Evaluation, numEq>;
338
339 using LocalEvalBlockVector = typename LocalResidual::LocalEvalBlockVector;
340
341public:
343 {
344 protected:
345 SolutionVector blockVector_;
346 public:
347 BlockVectorWrapper(const std::string&, const size_t size)
348 : blockVector_(size)
349 {}
350
351 BlockVectorWrapper() = default;
352
353 static BlockVectorWrapper serializationTestObject()
354 {
355 BlockVectorWrapper result("dummy", 3);
356 result.blockVector_[0] = 1.0;
357 result.blockVector_[1] = 2.0;
358 result.blockVector_[2] = 3.0;
359
360 return result;
361 }
362
363 SolutionVector& blockVector()
364 { return blockVector_; }
365 const SolutionVector& blockVector() const
366 { return blockVector_; }
367
368 bool operator==(const BlockVectorWrapper& wrapper) const
369 {
370 return std::equal(this->blockVector_.begin(), this->blockVector_.end(),
371 wrapper.blockVector_.begin(), wrapper.blockVector_.end());
372 }
373
374 template<class Serializer>
375 void serializeOp(Serializer& serializer)
376 {
377 serializer(blockVector_);
378 }
379 };
380
381private:
384
385 // copying a discretization object is not a good idea
387
388public:
389 // this constructor required to be explicitly specified because
390 // we've defined a constructor above which deletes all implicitly
391 // generated constructors in C++.
392 explicit FvBaseDiscretization(Simulator& simulator)
393 : simulator_(simulator)
394 , gridView_(simulator.gridView())
395 , elementMapper_(gridView_, Dune::mcmgElementLayout())
396 , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
397 , newtonMethod_(simulator)
398 , localLinearizer_(ThreadManager::maxThreads())
399 , linearizer_(new Linearizer())
400 , enableGridAdaptation_(Parameters::Get<Parameters::EnableGridAdaptation>() )
401 , enableIntensiveQuantityCache_(Parameters::Get<Parameters::EnableIntensiveQuantityCache>())
402 , enableStorageCache_(Parameters::Get<Parameters::EnableStorageCache>())
403 , enableThermodynamicHints_(Parameters::Get<Parameters::EnableThermodynamicHints>())
404 {
405 bool isEcfv = std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value;
406 if (enableGridAdaptation_ && !isEcfv)
407 throw std::invalid_argument("Grid adaptation currently only works for the "
408 "element-centered finite volume discretization (is: "
409 +Dune::className<Discretization>()+")");
410
411 PrimaryVariables::init();
412 size_t numDof = asImp_().numGridDof();
413 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
415 intensiveQuantityCache_[timeIdx].resize(numDof);
416 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
417 }
418
419 if (enableStorageCache_)
420 storageCache_[timeIdx].resize(numDof);
421 }
422
423 resizeAndResetIntensiveQuantitiesCache_();
424 asImp_().registerOutputModules_();
425 }
426
427 ~FvBaseDiscretization()
428 {
429 // delete all output modules
430 auto modIt = outputModules_.begin();
431 const auto& modEndIt = outputModules_.end();
432 for (; modIt != modEndIt; ++modIt)
433 delete *modIt;
434
435 delete linearizer_;
436 }
437
441 static void registerParameters()
442 {
443 Linearizer::registerParameters();
444 LocalLinearizer::registerParameters();
445 LocalResidual::registerParameters();
446 GradientCalculator::registerParameters();
447 IntensiveQuantities::registerParameters();
448 ExtensiveQuantities::registerParameters();
450 Linearizer::registerParameters();
451 PrimaryVariables::registerParameters();
452 // register runtime parameters of the output modules
454
455 Parameters::Register<Parameters::EnableGridAdaptation>
456 ("Enable adaptive grid refinement/coarsening");
457 Parameters::Register<Parameters::EnableVtkOutput>
458 ("Global switch for turning on writing VTK files");
459 Parameters::Register<Parameters::EnableThermodynamicHints>
460 ("Enable thermodynamic hints");
461 Parameters::Register<Parameters::EnableIntensiveQuantityCache>
462 ("Turn on caching of intensive quantities");
463 Parameters::Register<Parameters::EnableStorageCache>
464 ("Store previous storage terms and avoid re-calculating them.");
465 Parameters::Register<Parameters::OutputDir>
466 ("The directory to which result files are written");
467 }
468
473 {
474 // initialize the volume of the finite volumes to zero
475 size_t numDof = asImp_().numGridDof();
476 dofTotalVolume_.resize(numDof);
477 std::fill(dofTotalVolume_.begin(), dofTotalVolume_.end(), 0.0);
478
479 ElementContext elemCtx(simulator_);
480 gridTotalVolume_ = 0.0;
481
482 // iterate through the grid and evaluate the initial condition
483 for (const auto& elem : elements(gridView_)) {
484 const bool isInteriorElement = elem.partitionType() == Dune::InteriorEntity;
485 // ignore everything which is not in the interior if the
486 // current process' piece of the grid
488 continue;
489
490 // deal with the current element
491 elemCtx.updateStencil(elem);
492 const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
493
494 // loop over all element vertices, i.e. sub control volumes
495 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++) {
496 // map the local degree of freedom index to the global one
497 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
498
499 Scalar dofVolume = stencil.subControlVolume(dofIdx).volume();
500 dofTotalVolume_[globalIdx] += dofVolume;
501 gridTotalVolume_ += dofVolume;
502 }
503 }
504
505 // determine which DOFs should be considered to lie fully in the interior of the
506 // local process grid partition: those which do not have a non-zero volume
507 // before taking the peer processes into account...
508 isLocalDof_.resize(numDof);
509 for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx)
510 isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0);
511
512 // add the volumes of the DOFs on the process boundaries
513 const auto sumHandle =
514 GridCommHandleFactory::template sumHandle<Scalar>(dofTotalVolume_,
515 asImp_().dofMapper());
516 gridView_.communicate(*sumHandle,
517 Dune::InteriorBorder_All_Interface,
518 Dune::ForwardCommunication);
519
520 // sum up the volumes of the grid partitions
521 gridTotalVolume_ = gridView_.comm().sum(gridTotalVolume_);
522
523 linearizer_->init(simulator_);
524 for (unsigned threadId = 0; threadId < ThreadManager::maxThreads(); ++threadId)
525 localLinearizer_[threadId].init(simulator_);
526
527 resizeAndResetIntensiveQuantitiesCache_();
529 // invalidate all cached intensive quantities
530 for (unsigned timeIdx = 0; timeIdx < historySize; ++ timeIdx)
532 }
533
534 newtonMethod_.finishInit();
535 }
536
541 { return enableGridAdaptation_; }
542
548 {
549 // first set the whole domain to zero
550 SolutionVector& uCur = asImp_().solution(/*timeIdx=*/0);
551 uCur = Scalar(0.0);
552
553 ElementContext elemCtx(simulator_);
554
555 // iterate through the grid and evaluate the initial condition
556 for (const auto& elem : elements(gridView_)) {
557 // ignore everything which is not in the interior if the
558 // current process' piece of the grid
559 if (elem.partitionType() != Dune::InteriorEntity)
560 continue;
561
562 // deal with the current element
563 elemCtx.updateStencil(elem);
564
565 // loop over all element vertices, i.e. sub control volumes
566 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++)
567 {
568 // map the local degree of freedom index to the global one
569 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
570
571 // let the problem do the dirty work of nailing down
572 // the initial solution.
573 simulator_.problem().initial(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
574 asImp_().supplementInitialSolution_(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
575 uCur[globalIdx].checkDefined();
576 }
577 }
578
579 // synchronize the ghost DOFs (if necessary)
580 asImp_().syncOverlap();
581
582 // also set the solutions of the "previous" time steps to the initial solution.
583 for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx)
584 solution(timeIdx) = solution(/*timeIdx=*/0);
585
586 simulator_.problem().initialSolutionApplied();
587
588#ifndef NDEBUG
589 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
590 const auto& sol = solution(timeIdx);
591 for (unsigned dofIdx = 0; dofIdx < sol.size(); ++dofIdx)
592 sol[dofIdx].checkDefined();
593 }
594#endif // NDEBUG
595 }
596
601 void prefetch(const Element&) const
602 {
603 // do nothing by default
604 }
605
609 NewtonMethod& newtonMethod()
610 { return newtonMethod_; }
611
615 const NewtonMethod& newtonMethod() const
616 { return newtonMethod_; }
617
633 const IntensiveQuantities* thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
634 {
635 if (!enableThermodynamicHints_)
636 return 0;
637
638 // the intensive quantities cache doubles as thermodynamic hint
639 return cachedIntensiveQuantities(globalIdx, timeIdx);
640 }
641
653 const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
654 {
655 if (!enableIntensiveQuantityCache_ || !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
656 return nullptr;
657 }
658
659 // With the storage cache enabled, usually only the
660 // intensive quantities for the most recent time step are
661 // cached. However, this may be false for some Problem
662 // variants, so we should check if the cache exists for
663 // the timeIdx in question.
664 if (timeIdx > 0 && enableStorageCache_ && intensiveQuantityCache_[timeIdx].empty()) {
665 return nullptr;
666 }
667
668 return &intensiveQuantityCache_[timeIdx][globalIdx];
669 }
670
679 void updateCachedIntensiveQuantities(const IntensiveQuantities& intQuants,
680 unsigned globalIdx,
681 unsigned timeIdx) const
682 {
684 return;
685
686 intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
687 intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = 1;
688 }
689
698 unsigned timeIdx,
699 bool newValue) const
700 {
702 return;
703
704 intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue ? 1 : 0;
705 }
706
713 {
715 std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
716 intensiveQuantityCacheUpToDate_[timeIdx].end(),
717 /*value=*/0);
718 }
719 }
720
721 void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
722 {
724
725 // loop over all elements...
726 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
727#ifdef _OPENMP
728#pragma omp parallel
729#endif
730 {
731 ElementContext elemCtx(simulator_);
732 ElementIterator elemIt = threadedElemIt.beginParallel();
733 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
734 const Element& elem = *elemIt;
735 elemCtx.updatePrimaryStencil(elem);
736 elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
737 }
738 }
739 }
740
741 template <class GridViewType>
742 void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
743 {
744 // loop over all elements...
745 ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(gridView);
746#ifdef _OPENMP
747#pragma omp parallel
748#endif
749 {
750
751 ElementContext elemCtx(simulator_);
752 auto elemIt = threadedElemIt.beginParallel();
753 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
754 if (elemIt->partitionType() != Dune::InteriorEntity) {
755 continue;
756 }
757 const Element& elem = *elemIt;
758 elemCtx.updatePrimaryStencil(elem);
759 // Mark cache for this element as invalid.
760 const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
761 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
762 const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
764 }
765 // Update for this element.
766 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
767 }
768 }
769 }
770
780 {
782 return;
783
784 if (enableStorageCache() && simulator_.problem().recycleFirstIterationStorage()) {
785 // If the storage term is cached, the intensive quantities of the previous
786 // time steps do not need to be accessed, and we can thus spare ourselves to
787 // copy the objects for the intensive quantities.
788 // However, if the storage term at the start of the timestep cannot be deduced
789 // from the primary variables, we must calculate it from the old intensive
790 // quantities, and need to shift them.
791 return;
792 }
793
794 assert(numSlots > 0);
795
796 for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++ timeIdx) {
797 intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
798 intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
799 }
800
801 // the cache for the most recent time indices do not need to be invalidated
802 // because the solution for them did not change (TODO: that assumes that there is
803 // no post-processing of the solution after a time step! fix it?)
804 }
805
813 { return enableStorageCache_; }
814
822 { enableStorageCache_= enableStorageCache; }
823
834 const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
835 {
836 assert(enableStorageCache_);
837 return storageCache_[timeIdx][globalIdx];
838 }
839
851 void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
852 {
853 assert(enableStorageCache_);
854 storageCache_[timeIdx][globalIdx] = value;
855 }
856
864 Scalar globalResidual(GlobalEqVector& dest,
865 const SolutionVector& u) const
866 {
867 SolutionVector tmp(asImp_().solution(/*timeIdx=*/0));
868 mutableSolution(/*timeIdx=*/0) = u;
869 Scalar res = asImp_().globalResidual(dest);
870 mutableSolution(/*timeIdx=*/0) = tmp;
871 return res;
872 }
873
880 Scalar globalResidual(GlobalEqVector& dest) const
881 {
882 dest = 0;
883
884 std::mutex mutex;
885 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
886#ifdef _OPENMP
887#pragma omp parallel
888#endif
889 {
890 // Attention: the variables below are thread specific and thus cannot be
891 // moved in front of the #pragma!
892 unsigned threadId = ThreadManager::threadId();
893 ElementContext elemCtx(simulator_);
894 ElementIterator elemIt = threadedElemIt.beginParallel();
895 LocalEvalBlockVector residual, storageTerm;
896
897 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
898 const Element& elem = *elemIt;
899 if (elem.partitionType() != Dune::InteriorEntity)
900 continue;
901
902 elemCtx.updateAll(elem);
903 residual.resize(elemCtx.numDof(/*timeIdx=*/0));
904 storageTerm.resize(elemCtx.numPrimaryDof(/*timeIdx=*/0));
905 asImp_().localResidual(threadId).eval(residual, elemCtx);
906
907 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
908 mutex.lock();
909 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
910 unsigned globalI = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
911 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
912 dest[globalI][eqIdx] += Toolbox::value(residual[dofIdx][eqIdx]);
913 }
914 mutex.unlock();
915 }
916 }
917
918 // add up the residuals on the process borders
919 const auto sumHandle =
920 GridCommHandleFactory::template sumHandle<EqVector>(dest, asImp_().dofMapper());
921 gridView_.communicate(*sumHandle,
922 Dune::InteriorBorder_InteriorBorder_Interface,
923 Dune::ForwardCommunication);
924
925 // calculate the square norm of the residual. this is not
926 // entirely correct, since the residual for the finite volumes
927 // which are on the boundary are counted once for every
928 // process. As often in life: shit happens (, we don't care)...
929 Scalar result2 = dest.two_norm2();
930 result2 = asImp_().gridView().comm().sum(result2);
931
932 return std::sqrt(result2);
933 }
934
941 void globalStorage(EqVector& storage, unsigned timeIdx = 0) const
942 {
943 storage = 0;
944
945 std::mutex mutex;
946 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
947#ifdef _OPENMP
948#pragma omp parallel
949#endif
950 {
951 // Attention: the variables below are thread specific and thus cannot be
952 // moved in front of the #pragma!
953 unsigned threadId = ThreadManager::threadId();
954 ElementContext elemCtx(simulator_);
955 ElementIterator elemIt = threadedElemIt.beginParallel();
956 LocalEvalBlockVector elemStorage;
957
958 // in this method, we need to disable the storage cache because we want to
959 // evaluate the storage term for other time indices than the most recent one
960 elemCtx.setEnableStorageCache(false);
961
962 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
963 const Element& elem = *elemIt;
964 if (elem.partitionType() != Dune::InteriorEntity)
965 continue; // ignore ghost and overlap elements
966
967 elemCtx.updateStencil(elem);
968 elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
969
970 size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
971 elemStorage.resize(numPrimaryDof);
972
973 localResidual(threadId).evalStorage(elemStorage, elemCtx, timeIdx);
974
975 mutex.lock();
976 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx)
977 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
978 storage[eqIdx] += Toolbox::value(elemStorage[dofIdx][eqIdx]);
979 mutex.unlock();
980 }
981 }
982
983 storage = gridView_.comm().sum(storage);
984 }
985
993 void checkConservativeness([[maybe_unused]] Scalar tolerance = -1,
994 [[maybe_unused]] bool verbose=false) const
995 {
996#ifndef NDEBUG
997 Scalar totalBoundaryArea(0.0);
998 Scalar totalVolume(0.0);
999 EvalEqVector totalRate(0.0);
1000
1001 // take the newton tolerance times the total volume of the grid if we're not
1002 // given an explicit tolerance...
1003 if (tolerance <= 0) {
1004 tolerance =
1005 simulator_.model().newtonMethod().tolerance()
1006 * simulator_.model().gridTotalVolume()
1007 * 1000;
1008 }
1009
1010 // we assume the implicit Euler time discretization for now...
1011 assert(historySize == 2);
1012
1013 EqVector storageBeginTimeStep(0.0);
1014 globalStorage(storageBeginTimeStep, /*timeIdx=*/1);
1015
1016 EqVector storageEndTimeStep(0.0);
1017 globalStorage(storageEndTimeStep, /*timeIdx=*/0);
1018
1019 // calculate the rate at the boundary and the source rate
1020 ElementContext elemCtx(simulator_);
1021 elemCtx.setEnableStorageCache(false);
1022 auto eIt = simulator_.gridView().template begin</*codim=*/0>();
1023 const auto& elemEndIt = simulator_.gridView().template end</*codim=*/0>();
1024 for (; eIt != elemEndIt; ++eIt) {
1025 if (eIt->partitionType() != Dune::InteriorEntity)
1026 continue; // ignore ghost and overlap elements
1027
1028 elemCtx.updateAll(*eIt);
1029
1030 // handle the boundary terms
1031 if (elemCtx.onBoundary()) {
1032 BoundaryContext boundaryCtx(elemCtx);
1033
1034 for (unsigned faceIdx = 0; faceIdx < boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); ++faceIdx) {
1035 BoundaryRateVector values;
1036 simulator_.problem().boundary(values,
1038 faceIdx,
1039 /*timeIdx=*/0);
1040 Valgrind::CheckDefined(values);
1041
1042 unsigned dofIdx = boundaryCtx.interiorScvIndex(faceIdx, /*timeIdx=*/0);
1043 const auto& insideIntQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1044
1045 Scalar bfArea =
1046 boundaryCtx.boundarySegmentArea(faceIdx, /*timeIdx=*/0)
1047 * insideIntQuants.extrusionFactor();
1048
1049 for (unsigned i = 0; i < values.size(); ++i)
1050 values[i] *= bfArea;
1051
1053 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
1054 totalRate[eqIdx] += values[eqIdx];
1055 }
1056 }
1057
1058 // deal with the source terms
1059 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
1060 RateVector values;
1061 simulator_.problem().source(values,
1062 elemCtx,
1063 dofIdx,
1064 /*timeIdx=*/0);
1065 Valgrind::CheckDefined(values);
1066
1067 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1068 Scalar dofVolume =
1069 elemCtx.dofVolume(dofIdx, /*timeIdx=*/0)
1070 * intQuants.extrusionFactor();
1071 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
1072 totalRate[eqIdx] += -dofVolume*Toolbox::value(values[eqIdx]);
1073 totalVolume += dofVolume;
1074 }
1075 }
1076
1077 // summarize everything over all processes
1078 const auto& comm = simulator_.gridView().comm();
1079 totalRate = comm.sum(totalRate);
1081 totalVolume = comm.sum(totalVolume);
1082
1083 if (comm.rank() == 0) {
1086 storageRate /= simulator_.timeStepSize();
1087 if (verbose) {
1088 std::cout << "storage at beginning of time step: " << storageBeginTimeStep << "\n";
1089 std::cout << "storage at end of time step: " << storageEndTimeStep << "\n";
1090 std::cout << "rate based on storage terms: " << storageRate << "\n";
1091 std::cout << "rate based on source and boundary terms: " << totalRate << "\n";
1092 std::cout << "difference in rates: ";
1093 for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx)
1094 std::cout << (storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) << " ";
1095 std::cout << "\n";
1096 }
1097 for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
1098 Scalar eps =
1099 (std::abs(storageRate[eqIdx]) + Toolbox::value(totalRate[eqIdx]))*tolerance;
1100 eps = std::max(tolerance, eps);
1101 assert(std::abs(storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) <= eps);
1102 }
1103 }
1104#endif // NDEBUG
1105 }
1106
1112 Scalar dofTotalVolume(unsigned globalIdx) const
1113 { return dofTotalVolume_[globalIdx]; }
1114
1120 bool isLocalDof(unsigned globalIdx) const
1121 { return isLocalDof_[globalIdx]; }
1122
1127 Scalar gridTotalVolume() const
1128 { return gridTotalVolume_; }
1129
1135 const SolutionVector& solution(unsigned timeIdx) const
1136 { return solution_[timeIdx]->blockVector(); }
1137
1141 SolutionVector& solution(unsigned timeIdx)
1142 { return solution_[timeIdx]->blockVector(); }
1143
1144 protected:
1148 SolutionVector& mutableSolution(unsigned timeIdx) const
1149 { return solution_[timeIdx]->blockVector(); }
1150
1151 public:
1156 const Linearizer& linearizer() const
1157 { return *linearizer_; }
1158
1163 Linearizer& linearizer()
1164 { return *linearizer_; }
1165
1174 const LocalLinearizer& localLinearizer(unsigned openMpThreadId) const
1175 { return localLinearizer_[openMpThreadId]; }
1179 LocalLinearizer& localLinearizer(unsigned openMpThreadId)
1180 { return localLinearizer_[openMpThreadId]; }
1181
1185 const LocalResidual& localResidual(unsigned openMpThreadId) const
1186 { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1190 LocalResidual& localResidual(unsigned openMpThreadId)
1191 { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1192
1200 Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
1201 {
1202 Scalar absPv = std::abs(asImp_().solution(/*timeIdx=*/1)[globalDofIdx][pvIdx]);
1203 return 1.0/std::max(absPv, 1.0);
1204 }
1205
1212 Scalar eqWeight(unsigned, unsigned) const
1213 { return 1.0; }
1214
1225 const PrimaryVariables& pv1,
1226 const PrimaryVariables& pv2) const
1227 {
1228 Scalar result = 0.0;
1229 for (unsigned j = 0; j < numEq; ++j) {
1230 Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
1231 Scalar eqErr = std::abs((pv1[j] - pv2[j])*weight);
1232 //Scalar eqErr = std::abs(pv1[j] - pv2[j]);
1233 //eqErr *= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
1234
1235 result = std::max(result, eqErr);
1236 }
1237 return result;
1238 }
1239
1245 bool update()
1246 {
1247 TimerGuard prePostProcessGuard(prePostProcessTimer_);
1248
1249#ifndef NDEBUG
1250 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1251 // Make sure that the primary variables are defined. Note that because of padding
1252 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1253 // for definedness...
1254 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1255 asImp_().solution(timeIdx)[i].checkDefined();
1256 }
1257 }
1258#endif // NDEBUG
1259
1260 // make sure all timers are prestine
1261 prePostProcessTimer_.halt();
1262 linearizeTimer_.halt();
1263 solveTimer_.halt();
1264 updateTimer_.halt();
1265
1266 prePostProcessTimer_.start();
1267 asImp_().updateBegin();
1268 prePostProcessTimer_.stop();
1269
1270 bool converged = false;
1271
1272 try {
1273 converged = newtonMethod_.apply();
1274 }
1275 catch(...) {
1276 prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1277 linearizeTimer_ += newtonMethod_.linearizeTimer();
1278 solveTimer_ += newtonMethod_.solveTimer();
1279 updateTimer_ += newtonMethod_.updateTimer();
1280
1281 throw;
1282 }
1283
1284#ifndef NDEBUG
1285 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1286 // Make sure that the primary variables are defined. Note that because of padding
1287 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1288 // for definedness...
1289 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1290 asImp_().solution(timeIdx)[i].checkDefined();
1291 }
1292 }
1293#endif // NDEBUG
1294
1295 prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1296 linearizeTimer_ += newtonMethod_.linearizeTimer();
1297 solveTimer_ += newtonMethod_.solveTimer();
1298 updateTimer_ += newtonMethod_.updateTimer();
1299
1300 prePostProcessTimer_.start();
1301 if (converged)
1302 asImp_().updateSuccessful();
1303 else
1304 asImp_().updateFailed();
1305 prePostProcessTimer_.stop();
1306
1307#ifndef NDEBUG
1308 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1309 // Make sure that the primary variables are defined. Note that because of padding
1310 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1311 // for definedness...
1312 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1313 asImp_().solution(timeIdx)[i].checkDefined();
1314 }
1315 }
1316#endif // NDEBUG
1317
1318 return converged;
1319 }
1320
1329 { }
1330
1337 { }
1338
1344 { }
1345
1350 {
1351 throw std::invalid_argument("Grid adaptation need to be implemented for "
1352 "specific settings of grid and function spaces");
1353 }
1354
1361 {
1362 // Reset the current solution to the one of the
1363 // previous time step so that we can start the next
1364 // update at a physically meaningful solution.
1365 solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
1366 invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1367
1368#ifndef NDEBUG
1369 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1370 // Make sure that the primary variables are defined. Note that because of padding
1371 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1372 // for definedness...
1373 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i)
1374 asImp_().solution(timeIdx)[i].checkDefined();
1375 }
1376#endif // NDEBUG
1377 }
1378
1387 {
1388 // at this point we can adapt the grid
1389 if (this->enableGridAdaptation_) {
1390 asImp_().adaptGrid();
1391 }
1392
1393 // make the current solution the previous one.
1394 solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
1395
1396 // shift the intensive quantities cache by one position in the
1397 // history
1398 asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
1399 }
1400
1408 template <class Restarter>
1410 {
1411 throw std::runtime_error("Not implemented: The discretization chosen for this problem "
1412 "does not support restart files. (serialize() method unimplemented)");
1413 }
1414
1422 template <class Restarter>
1424 {
1425 throw std::runtime_error("Not implemented: The discretization chosen for this problem "
1426 "does not support restart files. (deserialize() method unimplemented)");
1427 }
1428
1437 template <class DofEntity>
1438 void serializeEntity(std::ostream& outstream,
1439 const DofEntity& dof)
1440 {
1441 unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1442
1443 // write phase state
1444 if (!outstream.good()) {
1445 throw std::runtime_error("Could not serialize degree of freedom "
1446 +std::to_string(dofIdx));
1447 }
1448
1449 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1450 outstream << solution(/*timeIdx=*/0)[dofIdx][eqIdx] << " ";
1451 }
1452 }
1453
1462 template <class DofEntity>
1463 void deserializeEntity(std::istream& instream,
1464 const DofEntity& dof)
1465 {
1466 unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1467
1468 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1469 if (!instream.good())
1470 throw std::runtime_error("Could not deserialize degree of freedom "
1471 +std::to_string(dofIdx));
1472 instream >> solution(/*timeIdx=*/0)[dofIdx][eqIdx];
1473 }
1474 }
1475
1479 size_t numGridDof() const
1480 { throw std::logic_error("The discretization class must implement the numGridDof() method!"); }
1481
1485 size_t numAuxiliaryDof() const
1486 {
1487 size_t result = 0;
1488 auto auxModIt = auxEqModules_.begin();
1489 const auto& auxModEndIt = auxEqModules_.end();
1490 for (; auxModIt != auxModEndIt; ++auxModIt)
1491 result += (*auxModIt)->numDofs();
1492
1493 return result;
1494 }
1495
1499 size_t numTotalDof() const
1500 { return asImp_().numGridDof() + numAuxiliaryDof(); }
1501
1506 const DofMapper& dofMapper() const
1507 { throw std::logic_error("The discretization class must implement the dofMapper() method!"); }
1508
1512 const VertexMapper& vertexMapper() const
1513 { return vertexMapper_; }
1514
1518 const ElementMapper& elementMapper() const
1519 { return elementMapper_; }
1520
1526 {
1527 delete linearizer_;
1528 linearizer_ = new Linearizer;
1529 linearizer_->init(simulator_);
1530 }
1531
1535 static std::string discretizationName()
1536 { return ""; }
1537
1543 std::string primaryVarName(unsigned pvIdx) const
1544 {
1545 std::ostringstream oss;
1546 oss << "primary variable_" << pvIdx;
1547 return oss.str();
1548 }
1549
1555 std::string eqName(unsigned eqIdx) const
1556 {
1557 std::ostringstream oss;
1558 oss << "equation_" << eqIdx;
1559 return oss.str();
1560 }
1561
1568 void updatePVWeights(const ElementContext&) const
1569 { }
1570
1575 { outputModules_.push_back(newModule); }
1576
1585 template <class VtkMultiWriter>
1587 const SolutionVector& u,
1588 const GlobalEqVector& deltaU) const
1589 {
1590 using ScalarBuffer = std::vector<double>;
1591
1592 GlobalEqVector globalResid(u.size());
1593 asImp_().globalResidual(globalResid, u);
1594
1595 // create the required scalar fields
1596 size_t numGridDof = asImp_().numGridDof();
1597
1598 // global defect of the two auxiliary equations
1599 ScalarBuffer* def[numEq];
1600 ScalarBuffer* delta[numEq];
1601 ScalarBuffer* priVars[numEq];
1602 ScalarBuffer* priVarWeight[numEq];
1603 ScalarBuffer* relError = writer.allocateManagedScalarBuffer(numGridDof);
1604 ScalarBuffer* normalizedRelError = writer.allocateManagedScalarBuffer(numGridDof);
1605 for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1606 priVars[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1607 priVarWeight[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1608 delta[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1609 def[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1610 }
1611
1612 Scalar minRelErr = 1e30;
1613 Scalar maxRelErr = -1e30;
1614 for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx) {
1615 for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1616 (*priVars[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
1617 (*priVarWeight[pvIdx])[globalIdx] = asImp_().primaryVarWeight(globalIdx, pvIdx);
1618 (*delta[pvIdx])[globalIdx] = - deltaU[globalIdx][pvIdx];
1619 (*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
1620 }
1621
1622 PrimaryVariables uOld(u[globalIdx]);
1623 PrimaryVariables uNew(uOld);
1624 uNew -= deltaU[globalIdx];
1625
1626 Scalar err = asImp_().relativeDofError(globalIdx, uOld, uNew);
1627 (*relError)[globalIdx] = err;
1628 (*normalizedRelError)[globalIdx] = err;
1629 minRelErr = std::min(err, minRelErr);
1630 maxRelErr = std::max(err, maxRelErr);
1631 }
1632
1633 // do the normalization of the relative error
1634 Scalar alpha = std::max(Scalar{1e-20},
1635 std::max(std::abs(maxRelErr),
1636 std::abs(minRelErr)));
1637 for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx)
1638 (*normalizedRelError)[globalIdx] /= alpha;
1639
1640 DiscBaseOutputModule::attachScalarDofData_(writer, *relError, "relative error");
1641 DiscBaseOutputModule::attachScalarDofData_(writer, *normalizedRelError, "normalized relative error");
1642
1643 for (unsigned i = 0; i < numEq; ++i) {
1644 std::ostringstream oss;
1645 oss.str(""); oss << "priVar_" << asImp_().primaryVarName(i);
1646 DiscBaseOutputModule::attachScalarDofData_(writer,
1647 *priVars[i],
1648 oss.str());
1649
1650 oss.str(""); oss << "delta_" << asImp_().primaryVarName(i);
1651 DiscBaseOutputModule::attachScalarDofData_(writer,
1652 *delta[i],
1653 oss.str());
1654
1655 oss.str(""); oss << "weight_" << asImp_().primaryVarName(i);
1656 DiscBaseOutputModule::attachScalarDofData_(writer,
1657 *priVarWeight[i],
1658 oss.str());
1659
1660 oss.str(""); oss << "defect_" << asImp_().eqName(i);
1661 DiscBaseOutputModule::attachScalarDofData_(writer,
1662 *def[i],
1663 oss.str());
1664 }
1665
1666 asImp_().prepareOutputFields();
1667 asImp_().appendOutputFields(writer);
1668 }
1669
1675 {
1676 bool needFullContextUpdate = false;
1677 auto modIt = outputModules_.begin();
1678 const auto& modEndIt = outputModules_.end();
1679 for (; modIt != modEndIt; ++modIt) {
1680 (*modIt)->allocBuffers();
1681 needFullContextUpdate = needFullContextUpdate || (*modIt)->needExtensiveQuantities();
1682 }
1683
1684 // iterate over grid
1685 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
1686#ifdef _OPENMP
1687#pragma omp parallel
1688#endif
1689 {
1690 ElementContext elemCtx(simulator_);
1691 ElementIterator elemIt = threadedElemIt.beginParallel();
1692 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
1693 const auto& elem = *elemIt;
1694 if (elem.partitionType() != Dune::InteriorEntity)
1695 // ignore non-interior entities
1696 continue;
1697
1699 elemCtx.updateAll(elem);
1700 else {
1701 elemCtx.updatePrimaryStencil(elem);
1702 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1703 }
1704
1705 // we cannot reuse the "modIt" variable here because the code here might
1706 // be threaded and "modIt" is is the same for all threads, i.e., if a
1707 // given thread modifies it, the changes affect all threads.
1708 auto modIt2 = outputModules_.begin();
1709 for (; modIt2 != modEndIt; ++modIt2)
1710 (*modIt2)->processElement(elemCtx);
1711 }
1712 }
1713 }
1714
1720 {
1721 auto modIt = outputModules_.begin();
1722 const auto& modEndIt = outputModules_.end();
1723 for (; modIt != modEndIt; ++modIt)
1724 (*modIt)->commitBuffers(writer);
1725 }
1726
1730 const GridView& gridView() const
1731 { return gridView_; }
1732
1745 {
1746 auxMod->setDofOffset(numTotalDof());
1747 auxEqModules_.push_back(auxMod);
1748
1749 // resize the solutions
1750 if (enableGridAdaptation_
1751 && !std::is_same<DiscreteFunction, BlockVectorWrapper>::value)
1752 {
1753 throw std::invalid_argument("Problems which require auxiliary modules cannot be used in"
1754 " conjunction with dune-fem");
1755 }
1756
1757 size_t numDof = numTotalDof();
1758 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx)
1759 solution(timeIdx).resize(numDof);
1760
1761 auxMod->applyInitial();
1762 }
1763
1770 {
1771 auxEqModules_.clear();
1772 linearizer_->eraseMatrix();
1773 newtonMethod_.eraseMatrix();
1774 }
1775
1779 size_t numAuxiliaryModules() const
1780 { return auxEqModules_.size(); }
1781
1786 { return auxEqModules_[auxEqModIdx]; }
1787
1792 { return auxEqModules_[auxEqModIdx]; }
1793
1798 { return enableIntensiveQuantityCache_ || enableThermodynamicHints_; }
1799
1800 const Timer& prePostProcessTimer() const
1801 { return prePostProcessTimer_; }
1802
1803 const Timer& linearizeTimer() const
1804 { return linearizeTimer_; }
1805
1806 const Timer& solveTimer() const
1807 { return solveTimer_; }
1808
1809 const Timer& updateTimer() const
1810 { return updateTimer_; }
1811
1812 template<class Serializer>
1813 void serializeOp(Serializer& serializer)
1814 {
1816 using Helper = typename BaseDiscretization::template SerializeHelper<Serializer>;
1817 Helper::serializeOp(serializer, solution_);
1818 }
1819
1820 bool operator==(const FvBaseDiscretization& rhs) const
1821 {
1822 return std::equal(this->solution_.begin(), this->solution_.end(),
1823 rhs.solution_.begin(), rhs.solution_.end(),
1824 [](const auto& x, const auto& y)
1825 {
1826 return *x == *y;
1827 });
1828 }
1829
1830protected:
1831 void resizeAndResetIntensiveQuantitiesCache_()
1832 {
1833 // allocate the storage cache
1834 if (enableStorageCache()) {
1835 size_t numDof = asImp_().numGridDof();
1836 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1837 storageCache_[timeIdx].resize(numDof);
1838 }
1839 }
1840
1841 // allocate the intensive quantities cache
1843 size_t numDof = asImp_().numGridDof();
1844 for(unsigned timeIdx=0; timeIdx<historySize; ++timeIdx) {
1845 intensiveQuantityCache_[timeIdx].resize(numDof);
1846 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
1848 }
1849 }
1850 }
1851 template <class Context>
1852 void supplementInitialSolution_(PrimaryVariables&,
1853 const Context&,
1854 unsigned,
1855 unsigned)
1856 { }
1857
1866 {
1867 // add the output modules available on all model
1868 auto *mod = new VtkPrimaryVarsModule<TypeTag>(simulator_);
1869 this->outputModules_.push_back(mod);
1870 }
1871
1875 LocalResidual& localResidual_()
1876 { return localLinearizer_.localResidual(); }
1877
1881 bool verbose_() const
1882 { return gridView_.comm().rank() == 0; }
1883
1884 Implementation& asImp_()
1885 { return *static_cast<Implementation*>(this); }
1886 const Implementation& asImp_() const
1887 { return *static_cast<const Implementation*>(this); }
1888
1889 // the problem we want to solve. defines the constitutive
1890 // relations, matxerial laws, etc.
1891 Simulator& simulator_;
1892
1893 // the representation of the spatial domain of the problem
1894 GridView gridView_;
1895
1896 // the mappers for element and vertex entities to global indices
1897 ElementMapper elementMapper_;
1898 VertexMapper vertexMapper_;
1899
1900 // a vector with all auxiliary equations to be considered
1901 std::vector<BaseAuxiliaryModule<TypeTag>*> auxEqModules_;
1902
1903 NewtonMethod newtonMethod_;
1904
1905 Timer prePostProcessTimer_;
1906 Timer linearizeTimer_;
1907 Timer solveTimer_;
1908 Timer updateTimer_;
1909
1910 // calculates the local jacobian matrix for a given element
1911 std::vector<LocalLinearizer> localLinearizer_;
1912 // Linearizes the problem at the current time step using the
1913 // local jacobian
1914 Linearizer *linearizer_;
1915
1916 // cur is the current iterative solution, prev the converged
1917 // solution of the previous time step
1918 mutable IntensiveQuantitiesVector intensiveQuantityCache_[historySize];
1919 // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
1920 mutable std::vector<unsigned char> intensiveQuantityCacheUpToDate_[historySize];
1921
1922 mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
1923
1924 std::list<BaseOutputModule<TypeTag>*> outputModules_;
1925
1926 Scalar gridTotalVolume_;
1927 std::vector<Scalar> dofTotalVolume_;
1928 std::vector<bool> isLocalDof_;
1929
1930 mutable GlobalEqVector storageCache_[historySize];
1931
1932 bool enableGridAdaptation_;
1933 bool enableIntensiveQuantityCache_;
1934 bool enableStorageCache_;
1935 bool enableThermodynamicHints_;
1936};
1937
1943template<class TypeTag>
1945{
1949
1950 static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
1951
1952public:
1953 template<class Serializer>
1955 template<class SolutionType>
1956 static void serializeOp(Serializer& serializer,
1958 {
1959 for (auto& sol : solution) {
1960 serializer(*sol);
1961 }
1962 }
1963 };
1964
1965 explicit FvBaseDiscretizationNoAdapt(Simulator& simulator)
1966 : ParentType(simulator)
1967 {
1968 if (this->enableGridAdaptation_) {
1969 throw std::invalid_argument("Grid adaptation need to use"
1970 " BaseDiscretization = FvBaseDiscretizationFemAdapt"
1971 " which currently requires the presence of the"
1972 " dune-fem module");
1973 }
1974 size_t numDof = this->asImp_().numGridDof();
1975 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1976 this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", numDof);
1977 }
1978 }
1979};
1980
1981} // namespace Opm
1982
1983#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.
Base class for specifying auxiliary equations.
Definition baseauxiliarymodule.hh:56
The base class for writer modules.
Definition baseoutputmodule.hh:67
The base class for all output writers.
Definition baseoutputwriter.hh:44
Represents all quantities which available on boundary segments.
Definition fvbaseboundarycontext.hh:44
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:46
The base class for the finite volume discretization schemes without adaptation.
Definition fvbasediscretization.hh:1945
Definition fvbasediscretization.hh:343
The base class for the finite volume discretization schemes.
Definition fvbasediscretization.hh:293
LocalLinearizer & localLinearizer(unsigned openMpThreadId)
Definition fvbasediscretization.hh:1179
size_t numAuxiliaryDof() const
Returns the number of degrees of freedom (DOFs) of the auxiliary equations.
Definition fvbasediscretization.hh:1485
void prepareOutputFields() const
Prepare the quantities relevant for the current solution to be appended to the output writers.
Definition fvbasediscretization.hh:1674
void shiftIntensiveQuantityCache(unsigned numSlots=1)
Move the intensive quantities for a given time index to the back.
Definition fvbasediscretization.hh:779
void adaptGrid()
Called by the update() method when the grid should be refined.
Definition fvbasediscretization.hh:1349
void addAuxiliaryModule(BaseAuxiliaryModule< TypeTag > *auxMod)
Add a module for an auxiliary equation.
Definition fvbasediscretization.hh:1744
void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx, unsigned timeIdx, bool newValue) const
Invalidate the cache for a given intensive quantities object.
Definition fvbasediscretization.hh:697
void finishInit()
Apply the initial conditions to the model.
Definition fvbasediscretization.hh:472
void prefetch(const Element &) const
Allows to improve the performance by prefetching all data which is associated with a given element.
Definition fvbasediscretization.hh:601
void updateSuccessful()
Called by the update() method if it was successful.
Definition fvbasediscretization.hh:1343
static std::string discretizationName()
Returns a string of discretization's human-readable name.
Definition fvbasediscretization.hh:1535
BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx)
Returns a given module for auxiliary equations.
Definition fvbasediscretization.hh:1785
bool isLocalDof(unsigned globalIdx) const
Returns if the overlap of the volume ofa degree of freedom is non-zero.
Definition fvbasediscretization.hh:1120
LocalResidual & localResidual_()
Reference to the local residal object.
Definition fvbasediscretization.hh:1875
void registerOutputModules_()
Register all output modules which make sense for the model.
Definition fvbasediscretization.hh:1865
bool enableGridAdaptation() const
Returns whether the grid ought to be adapted to the solution during the simulation.
Definition fvbasediscretization.hh:540
const NewtonMethod & newtonMethod() const
Returns the newton method object.
Definition fvbasediscretization.hh:615
SolutionVector & mutableSolution(unsigned timeIdx) const
Definition fvbasediscretization.hh:1148
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition fvbasediscretization.hh:1386
NewtonMethod & newtonMethod()
Returns the newton method object.
Definition fvbasediscretization.hh:609
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition fvbasediscretization.hh:1512
const EqVector & cachedStorage(unsigned globalIdx, unsigned timeIdx) const
Retrieve an entry of the cache for the storage term.
Definition fvbasediscretization.hh:834
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:1568
void updateFailed()
Called by the update() method if it was unsuccessful.
Definition fvbasediscretization.hh:1360
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition fvbasediscretization.hh:1336
Scalar globalResidual(GlobalEqVector &dest) const
Compute the global residual for the current solution vector.
Definition fvbasediscretization.hh:880
LocalResidual & localResidual(unsigned openMpThreadId)
Definition fvbasediscretization.hh:1190
void serializeEntity(std::ostream &outstream, const DofEntity &dof)
Write the current solution for a degree of freedom to a restart file.
Definition fvbasediscretization.hh:1438
void setEnableStorageCache(bool enableStorageCache)
Set the value of enable storage cache.
Definition fvbasediscretization.hh:821
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition fvbasediscretization.hh:1555
Scalar eqWeight(unsigned, unsigned) const
Returns the relative weight of an equation.
Definition fvbasediscretization.hh:1212
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition fvbasediscretization.hh:1200
void deserialize(Restarter &)
Deserializes the state of the model.
Definition fvbasediscretization.hh:1423
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:993
Scalar gridTotalVolume() const
Returns the volume of the whole grid which represents the spatial domain.
Definition fvbasediscretization.hh:1127
const IntensiveQuantities * thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
Return the thermodynamic hint for a entity on the grid at given time.
Definition fvbasediscretization.hh:633
const Linearizer & linearizer() const
Returns the operator linearizer for the global jacobian of the problem.
Definition fvbasediscretization.hh:1156
void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector &value) const
Set an entry of the cache for the storage term.
Definition fvbasediscretization.hh:851
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:1586
Scalar relativeDofError(unsigned vertexIdx, const PrimaryVariables &pv1, const PrimaryVariables &pv2) const
Returns the relative error between two vectors of primary variables.
Definition fvbasediscretization.hh:1224
Scalar globalResidual(GlobalEqVector &dest, const SolutionVector &u) const
Compute the global residual for an arbitrary solution vector.
Definition fvbasediscretization.hh:864
SolutionVector & solution(unsigned timeIdx)
Definition fvbasediscretization.hh:1141
void addOutputModule(BaseOutputModule< TypeTag > *newModule)
Add an module for writing visualization output after a timestep.
Definition fvbasediscretization.hh:1574
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:1463
const LocalLinearizer & localLinearizer(unsigned openMpThreadId) const
Returns the local jacobian which calculates the local stiffness matrix for an arbitrary element.
Definition fvbasediscretization.hh:1174
bool update()
Try to progress the model to the next timestep.
Definition fvbasediscretization.hh:1245
void clearAuxiliaryModules()
Causes the list of auxiliary equations to be cleared.
Definition fvbasediscretization.hh:1769
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition fvbasediscretization.hh:1518
size_t numGridDof() const
Returns the number of degrees of freedom (DOFs) for the computational grid.
Definition fvbasediscretization.hh:1479
void syncOverlap()
Syncronize the values of the primary variables on the degrees of freedom that overlap with the neighb...
Definition fvbasediscretization.hh:1328
void appendOutputFields(BaseOutputWriter &writer) const
Append the quantities relevant for the current solution to an output writer.
Definition fvbasediscretization.hh:1719
void applyInitialSolution()
Applies the initial solution for all degrees of freedom to which the model applies.
Definition fvbasediscretization.hh:547
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:679
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition fvbasediscretization.hh:1543
void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
Invalidate the whole intensive quantity cache for time index.
Definition fvbasediscretization.hh:712
Linearizer & linearizer()
Returns the object which linearizes the global system of equations at the current solution.
Definition fvbasediscretization.hh:1163
const BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx) const
Returns a given module for auxiliary equations.
Definition fvbasediscretization.hh:1791
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:941
Scalar dofTotalVolume(unsigned globalIdx) const
Returns the volume of a given control volume.
Definition fvbasediscretization.hh:1112
static void registerParameters()
Register all run-time parameters for the model.
Definition fvbasediscretization.hh:441
const SolutionVector & solution(unsigned timeIdx) const
Reference to the solution at a given history index as a block vector.
Definition fvbasediscretization.hh:1135
bool verbose_() const
Returns whether messages should be printed.
Definition fvbasediscretization.hh:1881
const LocalResidual & localResidual(unsigned openMpThreadId) const
Returns the object to calculate the local residual function.
Definition fvbasediscretization.hh:1185
const DofMapper & dofMapper() const
Mapper to convert the Dune entities of the discretization's degrees of freedoms are to indices.
Definition fvbasediscretization.hh:1506
bool storeIntensiveQuantities() const
Returns true if the cache for intensive quantities is enabled.
Definition fvbasediscretization.hh:1797
void serialize(Restarter &)
Serializes the current state of the model.
Definition fvbasediscretization.hh:1409
size_t numTotalDof() const
Returns the total number of degrees of freedom (i.e., grid plux auxiliary DOFs)
Definition fvbasediscretization.hh:1499
void resetLinearizer()
Resets the Jacobian matrix linearizer, so that the boundary types can be altered.
Definition fvbasediscretization.hh:1525
bool enableStorageCache() const
Returns true iff the storage term is cached.
Definition fvbasediscretization.hh:812
size_t numAuxiliaryModules() const
Returns the number of modules for auxiliary equations.
Definition fvbasediscretization.hh:1779
const GridView & gridView() const
Reference to the grid view of the spatial domain.
Definition fvbasediscretization.hh:1730
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:653
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
Definition fvbaseelementcontext.hh:52
Provide the properties at a face which make sense indepentently of the conserved quantities.
Definition fvbaseextensivequantities.hh:46
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition fvbasegradientcalculator.hh:47
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:71
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition fvbaselocalresidual.hh:58
Represents the primary variables used by the a model.
Definition fvbaseprimaryvariables.hh:52
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:129
Manages the initializing and running of time dependent problems.
Definition simulator.hh:97
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
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
Definition threadedentityiterator.hh:43
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition timerguard.hh:41
Provides an encapsulation to measure the system time.
Definition timer.hpp:46
void start()
Start counting the time resources used by the simulation.
Definition timer.cpp:46
void halt()
Stop the measurement reset all timing values.
Definition timer.cpp:75
double stop()
Stop counting the time resources.
Definition timer.cpp:52
Simplifies writing multi-file VTK datasets.
Definition vtkmultiwriter.hh:66
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkprimaryvarsmodule.hpp:73
Definition alignedallocator.hh:94
Calculates the local residual and its Jacobian for a single element of the grid.
Represents all quantities which available on boundary segments.
Class to specify constraints for a finite volume spatial discretization.
Represents all quantities which available for calculating constraints.
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
Provide the properties at a face which make sense indepentently of the conserved quantities.
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Base class for the model specific class which provides access to all intensive (i....
The common code for the linearizers of non-linear systems of equations.
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
A Newton method for models using a finite volume discretization.
Represents the primary variables used by the a model.
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.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
This is a grid manager which does not create any border list.
Manages the initializing and running of time dependent problems.
Definition fvbasediscretization.hh:1954
Definition fvbasediscretization.hh:267
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
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
Definition fvbaseproperties.hh:77
The secondary variables of all degrees of freedom in an element's stencil.
Definition fvbaseproperties.hh:140
A vector of holding a quantity for each equation for each DOF of an element.
Definition fvbaseproperties.hh:112
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
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
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
Vector containing volumetric or areal rates of quantities.
Definition fvbaseproperties.hh:116
Manages the simulation time.
Definition basicproperties.hh:116
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:40
use locking to prevent race conditions when linearizing the global system of equations in multi-threa...
Definition fvbaseproperties.hh:181
Specify whether to use volumetric residuals or not.
Definition fvbaseproperties.hh:237
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
Simplifies multi-threaded capabilities.
Provides an encapsulation to measure the system time.
A simple class which makes sure that a timer gets stopped if an exception is thrown.
VTK output module for the fluid composition.