Commit 2f250712 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[staggeredGrid] Add more required classes, some cleanup

parent 80fb8768
......@@ -121,40 +121,11 @@ public:
// hOutside -= rho*(gOutside*xOutside);
// }
// }
//
// const GlobalPosition vector = {1,0};
// const Scalar angle = std::acos(scvFace.unitOuterNormal() * vector / scvFace.unitOuterNormal().two_norm() / vector.two_norm())*180.0/M_PI;
// const Scalar threshold = 1e-6;
//
// if(std::abs(angle) < threshold)
return 1.0;
// // else if(std::abs(angle-180) < threshold)
// // return -1.0;
// else
// return 0;
//
// // if(std::abs(angle) < threshold || std::abs(angle-180) < threshold)
// // return tij*(hInside - hOutside);
// // std::cout << "normal " << scvFace.unitOuterNormal() << std::endl;
//
}
// static Stencil pressureStencil(const Problem& problem,
// const Element& element,
// const FVElementGeometry& fvGeometry,
// const SubControlVolumeFace& scvFace)
// {
// Stencil pressureStencil;
// if (!scvFace.boundary())
// {
// pressureStencil.push_back(scvFace.insideScvIdx());
// pressureStencil.push_back(scvFace.outsideScvIdx());
// }
// else
// pressureStencil.push_back(scvFace.insideScvIdx());
//
// return pressureStencil;
// }
}
// The flux variables cache has to be bound to an element prior to flux calculations
// During the binding, the transmissibilities will be computed and stored using the method below.
......@@ -164,75 +135,10 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvFace)
{
Scalar tij;
const auto insideScvIdx = scvFace.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = elemVolVars[insideScvIdx];
const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv, insideVolVars);
Scalar ti = calculateOmega_(problem, scvFace, insideK, element, insideScv);
if (!scvFace.boundary())
{
const auto outsideScvIdx = scvFace.outsideScvIdx();
// as we assemble fluxes from the neighbor to our element the outside index
// refers to the scv of our element, so we use the scv method
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto outsideElement = fvGeometry.globalFvGeometry().element(outsideScvIdx);
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
const auto outsideK = problem.spatialParams().intrinsicPermeability(outsideScv, outsideVolVars);
Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideElement, outsideScv);
// check for division by zero!
if (ti*tj <= 0.0)
tij = 0;
else
tij = scvFace.area()*(ti * tj)/(ti + tj);
}
else
{
tij = scvFace.area()*ti;
}
return tij;
return 0;
}
private:
static Scalar calculateOmega_(const Problem& problem,
const SubControlVolumeFace& scvFace,
const DimWorldMatrix &K,
const Element& element,
const SubControlVolume &scv)
{
GlobalPosition Knormal;
K.mv(scvFace.unitOuterNormal(), Knormal);
auto distanceVector = scvFace.center();
distanceVector -= scv.center();
distanceVector /= distanceVector.two_norm2();
Scalar omega = Knormal * distanceVector;
omega *= problem.boxExtrusionFactor(element, scv);
return omega;
}
static Scalar calculateOmega_(const Problem& problem,
const SubControlVolumeFace& scvFace,
const Scalar K,
const Element& element,
const SubControlVolume &scv)
{
auto distanceVector = scvFace.center();
distanceVector -= scv.center();
distanceVector /= distanceVector.two_norm2();
Scalar omega = K * (distanceVector * scvFace.unitOuterNormal());
omega *= problem.boxExtrusionFactor(element, scv);
return omega;
}
};
} // end namespace
......
// -*- 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_STAGGERED_ELEMENT_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_STAGGERED_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 StaggeredElementFluxVariablesCache;
/*!
* \ingroup ImplicitModel
* \brief Spezialization when caching globally
*/
template<class TypeTag>
class StaggeredElementFluxVariablesCache<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:
StaggeredElementFluxVariablesCache(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) {}
// aStaggeredess 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 StaggeredElementFluxVariablesCache<TypeTag, false>
{
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:
StaggeredElementFluxVariablesCache(const GlobalFluxVariablesCache& global)
: globalFluxVarsCachePtr_(&global) {}
// This function has to be called prior to flux calculations on the element.
// Prepares the transmissibilities of the scv faces in an element. The FvGeometry is assumed to be bound.
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
// resizing of the cache
const auto numScvf = fvGeometry.numScvf();
fluxVarsCache_.resize(numScvf);
globalScvfIndices_.resize(numScvf);
IndexType localScvfIdx = 0;
// fill the containers
for (auto&& scvf : scvfs(fvGeometry))
{
fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
globalScvfIndices_[localScvfIdx] = scvf.index();
localScvfIdx++;
}
}
// This function is called by the StaggeredLocalResidual before flux calculations during assembly.
// Prepares the transmissibilities of the scv faces in the stencil. The FvGeometries are assumed to be bound.
void bind(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
const auto globalI = globalFluxVarsCache().problem_().elementMapper().index(element);
const auto& neighborStencil = globalFluxVarsCache().problem_().model().stencils(element).neighborStencil();
const auto numNeighbors = neighborStencil.size();
// find the number of scv faces that need to be prepared
auto numScvf = fvGeometry.numScvf();
for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
{
const auto& fluxVarIndicesJ = globalFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
numScvf += fluxVarIndicesJ.size();
}
// fill the containers with the data on the scv faces inside the actual element
fluxVarsCache_.resize(numScvf);
globalScvfIndices_.resize(numScvf);
IndexType localScvfIdx = 0;
for (auto&& scvf : scvfs(fvGeometry))
{
fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
globalScvfIndices_[localScvfIdx] = scvf.index();
localScvfIdx++;
}
// add required data on the scv faces in the neighboring elements
for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
{
const auto& fluxVarIndicesJ = globalFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
const auto elementJ = fvGeometry.globalFvGeometry().element(neighborStencil[localIdxJ]);
for (auto fluxVarIdx : fluxVarIndicesJ)
{
auto&& scvfJ = fvGeometry.scvf(fluxVarIdx);
fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), elementJ, fvGeometry, elemVolVars, scvfJ);
globalScvfIndices_[localScvfIdx] = scvfJ.index();
localScvfIdx++;
}
}
}
void bindScvf(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
{
fluxVarsCache_.resize(1);
globalScvfIndices_.resize(1);
fluxVarsCache_[0].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
globalScvfIndices_[0] = scvf.index();
}
// access operators in the case of no caching
const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
{ return fluxVarsCache_[getLocalScvfIdx_(scvf.index())]; }
FluxVariablesCache& operator [](const SubControlVolumeFace& scvf)
{ return fluxVarsCache_[getLocalScvfIdx_(scvf.index())]; }
//! The global object we are a restriction of
const GlobalFluxVariablesCache& globalFluxVarsCache() const
{ return *globalFluxVarsCachePtr_; }
private:
const GlobalFluxVariablesCache* globalFluxVarsCachePtr_;
// get index of scvf in the local container
int getLocalScvfIdx_(const int scvfIdx) const
{
auto it = std::find(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
assert(it != globalScvfIndices_.end() && "Could not find the flux vars cache for scvfIdx");
return std::distance(globalScvfIndices_.begin(), it);
}
std::vector<FluxVariablesCache> fluxVarsCache_;
std::vector<IndexType> globalScvfIndices_;
};
} // 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_STAGGERED_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_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 StaggeredElementVolumeVariables
{};
// specialization in case of storing the volume variables globally
template<class TypeTag>
class StaggeredElementVolumeVariables<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
StaggeredElementVolumeVariables(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 Staggered 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 StaggeredElementVolumeVariables<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
StaggeredElementVolumeVariables(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;
// check if boundary is a pure dirichlet boundary
const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
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 Staggered 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
......@@ -20,11 +20,11 @@
* \file
* \brief The global object of flux var caches
*/
#ifndef DUMUX_DISCRETIZATION_CC_GLOBAL_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_CC_GLOBAL_FLUXVARSCACHE_HH
#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FLUXVARSCACHE_HH
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/cellcentered/elementfluxvariablescache.hh>
#include <dumux/discretization/staggered/elementfluxvariablescache.hh>
namespace Dumux
{
......@@ -34,17 +34,17 @@ namespace Dumux
* \brief Base class for the flux variables cache vector, we store one cache per face
*/
template<class TypeTag, bool EnableGlobalFluxVariablesCache>
class CCGlobalFluxVariablesCache;
class StaggeredGlobalFluxVariablesCache;
/*!
* \ingroup ImplicitModel
* \brief Spezialization when caching globally
*/
template<class TypeTag>
class CCGlobalFluxVariablesCache<TypeTag, true>
class StaggeredGlobalFluxVariablesCache<TypeTag, true>
{
// the local class needs access to the problem
friend CCElementFluxVariablesCache<TypeTag, true>;
friend StaggeredElementFluxVariablesCache<TypeTag, true>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
......@@ -80,7 +80,7 @@ public:
* The local object is only functional after calling its bind/bindElement method