Commit 89dc5853 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

add the M-phase, N-component implicit model

so far the "new" fluid systems have been copied but are not used by
any other model. some discussion still needed...

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6657 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent e323d1c4
......@@ -29,6 +29,7 @@ AC_CONFIG_FILES([dumux.pc
dumux/boxmodels/2p2cni/Makefile
dumux/boxmodels/2pni/Makefile
dumux/boxmodels/common/Makefile
dumux/boxmodels/MpNc/Makefile
dumux/boxmodels/richards/Makefile
dumux/common/Makefile
dumux/decoupled/Makefile
......@@ -50,9 +51,14 @@ AC_CONFIG_FILES([dumux.pc
dumux/material/binarycoefficients/Makefile
dumux/material/components/Makefile
dumux/material/components/iapws/Makefile
dumux/material/constraintsolvers/Makefile
dumux/material/eos/Makefile
dumux/material/fluidmatrixinteractions/Makefile
dumux/material/fluidmatrixinteractions/2p/Makefile
dumux/material/fluidmatrixinteractions/Mp/Makefile
dumux/material/fluidstates/Makefile
dumux/material/fluidsystems/Makefile
dumux/material/new_fluidsystems/Makefile
dumux/material/spatialparameters/Makefile
dumux/nonlinear/Makefile
m4/Makefile
......@@ -64,6 +70,7 @@ AC_CONFIG_FILES([dumux.pc
test/boxmodels/1p2c/Makefile
test/boxmodels/2p2c/Makefile
test/boxmodels/2p2cni/Makefile
test/boxmodels/MpNc/Makefile
test/boxmodels/richards/Makefile
test/common/Makefile
test/common/generalproblem/Makefile
......
SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni richards
SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni MpNc richards
boxmodelsdir = $(includedir)/dumux/boxmodels
......
MpNcdir = $(includedir)/dumux/boxmodels/MpNc
MpNc_HEADERS = *.hh
include $(top_srcdir)/am/global-rules
/*****************************************************************************
* Copyright (C) 2009-2010 by Andreas Lauser *
* Copyright (C) 2008-2009 by Klaus Mosthaf *
* Copyright (C) 2008 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* 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 data which is required to calculate
* all fluxes of components over a face of a finite volume.
*
* This means pressure, concentration and temperature gradients, phase
* densities at the integration point, etc.
*/
#ifndef DUMUX_MPNC_FLUX_VARIABLES_HH
#define DUMUX_MPNC_FLUX_VARIABLES_HH
#include <dumux/common/spline.hh>
#include "diffusion/fluxvariables.hh"
#include "energy/MpNcfluxvariablesenergy.hh"
namespace Dumux
{
/*!
* \brief This template class contains the data which is required to
* calculate all fluxes of components over a face of a finite
* volume for the two-phase, three-component model.
*
* This means pressure and concentration gradients, phase densities at
* the intergration point, etc.
*/
template <class TypeTag>
class MPNCFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
enum {
dimWorld = GridView::dimensionworld,
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
enableDiffusion = GET_PROP_VALUE(TypeTag, PTAG(EnableDiffusion)),
enableEnergy = GET_PROP_VALUE(TypeTag, PTAG(EnableEnergy)),
enableKinetic = GET_PROP_VALUE(TypeTag, PTAG(EnableKinetic)),
enableKineticEnergy = GET_PROP_VALUE(TypeTag, PTAG(EnableKineticEnergy)),
enableGravity = GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)),
};
typedef typename GridView::ctype CoordScalar;
typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolume SCV;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef MPNCFluxVariablesDiffusion<TypeTag, enableDiffusion> FluxVariablesDiffusion;
typedef MPNCFluxVariablesEnergy<TypeTag, enableEnergy, enableKineticEnergy> FluxVariablesEnergy;
public:
MPNCFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int scvfIdx,
const ElementVolumeVariables &elemVolVars)
: fvElemGeom_(elemGeom), volVars_(elemVolVars)
{
scvfIdx_ = scvfIdx;
// update the base module (i.e. advection)
calculateGradients_(problem, element, elemVolVars);
calculateVelocities_(problem, element, elemVolVars);
// update the flux data of the energy module (i.e. isothermal
// or non-isothermal)
energyDat_.update(problem, element, elemGeom, scvfIdx, *this, elemVolVars);
// update the flux data of the diffusion module (i.e. with or
// without diffusion)
diffusionDat_.update(problem, element, elemGeom, scvfIdx, elemVolVars);
extrusionFactor_ =
(elemVolVars[face().i].extrusionFactor()
+ elemVolVars[face().j].extrusionFactor()) / 2;
}
/*!
* \brief Calculate a phase's darcy velocity [m/s] at a
* sub-control volume face.
*
* So far, this method only exists in the Mp-Nc model!
*
* Of course, in the setting of a finite volume scheme, the velocities are
* on the faces rather than in the volume. Therefore, the velocity
*
* \param problem The problem object
* \param element The current element
* \param elemGeom The finite-volume geometry in the box scheme
* \param scvIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
void computeDarcy(Vector & vDarcy,
const ElementVolumeVariables &elemVolVars,
int phaseIdx) const
{
intrinsicPermeability().mv(potentialGrad(phaseIdx),
vDarcy);
// darcy velocity is along *negative* potential gradient
// (i.e. from high to low pressures), this means that we need
// to negate the product of the intrinsic permeability and the
// potential gradient!
vDarcy *= -1;
// JUST for upstream decision
Scalar normalFlux = vDarcy * face().normal;
// data attached to upstream and the downstream vertices
// of the current phase
int upIdx = face().i;
int dnIdx = face().j;
if (!std::isfinite(normalFlux))
DUNE_THROW(NumericalProblem, "Calculated non-finite normal flux");
if (normalFlux < 0)
std::swap(upIdx, dnIdx);
const VolumeVariables &up = elemVolVars[upIdx];
////////
// Jipie this is a velocity now, finally deserves the name
////////
vDarcy *= up.mobility(phaseIdx);
}
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
const VolumeVariables &volVars(int idx) const
{ return volVars_[idx]; }
/*!
* \brief Returns th extrusion factor for the sub-control volume face
*/
Scalar extrusionFactor() const
{ return extrusionFactor_; }
/*!
* \brief Return the intrinsic permeability.
*/
const Tensor &intrinsicPermeability() const
{ return K_; }
/*!
* \brief Return the pressure potential gradient.
*/
const Vector &potentialGrad(int phaseIdx) const
{ return potentialGrad_[phaseIdx]; }
////////////////////////////////////////////////
// forward calls to the diffusion module
Scalar porousDiffCoeffL(int compIdx) const
{ return diffusionDat_.porousDiffCoeffL(compIdx); };
Scalar porousDiffCoeffG(int compIIdx, int compJIdx) const
{ return diffusionDat_.porousDiffCoeffG(compIIdx, compJIdx); };
const Scalar moleFrac(int phaseIdx, int compIdx) const
{ return diffusionDat_.moleFrac(phaseIdx, compIdx); };
const Vector &moleFracGrad(int phaseIdx,
int compIdx) const
{ return diffusionDat_.moleFracGrad(phaseIdx, compIdx); };
// end of forward calls to the diffusion module
////////////////////////////////////////////////
////////////////////////////////////////////////
// forward calls to the temperature module
const Vector &temperatureGrad() const
{ return energyDat_.temperatureGrad(); };
const FluxVariablesEnergy &energyData() const
{ return energyDat_; }
// end of forward calls to the temperature module
////////////////////////////////////////////////
private:
void calculateGradients_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
potentialGrad_[phaseIdx] = Scalar(0);
}
// calculate pressure gradients using finite element gradients
Vector tmp(0.0);
for (int idx = 0;
idx < fvElemGeom_.numVertices;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
const Vector &feGrad = face().grad[idx];
// TODO: only calculate the gradients for the present
// phases.
//
// compute sum of pressure gradients for each phase
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
{
// the pressure gradient
tmp = feGrad;
tmp *= elemVolVars[idx].fluidState().pressure(phaseIdx);
potentialGrad_[phaseIdx] += tmp;
}
}
///////////////
// correct the pressure gradients by the gravitational acceleration
///////////////
if (enableGravity) {
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
Vector g(problem.boxGravity(element, fvElemGeom_, face().i));
g += problem.boxGravity(element, fvElemGeom_, face().j);
g /= 2;
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
// calculate the phase density at the integration point. we
// only do this if the wetting phase is present in both cells
Scalar SI = elemVolVars[face().i].fluidState().saturation(phaseIdx);
Scalar SJ = elemVolVars[face().j].fluidState().saturation(phaseIdx);
Scalar rhoI = elemVolVars[face().i].fluidState().density(phaseIdx);
Scalar rhoJ = elemVolVars[face().j].fluidState().density(phaseIdx);
Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
if (fI + fJ == 0)
// doesn't matter because no wetting phase is present in
// both cells!
fI = fJ = 0.5;
Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
// make gravity acceleration a force
Vector f(g);
f *= density;
// calculate the final potential gradient
potentialGrad_[phaseIdx] -= f;
}
}
}
void calculateVelocities_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// multiply the pressure potential with the intrinsic
// permeability
const SpatialParameters &sp = problem.spatialParameters();
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
sp.meanK(K_,
sp.intrinsicPermeability(element,
fvElemGeom_,
face().i),
sp.intrinsicPermeability(element,
fvElemGeom_,
face().j));
}
}
const FVElementGeometry &fvElemGeom_;
int scvfIdx_;
const ElementVolumeVariables &volVars_;
// The extrusion factor for the sub-control volume face
Scalar extrusionFactor_;
// pressure potential gradients
Vector potentialGrad_[numPhases];
// intrinsic permeability
Tensor K_;
FluxVariablesDiffusion diffusionDat_;
FluxVariablesEnergy energyDat_;
};
} // end namepace
#endif
/*****************************************************************************
* Copyright (C) 2008-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* 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/>. *
*****************************************************************************/
#ifndef DUMUX_MPNC_INDICES_HH
#define DUMUX_MPNC_INDICES_HH
#include "MpNcproperties.hh"
#include "mass/MpNcindicesmass.hh"
#include "energy/MpNcindicesenergy.hh"
namespace Dumux
{
/*!
* \brief The primary variable and equation indices for the MpNc
* model.
*/
template <class TypeTag, int BasePVOffset = 0>
struct MPNCIndices :
public MPNCMassIndices<BasePVOffset,
TypeTag,
GET_PROP_VALUE(TypeTag, PTAG(EnableKinetic)) >,
public MPNCEnergyIndices<BasePVOffset +
MPNCMassIndices<0, TypeTag, GET_PROP_VALUE(TypeTag, PTAG(EnableKinetic)) >::NumPrimaryVars,
GET_PROP_VALUE(TypeTag, PTAG(EnableEnergy)),
GET_PROP_VALUE(TypeTag, PTAG(EnableKineticEnergy))>
{
private:
enum { enableEnergy = GET_PROP_VALUE(TypeTag, PTAG(EnableEnergy)) };
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, PTAG(EnableDiffusion)) };
enum { enableKinetic = GET_PROP_VALUE(TypeTag, PTAG(EnableKinetic)) }; //mass transfer
enum { enableKineticEnergy = GET_PROP_VALUE(TypeTag, PTAG(EnableKineticEnergy)) }; // energy transfer
typedef MPNCMassIndices<BasePVOffset, TypeTag, enableKinetic> MassIndices;
typedef MPNCEnergyIndices<BasePVOffset + MassIndices::NumPrimaryVars, enableEnergy, enableKineticEnergy> EnergyIndices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
enum { numComponents = FluidSystem::numComponents };
enum { numPhases = FluidSystem::numPhases };
public:
/*!
* \brief The number of primary variables / equations.
*/
// temperature + Mass Balance + constraints for switch stuff
static const int NumPrimaryVars =
MassIndices::NumPrimaryVars +
EnergyIndices::NumPrimaryVars +
numPhases;
/*!
* \brief The number of primary variables / equations of the erngy module.
*/
static const int NumPrimaryEnergyVars =
EnergyIndices::NumPrimaryVars ;
/*!
* \brief Index of the saturation of the first phase in a vector
* of primary variables.
*
* The following (numPhases - 1) primary variables represent the
* saturations for the phases [1, ..., numPhases - 1]
*/
static const int S0Idx =
MassIndices::NumPrimaryVars +
EnergyIndices::NumPrimaryVars;
/*!
* \brief Index of the first phase' pressure in a vector of
* primary variables.
*/
static const int p0Idx =
MassIndices::NumPrimaryVars +
EnergyIndices::NumPrimaryVars +
numPhases - 1;
/*!
* \brief Index of the first phase NCP equation.
*
* The index for the remaining phases are consecutive.
*/
static const int phase0NcpIdx =
MassIndices::NumPrimaryVars +
EnergyIndices::NumPrimaryVars;
};
}
#endif
/*****************************************************************************
* Copyright (C) 2009-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* 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/>. *
*****************************************************************************/
#ifndef DUMUX_MPNC_LOCAL_RESIDUAL_HH
#define DUMUX_MPNC_LOCAL_RESIDUAL_HH
#include "MpNcfluxvariables.hh"
#include "diffusion/diffusion.hh"
#include "energy/MpNclocalresidualenergy.hh"
#include "mass/MpNclocalresidualmass.hh"
#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
namespace Dumux
{
/*!
* \brief two-phase, N-component specific details needed to
* approximately calculate the local defect in the BOX scheme.
*
* This class is used to fill the gaps in BoxLocalResidual for the
* two-phase, N-component twophase flow.
*/
template<class TypeTag>
class MPNCLocalResidual : public BoxLocalResidual<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCIndices)) Indices;
protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
typedef MPNCLocalResidual<TypeTag> ThisType;
typedef BoxLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
enableEnergy = GET_PROP_VALUE(TypeTag, PTAG(EnableEnergy)),
enableKineticEnergy = GET_PROP_VALUE(TypeTag, PTAG(EnableKineticEnergy)),
enableDiffusion = GET_PROP_VALUE(TypeTag, PTAG(EnableDiffusion)),
enableKinetic = GET_PROP_VALUE(TypeTag, PTAG(EnableKinetic)),
enableSmoothUpwinding = GET_PROP_VALUE(TypeTag, PTAG(EnableSmoothUpwinding)),
phase0NcpIdx = Indices::phase0NcpIdx
};
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GridView::CollectiveCommunication CollectiveCommunication;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyResid;
typedef MPNCLocalResidualMass<TypeTag, enableKinetic> MassResid;
public:
/*!
* \brief Evaluate the amount all conservation quantites
* (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)
*/
void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
{
// if flag usePrevSol is set, the solution from the previous
// time step is used, otherwise the current solution is
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit euler method.
const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
// compute mass and energy storage terms
MassResid::computeStorage(storage, volVars);
EnergyResid::computeStorage(storage, volVars);
Valgrind::CheckDefined(storage);
}
/*!
* \brief Evaluate the amount all conservation quantites
* (e.g. phase mass) within all sub-control volumes of an
* element.
*/
void addPhaseStorage(PrimaryVariables &<