Commit c427655a authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

add mpfa framework and implementation for the o method

parent 7e427364
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief This file contains the data which is required to calculate
* volume and mass fluxes of fluid phases over a face of a finite volume by means
* of the Darcy approximation. Specializations are provided for the different discretization methods.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
#include <memory>
#include <dune/common/float_cmp.hh>
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
namespace Dumux
{
namespace Properties
{
// forward declaration of properties
NEW_PROP_TAG(ProblemEnableGravity);
}
/*!
* \ingroup DarcysLaw
* \brief Specialization of Darcy's Law for the CCTpfa method.
*/
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FluxVarsCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
public:
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx,
const FluxVarsCache& fluxVarsCache)
{
const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil(phaseIdx);
const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions(phaseIdx);
const auto& tij = fluxVarsCache.advectionTij(phaseIdx);
// calculate Tij*pj
Scalar flux(0.0);
unsigned int localIdx = 0;
for (const auto volVarIdx : volVarsStencil)
{
const auto& volVars = elemVolVars[volVarIdx];
Scalar h = volVars.pressure(phaseIdx);
if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
{
// gravitational acceleration in the center of the actual element
const auto rho = volVars.density(phaseIdx);
const auto x = volVarsPositions[localIdx];
const auto g = problem.gravityAtPos(x);
h -= rho*(g*x);
}
flux += tij[localIdx++]*h;
}
return flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
}
static Stencil stencil(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
Stencil stencil;
if(problem.model().globalFvGeometry().scvfTouchesBoundary(scvf))
{
const auto& ivSeed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf, /*dummy*/0);
const auto& localScvSeeds = ivSeed.scvSeeds();
stencil.reserve(localScvSeeds.size());
for (const auto& localScvSeed : ivSeed.scvSeeds())
stencil.push_back(localScvSeed.globalIndex());
}
else
{
const auto& ivSeed = problem.model().globalFvGeometry().interactionVolumeSeed(scvf);
const auto& localScvSeeds = ivSeed.scvSeeds();
stencil.reserve(localScvSeeds.size());
for (const auto& localScvSeed : ivSeed.scvSeeds())
stencil.push_back(localScvSeed.globalIndex());
}
return stencil;
}
};
} // end namespace
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief The global object of flux var caches
*/
#ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_FLUXVARSCACHE_HH
#include <dumux/implicit/properties.hh>
namespace Dumux
{
/*!
* \ingroup ImplicitModel
* \brief Base class for the stencil local flux variables cache
*/
template<class TypeTag, bool EnableGlobalFluxVariablesCache>
class CCMpfaElementFluxVariablesCache;
/*!
* \ingroup ImplicitModel
* \brief Spezialization when caching globally
*/
template<class TypeTag>
class CCMpfaElementFluxVariablesCache<TypeTag, true>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using GlobalFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GlobalFluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
CCMpfaElementFluxVariablesCache(const GlobalFluxVariablesCache& global)
: globalFluxVarsCachePtr_(&global) {}
// Specialization for the global caching being enabled - do nothing here
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars) {}
// Specialization for the global caching being enabled - do nothing here
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars) {}
// Specialization for the global caching being enabled - do nothing here
void bindScvf(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) {}
// access operators in the case of caching
const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
{ return (*globalFluxVarsCachePtr_)[scvf.index()]; }
//! The global object we are a restriction of
const GlobalFluxVariablesCache& globalFluxVarsCache() const
{ return *globalFluxVarsCachePtr_; }
private:
const GlobalFluxVariablesCache* globalFluxVarsCachePtr_;
};
/*!
* \ingroup ImplicitModel
* \brief Spezialization when not using global caching
*/
template<class TypeTag>
class CCMpfaElementFluxVariablesCache<TypeTag, false>
{
// TODO: implement
};
} // end namespace
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief The local (stencil) volume variables class for cell centered models
*/
#ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH
#include <dumux/implicit/properties.hh>
namespace Dumux
{
/*!
* \ingroup ImplicitModel
* \brief Base class for the volume variables vector
*/
template<class TypeTag, bool enableGlobalVolVarsCache>
class CCMpfaElementVolumeVariables
{};
// specialization in case of storing the volume variables globally
template<class TypeTag>
class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
public:
//! Constructor
CCMpfaElementVolumeVariables(const GlobalVolumeVariables& globalVolVars)
: globalVolVarsPtr_(&globalVolVars) {}
const VolumeVariables& operator [](const SubControlVolume& scv) const
{ return globalVolVars().volVars(scv.index()); }
// operator for the access with an index
// needed for cc methods for the access to the boundary volume variables
const VolumeVariables& operator [](const IndexType scvIdx) const
{ return globalVolVars().volVars(scvIdx); }
// For compatibility reasons with the case of not storing the vol vars.
// function to be called before assembling an element, preparing the vol vars within the stencil
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{}
// function to prepare the vol vars within the element
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{}
//! The global volume variables object we are a restriction of
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
private:
const GlobalVolumeVariables* globalVolVarsPtr_;
};
// Specialization when the current volume variables are not stored
template<class TypeTag>
class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
public:
//! Constructor
CCMpfaElementVolumeVariables(const GlobalVolumeVariables& globalVolVars)
: globalVolVarsPtr_(&globalVolVars) {}
// Binding of an element, prepares the volume variables within the element stencil
// called by the local jacobian to prepare element assembly
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
auto eIdx = globalVolVars().problem_().elementMapper().index(element);
// stencil information
const auto& neighborStencil = globalVolVars().problem_().model().stencils(element).neighborStencil();
const auto numDofs = neighborStencil.size() + 1;
// resize local containers to the required size (for internal elements)
volumeVariables_.resize(numDofs);
volVarIndices_.resize(numDofs);
int localIdx = 0;
// update the volume variables of the element at hand
auto&& scvI = fvGeometry.scv(eIdx);
volumeVariables_[localIdx].update(sol[eIdx], globalVolVars().problem_(), element, scvI);
volVarIndices_[localIdx] = scvI.index();
++localIdx;
// Update the volume variables of the neighboring elements
for (auto globalJ : neighborStencil)
{
const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
volumeVariables_[localIdx].update(sol[globalJ], globalVolVars().problem_(), elementJ, scvJ);
volVarIndices_[localIdx] = scvJ.index();
++localIdx;
}
// Update boundary volume variables
for (auto&& scvf : scvfs(fvGeometry))
{
// if we are not on a boundary, skip to the next scvf
if (!scvf.boundary())
continue;
const auto dirichletPriVars = globalVolVars().problem_().dirichlet(element, scvf);
volumeVariables_.resize(localIdx+1);
volVarIndices_.resize(localIdx+1);
volumeVariables_[localIdx].update(dirichletPriVars, globalVolVars().problem_(), element, scvI);
volVarIndices_[localIdx] = scvf.outsideScvIdx();
++localIdx;
}
}
// Binding of an element, prepares only the volume variables of the element
// specialization for cc models
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
auto eIdx = globalVolVars().problem_().elementMapper().index(element);
volumeVariables_.resize(1);
volVarIndices_.resize(1);
// update the volume variables of the element
auto&& scv = fvGeometry.scv(eIdx);
volumeVariables_[0].update(sol[eIdx], globalVolVars().problem_(), element, scv);
volVarIndices_[0] = scv.index();
}
const VolumeVariables& operator [](const SubControlVolume& scv) const
{ return volumeVariables_[getLocalIdx_(scv.index())]; }
VolumeVariables& operator [](const SubControlVolume& scv)
{ return volumeVariables_[getLocalIdx_(scv.index())]; }
const VolumeVariables& operator [](IndexType scvIdx) const
{ return volumeVariables_[getLocalIdx_(scvIdx)]; }
VolumeVariables& operator [](IndexType scvIdx)
{ return volumeVariables_[getLocalIdx_(scvIdx)]; }
//! The global volume variables object we are a restriction of
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
private:
const GlobalVolumeVariables* globalVolVarsPtr_;
const int getLocalIdx_(const int volVarIdx) const
{
auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
assert(it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
return std::distance(volVarIndices_.begin(), it);
}
std::vector<IndexType> volVarIndices_;
std::vector<VolumeVariables> volumeVariables_;
};
} // end namespace
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Face types of the sub control volume faces in the mpfa method.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FACETYPES_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FACETYPES_HH
namespace Dumux
{
enum class MpfaFaceTypes : unsigned int
{
interior,
neumann,
dirichlet,
interiorNeumann,
interiorDirichlet
};
} // end namespace
#endif
\ No newline at end of file
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Base class for a local finite volume geometry for cell-centered TPFA models
* This builds up the sub control volumes and sub control volume faces
* for each element in the local scope we are restricting to, e.g. stencil or element.
*/
#ifndef DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH
#include <dune/common/iteratorrange.hh>
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/implicit/cellcentered/mpfa/properties.hh>
namespace Dumux
{
//! forward declaration of the global finite volume geometry
template<class TypeTag, bool EnableGlobalFVGeometryCache>
class CCMpfaGlobalFVGeometry;
/*!
* \ingroup ImplicitModel
* \brief Base class for the finite volume geometry vector for cell-centered TPFA models
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
template<class TypeTag, bool EnableGlobalFVGeometryCache>
class CCMpfaFVElementGeometry
{};
//! specialization in case the FVElementGeometries are stored globally
//! In this case we just forward internally to the global object
template<class TypeTag>
class CCMpfaFVElementGeometry<TypeTag, true>
{
using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
using ScvIterator = ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
using ScvfIterator = ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
public:
//! Constructor
CCMpfaFVElementGeometry(const GlobalFVGeometry& globalFvGeometry)
: globalFvGeometryPtr_(&globalFvGeometry) {}
//! Get an element sub control volume with a global scv index
const SubControlVolume& scv(IndexType scvIdx) const
{
return globalFvGeometry().scv(scvIdx);
}
//! Get an element sub control volume face with a global scvf index
const SubControlVolumeFace& scvf(IndexType scvfIdx) const
{
return globalFvGeometry().scvf(scvfIdx);
}
//! iterator range for sub control volumes. Iterates over
//! all scvs of the bound element (not including neighbor scvs)
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
friend inline Dune::IteratorRange<ScvIterator>
scvs(const CCMpfaFVElementGeometry& fvGeometry)
{
return Dune::IteratorRange<ScvIterator>(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry),
ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry));
}
//! iterator range for sub control volumes faces. Iterates over
//! all scvfs of the bound element (not including neighbor scvfs)
//! This is a free function found by means of ADL
//! To iterate over all sub control volume faces of this FVElementGeometry use