// -*- 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 Element-wise calculation of the local residual for problems
* using compositional fully implicit model.
*/
#ifndef DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
#define DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
namespace Dumux
{
/*!
* \ingroup Implicit
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the local residual for problems
* using compositional fully implicit model.
*
*/
template<class TypeTag>
class MyCompositionalLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
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;
using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
static constexpr int numPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases();
static constexpr int numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
enum { conti0EqIdx = Indices::conti0EqIdx };
public:
using ParentType::ParentType;
/*!
* \brief Evaluate the amount of all conservation quantities
* (e.g. phase mass) within a sub-control volume.
*
* The result should be averaged over the volume (e.g. phase mass
* inside a sub control volume divided by the volume)
*
* \param storage The mass of the component within the sub-control volume
* \param scvIdx The SCV (sub-control-volume) index
* \param usePrevSol Evaluate function with solution of current or previous time step
*/
ResidualVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
ResidualVector storage(0.0);
// compute storage term of all components within all phases
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
auto eqIdx = conti0EqIdx + compIdx;
storage[eqIdx] += volVars.porosity()
* volVars.saturation(phaseIdx)
* volVars.molarDensity(phaseIdx)
* volVars.moleFraction(phaseIdx, compIdx);
}
}
return storage;
}
/*!
* \brief Evaluates the total flux of all conservation quantities
* over a face of a sub-control volume.
*
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param fIdx The index of the SCV face
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
ResidualVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get upwind weights into local scope
ResidualVector flux(0.0);
auto enableDiffusion = getParam<bool>("Problem.EnableDiffusion", true);
// advective fluxes
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
// get equation index
const auto eqIdx = conti0EqIdx + compIdx;
// the physical quantities for which we perform upwinding
const auto upwindTerm = [phaseIdx,compIdx] (const auto& volVars)
{ return volVars.molarDensity(phaseIdx)*volVars.moleFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
// TODO: dumux-course-task
// same here: only add diffusive fluxes if diffusion is enabled
if (enableDiffusion)
flux[eqIdx] += diffusiveFluxes[compIdx];
}
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
}
//! Add diffusive energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
return flux;
}
protected:
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
const Implementation *asImp_() const
{ return static_cast<const Implementation *> (this); }
};
} // end namespace Dumux
#endif