Loading...
Searching...
No Matches
fvbasediscretization.hh
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_FV_BASE_DISCRETIZATION_HH
#define EWOMS_FV_BASE_DISCRETIZATION_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/version.hh>
#include <dune/istl/bvector.hh>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/models/discretization/common/fvbaseadlocallinearizer.hh>
#include <opm/models/discretization/common/fvbaseboundarycontext.hh>
#include <opm/models/discretization/common/fvbaseconstraints.hh>
#include <opm/models/discretization/common/fvbaseconstraintscontext.hh>
#include <opm/models/discretization/common/fvbaseelementcontext.hh>
#include <opm/models/discretization/common/fvbaseextensivequantities.hh>
#include <opm/models/discretization/common/fvbasefdlocallinearizer.hh>
#include <opm/models/discretization/common/fvbasegradientcalculator.hh>
#include <opm/models/discretization/common/fvbaseintensivequantities.hh>
#include <opm/models/discretization/common/fvbaselinearizer.hh>
#include <opm/models/discretization/common/fvbaselocalresidual.hh>
#include <opm/models/discretization/common/fvbasenewtonmethod.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
#include <opm/models/io/vtkprimaryvarsmodule.hpp>
#include <opm/models/parallel/gridcommhandles.hh>
#include <opm/models/parallel/threadmanager.hpp>
#include <opm/models/utils/alignedallocator.hh>
#include <opm/models/utils/simulator.hh>
#include <opm/models/utils/timer.hpp>
#include <opm/models/utils/timerguard.hh>
#include <opm/simulators/linalg/linalgparameters.hh>
#include <opm/simulators/linalg/nullborderlistmanager.hh>
#include <algorithm>
#include <array>
#include <cstddef>
#include <list>
#include <memory>
#include <stdexcept>
#include <sstream>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
namespace Opm {
template<class TypeTag>
class FvBaseDiscretizationNoAdapt;
template<class TypeTag>
class FvBaseDiscretization;
} // namespace Opm
namespace Opm::Properties {
template<class TypeTag>
struct Simulator<TypeTag, TTag::FvBaseDiscretization>
{ using type = ::Opm::Simulator<TypeTag>; };
template<class TypeTag>
struct VertexMapper<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
template<class TypeTag>
struct ElementMapper<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
template<class TypeTag>
struct BorderListCreator<TypeTag, TTag::FvBaseDiscretization>
{
using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
public:
using type = Linear::NullBorderListCreator<GridView, DofMapper>;
};
template<class TypeTag>
struct DiscLocalResidual<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseLocalResidual<TypeTag>; };
template<class TypeTag>
struct DiscIntensiveQuantities<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseIntensiveQuantities<TypeTag>; };
template<class TypeTag>
struct DiscExtensiveQuantities<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseExtensiveQuantities<TypeTag>; };
template<class TypeTag>
struct GradientCalculator<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseGradientCalculator<TypeTag>; };
template<class TypeTag>
struct EqVector<TypeTag, TTag::FvBaseDiscretization>
{
using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
getPropValue<TypeTag, Properties::NumEq>()>;
};
template<class TypeTag>
struct RateVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = GetPropType<TypeTag, Properties::EqVector>; };
template<class TypeTag>
struct BoundaryRateVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = GetPropType<TypeTag, Properties::RateVector>; };
template<class TypeTag>
struct Constraints<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseConstraints<TypeTag>; };
template<class TypeTag>
struct ElementEqVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
template<class TypeTag>
struct GlobalEqVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
template<class TypeTag>
struct PrimaryVariables<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBasePrimaryVariables<TypeTag>; };
template<class TypeTag>
struct SolutionVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; };
template<class TypeTag>
struct IntensiveQuantities<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseIntensiveQuantities<TypeTag>; };
template<class TypeTag>
struct ElementContext<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseElementContext<TypeTag>; };
template<class TypeTag>
struct BoundaryContext<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseBoundaryContext<TypeTag>; };
template<class TypeTag>
struct ConstraintsContext<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseConstraintsContext<TypeTag>; };
template<class TypeTag>
struct ThreadManager<TypeTag, TTag::FvBaseDiscretization>
{ using type = ::Opm::ThreadManager; };
template<class TypeTag>
struct UseLinearizationLock<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = true; };
template<class TypeTag>
struct Linearizer<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseLinearizer<TypeTag>; };
template<class TypeTag>
struct VtkOutputFormat<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr auto value = Dune::VTK::ascii; };
// disable constraints by default
template<class TypeTag>
struct EnableConstraints<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = false; };
template<class TypeTag>
struct TimeDiscHistorySize<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr int value = 2; };
template<class TypeTag>
struct ExtensiveStorageTerm<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = false; };
// use volumetric residuals is default
template<class TypeTag>
struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = true; };
template<class TypeTag>
struct EnableExperiments<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = true; };
template <class TypeTag, class MyTypeTag>
struct BaseDiscretizationType { using type = UndefinedProperty; };
#if !HAVE_DUNE_FEM
template<class TypeTag>
struct BaseDiscretizationType<TypeTag,TTag::FvBaseDiscretization>
{ using type = FvBaseDiscretizationNoAdapt<TypeTag>; };
template<class TypeTag>
struct DiscreteFunction<TypeTag, TTag::FvBaseDiscretization>
{
using BaseDiscretization = FvBaseDiscretization<TypeTag>;
using type = typename BaseDiscretization::BlockVectorWrapper;
};
#endif
} // namespace Opm::Properties
namespace Opm {
template<class TypeTag>
class FvBaseDiscretization
{
using Implementation = GetPropType<TypeTag, Properties::Model>;
using Discretization = GetPropType<TypeTag, Properties::Discretization>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Grid = GetPropType<TypeTag, Properties::Grid>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using VertexMapper = GetPropType<TypeTag, Properties::VertexMapper>;
using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Linearizer = GetPropType<TypeTag, Properties::Linearizer>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using BoundaryContext = GetPropType<TypeTag, Properties::BoundaryContext>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using GradientCalculator = GetPropType<TypeTag, Properties::GradientCalculator>;
using Stencil = GetPropType<TypeTag, Properties::Stencil>;
using DiscBaseOutputModule = GetPropType<TypeTag, Properties::DiscBaseOutputModule>;
using GridCommHandleFactory = GetPropType<TypeTag, Properties::GridCommHandleFactory>;
using NewtonMethod = GetPropType<TypeTag, Properties::NewtonMethod>;
using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
using LocalLinearizer = GetPropType<TypeTag, Properties::LocalLinearizer>;
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
enum {
numEq = getPropValue<TypeTag, Properties::NumEq>(),
historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
};
using IntensiveQuantitiesVector = std::vector<IntensiveQuantities,
aligned_allocator<IntensiveQuantities,
alignof(IntensiveQuantities)>>;
using Element = typename GridView::template Codim<0>::Entity;
using ElementIterator = typename GridView::template Codim<0>::Iterator;
using Toolbox = MathToolbox<Evaluation>;
using VectorBlock = Dune::FieldVector<Evaluation, numEq>;
using EvalEqVector = Dune::FieldVector<Evaluation, numEq>;
using LocalEvalBlockVector = typename LocalResidual::LocalEvalBlockVector;
public:
class BlockVectorWrapper
{
protected:
SolutionVector blockVector_;
public:
BlockVectorWrapper(const std::string&, const std::size_t size)
: blockVector_(size)
{}
BlockVectorWrapper() = default;
static BlockVectorWrapper serializationTestObject()
{
BlockVectorWrapper result("dummy", 3);
result.blockVector_[0] = 1.0;
result.blockVector_[1] = 2.0;
result.blockVector_[2] = 3.0;
return result;
}
SolutionVector& blockVector()
{ return blockVector_; }
const SolutionVector& blockVector() const
{ return blockVector_; }
bool operator==(const BlockVectorWrapper& wrapper) const
{
return std::equal(this->blockVector_.begin(), this->blockVector_.end(),
wrapper.blockVector_.begin(), wrapper.blockVector_.end());
}
template<class Serializer>
void serializeOp(Serializer& serializer)
{
serializer(blockVector_);
}
};
private:
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
using DiscreteFunction = GetPropType<TypeTag, Properties::DiscreteFunction>;
public:
explicit FvBaseDiscretization(Simulator& simulator)
: simulator_(simulator)
, gridView_(simulator.gridView())
, elementMapper_(gridView_, Dune::mcmgElementLayout())
, vertexMapper_(gridView_, Dune::mcmgVertexLayout())
, newtonMethod_(simulator)
, localLinearizer_(ThreadManager::maxThreads())
, linearizer_(std::make_unique<Linearizer>())
, enableGridAdaptation_(Parameters::Get<Parameters::EnableGridAdaptation>() )
, enableIntensiveQuantityCache_(Parameters::Get<Parameters::EnableIntensiveQuantityCache>())
, enableStorageCache_(Parameters::Get<Parameters::EnableStorageCache>())
, enableThermodynamicHints_(Parameters::Get<Parameters::EnableThermodynamicHints>())
, cachedIntensiveQuantityHistorySize_(-1)
{
const bool isEcfv = std::is_same_v<Discretization, EcfvDiscretization<TypeTag>>;
if (enableGridAdaptation_ && !isEcfv) {
throw std::invalid_argument("Grid adaptation currently only works for the "
"element-centered finite volume discretization (is: " +
Dune::className<Discretization>() + ")");
}
PrimaryVariables::init();
// Setting up the intensive quantities cache and storage cache is done in finishInit()
// and applyInitialSolution() to ensure that the history size is correct.
asImp_().registerOutputModules_();
}
// copying a discretization object is not a good idea
FvBaseDiscretization(const FvBaseDiscretization&) = delete;
static void registerParameters()
{
Linearizer::registerParameters();
LocalLinearizer::registerParameters();
LocalResidual::registerParameters();
GradientCalculator::registerParameters();
IntensiveQuantities::registerParameters();
ExtensiveQuantities::registerParameters();
NewtonMethod::registerParameters();
Linearizer::registerParameters();
PrimaryVariables::registerParameters();
// register runtime parameters of the output modules
VtkPrimaryVarsModule<TypeTag>::registerParameters();
Parameters::Register<Parameters::EnableGridAdaptation>
("Enable adaptive grid refinement/coarsening");
Parameters::Register<Parameters::EnableVtkOutput>
("Global switch for turning on writing VTK files");
Parameters::Register<Parameters::EnableThermodynamicHints>
("Enable thermodynamic hints");
Parameters::Register<Parameters::EnableIntensiveQuantityCache>
("Turn on caching of intensive quantities");
Parameters::Register<Parameters::EnableStorageCache>
("Store previous storage terms and avoid re-calculating them.");
Parameters::Register<Parameters::OutputDir>
("The directory to which result files are written");
}
void finishInit()
{
// initialize the volume of the finite volumes to zero
const std::size_t numDof = asImp_().numGridDof();
dofTotalVolume_.resize(numDof);
std::fill(dofTotalVolume_.begin(), dofTotalVolume_.end(), 0.0);
ElementContext elemCtx(simulator_);
gridTotalVolume_ = 0.0;
// iterate through the grid and evaluate the initial condition
for (const auto& elem : elements(gridView_)) {
// ignore everything which is not in the interior if the
// current process' piece of the grid
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
// deal with the current element
elemCtx.updateStencil(elem);
const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
// loop over all element vertices, i.e. sub control volumes
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++) {
// map the local degree of freedom index to the global one
const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
const Scalar dofVolume = stencil.subControlVolume(dofIdx).volume();
dofTotalVolume_[globalIdx] += dofVolume;
gridTotalVolume_ += dofVolume;
}
}
// determine which DOFs should be considered to lie fully in the interior of the
// local process grid partition: those which do not have a non-zero volume
// before taking the peer processes into account...
isLocalDof_.resize(numDof);
for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0);
}
// add the volumes of the DOFs on the process boundaries
const auto sumHandle =
GridCommHandleFactory::template sumHandle<Scalar>(dofTotalVolume_,
asImp_().dofMapper());
gridView_.communicate(*sumHandle,
Dune::InteriorBorder_All_Interface,
Dune::ForwardCommunication);
// sum up the volumes of the grid partitions
gridTotalVolume_ = gridView_.comm().sum(gridTotalVolume_);
linearizer_->init(simulator_);
for (unsigned threadId = 0; threadId < ThreadManager::maxThreads(); ++threadId) {
localLinearizer_[threadId].init(simulator_);
}
resizeAndResetIntensiveQuantitiesCache_();
newtonMethod_.finishInit();
}
bool enableGridAdaptation() const
{ return enableGridAdaptation_; }
void applyInitialSolution()
{
// first set the whole domain to zero
SolutionVector& uCur = asImp_().solution(/*timeIdx=*/0);
uCur = Scalar(0.0);
ElementContext elemCtx(simulator_);
// iterate through the grid and evaluate the initial condition
for (const auto& elem : elements(gridView_)) {
// ignore everything which is not in the interior if the
// current process' piece of the grid
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
// deal with the current element
elemCtx.updateStencil(elem);
// loop over all element vertices, i.e. sub control volumes
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
// map the local degree of freedom index to the global one
const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
// let the problem do the dirty work of nailing down
// the initial solution.
simulator_.problem().initial(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
asImp_().supplementInitialSolution_(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
uCur[globalIdx].checkDefined();
}
}
// synchronize the ghost DOFs (if necessary)
asImp_().syncOverlap();
// also set the solutions of the "previous" time steps to the initial solution.
for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx) {
solution(timeIdx) = solution(/*timeIdx=*/0);
}
// Initialize intensive quantities cache now that all problem-specific parameters are available.
// This ensures intensiveQuantityHistorySize is correct based on recycleFirstIterationStorage().
// TODO: Where this is done should perhaps be changed once finishInit() is refactored.
resizeAndResetIntensiveQuantitiesCache_();
simulator_.problem().initialSolutionApplied();
#ifndef NDEBUG
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
const auto& sol = solution(timeIdx);
for (unsigned dofIdx = 0; dofIdx < sol.size(); ++dofIdx) {
sol[dofIdx].checkDefined();
}
}
#endif // NDEBUG
}
void prefetch(const Element&) const
{
// do nothing by default
}
NewtonMethod& newtonMethod()
{ return newtonMethod_; }
const NewtonMethod& newtonMethod() const
{ return newtonMethod_; }
const IntensiveQuantities* thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
{
if (!enableThermodynamicHints_) {
return 0;
}
// the intensive quantities cache doubles as thermodynamic hint
return cachedIntensiveQuantities(globalIdx, timeIdx);
}
const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
{
if (!enableIntensiveQuantityCache_ ||
timeIdx >= cachedIntensiveQuantityHistorySize_ ||
!intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
return nullptr;
}
// With the storage cache enabled, usually only the
// intensive quantities for the most recent time step are
// cached. However, this may be false for some Problem
// variants, so we should check if the cache exists for
// the timeIdx in question.
if (timeIdx > 0 && enableStorageCache_ && intensiveQuantityCache_[timeIdx].empty()) {
return nullptr;
}
return &intensiveQuantityCache_[timeIdx][globalIdx];
}
void updateCachedIntensiveQuantities(const IntensiveQuantities& intQuants,
unsigned globalIdx,
unsigned timeIdx) const
{
if (!storeIntensiveQuantities() || timeIdx >= cachedIntensiveQuantityHistorySize_) {
return;
}
intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
}
void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx,
unsigned timeIdx,
bool newValue) const
{
if (!storeIntensiveQuantities()) {
return;
}
intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue ? 1 : 0;
}
unsigned cachedIntensiveQuantityHistorySize() const
{ return cachedIntensiveQuantityHistorySize_; }
void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
{
if (timeIdx >= cachedIntensiveQuantityHistorySize_) {
return;
}
if (storeIntensiveQuantities()) {
std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
intensiveQuantityCacheUpToDate_[timeIdx].end(),
/*value=*/0);
}
}
void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
{
invalidateIntensiveQuantitiesCache(timeIdx);
// loop over all elements...
ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
#ifdef _OPENMP
#pragma omp parallel
#endif
{
ElementContext elemCtx(simulator_);
ElementIterator elemIt = threadedElemIt.beginParallel();
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
}
}
}
template <class GridViewType>
void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
{
// loop over all elements...
ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(gridView);
#ifdef _OPENMP
#pragma omp parallel
#endif
{
ElementContext elemCtx(simulator_);
auto elemIt = threadedElemIt.beginParallel();
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
if (elemIt->partitionType() != Dune::InteriorEntity) {
continue;
}
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
// Mark cache for this element as invalid.
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false);
}
// Update for this element.
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
}
}
}
void shiftIntensiveQuantityCache(unsigned numSlots = 1)
{
if (!storeIntensiveQuantities() || numSlots <= 0) {
return;
}
if (enableStorageCache() && simulator_.problem().recycleFirstIterationStorage()) {
// If the storage term is cached, the intensive quantities of the previous
// time steps do not need to be accessed, and we can thus spare ourselves to
// copy the objects for the intensive quantities.
// However, if the storage term at the start of the timestep cannot be deduced
// from the primary variables, we must calculate it from the old intensive
// quantities, and need to shift them.
return;
}
const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
for (unsigned timeIdx = 0; timeIdx < intensiveHistorySize - numSlots; ++timeIdx) {
intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
}
// the cache for the most recent time indices do not need to be invalidated
// because the solution for them did not change (TODO: that assumes that there is
// no post-processing of the solution after a time step! fix it?)
}
bool enableStorageCache() const
{ return enableStorageCache_; }
void setEnableStorageCache(bool enableStorageCache)
{ enableStorageCache_= enableStorageCache; }
const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
{
if (!enableStorageCache_ ||
timeIdx >= historySize ||
!storageCacheUpToDate_[timeIdx][globalIdx]) {
throw std::logic_error("Cached storage is not available or up to date for the requested "
"global index and time index. Make sure storage cache is enabled "
"and the entry is valid before calling this method.");
}
return storageCache_[timeIdx][globalIdx];
}
void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
{
if (!enableStorageCache_ || timeIdx >= historySize) {
return;
}
storageCache_[timeIdx][globalIdx] = value;
storageCacheUpToDate_[timeIdx][globalIdx] = 1;
}
bool storageCacheIsUpToDate(unsigned globalIdx, unsigned timeIdx) const
{
if (!enableStorageCache_ || timeIdx >= historySize) {
return false;
}
return storageCacheUpToDate_[timeIdx][globalIdx] != 0;
}
void invalidateStorageCacheEntry(unsigned globalIdx, unsigned timeIdx) const
{
if (enableStorageCache_ && timeIdx < historySize) {
storageCacheUpToDate_[timeIdx][globalIdx] = 0;
}
}
void invalidateStorageCache(unsigned timeIdx) const
{
if (enableStorageCache_ && timeIdx < historySize) {
std::fill(storageCacheUpToDate_[timeIdx].begin(),
storageCacheUpToDate_[timeIdx].end(),
/*value=*/0);
}
}
void shiftStorageCache(unsigned numSlots = 1) const
{
// If we cannot recycle first iteration storage, it does not make sense to shift the storage cache.
if (enableStorageCache_ && !simulator_.problem().recycleFirstIterationStorage()) {
for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++timeIdx) {
storageCache_[timeIdx + numSlots] = storageCache_[timeIdx];
storageCacheUpToDate_[timeIdx + numSlots] = storageCacheUpToDate_[timeIdx];
}
// should we invalidate the cache for the most recent time indices? (see shiftIntensiveQuantityCache)
}
}
Scalar globalResidual(GlobalEqVector& dest,
const SolutionVector& u) const
{
mutableSolution(/*timeIdx=*/0) = u;
const Scalar res = asImp_().globalResidual(dest);
mutableSolution(/*timeIdx=*/0) = asImp_().solution(/*timeIdx=*/0);
return res;
}
Scalar globalResidual(GlobalEqVector& dest) const
{
dest = 0;
std::mutex mutex;
ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
#ifdef _OPENMP
#pragma omp parallel
#endif
{
// Attention: the variables below are thread specific and thus cannot be
// moved in front of the #pragma!
const unsigned threadId = ThreadManager::threadId();
ElementContext elemCtx(simulator_);
ElementIterator elemIt = threadedElemIt.beginParallel();
LocalEvalBlockVector residual, storageTerm;
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
const Element& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updateAll(elem);
residual.resize(elemCtx.numDof(/*timeIdx=*/0));
storageTerm.resize(elemCtx.numPrimaryDof(/*timeIdx=*/0));
asImp_().localResidual(threadId).eval(residual, elemCtx);
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
mutex.lock();
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
const unsigned globalI = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
dest[globalI][eqIdx] += Toolbox::value(residual[dofIdx][eqIdx]);
}
}
mutex.unlock();
}
}
// add up the residuals on the process borders
const auto sumHandle =
GridCommHandleFactory::template sumHandle<EqVector>(dest, asImp_().dofMapper());
gridView_.communicate(*sumHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
// calculate the square norm of the residual. this is not
// entirely correct, since the residual for the finite volumes
// which are on the boundary are counted once for every
// process. As often in life: shit happens (, we don't care)...
return std::sqrt(asImp_().gridView().comm().sum(dest.two_norm2()));
}
void globalStorage(EqVector& storage, unsigned timeIdx = 0) const
{
storage = 0;
std::mutex mutex;
ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
#ifdef _OPENMP
#pragma omp parallel
#endif
{
// Attention: the variables below are thread specific and thus cannot be
// moved in front of the #pragma!
const unsigned threadId = ThreadManager::threadId();
ElementContext elemCtx(simulator_);
ElementIterator elemIt = threadedElemIt.beginParallel();
LocalEvalBlockVector elemStorage;
// in this method, we need to disable the storage cache because we want to
// evaluate the storage term for other time indices than the most recent one
elemCtx.setEnableStorageCache(false);
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
const Element& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue; // ignore ghost and overlap elements
}
elemCtx.updateStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
elemStorage.resize(numPrimaryDof);
localResidual(threadId).evalStorage(elemStorage, elemCtx, timeIdx);
mutex.lock();
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
storage[eqIdx] += Toolbox::value(elemStorage[dofIdx][eqIdx]);
}
}
mutex.unlock();
}
}
storage = gridView_.comm().sum(storage);
}
void checkConservativeness([[maybe_unused]] Scalar tolerance = -1,
[[maybe_unused]] bool verbose = false) const
{
#ifndef NDEBUG
Scalar totalBoundaryArea(0.0);
Scalar totalVolume(0.0);
EvalEqVector totalRate(0.0);
// take the newton tolerance times the total volume of the grid if we're not
// given an explicit tolerance...
if (tolerance <= 0) {
tolerance =
simulator_.model().newtonMethod().tolerance() *
simulator_.model().gridTotalVolume() *
1000;
}
// we assume the implicit Euler time discretization for now...
assert(historySize == 2);
EqVector storageBeginTimeStep(0.0);
globalStorage(storageBeginTimeStep, /*timeIdx=*/1);
EqVector storageEndTimeStep(0.0);
globalStorage(storageEndTimeStep, /*timeIdx=*/0);
// calculate the rate at the boundary and the source rate
ElementContext elemCtx(simulator_);
elemCtx.setEnableStorageCache(false);
for (const auto& elem : elements(simulator_.gridView())) {
if (elem.partitionType() != Dune::InteriorEntity) {
continue; // ignore ghost and overlap elements
}
elemCtx.updateAll(elem);
// handle the boundary terms
if (elemCtx.onBoundary()) {
BoundaryContext boundaryCtx(elemCtx);
for (unsigned faceIdx = 0; faceIdx < boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); ++faceIdx) {
BoundaryRateVector values;
simulator_.problem().boundary(values,
boundaryCtx,
faceIdx,
/*timeIdx=*/0);
Valgrind::CheckDefined(values);
const unsigned dofIdx = boundaryCtx.interiorScvIndex(faceIdx, /*timeIdx=*/0);
const auto& insideIntQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
const Scalar bfArea =
boundaryCtx.boundarySegmentArea(faceIdx, /*timeIdx=*/0) *
insideIntQuants.extrusionFactor();
for (unsigned i = 0; i < values.size(); ++i) {
values[i] *= bfArea;
}
totalBoundaryArea += bfArea;
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
totalRate[eqIdx] += values[eqIdx];
}
}
}
// deal with the source terms
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
RateVector values;
simulator_.problem().source(values,
elemCtx,
dofIdx,
/*timeIdx=*/0);
Valgrind::CheckDefined(values);
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
Scalar dofVolume =
elemCtx.dofVolume(dofIdx, /*timeIdx=*/0) *
intQuants.extrusionFactor();
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
totalRate[eqIdx] += -dofVolume*Toolbox::value(values[eqIdx]);
}
totalVolume += dofVolume;
}
}
// summarize everything over all processes
const auto& comm = simulator_.gridView().comm();
totalRate = comm.sum(totalRate);
totalBoundaryArea = comm.sum(totalBoundaryArea);
totalVolume = comm.sum(totalVolume);
if (comm.rank() == 0) {
EqVector storageRate = storageBeginTimeStep;
storageRate -= storageEndTimeStep;
storageRate /= simulator_.timeStepSize();
if (verbose) {
std::cout << "storage at beginning of time step: " << storageBeginTimeStep << "\n";
std::cout << "storage at end of time step: " << storageEndTimeStep << "\n";
std::cout << "rate based on storage terms: " << storageRate << "\n";
std::cout << "rate based on source and boundary terms: " << totalRate << "\n";
std::cout << "difference in rates: ";
for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
std::cout << (storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) << " ";
}
std::cout << "\n";
}
for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
Scalar eps =
(std::abs(storageRate[eqIdx]) + Toolbox::value(totalRate[eqIdx])) * tolerance;
eps = std::max(tolerance, eps);
assert(std::abs(storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) <= eps);
}
}
#endif // NDEBUG
}
Scalar dofTotalVolume(unsigned globalIdx) const
{ return dofTotalVolume_[globalIdx]; }
bool isLocalDof(unsigned globalIdx) const
{ return isLocalDof_[globalIdx]; }
Scalar gridTotalVolume() const
{ return gridTotalVolume_; }
const SolutionVector& solution(unsigned timeIdx) const
{ return solution_[timeIdx]->blockVector(); }
SolutionVector& solution(unsigned timeIdx)
{ return solution_[timeIdx]->blockVector(); }
protected:
SolutionVector& mutableSolution(unsigned timeIdx) const
{ return solution_[timeIdx]->blockVector(); }
public:
const Linearizer& linearizer() const
{ return *linearizer_; }
Linearizer& linearizer()
{ return *linearizer_; }
const LocalLinearizer& localLinearizer(unsigned openMpThreadId) const
{ return localLinearizer_[openMpThreadId]; }
LocalLinearizer& localLinearizer(unsigned openMpThreadId)
{ return localLinearizer_[openMpThreadId]; }
const LocalResidual& localResidual(unsigned openMpThreadId) const
{ return asImp_().localLinearizer(openMpThreadId).localResidual(); }
LocalResidual& localResidual(unsigned openMpThreadId)
{ return asImp_().localLinearizer(openMpThreadId).localResidual(); }
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
{
const Scalar absPv = std::abs(asImp_().solution(/*timeIdx=*/1)[globalDofIdx][pvIdx]);
return 1.0 / std::max(absPv, 1.0);
}
Scalar eqWeight(unsigned, unsigned) const
{ return 1.0; }
Scalar relativeDofError(unsigned vertexIdx,
const PrimaryVariables& pv1,
const PrimaryVariables& pv2) const
{
Scalar result = 0.0;
for (unsigned j = 0; j < numEq; ++j) {
const Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
const Scalar eqErr = std::abs((pv1[j] - pv2[j])*weight);
//Scalar eqErr = std::abs(pv1[j] - pv2[j]);
//eqErr *= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
result = std::max(result, eqErr);
}
return result;
}
bool update()
{
const TimerGuard prePostProcessGuard(prePostProcessTimer_);
#ifndef NDEBUG
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
// Make sure that the primary variables are defined. Note that because of padding
// bytes, we can't just simply ask valgrind to check the whole solution vectors
// for definedness...
for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
asImp_().solution(timeIdx)[i].checkDefined();
}
}
#endif // NDEBUG
// make sure all timers are prestine
prePostProcessTimer_.halt();
linearizeTimer_.halt();
solveTimer_.halt();
updateTimer_.halt();
prePostProcessTimer_.start();
asImp_().updateBegin();
prePostProcessTimer_.stop();
bool converged = false;
try {
converged = newtonMethod_.apply();
}
catch(...) {
prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
linearizeTimer_ += newtonMethod_.linearizeTimer();
solveTimer_ += newtonMethod_.solveTimer();
updateTimer_ += newtonMethod_.updateTimer();
throw;
}
#ifndef NDEBUG
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
// Make sure that the primary variables are defined. Note that because of padding
// bytes, we can't just simply ask valgrind to check the whole solution vectors
// for definedness...
for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
asImp_().solution(timeIdx)[i].checkDefined();
}
}
#endif // NDEBUG
prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
linearizeTimer_ += newtonMethod_.linearizeTimer();
solveTimer_ += newtonMethod_.solveTimer();
updateTimer_ += newtonMethod_.updateTimer();
prePostProcessTimer_.start();
if (converged) {
asImp_().updateSuccessful();
}
else {
asImp_().updateFailed();
}
prePostProcessTimer_.stop();
#ifndef NDEBUG
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
// Make sure that the primary variables are defined. Note that because of padding
// bytes, we can't just simply ask valgrind to check the whole solution vectors
// for definedness...
for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
asImp_().solution(timeIdx)[i].checkDefined();
}
}
#endif // NDEBUG
return converged;
}
void syncOverlap()
{}
void updateBegin()
{}
void updateSuccessful()
{}
void adaptGrid()
{
throw std::invalid_argument("Grid adaptation need to be implemented for "
"specific settings of grid and function spaces");
}
void updateFailed()
{
// Reset the current solution to the one of the
// previous time step so that we can start the next
// update at a physically meaningful solution.
solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
#ifndef NDEBUG
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
// Make sure that the primary variables are defined. Note that because of padding
// bytes, we can't just simply ask valgrind to check the whole solution vectors
// for definedness...
for (std::size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
asImp_().solution(timeIdx)[i].checkDefined();
}
}
#endif // NDEBUG
}
void advanceTimeLevel()
{
// at this point we can adapt the grid
if (this->enableGridAdaptation_) {
asImp_().adaptGrid();
}
// make the current solution the previous one.
solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
// shift the storage cache by one position in the history
asImp_().shiftStorageCache(/*numSlots=*/1);
// shift the intensive quantities cache by one position in the
// history
asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
}
template <class Restarter>
void serialize(Restarter&)
{
throw std::runtime_error("Not implemented: The discretization chosen for this problem "
"does not support restart files. (serialize() method unimplemented)");
}
template <class Restarter>
void deserialize(Restarter&)
{
throw std::runtime_error("Not implemented: The discretization chosen for this problem "
"does not support restart files. (deserialize() method unimplemented)");
}
template <class DofEntity>
void serializeEntity(std::ostream& outstream,
const DofEntity& dof)
{
const unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
// write phase state
if (!outstream.good()) {
throw std::runtime_error("Could not serialize degree of freedom " +
std::to_string(dofIdx));
}
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
outstream << solution(/*timeIdx=*/0)[dofIdx][eqIdx] << " ";
}
}
template <class DofEntity>
void deserializeEntity(std::istream& instream,
const DofEntity& dof)
{
const unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!instream.good()) {
throw std::runtime_error("Could not deserialize degree of freedom " +
std::to_string(dofIdx));
}
instream >> solution(/*timeIdx=*/0)[dofIdx][eqIdx];
}
}
std::size_t numGridDof() const
{ throw std::logic_error("The discretization class must implement the numGridDof() method!"); }
std::size_t numAuxiliaryDof() const
{
return std::accumulate(auxEqModules_.begin(), auxEqModules_.end(),
std::size_t{0},
[](const auto acc, const auto& mod)
{ return acc + mod->numDofs(); });
}
std::size_t numTotalDof() const
{ return asImp_().numGridDof() + numAuxiliaryDof(); }
const DofMapper& dofMapper() const
{ throw std::logic_error("The discretization class must implement the dofMapper() method!"); }
const VertexMapper& vertexMapper() const
{ return vertexMapper_; }
const ElementMapper& elementMapper() const
{ return elementMapper_; }
void resetLinearizer ()
{
linearizer_ = std::make_unique<Linearizer>();
linearizer_->init(simulator_);
}
static std::string discretizationName()
{ return ""; }
std::string primaryVarName(unsigned pvIdx) const
{
std::ostringstream oss;
oss << "primary variable_" << pvIdx;
return oss.str();
}
std::string eqName(unsigned eqIdx) const
{
std::ostringstream oss;
oss << "equation_" << eqIdx;
return oss.str();
}
void updatePVWeights(const ElementContext&) const
{}
void addOutputModule(std::unique_ptr<BaseOutputModule<TypeTag>> newModule)
{ outputModules_.push_back(std::move(newModule)); }
template <class VtkMultiWriter>
void addConvergenceVtkFields(VtkMultiWriter& writer,
const SolutionVector& u,
const GlobalEqVector& deltaU) const
{
using ScalarBuffer = std::vector<double>;
GlobalEqVector globalResid(u.size());
asImp_().globalResidual(globalResid, u);
// create the required scalar fields
const std::size_t numDof = asImp_().numGridDof();
// global defect of the two auxiliary equations
std::array<ScalarBuffer*, numEq> def;
std::array<ScalarBuffer*, numEq> delta;
std::array<ScalarBuffer*, numEq> priVars;
std::array<ScalarBuffer*, numEq> priVarWeight;
ScalarBuffer* relError = writer.allocateManagedScalarBuffer(numDof);
ScalarBuffer* normalizedRelError = writer.allocateManagedScalarBuffer(numDof);
for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
priVars[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
priVarWeight[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
delta[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
def[pvIdx] = writer.allocateManagedScalarBuffer(numDof);
}
Scalar minRelErr = 1e30;
Scalar maxRelErr = -1e30;
for (unsigned globalIdx = 0; globalIdx < numDof; ++ globalIdx) {
for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
(*priVars[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
(*priVarWeight[pvIdx])[globalIdx] = asImp_().primaryVarWeight(globalIdx, pvIdx);
(*delta[pvIdx])[globalIdx] = - deltaU[globalIdx][pvIdx];
(*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
}
PrimaryVariables uOld(u[globalIdx]);
PrimaryVariables uNew(uOld);
uNew -= deltaU[globalIdx];
const Scalar err = asImp_().relativeDofError(globalIdx, uOld, uNew);
(*relError)[globalIdx] = err;
(*normalizedRelError)[globalIdx] = err;
minRelErr = std::min(err, minRelErr);
maxRelErr = std::max(err, maxRelErr);
}
// do the normalization of the relative error
const Scalar alpha = std::max(Scalar{1e-20},
std::max(std::abs(maxRelErr),
std::abs(minRelErr)));
for (unsigned globalIdx = 0; globalIdx < numDof; ++globalIdx) {
(*normalizedRelError)[globalIdx] /= alpha;
}
DiscBaseOutputModule::attachScalarDofData_(writer, *relError, "relative error");
DiscBaseOutputModule::attachScalarDofData_(writer, *normalizedRelError, "normalized relative error");
for (unsigned i = 0; i < numEq; ++i) {
std::ostringstream oss;
oss.str(""); oss << "priVar_" << asImp_().primaryVarName(i);
DiscBaseOutputModule::attachScalarDofData_(writer,
*priVars[i],
oss.str());
oss.str(""); oss << "delta_" << asImp_().primaryVarName(i);
DiscBaseOutputModule::attachScalarDofData_(writer,
*delta[i],
oss.str());
oss.str(""); oss << "weight_" << asImp_().primaryVarName(i);
DiscBaseOutputModule::attachScalarDofData_(writer,
*priVarWeight[i],
oss.str());
oss.str(""); oss << "defect_" << asImp_().eqName(i);
DiscBaseOutputModule::attachScalarDofData_(writer,
*def[i],
oss.str());
}
asImp_().prepareOutputFields();
asImp_().appendOutputFields(writer);
}
void prepareOutputFields() const
{
const bool needFullContextUpdate =
std::any_of(outputModules_.begin(), outputModules_.end(),
[](const auto& mod) { return mod->needExtensiveQuantities(); });
std::for_each(outputModules_.begin(), outputModules_.end(),
[](auto& mod) { mod->allocBuffers(); });
// iterate over grid
ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
#ifdef _OPENMP
#pragma omp parallel
#endif
{
ElementContext elemCtx(simulator_);
ElementIterator elemIt = threadedElemIt.beginParallel();
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
// ignore non-interior entities
continue;
}
if (needFullContextUpdate) {
elemCtx.updateAll(elem);
}
else {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
}
std::for_each(outputModules_.begin(), outputModules_.end(),
[&elemCtx](auto& mod) { mod->processElement(elemCtx); });
}
}
}
void appendOutputFields(BaseOutputWriter& writer) const
{
std::for_each(outputModules_.begin(), outputModules_.end(),
[&writer](auto& mod) { mod->commitBuffers(writer); });
}
const GridView& gridView() const
{ return gridView_; }
void addAuxiliaryModule(BaseAuxiliaryModule<TypeTag>* auxMod)
{
auxMod->setDofOffset(numTotalDof());
auxEqModules_.push_back(auxMod);
// resize the solutions
if (enableGridAdaptation_ && !std::is_same_v<DiscreteFunction, BlockVectorWrapper>) {
throw std::invalid_argument("Problems which require auxiliary modules cannot be used in"
" conjunction with dune-fem");
}
const std::size_t numDof = numTotalDof();
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
solution(timeIdx).resize(numDof);
}
auxMod->applyInitial();
}
void clearAuxiliaryModules()
{
auxEqModules_.clear();
linearizer_->eraseMatrix();
newtonMethod_.eraseMatrix();
}
std::size_t numAuxiliaryModules() const
{ return auxEqModules_.size(); }
BaseAuxiliaryModule<TypeTag>* auxiliaryModule(unsigned auxEqModIdx)
{ return auxEqModules_[auxEqModIdx]; }
const BaseAuxiliaryModule<TypeTag>* auxiliaryModule(unsigned auxEqModIdx) const
{ return auxEqModules_[auxEqModIdx]; }
bool storeIntensiveQuantities() const
{ return enableIntensiveQuantityCache_ || enableThermodynamicHints_; }
const Timer& prePostProcessTimer() const
{ return prePostProcessTimer_; }
const Timer& linearizeTimer() const
{ return linearizeTimer_; }
const Timer& solveTimer() const
{ return solveTimer_; }
const Timer& updateTimer() const
{ return updateTimer_; }
template<class Serializer>
void serializeOp(Serializer& serializer)
{
using BaseDiscretization = GetPropType<TypeTag, Properties::BaseDiscretizationType>;
using Helper = typename BaseDiscretization::template SerializeHelper<Serializer>;
Helper::serializeOp(serializer, solution_);
}
bool operator==(const FvBaseDiscretization& rhs) const
{
return std::equal(this->solution_.begin(), this->solution_.end(),
rhs.solution_.begin(), rhs.solution_.end(),
[](const auto& x, const auto& y)
{
return *x == *y;
});
}
protected:
void resizeAndResetIntensiveQuantitiesCache_()
{
// allocate the storage cache
if (enableStorageCache()) {
const std::size_t numDof = asImp_().numGridDof();
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
storageCache_[timeIdx].resize(numDof);
storageCacheUpToDate_[timeIdx].resize(numDof, /*value=*/0);
}
}
// allocate the intensive quantities cache
if (storeIntensiveQuantities()) {
const std::size_t numDof = asImp_().numGridDof();
cachedIntensiveQuantityHistorySize_ = simulator_.problem().intensiveQuantityHistorySize();
const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_;
// resize the vectors based on runtime history size
intensiveQuantityCache_.resize(intensiveHistorySize);
intensiveQuantityCacheUpToDate_.resize(intensiveHistorySize);
for(unsigned timeIdx = 0; timeIdx < intensiveHistorySize; ++timeIdx) {
intensiveQuantityCache_[timeIdx].resize(numDof);
intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
invalidateIntensiveQuantitiesCache(timeIdx);
}
}
}
template <class Context>
void supplementInitialSolution_(PrimaryVariables&,
const Context&,
unsigned,
unsigned)
{}
void registerOutputModules_()
{
// add the output modules available on all model
this->outputModules_.push_back(std::make_unique<VtkPrimaryVarsModule<TypeTag>>(simulator_));
}
LocalResidual& localResidual_()
{ return localLinearizer_.localResidual(); }
bool verbose_() const
{ return gridView_.comm().rank() == 0; }
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
// the problem we want to solve. defines the constitutive
// relations, matxerial laws, etc.
Simulator& simulator_;
// the representation of the spatial domain of the problem
GridView gridView_;
// the mappers for element and vertex entities to global indices
ElementMapper elementMapper_;
VertexMapper vertexMapper_;
// a vector with all auxiliary equations to be considered
std::vector<BaseAuxiliaryModule<TypeTag>*> auxEqModules_;
NewtonMethod newtonMethod_;
Timer prePostProcessTimer_;
Timer linearizeTimer_;
Timer solveTimer_;
Timer updateTimer_;
// calculates the local jacobian matrix for a given element
std::vector<LocalLinearizer> localLinearizer_;
// Linearizes the problem at the current time step using the
// local jacobian
std::unique_ptr<Linearizer> linearizer_;
// cur is the current iterative solution, prev the converged
// solution of the previous time step
mutable std::vector<IntensiveQuantitiesVector> intensiveQuantityCache_;
// while these are logically bools, concurrent writes to vector<bool> are not thread safe.
mutable std::vector<std::vector<unsigned char>> intensiveQuantityCacheUpToDate_;
std::array<std::unique_ptr<DiscreteFunction>, historySize> solution_;
std::list<std::unique_ptr<BaseOutputModule<TypeTag>>> outputModules_;
Scalar gridTotalVolume_;
std::vector<Scalar> dofTotalVolume_;
std::vector<bool> isLocalDof_;
mutable std::array<GlobalEqVector, historySize> storageCache_;
// while these are logically bools, concurrent writes to vector<bool> are not thread safe.
mutable std::array<std::vector<unsigned char>, historySize> storageCacheUpToDate_;
bool enableGridAdaptation_;
bool enableIntensiveQuantityCache_;
bool enableStorageCache_;
bool enableThermodynamicHints_;
mutable unsigned cachedIntensiveQuantityHistorySize_;
};
template<class TypeTag>
class FvBaseDiscretizationNoAdapt : public FvBaseDiscretization<TypeTag>
{
using ParentType = FvBaseDiscretization<TypeTag>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using DiscreteFunction = GetPropType<TypeTag, Properties::DiscreteFunction>;
static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
public:
template<class Serializer>
struct SerializeHelper
{
template<class SolutionType>
static void serializeOp(Serializer& serializer,
SolutionType& solution)
{
for (auto& sol : solution) {
serializer(*sol);
}
}
};
explicit FvBaseDiscretizationNoAdapt(Simulator& simulator)
: ParentType(simulator)
{
if (this->enableGridAdaptation_) {
throw std::invalid_argument("Grid adaptation need to use"
" BaseDiscretization = FvBaseDiscretizationFemAdapt"
" which currently requires the presence of the"
" dune-fem module");
}
const std::size_t numDof = this->asImp_().numGridDof();
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", numDof);
}
}
};
} // namespace Opm
#endif // EWOMS_FV_BASE_DISCRETIZATION_HH
alignedallocator.hh
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1....
baseauxiliarymodule.hh
Base class for specifying auxiliary equations.
Opm::BaseAuxiliaryModule
Base class for specifying auxiliary equations.
Definition baseauxiliarymodule.hh:56
Opm::BaseAuxiliaryModule::applyInitial
virtual void applyInitial()=0
Set the initial condition of the auxiliary module in the solution vector.
Opm::BaseAuxiliaryModule::setDofOffset
void setDofOffset(int value)
Set the offset in the global system of equations for the first degree of freedom of this auxiliary mo...
Definition baseauxiliarymodule.hh:78
Opm::BaseOutputModule
The base class for writer modules.
Definition baseoutputmodule.hh:68
Opm::BaseOutputWriter
The base class for all output writers.
Definition baseoutputwriter.hh:46
Opm::FvBaseBoundaryContext
Represents all quantities which available on boundary segments.
Definition fvbaseboundarycontext.hh:46
Opm::FvBaseConstraintsContext
Represents all quantities which available for calculating constraints.
Definition fvbaseconstraintscontext.hh:44
Opm::FvBaseConstraints
Class to specify constraints for a finite volume spatial discretization.
Definition fvbaseconstraints.hh:48
Opm::FvBaseDiscretizationNoAdapt
The base class for the finite volume discretization schemes without adaptation.
Definition fvbasediscretization.hh:2048
Opm::FvBaseDiscretization::BlockVectorWrapper
Definition fvbasediscretization.hh:350
Opm::FvBaseDiscretization
The base class for the finite volume discretization schemes.
Definition fvbasediscretization.hh:298
Opm::FvBaseDiscretization::localLinearizer
LocalLinearizer & localLinearizer(unsigned openMpThreadId)
Definition fvbasediscretization.hh:1269
Opm::FvBaseDiscretization::prepareOutputFields
void prepareOutputFields() const
Prepare the quantities relevant for the current solution to be appended to the output writers.
Definition fvbasediscretization.hh:1769
Opm::FvBaseDiscretization::shiftIntensiveQuantityCache
void shiftIntensiveQuantityCache(unsigned numSlots=1)
Move the intensive quantities for a given time index to the back.
Definition fvbasediscretization.hh:785
Opm::FvBaseDiscretization::adaptGrid
void adaptGrid()
Called by the update() method when the grid should be refined.
Definition fvbasediscretization.hh:1442
Opm::FvBaseDiscretization::addAuxiliaryModule
void addAuxiliaryModule(BaseAuxiliaryModule< TypeTag > *auxMod)
Add a module for an auxiliary equation.
Definition fvbasediscretization.hh:1833
Opm::FvBaseDiscretization::setIntensiveQuantitiesCacheEntryValidity
void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx, unsigned timeIdx, bool newValue) const
Invalidate the cache for a given intensive quantities object.
Definition fvbasediscretization.hh:692
Opm::FvBaseDiscretization::finishInit
void finishInit()
Apply the initial conditions to the model.
Definition fvbasediscretization.hh:458
Opm::FvBaseDiscretization::prefetch
void prefetch(const Element &) const
Allows to improve the performance by prefetching all data which is associated with a given element.
Definition fvbasediscretization.hh:591
Opm::FvBaseDiscretization::updateSuccessful
void updateSuccessful()
Called by the update() method if it was successful.
Definition fvbasediscretization.hh:1436
Opm::FvBaseDiscretization::discretizationName
static std::string discretizationName()
Returns a string of discretization's human-readable name.
Definition fvbasediscretization.hh:1629
Opm::FvBaseDiscretization::cachedIntensiveQuantityHistorySize
unsigned cachedIntensiveQuantityHistorySize() const
Get the cached intensive quantity history size.
Definition fvbasediscretization.hh:706
Opm::FvBaseDiscretization::auxiliaryModule
BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx)
Returns a given module for auxiliary equations.
Definition fvbasediscretization.hh:1873
Opm::FvBaseDiscretization::isLocalDof
bool isLocalDof(unsigned globalIdx) const
Returns if the overlap of the volume ofa degree of freedom is non-zero.
Definition fvbasediscretization.hh:1209
Opm::FvBaseDiscretization::numAuxiliaryModules
std::size_t numAuxiliaryModules() const
Returns the number of modules for auxiliary equations.
Definition fvbasediscretization.hh:1867
Opm::FvBaseDiscretization::localResidual_
LocalResidual & localResidual_()
Reference to the local residal object.
Definition fvbasediscretization.hh:1971
Opm::FvBaseDiscretization::registerOutputModules_
void registerOutputModules_()
Register all output modules which make sense for the model.
Definition fvbasediscretization.hh:1962
Opm::FvBaseDiscretization::enableGridAdaptation
bool enableGridAdaptation() const
Returns whether the grid ought to be adapted to the solution during the simulation.
Definition fvbasediscretization.hh:523
Opm::FvBaseDiscretization::newtonMethod
const NewtonMethod & newtonMethod() const
Returns the newton method object.
Definition fvbasediscretization.hh:605
Opm::FvBaseDiscretization::mutableSolution
SolutionVector & mutableSolution(unsigned timeIdx) const
Definition fvbasediscretization.hh:1237
Opm::FvBaseDiscretization::advanceTimeLevel
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition fvbasediscretization.hh:1480
Opm::FvBaseDiscretization::newtonMethod
NewtonMethod & newtonMethod()
Returns the newton method object.
Definition fvbasediscretization.hh:599
Opm::FvBaseDiscretization::vertexMapper
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition fvbasediscretization.hh:1607
Opm::FvBaseDiscretization::cachedStorage
const EqVector & cachedStorage(unsigned globalIdx, unsigned timeIdx) const
Retrieve an entry of the cache for the storage term.
Definition fvbasediscretization.hh:843
Opm::FvBaseDiscretization::updatePVWeights
void updatePVWeights(const ElementContext &) const
Update the weights of all primary variables within an element given the complete set of intensive qua...
Definition fvbasediscretization.hh:1662
Opm::FvBaseDiscretization::updateFailed
void updateFailed()
Called by the update() method if it was unsuccessful.
Definition fvbasediscretization.hh:1453
Opm::FvBaseDiscretization::updateBegin
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition fvbasediscretization.hh:1429
Opm::FvBaseDiscretization::globalResidual
Scalar globalResidual(GlobalEqVector &dest) const
Compute the global residual for the current solution vector.
Definition fvbasediscretization.hh:964
Opm::FvBaseDiscretization::localResidual
LocalResidual & localResidual(unsigned openMpThreadId)
Definition fvbasediscretization.hh:1281
Opm::FvBaseDiscretization::serializeEntity
void serializeEntity(std::ostream &outstream, const DofEntity &dof)
Write the current solution for a degree of freedom to a restart file.
Definition fvbasediscretization.hh:1535
Opm::FvBaseDiscretization::setEnableStorageCache
void setEnableStorageCache(bool enableStorageCache)
Set the value of enable storage cache.
Definition fvbasediscretization.hh:827
Opm::FvBaseDiscretization::eqName
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition fvbasediscretization.hh:1649
Opm::FvBaseDiscretization::eqWeight
Scalar eqWeight(unsigned, unsigned) const
Returns the relative weight of an equation.
Definition fvbasediscretization.hh:1303
Opm::FvBaseDiscretization::primaryVarWeight
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition fvbasediscretization.hh:1291
Opm::FvBaseDiscretization::deserialize
void deserialize(Restarter &)
Deserializes the state of the model.
Definition fvbasediscretization.hh:1520
Opm::FvBaseDiscretization::checkConservativeness
void checkConservativeness(Scalar tolerance=-1, bool verbose=false) const
Ensure that the difference between the storage terms of the last and of the current time step is cons...
Definition fvbasediscretization.hh:1079
Opm::FvBaseDiscretization::gridTotalVolume
Scalar gridTotalVolume() const
Returns the volume of the whole grid which represents the spatial domain.
Definition fvbasediscretization.hh:1216
Opm::FvBaseDiscretization::thermodynamicHint
const IntensiveQuantities * thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
Return the thermodynamic hint for a entity on the grid at given time.
Definition fvbasediscretization.hh:623
Opm::FvBaseDiscretization::storageCacheIsUpToDate
bool storageCacheIsUpToDate(unsigned globalIdx, unsigned timeIdx) const
Returns true if the storage cache entry for a given DOF and time index is up to date.
Definition fvbasediscretization.hh:883
Opm::FvBaseDiscretization::linearizer
const Linearizer & linearizer() const
Returns the operator linearizer for the global jacobian of the problem.
Definition fvbasediscretization.hh:1245
Opm::FvBaseDiscretization::updateCachedStorage
void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector &value) const
Set an entry of the cache for the storage term.
Definition fvbasediscretization.hh:867
Opm::FvBaseDiscretization::addConvergenceVtkFields
void addConvergenceVtkFields(VtkMultiWriter &writer, const SolutionVector &u, const GlobalEqVector &deltaU) const
Add the vector fields for analysing the convergence of the newton method to the a VTK writer.
Definition fvbasediscretization.hh:1680
Opm::FvBaseDiscretization::relativeDofError
Scalar relativeDofError(unsigned vertexIdx, const PrimaryVariables &pv1, const PrimaryVariables &pv2) const
Returns the relative error between two vectors of primary variables.
Definition fvbasediscretization.hh:1315
Opm::FvBaseDiscretization::globalResidual
Scalar globalResidual(GlobalEqVector &dest, const SolutionVector &u) const
Compute the global residual for an arbitrary solution vector.
Definition fvbasediscretization.hh:949
Opm::FvBaseDiscretization::solution
SolutionVector & solution(unsigned timeIdx)
Definition fvbasediscretization.hh:1230
Opm::FvBaseDiscretization::deserializeEntity
void deserializeEntity(std::istream &instream, const DofEntity &dof)
Reads the current solution variables for a degree of freedom from a restart file.
Definition fvbasediscretization.hh:1560
Opm::FvBaseDiscretization::localLinearizer
const LocalLinearizer & localLinearizer(unsigned openMpThreadId) const
Returns the local jacobian which calculates the local stiffness matrix for an arbitrary element.
Definition fvbasediscretization.hh:1263
Opm::FvBaseDiscretization::update
bool update()
Try to progress the model to the next timestep.
Definition fvbasediscretization.hh:1336
Opm::FvBaseDiscretization::numAuxiliaryDof
std::size_t numAuxiliaryDof() const
Returns the number of degrees of freedom (DOFs) of the auxiliary equations.
Definition fvbasediscretization.hh:1583
Opm::FvBaseDiscretization::clearAuxiliaryModules
void clearAuxiliaryModules()
Causes the list of auxiliary equations to be cleared.
Definition fvbasediscretization.hh:1857
Opm::FvBaseDiscretization::elementMapper
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition fvbasediscretization.hh:1613
Opm::FvBaseDiscretization::syncOverlap
void syncOverlap()
Syncronize the values of the primary variables on the degrees of freedom that overlap with the neighb...
Definition fvbasediscretization.hh:1421
Opm::FvBaseDiscretization::appendOutputFields
void appendOutputFields(BaseOutputWriter &writer) const
Append the quantities relevant for the current solution to an output writer.
Definition fvbasediscretization.hh:1810
Opm::FvBaseDiscretization::numTotalDof
std::size_t numTotalDof() const
Returns the total number of degrees of freedom (i.e., grid plux auxiliary DOFs).
Definition fvbasediscretization.hh:1594
Opm::FvBaseDiscretization::applyInitialSolution
void applyInitialSolution()
Applies the initial solution for all degrees of freedom to which the model applies.
Definition fvbasediscretization.hh:530
Opm::FvBaseDiscretization::updateCachedIntensiveQuantities
void updateCachedIntensiveQuantities(const IntensiveQuantities &intQuants, unsigned globalIdx, unsigned timeIdx) const
Update the intensive quantity cache for a entity on the grid at given time.
Definition fvbasediscretization.hh:672
Opm::FvBaseDiscretization::primaryVarName
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition fvbasediscretization.hh:1637
Opm::FvBaseDiscretization::invalidateIntensiveQuantitiesCache
void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
Invalidate the whole intensive quantity cache for time index.
Definition fvbasediscretization.hh:714
Opm::FvBaseDiscretization::linearizer
Linearizer & linearizer()
Returns the object which linearizes the global system of equations at the current solution.
Definition fvbasediscretization.hh:1252
Opm::FvBaseDiscretization::auxiliaryModule
const BaseAuxiliaryModule< TypeTag > * auxiliaryModule(unsigned auxEqModIdx) const
Returns a given module for auxiliary equations.
Definition fvbasediscretization.hh:1879
Opm::FvBaseDiscretization::globalStorage
void globalStorage(EqVector &storage, unsigned timeIdx=0) const
Compute the integral over the domain of the storage terms of all conservation quantities.
Definition fvbasediscretization.hh:1024
Opm::FvBaseDiscretization::dofTotalVolume
Scalar dofTotalVolume(unsigned globalIdx) const
Returns the volume of a given control volume.
Definition fvbasediscretization.hh:1201
Opm::FvBaseDiscretization::registerParameters
static void registerParameters()
Register all run-time parameters for the model.
Definition fvbasediscretization.hh:427
Opm::FvBaseDiscretization::solution
const SolutionVector & solution(unsigned timeIdx) const
Reference to the solution at a given history index as a block vector.
Definition fvbasediscretization.hh:1224
Opm::FvBaseDiscretization::verbose_
bool verbose_() const
Returns whether messages should be printed.
Definition fvbasediscretization.hh:1977
Opm::FvBaseDiscretization::shiftStorageCache
void shiftStorageCache(unsigned numSlots=1) const
Shift storage cache by a given number of time step slots.
Definition fvbasediscretization.hh:929
Opm::FvBaseDiscretization::localResidual
const LocalResidual & localResidual(unsigned openMpThreadId) const
Returns the object to calculate the local residual function.
Definition fvbasediscretization.hh:1275
Opm::FvBaseDiscretization::numGridDof
std::size_t numGridDof() const
Returns the number of degrees of freedom (DOFs) for the computational grid.
Definition fvbasediscretization.hh:1577
Opm::FvBaseDiscretization::dofMapper
const DofMapper & dofMapper() const
Mapper to convert the Dune entities of the discretization's degrees of freedoms are to indices.
Definition fvbasediscretization.hh:1601
Opm::FvBaseDiscretization::storeIntensiveQuantities
bool storeIntensiveQuantities() const
Returns true if the cache for intensive quantities is enabled.
Definition fvbasediscretization.hh:1885
Opm::FvBaseDiscretization::invalidateStorageCache
void invalidateStorageCache(unsigned timeIdx) const
Invalidate the whole storage cache for a given time index.
Definition fvbasediscretization.hh:909
Opm::FvBaseDiscretization::serialize
void serialize(Restarter &)
Serializes the current state of the model.
Definition fvbasediscretization.hh:1506
Opm::FvBaseDiscretization::invalidateStorageCacheEntry
void invalidateStorageCacheEntry(unsigned globalIdx, unsigned timeIdx) const
Invalidate the storage cache for a given DOF and time index.
Definition fvbasediscretization.hh:897
Opm::FvBaseDiscretization::resetLinearizer
void resetLinearizer()
Resets the Jacobian matrix linearizer, so that the boundary types can be altered.
Definition fvbasediscretization.hh:1620
Opm::FvBaseDiscretization::enableStorageCache
bool enableStorageCache() const
Returns true iff the storage term is cached.
Definition fvbasediscretization.hh:818
Opm::FvBaseDiscretization::gridView
const GridView & gridView() const
Reference to the grid view of the spatial domain.
Definition fvbasediscretization.hh:1819
Opm::FvBaseDiscretization::addOutputModule
void addOutputModule(std::unique_ptr< BaseOutputModule< TypeTag > > newModule)
Add an module for writing visualization output after a timestep.
Definition fvbasediscretization.hh:1668
Opm::FvBaseDiscretization::cachedIntensiveQuantities
const IntensiveQuantities * cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
Return the cached intensive quantities for a entity on the grid at given time.
Definition fvbasediscretization.hh:644
Opm::FvBaseElementContext
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
Definition fvbaseelementcontext.hh:55
Opm::FvBaseExtensiveQuantities
Provide the properties at a face which make sense independently of the conserved quantities.
Definition fvbaseextensivequantities.hh:48
Opm::FvBaseGradientCalculator
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition fvbasegradientcalculator.hh:52
Opm::FvBaseIntensiveQuantities
Base class for the model specific class which provides access to all intensive (i....
Definition fvbaseintensivequantities.hh:45
Opm::FvBaseLinearizer
The common code for the linearizers of non-linear systems of equations.
Definition fvbaselinearizer.hh:78
Opm::FvBaseLocalResidual
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition fvbaselocalresidual.hh:63
Opm::FvBasePrimaryVariables
Represents the primary variables used by the a model.
Definition fvbaseprimaryvariables.hh:53
Opm::Linear::NullBorderListCreator
This is a grid manager which does not create any border list.
Definition nullborderlistmanager.hh:44
Opm::NewtonMethod::registerParameters
static void registerParameters()
Register all run-time parameters for the Newton method.
Definition newtonmethod.hh:136
Opm::Simulator
Manages the initializing and running of time dependent problems.
Definition simulator.hh:84
Opm::ThreadManager
Simplifies multi-threaded capabilities.
Definition threadmanager.hpp:36
Opm::ThreadManager::maxThreads
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition threadmanager.hpp:66
Opm::ThreadManager::threadId
static unsigned threadId()
Return the index of the current OpenMP thread.
Definition threadmanager.cpp:84
Opm::ThreadedEntityIterator
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
Definition threadedentityiterator.hh:42
Opm::TimerGuard
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition timerguard.hh:42
Opm::Timer
Provides an encapsulation to measure the system time.
Definition timer.hpp:46
Opm::VtkMultiWriter
Simplifies writing multi-file VTK datasets.
Definition vtkmultiwriter.hh:65
Opm::VtkMultiWriter::allocateManagedScalarBuffer
ScalarBuffer * allocateManagedScalarBuffer(std::size_t numEntities)
Allocate a managed buffer for a scalar field.
Definition vtkmultiwriter.hh:192
Opm::VtkPrimaryVarsModule
VTK output module for the fluid composition.
Definition vtkprimaryvarsmodule.hpp:48
Opm::VtkPrimaryVarsModule::registerParameters
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkprimaryvarsmodule.hpp:74
Opm::aligned_allocator
Definition alignedallocator.hh:97
fvbaseadlocallinearizer.hh
Calculates the local residual and its Jacobian for a single element of the grid.
fvbaseboundarycontext.hh
Represents all quantities which available on boundary segments.
fvbaseconstraints.hh
Class to specify constraints for a finite volume spatial discretization.
fvbaseconstraintscontext.hh
Represents all quantities which available for calculating constraints.
fvbaseelementcontext.hh
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
fvbaseextensivequantities.hh
Provide the properties at a face which make sense independently of the conserved quantities.
fvbasefdlocallinearizer.hh
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
fvbasegradientcalculator.hh
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
fvbaseintensivequantities.hh
Base class for the model specific class which provides access to all intensive (i....
fvbaselinearizer.hh
The common code for the linearizers of non-linear systems of equations.
fvbaselocalresidual.hh
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
fvbasenewtonmethod.hh
A Newton method for models using a finite volume discretization.
fvbaseprimaryvariables.hh
Represents the primary variables used by the a model.
fvbaseproperties.hh
Declare the properties used by the infrastructure code of the finite volume discretizations.
gridcommhandles.hh
Provides data handles for parallel communication which operate on DOFs.
linalgparameters.hh
Declares the parameters for the black oil model.
Opm::Properties::TTag
The generic type tag for problems using the immiscible multi-phase model.
Definition blackoilmodel.hh:81
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:233
Opm::getPropValue
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:240
nullborderlistmanager.hh
This is a grid manager which does not create any border list.
Opm::Parameters::Register
void Register(const char *usageString)
Register a run-time parameter.
Definition parametersystem.hpp:292
simulator.hh
Manages the initializing and running of time dependent problems.
Opm::FvBaseDiscretizationNoAdapt::SerializeHelper
Definition fvbasediscretization.hh:2058
Opm::Properties::BaseDiscretizationType
Definition fvbasediscretization.hh:272
Opm::Properties::BorderListCreator
The class which marks the border indices associated with the degrees of freedom on a process boundary...
Definition basicproperties.hh:125
Opm::Properties::BoundaryContext
The secondary variables of a boundary segment.
Definition fvbaseproperties.hh:143
Opm::Properties::BoundaryRateVector
Type of object for specifying boundary conditions.
Definition fvbaseproperties.hh:119
Opm::Properties::ConstraintsContext
The secondary variables of a constraint degree of freedom.
Definition fvbaseproperties.hh:146
Opm::Properties::Constraints
The class which represents a constraint degree of freedom.
Definition fvbaseproperties.hh:122
Opm::Properties::DiscExtensiveQuantities
The part of the extensive quantities which is specific to the spatial discretization.
Definition fvbaseproperties.hh:160
Opm::Properties::DiscIntensiveQuantities
The discretization specific part of the intensive quantities.
Definition fvbaseproperties.hh:136
Opm::Properties::DiscLocalResidual
The discretization specific part of the local residual.
Definition fvbaseproperties.hh:91
Opm::Properties::DiscreteFunction
Definition fvbaseproperties.hh:77
Opm::Properties::ElementContext
The secondary variables of all degrees of freedom in an element's stencil.
Definition fvbaseproperties.hh:140
Opm::Properties::ElementEqVector
A vector of holding a quantity for each equation for each DOF of an element.
Definition fvbaseproperties.hh:112
Opm::Properties::ElementMapper
The mapper to find the global index of an element.
Definition fvbaseproperties.hh:213
Opm::Properties::EnableConstraints
Specify whether the some degrees of fredom can be constraint.
Definition fvbaseproperties.hh:199
Opm::Properties::EnableExperiments
Specify if experimental features should be enabled or not.
Definition fvbaseproperties.hh:241
Opm::Properties::EqVector
A vector of holding a quantity for each equation (usually at a given spatial location).
Definition fvbaseproperties.hh:109
Opm::Properties::ExtensiveStorageTerm
Specify whether the storage terms use extensive quantities or not.
Definition fvbaseproperties.hh:233
Opm::Properties::GlobalEqVector
Vector containing a quantity of for equation for each DOF of the whole grid.
Definition linalgproperties.hh:54
Opm::Properties::GradientCalculator
Calculates gradients of arbitrary quantities at flux integration points.
Definition fvbaseproperties.hh:152
Opm::Properties::IntensiveQuantities
The secondary variables within a sub-control volume.
Definition fvbaseproperties.hh:133
Opm::Properties::Linearizer
The class which linearizes the non-linear system of equations.
Definition newtonmethodproperties.hh:36
Opm::Properties::PrimaryVariables
A vector of primary variables within a sub-control volume.
Definition fvbaseproperties.hh:130
Opm::Properties::RateVector
Vector containing volumetric or areal rates of quantities.
Definition fvbaseproperties.hh:116
Opm::Properties::Simulator
Manages the simulation time.
Definition basicproperties.hh:116
Opm::Properties::SolutionVector
Vector containing all primary variables of the grid.
Definition fvbaseproperties.hh:126
Opm::Properties::ThreadManager
The OpenMP threads manager.
Definition fvbaseproperties.hh:174
Opm::Properties::TimeDiscHistorySize
The history size required by the time discretization.
Definition fvbaseproperties.hh:225
Opm::Properties::UndefinedProperty
a tag to mark properties as undefined
Definition propertysystem.hh:38
Opm::Properties::UseLinearizationLock
use locking to prevent race conditions when linearizing the global system of equations in multi-threa...
Definition fvbaseproperties.hh:181
Opm::Properties::UseVolumetricResidual
Specify whether to use volumetric residuals or not.
Definition fvbaseproperties.hh:237
Opm::Properties::VertexMapper
The mapper to find the global index of a vertex.
Definition fvbaseproperties.hh:207
Opm::Properties::VtkOutputFormat
Specify the format the VTK output is written to disk.
Definition fvbaseproperties.hh:195
threadmanager.hpp
Simplifies multi-threaded capabilities.
timer.hpp
Provides an encapsulation to measure the system time.
timerguard.hh
A simple class which makes sure that a timer gets stopped if an exception is thrown.
vtkprimaryvarsmodule.hpp
VTK output module for the fluid composition.