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 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;
502 gridTotalVolume_ += dofVolume;
503 }
504 }
505
506 // determine which DOFs should be considered to lie fully in the interior of the
507 // local process grid partition: those which do not have a non-zero volume
508 // before taking the peer processes into account...
509 isLocalDof_.resize(numDof);
510 for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx)
511 isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0);
512
513 // add the volumes of the DOFs on the process boundaries
514 const auto sumHandle =
515 GridCommHandleFactory::template sumHandle<Scalar>(dofTotalVolume_,
516 asImp_().dofMapper());
517 gridView_.communicate(*sumHandle,
518 Dune::InteriorBorder_All_Interface,
519 Dune::ForwardCommunication);
520
521 // sum up the volumes of the grid partitions
522 gridTotalVolume_ = gridView_.comm().sum(gridTotalVolume_);
523
524 linearizer_->init(simulator_);
525 for (unsigned threadId = 0; threadId < ThreadManager::maxThreads(); ++threadId)
526 localLinearizer_[threadId].init(simulator_);
527
528 resizeAndResetIntensiveQuantitiesCache_();
530 // invalidate all cached intensive quantities
531 for (unsigned timeIdx = 0; timeIdx < historySize; ++ timeIdx)
533 }
534
535 newtonMethod_.finishInit();
536 }
537
542 { return enableGridAdaptation_; }
543
549 {
550 // first set the whole domain to zero
551 SolutionVector& uCur = asImp_().solution(/*timeIdx=*/0);
552 uCur = Scalar(0.0);
553
554 ElementContext elemCtx(simulator_);
555
556 // iterate through the grid and evaluate the initial condition
557 for (const auto& elem : elements(gridView_)) {
558 // ignore everything which is not in the interior if the
559 // current process' piece of the grid
560 if (elem.partitionType() != Dune::InteriorEntity)
561 continue;
562
563 // deal with the current element
564 elemCtx.updateStencil(elem);
565
566 // loop over all element vertices, i.e. sub control volumes
567 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++)
568 {
569 // map the local degree of freedom index to the global one
570 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
571
572 // let the problem do the dirty work of nailing down
573 // the initial solution.
574 simulator_.problem().initial(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
575 asImp_().supplementInitialSolution_(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
576 uCur[globalIdx].checkDefined();
577 }
578 }
579
580 // synchronize the ghost DOFs (if necessary)
581 asImp_().syncOverlap();
582
583 // also set the solutions of the "previous" time steps to the initial solution.
584 for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx)
585 solution(timeIdx) = solution(/*timeIdx=*/0);
586
587 simulator_.problem().initialSolutionApplied();
588
589#ifndef NDEBUG
590 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
591 const auto& sol = solution(timeIdx);
592 for (unsigned dofIdx = 0; dofIdx < sol.size(); ++dofIdx)
593 sol[dofIdx].checkDefined();
594 }
595#endif // NDEBUG
596 }
597
602 void prefetch(const Element&) const
603 {
604 // do nothing by default
605 }
606
610 NewtonMethod& newtonMethod()
611 { return newtonMethod_; }
612
616 const NewtonMethod& newtonMethod() const
617 { return newtonMethod_; }
618
634 const IntensiveQuantities* thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
635 {
636 if (!enableThermodynamicHints_)
637 return 0;
638
639 // the intensive quantities cache doubles as thermodynamic hint
640 return cachedIntensiveQuantities(globalIdx, timeIdx);
641 }
642
654 const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
655 {
656 if (!enableIntensiveQuantityCache_ || !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
657 return nullptr;
658 }
659
660 // With the storage cache enabled, usually only the
661 // intensive quantities for the most recent time step are
662 // cached. However, this may be false for some Problem
663 // variants, so we should check if the cache exists for
664 // the timeIdx in question.
665 if (timeIdx > 0 && enableStorageCache_ && intensiveQuantityCache_[timeIdx].empty()) {
666 return nullptr;
667 }
668
669 return &intensiveQuantityCache_[timeIdx][globalIdx];
670 }
671
680 void updateCachedIntensiveQuantities(const IntensiveQuantities& intQuants,
681 unsigned globalIdx,
682 unsigned timeIdx) const
683 {
685 return;
686
687 intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
688 intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = 1;
689 }
690
699 unsigned timeIdx,
700 bool newValue) const
701 {
703 return;
704
705 intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue ? 1 : 0;
706 }
707
714 {
716 std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
717 intensiveQuantityCacheUpToDate_[timeIdx].end(),
718 /*value=*/0);
719 }
720 }
721
722 void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
723 {
725
726 // loop over all elements...
727 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
728#ifdef _OPENMP
729#pragma omp parallel
730#endif
731 {
732 ElementContext elemCtx(simulator_);
733 ElementIterator elemIt = threadedElemIt.beginParallel();
734 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
735 const Element& elem = *elemIt;
736 elemCtx.updatePrimaryStencil(elem);
737 elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
738 }
739 }
740 }
741
742 template <class GridViewType>
743 void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
744 {
745 // loop over all elements...
746 ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(gridView);
747#ifdef _OPENMP
748#pragma omp parallel
749#endif
750 {
751
752 ElementContext elemCtx(simulator_);
753 auto elemIt = threadedElemIt.beginParallel();
754 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
755 if (elemIt->partitionType() != Dune::InteriorEntity) {
756 continue;
757 }
758 const Element& elem = *elemIt;
759 elemCtx.updatePrimaryStencil(elem);
760 // Mark cache for this element as invalid.
761 const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
762 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
763 const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
765 }
766 // Update for this element.
767 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
768 }
769 }
770 }
771
781 {
783 return;
784
785 if (enableStorageCache() && simulator_.problem().recycleFirstIterationStorage()) {
786 // If the storage term is cached, the intensive quantities of the previous
787 // time steps do not need to be accessed, and we can thus spare ourselves to
788 // copy the objects for the intensive quantities.
789 // However, if the storage term at the start of the timestep cannot be deduced
790 // from the primary variables, we must calculate it from the old intensive
791 // quantities, and need to shift them.
792 return;
793 }
794
795 assert(numSlots > 0);
796
797 for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++ timeIdx) {
798 intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
799 intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
800 }
801
802 // the cache for the most recent time indices do not need to be invalidated
803 // because the solution for them did not change (TODO: that assumes that there is
804 // no post-processing of the solution after a time step! fix it?)
805 }
806
814 { return enableStorageCache_; }
815
823 { enableStorageCache_= enableStorageCache; }
824
835 const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
836 {
837 assert(enableStorageCache_);
838 return storageCache_[timeIdx][globalIdx];
839 }
840
852 void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
853 {
854 assert(enableStorageCache_);
855 storageCache_[timeIdx][globalIdx] = value;
856 }
857
865 Scalar globalResidual(GlobalEqVector& dest,
866 const SolutionVector& u) const
867 {
868 SolutionVector tmp(asImp_().solution(/*timeIdx=*/0));
869 mutableSolution(/*timeIdx=*/0) = u;
870 Scalar res = asImp_().globalResidual(dest);
871 mutableSolution(/*timeIdx=*/0) = tmp;
872 return res;
873 }
874
881 Scalar globalResidual(GlobalEqVector& dest) const
882 {
883 dest = 0;
884
885 std::mutex mutex;
886 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
887#ifdef _OPENMP
888#pragma omp parallel
889#endif
890 {
891 // Attention: the variables below are thread specific and thus cannot be
892 // moved in front of the #pragma!
893 unsigned threadId = ThreadManager::threadId();
894 ElementContext elemCtx(simulator_);
895 ElementIterator elemIt = threadedElemIt.beginParallel();
896 LocalEvalBlockVector residual, storageTerm;
897
898 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
899 const Element& elem = *elemIt;
900 if (elem.partitionType() != Dune::InteriorEntity)
901 continue;
902
903 elemCtx.updateAll(elem);
904 residual.resize(elemCtx.numDof(/*timeIdx=*/0));
905 storageTerm.resize(elemCtx.numPrimaryDof(/*timeIdx=*/0));
906 asImp_().localResidual(threadId).eval(residual, elemCtx);
907
908 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
909 mutex.lock();
910 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
911 unsigned globalI = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
912 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
913 dest[globalI][eqIdx] += Toolbox::value(residual[dofIdx][eqIdx]);
914 }
915 mutex.unlock();
916 }
917 }
918
919 // add up the residuals on the process borders
920 const auto sumHandle =
921 GridCommHandleFactory::template sumHandle<EqVector>(dest, asImp_().dofMapper());
922 gridView_.communicate(*sumHandle,
923 Dune::InteriorBorder_InteriorBorder_Interface,
924 Dune::ForwardCommunication);
925
926 // calculate the square norm of the residual. this is not
927 // entirely correct, since the residual for the finite volumes
928 // which are on the boundary are counted once for every
929 // process. As often in life: shit happens (, we don't care)...
930 Scalar result2 = dest.two_norm2();
931 result2 = asImp_().gridView().comm().sum(result2);
932
933 return std::sqrt(result2);
934 }
935
942 void globalStorage(EqVector& storage, unsigned timeIdx = 0) const
943 {
944 storage = 0;
945
946 std::mutex mutex;
947 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
948#ifdef _OPENMP
949#pragma omp parallel
950#endif
951 {
952 // Attention: the variables below are thread specific and thus cannot be
953 // moved in front of the #pragma!
954 unsigned threadId = ThreadManager::threadId();
955 ElementContext elemCtx(simulator_);
956 ElementIterator elemIt = threadedElemIt.beginParallel();
957 LocalEvalBlockVector elemStorage;
958
959 // in this method, we need to disable the storage cache because we want to
960 // evaluate the storage term for other time indices than the most recent one
961 elemCtx.setEnableStorageCache(false);
962
963 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
964 const Element& elem = *elemIt;
965 if (elem.partitionType() != Dune::InteriorEntity)
966 continue; // ignore ghost and overlap elements
967
968 elemCtx.updateStencil(elem);
969 elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
970
971 size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
972 elemStorage.resize(numPrimaryDof);
973
974 localResidual(threadId).evalStorage(elemStorage, elemCtx, timeIdx);
975
976 mutex.lock();
977 for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx)
978 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
979 storage[eqIdx] += Toolbox::value(elemStorage[dofIdx][eqIdx]);
980 mutex.unlock();
981 }
982 }
983
984 storage = gridView_.comm().sum(storage);
985 }
986
994 void checkConservativeness([[maybe_unused]] Scalar tolerance = -1,
995 [[maybe_unused]] bool verbose=false) const
996 {
997#ifndef NDEBUG
998 Scalar totalBoundaryArea(0.0);
999 Scalar totalVolume(0.0);
1000 EvalEqVector totalRate(0.0);
1001
1002 // take the newton tolerance times the total volume of the grid if we're not
1003 // given an explicit tolerance...
1004 if (tolerance <= 0) {
1005 tolerance =
1006 simulator_.model().newtonMethod().tolerance()
1007 * simulator_.model().gridTotalVolume()
1008 * 1000;
1009 }
1010
1011 // we assume the implicit Euler time discretization for now...
1012 assert(historySize == 2);
1013
1014 EqVector storageBeginTimeStep(0.0);
1015 globalStorage(storageBeginTimeStep, /*timeIdx=*/1);
1016
1017 EqVector storageEndTimeStep(0.0);
1018 globalStorage(storageEndTimeStep, /*timeIdx=*/0);
1019
1020 // calculate the rate at the boundary and the source rate
1021 ElementContext elemCtx(simulator_);
1022 elemCtx.setEnableStorageCache(false);
1023 auto eIt = simulator_.gridView().template begin</*codim=*/0>();
1024 const auto& elemEndIt = simulator_.gridView().template end</*codim=*/0>();
1025 for (; eIt != elemEndIt; ++eIt) {
1026 if (eIt->partitionType() != Dune::InteriorEntity)
1027 continue; // ignore ghost and overlap elements
1028
1029 elemCtx.updateAll(*eIt);
1030
1031 // handle the boundary terms
1032 if (elemCtx.onBoundary()) {
1033 BoundaryContext boundaryCtx(elemCtx);
1034
1035 for (unsigned faceIdx = 0; faceIdx < boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); ++faceIdx) {
1036 BoundaryRateVector values;
1037 simulator_.problem().boundary(values,
1039 faceIdx,
1040 /*timeIdx=*/0);
1041 Valgrind::CheckDefined(values);
1042
1043 unsigned dofIdx = boundaryCtx.interiorScvIndex(faceIdx, /*timeIdx=*/0);
1044 const auto& insideIntQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1045
1046 Scalar bfArea =
1047 boundaryCtx.boundarySegmentArea(faceIdx, /*timeIdx=*/0)
1048 * insideIntQuants.extrusionFactor();
1049
1050 for (unsigned i = 0; i < values.size(); ++i)
1051 values[i] *= bfArea;
1052
1054 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
1055 totalRate[eqIdx] += values[eqIdx];
1056 }
1057 }
1058
1059 // deal with the source terms
1060 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
1061 RateVector values;
1062 simulator_.problem().source(values,
1063 elemCtx,
1064 dofIdx,
1065 /*timeIdx=*/0);
1066 Valgrind::CheckDefined(values);
1067
1068 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1069 Scalar dofVolume =
1070 elemCtx.dofVolume(dofIdx, /*timeIdx=*/0)
1071 * intQuants.extrusionFactor();
1072 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
1073 totalRate[eqIdx] += -dofVolume*Toolbox::value(values[eqIdx]);
1074 totalVolume += dofVolume;
1075 }
1076 }
1077
1078 // summarize everything over all processes
1079 const auto& comm = simulator_.gridView().comm();
1080 totalRate = comm.sum(totalRate);
1082 totalVolume = comm.sum(totalVolume);
1083
1084 if (comm.rank() == 0) {
1087 storageRate /= simulator_.timeStepSize();
1088 if (verbose) {
1089 std::cout << "storage at beginning of time step: " << storageBeginTimeStep << "\n";
1090 std::cout << "storage at end of time step: " << storageEndTimeStep << "\n";
1091 std::cout << "rate based on storage terms: " << storageRate << "\n";
1092 std::cout << "rate based on source and boundary terms: " << totalRate << "\n";
1093 std::cout << "difference in rates: ";
1094 for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx)
1095 std::cout << (storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) << " ";
1096 std::cout << "\n";
1097 }
1098 for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
1099 Scalar eps =
1100 (std::abs(storageRate[eqIdx]) + Toolbox::value(totalRate[eqIdx]))*tolerance;
1101 eps = std::max(tolerance, eps);
1102 assert(std::abs(storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) <= eps);
1103 }
1104 }
1105#endif // NDEBUG
1106 }
1107
1113 Scalar dofTotalVolume(unsigned globalIdx) const
1114 { return dofTotalVolume_[globalIdx]; }
1115
1121 bool isLocalDof(unsigned globalIdx) const
1122 { return isLocalDof_[globalIdx]; }
1123
1128 Scalar gridTotalVolume() const
1129 { return gridTotalVolume_; }
1130
1136 const SolutionVector& solution(unsigned timeIdx) const
1137 { return solution_[timeIdx]->blockVector(); }
1138
1142 SolutionVector& solution(unsigned timeIdx)
1143 { return solution_[timeIdx]->blockVector(); }
1144
1145 protected:
1149 SolutionVector& mutableSolution(unsigned timeIdx) const
1150 { return solution_[timeIdx]->blockVector(); }
1151
1152 public:
1157 const Linearizer& linearizer() const
1158 { return *linearizer_; }
1159
1164 Linearizer& linearizer()
1165 { return *linearizer_; }
1166
1175 const LocalLinearizer& localLinearizer(unsigned openMpThreadId) const
1176 { return localLinearizer_[openMpThreadId]; }
1180 LocalLinearizer& localLinearizer(unsigned openMpThreadId)
1181 { return localLinearizer_[openMpThreadId]; }
1182
1186 const LocalResidual& localResidual(unsigned openMpThreadId) const
1187 { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1191 LocalResidual& localResidual(unsigned openMpThreadId)
1192 { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1193
1201 Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
1202 {
1203 Scalar absPv = std::abs(asImp_().solution(/*timeIdx=*/1)[globalDofIdx][pvIdx]);
1204 return 1.0/std::max(absPv, 1.0);
1205 }
1206
1213 Scalar eqWeight(unsigned, unsigned) const
1214 { return 1.0; }
1215
1226 const PrimaryVariables& pv1,
1227 const PrimaryVariables& pv2) const
1228 {
1229 Scalar result = 0.0;
1230 for (unsigned j = 0; j < numEq; ++j) {
1231 Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
1232 Scalar eqErr = std::abs((pv1[j] - pv2[j])*weight);
1233 //Scalar eqErr = std::abs(pv1[j] - pv2[j]);
1234 //eqErr *= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
1235
1236 result = std::max(result, eqErr);
1237 }
1238 return result;
1239 }
1240
1246 bool update()
1247 {
1248 TimerGuard prePostProcessGuard(prePostProcessTimer_);
1249
1250#ifndef NDEBUG
1251 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1252 // Make sure that the primary variables are defined. Note that because of padding
1253 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1254 // for definedness...
1255 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1256 asImp_().solution(timeIdx)[i].checkDefined();
1257 }
1258 }
1259#endif // NDEBUG
1260
1261 // make sure all timers are prestine
1262 prePostProcessTimer_.halt();
1263 linearizeTimer_.halt();
1264 solveTimer_.halt();
1265 updateTimer_.halt();
1266
1267 prePostProcessTimer_.start();
1268 asImp_().updateBegin();
1269 prePostProcessTimer_.stop();
1270
1271 bool converged = false;
1272
1273 try {
1274 converged = newtonMethod_.apply();
1275 }
1276 catch(...) {
1277 prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1278 linearizeTimer_ += newtonMethod_.linearizeTimer();
1279 solveTimer_ += newtonMethod_.solveTimer();
1280 updateTimer_ += newtonMethod_.updateTimer();
1281
1282 throw;
1283 }
1284
1285#ifndef NDEBUG
1286 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1287 // Make sure that the primary variables are defined. Note that because of padding
1288 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1289 // for definedness...
1290 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1291 asImp_().solution(timeIdx)[i].checkDefined();
1292 }
1293 }
1294#endif // NDEBUG
1295
1296 prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1297 linearizeTimer_ += newtonMethod_.linearizeTimer();
1298 solveTimer_ += newtonMethod_.solveTimer();
1299 updateTimer_ += newtonMethod_.updateTimer();
1300
1301 prePostProcessTimer_.start();
1302 if (converged)
1303 asImp_().updateSuccessful();
1304 else
1305 asImp_().updateFailed();
1306 prePostProcessTimer_.stop();
1307
1308#ifndef NDEBUG
1309 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1310 // Make sure that the primary variables are defined. Note that because of padding
1311 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1312 // for definedness...
1313 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1314 asImp_().solution(timeIdx)[i].checkDefined();
1315 }
1316 }
1317#endif // NDEBUG
1318
1319 return converged;
1320 }
1321
1330 { }
1331
1338 { }
1339
1345 { }
1346
1351 {
1352 throw std::invalid_argument("Grid adaptation need to be implemented for "
1353 "specific settings of grid and function spaces");
1354 }
1355
1362 {
1363 // Reset the current solution to the one of the
1364 // previous time step so that we can start the next
1365 // update at a physically meaningful solution.
1366 solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
1367 invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1368
1369#ifndef NDEBUG
1370 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1371 // Make sure that the primary variables are defined. Note that because of padding
1372 // bytes, we can't just simply ask valgrind to check the whole solution vectors
1373 // for definedness...
1374 for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i)
1375 asImp_().solution(timeIdx)[i].checkDefined();
1376 }
1377#endif // NDEBUG
1378 }
1379
1388 {
1389 // at this point we can adapt the grid
1390 if (this->enableGridAdaptation_) {
1391 asImp_().adaptGrid();
1392 }
1393
1394 // make the current solution the previous one.
1395 solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
1396
1397 // shift the intensive quantities cache by one position in the
1398 // history
1399 asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
1400 }
1401
1409 template <class Restarter>
1411 {
1412 throw std::runtime_error("Not implemented: The discretization chosen for this problem "
1413 "does not support restart files. (serialize() method unimplemented)");
1414 }
1415
1423 template <class Restarter>
1425 {
1426 throw std::runtime_error("Not implemented: The discretization chosen for this problem "
1427 "does not support restart files. (deserialize() method unimplemented)");
1428 }
1429
1438 template <class DofEntity>
1439 void serializeEntity(std::ostream& outstream,
1440 const DofEntity& dof)
1441 {
1442 unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1443
1444 // write phase state
1445 if (!outstream.good()) {
1446 throw std::runtime_error("Could not serialize degree of freedom "
1447 +std::to_string(dofIdx));
1448 }
1449
1450 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1451 outstream << solution(/*timeIdx=*/0)[dofIdx][eqIdx] << " ";
1452 }
1453 }
1454
1463 template <class DofEntity>
1464 void deserializeEntity(std::istream& instream,
1465 const DofEntity& dof)
1466 {
1467 unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1468
1469 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1470 if (!instream.good())
1471 throw std::runtime_error("Could not deserialize degree of freedom "
1472 +std::to_string(dofIdx));
1473 instream >> solution(/*timeIdx=*/0)[dofIdx][eqIdx];
1474 }
1475 }
1476
1480 size_t numGridDof() const
1481 { throw std::logic_error("The discretization class must implement the numGridDof() method!"); }
1482
1486 size_t numAuxiliaryDof() const
1487 {
1488 size_t result = 0;
1489 auto auxModIt = auxEqModules_.begin();
1490 const auto& auxModEndIt = auxEqModules_.end();
1491 for (; auxModIt != auxModEndIt; ++auxModIt)
1492 result += (*auxModIt)->numDofs();
1493
1494 return result;
1495 }
1496
1500 size_t numTotalDof() const
1501 { return asImp_().numGridDof() + numAuxiliaryDof(); }
1502
1507 const DofMapper& dofMapper() const
1508 { throw std::logic_error("The discretization class must implement the dofMapper() method!"); }
1509
1513 const VertexMapper& vertexMapper() const
1514 { return vertexMapper_; }
1515
1519 const ElementMapper& elementMapper() const
1520 { return elementMapper_; }
1521
1527 {
1528 delete linearizer_;
1529 linearizer_ = new Linearizer;
1530 linearizer_->init(simulator_);
1531 }
1532
1536 static std::string discretizationName()
1537 { return ""; }
1538
1544 std::string primaryVarName(unsigned pvIdx) const
1545 {
1546 std::ostringstream oss;
1547 oss << "primary variable_" << pvIdx;
1548 return oss.str();
1549 }
1550
1556 std::string eqName(unsigned eqIdx) const
1557 {
1558 std::ostringstream oss;
1559 oss << "equation_" << eqIdx;
1560 return oss.str();
1561 }
1562
1569 void updatePVWeights(const ElementContext&) const
1570 { }
1571
1576 { outputModules_.push_back(newModule); }
1577
1586 template <class VtkMultiWriter>
1588 const SolutionVector& u,
1589 const GlobalEqVector& deltaU) const
1590 {
1591 using ScalarBuffer = std::vector<double>;
1592
1593 GlobalEqVector globalResid(u.size());
1594 asImp_().globalResidual(globalResid, u);
1595
1596 // create the required scalar fields
1597 size_t numGridDof = asImp_().numGridDof();
1598
1599 // global defect of the two auxiliary equations
1600 ScalarBuffer* def[numEq];
1601 ScalarBuffer* delta[numEq];
1602 ScalarBuffer* priVars[numEq];
1603 ScalarBuffer* priVarWeight[numEq];
1604 ScalarBuffer* relError = writer.allocateManagedScalarBuffer(numGridDof);
1605 ScalarBuffer* normalizedRelError = writer.allocateManagedScalarBuffer(numGridDof);
1606 for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1607 priVars[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1608 priVarWeight[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1609 delta[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1610 def[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1611 }
1612
1613 Scalar minRelErr = 1e30;
1614 Scalar maxRelErr = -1e30;
1615 for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx) {
1616 for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1617 (*priVars[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
1618 (*priVarWeight[pvIdx])[globalIdx] = asImp_().primaryVarWeight(globalIdx, pvIdx);
1619 (*delta[pvIdx])[globalIdx] = - deltaU[globalIdx][pvIdx];
1620 (*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
1621 }
1622
1623 PrimaryVariables uOld(u[globalIdx]);
1624 PrimaryVariables uNew(uOld);
1625 uNew -= deltaU[globalIdx];
1626
1627 Scalar err = asImp_().relativeDofError(globalIdx, uOld, uNew);
1628 (*relError)[globalIdx] = err;
1629 (*normalizedRelError)[globalIdx] = err;
1630 minRelErr = std::min(err, minRelErr);
1631 maxRelErr = std::max(err, maxRelErr);
1632 }
1633
1634 // do the normalization of the relative error
1635 Scalar alpha = std::max(Scalar{1e-20},
1636 std::max(std::abs(maxRelErr),
1637 std::abs(minRelErr)));
1638 for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx)
1639 (*normalizedRelError)[globalIdx] /= alpha;
1640
1641 DiscBaseOutputModule::attachScalarDofData_(writer, *relError, "relative error");
1642 DiscBaseOutputModule::attachScalarDofData_(writer, *normalizedRelError, "normalized relative error");
1643
1644 for (unsigned i = 0; i < numEq; ++i) {
1645 std::ostringstream oss;
1646 oss.str(""); oss << "priVar_" << asImp_().primaryVarName(i);
1647 DiscBaseOutputModule::attachScalarDofData_(writer,
1648 *priVars[i],
1649 oss.str());
1650
1651 oss.str(""); oss << "delta_" << asImp_().primaryVarName(i);
1652 DiscBaseOutputModule::attachScalarDofData_(writer,
1653 *delta[i],
1654 oss.str());
1655
1656 oss.str(""); oss << "weight_" << asImp_().primaryVarName(i);
1657 DiscBaseOutputModule::attachScalarDofData_(writer,
1658 *priVarWeight[i],
1659 oss.str());
1660
1661 oss.str(""); oss << "defect_" << asImp_().eqName(i);
1662 DiscBaseOutputModule::attachScalarDofData_(writer,
1663 *def[i],
1664 oss.str());
1665 }
1666
1667 asImp_().prepareOutputFields();
1668 asImp_().appendOutputFields(writer);
1669 }
1670
1676 {
1677 bool needFullContextUpdate = false;
1678 auto modIt = outputModules_.begin();
1679 const auto& modEndIt = outputModules_.end();
1680 for (; modIt != modEndIt; ++modIt) {
1681 (*modIt)->allocBuffers();
1682 needFullContextUpdate = needFullContextUpdate || (*modIt)->needExtensiveQuantities();
1683 }
1684
1685 // iterate over grid
1686 ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
1687#ifdef _OPENMP
1688#pragma omp parallel
1689#endif
1690 {
1691 ElementContext elemCtx(simulator_);
1692 ElementIterator elemIt = threadedElemIt.beginParallel();
1693 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
1694 const auto& elem = *elemIt;
1695 if (elem.partitionType() != Dune::InteriorEntity)
1696 // ignore non-interior entities
1697 continue;
1698
1700 elemCtx.updateAll(elem);
1701 else {
1702 elemCtx.updatePrimaryStencil(elem);
1703 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1704 }
1705
1706 // we cannot reuse the "modIt" variable here because the code here might
1707 // be threaded and "modIt" is is the same for all threads, i.e., if a
1708 // given thread modifies it, the changes affect all threads.
1709 auto modIt2 = outputModules_.begin();
1710 for (; modIt2 != modEndIt; ++modIt2)
1711 (*modIt2)->processElement(elemCtx);
1712 }
1713 }
1714 }
1715
1721 {
1722 auto modIt = outputModules_.begin();
1723 const auto& modEndIt = outputModules_.end();
1724 for (; modIt != modEndIt; ++modIt)
1725 (*modIt)->commitBuffers(writer);
1726 }
1727
1731 const GridView& gridView() const
1732 { return gridView_; }
1733
1746 {
1747 auxMod->setDofOffset(numTotalDof());
1748 auxEqModules_.push_back(auxMod);
1749
1750 // resize the solutions
1751 if (enableGridAdaptation_
1752 && !std::is_same<DiscreteFunction, BlockVectorWrapper>::value)
1753 {
1754 throw std::invalid_argument("Problems which require auxiliary modules cannot be used in"
1755 " conjunction with dune-fem");
1756 }
1757
1758 size_t numDof = numTotalDof();
1759 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx)
1760 solution(timeIdx).resize(numDof);
1761
1762 auxMod->applyInitial();
1763 }
1764
1771 {
1772 auxEqModules_.clear();
1773 linearizer_->eraseMatrix();
1774 newtonMethod_.eraseMatrix();
1775 }
1776
1780 size_t numAuxiliaryModules() const
1781 { return auxEqModules_.size(); }
1782
1787 { return auxEqModules_[auxEqModIdx]; }
1788
1793 { return auxEqModules_[auxEqModIdx]; }
1794
1799 { return enableIntensiveQuantityCache_ || enableThermodynamicHints_; }
1800
1801 const Timer& prePostProcessTimer() const
1802 { return prePostProcessTimer_; }
1803
1804 const Timer& linearizeTimer() const
1805 { return linearizeTimer_; }
1806
1807 const Timer& solveTimer() const
1808 { return solveTimer_; }
1809
1810 const Timer& updateTimer() const
1811 { return updateTimer_; }
1812
1813 template<class Serializer>
1814 void serializeOp(Serializer& serializer)
1815 {
1817 using Helper = typename BaseDiscretization::template SerializeHelper<Serializer>;
1818 Helper::serializeOp(serializer, solution_);
1819 }
1820
1821 bool operator==(const FvBaseDiscretization& rhs) const
1822 {
1823 return std::equal(this->solution_.begin(), this->solution_.end(),
1824 rhs.solution_.begin(), rhs.solution_.end(),
1825 [](const auto& x, const auto& y)
1826 {
1827 return *x == *y;
1828 });
1829 }
1830
1831protected:
1832 void resizeAndResetIntensiveQuantitiesCache_()
1833 {
1834 // allocate the storage cache
1835 if (enableStorageCache()) {
1836 size_t numDof = asImp_().numGridDof();
1837 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1838 storageCache_[timeIdx].resize(numDof);
1839 }
1840 }
1841
1842 // allocate the intensive quantities cache
1844 size_t numDof = asImp_().numGridDof();
1845 for(unsigned timeIdx=0; timeIdx<historySize; ++timeIdx) {
1846 intensiveQuantityCache_[timeIdx].resize(numDof);
1847 intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
1849 }
1850 }
1851 }
1852 template <class Context>
1853 void supplementInitialSolution_(PrimaryVariables&,
1854 const Context&,
1855 unsigned,
1856 unsigned)
1857 { }
1858
1867 {
1868 // add the output modules available on all model
1869 auto *mod = new VtkPrimaryVarsModule<TypeTag>(simulator_);
1870 this->outputModules_.push_back(mod);
1871 }
1872
1876 LocalResidual& localResidual_()
1877 { return localLinearizer_.localResidual(); }
1878
1882 bool verbose_() const
1883 { return gridView_.comm().rank() == 0; }
1884
1885 Implementation& asImp_()
1886 { return *static_cast<Implementation*>(this); }
1887 const Implementation& asImp_() const
1888 { return *static_cast<const Implementation*>(this); }
1889
1890 // the problem we want to solve. defines the constitutive
1891 // relations, matxerial laws, etc.
1892 Simulator& simulator_;
1893
1894 // the representation of the spatial domain of the problem
1895 GridView gridView_;
1896
1897 // the mappers for element and vertex entities to global indices
1898 ElementMapper elementMapper_;
1899 VertexMapper vertexMapper_;
1900
1901 // a vector with all auxiliary equations to be considered
1902 std::vector<BaseAuxiliaryModule<TypeTag>*> auxEqModules_;
1903
1904 NewtonMethod newtonMethod_;
1905
1906 Timer prePostProcessTimer_;
1907 Timer linearizeTimer_;
1908 Timer solveTimer_;
1909 Timer updateTimer_;
1910
1911 // calculates the local jacobian matrix for a given element
1912 std::vector<LocalLinearizer> localLinearizer_;
1913 // Linearizes the problem at the current time step using the
1914 // local jacobian
1915 Linearizer *linearizer_;
1916
1917 // cur is the current iterative solution, prev the converged
1918 // solution of the previous time step
1919 mutable IntensiveQuantitiesVector intensiveQuantityCache_[historySize];
1920 // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
1921 mutable std::vector<unsigned char> intensiveQuantityCacheUpToDate_[historySize];
1922
1923 mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
1924
1925 std::list<BaseOutputModule<TypeTag>*> outputModules_;
1926
1927 Scalar gridTotalVolume_;
1928 std::vector<Scalar> dofTotalVolume_;
1929 std::vector<bool> isLocalDof_;
1930
1931 mutable GlobalEqVector storageCache_[historySize];
1932
1933 bool enableGridAdaptation_;
1934 bool enableIntensiveQuantityCache_;
1935 bool enableStorageCache_;
1936 bool enableThermodynamicHints_;
1937};
1938
1944template<class TypeTag>
1946{
1950
1951 static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
1952
1953public:
1954 template<class Serializer>
1956 template<class SolutionType>
1957 static void serializeOp(Serializer& serializer,
1959 {
1960 for (auto& sol : solution) {
1961 serializer(*sol);
1962 }
1963 }
1964 };
1965
1966 FvBaseDiscretizationNoAdapt(Simulator& simulator)
1967 : ParentType(simulator)
1968 {
1969 if (this->enableGridAdaptation_) {
1970 throw std::invalid_argument("Grid adaptation need to use"
1971 " BaseDiscretization = FvBaseDiscretizationFemAdapt"
1972 " which currently requires the presence of the"
1973 " dune-fem module");
1974 }
1975 size_t numDof = this->asImp_().numGridDof();
1976 for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1977 this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", numDof);
1978 }
1979 }
1980};
1981
1982} // namespace Opm
1983
1984#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:1946
Definition fvbasediscretization.hh:343
The base class for the finite volume discretization schemes.
Definition fvbasediscretization.hh:293
LocalLinearizer & localLinearizer(unsigned openMpThreadId)
Definition fvbasediscretization.hh:1180
size_t numAuxiliaryDof() const
Returns the number of degrees of freedom (DOFs) of the auxiliary equations.
Definition fvbasediscretization.hh:1486
void prepareOutputFields() const
Prepare the quantities relevant for the current solution to be appended to the output writers.
Definition fvbasediscretization.hh:1675
void shiftIntensiveQuantityCache(unsigned numSlots=1)
Move the intensive quantities for a given time index to the back.
Definition fvbasediscretization.hh:780
void adaptGrid()
Called by the update() method when the grid should be refined.
Definition fvbasediscretization.hh:1350
void addAuxiliaryModule(BaseAuxiliaryModule< TypeTag > *auxMod)
Add a module for an auxiliary equation.
Definition fvbasediscretization.hh:1745
void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx, unsigned timeIdx, bool newValue) const
Invalidate the cache for a given intensive quantities object.
Definition fvbasediscretization.hh:698
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:602
void updateSuccessful()
Called by the update() method if it was successful.
Definition fvbasediscretization.hh:1344
static std::string discretizationName()
Returns a string of discretization's human-readable name.
Definition fvbasediscretization.hh:1536
BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx)
Returns a given module for auxiliary equations.
Definition fvbasediscretization.hh:1786
bool isLocalDof(unsigned globalIdx) const
Returns if the overlap of the volume ofa degree of freedom is non-zero.
Definition fvbasediscretization.hh:1121
LocalResidual & localResidual_()
Reference to the local residal object.
Definition fvbasediscretization.hh:1876
void registerOutputModules_()
Register all output modules which make sense for the model.
Definition fvbasediscretization.hh:1866
bool enableGridAdaptation() const
Returns whether the grid ought to be adapted to the solution during the simulation.
Definition fvbasediscretization.hh:541
const NewtonMethod & newtonMethod() const
Returns the newton method object.
Definition fvbasediscretization.hh:616
SolutionVector & mutableSolution(unsigned timeIdx) const
Definition fvbasediscretization.hh:1149
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition fvbasediscretization.hh:1387
NewtonMethod & newtonMethod()
Returns the newton method object.
Definition fvbasediscretization.hh:610
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition fvbasediscretization.hh:1513
const EqVector & cachedStorage(unsigned globalIdx, unsigned timeIdx) const
Retrieve an entry of the cache for the storage term.
Definition fvbasediscretization.hh:835
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:1569
void updateFailed()
Called by the update() method if it was unsuccessful.
Definition fvbasediscretization.hh:1361
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition fvbasediscretization.hh:1337
Scalar globalResidual(GlobalEqVector &dest) const
Compute the global residual for the current solution vector.
Definition fvbasediscretization.hh:881
LocalResidual & localResidual(unsigned openMpThreadId)
Definition fvbasediscretization.hh:1191
void serializeEntity(std::ostream &outstream, const DofEntity &dof)
Write the current solution for a degree of freedom to a restart file.
Definition fvbasediscretization.hh:1439
void setEnableStorageCache(bool enableStorageCache)
Set the value of enable storage cache.
Definition fvbasediscretization.hh:822
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition fvbasediscretization.hh:1556
Scalar eqWeight(unsigned, unsigned) const
Returns the relative weight of an equation.
Definition fvbasediscretization.hh:1213
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition fvbasediscretization.hh:1201
void deserialize(Restarter &)
Deserializes the state of the model.
Definition fvbasediscretization.hh:1424
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:994
Scalar gridTotalVolume() const
Returns the volume of the whole grid which represents the spatial domain.
Definition fvbasediscretization.hh:1128
const IntensiveQuantities * thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
Return the thermodynamic hint for a entity on the grid at given time.
Definition fvbasediscretization.hh:634
const Linearizer & linearizer() const
Returns the operator linearizer for the global jacobian of the problem.
Definition fvbasediscretization.hh:1157
void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector &value) const
Set an entry of the cache for the storage term.
Definition fvbasediscretization.hh:852
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:1587
Scalar relativeDofError(unsigned vertexIdx, const PrimaryVariables &pv1, const PrimaryVariables &pv2) const
Returns the relative error between two vectors of primary variables.
Definition fvbasediscretization.hh:1225
Scalar globalResidual(GlobalEqVector &dest, const SolutionVector &u) const
Compute the global residual for an arbitrary solution vector.
Definition fvbasediscretization.hh:865
SolutionVector & solution(unsigned timeIdx)
Definition fvbasediscretization.hh:1142
void addOutputModule(BaseOutputModule< TypeTag > *newModule)
Add an module for writing visualization output after a timestep.
Definition fvbasediscretization.hh:1575
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:1464
const LocalLinearizer & localLinearizer(unsigned openMpThreadId) const
Returns the local jacobian which calculates the local stiffness matrix for an arbitrary element.
Definition fvbasediscretization.hh:1175
bool update()
Try to progress the model to the next timestep.
Definition fvbasediscretization.hh:1246
void clearAuxiliaryModules()
Causes the list of auxiliary equations to be cleared.
Definition fvbasediscretization.hh:1770
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition fvbasediscretization.hh:1519
size_t numGridDof() const
Returns the number of degrees of freedom (DOFs) for the computational grid.
Definition fvbasediscretization.hh:1480
void syncOverlap()
Syncronize the values of the primary variables on the degrees of freedom that overlap with the neighb...
Definition fvbasediscretization.hh:1329
void appendOutputFields(BaseOutputWriter &writer) const
Append the quantities relevant for the current solution to an output writer.
Definition fvbasediscretization.hh:1720
void applyInitialSolution()
Applies the initial solution for all degrees of freedom to which the model applies.
Definition fvbasediscretization.hh:548
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:680
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition fvbasediscretization.hh:1544
void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
Invalidate the whole intensive quantity cache for time index.
Definition fvbasediscretization.hh:713
Linearizer & linearizer()
Returns the object which linearizes the global system of equations at the current solution.
Definition fvbasediscretization.hh:1164
const BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx) const
Returns a given module for auxiliary equations.
Definition fvbasediscretization.hh:1792
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:942
Scalar dofTotalVolume(unsigned globalIdx) const
Returns the volume of a given control volume.
Definition fvbasediscretization.hh:1113
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:1136
bool verbose_() const
Returns whether messages should be printed.
Definition fvbasediscretization.hh:1882
const LocalResidual & localResidual(unsigned openMpThreadId) const
Returns the object to calculate the local residual function.
Definition fvbasediscretization.hh:1186
const DofMapper & dofMapper() const
Mapper to convert the Dune entities of the discretization's degrees of freedoms are to indices.
Definition fvbasediscretization.hh:1507
bool storeIntensiveQuantities() const
Returns true if the cache for intensive quantities is enabled.
Definition fvbasediscretization.hh:1798
void serialize(Restarter &)
Serializes the current state of the model.
Definition fvbasediscretization.hh:1410
size_t numTotalDof() const
Returns the total number of degrees of freedom (i.e., grid plux auxiliary DOFs)
Definition fvbasediscretization.hh:1500
void resetLinearizer()
Resets the Jacobian matrix linearizer, so that the boundary types can be altered.
Definition fvbasediscretization.hh:1526
bool enableStorageCache() const
Returns true iff the storage term is cached.
Definition fvbasediscretization.hh:813
size_t numAuxiliaryModules() const
Returns the number of modules for auxiliary equations.
Definition fvbasediscretization.hh:1780
const GridView & gridView() const
Reference to the grid view of the spatial domain.
Definition fvbasediscretization.hh:1731
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:654
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:92
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:1955
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.