Skip to content
GitLab
Commit f4c85a46 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[staggered] Introduce elementFaceVars 

* can be used like elemVolVars: elemFaceVars[scvf]
parent f72e5874
Hide whitespace changes
Inline Side-by-side
dumux/assembly/staggeredlocalassembler.hh
View file @ f4c85a46
......@@ -78,6 +78,9 @@ class StaggeredLocalAssembler<TypeTag,
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
static constexpr bool enableGlobalFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache);
......@@ -111,8 +114,10 @@ public:
auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
auto& curGlobalFaceVars = gridVariables.curGridFaceVars();
auto& prevGlobalFaceVars = gridVariables.prevGridFaceVars();
auto curElemeFaceVars = localView(gridVariables.curGridFaceVars());
auto prevElemeFaceVars = localView(gridVariables.prevGridFaceVars());
curElemeFaceVars.bind(element, fvGeometry, curSol);
prevElemeFaceVars.bind(element, fvGeometry, curSol); //TODO: stationary
const bool isStationary = localResidual.isStationary();
auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
......@@ -128,8 +133,8 @@ public:
fvGeometry,
prevElemVolVars,
curElemVolVars,
prevGlobalFaceVars,
curGlobalFaceVars,
prevElemeFaceVars,
curElemeFaceVars,
elemBcTypes,
elemFluxVarsCache)[0];
......@@ -148,8 +153,8 @@ public:
scvf,
prevElemVolVars,
curElemVolVars,
prevGlobalFaceVars,
curGlobalFaceVars,
prevElemeFaceVars,
curElemeFaceVars,
elemBcTypes,
elemFluxVarsCache)[0];
......@@ -163,8 +168,8 @@ public:
curSol,
prevElemVolVars,
curElemVolVars,
prevGlobalFaceVars,
curGlobalFaceVars,
prevElemeFaceVars,
curElemeFaceVars,
elemFluxVarsCache,
elemBcTypes,
jac,
......@@ -261,8 +266,8 @@ protected:
const SolutionVector& curSol,
const ElementVolumeVariables& prevElemVolVars,
ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevGlobalFaceVars,
GlobalFaceVars& curGlobalFaceVars,
const ElementFaceVariables& prevElemFaceVars,
ElementFaceVariables& curElemFaceVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
......@@ -270,16 +275,16 @@ protected:
const FaceSolutionVector& faceResidualCache)
{
// compute the derivatives of the cell center residuals with respect to cell center dofs
dCCdCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
dCCdCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
// compute the derivatives of the cell center residuals with respect to face dofs
dCCdFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
dCCdFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
// compute the derivatives of the face residuals with respect to cell center dofs
dFacedCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
dFacedCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
// compute the derivatives of the face residuals with respect to face dofs
dFacedFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
dFacedFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
}
/*!
......@@ -292,8 +297,8 @@ static void dCCdCC_(Assembler& assembler,
const SolutionVector& curSol,
const ElementVolumeVariables& prevElemVolVars,
ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevGlobalFaceVars,
const GlobalFaceVars& curGlobalFaceVars,
const ElementFaceVariables& prevElemFaceVars,
const ElementFaceVariables& curElemFaceVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
......@@ -330,7 +335,7 @@ static void dCCdCC_(Assembler& assembler,
curVolVars.update(elemSol, problem, elementJ, scvJ);
const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry, prevElemVolVars, curElemVolVars,
prevGlobalFaceVars, curGlobalFaceVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
const auto partialDeriv = (deflectedResidual - ccResidual) / eps;
......@@ -354,8 +359,8 @@ static void dCCdFace_(Assembler& assembler,
const SolutionVector& curSol,
const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevGlobalFaceVars,
GlobalFaceVars& curGlobalFaceVars,
const ElementFaceVariables& prevElemFaceVars,
ElementFaceVariables& curElemFaceVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
......@@ -368,6 +373,7 @@ static void dCCdFace_(Assembler& assembler,
const auto& problem = assembler.problem();
auto& localResidual = assembler.localResidual();
auto& gridVariables = assembler.gridVariables();
// build derivatives with for cell center dofs w.r.t. cell center dofs
const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
......@@ -377,8 +383,8 @@ static void dCCdFace_(Assembler& assembler,
const auto globalJ = scvfJ.dofIndex();
// get the faceVars of the face with respect to which we are going to build the derivative
auto origFaceVars = curGlobalFaceVars.faceVars(scvfJ.index());
auto& curFaceVars = curGlobalFaceVars.faceVars(scvfJ.index());
auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvfJ);
const auto origFaceVars(faceVars);
for(auto pvIdx : PriVarIndices(faceIdx))
{
......@@ -391,11 +397,11 @@ static void dCCdFace_(Assembler& assembler,
priVars[pvIdx] += eps;
faceSolution[globalJ][pvIdx] += eps;
curFaceVars.update(scvfJ, faceSolution);
faceVars.update(scvfJ, faceSolution);
const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry,
prevElemVolVars, curElemVolVars,
prevGlobalFaceVars, curGlobalFaceVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
const auto partialDeriv = (deflectedResidual - ccResidual) / eps;
......@@ -404,7 +410,7 @@ static void dCCdFace_(Assembler& assembler,
updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], cellCenterGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
// restore the original faceVars
curFaceVars = origFaceVars;
faceVars = origFaceVars;
}
}
}
......@@ -419,8 +425,8 @@ static void dFacedCC_(Assembler& assembler,
const SolutionVector& curSol,
const ElementVolumeVariables& prevElemVolVars,
ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevGlobalFaceVars,
const GlobalFaceVars& curGlobalFaceVars,
const ElementFaceVariables& prevElemFaceVars,
const ElementFaceVariables& curElemFaceVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
......@@ -461,7 +467,7 @@ static void dFacedCC_(Assembler& assembler,
const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
prevElemVolVars, curElemVolVars,
prevGlobalFaceVars, curGlobalFaceVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
const auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]) / eps;
......@@ -486,8 +492,8 @@ static void dFacedFace_(Assembler& assembler,
const SolutionVector& curSol,
const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevGlobalFaceVars,
GlobalFaceVars& curGlobalFaceVars,
const ElementFaceVariables& prevElemFaceVars,
ElementFaceVariables& curElemFaceVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
......@@ -499,6 +505,7 @@ static void dFacedFace_(Assembler& assembler,
const auto& problem = assembler.problem();
auto& localResidual = assembler.localResidual();
const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
auto& gridVariables = assembler.gridVariables();
for(auto&& scvf : scvfs(fvGeometry))
{
......@@ -509,8 +516,8 @@ static void dFacedFace_(Assembler& assembler,
for(const auto& globalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
{
// get the faceVars of the face with respect to which we are going to build the derivative
auto origFaceVars = curGlobalFaceVars.faceVars(scvf.index());
auto& curFaceVars = curGlobalFaceVars.faceVars(scvf.index());
auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvf);
const auto origFaceVars(faceVars);
for(auto pvIdx : PriVarIndices(faceIdx))
{
......@@ -519,11 +526,11 @@ static void dFacedFace_(Assembler& assembler,
const Scalar eps = numericEpsilon(faceSolution[globalJ][pvIdx], faceIdx, faceIdx);
faceSolution[globalJ][pvIdx] += eps;
curFaceVars.update(scvf, faceSolution);
faceVars.update(scvf, faceSolution);
const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
prevElemVolVars, curElemVolVars,
prevGlobalFaceVars, curGlobalFaceVars,
prevElemFaceVars, curElemFaceVars,
elemBcTypes, elemFluxVarsCache);
const auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]) / eps;
......@@ -532,7 +539,7 @@ static void dFacedFace_(Assembler& assembler,
updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
// restore the original faceVars
curFaceVars = origFaceVars;
faceVars = origFaceVars;
}
}
}
......@@ -603,6 +610,16 @@ private:
static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
{ return gridVolVars.volVars(scv); }
template<class T = TypeTag>
static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), FaceVariables&>::type //TODO: introduce correct property
getFaceVarAccess(GlobalFaceVars& gridFaceVars, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
{ return elemFaceVars[scvf]; }
template<class T = TypeTag>
static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), FaceVariables&>::type
getFaceVarAccess(GlobalFaceVars& gridFaceVars, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
{ return gridFaceVars.faceVars(scvf.index()); }
};
} // end namespace Dumux
......
dumux/discretization/staggered/elementfacevariables.hh 0 → 100644
View file @ f4c85a46
// -*- 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 face variables class for free flow staggered grid models
*/
#ifndef DUMUX_DISCRETIZATION_STAGGERED_ELEMENTFACEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_ELEMENTFACEVARIABLES_HH
#include <dumux/common/basicproperties.hh>
namespace Dumux
{
template<class TypeTag>
class StaggeredElementFaceVariables
{
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
public:
StaggeredElementFaceVariables(const GlobalFaceVars& globalFacesVars) : globalFaceVarsPtr_(&globalFacesVars) {}
const FaceVariables& operator [](const SubControlVolumeFace& scvf) const
{ return globalFaceVars().faceVars(scvf.index()); }
// operator for the access with an index
// needed for cc methods for the access to the boundary volume variables
const FaceVariables& operator [](const IndexType scvfIdx) const
{ return globalFaceVars().faceVars(scvfIdx); }
//! 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)
{}
//! The global volume variables object we are a restriction of
const GlobalFaceVars& globalFaceVars() const
{ return *globalFaceVarsPtr_; }
private:
const GlobalFaceVars* globalFaceVarsPtr_;
};
} // end namespace
#endif
dumux/discretization/staggered/freeflow/facevariables.hh
View file @ f4c85a46
......@@ -40,6 +40,7 @@ class StaggeredFaceVariables
using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
using FaceSolution = typename GET_PROP_TYPE(TypeTag, StaggeredFaceSolution);
struct SubFaceData
......@@ -51,6 +52,19 @@ class StaggeredFaceVariables
};
public:
StaggeredFaceVariables() = default;
StaggeredFaceVariables(const GlobalFaceVars& globalFacesVars) : globalFaceVarsPtr_(&globalFacesVars) {}
const StaggeredFaceVariables& operator [](const SubControlVolumeFace& scvf) const
{ return globalFaceVars().faceVars(scvf.index()); }
// // operator for the access with an index
// // needed for cc methods for the access to the boundary volume variables
// const StaggeredFaceVariables& operator [](const IndexType scvIdx) const
// { return globalVolVars().volVars(scvIdx); }
void update(const FacePrimaryVariables &facePrivars)
{
velocitySelf_ = facePrivars;
......@@ -120,8 +134,13 @@ public:
return subFaceVelocities_[localSubFaceIdx];
}
//! The global volume variables object we are a restriction of
const GlobalFaceVars& globalFaceVars() const
{ return *globalFaceVarsPtr_; }
private:
const GlobalFaceVars* globalFaceVarsPtr_;
Scalar velocity_;
Scalar velocitySelf_;
Scalar velocityOpposite_;
......
dumux/discretization/staggered/globalfacevariables.hh
View file @ f4c85a46
......@@ -29,6 +29,11 @@
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(ElementFaceVariables);
}
template<class TypeTag>
class StaggeredGlobalFaceVariables
{
......@@ -37,6 +42,7 @@ class StaggeredGlobalFaceVariables
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
......@@ -49,14 +55,7 @@ public:
void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
{
const auto& faceSol = sol[faceIdx];
// const auto numFaceDofs = fvGridGeometry.gridView().size(1);
// faceVariables_.resize(numFaceDofs);
faceVariables_.resize(fvGridGeometry.numScvf());
for(auto&& element : elements(fvGridGeometry.gridView()))
......@@ -67,21 +66,8 @@ public:
for(auto&& scvf : scvfs(fvGeometry))
{
faceVariables_[scvf.index()].update(scvf, faceSol);
auto test = StaggeredFaceSolution<TypeTag>(scvf, faceSol, fvGridGeometry);
// std::cout << "test" << test[scvf.dofIndex()] << std::endl;
}
}
// for(auto& faceVariable : faceVariables_)
// {
// faceVariable.update()
// }
// // assert(faceVariables_.size() == faceSol.size());
//
// for(int i = 0; i < numFaceDofs; ++i)
// {
// faceVariables_[i].update(faceSol[i]);
// }
}
const FaceVariables& faceVars(const IndexType facetIdx) const
......@@ -89,6 +75,16 @@ public:
FaceVariables& faceVars(const IndexType facetIdx)
{ return faceVariables_[facetIdx]; }
/*!
* \brief Return a local restriction of this global object
* The local object is only functional after calling its bind/bindElement method
* This is a free function that will be found by means of ADL
*/
friend inline ElementFaceVariables localView(const StaggeredGlobalFaceVariables& global)
{ return ElementFaceVariables(global); }
private:
const Problem& problem_() const
{ return *problemPtr_; }
......
dumux/discretization/staggered/properties.hh
View file @ f4c85a46
......@@ -45,6 +45,7 @@
#include <dumux/discretization/staggered/fvelementgeometry.hh>
#include <dumux/discretization/staggered/globalfacevariables.hh>
#include <dumux/discretization/staggered/facesolution.hh>
#include <dumux/discretization/staggered/elementfacevariables.hh>
#include <dumux/common/intersectionmapper.hh>
#include <dune/istl/multitypeblockvector.hh>
......@@ -62,6 +63,7 @@ namespace Properties
NEW_PROP_TAG(CellCenterSolutionVector);
NEW_PROP_TAG(FaceSolutionVector);
NEW_PROP_TAG(StaggeredFaceSolution);
NEW_PROP_TAG(ElementFaceVariables);
//! Type tag for the box scheme.
NEW_TYPE_TAG(StaggeredModel, INHERITS_FROM(FiniteVolumeModel, NumericModel));
......@@ -113,6 +115,8 @@ SET_TYPE_PROP(StaggeredModel, IntersectionMapper, Dumux::ConformingGridIntersect
SET_TYPE_PROP(StaggeredModel, StaggeredFaceSolution, StaggeredFaceSolution<TypeTag>);
SET_TYPE_PROP(StaggeredModel, ElementFaceVariables, StaggeredElementFaceVariables<TypeTag>);
//! Definition of the indices for cell center and face dofs in the global solution vector
SET_PROP(StaggeredModel, DofTypeIndices)
{
......
dumux/freeflow/staggered/fluxvariables.hh
View file @ f4c85a46
......@@ -66,7 +66,6 @@ class FreeFlowFluxVariablesImpl<TypeTag, false>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
......@@ -74,6 +73,8 @@ class FreeFlowFluxVariablesImpl<TypeTag, false>
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
enum {
......@@ -100,7 +101,7 @@ public:
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf,
const FluxVariablesCache& fluxVarsCache)
{
......@@ -111,8 +112,7 @@ public:
const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
CellCenterPrimaryVariables flux(0.0);
// const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
const Scalar velocity = globalFaceVars.faceVars(scvf.index()).velocitySelf();
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
......@@ -191,19 +191,19 @@ public:
* \param scvf The sub control volume face
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param globalFaceVars The face variables
* \param elementFaceVars The face variables
*/
FacePrimaryVariables computeNormalMomentumFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars)
const ElementFaceVariables& elementFaceVars)
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
const Scalar velocitySelf = globalFaceVars.faceVars(scvf.index()).velocitySelf() ;
const Scalar velocityOpposite = globalFaceVars.faceVars(scvf.index()).velocityOpposite();
const Scalar velocitySelf = elementFaceVars[scvf].velocitySelf() ;
const Scalar velocityOpposite = elementFaceVars[scvf].velocityOpposite();
FacePrimaryVariables normalFlux(0.0);
if(navierStokes)
......@@ -243,28 +243,20 @@ public:
* \param scvf The sub control volume face
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param globalFaceVars The face variables
* \param elementFaceVars The face variables
*/
FacePrimaryVariables computeTangetialMomentumFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars)
const ElementFaceVariables& elementFaceVars)
{
FacePrimaryVariables tangentialFlux(0.0);
// // convenience function to get the velocity on a face
// auto velocity = [&globalFaceVars](const int dofIdx)
// {
// return globalFaceVars.faceVars(scvf.dofIdx()).velocity();
// };
auto& faceVars = globalFaceVars.faceVars(scvf.index());
const int numSubFaces = scvf.pairData().size();
auto& faceVars = elementFaceVars[scvf];
const int numSubFaces = scvf.pairData().size();
// account for all sub-faces
// for(const auto& subFaceData : scvf.pairData())
for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
......@@ -305,7 +297,6 @@ private:
const FaceVars& faceVars,
const int localSubFaceIdx)
{
// const Scalar transportingVelocity = velocity(subFaceData.normalPair.first);
const Scalar transportingVelocity = faceVars.subFaceData(localSubFaceIdx).velocityNormalInside;
const auto insideScvIdx = normalFace.insideScvIdx();
const auto outsideScvIdx = normalFace.outsideScvIdx();
......@@ -324,12 +315,10 @@ private:
if(innerElementIsUpstream)
transportedVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelInside;
// transportedVelocity = velocity(scvf.dofIndex());
else
{
const int outerDofIdx = scvf.pairData(localSubFaceIdx).outerParallelFaceDofIdx;
if(outerDofIdx >= 0)
// transportedVelocity = velocity(outerDofIdx);
transportedVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelOutside;
else // this is the case when the outer parallal dof would lie outside the domain TODO: discuss which one is better
// transportedVelocity = problem.dirichlet(makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
......@@ -370,10 +359,8 @@ private:
const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;
// the normal derivative
// const int innerNormalVelocityIdx = subFaceData.normalPair.first;
const int outerNormalVelocityIdx = scvf.pairData(localSubFaceIdx).normalPair.second;
// const Scalar innerNormalVelocity = velocity(innerNormalVelocityIdx);
const Scalar innerNormalVelocity = faceVars.subFaceData(localSubFaceIdx).velocityNormalInside;
const Scalar outerNormalVelocity = outerNormalVelocityIdx >= 0 ?
......
dumux/freeflow/staggered/localresidual.hh
View file @ f4c85a46
......@@ -35,6 +35,7 @@ NEW_PROP_TAG(EnableInertiaTerms);
NEW_PROP_TAG(ReplaceCompEqIdx);
NEW_PROP_TAG(EnergyFluxVariables);
NEW_PROP_TAG(NormalizePressure);
NEW_PROP_TAG(ElementFaceVariables);
}
/*!
......@@ -86,6 +87,7 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
using EnergyFluxVariables = typename GET_PROP_TYPE(TypeTag, EnergyFluxVariables);
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
......@@ -109,7 +111,6 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere