Commit bff3aa6c authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[staggered][freeflow] First step to improve free-flow structure

parent 62edaf04
// -*- 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/>. *
*****************************************************************************/
/*!
* \ingroup Properties
* \file
*
* \brief Defines a type tag and some properties for free-flow models using the staggered scheme.
*/
#ifndef DUMUX_STAGGERD_FREE_FLOW_PROPERTIES_HH
#define DUMUX_STAGGERD_FREE_FLOW_PROPERTIES_HH
#include <dumux/common/properties.hh>
#include <dumux/discretization/staggered/properties.hh>
#include <dumux/freeflow/navierstokes/staggered/boundarytypes.hh>
#include "subcontrolvolumeface.hh"
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(NumEq);
//! Type tag for the staggered scheme.
NEW_TYPE_TAG(StaggeredFreeFlowModel, INHERITS_FROM(StaggeredModel));
// UNSET_PROP(FreeFlow, PrimaryVariables);
UNSET_PROP(FreeFlow, NumEqVector);
SET_INT_PROP(StaggeredFreeFlowModel, NumEqFace, 1); //!< set the number of equations to 1
SET_PROP(StaggeredFreeFlowModel, NumEqCellCenter)
{
private:
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
static constexpr auto dim = GridView::dimension;
static constexpr auto numEq = GET_PROP_VALUE(TypeTag, NumEq);
public:
// static_assert(true, "fuuu");
// static_assert(dim == 2, "fuuuDim");
// static constexpr int value = GET_PROP_VALUE(TypeTag, NumEq) - dim;
static constexpr int value = 1;
// static_assert(value == 1, "fuuuVal");
};
//! The default sub-controlvolume face
SET_PROP(StaggeredFreeFlowModel, SubControlVolumeFace)
{
private:
using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
static constexpr int dim = Grid::dimension;
static constexpr int dimWorld = Grid::dimensionworld;
struct ScvfGeometryTraits
{
using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
using LocalIndexType = unsigned int;
using Scalar = typename Grid::ctype;
using Geometry = typename Grid::template Codim<1>::Geometry;
using GlobalPosition = Dune::FieldVector<Scalar, dim>;
};
public:
using type = FreeFlowStaggeredSubControlVolumeFace<ScvfGeometryTraits>;
};
//! The default geometry helper required for the stencils, etc.
SET_PROP(StaggeredFreeFlowModel, StaggeredGeometryHelper)
{
private:
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
public:
using type = FreeFlowStaggeredGeometryHelper<GridView>;
};
//! The variables living on the faces
SET_TYPE_PROP(StaggeredFreeFlowModel, FaceVariables, StaggeredFaceVariables<TypeTag>);
//! A container class used to specify values for boundary/initial conditions
SET_PROP(StaggeredFreeFlowModel, PrimaryVariables)
{
private:
using CellCenterBoundaryValues = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FaceBoundaryValues = Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
GridView::dimension>;
public:
using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
};
//! A container class used to specify values for sources and Neumann BCs
SET_PROP(StaggeredFreeFlowModel, NumEqVector)
{
private:
using CellCenterBoundaryValues = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FaceBoundaryValues = Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
GridView::dimension>;
public:
using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
};
//! Boundary types at a single degree of freedom
SET_PROP(StaggeredFreeFlowModel, BoundaryTypes)
{
private:
static constexpr auto size = GET_PROP_VALUE(TypeTag, NumEqCellCenter) + GET_PROP_VALUE(TypeTag, NumEqFace);
public:
using type = StaggeredFreeFlowBoundaryTypes<size>;
};
} // namespace Properties
} // namespace Dumux
#endif
add_subdirectory("staggered")
add_subdirectory("navierstokes")
add_subdirectory("staggerednc")
add_subdirectory("staggeredni")
// -*- 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 the flux variables
*/
#ifndef DUMUX_FREELOW_IMPLICIT_FLUXVARIABLES_HH
#define DUMUX_FREELOW_IMPLICIT_FLUXVARIABLES_HH
#include <dumux/common/properties.hh>
#include <dumux/discretization/fluxvariablesbase.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/freeflow/navierstokes/staggered/fluxvariables.hh>
namespace Dumux
{
namespace Properties
{
// forward declaration
NEW_PROP_TAG(EnableInertiaTerms);
}
// forward declaration
template<class TypeTag, DiscretizationMethods Method>
class NavierStokesFluxVariablesImpl;
/*!
* \ingroup ImplicitModel
* \brief The flux variables class
* specializations are provided for combinations of physical processes
* \note Not all specializations are currently implemented
*/
template<class TypeTag>
using FreeFlowFluxVariables = NavierStokesFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, DiscretizationMethod)>;
} // end namespace
#endif
......@@ -79,7 +79,7 @@ NEW_PROP_TAG(EnergyFluxVariables); //!< The energy flux variables
///////////////////////////////////////////////////////////////////////////
// default property values for the isothermal single phase model
///////////////////////////////////////////////////////////////////////////
SET_INT_PROP(NavierStokes, NumEqCellCenter, 1); //! set the number of equations to 1
// SET_INT_PROP(NavierStokes, NumEqCellCenter, 1); //! set the number of equations to 1
SET_INT_PROP(NavierStokes, NumPhases, 1); //! The number of phases in the 1p model is 1
SET_INT_PROP(NavierStokes, NumComponents, 1); //! The number of components in the 1p model is 1
......
// -*- 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 the flux variables
*/
#ifndef DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
#define DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
#include <dumux/common/properties.hh>
#include <dumux/discretization/fluxvariablesbase.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/freeflow/navierstokes/fluxvariables.hh>
namespace Dumux
{
namespace Properties
{
// forward declaration
NEW_PROP_TAG(EnableInertiaTerms);
}
// forward declaration
template<class TypeTag, DiscretizationMethods Method>
class NavierStokesFluxVariablesImpl;
/*!
* \ingroup Discretization
* \brief Base class for the flux variables
* Actual flux variables inherit from this class
*/
// specialization for immiscible, isothermal flow
template<class TypeTag>
class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
: public FluxVariablesBase<TypeTag>
{
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
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 SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
enum {
// grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
pressureIdx = Indices::pressureIdx,
velocityIdx = Indices::velocityIdx,
massBalanceIdx = Indices::massBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx
};
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
public:
CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf,
const FluxVariablesCache& fluxVarsCache)
{
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
// if we are on an inflow/outflow boundary, use the volVars of the element itself
const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
CellCenterPrimaryVariables flux(0.0);
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
static const Scalar upWindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
flux = (upWindWeight * upstreamVolVars.density() +
(1.0 - upWindWeight) * downstreamVolVars.density()) * velocity;
flux *= scvf.area() * sign(scvf.outerNormalScalar());
return flux;
}
void computeCellCenterToCellCenterStencil(Stencil& stencil,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
// the first entry is always the cc dofIdx itself
if(stencil.empty())
stencil.push_back(scvf.insideScvIdx());
if(!scvf.boundary())
stencil.push_back(scvf.outsideScvIdx());
}
void computeCellCenterToFaceStencil(Stencil& stencil,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
stencil.push_back(scvf.dofIndex());
}
void computeFaceToCellCenterStencil(Stencil& stencil,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
const int eIdx = scvf.insideScvIdx();
stencil.push_back(scvf.insideScvIdx());
for(const auto& data : scvf.pairData())
{
auto& normalFace = fvGeometry.scvf(eIdx, data.localNormalFaceIdx);
const auto outerParallelElementDofIdx = normalFace.outsideScvIdx();
if(!normalFace.boundary())
stencil.push_back(outerParallelElementDofIdx);
}
}
void computeFaceToFaceStencil(Stencil& stencil,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
// the first entries are always the face dofIdx itself and the one of the opposing face
if(stencil.empty())
{
stencil.push_back(scvf.dofIndex());
stencil.push_back(scvf.dofIndexOpposingFace());
}
for(const auto& data : scvf.pairData())
{
stencil.push_back(data.normalPair.first);
const auto outerParallelFaceDofIdx = data.outerParallelFaceDofIdx;
if(outerParallelFaceDofIdx >= 0)
stencil.push_back(outerParallelFaceDofIdx);
if(!scvf.boundary())
stencil.push_back(data.normalPair.second);
}
}
/*!
* \brief Returns the normal part of the momentum flux
* \param scvf The sub control volume face
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param elementFaceVars The face variables
*/
FacePrimaryVariables computeNormalMomentumFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elementFaceVars)
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
const Scalar velocitySelf = elementFaceVars[scvf].velocitySelf() ;
const Scalar velocityOpposite = elementFaceVars[scvf].velocityOpposite();
FacePrimaryVariables normalFlux(0.0);
if(navierStokes)
{
// advective part
const Scalar vAvg = (velocitySelf + velocityOpposite) * 0.5;
const Scalar vUp = (sign(scvf.outerNormalScalar()) == sign(vAvg)) ? velocityOpposite : velocitySelf;
normalFlux += vAvg * vUp * insideVolVars.density();
}
// diffusive part
const Scalar deltaV = scvf.normalInPosCoordDir() ?
(velocitySelf - velocityOpposite) :
(velocityOpposite - velocitySelf);
const Scalar deltaX = scvf.selfToOppositeDistance();
normalFlux -= insideVolVars.viscosity() * 2.0 * deltaV/deltaX;
// account for the orientation of the face
const Scalar sgn = -1.0 * sign(scvf.outerNormalScalar());
Scalar result = normalFlux * sgn * scvf.area();
// treat outflow conditions
if(navierStokes && scvf.boundary())
{
const auto& upVolVars = (sign(scvf.outerNormalScalar()) == sign(velocitySelf)) ?
elemVolVars[insideScvIdx] : elemVolVars[scvf.outsideScvIdx()] ;
result += velocitySelf * velocitySelf * upVolVars.density() * sign(scvf.outerNormalScalar()) * scvf.area() ;
}
return result;
}
/*!
* \brief Returns the tangential part of the momentum flux
* \param scvf The sub control volume face
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param elementFaceVars The face variables
*/
FacePrimaryVariables computeTangetialMomentumFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elementFaceVars)
{
FacePrimaryVariables tangentialFlux(0.0);
auto& faceVars = elementFaceVars[scvf];
const int numSubFaces = scvf.pairData().size();
// account for all sub-faces
for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
// Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
if(scvf.pairData()[localSubFaceIdx].outerParallelFaceDofIdx < 0)
{
// lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
auto makeGhostFace = [eIdx] (const GlobalPosition& pos)
{
return SubControlVolumeFace(pos, std::vector<unsigned int>{eIdx,eIdx});
};
// use the ghost face to check if there is a symmetry boundary condition and skip any further steps if yes
const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(scvf.pairData()[localSubFaceIdx].virtualOuterParallelFaceDofPos));
if(bcTypes.isSymmetry())
continue;
}
// if there is no symmetry boundary condition, proceed to calculate the tangential momentum flux
if(navierStokes)
tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
tangentialFlux += computeDiffusivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
}
return tangentialFlux;
}
private:
FacePrimaryVariables computeAdvectivePartOfTangentialMomentumFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
const ElementVolumeVariables& elemVolVars,
const FaceVariables& faceVars,
const int localSubFaceIdx)
{
const Scalar transportingVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
const auto insideScvIdx = normalFace.insideScvIdx();
const auto outsideScvIdx = normalFace.outsideScvIdx();
const bool innerElementIsUpstream = ( sign(normalFace.outerNormalScalar()) == sign(transportingVelocity) );
const auto& upVolVars = innerElementIsUpstream ? elemVolVars[insideScvIdx] : elemVolVars[outsideScvIdx];
const Scalar transportedVelocity = innerElementIsUpstream ?
faceVars.velocitySelf() :
faceVars.velocityParallel(localSubFaceIdx);
const Scalar momentum = upVolVars.density() * transportedVelocity;
const int sgn = sign(normalFace.outerNormalScalar());
return transportingVelocity * momentum * sgn * normalFace.area() * 0.5;
}
FacePrimaryVariables computeDiffusivePartOfTangentialMomentumFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
const ElementVolumeVariables& elemVolVars,
const FaceVariables& faceVars,
const int localSubFaceIdx)
{
FacePrimaryVariables tangentialDiffusiveFlux(0.0);
const auto insideScvIdx = normalFace.insideScvIdx();
const auto outsideScvIdx = normalFace.outsideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
// the averaged viscosity at the face normal to our face of interest (where we assemble the face residual)
const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;
// the normal derivative
const Scalar innerNormalVelocity = faceVars.velocityNormalInside(localSubFaceIdx);
const Scalar outerNormalVelocity = faceVars.velocityNormalOutside(localSubFaceIdx);
const Scalar normalDeltaV = scvf.normalInPosCoordDir() ?
(outerNormalVelocity - innerNormalVelocity) :
(innerNormalVelocity - outerNormalVelocity);
const Scalar normalDerivative = normalDeltaV / scvf.pairData(localSubFaceIdx).normalDistance;
tangentialDiffusiveFlux -= muAvg * normalDerivative;
// the parallel derivative
const Scalar innerParallelVelocity = faceVars.velocitySelf();
const Scalar outerParallelVelocity = faceVars.velocityParallel(localSubFaceIdx);
const Scalar parallelDeltaV = normalFace.normalInPosCoordDir() ?
(outerParallelVelocity - innerParallelVelocity) :
(innerParallelVelocity - outerParallelVelocity);
const Scalar parallelDerivative = parallelDeltaV / scvf.pairData(localSubFaceIdx).parallelDistance;
tangentialDiffusiveFlux -= muAvg * parallelDerivative;
const Scalar sgn = sign(normalFace.outerNormalScalar());
return tangentialDiffusiveFlux * sgn * normalFace.area() * 0.5;
}
};
} // end namespace
#endif