Skip to content
GitLab
Commit 9f0f7bd5 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[folder structure] rename and move files from implicit/mpnc, adapt includes

parent 4c3fde88
Expand all Hide whitespace changes
Inline Side-by-side
dumux/implicit/mpnc/diffusion/diffusion.hh
View file @ 9f0f7bd5
// -*- 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 This file contains parts to calculate the diffusive flux in
* the fully coupled MpNc model
*/
#ifndef DUMUX_MPNC_DIFFUSION_HH
#define DUMUX_MPNC_DIFFUSION_HH
#ifndef DUMUX_MPNC_DIFFUSION_HH_OLD
#define DUMUX_MPNC_DIFFUSION_HH_OLD
#include <dune/common/float_cmp.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#warning this header is deprecated, use dumux/porousmediumflow/mpnc/implicit/diffusion/diffusion.hh instead
#include <dumux/porousmediumflow/mpnc/implicit/diffusion/diffusion.hh>
#include <dumux/implicit/mpnc/mpncproperties.hh>
namespace Dumux {
/*!
* \brief Calculates the diffusive flux in the fully coupled MpNc model.
* Called from the mass module.
*/
template <class TypeTag, bool enableDiffusion>
class MPNCDiffusion
{
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum { nPhaseIdx = FluidSystem::nPhaseIdx };
enum { wPhaseIdx = FluidSystem::wPhaseIdx };
enum { nCompIdx = FluidSystem::nCompIdx };
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef Dune::FieldMatrix<Scalar, numComponents, numComponents> ComponentMatrix;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
public:
/*!
* \brief Calls the function for the diffusion in the gas and liquid phases, respectively.
*
* In the gas phase, the mutual influence of mole fractions can be considered.
*
* \param fluxes The Diffusive flux over the sub-control-volume face for each component
* \param phaseIdx The index of the phase we are calculating the diffusive flux for
* \param fluxVars The flux variables at the current sub control volume face
* \param molarDensity The (molar) density of the phase
*/
static void flux(ComponentVector &fluxes,
const unsigned int phaseIdx,
const FluxVariables &fluxVars,
const Scalar molarDensity )
{
if ( not FluidSystem::isLiquid(phaseIdx) )
gasFlux_(fluxes, fluxVars, molarDensity);
else if ( FluidSystem::isLiquid(phaseIdx) ){
#if MACROSCALE_DIFFUSION_ONLY_GAS
return ; // in the case that only the diffusion in the gas phase is considered,
// the liquidFlux should not be called
#endif
liquidFlux_(fluxes, fluxVars, molarDensity);
}
else
DUNE_THROW(Dune::InvalidStateException,
"Invalid phase index: " << phaseIdx);
}
protected:
/*!
* \brief Calculates the diffusive flux in the liquid phase: only Fick diffusion (no mutual influence) is considered.
*
* \param fluxes The Diffusive flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current sub control volume face
* \param molarDensity The (molar) density of the phase
*/
static void liquidFlux_(ComponentVector &fluxes,
const FluxVariables &fluxVars,
const Scalar molarDensity)
{
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
// TODO: tensorial diffusion coefficients
const Scalar xGrad = fluxVars.moleFractionGrad(wPhaseIdx, compIdx)*fluxVars.face().normal;
fluxes[compIdx] =
- xGrad *
molarDensity *
fluxVars.porousDiffCoeffL(compIdx) ;
}
}
/*!
* \brief Calculates the diffusive flux in the gas phase:
* The property UseMaxwellDiffusion selects whether Maxwell or Fick diffusion is applied.
* Has to be the same for a 2-component system. However, Maxwells does not work always.
*
* \param fluxes The Diffusive flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current sub control volume face
* \param molarDensity The (molar) density of the phase
*
* Reid et al. (1987, p. 596) \cite reid1987 <BR>
*/
static void gasFlux_(ComponentVector &fluxes,
const FluxVariables &fluxVars,
const Scalar molarDensity)
{
// Alternative: use Fick Diffusion: no mutual influence of diffusion
if (GET_PROP_VALUE(TypeTag, UseMaxwellDiffusion) ){
// Stefan-Maxwell equation
//
// See: R. Reid, et al.: "The Properties of Liquids and
// Gases", 4th edition, 1987, McGraw-Hill, p 596
// TODO: tensorial diffusion coefficients
ComponentMatrix M(0);
for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) {
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
Scalar Dij = fluxVars.porousDiffCoeffG(compIIdx, compJIdx);
if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(Dij, 0.0, 1.0e-30)) {
M[compIIdx][compJIdx] += fluxVars.moleFraction(nPhaseIdx, compIIdx) / Dij;
M[compIIdx][compIIdx] -= fluxVars.moleFraction(nPhaseIdx, compJIdx) / Dij;
}
}
}
for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
M[numComponents - 1][compIIdx] = 1.0;
}
ComponentVector rightHandSide ; // see source cited above
for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) {
rightHandSide[compIIdx] = molarDensity*(fluxVars.moleFractionGrad(nPhaseIdx, compIIdx)*fluxVars.face().normal);
}
rightHandSide[numComponents - 1] = 0.0;
M.solve(fluxes, rightHandSide);
}
else{// Fick Diffusion
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
// TODO: tensorial diffusion coefficients
const Scalar xGrad = fluxVars.moleFractionGrad(nPhaseIdx, compIdx)*fluxVars.face().normal;
fluxes[compIdx] =
- xGrad *
molarDensity
* fluxVars.porousDiffCoeffG(compIdx, nCompIdx) ; // this is == 0 for nComp==comp,
// i.e. no diffusion of the main component of the phase
}
}
}
// return whether a concentration can be assumed to be a trace
// component in the context of diffusion
static Scalar isTraceComp_(Scalar x)
{ return x < 0.5/numComponents; }
};
/*!
* \brief Specialization of the diffusion module for the case where
* diffusion is disabled.
*
* This class just does nothing.
*/
template <class TypeTag>
class MPNCDiffusion<TypeTag, false>
{
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
public:
static void flux(ComponentVector &fluxes,
const unsigned int phaseIdx,
const FluxVariables &fluxVars,
const Scalar molarDensity)
{ fluxes = 0; }
};
}
#endif // DUMUX_MPNC_DIFFUSION_HH
#endif
dumux/implicit/mpnc/diffusion/fluxvariables.hh
View file @ 9f0f7bd5
// -*- 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 This file contains the diffusion module for the flux data of
* the fully coupled MpNc model
*/
#ifndef DUMUX_MPNC_DIFFUSION_FLUX_VARIABLES_HH
#define DUMUX_MPNC_DIFFUSION_FLUX_VARIABLES_HH
#ifndef DUMUX_MPNC_DIFFUSION_FLUX_VARIABLES_HH_OLD
#define DUMUX_MPNC_DIFFUSION_FLUX_VARIABLES_HH_OLD
#include <dune/common/fvector.hh>
#warning this header is deprecated, use dumux/porousmediumflow/mpnc/implicit/diffusion/fluxvariables.hh instead
#include "../mpncproperties.hh"
#include <dumux/porousmediumflow/mpnc/implicit/diffusion/fluxvariables.hh>
namespace Dumux {
/*!
* \ingroup MPNCModel
* \ingroup ImplicitFluxVariables
* \brief Variables for the diffusive fluxes in the MpNc model
*/
template<class TypeTag, bool enableDiffusion>
class MPNCFluxVariablesDiffusion
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
enum{dim = GridView::dimension};
enum{dimWorld = GridView::dimensionworld};
enum{numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
enum{numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum{wPhaseIdx = FluidSystem::wPhaseIdx};
enum{nPhaseIdx = FluidSystem::nPhaseIdx};
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
/*!
* \brief The constructor
*/
MPNCFluxVariablesDiffusion()
{}
/*!
* \brief update
*
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param face The SCV (sub-control-volume) face
* \param elemVolVars The volume variables of the current element
*/
void update(const Problem & problem,
const Element & element,
const FVElementGeometry & fvGeometry,
const SCVFace & face,
const ElementVolumeVariables & elemVolVars)
{
const unsigned int i = face.i;
const unsigned int j = face.j;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx){
for (int compIdx = 0; compIdx < numComponents; ++compIdx){
moleFraction_[phaseIdx][compIdx] = 0. ;
moleFractionGrad_[phaseIdx][compIdx] = 0. ;
}
}
GlobalPosition tmp ;
for (unsigned int idx = 0;
idx < face.numFap;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
const GlobalPosition & feGrad = face.grad[idx];
// index for the element volume variables
int volVarsIdx = face.fapIndices[idx];
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++){
for (int compIdx = 0; compIdx < numComponents; ++compIdx){
// calculate mole fractions at the integration points of the face
moleFraction_[phaseIdx][compIdx] += elemVolVars[volVarsIdx].fluidState().moleFraction(phaseIdx, compIdx)*
face.shapeValue[idx];
// calculate mole fraction gradients
tmp = feGrad;
tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(phaseIdx, compIdx);
moleFractionGrad_[phaseIdx][compIdx] += tmp;
}
}
}
// initialize the diffusion coefficients to zero
for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
porousDiffCoeffL_[compIIdx] = 0.0;
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
porousDiffCoeffG_[compIIdx][compJIdx] = 0.0;
}
// calculate the diffusion coefficients at the integration
// point in the porous medium
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// make sure to only calculate diffusion coefficents
// for phases which exist in both finite volumes
if (elemVolVars[i].fluidState().saturation(phaseIdx) <= 1e-4 ||
elemVolVars[j].fluidState().saturation(phaseIdx) <= 1e-4)
{
continue;
}
// reduction factor for the diffusion coefficients in the
// porous medium in nodes i and j. this is the tortuosity
// times porosity times phase saturation at the nodes i
// and j
//
Scalar red_i =
elemVolVars[i].fluidState().saturation(phaseIdx)/elemVolVars[i].porosity() *
pow(elemVolVars[i].porosity() * elemVolVars[i].fluidState().saturation(phaseIdx), 7.0/3);
Scalar red_j =
elemVolVars[j].fluidState().saturation(phaseIdx)/elemVolVars[j].porosity() *
pow(elemVolVars[j].porosity() * elemVolVars[j].fluidState().saturation(phaseIdx), 7.0/3);
if (phaseIdx == wPhaseIdx) {
// Liquid phase diffusion coefficients in the porous medium
for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
// -> arithmetic mean
porousDiffCoeffL_[compIIdx]
= 1./2*(red_i * elemVolVars[i].diffCoeff(wPhaseIdx, 0, compIIdx) +
red_j * elemVolVars[j].diffCoeff(wPhaseIdx, 0, compIIdx));
}
}
else {
// Gas phase diffusion coefficients in the porous medium
for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
// -> arithmetic mean
porousDiffCoeffG_[compIIdx][compJIdx]
= 1./2*(red_i * elemVolVars[i].diffCoeff(nPhaseIdx, compIIdx, compJIdx) +
red_j * elemVolVars[j].diffCoeff(nPhaseIdx, compIIdx, compJIdx));
}
}
}
}
}
/*!
* \brief The binary diffusion coefficient for each component in the fluid phase.
* \param compIdx The local index of the components
*/
Scalar porousDiffCoeffL(const unsigned int compIdx) const
{
// TODO: tensorial diffusion coefficients
return porousDiffCoeffL_[compIdx];
}
/*!
* \brief The binary diffusion coefficient for each component in the gas phase.
* \param compIIdx The local index of the first component in the phase
* \param compJIdx The local index of the second component in the phase
*/
Scalar porousDiffCoeffG(const unsigned int compIIdx,
const unsigned int compJIdx) const
{
// TODO: tensorial diffusion coefficients
return porousDiffCoeffG_[compIIdx][compJIdx];
}
/*!
* \brief The mole fraction and concentration gradient for all phases and components
* \param phaseIdx The local index of the phases
* \param compIdx The local index of the component
*/
Scalar moleFraction(const unsigned int phaseIdx,
const unsigned int compIdx) const
{ return moleFraction_[phaseIdx][compIdx]; }
const GlobalPosition &moleFractionGrad(const unsigned int phaseIdx,
const unsigned int compIdx) const
{ return moleFractionGrad_[phaseIdx][compIdx];}
protected:
// the diffusion coefficients for the porous medium for the
// liquid phase
Scalar porousDiffCoeffL_[numComponents];
// the diffusion coefficients for the porous medium for the
// gas phase
Scalar porousDiffCoeffG_[numComponents][numComponents];
// the concentration gradients of all components in all phases
GlobalPosition moleFractionGrad_[numPhases][numComponents];
// the mole fractions of each component at the integration point
Scalar moleFraction_[numPhases][numComponents];
};
template<class TypeTag>
class MPNCFluxVariablesDiffusion<TypeTag, false>
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
public:
/*!
* \brief The constructor
*/
MPNCFluxVariablesDiffusion()
{}
/*!
* \brief update
*
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param face The SCV (sub-control-volume) face
* \param elemVolVars The volume variables of the current element
*/
void update(const Problem & problem,
const Element & element,
const FVElementGeometry & fvGeometry,
const SCVFace & face,
const ElementVolumeVariables & elemVolVars)
{
}
};
}
#endif // DUMUX_MPNC_DIFFUSION_FLUX_VARIABLES_HH
#endif
dumux/implicit/mpnc/diffusion/volumevariables.hh
View file @ 9f0f7bd5
// -*- 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 This file contains the diffusion module for the vertex data
* of the fully coupled MpNc model
*/
#ifndef DUMUX_MPNC_DIFFUSION_VOLUME_VARIABLES_HH
#define DUMUX_MPNC_DIFFUSION_VOLUME_VARIABLES_HH
#ifndef DUMUX_MPNC_DIFFUSION_VOLUME_VARIABLES_HH_OLD
#define DUMUX_MPNC_DIFFUSION_VOLUME_VARIABLES_HH_OLD
#include <dumux/common/valgrind.hh>
#include <dumux/implicit/mpnc/mpncproperties.hh>
#warning this header is deprecated, use dumux/porousmediumflow/mpnc/implicit/diffusion/volumevariables.hh instead
namespace Dumux {
#include <dumux/porousmediumflow/mpnc/implicit/diffusion/volumevariables.hh>
/*!
* \brief Variables for the diffusive fluxes in the MpNc model within
* a finite volume.
*/
template<class TypeTag, bool enableDiffusion>
class MPNCVolumeVariablesDiffusion
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename FluidSystem::ParameterCache ParameterCache;
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { wPhaseIdx = FluidSystem::wPhaseIdx };
enum { nPhaseIdx = FluidSystem::nPhaseIdx };
public:
/*!
* \brief The constructor
*/
MPNCVolumeVariablesDiffusion()
{}
/*!
* \brief update
*
* \param fluidState An arbitrary fluid state
* \param paramCache Container for cache parameters
* \param volVars The volume variables
* \param problem The problem
*/
void update(FluidState &fluidState,
ParameterCache &paramCache,
const VolumeVariables &volVars,
const Problem &problem)
{
Valgrind::SetUndefined(*this);
// diffusion coefficents in liquid
diffCoeffL_[0] = 0.0;
for (int compIdx = 1; compIdx < numComponents; ++compIdx) {
diffCoeffL_[compIdx] =
FluidSystem::binaryDiffusionCoefficient(fluidState,
paramCache,
wPhaseIdx,
0,
compIdx);
}
Valgrind::CheckDefined(diffCoeffL_);
// diffusion coefficents in gas
for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
diffCoeffG_[compIIdx][compIIdx] = 0;
for (int compJIdx = compIIdx + 1; compJIdx < numComponents; ++compJIdx) {
diffCoeffG_[compIIdx][compJIdx] =
FluidSystem::binaryDiffusionCoefficient(fluidState,
paramCache,
nPhaseIdx,
compIIdx,
compJIdx);
// fill the symmetric part of the diffusion coefficent
// matrix
diffCoeffG_[compJIdx][compIIdx] = diffCoeffG_[compIIdx][compJIdx];
}
}
Valgrind::CheckDefined(diffCoeffG_);
}
/*!
* \brief The binary diffusion coefficient for each fluid phase.
* \param phaseIdx The local index of the phases
* \param compIIdx The local index of the first component in the phase
* \param compJIdx The local index of the second component in the phase
*/
Scalar diffCoeff(const unsigned int phaseIdx,
const unsigned int compIIdx,
const unsigned int compJIdx) const
{
if (phaseIdx == nPhaseIdx)
// TODO: tensorial diffusion coefficients
return diffCoeffG_[compIIdx][compJIdx];
const unsigned int i = std::min(compIIdx, compJIdx);
const unsigned int j = std::max(compIIdx, compJIdx);
if (i != 0)
return 0;
return diffCoeffL_[j];