mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
changed: ewoms/models/ncp -> opm/models/ncp
This commit is contained in:
Arne Morten Kvarving
parent
d302771e6c
commit
09cdd3aadd
@ -29,7 +29,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
|
||||
#include "problems/co2injectionproblem.hh"
|
||||
|
||||
|
||||
@ -29,7 +29,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
|
||||
#include "problems/co2injectionproblem.hh"
|
||||
|
||||
|
||||
@ -29,7 +29,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
|
||||
#include "problems/co2injectionproblem.hh"
|
||||
|
||||
|
||||
@ -29,7 +29,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
|
||||
|
||||
#include "problems/co2injectionproblem.hh"
|
||||
|
||||
@ -28,7 +28,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include "problems/diffusionproblem.hh"
|
||||
|
||||
BEGIN_PROPERTIES
|
||||
|
||||
@ -28,7 +28,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
|
||||
#include "problems/obstacleproblem.hh"
|
||||
|
||||
|
||||
@ -28,7 +28,7 @@
|
||||
#ifndef EWOMS_POWER_INJECTION_PROBLEM_HH
|
||||
#define EWOMS_POWER_INJECTION_PROBLEM_HH
|
||||
|
||||
#include <ewoms/models/ncp/ncpproperties.hh>
|
||||
#include <opm/models/ncp/ncpproperties.hh>
|
||||
|
||||
#include <opm/models/io/cubegridvanguard.hh>
|
||||
|
||||
|
||||
@ -28,7 +28,7 @@
|
||||
#ifndef EWOMS_OBSTACLE_PROBLEM_HH
|
||||
#define EWOMS_OBSTACLE_PROBLEM_HH
|
||||
|
||||
#include <ewoms/models/ncp/ncpproperties.hh>
|
||||
#include <opm/models/ncp/ncpproperties.hh>
|
||||
|
||||
#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
|
||||
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
|
||||
|
||||
@ -28,7 +28,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
|
||||
#include "problems/reservoirproblem.hh"
|
||||
|
||||
|
||||
@ -29,7 +29,7 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/models/utils/start.hh>
|
||||
#include <ewoms/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/ncp/ncpmodel.hh>
|
||||
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
|
||||
#include "problems/reservoirproblem.hh"
|
||||
|
||||
|
||||
203
opm/models/ncp/ncpboundaryratevector.hh
Normal file
203
opm/models/ncp/ncpboundaryratevector.hh
Normal file
@ -0,0 +1,203 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpBoundaryRateVector
|
||||
*/
|
||||
#ifndef EWOMS_NCP_BOUNDARY_RATE_VECTOR_HH
|
||||
#define EWOMS_NCP_BOUNDARY_RATE_VECTOR_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief Implements a boundary vector for the fully implicit
|
||||
* compositional multi-phase NCP model.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpBoundaryRateVector : public GET_PROP_TYPE(TypeTag, RateVector)
|
||||
{
|
||||
typedef typename GET_PROP_TYPE(TypeTag, RateVector) ParentType;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
|
||||
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
|
||||
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||
enum { conti0EqIdx = Indices::conti0EqIdx };
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
|
||||
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
|
||||
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
public:
|
||||
NcpBoundaryRateVector() : ParentType()
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \copydoc
|
||||
* ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(Scalar)
|
||||
*/
|
||||
NcpBoundaryRateVector(const Evaluation& value) : ParentType(value)
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(const
|
||||
* ImmiscibleBoundaryRateVector&)
|
||||
*/
|
||||
NcpBoundaryRateVector(const NcpBoundaryRateVector& value) = default;
|
||||
NcpBoundaryRateVector& operator=(const NcpBoundaryRateVector& value) = default;
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleBoundaryRateVector::setFreeFlow
|
||||
*/
|
||||
template <class Context, class FluidState>
|
||||
void setFreeFlow(const Context& context,
|
||||
unsigned bfIdx,
|
||||
unsigned timeIdx,
|
||||
const FluidState& fluidState)
|
||||
{
|
||||
ExtensiveQuantities extQuants;
|
||||
extQuants.updateBoundary(context, bfIdx, timeIdx, fluidState);
|
||||
const auto& insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx);
|
||||
unsigned focusDofIdx = context.focusDofIndex();
|
||||
unsigned interiorDofIdx = context.interiorScvIndex(bfIdx, timeIdx);
|
||||
|
||||
////////
|
||||
// advective fluxes of all components in all phases
|
||||
////////
|
||||
(*this) = Evaluation(0.0);
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
Evaluation density;
|
||||
if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
|
||||
if (focusDofIdx == interiorDofIdx)
|
||||
density = fluidState.density(phaseIdx);
|
||||
else
|
||||
density = Opm::getValue(fluidState.density(phaseIdx));
|
||||
}
|
||||
else if (focusDofIdx == interiorDofIdx)
|
||||
density = insideIntQuants.fluidState().density(phaseIdx);
|
||||
else
|
||||
density = Opm::getValue(insideIntQuants.fluidState().density(phaseIdx));
|
||||
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
Evaluation molarity;
|
||||
if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
|
||||
if (focusDofIdx == interiorDofIdx)
|
||||
molarity = fluidState.molarity(phaseIdx, compIdx);
|
||||
else
|
||||
molarity = Opm::getValue(fluidState.molarity(phaseIdx, compIdx));
|
||||
}
|
||||
else if (focusDofIdx == interiorDofIdx)
|
||||
molarity = insideIntQuants.fluidState().molarity(phaseIdx, compIdx);
|
||||
else
|
||||
molarity = Opm::getValue(insideIntQuants.fluidState().molarity(phaseIdx, compIdx));
|
||||
|
||||
// add advective flux of current component in current
|
||||
// phase
|
||||
(*this)[conti0EqIdx + compIdx] += extQuants.volumeFlux(phaseIdx)*molarity;
|
||||
}
|
||||
|
||||
if (enableEnergy) {
|
||||
Evaluation specificEnthalpy;
|
||||
if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
|
||||
if (focusDofIdx == interiorDofIdx)
|
||||
specificEnthalpy = fluidState.enthalpy(phaseIdx);
|
||||
else
|
||||
specificEnthalpy = Opm::getValue(fluidState.enthalpy(phaseIdx));
|
||||
}
|
||||
else if (focusDofIdx == interiorDofIdx)
|
||||
specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx);
|
||||
else
|
||||
specificEnthalpy = Opm::getValue(insideIntQuants.fluidState().enthalpy(phaseIdx));
|
||||
|
||||
Evaluation enthalpyRate = density*extQuants.volumeFlux(phaseIdx)*specificEnthalpy;
|
||||
EnergyModule::addToEnthalpyRate(*this, enthalpyRate);
|
||||
}
|
||||
}
|
||||
|
||||
// thermal conduction
|
||||
EnergyModule::addToEnthalpyRate(*this, EnergyModule::thermalConductionRate(extQuants));
|
||||
|
||||
#ifndef NDEBUG
|
||||
for (unsigned i = 0; i < numEq; ++i) {
|
||||
Opm::Valgrind::CheckDefined((*this)[i]);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleBoundaryRateVector::setInFlow
|
||||
*/
|
||||
template <class Context, class FluidState>
|
||||
void setInFlow(const Context& context,
|
||||
unsigned bfIdx,
|
||||
unsigned timeIdx,
|
||||
const FluidState& fluidState)
|
||||
{
|
||||
this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
|
||||
|
||||
// we only allow fluxes in the direction opposite to the outer unit normal
|
||||
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
|
||||
this->operator[](eqIdx) = Toolbox::min(0.0, this->operator[](eqIdx));
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleBoundaryRateVector::setOutFlow
|
||||
*/
|
||||
template <class Context, class FluidState>
|
||||
void setOutFlow(const Context& context,
|
||||
unsigned bfIdx,
|
||||
unsigned timeIdx,
|
||||
const FluidState& fluidState)
|
||||
{
|
||||
this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
|
||||
|
||||
// we only allow fluxes in the same direction as the outer unit normal
|
||||
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
|
||||
this->operator[](eqIdx) = Toolbox::max(0.0, this->operator[](eqIdx));
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleBoundaryRateVector::setNoFlow
|
||||
*/
|
||||
void setNoFlow()
|
||||
{ (*this) = Evaluation(0.0); }
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
91
opm/models/ncp/ncpextensivequantities.hh
Normal file
91
opm/models/ncp/ncpextensivequantities.hh
Normal file
@ -0,0 +1,91 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpExtensiveQuantities
|
||||
*/
|
||||
#ifndef EWOMS_NCP_EXTENSIVE_QUANTITIES_HH
|
||||
#define EWOMS_NCP_EXTENSIVE_QUANTITIES_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
#include <opm/models/common/diffusionmodule.hh>
|
||||
#include <opm/models/common/multiphasebaseextensivequantities.hh>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
* \ingroup ExtensiveQuantities
|
||||
*
|
||||
* \brief This template class represents the extensive quantities of the compositional
|
||||
* NCP model.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpExtensiveQuantities
|
||||
: public MultiPhaseBaseExtensiveQuantities<TypeTag>
|
||||
, public EnergyExtensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy)>
|
||||
, public DiffusionExtensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableDiffusion)>
|
||||
{
|
||||
typedef MultiPhaseBaseExtensiveQuantities<TypeTag> ParentType;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
|
||||
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
|
||||
typedef Opm::DiffusionExtensiveQuantities<TypeTag, enableDiffusion> DiffusionExtensiveQuantities;
|
||||
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
typedef Opm::EnergyExtensiveQuantities<TypeTag, enableEnergy> EnergyExtensiveQuantities;
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \copydoc MultiPhaseBaseExtensiveQuantities::update
|
||||
*/
|
||||
void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
|
||||
{
|
||||
ParentType::update(elemCtx, scvfIdx, timeIdx);
|
||||
DiffusionExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);
|
||||
EnergyExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc MultiPhaseBaseExtensiveQuantities::updateBoundary
|
||||
*/
|
||||
template <class Context, class FluidState>
|
||||
void updateBoundary(const Context& context,
|
||||
unsigned bfIdx,
|
||||
unsigned timeIdx,
|
||||
const FluidState& fluidState)
|
||||
{
|
||||
ParentType::updateBoundary(context, bfIdx, timeIdx, fluidState);
|
||||
DiffusionExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState);
|
||||
EnergyExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState);
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
104
opm/models/ncp/ncpindices.hh
Normal file
104
opm/models/ncp/ncpindices.hh
Normal file
@ -0,0 +1,104 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpIndices
|
||||
*/
|
||||
#ifndef EWOMS_NCP_INDICES_HH
|
||||
#define EWOMS_NCP_INDICES_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief The primary variable and equation indices for the
|
||||
* compositional multi-phase NCP model.
|
||||
*/
|
||||
template <class TypeTag, int PVOffset = 0>
|
||||
struct NcpIndices
|
||||
: public EnergyIndices<PVOffset
|
||||
+ GET_PROP_VALUE(TypeTag, NumComponents)
|
||||
+ GET_PROP_VALUE(TypeTag, NumPhases),
|
||||
GET_PROP_VALUE(TypeTag, EnableEnergy)>
|
||||
{
|
||||
private:
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
enum { numComponents = FluidSystem::numComponents };
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
|
||||
typedef Opm::EnergyIndices<PVOffset + numComponents + numPhases, enableEnergy> EnergyIndices;
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \brief The number of primary variables / equations.
|
||||
*/
|
||||
static const int numEq = numComponents + numPhases + EnergyIndices::numEq_;
|
||||
|
||||
/*!
|
||||
* \brief Index of the equation for the continuity of mass of the
|
||||
* first component.
|
||||
*
|
||||
* numComponents equations follow...
|
||||
*/
|
||||
static const int conti0EqIdx = PVOffset;
|
||||
|
||||
/*!
|
||||
* \brief Index of the first phase NCP equation.
|
||||
*
|
||||
* The index for the remaining phases are consecutive.
|
||||
*/
|
||||
static const int ncp0EqIdx = conti0EqIdx + numComponents;
|
||||
|
||||
/*!
|
||||
* \brief Index of the primary variable for the fugacity of the
|
||||
* first component.
|
||||
*
|
||||
* numComponents primary variables follow...
|
||||
*/
|
||||
static const int fugacity0Idx = PVOffset;
|
||||
|
||||
/*!
|
||||
* \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 saturation0Idx = fugacity0Idx + numComponents;
|
||||
|
||||
/*!
|
||||
* \brief Index of the first phase' pressure in a vector of
|
||||
* primary variables.
|
||||
*/
|
||||
static const int pressure0Idx = saturation0Idx + numPhases - 1;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
247
opm/models/ncp/ncpintensivequantities.hh
Normal file
247
opm/models/ncp/ncpintensivequantities.hh
Normal file
@ -0,0 +1,247 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpIntensiveQuantities
|
||||
*/
|
||||
#ifndef EWOMS_NCP_INTENSIVE_QUANTITIES_HH
|
||||
#define EWOMS_NCP_INTENSIVE_QUANTITIES_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
#include <opm/models/common/diffusionmodule.hh>
|
||||
|
||||
#include <opm/material/constraintsolvers/NcpFlash.hpp>
|
||||
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
|
||||
#include <opm/material/constraintsolvers/CompositionFromFugacities.hpp>
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
* \ingroup IntensiveQuantities
|
||||
*
|
||||
* \brief Contains the quantities which are are constant within a
|
||||
* finite volume in the compositional multi-phase NCP model.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpIntensiveQuantities
|
||||
: public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
|
||||
, public DiffusionIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableDiffusion) >
|
||||
, public EnergyIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy) >
|
||||
, public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
|
||||
{
|
||||
typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
|
||||
|
||||
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
|
||||
enum { fugacity0Idx = Indices::fugacity0Idx };
|
||||
enum { saturation0Idx = Indices::saturation0Idx };
|
||||
enum { pressure0Idx = Indices::pressure0Idx };
|
||||
enum { dimWorld = GridView::dimensionworld };
|
||||
|
||||
typedef Opm::CompositionFromFugacities<Scalar, FluidSystem, Evaluation>
|
||||
CompositionFromFugacitiesSolver;
|
||||
typedef Opm::CompositionalFluidState<Evaluation, FluidSystem,
|
||||
/*storeEnthalpy=*/enableEnergy> FluidState;
|
||||
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
|
||||
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
|
||||
typedef Opm::DiffusionIntensiveQuantities<TypeTag, enableDiffusion> DiffusionIntensiveQuantities;
|
||||
typedef Opm::EnergyIntensiveQuantities<TypeTag, enableEnergy> EnergyIntensiveQuantities;
|
||||
typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
|
||||
|
||||
public:
|
||||
NcpIntensiveQuantities()
|
||||
{}
|
||||
|
||||
NcpIntensiveQuantities(const NcpIntensiveQuantities& other) = default;
|
||||
|
||||
NcpIntensiveQuantities& operator=(const NcpIntensiveQuantities& other) = default;
|
||||
|
||||
/*!
|
||||
* \brief IntensiveQuantities::update
|
||||
*/
|
||||
void update(const ElementContext& elemCtx,
|
||||
unsigned dofIdx,
|
||||
unsigned timeIdx)
|
||||
{
|
||||
ParentType::update(elemCtx, dofIdx, timeIdx);
|
||||
ParentType::checkDefined();
|
||||
|
||||
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
|
||||
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
|
||||
|
||||
// set the phase saturations
|
||||
Evaluation sumSat = 0;
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
|
||||
const Evaluation& val = priVars.makeEvaluation(saturation0Idx + phaseIdx, timeIdx);
|
||||
fluidState_.setSaturation(phaseIdx, val);
|
||||
sumSat += val;
|
||||
}
|
||||
fluidState_.setSaturation(numPhases - 1, 1.0 - sumSat);
|
||||
Opm::Valgrind::CheckDefined(sumSat);
|
||||
|
||||
// set the fluid phase temperature
|
||||
EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
|
||||
|
||||
// retrieve capillary pressure parameters
|
||||
const auto& problem = elemCtx.problem();
|
||||
const MaterialLawParams& materialParams =
|
||||
problem.materialLawParams(elemCtx, dofIdx, timeIdx);
|
||||
// calculate capillary pressures
|
||||
Evaluation capPress[numPhases];
|
||||
MaterialLaw::capillaryPressures(capPress, materialParams, fluidState_);
|
||||
// add to the pressure of the first fluid phase
|
||||
const Evaluation& pressure0 = priVars.makeEvaluation(pressure0Idx, timeIdx);
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
fluidState_.setPressure(phaseIdx, pressure0 + (capPress[phaseIdx] - capPress[0]));
|
||||
|
||||
ComponentVector fug;
|
||||
// retrieve component fugacities
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
fug[compIdx] = priVars.makeEvaluation(fugacity0Idx + compIdx, timeIdx);
|
||||
|
||||
// calculate phase compositions
|
||||
const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx);
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
// initial guess
|
||||
if (hint) {
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
// use the hint for the initial mole fraction!
|
||||
const Evaluation& moleFracIJ = hint->fluidState().moleFraction(phaseIdx, compIdx);
|
||||
fluidState_.setMoleFraction(phaseIdx, compIdx, moleFracIJ);
|
||||
}
|
||||
}
|
||||
else // !hint
|
||||
CompositionFromFugacitiesSolver::guessInitial(fluidState_, phaseIdx, fug);
|
||||
|
||||
// calculate the phase composition from the component
|
||||
// fugacities
|
||||
CompositionFromFugacitiesSolver::solve(fluidState_, paramCache, phaseIdx, fug);
|
||||
}
|
||||
|
||||
// porosity
|
||||
porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
|
||||
Opm::Valgrind::CheckDefined(porosity_);
|
||||
|
||||
// relative permeabilities
|
||||
MaterialLaw::relativePermeabilities(relativePermeability_, materialParams, fluidState_);
|
||||
|
||||
// dynamic viscosities
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
// viscosities
|
||||
const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
|
||||
fluidState_.setViscosity(phaseIdx, mu);
|
||||
|
||||
mobility_[phaseIdx] = relativePermeability_[phaseIdx]/mu;
|
||||
}
|
||||
|
||||
// intrinsic permeability
|
||||
intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
|
||||
|
||||
// update the quantities specific for the velocity model
|
||||
FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
|
||||
|
||||
// energy related quantities
|
||||
EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
|
||||
|
||||
// update the diffusion specific quantities of the intensive quantities
|
||||
DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
|
||||
|
||||
checkDefined();
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief ImmiscibleIntensiveQuantities::fluidState
|
||||
*/
|
||||
const FluidState& fluidState() const
|
||||
{ return fluidState_; }
|
||||
|
||||
/*!
|
||||
* \brief ImmiscibleIntensiveQuantities::intrinsicPermeability
|
||||
*/
|
||||
const DimMatrix& intrinsicPermeability() const
|
||||
{ return intrinsicPerm_; }
|
||||
|
||||
/*!
|
||||
* \brief ImmiscibleIntensiveQuantities::relativePermeability
|
||||
*/
|
||||
const Evaluation& relativePermeability(unsigned phaseIdx) const
|
||||
{ return relativePermeability_[phaseIdx]; }
|
||||
|
||||
/*!
|
||||
* \brief ImmiscibleIntensiveQuantities::mobility
|
||||
*/
|
||||
const Evaluation& mobility(unsigned phaseIdx) const
|
||||
{ return mobility_[phaseIdx]; }
|
||||
|
||||
/*!
|
||||
* \brief ImmiscibleIntensiveQuantities::porosity
|
||||
*/
|
||||
const Evaluation& porosity() const
|
||||
{ return porosity_; }
|
||||
|
||||
/*!
|
||||
* \brief IntensiveQuantities::checkDefined
|
||||
*/
|
||||
void checkDefined() const
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
ParentType::checkDefined();
|
||||
|
||||
Opm::Valgrind::CheckDefined(porosity_);
|
||||
Opm::Valgrind::CheckDefined(relativePermeability_);
|
||||
|
||||
fluidState_.checkDefined();
|
||||
#endif
|
||||
}
|
||||
|
||||
private:
|
||||
DimMatrix intrinsicPerm_;
|
||||
FluidState fluidState_;
|
||||
Evaluation porosity_;
|
||||
Evaluation relativePermeability_[numPhases];
|
||||
Evaluation mobility_[numPhases];
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
261
opm/models/ncp/ncplocalresidual.hh
Normal file
261
opm/models/ncp/ncplocalresidual.hh
Normal file
@ -0,0 +1,261 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpLocalResidual
|
||||
*/
|
||||
#ifndef EWOMS_NCP_LOCAL_RESIDUAL_HH
|
||||
#define EWOMS_NCP_LOCAL_RESIDUAL_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
|
||||
#include <opm/models/common/diffusionmodule.hh>
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief Details needed to calculate the local residual in the
|
||||
* compositional multi-phase NCP-model .
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpLocalResidual : public GET_PROP_TYPE(TypeTag, DiscLocalResidual)
|
||||
{
|
||||
typedef typename GET_PROP_TYPE(TypeTag, DiscLocalResidual) ParentType;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
|
||||
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
|
||||
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||
enum { ncp0EqIdx = Indices::ncp0EqIdx };
|
||||
enum { conti0EqIdx = Indices::conti0EqIdx };
|
||||
|
||||
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
|
||||
typedef Opm::DiffusionModule<TypeTag, enableDiffusion> DiffusionModule;
|
||||
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
|
||||
|
||||
typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
|
||||
typedef Dune::BlockVector<EvalEqVector> ElemEvalEqVector;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \copydoc ImmiscibleLocalResidual::addPhaseStorage
|
||||
*/
|
||||
template <class LhsEval>
|
||||
void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
|
||||
const ElementContext& elemCtx,
|
||||
unsigned dofIdx,
|
||||
unsigned timeIdx,
|
||||
unsigned phaseIdx) const
|
||||
{
|
||||
const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
|
||||
const auto& fluidState = intQuants.fluidState();
|
||||
|
||||
// compute storage term of all components within all phases
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
unsigned eqIdx = conti0EqIdx + compIdx;
|
||||
storage[eqIdx] +=
|
||||
Toolbox::template decay<LhsEval>(fluidState.molarity(phaseIdx, compIdx))
|
||||
* Toolbox::template decay<LhsEval>(fluidState.saturation(phaseIdx))
|
||||
* Toolbox::template decay<LhsEval>(intQuants.porosity());
|
||||
}
|
||||
|
||||
EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleLocalResidual::computeStorage
|
||||
*/
|
||||
template <class LhsEval>
|
||||
void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
|
||||
const ElementContext& elemCtx,
|
||||
unsigned dofIdx,
|
||||
unsigned timeIdx) const
|
||||
{
|
||||
storage = 0;
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
addPhaseStorage(storage, elemCtx, dofIdx, timeIdx, phaseIdx);
|
||||
|
||||
EnergyModule::addSolidEnergyStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleLocalResidual::computeFlux
|
||||
*/
|
||||
void computeFlux(RateVector& flux,
|
||||
const ElementContext& elemCtx,
|
||||
unsigned scvfIdx,
|
||||
unsigned timeIdx) const
|
||||
{
|
||||
flux = 0.0;
|
||||
addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
|
||||
Opm::Valgrind::CheckDefined(flux);
|
||||
|
||||
addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
|
||||
Opm::Valgrind::CheckDefined(flux);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleLocalResidual::addAdvectiveFlux
|
||||
*/
|
||||
void addAdvectiveFlux(RateVector& flux,
|
||||
const ElementContext& elemCtx,
|
||||
unsigned scvfIdx,
|
||||
unsigned timeIdx) const
|
||||
{
|
||||
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
|
||||
|
||||
unsigned focusDofIdx = elemCtx.focusDofIndex();
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
// data attached to upstream and the downstream DOFs
|
||||
// of the current phase
|
||||
unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
|
||||
const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
|
||||
|
||||
// this is a bit hacky because it is specific to the element-centered
|
||||
// finite volume scheme. (N.B. that if finite differences are used to
|
||||
// linearize the system of equations, it does not matter.)
|
||||
if (upIdx == focusDofIdx) {
|
||||
Evaluation tmp =
|
||||
up.fluidState().molarDensity(phaseIdx)
|
||||
* extQuants.volumeFlux(phaseIdx);
|
||||
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
flux[conti0EqIdx + compIdx] +=
|
||||
tmp*up.fluidState().moleFraction(phaseIdx, compIdx);
|
||||
}
|
||||
}
|
||||
else {
|
||||
Evaluation tmp =
|
||||
Toolbox::value(up.fluidState().molarDensity(phaseIdx))
|
||||
* extQuants.volumeFlux(phaseIdx);
|
||||
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
flux[conti0EqIdx + compIdx] +=
|
||||
tmp*Toolbox::value(up.fluidState().moleFraction(phaseIdx, compIdx));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleLocalResidual::addDiffusiveFlux
|
||||
*/
|
||||
void addDiffusiveFlux(RateVector& flux,
|
||||
const ElementContext& elemCtx,
|
||||
unsigned scvfIdx,
|
||||
unsigned timeIdx) const
|
||||
{
|
||||
DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
|
||||
EnergyModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseLocalResidual::computeSource
|
||||
*
|
||||
* By default, this method only asks the problem to specify a
|
||||
* source term.
|
||||
*/
|
||||
void computeSource(RateVector& source,
|
||||
const ElementContext& elemCtx,
|
||||
unsigned dofIdx,
|
||||
unsigned timeIdx) const
|
||||
{
|
||||
Opm::Valgrind::SetUndefined(source);
|
||||
elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
|
||||
Opm::Valgrind::CheckDefined(source);
|
||||
|
||||
// evaluate the NCPs (i.e., the "phase presence" equations)
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
source[ncp0EqIdx + phaseIdx] =
|
||||
phaseNcp(elemCtx, dofIdx, timeIdx, phaseIdx);
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the value of the NCP-function for a phase.
|
||||
*/
|
||||
template <class LhsEval = Evaluation>
|
||||
LhsEval phaseNcp(const ElementContext& elemCtx,
|
||||
unsigned dofIdx,
|
||||
unsigned timeIdx,
|
||||
unsigned phaseIdx) const
|
||||
{
|
||||
const auto& fluidState = elemCtx.intensiveQuantities(dofIdx, timeIdx).fluidState();
|
||||
typedef typename std::remove_const<typename std::remove_reference<decltype(fluidState)>::type>::type FluidState;
|
||||
|
||||
typedef Opm::MathToolbox<LhsEval> LhsToolbox;
|
||||
|
||||
const LhsEval& a = phaseNotPresentIneq_<FluidState, LhsEval>(fluidState, phaseIdx);
|
||||
const LhsEval& b = phasePresentIneq_<FluidState, LhsEval>(fluidState, phaseIdx);
|
||||
return LhsToolbox::min(a, b);
|
||||
}
|
||||
|
||||
private:
|
||||
/*!
|
||||
* \brief Returns the value of the inequality where a phase is
|
||||
* present.
|
||||
*/
|
||||
template <class FluidState, class LhsEval>
|
||||
LhsEval phasePresentIneq_(const FluidState& fluidState, unsigned phaseIdx) const
|
||||
{
|
||||
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
return FsToolbox::template decay<LhsEval>(fluidState.saturation(phaseIdx));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the value of the inequality where a phase is not
|
||||
* present.
|
||||
*/
|
||||
template <class FluidState, class LhsEval>
|
||||
LhsEval phaseNotPresentIneq_(const FluidState& fluidState, unsigned phaseIdx) const
|
||||
{
|
||||
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
// difference of sum of mole fractions in the phase from 100%
|
||||
LhsEval a = 1.0;
|
||||
for (unsigned i = 0; i < numComponents; ++i)
|
||||
a -= FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
|
||||
return a;
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
455
opm/models/ncp/ncpmodel.hh
Normal file
455
opm/models/ncp/ncpmodel.hh
Normal file
@ -0,0 +1,455 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpModel
|
||||
*/
|
||||
#ifndef EWOMS_NCP_MODEL_HH
|
||||
#define EWOMS_NCP_MODEL_HH
|
||||
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
#include "ncplocalresidual.hh"
|
||||
#include "ncpextensivequantities.hh"
|
||||
#include "ncpprimaryvariables.hh"
|
||||
#include "ncpboundaryratevector.hh"
|
||||
#include "ncpratevector.hh"
|
||||
#include "ncpintensivequantities.hh"
|
||||
#include "ncpnewtonmethod.hh"
|
||||
#include "ncpindices.hh"
|
||||
|
||||
#include <opm/models/common/multiphasebasemodel.hh>
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
#include <opm/models/common/diffusionmodule.hh>
|
||||
#include <opm/models/io/vtkcompositionmodule.hh>
|
||||
#include <opm/models/io/vtkenergymodule.hh>
|
||||
#include <opm/models/io/vtkdiffusionmodule.hh>
|
||||
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
#include <opm/material/common/Exceptions.hpp>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <array>
|
||||
|
||||
namespace Opm {
|
||||
template <class TypeTag>
|
||||
class NcpModel;
|
||||
}
|
||||
|
||||
BEGIN_PROPERTIES
|
||||
|
||||
/*!
|
||||
* \brief Define the type tag for the compositional NCP model.
|
||||
*/
|
||||
NEW_TYPE_TAG(NcpModel, INHERITS_FROM(MultiPhaseBaseModel,
|
||||
VtkComposition,
|
||||
VtkEnergy,
|
||||
VtkDiffusion));
|
||||
|
||||
//! Use the Ncp local jacobian operator for the compositional NCP model
|
||||
SET_TYPE_PROP(NcpModel,
|
||||
LocalResidual,
|
||||
Opm::NcpLocalResidual<TypeTag>);
|
||||
|
||||
//! Use the Ncp specific newton method for the compositional NCP model
|
||||
SET_TYPE_PROP(NcpModel, NewtonMethod, Opm::NcpNewtonMethod<TypeTag>);
|
||||
|
||||
//! the Model property
|
||||
SET_TYPE_PROP(NcpModel, Model, Opm::NcpModel<TypeTag>);
|
||||
|
||||
//! The type of the base base class for actual problems
|
||||
SET_TYPE_PROP(NcpModel, BaseProblem, Opm::MultiPhaseBaseProblem<TypeTag>);
|
||||
|
||||
//! Disable the energy equation by default
|
||||
SET_BOOL_PROP(NcpModel, EnableEnergy, false);
|
||||
|
||||
//! disable diffusion by default
|
||||
SET_BOOL_PROP(NcpModel, EnableDiffusion, false);
|
||||
|
||||
//! the RateVector property
|
||||
SET_TYPE_PROP(NcpModel, RateVector, Opm::NcpRateVector<TypeTag>);
|
||||
|
||||
//! the BoundaryRateVector property
|
||||
SET_TYPE_PROP(NcpModel, BoundaryRateVector, Opm::NcpBoundaryRateVector<TypeTag>);
|
||||
|
||||
//! the PrimaryVariables property
|
||||
SET_TYPE_PROP(NcpModel, PrimaryVariables, Opm::NcpPrimaryVariables<TypeTag>);
|
||||
|
||||
//! the IntensiveQuantities property
|
||||
SET_TYPE_PROP(NcpModel, IntensiveQuantities, Opm::NcpIntensiveQuantities<TypeTag>);
|
||||
|
||||
//! the ExtensiveQuantities property
|
||||
SET_TYPE_PROP(NcpModel, ExtensiveQuantities, Opm::NcpExtensiveQuantities<TypeTag>);
|
||||
|
||||
//! The indices required by the compositional NCP model
|
||||
SET_TYPE_PROP(NcpModel, Indices, Opm::NcpIndices<TypeTag, 0>);
|
||||
|
||||
//! The unmodified weight for the pressure primary variable
|
||||
SET_SCALAR_PROP(NcpModel, NcpPressureBaseWeight, 1.0);
|
||||
//! The weight for the saturation primary variables
|
||||
SET_SCALAR_PROP(NcpModel, NcpSaturationsBaseWeight, 1.0);
|
||||
//! The unmodified weight for the fugacity primary variables
|
||||
SET_SCALAR_PROP(NcpModel, NcpFugacitiesBaseWeight, 1.0e-6);
|
||||
|
||||
END_PROPERTIES
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief A compositional multi-phase model based on non-linear
|
||||
* complementarity functions.
|
||||
*
|
||||
* This model implements a \f$M\f$-phase flow of a fluid mixture
|
||||
* composed of \f$N\f$ chemical species. The phases are denoted by
|
||||
* lower index \f$\alpha \in \{ 1, \dots, M \}\f$. All fluid phases
|
||||
* are mixtures of \f$N \geq M - 1\f$ chemical species which are
|
||||
* denoted by the upper index \f$\kappa \in \{ 1, \dots, N \} \f$.
|
||||
*
|
||||
*
|
||||
* By default, the standard multi-phase Darcy approach is used to determine
|
||||
* the velocity, i.e.
|
||||
* \f[
|
||||
* \mathbf{v}_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
|
||||
* \left(\mathbf{grad}\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right) \;,
|
||||
* \f]
|
||||
* although the actual approach which is used can be specified via the
|
||||
* \c FluxModule property. For example, the velocity model can by
|
||||
* changed to the Forchheimer approach by
|
||||
* \code
|
||||
* SET_TYPE_PROP(MyProblemTypeTag, FluxModule, Opm::ForchheimerFluxModule<TypeTag>);
|
||||
* \endcode
|
||||
*
|
||||
* The core of the model is the conservation mass of each component by
|
||||
* means of the equation
|
||||
* \f[
|
||||
* \sum_\alpha \frac{\partial\;\phi c_\alpha^\kappa S_\alpha }{\partial t}
|
||||
* - \sum_\alpha \mathrm{div} \left\{ c_\alpha^\kappa \mathbf{v}_\alpha \right\}
|
||||
* - q^\kappa = 0 \;.
|
||||
* \f]
|
||||
*
|
||||
* For the missing \f$M\f$ model assumptions, the model uses
|
||||
* non-linear complementarity functions. These are based on the
|
||||
* observation that if a fluid phase is not present, the sum of the
|
||||
* mole fractions of this fluid phase is smaller than \f$1\f$, i.e.
|
||||
* \f[ \forall \alpha: S_\alpha = 0 \implies \sum_\kappa
|
||||
* x_\alpha^\kappa \leq 1 \f]
|
||||
*
|
||||
* Also, if a fluid phase may be present at a given spatial location
|
||||
* its saturation must be non-negative:
|
||||
* \f[ \forall \alpha: \sum_\kappa x_\alpha^\kappa = 1 \implies S_\alpha \geq 0
|
||||
*\f]
|
||||
*
|
||||
* Since at any given spatial location, a phase is always either
|
||||
* present or not present, one of the strict equalities on the
|
||||
* right hand side is always true, i.e.
|
||||
* \f[
|
||||
* \forall \alpha: S_\alpha \left( \sum_\kappa x_\alpha^\kappa - 1 \right) = 0
|
||||
* \f]
|
||||
* always holds.
|
||||
*
|
||||
* These three equations constitute a non-linear complementarity
|
||||
* problem, which can be solved using so-called non-linear
|
||||
* complementarity functions \f$\Phi(a, b)\f$. Such functions have the property
|
||||
* \f[\Phi(a,b) = 0 \iff a \geq0 \land b \geq0 \land a \cdot b = 0 \f]
|
||||
*
|
||||
* Several non-linear complementarity functions have been suggested,
|
||||
* e.g. the Fischer-Burmeister function
|
||||
* \f[ \Phi(a,b) = a + b - \sqrt{a^2 + b^2} \;. \f]
|
||||
* This model uses
|
||||
* \f[ \Phi(a,b) = \min \{a, b \}\;, \f]
|
||||
* because of its piecewise linearity.
|
||||
*
|
||||
* The model assumes local thermodynamic equilibrium and uses the
|
||||
* following primary variables:
|
||||
* - The pressure of the first phase \f$p_1\f$
|
||||
* - The component fugacities \f$f^1, \dots, f^{N}\f$
|
||||
* - The saturations of the first \f$M-1\f$ phases \f$S_1, \dots, S_{M-1}\f$
|
||||
* - Temperature \f$T\f$ if the energy equation is enabled
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpModel
|
||||
: public MultiPhaseBaseModel<TypeTag>
|
||||
{
|
||||
typedef MultiPhaseBaseModel<TypeTag> ParentType;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
enum { numComponents = FluidSystem::numComponents };
|
||||
enum { fugacity0Idx = Indices::fugacity0Idx };
|
||||
enum { pressure0Idx = Indices::pressure0Idx };
|
||||
enum { saturation0Idx = Indices::saturation0Idx };
|
||||
enum { conti0EqIdx = Indices::conti0EqIdx };
|
||||
enum { ncp0EqIdx = Indices::ncp0EqIdx };
|
||||
enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
|
||||
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
|
||||
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
|
||||
typedef Opm::DiffusionModule<TypeTag, enableDiffusion> DiffusionModule;
|
||||
|
||||
public:
|
||||
NcpModel(Simulator& simulator)
|
||||
: ParentType(simulator)
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \brief Register all run-time parameters for the immiscible model.
|
||||
*/
|
||||
static void registerParameters()
|
||||
{
|
||||
ParentType::registerParameters();
|
||||
|
||||
DiffusionModule::registerParameters();
|
||||
EnergyModule::registerParameters();
|
||||
|
||||
// register runtime parameters of the VTK output modules
|
||||
Opm::VtkCompositionModule<TypeTag>::registerParameters();
|
||||
|
||||
if (enableDiffusion)
|
||||
Opm::VtkDiffusionModule<TypeTag>::registerParameters();
|
||||
|
||||
if (enableEnergy)
|
||||
Opm::VtkEnergyModule<TypeTag>::registerParameters();
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::finishInit()
|
||||
*/
|
||||
void finishInit()
|
||||
{
|
||||
ParentType::finishInit();
|
||||
|
||||
minActivityCoeff_.resize(this->numGridDof());
|
||||
std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
|
||||
}
|
||||
|
||||
void adaptGrid()
|
||||
{
|
||||
ParentType::adaptGrid();
|
||||
minActivityCoeff_.resize(this->numGridDof());
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::name
|
||||
*/
|
||||
static std::string name()
|
||||
{ return "ncp"; }
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::primaryVarName
|
||||
*/
|
||||
std::string primaryVarName(unsigned pvIdx) const
|
||||
{
|
||||
std::string s;
|
||||
if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
|
||||
return s;
|
||||
|
||||
std::ostringstream oss;
|
||||
if (pvIdx == pressure0Idx)
|
||||
oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
|
||||
else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1))
|
||||
oss << "saturation_" << FluidSystem::phaseName(/*phaseIdx=*/pvIdx - saturation0Idx);
|
||||
else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents)
|
||||
oss << "fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
|
||||
else
|
||||
assert(false);
|
||||
|
||||
return oss.str();
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::eqName
|
||||
*/
|
||||
std::string eqName(unsigned eqIdx) const
|
||||
{
|
||||
std::string s;
|
||||
if (!(s = EnergyModule::eqName(eqIdx)).empty())
|
||||
return s;
|
||||
|
||||
std::ostringstream oss;
|
||||
if (conti0EqIdx <= eqIdx && eqIdx < conti0EqIdx + numComponents)
|
||||
oss << "continuity^" << FluidSystem::componentName(eqIdx - conti0EqIdx);
|
||||
else if (ncp0EqIdx <= eqIdx && eqIdx < ncp0EqIdx + numPhases)
|
||||
oss << "ncp_" << FluidSystem::phaseName(/*phaseIdx=*/eqIdx - ncp0EqIdx);
|
||||
else
|
||||
assert(false);
|
||||
|
||||
return oss.str();
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::updateBegin
|
||||
*/
|
||||
void updateBegin()
|
||||
{
|
||||
ParentType::updateBegin();
|
||||
|
||||
// find the a reference pressure. The first degree of freedom
|
||||
// might correspond to non-interior entities which would lead
|
||||
// to an undefined value, so we have to iterate...
|
||||
for (unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++ dofIdx) {
|
||||
if (this->isLocalDof(dofIdx)) {
|
||||
referencePressure_ =
|
||||
this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::updatePVWeights
|
||||
*/
|
||||
void updatePVWeights(const ElementContext& elemCtx) const
|
||||
{
|
||||
for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
|
||||
unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
|
||||
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
minActivityCoeff_[globalIdx][compIdx] = 1e100;
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
|
||||
|
||||
minActivityCoeff_[globalIdx][compIdx] =
|
||||
std::min(minActivityCoeff_[globalIdx][compIdx],
|
||||
Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx))
|
||||
* Toolbox::value(fs.pressure(phaseIdx)));
|
||||
Opm::Valgrind::CheckDefined(minActivityCoeff_[globalIdx][compIdx]);
|
||||
}
|
||||
if (minActivityCoeff_[globalIdx][compIdx] <= 0)
|
||||
throw Opm::NumericalIssue("The minumum activity coefficient for component "+std::to_string(compIdx)
|
||||
+" on DOF "+std::to_string(globalIdx)+" is negative or zero!");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::primaryVarWeight
|
||||
*/
|
||||
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
|
||||
{
|
||||
Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
|
||||
Scalar result;
|
||||
if (tmp > 0)
|
||||
// energy related quantity
|
||||
result = tmp;
|
||||
else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
|
||||
// component fugacity
|
||||
unsigned compIdx = pvIdx - fugacity0Idx;
|
||||
assert(0 <= compIdx && compIdx <= numComponents);
|
||||
|
||||
Opm::Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
|
||||
static const Scalar fugacityBaseWeight =
|
||||
GET_PROP_VALUE(TypeTag, NcpFugacitiesBaseWeight);
|
||||
result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
|
||||
}
|
||||
else if (Indices::pressure0Idx == pvIdx) {
|
||||
static const Scalar pressureBaseWeight = GET_PROP_VALUE(TypeTag, NcpPressureBaseWeight);
|
||||
result = pressureBaseWeight / referencePressure_;
|
||||
}
|
||||
else {
|
||||
#ifndef NDEBUG
|
||||
unsigned phaseIdx = pvIdx - saturation0Idx;
|
||||
assert(0 <= phaseIdx && phaseIdx < numPhases - 1);
|
||||
#endif
|
||||
|
||||
// saturation
|
||||
static const Scalar saturationsBaseWeight =
|
||||
GET_PROP_VALUE(TypeTag, NcpSaturationsBaseWeight);
|
||||
result = saturationsBaseWeight;
|
||||
}
|
||||
|
||||
assert(std::isfinite(result));
|
||||
assert(result > 0);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseDiscretization::eqWeight
|
||||
*/
|
||||
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
|
||||
{
|
||||
Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
|
||||
if (tmp > 0)
|
||||
// an energy related equation
|
||||
return tmp;
|
||||
// an NCP
|
||||
else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
|
||||
return 1.0;
|
||||
|
||||
// a mass conservation equation
|
||||
unsigned compIdx = eqIdx - Indices::conti0EqIdx;
|
||||
assert(0 <= compIdx && compIdx <= numComponents);
|
||||
|
||||
// make all kg equal
|
||||
return FluidSystem::molarMass(compIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the smallest activity coefficient of a component for the
|
||||
* most current solution at a vertex.
|
||||
*
|
||||
* \param globalDofIdx The global index of the vertex (i.e. finite volume) of interest.
|
||||
* \param compIdx The index of the component of interest.
|
||||
*/
|
||||
Scalar minActivityCoeff(unsigned globalDofIdx, unsigned compIdx) const
|
||||
{ return minActivityCoeff_[globalDofIdx][compIdx]; }
|
||||
|
||||
/*!
|
||||
* \internal
|
||||
*/
|
||||
void registerOutputModules_()
|
||||
{
|
||||
ParentType::registerOutputModules_();
|
||||
|
||||
this->addOutputModule(new Opm::VtkCompositionModule<TypeTag>(this->simulator_));
|
||||
if (enableDiffusion)
|
||||
this->addOutputModule(new Opm::VtkDiffusionModule<TypeTag>(this->simulator_));
|
||||
if (enableEnergy)
|
||||
this->addOutputModule(new Opm::VtkEnergyModule<TypeTag>(this->simulator_));
|
||||
}
|
||||
|
||||
mutable Scalar referencePressure_;
|
||||
mutable std::vector<ComponentVector> minActivityCoeff_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
211
opm/models/ncp/ncpnewtonmethod.hh
Normal file
211
opm/models/ncp/ncpnewtonmethod.hh
Normal file
@ -0,0 +1,211 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpNewtonMethod
|
||||
*/
|
||||
#ifndef EWOMS_NCP_NEWTON_METHOD_HH
|
||||
#define EWOMS_NCP_NEWTON_METHOD_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
|
||||
#include <opm/material/common/Unused.hpp>
|
||||
#include <opm/material/common/Exceptions.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief A Newton solver specific to the NCP model.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpNewtonMethod : public GET_PROP_TYPE(TypeTag, DiscNewtonMethod)
|
||||
{
|
||||
typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
|
||||
|
||||
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
|
||||
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||
enum { fugacity0Idx = Indices::fugacity0Idx };
|
||||
enum { saturation0Idx = Indices::saturation0Idx };
|
||||
enum { pressure0Idx = Indices::pressure0Idx };
|
||||
enum { ncp0EqIdx = Indices::ncp0EqIdx };
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& )
|
||||
*/
|
||||
NcpNewtonMethod(Simulator& simulator) : ParentType(simulator)
|
||||
{}
|
||||
|
||||
protected:
|
||||
friend ParentType;
|
||||
friend NewtonMethod<TypeTag>;
|
||||
|
||||
void preSolve_(const SolutionVector& currentSolution OPM_UNUSED,
|
||||
const GlobalEqVector& currentResidual)
|
||||
{
|
||||
const auto& constraintsMap = this->model().linearizer().constraintsMap();
|
||||
this->lastError_ = this->error_;
|
||||
|
||||
// calculate the error as the maximum weighted tolerance of
|
||||
// the solution's residual
|
||||
this->error_ = 0;
|
||||
for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
|
||||
// do not consider auxiliary DOFs for the error
|
||||
if (dofIdx >= this->model().numGridDof() || this->model().dofTotalVolume(dofIdx) <= 0.0)
|
||||
continue;
|
||||
|
||||
// also do not consider DOFs which are constraint
|
||||
if (this->enableConstraints_()) {
|
||||
if (constraintsMap.count(dofIdx) > 0)
|
||||
continue;
|
||||
}
|
||||
|
||||
const auto& r = currentResidual[dofIdx];
|
||||
for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) {
|
||||
if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
|
||||
continue;
|
||||
this->error_ =
|
||||
std::max(std::abs(r[eqIdx]*this->model().eqWeight(dofIdx, eqIdx)),
|
||||
this->error_);
|
||||
}
|
||||
}
|
||||
|
||||
// take the other processes into account
|
||||
this->error_ = this->comm_.max(this->error_);
|
||||
|
||||
// make sure that the error never grows beyond the maximum
|
||||
// allowed one
|
||||
if (this->error_ > EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError))
|
||||
throw Opm::NumericalIssue("Newton: Error "+std::to_string(double(this->error_))+
|
||||
+" is larger than maximum allowed error of "
|
||||
+std::to_string(double(EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError))));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc FvBaseNewtonMethod::updatePrimaryVariables_
|
||||
*/
|
||||
void updatePrimaryVariables_(unsigned globalDofIdx,
|
||||
PrimaryVariables& nextValue,
|
||||
const PrimaryVariables& currentValue,
|
||||
const EqVector& update,
|
||||
const EqVector& currentResidual OPM_UNUSED)
|
||||
{
|
||||
// normal Newton-Raphson update
|
||||
nextValue = currentValue;
|
||||
nextValue -= update;
|
||||
|
||||
////
|
||||
// put crash barriers along the update path
|
||||
////
|
||||
|
||||
// saturations: limit the change of any saturation to at most 20%
|
||||
Scalar sumSatDelta = 0.0;
|
||||
Scalar maxSatDelta = 0.0;
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
|
||||
maxSatDelta = std::max(std::abs(update[saturation0Idx + phaseIdx]),
|
||||
maxSatDelta);
|
||||
sumSatDelta += update[saturation0Idx + phaseIdx];
|
||||
}
|
||||
maxSatDelta = std::max(std::abs(- sumSatDelta), maxSatDelta);
|
||||
|
||||
if (maxSatDelta > 0.2) {
|
||||
Scalar alpha = 0.2/maxSatDelta;
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
|
||||
nextValue[saturation0Idx + phaseIdx] =
|
||||
currentValue[saturation0Idx + phaseIdx]
|
||||
- alpha*update[saturation0Idx + phaseIdx];
|
||||
}
|
||||
}
|
||||
|
||||
// limit pressure reference change to 20% of the total value per iteration
|
||||
clampValue_(nextValue[pressure0Idx],
|
||||
currentValue[pressure0Idx]*0.8,
|
||||
currentValue[pressure0Idx]*1.2);
|
||||
|
||||
// fugacities
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
Scalar& val = nextValue[fugacity0Idx + compIdx];
|
||||
Scalar oldVal = currentValue[fugacity0Idx + compIdx];
|
||||
|
||||
// get the minimum activity coefficient for the component (i.e., the activity
|
||||
// coefficient of the phase for which the component has the highest affinity)
|
||||
Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
|
||||
// Make sure that the activity coefficient does not get too small.
|
||||
minPhi = std::max(0.001*currentValue[pressure0Idx], minPhi);
|
||||
|
||||
// allow the mole fraction of the component to change at most 70% in any
|
||||
// phase (assuming composition independent fugacity coefficients).
|
||||
Scalar maxDelta = 0.7 * minPhi;
|
||||
clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
|
||||
|
||||
// make sure that fugacities do not become negative
|
||||
val = std::max(val, 0.0);
|
||||
}
|
||||
|
||||
// do not become grossly unphysical in a single iteration for the first few
|
||||
// iterations of a time step
|
||||
if (this->numIterations_ < 3) {
|
||||
// fugacities
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
Scalar& val = nextValue[fugacity0Idx + compIdx];
|
||||
Scalar oldVal = currentValue[fugacity0Idx + compIdx];
|
||||
Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
|
||||
if (oldVal < 1.0*minPhi && val > 1.0*minPhi)
|
||||
val = 1.0*minPhi;
|
||||
else if (oldVal > 0.0 && val < 0.0)
|
||||
val = 0.0;
|
||||
}
|
||||
|
||||
// saturations
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
|
||||
Scalar& val = nextValue[saturation0Idx + phaseIdx];
|
||||
Scalar oldVal = currentValue[saturation0Idx + phaseIdx];
|
||||
if (oldVal < 1.0 && val > 1.0)
|
||||
val = 1.0;
|
||||
else if (oldVal > 0.0 && val < 0.0)
|
||||
val = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const
|
||||
{ val = std::max(minVal, std::min(val, maxVal)); }
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
192
opm/models/ncp/ncpprimaryvariables.hh
Normal file
192
opm/models/ncp/ncpprimaryvariables.hh
Normal file
@ -0,0 +1,192 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpPrimaryVariables
|
||||
*/
|
||||
#ifndef EWOMS_NCP_PRIMARY_VARIABLES_HH
|
||||
#define EWOMS_NCP_PRIMARY_VARIABLES_HH
|
||||
|
||||
#include "ncpproperties.hh"
|
||||
|
||||
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
|
||||
#include <opm/models/common/energymodule.hh>
|
||||
|
||||
#include <opm/material/constraintsolvers/NcpFlash.hpp>
|
||||
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief Represents the primary variables used by the compositional
|
||||
* multi-phase NCP model.
|
||||
*
|
||||
* This class is basically a Dune::FieldVector which can retrieve its
|
||||
* contents from an aribitatry fluid state.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
|
||||
{
|
||||
typedef FvBasePrimaryVariables<TypeTag> ParentType;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
enum { pressure0Idx = Indices::pressure0Idx };
|
||||
enum { saturation0Idx = Indices::saturation0Idx };
|
||||
enum { fugacity0Idx = Indices::fugacity0Idx };
|
||||
|
||||
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
|
||||
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
|
||||
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
|
||||
|
||||
typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
public:
|
||||
NcpPrimaryVariables() : ParentType()
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(Scalar)
|
||||
*/
|
||||
NcpPrimaryVariables(Scalar value) : ParentType(value)
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(const
|
||||
* ImmisciblePrimaryVariables& )
|
||||
*/
|
||||
NcpPrimaryVariables(const NcpPrimaryVariables& value) = default;
|
||||
NcpPrimaryVariables& operator=(const NcpPrimaryVariables& value) = default;
|
||||
|
||||
/*!
|
||||
* \copydoc ImmisciblePrimaryVariables::assignMassConservative
|
||||
*/
|
||||
template <class FluidState>
|
||||
void assignMassConservative(const FluidState& fluidState,
|
||||
const MaterialLawParams& matParams,
|
||||
bool isInEquilibrium = false)
|
||||
{
|
||||
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
#ifndef NDEBUG
|
||||
// make sure the temperature is the same in all fluid phases
|
||||
for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
|
||||
assert(fluidState.temperature(0) == fluidState.temperature(phaseIdx));
|
||||
}
|
||||
#endif // NDEBUG
|
||||
|
||||
// for the equilibrium case, we don't need complicated
|
||||
// computations.
|
||||
if (isInEquilibrium) {
|
||||
assignNaive(fluidState);
|
||||
return;
|
||||
}
|
||||
|
||||
// use a flash calculation to calculate a fluid state in
|
||||
// thermodynamic equilibrium
|
||||
typename FluidSystem::template ParameterCache<Scalar> paramCache;
|
||||
Opm::CompositionalFluidState<Scalar, FluidSystem> fsFlash;
|
||||
|
||||
// use the externally given fluid state as initial value for
|
||||
// the flash calculation
|
||||
fsFlash.assign(fluidState);
|
||||
|
||||
// calculate the phase densities
|
||||
paramCache.updateAll(fsFlash);
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
Scalar rho = FluidSystem::density(fsFlash, paramCache, phaseIdx);
|
||||
fsFlash.setDensity(phaseIdx, rho);
|
||||
}
|
||||
|
||||
// calculate the "global molarities"
|
||||
ComponentVector globalMolarities(0.0);
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
globalMolarities[compIdx] +=
|
||||
FsToolbox::value(fsFlash.saturation(phaseIdx))
|
||||
* FsToolbox::value(fsFlash.molarity(phaseIdx, compIdx));
|
||||
}
|
||||
}
|
||||
|
||||
// run the flash calculation
|
||||
NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
|
||||
|
||||
// use the result to assign the primary variables
|
||||
assignNaive(fsFlash);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmisciblePrimaryVariables::assignNaive
|
||||
*/
|
||||
template <class FluidState>
|
||||
void assignNaive(const FluidState& fluidState, unsigned refPhaseIdx = 0)
|
||||
{
|
||||
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
// assign the phase temperatures. this is out-sourced to
|
||||
// the energy module
|
||||
EnergyModule::setPriVarTemperatures(*this, fluidState);
|
||||
|
||||
// assign fugacities.
|
||||
typename FluidSystem::template ParameterCache<Scalar> paramCache;
|
||||
paramCache.updatePhase(fluidState, refPhaseIdx);
|
||||
Scalar pRef = FsToolbox::value(fluidState.pressure(refPhaseIdx));
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
// we always compute the fugacities because they are quite exotic quantities
|
||||
// and this easily forgotten to be specified
|
||||
Scalar fugCoeff =
|
||||
FluidSystem::template fugacityCoefficient<FluidState, Scalar>(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx,
|
||||
compIdx);
|
||||
(*this)[fugacity0Idx + compIdx] =
|
||||
fugCoeff*fluidState.moleFraction(refPhaseIdx, compIdx)*pRef;
|
||||
}
|
||||
|
||||
// assign pressure of first phase
|
||||
(*this)[pressure0Idx] = FsToolbox::value(fluidState.pressure(/*phaseIdx=*/0));
|
||||
|
||||
// assign first M - 1 saturations
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
|
||||
(*this)[saturation0Idx + phaseIdx] = FsToolbox::value(fluidState.saturation(phaseIdx));
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
59
opm/models/ncp/ncpproperties.hh
Normal file
59
opm/models/ncp/ncpproperties.hh
Normal file
@ -0,0 +1,59 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief Declares the properties required for the NCP compositional
|
||||
* multi-phase model.
|
||||
*/
|
||||
#ifndef EWOMS_NCP_PROPERTIES_HH
|
||||
#define EWOMS_NCP_PROPERTIES_HH
|
||||
|
||||
#include <opm/models/common/multiphasebaseproperties.hh>
|
||||
#include <opm/models/io/vtkcompositionmodule.hh>
|
||||
#include <opm/models/io/vtkenergymodule.hh>
|
||||
#include <opm/models/io/vtkdiffusionmodule.hh>
|
||||
|
||||
BEGIN_PROPERTIES
|
||||
|
||||
//! Enable the energy equation?
|
||||
NEW_PROP_TAG(EnableEnergy);
|
||||
|
||||
//! Enable diffusive fluxes?
|
||||
NEW_PROP_TAG(EnableDiffusion);
|
||||
|
||||
//! The unmodified weight for the pressure primary variable
|
||||
NEW_PROP_TAG(NcpPressureBaseWeight);
|
||||
//! The weight for the saturation primary variables
|
||||
NEW_PROP_TAG(NcpSaturationsBaseWeight);
|
||||
//! The unmodified weight for the fugacity primary variables
|
||||
NEW_PROP_TAG(NcpFugacitiesBaseWeight);
|
||||
|
||||
//! The themodynamic constraint solver which calculates the
|
||||
//! composition of any phase given all component fugacities.
|
||||
NEW_PROP_TAG(NcpCompositionFromFugacitiesSolver);
|
||||
|
||||
END_PROPERTIES
|
||||
|
||||
#endif
|
||||
151
opm/models/ncp/ncpratevector.hh
Normal file
151
opm/models/ncp/ncpratevector.hh
Normal file
@ -0,0 +1,151 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM 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.
|
||||
|
||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::NcpRateVector
|
||||
*/
|
||||
#ifndef EWOMS_NCP_RATE_VECTOR_HH
|
||||
#define EWOMS_NCP_RATE_VECTOR_HH
|
||||
|
||||
#include "ncpindices.hh"
|
||||
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
#include <opm/material/constraintsolvers/NcpFlash.hpp>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \ingroup NcpModel
|
||||
*
|
||||
* \brief Implements a vector representing mass, molar or volumetric rates.
|
||||
*
|
||||
* This class is basically a Dune::FieldVector which can be set using
|
||||
* either mass, molar or volumetric rates.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class NcpRateVector
|
||||
: public Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Evaluation),
|
||||
GET_PROP_VALUE(TypeTag, NumEq)>
|
||||
{
|
||||
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
typedef Dune::FieldVector<Evaluation, numEq> ParentType;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
|
||||
enum { conti0EqIdx = Indices::conti0EqIdx };
|
||||
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
|
||||
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
|
||||
|
||||
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
|
||||
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
public:
|
||||
NcpRateVector() : ParentType()
|
||||
{ Opm::Valgrind::SetUndefined(*this); }
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleRateVector::ImmiscibleRateVector(Scalar)
|
||||
*/
|
||||
NcpRateVector(const Evaluation& value)
|
||||
: ParentType(value)
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const
|
||||
* ImmiscibleRateVector& )
|
||||
*/
|
||||
NcpRateVector(const NcpRateVector& value)
|
||||
: ParentType(value)
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleRateVector::setMassRate
|
||||
*/
|
||||
void setMassRate(const ParentType& value)
|
||||
{
|
||||
// convert to molar rates
|
||||
ParentType molarRate(value);
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
molarRate[conti0EqIdx + compIdx] /= FluidSystem::molarMass(compIdx);
|
||||
|
||||
// set the molar rate
|
||||
setMolarRate(molarRate);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleRateVector::setMolarRate
|
||||
*/
|
||||
void setMolarRate(const ParentType& value)
|
||||
{ ParentType::operator=(value); }
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleRateVector::setEnthalpyRate
|
||||
*/
|
||||
template <class RhsEval>
|
||||
void setEnthalpyRate(const RhsEval& rate)
|
||||
{ EnergyModule::setEnthalpyRate(*this, rate); }
|
||||
|
||||
/*!
|
||||
* \copydoc ImmiscibleRateVector::setVolumetricRate
|
||||
*/
|
||||
template <class FluidState, class RhsEval>
|
||||
void setVolumetricRate(const FluidState& fluidState, unsigned phaseIdx, const RhsEval& volume)
|
||||
{
|
||||
*this = 0.0;
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
(*this)[conti0EqIdx + compIdx] = fluidState.molarity(phaseIdx, compIdx) * volume;
|
||||
|
||||
EnergyModule::setEnthalpyRate(*this, fluidState, phaseIdx, volume);
|
||||
Opm::Valgrind::CheckDefined(*this);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Assignment operator from a scalar or a function evaluation
|
||||
*/
|
||||
template <class RhsEval>
|
||||
NcpRateVector& operator=(const RhsEval& value)
|
||||
{
|
||||
for (unsigned i=0; i < this->size(); ++i)
|
||||
(*this)[i] = value;
|
||||
return *this;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Assignment operator from another rate vector
|
||||
*/
|
||||
NcpRateVector& operator=(const NcpRateVector& other)
|
||||
{
|
||||
for (unsigned i=0; i < this->size(); ++i)
|
||||
(*this)[i] = other[i];
|
||||
return *this;
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue
Block a user