diff --git a/examples/co2injection_flash_ecfv.cpp b/examples/co2injection_flash_ecfv.cpp index 8574db45c..3216811f8 100644 --- a/examples/co2injection_flash_ecfv.cpp +++ b/examples/co2injection_flash_ecfv.cpp @@ -33,7 +33,7 @@ #endif #include -#include +#include #include #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_flash_ni_ecfv.cpp b/examples/co2injection_flash_ni_ecfv.cpp index 09c367f3d..b4964d290 100644 --- a/examples/co2injection_flash_ni_ecfv.cpp +++ b/examples/co2injection_flash_ni_ecfv.cpp @@ -30,7 +30,7 @@ #include #include -#include +#include #include #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_flash_ni_vcfv.cpp b/examples/co2injection_flash_ni_vcfv.cpp index 86c38982a..096300549 100644 --- a/examples/co2injection_flash_ni_vcfv.cpp +++ b/examples/co2injection_flash_ni_vcfv.cpp @@ -30,7 +30,7 @@ #include #include -#include +#include #include #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_flash_vcfv.cpp b/examples/co2injection_flash_vcfv.cpp index 590d4da8c..0f6061abc 100644 --- a/examples/co2injection_flash_vcfv.cpp +++ b/examples/co2injection_flash_vcfv.cpp @@ -33,7 +33,7 @@ #endif #include -#include +#include #include #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh" diff --git a/examples/diffusion_flash.cpp b/examples/diffusion_flash.cpp index aaae2da68..afcb9d658 100644 --- a/examples/diffusion_flash.cpp +++ b/examples/diffusion_flash.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include "problems/diffusionproblem.hh" BEGIN_PROPERTIES diff --git a/opm/models/flash/flashboundaryratevector.hh b/opm/models/flash/flashboundaryratevector.hh new file mode 100644 index 000000000..d15781528 --- /dev/null +++ b/opm/models/flash/flashboundaryratevector.hh @@ -0,0 +1,205 @@ +// -*- 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 . + + 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::FlashBoundaryRateVector + */ +#ifndef EWOMS_FLASH_BOUNDARY_RATE_VECTOR_HH +#define EWOMS_FLASH_BOUNDARY_RATE_VECTOR_HH + +#include "flashproperties.hh" + +#include +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief Implements a boundary vector for the fully implicit + * compositional multi-phase model which is based on flash + * calculations. + */ +template +class FlashBoundaryRateVector : 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 EnergyModule; + typedef Opm::MathToolbox Toolbox; + +public: + FlashBoundaryRateVector() : ParentType() + {} + + /*! + * \copydoc + * ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(Scalar) + */ + FlashBoundaryRateVector(const Evaluation& value) : ParentType(value) + {} + + /*! + * \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(const + * ImmiscibleBoundaryRateVector& ) + */ + FlashBoundaryRateVector(const FlashBoundaryRateVector& value) = default; + FlashBoundaryRateVector& operator=(const FlashBoundaryRateVector& value) = default; + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setFreeFlow + */ + template + 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 + 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) { + Evaluation& val = this->operator[](eqIdx); + val = Toolbox::min(0.0, val); + } + } + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setOutFlow + */ + template + 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) { + Evaluation& val = this->operator[](eqIdx); + val = Toolbox::max(0.0, val); + } + } + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setNoFlow + */ + void setNoFlow() + { (*this) = Evaluation(0.0); } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/flash/flashextensivequantities.hh b/opm/models/flash/flashextensivequantities.hh new file mode 100644 index 000000000..30fa5fd42 --- /dev/null +++ b/opm/models/flash/flashextensivequantities.hh @@ -0,0 +1,96 @@ +// -*- 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 . + + 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::FlashExtensiveQuantities + */ +#ifndef EWOMS_FLASH_EXTENSIVE_QUANTITIES_HH +#define EWOMS_FLASH_EXTENSIVE_QUANTITIES_HH + +#include "flashproperties.hh" + +#include +#include +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * \ingroup ExtensiveQuantities + * + * \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 compositional multi-phase model based on + * flash calculations. + * + * This means pressure and concentration gradients, phase densities at + * the integration point, etc. + */ +template +class FlashExtensiveQuantities + : public MultiPhaseBaseExtensiveQuantities + , public EnergyExtensiveQuantities + , public DiffusionExtensiveQuantities +{ + typedef MultiPhaseBaseExtensiveQuantities ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) }; + typedef Opm::DiffusionExtensiveQuantities DiffusionExtensiveQuantities; + + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + typedef Opm::EnergyExtensiveQuantities 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 + 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 diff --git a/opm/models/flash/flashindices.hh b/opm/models/flash/flashindices.hh new file mode 100644 index 000000000..048e7593d --- /dev/null +++ b/opm/models/flash/flashindices.hh @@ -0,0 +1,72 @@ +// -*- 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 . + + 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::FlashIndices + */ +#ifndef EWOMS_FLASH_INDICES_HH +#define EWOMS_FLASH_INDICES_HH + +#include "flashproperties.hh" +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief Defines the primary variable and equation indices for the + * compositional multi-phase model based on flash calculations. + * + * \tparam PVOffset The first index in a primary variable vector. + */ +template +class FlashIndices + : public EnergyIndices +{ + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + typedef Opm::EnergyIndices EnergyIndices; + +public: + //! number of equations/primary variables + static const int numEq = numComponents + EnergyIndices::numEq_; + + // Primary variable indices + + //! Index of the total concentration of the first component in the + //! pore space. + static const int cTot0Idx = PVOffset; + + // equation indices + + //! Index of the mass conservation equation for the first + //! component. + static const int conti0EqIdx = PVOffset; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/flash/flashintensivequantities.hh b/opm/models/flash/flashintensivequantities.hh new file mode 100644 index 000000000..ffc110ebc --- /dev/null +++ b/opm/models/flash/flashintensivequantities.hh @@ -0,0 +1,217 @@ +// -*- 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 . + + 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::FlashIntensiveQuantities + */ +#ifndef EWOMS_FLASH_INTENSIVE_QUANTITIES_HH +#define EWOMS_FLASH_INTENSIVE_QUANTITIES_HH + +#include "flashproperties.hh" +#include "flashindices.hh" + +#include +#include + +#include +#include + +#include +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * \ingroup IntensiveQuantities + * + * \brief Contains the intensive quantities of the flash-based compositional multi-phase model + */ +template +class FlashIntensiveQuantities + : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) + , public DiffusionIntensiveQuantities + , public EnergyIntensiveQuantities + , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + + // primary variable indices + enum { cTot0Idx = Indices::cTot0Idx }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { dimWorld = GridView::dimensionworld }; + + 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, FlashSolver) FlashSolver; + + typedef Dune::FieldVector ComponentVector; + typedef Dune::FieldMatrix DimMatrix; + + typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities; + typedef Opm::DiffusionIntensiveQuantities DiffusionIntensiveQuantities; + typedef Opm::EnergyIntensiveQuantities EnergyIntensiveQuantities; + +public: + //! The type of the object returned by the fluidState() method + typedef Opm::CompositionalFluidState FluidState; + + FlashIntensiveQuantities() + { } + + FlashIntensiveQuantities(const FlashIntensiveQuantities& other) = default; + + FlashIntensiveQuantities& operator=(const FlashIntensiveQuantities& other) = default; + + /*! + * \copydoc IntensiveQuantities::update + */ + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + { + ParentType::update(elemCtx, dofIdx, timeIdx); + EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx); + + const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); + const auto& problem = elemCtx.problem(); + Scalar flashTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, FlashTolerance); + + // extract the total molar densities of the components + ComponentVector cTotal; + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) + cTotal[compIdx] = priVars.makeEvaluation(cTot0Idx + compIdx, timeIdx); + + const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx); + if (hint) { + // use the same fluid state as the one of the hint, but + // make sure that we don't overwrite the temperature + // specified by the primary variables + Evaluation T = fluidState_.temperature(/*phaseIdx=*/0); + fluidState_.assign(hint->fluidState()); + fluidState_.setTemperature(T); + } + else + FlashSolver::guessInitial(fluidState_, cTotal); + + // compute the phase compositions, densities and pressures + typename FluidSystem::template ParameterCache paramCache; + const MaterialLawParams& materialParams = + problem.materialLawParams(elemCtx, dofIdx, timeIdx); + FlashSolver::template solve(fluidState_, + materialParams, + paramCache, + cTotal, + flashTolerance); + + // calculate relative permeabilities + MaterialLaw::relativePermeabilities(relativePermeability_, + materialParams, fluidState_); + Opm::Valgrind::CheckDefined(relativePermeability_); + + // set the phase viscosities + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + paramCache.updatePhase(fluidState_, phaseIdx); + + const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx); + fluidState_.setViscosity(phaseIdx, mu); + + mobility_[phaseIdx] = relativePermeability_[phaseIdx] / mu; + Opm::Valgrind::CheckDefined(mobility_[phaseIdx]); + } + + ///////////// + // calculate the remaining quantities + ///////////// + + // porosity + porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx); + Opm::Valgrind::CheckDefined(porosity_); + + // 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); + } + + /*! + * \copydoc ImmiscibleIntensiveQuantities::fluidState + */ + const FluidState& fluidState() const + { return fluidState_; } + + /*! + * \copydoc ImmiscibleIntensiveQuantities::intrinsicPermeability + */ + const DimMatrix& intrinsicPermeability() const + { return intrinsicPerm_; } + + /*! + * \copydoc ImmiscibleIntensiveQuantities::relativePermeability + */ + const Evaluation& relativePermeability(unsigned phaseIdx) const + { return relativePermeability_[phaseIdx]; } + + /*! + * \copydoc ImmiscibleIntensiveQuantities::mobility + */ + const Evaluation& mobility(unsigned phaseIdx) const + { + return mobility_[phaseIdx]; + } + + /*! + * \copydoc ImmiscibleIntensiveQuantities::porosity + */ + const Evaluation& porosity() const + { return porosity_; } + +private: + DimMatrix intrinsicPerm_; + FluidState fluidState_; + Evaluation porosity_; + Evaluation relativePermeability_[numPhases]; + Evaluation mobility_[numPhases]; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/flash/flashlocalresidual.hh b/opm/models/flash/flashlocalresidual.hh new file mode 100644 index 000000000..5425ca572 --- /dev/null +++ b/opm/models/flash/flashlocalresidual.hh @@ -0,0 +1,198 @@ +// -*- 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 . + + 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::FlashLocalResidual + */ +#ifndef EWOMS_FLASH_LOCAL_RESIDUAL_HH +#define EWOMS_FLASH_LOCAL_RESIDUAL_HH + +#include "flashproperties.hh" + +#include +#include + +#include + +namespace Opm { +/*! + * \ingroup FlashModel + * + * \brief Calculates the local residual of the compositional multi-phase + * model based on flash calculations. + */ +template +class FlashLocalResidual: public GET_PROP_TYPE(TypeTag, DiscLocalResidual) +{ + 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, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + 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 { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) }; + typedef Opm::DiffusionModule DiffusionModule; + + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + typedef Opm::EnergyModule EnergyModule; + + typedef Opm::MathToolbox Toolbox; + +public: + /*! + * \copydoc ImmiscibleLocalResidual::addPhaseStorage + */ + template + void addPhaseStorage(Dune::FieldVector& storage, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx, + unsigned phaseIdx) const + { + const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + const auto& fs = 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(fs.molarity(phaseIdx, compIdx)) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + } + + EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx); + } + + /*! + * \copydoc FvBaseLocalResidual::computeStorage + */ + template + void computeStorage(Dune::FieldVector& 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 FvBaseLocalResidual::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 finite volume of the current phase + unsigned upIdx = static_cast(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 ImmiscibleLocalResidual::computeSource + */ + 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); + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/flash/flashmodel.hh b/opm/models/flash/flashmodel.hh new file mode 100644 index 000000000..14fa5bda1 --- /dev/null +++ b/opm/models/flash/flashmodel.hh @@ -0,0 +1,318 @@ +// -*- 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 . + + 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::FlashModel + */ +#ifndef EWOMS_FLASH_MODEL_HH +#define EWOMS_FLASH_MODEL_HH + +#include + +#include "flashproperties.hh" +#include "flashprimaryvariables.hh" +#include "flashlocalresidual.hh" +#include "flashratevector.hh" +#include "flashboundaryratevector.hh" +#include "flashintensivequantities.hh" +#include "flashextensivequantities.hh" +#include "flashindices.hh" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Opm { +template +class FlashModel; +} + +BEGIN_PROPERTIES + +//! The type tag for the isothermal single phase problems +NEW_TYPE_TAG(FlashModel, INHERITS_FROM(MultiPhaseBaseModel, + VtkComposition, + VtkEnergy, + VtkDiffusion)); + +//! Use the FlashLocalResidual function for the flash model +SET_TYPE_PROP(FlashModel, LocalResidual, + Opm::FlashLocalResidual); + +//! Use the NCP flash solver by default +SET_TYPE_PROP(FlashModel, FlashSolver, + Opm::NcpFlash); + +//! Let the flash solver choose its tolerance by default +SET_SCALAR_PROP(FlashModel, FlashTolerance, -1.0); + +//! the Model property +SET_TYPE_PROP(FlashModel, Model, Opm::FlashModel); + +//! the PrimaryVariables property +SET_TYPE_PROP(FlashModel, PrimaryVariables, Opm::FlashPrimaryVariables); + +//! the RateVector property +SET_TYPE_PROP(FlashModel, RateVector, Opm::FlashRateVector); + +//! the BoundaryRateVector property +SET_TYPE_PROP(FlashModel, BoundaryRateVector, Opm::FlashBoundaryRateVector); + +//! the IntensiveQuantities property +SET_TYPE_PROP(FlashModel, IntensiveQuantities, Opm::FlashIntensiveQuantities); + +//! the ExtensiveQuantities property +SET_TYPE_PROP(FlashModel, ExtensiveQuantities, Opm::FlashExtensiveQuantities); + +//! The indices required by the flash-baseed isothermal compositional model +SET_TYPE_PROP(FlashModel, Indices, Opm::FlashIndices); + +// The updates of intensive quantities tend to be _very_ expensive for this +// model, so let's try to minimize the number of required ones +SET_BOOL_PROP(FlashModel, EnableIntensiveQuantityCache, true); + +// since thermodynamic hints are basically free if the cache for intensive quantities is +// enabled, and this model usually shows quite a performance improvment if they are +// enabled, let's enable them by default. +SET_BOOL_PROP(FlashModel, EnableThermodynamicHints, true); + +// disable molecular diffusion by default +SET_BOOL_PROP(FlashModel, EnableDiffusion, false); + +//! Disable the energy equation by default +SET_BOOL_PROP(FlashModel, EnableEnergy, false); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief A compositional multi-phase model based on flash-calculations + * + * This model assumes a flow of \f$M \geq 1\f$ fluid phases + * \f$\alpha\f$, each of which is assumed to be a mixture \f$N \geq + * M\f$ chemical species (denoted by the upper index \f$\kappa\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); + * \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] + * + * To determine the quanties that occur in the equations above, this + * model uses flash calculations. A flash solver starts with + * the total mass or molar mass per volume for each component and, + * calculates the compositions, saturations and pressures of all + * phases at a given temperature. For this the flash solver has to use + * some model assumptions internally. (Often these are the same + * primary variable switching or NCP assumptions as used by the other + * fully implicit compositional multi-phase models provided by eWoms.) + * + * Using flash calculations for the flow model has some disadvantages: + * - The accuracy of the flash solver needs to be sufficient to + * calculate the parital derivatives using numerical differentiation + * which are required for the Newton scheme. + * - Flash calculations tend to be quite computationally expensive and + * are often numerically unstable. + * + * It is thus adviced to increase the target tolerance of the Newton + * scheme or a to use type for scalar values which exhibits higher + * precision than the standard \c double (e.g. \c quad) if this model + * ought to be used. + * + * The model uses the following primary variables: + * - The total molar concentration of each component: + * \f$c^\kappa = \sum_\alpha S_\alpha x_\alpha^\kappa \rho_{mol, \alpha}\f$ + * - The absolute temperature $T$ in Kelvins if the energy equation enabled. + */ +template +class FlashModel + : public MultiPhaseBaseModel +{ + typedef MultiPhaseBaseModel ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + + + typedef Opm::EnergyModule EnergyModule; + +public: + FlashModel(Simulator& simulator) + : ParentType(simulator) + {} + + /*! + * \brief Register all run-time parameters for the immiscible model. + */ + static void registerParameters() + { + ParentType::registerParameters(); + + // register runtime parameters of the VTK output modules + Opm::VtkCompositionModule::registerParameters(); + + if (enableDiffusion) + Opm::VtkDiffusionModule::registerParameters(); + + if (enableEnergy) + Opm::VtkEnergyModule::registerParameters(); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, FlashTolerance, + "The maximum tolerance for the flash solver to " + "consider the solution converged"); + } + + /*! + * \copydoc FvBaseDiscretization::name + */ + static std::string name() + { return "flash"; } + + /*! + * \copydoc FvBaseDiscretization::primaryVarName + */ + std::string primaryVarName(unsigned pvIdx) const + { + const std::string& tmp = EnergyModule::primaryVarName(pvIdx); + if (tmp != "") + return tmp; + + std::ostringstream oss; + if (Indices::cTot0Idx <= pvIdx && pvIdx < Indices::cTot0Idx + + numComponents) + oss << "c_tot," << FluidSystem::componentName(/*compIdx=*/pvIdx + - Indices::cTot0Idx); + else + assert(false); + + return oss.str(); + } + + /*! + * \copydoc FvBaseDiscretization::eqName + */ + std::string eqName(unsigned eqIdx) const + { + const std::string& tmp = EnergyModule::eqName(eqIdx); + if (tmp != "") + return tmp; + + std::ostringstream oss; + if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx + + numComponents) { + unsigned compIdx = eqIdx - Indices::conti0EqIdx; + oss << "continuity^" << FluidSystem::componentName(compIdx); + } + else + assert(false); + + return oss.str(); + } + + /*! + * \copydoc FvBaseDiscretization::primaryVarWeight + */ + Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const + { + Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx); + if (tmp > 0) + return tmp; + + unsigned compIdx = pvIdx - Indices::cTot0Idx; + + // make all kg equal. also, divide the weight of all total + // compositions by 100 to make the relative errors more + // comparable to the ones of the other models (at 10% porosity + // the medium is fully saturated with water at atmospheric + // conditions if 100 kg/m^3 are present!) + return FluidSystem::molarMass(compIdx) / 100.0; + } + + /*! + * \copydoc FvBaseDiscretization::eqWeight + */ + Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const + { + Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx); + if (tmp > 0) + return tmp; + + unsigned compIdx = eqIdx - Indices::conti0EqIdx; + + // make all kg equal + return FluidSystem::molarMass(compIdx); + } + + void registerOutputModules_() + { + ParentType::registerOutputModules_(); + + // add the VTK output modules which are meaningful for the model + this->addOutputModule(new Opm::VtkCompositionModule(this->simulator_)); + if (enableDiffusion) + this->addOutputModule(new Opm::VtkDiffusionModule(this->simulator_)); + if (enableEnergy) + this->addOutputModule(new Opm::VtkEnergyModule(this->simulator_)); + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/flash/flashprimaryvariables.hh b/opm/models/flash/flashprimaryvariables.hh new file mode 100644 index 000000000..29aa95b1a --- /dev/null +++ b/opm/models/flash/flashprimaryvariables.hh @@ -0,0 +1,150 @@ +// -*- 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 . + + 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::FlashPrimaryVariables + */ +#ifndef EWOMS_FLASH_PRIMARY_VARIABLES_HH +#define EWOMS_FLASH_PRIMARY_VARIABLES_HH + +#include "flashindices.hh" +#include "flashproperties.hh" + +#include +#include + +#include +#include +#include +#include + +#include + +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief Represents the primary variables used by the compositional + * flow model based on flash calculations. + * + * This class is basically a Dune::FieldVector which can retrieve its + * contents from an aribitatry fluid state. + */ +template +class FlashPrimaryVariables : public FvBasePrimaryVariables +{ + typedef FvBasePrimaryVariables ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + // primary variable indices + enum { cTot0Idx = Indices::cTot0Idx }; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + + typedef typename Opm::MathToolbox Toolbox; + typedef Opm::EnergyModule EnergyModule; + +public: + FlashPrimaryVariables() : ParentType() + { Opm::Valgrind::SetDefined(*this); } + + /*! + * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(Scalar) + */ + FlashPrimaryVariables(Scalar value) : ParentType(value) + { + Opm::Valgrind::CheckDefined(value); + Opm::Valgrind::SetDefined(*this); + } + + /*! + * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(const + * ImmisciblePrimaryVariables& ) + */ + FlashPrimaryVariables(const FlashPrimaryVariables& value) = default; + FlashPrimaryVariables& operator=(const FlashPrimaryVariables& value) = default; + + /*! + * \copydoc ImmisciblePrimaryVariables::assignMassConservative + */ + template + void assignMassConservative(const FluidState& fluidState, + const MaterialLawParams& matParams OPM_UNUSED, + bool isInEquilibrium OPM_UNUSED= false) + { + // there is no difference between naive and mass conservative + // assignment in the flash model. (we only need the total + // concentrations.) + assignNaive(fluidState); + } + + /*! + * \copydoc ImmisciblePrimaryVariables::assignNaive + */ + template + void assignNaive(const FluidState& fluidState) + { + // reset everything + (*this) = 0.0; + + // assign the phase temperatures. this is out-sourced to + // the energy module + EnergyModule::setPriVarTemperatures(*this, fluidState); + + // determine the phase presence. + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + this->operator[](cTot0Idx + compIdx) += + fluidState.molarity(phaseIdx, compIdx) * fluidState.saturation(phaseIdx); + } + } + } + + /*! + * \brief Prints the names of the primary variables and their values. + * + * \param os The \c std::ostream which should be used for the output. + */ + void print(std::ostream& os = std::cout) const + { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + os << "(c_tot," << FluidSystem::componentName(compIdx) << " = " + << this->operator[](cTot0Idx + compIdx); + } + os << ")" << std::flush; + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/flash/flashproperties.hh b/opm/models/flash/flashproperties.hh new file mode 100644 index 000000000..a679a3f3c --- /dev/null +++ b/opm/models/flash/flashproperties.hh @@ -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 . + + 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 FlashModel + * + * \brief Declares the properties required by the compositional + * multi-phase model based on flash calculations. + */ +#ifndef EWOMS_FLASH_PROPERTIES_HH +#define EWOMS_FLASH_PROPERTIES_HH + +#include +#include +#include +#include + +BEGIN_PROPERTIES + +//! Provides the thermodynamic relations +NEW_PROP_TAG(FluidSystem); +//! The type of the flash constraint solver +NEW_PROP_TAG(FlashSolver); +//! The maximum accepted error of the flash solver +NEW_PROP_TAG(FlashTolerance); + +//! The thermal conduction law which ought to be used +NEW_PROP_TAG(ThermalConductionLaw); +//! The parameters of the thermal conduction law +NEW_PROP_TAG(ThermalConductionLawParams); + +//! Specifies whether energy should be considered as a conservation quantity or not +NEW_PROP_TAG(EnableEnergy); +//! Enable diffusive fluxes? +NEW_PROP_TAG(EnableDiffusion); + +END_PROPERTIES + +#endif diff --git a/opm/models/flash/flashratevector.hh b/opm/models/flash/flashratevector.hh new file mode 100644 index 000000000..c18d99958 --- /dev/null +++ b/opm/models/flash/flashratevector.hh @@ -0,0 +1,144 @@ +// -*- 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 . + + 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::FlashRateVector + */ +#ifndef EWOMS_FLASH_RATE_VECTOR_HH +#define EWOMS_FLASH_RATE_VECTOR_HH + +#include + +#include +#include +#include + +#include "flashintensivequantities.hh" + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \copydoc Opm::ImmiscibleRateVector + */ +template +class FlashRateVector + : public Dune::FieldVector +{ + 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, Indices) Indices; + + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + + typedef Dune::FieldVector ParentType; + typedef Opm::EnergyModule EnergyModule; + +public: + FlashRateVector() : ParentType() + { Opm::Valgrind::SetUndefined(*this); } + + /*! + * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(Scalar) + */ + FlashRateVector(const Evaluation& value) : ParentType(value) + {} + + /*! + * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const + * ImmiscibleRateVector& ) + */ + FlashRateVector(const FlashRateVector& 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); + + setMolarRate(molarRate); + } + + /*! + * \copydoc ImmiscibleRateVector::setMolarRate + */ + void setMolarRate(const ParentType& value) + { ParentType::operator=(value); } + + /*! + * \copydoc ImmiscibleRateVector::setEnthalpyRate + */ + void setEnthalpyRate(const Evaluation& rate) + { EnergyModule::setEnthalpyRate(*this, rate); } + + /*! + * \copydoc ImmiscibleRateVector::setVolumetricRate + */ + template + void setVolumetricRate(const FluidState& fluidState, unsigned phaseIdx, const RhsEval& volume) + { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) + (*this)[conti0EqIdx + compIdx] = + fluidState.density(phaseIdx, compIdx) + * fluidState.moleFraction(phaseIdx, compIdx) + * volume; + + EnergyModule::setEnthalpyRate(*this, fluidState, phaseIdx, volume); + } + + /*! + * \brief Assignment operator from a scalar or a function evaluation + */ + template + FlashRateVector& operator=(const RhsEval& value) + { + for (unsigned i=0; i < this->size(); ++i) + (*this)[i] = value; + return *this; + } + + /*! + * \brief Assignment operator from another rate vector + */ + FlashRateVector& operator=(const FlashRateVector& other) + { + for (unsigned i=0; i < this->size(); ++i) + (*this)[i] = other[i]; + return *this; + } +}; + +} // namespace Opm + +#endif