Commit 7bbcf3b8 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[freeflow] Introduce EnergyLocalResidual

parent cb60ba16
......@@ -44,7 +44,6 @@ template <class TypeTag>
class FouriersLawImplementation<TypeTag, DiscretizationMethod::staggered >
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
......@@ -63,22 +62,14 @@ public:
//! We don't cache anything for this law
using Cache = FluxVariablesCaching::EmptyDiffusionCache;
//! calculate the molecular diffusive fluxes
static Scalar diffusiveFluxForCellCenter(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf)
//! calculate the diffusive energy fluxes
static Scalar flux(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf)
{
Scalar flux(0.0);
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
if(bcTypes.isOutflow(energyBalanceIdx) || bcTypes.isNeumann(energyBalanceIdx))
return flux;
}
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
......
......@@ -59,6 +59,8 @@ class FreeflowNCResidualImpl<TypeTag, BaseLocalResidual, DiscretizationMethod::s
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
static constexpr auto cellCenterOffset = ParentType::cellCenterOffset;
using EnergyLocalResidual = typename ParentType::EnergyLocalResidual;
public:
using ParentType::ParentType;
......@@ -88,8 +90,7 @@ public:
if(Indices::replaceCompEqIdx < numComponents)
storage[Indices::replaceCompEqIdx] = density;
this->computeStorageForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance() >(),
problem, scv, volVars, storage);
EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
return storage;
}
......
......@@ -42,9 +42,8 @@ class FreeflowNCVolumeVariables : public FreeFlowVolumeVariables< Traits, Freefl
using ParentType = FreeFlowVolumeVariables<Traits, ThisType>;
using Scalar = typename Traits::PrimaryVariables::value_type;
using Indices = typename Traits::ModelTraits::Indices;
static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
static constexpr int numComponents = Traits::ModelTraits::numComponents();
static constexpr bool useMoles = Traits::ModelTraits::useMoles();
......@@ -54,6 +53,8 @@ public:
using FluidSystem = typename Traits::FluidSystem;
//! export the fluid state type
using FluidState = typename Traits::FluidState;
//! export the indices type
using Indices = typename Traits::ModelTraits::Indices;
/*!
* \brief Update all quantities for a given control volume
......
......@@ -232,9 +232,6 @@ public:
using type = FreeflowNonIsothermalVtkOutputFields<IsothermalFields, ModelTraits>;
};
//! Use Fourier's Law as default heat conduction type
SET_TYPE_PROP(NavierStokesNI, HeatConductionType, FouriersLaw<TypeTag>);
// \}
}
......
......@@ -84,6 +84,8 @@ class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
public:
using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
/*!
* \brief Returns the advective flux over a sub control volume face.
* \param elemVolVars All volume variables for the element
......
......@@ -28,6 +28,7 @@
#include <dumux/discretization/methods.hh>
#include <dumux/assembly/staggeredlocalresidual.hh>
#include <dune/common/hybridutilities.hh>
#include <dumux/freeflow/nonisothermal/localresidual.hh>
namespace Dumux
{
......@@ -79,6 +80,7 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
public:
using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance()>;
// account for the offset of the cell center privars within the PrimaryVariables container
static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
......@@ -100,8 +102,7 @@ public:
CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars,
elemFaceVars, scvf, elemFluxVarsCache[scvf]);
computeFluxForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance()>(),
problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache, flux);
EnergyLocalResidual::heatFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
return flux;
}
......@@ -135,8 +136,7 @@ public:
CellCenterPrimaryVariables storage;
storage[Indices::conti0EqIdx - ModelTraits::dim()] = volVars.density();
computeStorageForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance() >(),
problem, scv, volVars, storage);
EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
return storage;
}
......@@ -204,24 +204,6 @@ public:
protected:
// /*!
// * \brief Evaluate boundary conditions
// */
// template<class ElementBoundaryTypes>
// void evalBoundary_(const Element& element,
// const FVElementGeometry& fvGeometry,
// const ElementVolumeVariables& elemVolVars,
// const ElementFaceVariables& elemFaceVars,
// const ElementBoundaryTypes& elemBcTypes,
// const ElementFluxVariablesCache& elemFluxVarsCache)
// {
// evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
// for (auto&& scvf : scvfs(fvGeometry))
// {
// evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
// }
// }
/*!
* \brief Evaluate boundary conditions for a cell center dof
*/
......@@ -245,6 +227,7 @@ protected:
// handle the actual boundary conditions:
const auto bcTypes = problem.boundaryTypes(element, scvf);
EnergyLocalResidual::heatFlux(boundaryFlux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
if(bcTypes.hasNeumann())
{
static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
......@@ -327,52 +310,6 @@ protected:
}
}
//! Evaluate energy fluxes entering or leaving the cell center control volume for non isothermal models
void computeFluxForCellCenterNonIsothermal_(std::true_type,
const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf,
const ElementFluxVariablesCache& elemFluxVarsCache,
CellCenterPrimaryVariables& flux) const
{
// if we are on an inflow/outflow boundary, use the volVars of the element itself
// TODO: catch neumann and outflow in localResidual's evalBoundary_()
bool isOutflow = false;
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
if(bcTypes.isOutflow(Indices::energyBalanceIdx))
isOutflow = true;
}
auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); };
using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
flux[Indices::energyBalanceIdx - cellCenterOffset] = FluxVariables::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
flux[Indices::energyBalanceIdx - cellCenterOffset] += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf);
}
//! Evaluate energy fluxes entering or leaving the cell center control volume for non isothermal models
template <typename... Args>
void computeFluxForCellCenterNonIsothermal_(std::false_type, Args&&... args) const {}
//! Evaluate energy storage for non isothermal models
void computeStorageForCellCenterNonIsothermal_(std::true_type,
const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars,
CellCenterPrimaryVariables& storage) const
{
storage[Indices::energyBalanceIdx - cellCenterOffset] = volVars.density() * volVars.internalEnergy();
}
//! Evaluate energy storage for isothermal models
template <typename... Args>
void computeStorageForCellCenterNonIsothermal_(std::false_type, Args&&... args) const {}
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
......
......@@ -40,15 +40,16 @@ class NavierStokesVolumeVariables : public FreeFlowVolumeVariables< Traits, Navi
using ParentType = FreeFlowVolumeVariables<Traits, ThisType>;
using Scalar = typename Traits::PrimaryVariables::value_type;
using Indices = typename Traits::ModelTraits::Indices;
static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
public:
//! export the underlying fluid system
using FluidSystem = typename Traits::FluidSystem;
//! export the fluid state type
using FluidState = typename Traits::FluidState;
//! export the indices type
using Indices = typename Traits::ModelTraits::Indices;
/*!
* \brief Update all quantities for a given control volume
......
// -*- 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
* \ingroup FreeflowNIModel
* \copydoc Dumux::FreeFlowEnergyLocalResidual
*/
#ifndef DUMUX_FREE_FLOW_ENERGY_LOCAL_RESIDUAL_HH
#define DUMUX_FREE_FLOW_ENERGY_LOCAL_RESIDUAL_HH
#include <dumux/discretization/methods.hh>
namespace Dumux {
// forward declaration
template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod, bool enableEneryBalance>
class FreeFlowEnergyLocalResidualImplementation;
/*!
* \ingroup FreeflowNIModel
* \brief Element-wise calculation of the local residual for non-isothermal
* free-flow models
*/
template<class FVGridGeometry, class FluxVariables, bool enableEneryBalance>
using FreeFlowEnergyLocalResidual =
FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
FluxVariables,
FVGridGeometry::discMethod,
enableEneryBalance>;
/*!
* \ingroup FreeflowNIModel
* \brief Specialization for isothermal models, does nothing
*/
template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod>
class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry, FluxVariables, discMethod, false>
{
public:
//! do nothing for the isothermal case
template <typename... Args>
static void fluidPhaseStorage(Args&&... args)
{}
//! do nothing for the isothermal case
template <typename... Args>
static void heatConvectionFlux(Args&&... args)
{}
//! do nothing for the isothermal case
template <typename... Args>
static void heatFlux(Args&&... args)
{}
};
/*!
* \ingroup FreeflowNIModel
* \brief Specialization for staggered non-isothermal models
*/
template<class FVGridGeometry, class FluxVariables>
class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
FluxVariables,
DiscretizationMethod::staggered,
true>
{
using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using HeatConductionType = typename FluxVariables::HeatConductionType;
public:
//! The energy storage in the fluid phase
template<class NumEqVector, class VolumeVariables>
static void fluidPhaseStorage(NumEqVector& storage,
const VolumeVariables& volVars)
{
static constexpr auto localEnergyBalanceIdx = NumEqVector::dimension - 1;
storage[localEnergyBalanceIdx] += volVars.density() * volVars.internalEnergy();
}
//! The convective and conductive heat fluxes in the fluid phase
template<class NumEqVector, class Problem, class ElementVolumeVariables, class ElementFaceVariables>
static void heatFlux(NumEqVector& flux,
const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf)
{
static constexpr auto localEnergyBalanceIdx = NumEqVector::dimension - 1;
using Indices = typename ElementVolumeVariables::VolumeVariables::Indices;
bool isOutflow = false;
if(scvf.boundary())
{
const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
if(bcTypes.isOutflow(Indices::energyBalanceIdx))
isOutflow = true;
}
auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); };
flux[localEnergyBalanceIdx] += FluxVariables::advectiveFluxForCellCenter(elemVolVars,
elemFaceVars,
scvf,
upwindTerm,
isOutflow);
if(!isOutflow)
flux[localEnergyBalanceIdx] += HeatConductionType::flux(element,
fvGeometry,
elemVolVars,
scvf);
}
};
} // end namespace Dumux
#endif
......@@ -28,6 +28,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/properties/model.hh>
#include <dumux/discretization/fourierslaw.hh>
namespace Dumux
{
......@@ -36,6 +37,9 @@ namespace Properties
//! Type tag for free-flow models
NEW_TYPE_TAG(FreeFlow, INHERITS_FROM(ModelProperties));
//! Use Fourier's Law as default heat conduction type
SET_TYPE_PROP(FreeFlow, HeatConductionType, FouriersLaw<TypeTag>);
} // namespace Properties
} // namespace Dumux
......
......@@ -47,14 +47,15 @@ class LowReKEpsilonVolumeVariables
using NavierStokesParentType = NSVolumeVariables;
using Scalar = typename Traits::PrimaryVariables::value_type;
using Indices = typename Traits::ModelTraits::Indices;
static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
public:
//! export the underlying fluid system
using FluidSystem = typename Traits::FluidSystem;
//! export the indices type
using Indices = typename Traits::ModelTraits::Indices;
/*!
* \brief Update all quantities for a given control volume
......
......@@ -46,16 +46,17 @@ class ZeroEqVolumeVariables
using NavierStokesParentType = NSVolumeVariables;
using Scalar = typename Traits::PrimaryVariables::value_type;
using Indices = typename Traits::ModelTraits::Indices;
static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
static constexpr int fluidSystemPhaseIdx = Indices::fluidSystemPhaseIdx;
static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
public:
//! export the underlying fluid system
using FluidSystem = typename Traits::FluidSystem;
//! export the fluid state type
using FluidState = typename Traits::FluidState;
//! export the indices type
using Indices = typename Traits::ModelTraits::Indices;
/*!
* \brief Update all quantities for a given control volume
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment