From e4c9d81e6546a2ef858c6b0f9d0f05318ac3245c Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 7 Sep 2021 13:09:59 +0200
Subject: [PATCH] [ff][staggered][navierstokes] Add helper classes for new
staggered models
---
.../navierstokes/momentum/fluxhelper.hh | 29 ++-
.../freeflow/navierstokes/scalarfluxhelper.hh | 183 ++++++++++++++++++
.../navierstokes/scalarfluxvariables.hh | 157 +++++++++++++++
.../scalarfluxvariablescachefiller.hh | 163 ++++++++++++++++
.../navierstokes/scalarvolumevariables.hh | 106 ++++++++++
dumux/freeflow/navierstokes/velocityoutput.hh | 103 ++++++++++
6 files changed, 723 insertions(+), 18 deletions(-)
create mode 100644 dumux/freeflow/navierstokes/scalarfluxhelper.hh
create mode 100644 dumux/freeflow/navierstokes/scalarfluxvariables.hh
create mode 100644 dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh
create mode 100644 dumux/freeflow/navierstokes/scalarvolumevariables.hh
create mode 100644 dumux/freeflow/navierstokes/velocityoutput.hh
diff --git a/dumux/freeflow/navierstokes/momentum/fluxhelper.hh b/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
index 34fdb81355..46a95d283a 100644
--- a/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
+++ b/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
@@ -88,6 +88,10 @@ struct NavierStokesMomentumBoundaryFluxHelper
if (scvf.isLateral())
{
+ // if the lateral scvf is adjacent to a slip boundary, call the helper function for this special case
+ if (scv.boundary() && problem.onSlipBoundary(fvGeometry, fvGeometry.frontalScvfOnBoundary(scv)))
+ return slipVelocityMomentumFlux(problem, fvGeometry, scvf, elemVolVars, elemFluxVarsCache);
+
// viscous terms
const Scalar mu = problem.effectiveViscosity(element, fvGeometry, scvf);
@@ -101,7 +105,7 @@ struct NavierStokesMomentumBoundaryFluxHelper
* scvf.directionSign();
// advective terms
- if (problem.enableInertiaTerms())
+ if (problem.enableInertiaTerms()) // TODO revise advection terms! Take care with upwinding!
{
const auto transportingVelocity = [&]()
{
@@ -164,11 +168,9 @@ struct NavierStokesMomentumBoundaryFluxHelper
// lateral face coinciding with boundary
else if (scvf.boundary())
{
- assert(false);
- // this should not happen? TODO revise
- // const auto insideDensity = problem.density(element, fvGeometry.scv(scvf.insideScvIdx()));
- // const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity();
- // flux[scv.dofAxis()] += innerVelocity * transportingVelocity * insideDensity * scvf.directionSign();
+ const auto insideDensity = problem.density(element, fvGeometry.scv(scvf.insideScvIdx()));
+ const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity();
+ flux[scv.dofAxis()] += innerVelocity * transportingVelocity * insideDensity * scvf.directionSign();
}
}
}
@@ -273,19 +275,10 @@ struct NavierStokesMomentumBoundaryFluxHelper
}
}();
- const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity();
- const auto outerVelocity = slipVelocity;
-
+ // Do not use upwinding here but directly take the slip velocity located on the boundary. Upwinding with a weight of 0.5
+ // would actually prevent second order grid convergence.
const auto rho = problem.getInsideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf);
- const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity);
-
- const auto insideMomentum = innerVelocity * rho.first;
- const auto outsideMomentum = outerVelocity * rho.second;
-
- static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
-
- const auto transportedMomentum = selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum)
- : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum);
+ const auto transportedMomentum = slipVelocity * rho.second;
flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign();
}
diff --git a/dumux/freeflow/navierstokes/scalarfluxhelper.hh b/dumux/freeflow/navierstokes/scalarfluxhelper.hh
new file mode 100644
index 0000000000..9e35c8b5f5
--- /dev/null
+++ b/dumux/freeflow/navierstokes/scalarfluxhelper.hh
@@ -0,0 +1,183 @@
+// -*- 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 3 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/>. *
+ *****************************************************************************/
+
+#ifndef DUMUX_NAVIERSTOKES_SCALAR_BOUNDARY_FLUXHELPER_HH
+#define DUMUX_NAVIERSTOKES_SCALAR_BOUNDARY_FLUXHELPER_HH
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/std/type_traits.hh>
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/discretization/elementsolution.hh>
+
+namespace Dumux {
+
+#ifndef DOXYGEN
+namespace Detail {
+// helper structs and functions detecting if the VolumeVariables belong to a non-isothermal model
+template <class Indices>
+using NonisothermalDetector = decltype(std::declval<Indices>().energyEqIdx);
+
+template<class Indices>
+static constexpr bool isNonIsothermal()
+{ return Dune::Std::is_detected<NonisothermalDetector, Indices>::value; }
+
+} // end namespace Detail
+#endif
+
+template<class AdvectiveFlux>
+struct NavierStokesScalarBoundaryFluxHelper
+{
+ /*!
+ * \brief Return the area-specific, weighted advective flux of a scalar quantity.
+ */
+ template<class VolumeVariables, class SubControlVolumeFace, class Scalar, class UpwindTerm>
+ static Scalar advectiveScalarUpwindFlux(const VolumeVariables& insideVolVars,
+ const VolumeVariables& outsideVolVars,
+ const SubControlVolumeFace& scvf,
+ const Scalar volumeFlux,
+ const Scalar upwindWeight,
+ UpwindTerm upwindTerm)
+ {
+ using std::signbit;
+ const bool insideIsUpstream = !signbit(volumeFlux);
+
+ const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+ const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars;
+
+ return (upwindWeight * upwindTerm(upstreamVolVars) +
+ (1.0 - upwindWeight) * upwindTerm(downstreamVolVars))
+ * volumeFlux;
+ }
+
+ template<class Indices, class NumEqVector, class UpwindFunction>
+ static void addModelSpecificAdvectiveFlux(NumEqVector& flux,
+ const UpwindFunction& upwind)
+ {
+ // add advective fluxes based on physical type of model
+ AdvectiveFlux::addAdvectiveFlux(flux, upwind);
+
+ // for non-isothermal models, add the energy flux
+ if constexpr (Detail::isNonIsothermal<Indices>())
+ {
+ auto upwindTerm = [](const auto& volVars) { return volVars.density()*volVars.enthalpy(); };
+ flux[Indices::energyEqIdx] = upwind(upwindTerm);
+ }
+ }
+
+ /*!
+ * \brief Return the area-specific outflow fluxes for all scalar balance equations.
+ * The values specified in outsideBoundaryPriVars are used in case of flow reversal.
+ */
+ template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables>
+ static auto scalarOutflowFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars,
+ typename ElementVolumeVariables::VolumeVariables::PrimaryVariables&& outsideBoundaryPriVars,
+ const typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type upwindWeight = 1.0)
+ {
+ using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
+ using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
+ using NumEqVector = PrimaryVariables;
+ NumEqVector flux;
+ const auto velocity = problem.faceVelocity(element,fvGeometry, scvf);
+ const auto volumeFlux = velocity * scvf.unitOuterNormal();
+ using std::signbit;
+ const bool insideIsUpstream = !signbit(volumeFlux);
+ const VolumeVariables& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const VolumeVariables& outsideVolVars = [&]()
+ {
+ // only use the inside volVars for "true" outflow conditions (avoid constructing the outside volVars)
+ if (insideIsUpstream && Dune::FloatCmp::eq(upwindWeight, 1.0, 1e-6))
+ return insideVolVars;
+ else
+ {
+ // construct outside volVars from the given priVars for situations of flow reversal
+ VolumeVariables boundaryVolVars;
+ boundaryVolVars.update(elementSolution<FVElementGeometry>(std::forward<PrimaryVariables>(outsideBoundaryPriVars)),
+ problem,
+ element,
+ fvGeometry.scv(scvf.insideScvIdx()));
+ return boundaryVolVars;
+ }
+ }();
+
+ auto upwindFuntion = [&](const auto& upwindTerm)
+ {
+ return advectiveScalarUpwindFlux(insideVolVars, outsideVolVars, scvf, volumeFlux, upwindWeight, upwindTerm);
+ };
+
+ addModelSpecificAdvectiveFlux<typename VolumeVariables::Indices>(flux, upwindFuntion);
+
+ return flux;
+ }
+
+ /*!
+ * \brief Return the area-specific outflow fluxes for all scalar balance equations.
+ * This should only be used of flow reversal does never occur.
+ * A (deactivable) warning is emitted otherwise.
+ */
+ template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables>
+ static auto scalarOutflowFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars)
+ {
+ using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
+ using NumEqVector = typename VolumeVariables::PrimaryVariables;
+ NumEqVector flux;
+ const auto velocity = problem.faceVelocity(element,fvGeometry, scvf);
+ const auto volumeFlux = velocity * scvf.unitOuterNormal();
+ using std::signbit;
+ const bool insideIsUpstream = !signbit(volumeFlux);
+ const VolumeVariables& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+
+ if constexpr (VolumeVariables::FluidSystem::isCompressible(0/*phaseIdx*/) /*TODO viscosityIsConstant*/ || NumEqVector::size() > 1)
+ {
+ static const bool verbose = getParamFromGroup<bool>(problem.paramGroup(), "Flux.EnableOutflowReversalWarning", true);
+ using std::abs;
+ if (verbose && !insideIsUpstream && abs(volumeFlux) > 1e-10)
+ {
+ std::cout << "velo " << velocity << ", flux " << volumeFlux << std::endl;
+ std::cout << "\n ********** WARNING ********** \n\n"
+ "Outflow condition set at " << scvf.center() << " might be invalid due to flow reversal. "
+ "Consider using \n"
+ "outflowFlux(problem, element, fvGeometry, scvf, elemVolVars, outsideBoundaryPriVars, upwindWeight) \n"
+ "instead where you can specify primary variables for inflow situations.\n"
+ "\n ***************************** \n" << std::endl;
+ }
+ }
+
+ auto upwindFuntion = [&](const auto& upwindTerm)
+ {
+ return advectiveScalarUpwindFlux(insideVolVars, insideVolVars, scvf, volumeFlux, 1.0 /*upwindWeight*/, upwindTerm);
+ };
+
+ addModelSpecificAdvectiveFlux<typename VolumeVariables::Indices>(flux, upwindFuntion);
+
+ return flux;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/scalarfluxvariables.hh b/dumux/freeflow/navierstokes/scalarfluxvariables.hh
new file mode 100644
index 0000000000..3a07b5e9bb
--- /dev/null
+++ b/dumux/freeflow/navierstokes/scalarfluxvariables.hh
@@ -0,0 +1,157 @@
+// -*- 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 3 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
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesFluxVariablesImpl
+ */
+#ifndef DUMUX_NAVIERSTOKES_SCALAR_CONSERVATION_MODEL_FLUXVARIABLES_HH
+#define DUMUX_NAVIERSTOKES_SCALAR_CONSERVATION_MODEL_FLUXVARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/typetraits/problem.hh>
+#include <dumux/flux/fluxvariablesbase.hh>
+#include <dumux/discretization/extrusion.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/flux/upwindscheme.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief The flux variables base class for scalar quantities balanced in the Navier-Stokes model.
+ */
+template<class Problem,
+ class ModelTraits,
+ class FluxTypes,
+ class ElementVolumeVariables,
+ class ElementFluxVariablesCache,
+ class UpwindScheme = UpwindScheme<typename ProblemTraits<Problem>::GridGeometry>>
+class NavierStokesScalarConservationModelFluxVariables
+: public FluxVariablesBase<Problem,
+ typename ProblemTraits<Problem>::GridGeometry::LocalView,
+ ElementVolumeVariables,
+ ElementFluxVariablesCache>
+{
+ using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
+ using GridGeometry = typename ProblemTraits<Problem>::GridGeometry;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using Extrusion = Extrusion_t<typename ProblemTraits<Problem>::GridGeometry>;
+
+public:
+
+ struct AdvectionType
+ {
+ static Scalar flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+ const auto velocity = problem.faceVelocity(element, fvGeometry, scvf);
+ const Scalar volumeFlux = velocity*scvf.unitOuterNormal()*Extrusion::area(scvf)*extrusionFactor_(elemVolVars, scvf);
+ return volumeFlux;
+ }
+ };
+
+ /*!
+ * \brief Returns the advective flux computed by the respective law.
+ */
+ template<typename FunctionType>
+ Scalar getAdvectiveFlux(const FunctionType& upwindTerm) const
+ {
+ if constexpr (ModelTraits::enableAdvection())
+ {
+ const auto& scvf = this->scvFace();
+ const auto velocity = this->problem().faceVelocity(this->element(), this->fvGeometry(), scvf);
+ const Scalar volumeFlux = velocity*scvf.unitOuterNormal()*Extrusion::area(scvf)*extrusionFactor_(this->elemVolVars(), scvf);
+ return UpwindScheme::apply(*this, upwindTerm, volumeFlux, 0/*phaseIdx*/);
+ }
+ else
+ return 0.0;
+ }
+
+ /*!
+ * \brief Returns the conductive energy flux computed by the respective law.
+ */
+ Scalar heatConductionFlux() const
+ {
+ if constexpr (ModelTraits::enableEnergyBalance())
+ {
+ using HeatConductionType = typename FluxTypes::HeatConductionType;
+ return HeatConductionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->elemVolVars(),
+ this->scvFace(),
+ this->elemFluxVarsCache());
+ }
+ else
+ return 0.0;
+ }
+
+ /*!
+ * \brief Returns the advective energy flux.
+ */
+ Scalar heatAdvectionFlux() const
+ {
+ if constexpr (ModelTraits::enableEnergyBalance())
+ {
+ const auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); };
+ return getAdvectiveFlux(upwindTerm);
+ }
+ else
+ return 0.0;
+ }
+
+ /*!
+ * \brief Returns the total energy flux.
+ */
+ Scalar heatFlux() const
+ {
+ return heatConductionFlux() + heatAdvectionFlux();
+ }
+
+ /*!
+ * \brief Adds the energy flux to a given flux vector.
+ */
+ template<class NumEqVector>
+ void addHeatFlux(NumEqVector& flux) const
+ {
+ if constexpr (ModelTraits::enableEnergyBalance())
+ flux[ModelTraits::Indices::energyEqIdx] = heatFlux();
+ }
+
+private:
+
+ static Scalar extrusionFactor_(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf)
+ {
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+ return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh b/dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh
new file mode 100644
index 0000000000..d764a3d139
--- /dev/null
+++ b/dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh
@@ -0,0 +1,163 @@
+// -*- 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 3 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
+ * \ingroup FreeflowModels
+ * \brief A helper class to fill the flux variables cache
+ */
+#ifndef DUMUX_FREEFLOW_SCALAR_FLUXVARIABLESCACHE_FILLER_HH
+#define DUMUX_FREEFLOW_SCALAR_FLUXVARIABLESCACHE_FILLER_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/discretization/method.hh>
+#include <dumux/discretization/extrusion.hh>
+#include <dumux/flux/referencesystemformulation.hh>
+#include <dumux/common/typetraits/problem.hh>
+
+namespace Dumux {
+
+// forward declaration
+template<class Problem, class ModelTraits, bool diffusionIsSolDependent, bool heatConductionIsSolDependent, DiscretizationMethod discMethod>
+class FreeFlowScalarFluxVariablesCacheFillerImplementation;
+
+/*!
+ * \ingroup PorousmediumflowModels
+ * \brief The flux variables cache filler class for free flow
+ *
+ * Helps filling the flux variables cache depending several policies
+ */
+template<class Problem, class ModelTraits, bool diffusionIsSolDependent, bool heatConductionIsSolDependent>
+using FreeFlowScalarFluxVariablesCacheFiller = FreeFlowScalarFluxVariablesCacheFillerImplementation<Problem, ModelTraits, diffusionIsSolDependent, heatConductionIsSolDependent, ProblemTraits<Problem>::GridGeometry::discMethod>;
+
+//! Specialization of the flux variables cache filler for the cell centered tpfa method
+template<class Problem, class ModelTraits, bool diffusionIsSolDependent, bool heatConductionIsSolDependent>
+class FreeFlowScalarFluxVariablesCacheFillerImplementation<Problem, ModelTraits, diffusionIsSolDependent, heatConductionIsSolDependent, DiscretizationMethod::cctpfa>
+{
+ using GridGeometry = typename ProblemTraits<Problem>::GridGeometry;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ static constexpr bool diffusionEnabled = ModelTraits::enableMolecularDiffusion();
+ static constexpr bool heatConductionEnabled = ModelTraits::enableEnergyBalance();
+
+public:
+ static constexpr bool isSolDependent = (diffusionEnabled && diffusionIsSolDependent) ||
+ (heatConductionEnabled && heatConductionIsSolDependent);
+
+ //! The constructor. Sets the problem pointer
+ FreeFlowScalarFluxVariablesCacheFillerImplementation(const Problem& problem)
+ : problemPtr_(&problem) {}
+
+ /*!
+ * \brief function to fill the flux variables caches
+ *
+ * \param fluxVarsCacheContainer Either the element or global flux variables cache
+ * \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The corresponding sub-control volume face
+ * \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones)
+ */
+ template<class FluxVariablesCacheContainer, class FluxVariablesCache, class ElementVolumeVariables>
+ void fill(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const bool forceUpdateAll = false)
+ {
+ // fill the physics-related quantities of the caches
+ if (forceUpdateAll)
+ {
+ if constexpr (diffusionEnabled)
+ fillDiffusion_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ if constexpr (heatConductionEnabled)
+ fillHeatConduction_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ }
+ else
+ {
+ if constexpr (diffusionEnabled && diffusionIsSolDependent)
+ fillDiffusion_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ if constexpr (heatConductionEnabled && heatConductionIsSolDependent)
+ fillHeatConduction_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ }
+ }
+
+private:
+
+ const Problem& problem() const
+ { return *problemPtr_; }
+
+
+ //! method to fill the diffusive quantities
+ template<class FluxVariablesCache, class ElementVolumeVariables>
+ void fillDiffusion_(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using DiffusionType = typename ElementVolumeVariables::VolumeVariables::MolecularDiffusionType;
+ using DiffusionFiller = typename DiffusionType::Cache::Filler;
+ using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem;
+
+ static constexpr int numPhases = ModelTraits::numFluidPhases();
+ static constexpr int numComponents = ModelTraits::numFluidComponents();
+
+ // forward to the filler of the diffusive quantities
+ if constexpr (FluidSystem::isTracerFluidSystem())
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
+ DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+ else
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
+ if (compIdx != FluidSystem::getMainComponent(phaseIdx))
+ DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+ }
+
+ //! method to fill the quantities related to heat conduction
+ template<class FluxVariablesCache, class ElementVolumeVariables>
+ void fillHeatConduction_(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using HeatConductionType = typename ElementVolumeVariables::VolumeVariables::HeatConductionType;
+ using HeatConductionFiller = typename HeatConductionType::Cache::Filler;
+
+ // forward to the filler of the diffusive quantities
+ HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+ }
+
+ const Problem* problemPtr_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/navierstokes/scalarvolumevariables.hh b/dumux/freeflow/navierstokes/scalarvolumevariables.hh
new file mode 100644
index 0000000000..ebbaa231f9
--- /dev/null
+++ b/dumux/freeflow/navierstokes/scalarvolumevariables.hh
@@ -0,0 +1,106 @@
+// -*- 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 3 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
+ * \ingroup NavierStokesModel
+ * \copydoc Dumux::NavierStokesVolumeVariables
+ */
+#ifndef DUMUX_NAVIERSTOKES_SCALAR_CONSERVATION_MODEL_VOLUME_VARIABLES_HH
+#define DUMUX_NAVIERSTOKES_SCALAR_CONSERVATION_MODEL_VOLUME_VARIABLES_HH
+
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Volume variables for the single-phase Navier-Stokes model.
+ */
+template <class Traits>
+class NavierStokesScalarConservationModelVolumeVariables
+{
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+
+public:
+ //! export the type used for the primary variables
+ using PrimaryVariables = typename Traits::PrimaryVariables;
+ //! export the indices type
+ using Indices = typename Traits::ModelTraits::Indices;
+ //! Export the underlying fluid system
+ using FluidSystem = typename Traits::FluidSystem;
+ //! Export the fluid state type
+ using FluidState = typename Traits::FluidState;
+
+ //! Return number of phases considered by the model
+ static constexpr int numFluidPhases() { return Traits::ModelTraits::numFluidPhases(); }
+ //! Return number of components considered by the model
+ static constexpr int numFluidComponents() { return Traits::ModelTraits::numFluidComponents(); }
+
+ /*!
+ * \brief Update all quantities for a given control volume
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to
+ * be simulated
+ * \param element An element which contains part of the control volume
+ * \param scv The sub-control volume
+ */
+ template<class ElementSolution, class Problem, class Element, class SubControlVolume>
+ void update(const ElementSolution& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ priVars_ = elemSol[scv.indexInElement()];
+ extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol);
+ }
+
+ /*!
+ * \brief Return how much the sub-control volume is extruded.
+ *
+ * This means the factor by which a lower-dimensional (1D or 2D)
+ * entity needs to be expanded to get a full dimensional cell. The
+ * default is 1.0 which means that 1D problems are actually
+ * thought as pipes with a cross section of 1 m^2 and 2D problems
+ * are assumed to extend 1 m to the back.
+ */
+ Scalar extrusionFactor() const
+ { return extrusionFactor_; }
+
+ /*!
+ * \brief Return a component of primary variable vector
+ *
+ * \param pvIdx The index of the primary variable of interest
+ */
+ Scalar priVar(const int pvIdx) const
+ { return priVars_[pvIdx]; }
+
+ /*!
+ * \brief Return the primary variable vector
+ */
+ const PrimaryVariables& priVars() const
+ { return priVars_; }
+
+protected:
+ PrimaryVariables priVars_;
+ Scalar extrusionFactor_;
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_NAVIERSTOKES_MOMENTUM_VOLUME_VARIABLES_HH
diff --git a/dumux/freeflow/navierstokes/velocityoutput.hh b/dumux/freeflow/navierstokes/velocityoutput.hh
new file mode 100644
index 0000000000..eb504effe3
--- /dev/null
+++ b/dumux/freeflow/navierstokes/velocityoutput.hh
@@ -0,0 +1,103 @@
+// -*- 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 3 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
+ * \ingroup StaggeredDiscretization
+ * \copydoc Dumux::StaggeredFreeFlowVelocityOutput
+ */
+#ifndef DUMUX_FREEFLOW_NAVIERSTOKES_VELOCITYOUTPUT_HH
+#define DUMUX_FREEFLOW_NAVIERSTOKES_VELOCITYOUTPUT_HH
+
+#include <dumux/io/velocityoutput.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/freeflow/navierstokes/momentum/velocityreconstruction.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup StaggeredDiscretization
+ * \brief Velocity output for staggered free-flow models
+ */
+template<class GridVariables>
+class NavierStokesVelocityOutput : public VelocityOutput<GridVariables>
+{
+ using ParentType = VelocityOutput<GridVariables>;
+ using GridGeometry = typename GridVariables::GridGeometry;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+ using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+ using ElementFluxVarsCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+ using VolumeVariables = typename GridVariables::VolumeVariables;
+ using FluidSystem = typename VolumeVariables::FluidSystem;
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+ using VelocityVector = typename ParentType::VelocityVector;
+
+ NavierStokesVelocityOutput(const std::string& paramGroup = "")
+ {
+ enableOutput_ = getParamFromGroup<bool>(paramGroup, "Vtk.AddVelocity", true);
+ }
+
+ //! Returns whether to enable the velocity output or not
+ bool enableOutput() const override { return enableOutput_; }
+
+ //! returns the phase name of a given phase index
+ std::string phaseName(int phaseIdx) const override { return FluidSystem::phaseName(phaseIdx); }
+
+ //! returns the number of phases
+ int numFluidPhases() const override { return VolumeVariables::numFluidPhases(); }
+
+ //! Calculate the velocities for the scvs in the element
+ //! We assume the local containers to be bound to the complete stencil
+ void calculateVelocity(VelocityVector& velocity,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVarsCache& elemFluxVarsCache,
+ int phaseIdx) const override
+ {
+ using Problem = std::decay_t<decltype(elemVolVars.gridVolVars().problem())>;
+ if constexpr (Problem::momentumDiscretizationMethod == DiscretizationMethod::fcstaggered)
+ calculateVelocityForStaggeredGrid_(velocity, element, fvGeometry, elemVolVars);
+ }
+
+private:
+ void calculateVelocityForStaggeredGrid_(VelocityVector& velocity,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars) const
+ {
+ const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element);
+ const auto getFaceVelocity = [&](const FVElementGeometry& fvG, const auto& scvf)
+ {
+ return elemVolVars.gridVolVars().problem().faceVelocity(element, fvGeometry, scvf);
+ };
+
+ velocity[eIdx] = StaggeredVelocityReconstruction::cellCenterVelocity(getFaceVelocity, fvGeometry);
+ }
+
+bool enableOutput_;
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab