OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

a simple 1D test demonstrating the usage of the PTFlash compostional modeling

Browse Source
This commit is contained in:
Trine Mykkeltvedt 2022-05-02 14:01:37 +02:00 committed by Kai Bao
parent 6c28415a7e
commit 83cb5670b7
23 changed files with 2661 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
examples/problems/co2injectiontestflash.hh Normal file
View File

@ -0,0 +1,66 @@
// -*- 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::Co2InjectionTestFlash
*/
#ifndef EWOMS_CO2_INJECTION_TEST_FLASH_HH
#define EWOMS_CO2_INJECTION_TEST_FLASH_HH
#include <opm/material/constraintsolvers/NcpFlash.hpp>
namespace Opm {
/*!
* \brief Flash solver used by the CO2 injection test problem.
*
* This class is just the NCP flash solver with the guessInitial()
* method that is adapted to the pressure regime of the CO2 injection
* problem.
*/
template <class Scalar, class FluidSystem>
class Co2InjectionTestFlash : public Opm::NcpFlash<Scalar, FluidSystem>
{
using ParentType = Opm::NcpFlash<Scalar, FluidSystem>;
enum { numPhases = FluidSystem::numPhases };
public:
/*!
* \brief Guess initial values for all quantities.
*/
template <class FluidState, class ComponentVector>
static void guessInitial(FluidState& fluidState, const ComponentVector& globalMolarities)
{
ParentType::guessInitial(fluidState, globalMolarities);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// pressure. use something close to the reservoir pressure as initial guess
fluidState.setPressure(phaseIdx, 100e5);
}
}
};
} // namespace Opm
#endif // EWOMS_CO2_INJECTION_TEST_FLASH_HH

667
examples/problems/simpletestproblem.hh Normal file
View File

@ -0,0 +1,667 @@
// -*- 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::simpletestproblem
*/
#ifndef EWOMS_SIMPLETEST_PROBLEM_HH
#define EWOMS_SIMPLETEST_PROBLEM_HH
#include <opm/common/Exceptions.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/BrooksCorey.hpp>
#include <opm/material/constraintsolvers/PTFlash.hpp>
#include <opm/material/fluidsystems/ThreeComponentFluidSystem.hh>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include <opm/models/ptflash/flashmodel.hh>
#include <opm/models/io/structuredgridvanguard.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/start.hh>
#include <opm/simulators/linalg/parallelistlbackend.hh>
#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <sstream>
#include <iostream>
#include <string>
namespace Opm {
template <class TypeTag>
class SimpleTest;
} // namespace Opm
namespace Opm::Properties {
namespace TTag {
struct SimpleTest {};
} // end namespace TTag
// declare the "simpletest" problem specify property tags
template <class TypeTag, class MyTypeTag>
struct Temperature { using type = UndefinedProperty; };
template <class TypeTag, class MyTypeTag>
struct SimulationName { using type = UndefinedProperty; };
template <class TypeTag, class MyTypeTag>
struct EpisodeLength { using type = UndefinedProperty;};
template <class TypeTag, class MyTypeTag>
struct Initialpressure { using type = UndefinedProperty;};
// Set the grid type: --->1D
template <class TypeTag>
struct Grid<TypeTag, TTag::SimpleTest> { using type = Dune::YaspGrid</*dim=*/2>; };
// Set the problem property
template <class TypeTag>
struct Problem<TypeTag, TTag::SimpleTest>
{ using type = Opm::SimpleTest<TypeTag>; };
// Set flash solver
template <class TypeTag>
struct FlashSolver<TypeTag, TTag::SimpleTest> {
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
public:
using type = Opm::PTFlash<Scalar, FluidSystem>;
};
// Set fluid configuration
template <class TypeTag>
struct FluidSystem<TypeTag, TTag::SimpleTest>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = Opm::ThreeComponentFluidSystem<Scalar>;
};
// Set the material Law
template <class TypeTag>
struct MaterialLaw<TypeTag, TTag::SimpleTest>
{
private:
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
// /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, TODO
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
// define the material law which is parameterized by effective saturations
using EffMaterialLaw = Opm::NullMaterial<Traits>;
//using EffMaterialLaw = Opm::BrooksCorey<Traits>;
public:
using type = EffMaterialLaw;
};
// Write the Newton convergence behavior to disk?
template <class TypeTag>
struct NewtonWriteConvergence<TypeTag, TTag::SimpleTest> {
static constexpr bool value = false; };
// Enable gravity false
template <class TypeTag>
struct EnableGravity<TypeTag, TTag::SimpleTest> { static constexpr bool value = false;
};
// set the defaults for the problem specific properties
template <class TypeTag>
struct Temperature<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 423.25;//TODO
};
template <class TypeTag>
struct Initialpressure<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e5;
};
template <class TypeTag>
struct SimulationName<TypeTag, TTag::SimpleTest> {
static constexpr auto value = "simpletest";
};
// The default for the end time of the simulation
template <class TypeTag>
struct EndTime<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 3 * 24. * 60. * 60.;//3 days
};
// convergence control
template <class TypeTag>
struct InitialTimeStepSize<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
//static constexpr type value = 30;
static constexpr type value = 1 * 24. * 60. * 60.;
};
template <class TypeTag>
struct LinearSolverTolerance<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e-3;
};
template <class TypeTag>
struct LinearSolverAbsTolerance<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 0.;
};
template <class TypeTag>
struct NewtonTolerance<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e-3;
};
template <class TypeTag>
struct MaxTimeStepSize<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 60 * 60;
};
template <class TypeTag>
struct NewtonMaxIterations<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 30;
};
template <class TypeTag>
struct NewtonTargetIterations<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 6;
};
// output
template <class TypeTag>
struct VtkWriteFilterVelocities<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWritePotentialGradients<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteTotalMassFractions<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteTotalMoleFractions<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteFugacityCoeffs<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteLiquidMoleFractions<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteEquilibriumConstants<TypeTag, TTag::SimpleTest> {
static constexpr bool value = true;
};
// write restart for every hour
template <class TypeTag>
struct EpisodeLength<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 60. * 60.;
};
// mesh grid
template <class TypeTag>
struct Vanguard<TypeTag, TTag::SimpleTest> {
using type = Opm::StructuredGridVanguard<TypeTag>;
};
//\Note: from the Julia code, the problem is a 1D problem with 3X1 cell.
//\Note: DomainSizeX is 3.0 meters
//\Note: DomainSizeY is 1.0 meters
template <class TypeTag>
struct DomainSizeX<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 3; // meter
};
template <class TypeTag>
struct DomainSizeY<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0;
};
template <class TypeTag>
struct DomainSizeZ<TypeTag, TTag::SimpleTest> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0;
};
template<class TypeTag>
struct CellsX<TypeTag, TTag::SimpleTest> { static constexpr int value = 3; };
template<class TypeTag>
struct CellsY<TypeTag, TTag::SimpleTest> { static constexpr int value = 1; };
template<class TypeTag>
struct CellsZ<TypeTag, TTag::SimpleTest> { static constexpr int value = 1; };
// compositional, with diffusion
template <class TypeTag>
struct EnableEnergy<TypeTag, TTag::SimpleTest> {
static constexpr bool value = false;
};
} // namespace Opm::Properties
namespace Opm {
/*!
* \ingroup TestProblems
*
* \brief 3 component simple testproblem with ["CO2", "C1", "C10"]
*
*/
template <class TypeTag>
class SimpleTest : public GetPropType<TypeTag, Properties::BaseProblem>
{
using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
enum { dim = GridView::dimension };
enum { dimWorld = GridView::dimensionworld };
using Indices = GetPropType<TypeTag, Properties::Indices>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Model = GetPropType<TypeTag, Properties::Model>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
using Toolbox = Opm::MathToolbox<Evaluation>;
using CoordScalar = typename GridView::ctype;
enum { numPhases = FluidSystem::numPhases };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { Comp2Idx = FluidSystem::Comp2Idx };
enum { Comp1Idx = FluidSystem::Comp1Idx };
enum { Comp0Idx = FluidSystem::Comp0Idx };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { contiCO2EqIdx = conti0EqIdx + Comp1Idx };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
enum { enableGravity = getPropValue<TypeTag, Properties::EnableGravity>() };
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using DimVector = Dune::FieldVector<Scalar, dimWorld>;
using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
using FlashSolver = GetPropType<TypeTag, Properties::FlashSolver>;
public:
using FluidState = Opm::CompositionalFluidState<Evaluation, FluidSystem, enableEnergy>;
/*!
* \copydoc Doxygen::defaultProblemConstructor
*/
SimpleTest(Simulator& simulator)
: ParentType(simulator)
{
Scalar epi_len = EWOMS_GET_PARAM(TypeTag, Scalar, EpisodeLength);
simulator.setEpisodeLength(epi_len);
}
void initPetrophysics()
{
temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
K_ = this->toDimMatrix_(9.869232667160131e-14);
porosity_ = 0.1;
}
template <class Context>
const DimVector&
gravity([[maybe_unused]]const Context& context, [[maybe_unused]] unsigned spaceIdx, [[maybe_unused]] unsigned timeIdx) const
{
return gravity();
}
const DimVector& gravity() const
{
return gravity_;
}
/*!
* \copydoc FvBaseProblem::finishInit
*/
void finishInit()
{
ParentType::finishInit();
// initialize fixed parameters; temperature, permeability, porosity
initPetrophysics();
}
/*!
* \copydoc FvBaseMultiPhaseProblem::registerParameters
*/
static void registerParameters()
{
ParentType::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
"The temperature [K] in the reservoir");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, Initialpressure,
"The initial pressure [Pa s] in the reservoir");
EWOMS_REGISTER_PARAM(TypeTag,
std::string,
SimulationName,
"The name of the simulation used for the output "
"files");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, EpisodeLength,
"Time interval [s] for episode length");
}
/*!
* \copydoc FvBaseProblem::name
*/
std::string name() const
{
std::ostringstream oss;
oss << EWOMS_GET_PARAM(TypeTag, std::string, SimulationName);
return oss.str();
}
// This method must be overridden for the simulator to continue with
// a new episode. We just start a new episode with the same length as
// the old one.
void endEpisode()
{
Scalar epi_len = EWOMS_GET_PARAM(TypeTag, Scalar, EpisodeLength);
this->simulator().startNextEpisode(epi_len);
}
// only write output when episodes change, aka. report steps, and
// include the initial timestep too
bool shouldWriteOutput()
{
return this->simulator().episodeWillBeOver() || (this->simulator().timeStepIndex() == -1);
}
// we don't care about doing restarts from every fifth timestep, it
// will just slow us down
bool shouldWriteRestartFile()
{
return false;
}
/*!
* \copydoc FvBaseProblem::endTimeStep
*/
void endTimeStep()
{
Scalar tol = this->model().newtonMethod().tolerance() * 1e5;
this->model().checkConservativeness(tol);
// Calculate storage terms
PrimaryVariables storageO, storageW;
this->model().globalPhaseStorage(storageO, oilPhaseIdx);
// Calculate storage terms
PrimaryVariables storageL, storageG;
this->model().globalPhaseStorage(storageL, /*phaseIdx=*/0);
this->model().globalPhaseStorage(storageG, /*phaseIdx=*/1);
// Write mass balance information for rank 0
// if (this->gridView().comm().rank() == 0) {
// std::cout << "Storage: liquid=[" << storageL << "]"
// << " gas=[" << storageG << "]\n" << std::flush;
// }
}
/*!
* \copydoc FvBaseProblem::initial
*/
template <class Context>
void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
Opm::CompositionalFluidState<Evaluation, FluidSystem> fs;
initialFluidState(fs, context, spaceIdx, timeIdx);
values.assignNaive(fs);
std::cout << "primary variables for cell " << context.globalSpaceIndex(spaceIdx, timeIdx) << ": " << values << "\n";
}
// Constant temperature
template <class Context>
Scalar temperature([[maybe_unused]] const Context& context, [[maybe_unused]] unsigned spaceIdx, [[maybe_unused]] unsigned timeIdx) const
{
return temperature_;
}
// Constant permeability
template <class Context>
const DimMatrix& intrinsicPermeability([[maybe_unused]] const Context& context,
[[maybe_unused]] unsigned spaceIdx,
[[maybe_unused]] unsigned timeIdx) const
{
return K_;
}
// Constant porosity
template <class Context>
Scalar porosity([[maybe_unused]] const Context& context, [[maybe_unused]] unsigned spaceIdx, [[maybe_unused]] unsigned timeIdx) const
{
int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
int inj = 0;
int prod = 2;
if (spatialIdx == inj || spatialIdx == prod)
return 1.0;
else
return porosity_;
}
/*!
* \copydoc FvBaseMultiPhaseProblem::materialLawParams
*/
template <class Context>
const MaterialLawParams& materialLawParams([[maybe_unused]] const Context& context,
[[maybe_unused]] unsigned spaceIdx,
[[maybe_unused]] unsigned timeIdx) const
{
return this->mat_;
}
// No flow (introduce fake wells instead)
template <class Context>
void boundary(BoundaryRateVector& values,
const Context& /*context*/,
unsigned /*spaceIdx*/,
unsigned /*timeIdx*/) const
{ values.setNoFlow(); }
// No source terms
template <class Context>
void source(RateVector& source_rate,
[[maybe_unused]] const Context& context,
[[maybe_unused]] unsigned spaceIdx,
[[maybe_unused]] unsigned timeIdx) const
{
source_rate = Scalar(0.0);
}
private:
/*!
* \copydoc FvBaseProblem::initial
*/
template <class FluidState, class Context>
void initialFluidState(FluidState& fs, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
using Scalar = double;
using FluidSystem = Opm::ThreeComponentFluidSystem<Scalar>;
constexpr auto numComponents = FluidSystem::numComponents;
using Evaluation = Opm::DenseAd::Evaluation<double, numComponents>;
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
// input from Olav
//z0 = [0.5, 0.3, 0.2]
//zi = [0.99, 0.01-1e-3, 1e-3]
//p0 = 75e5
//T0 = 423.25
int inj = 0;
int prod = 2;
int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
ComponentVector comp;
comp[0] = Evaluation::createVariable(0.5, 1);
comp[1] = Evaluation::createVariable(0.3, 2);
comp[2] = 1. - comp[0] - comp[1];
if (spatialIdx == inj){
comp[0] = Evaluation::createVariable(0.99, 1);
comp[1] = Evaluation::createVariable(0.01-1e-3, 2);
comp[2] = 1. - comp[0] - comp[1];
}
ComponentVector sat;
sat[0] = 1.0; sat[1] = 1.0-sat[0];
// TODO: should we put the derivative against the temperature here?
const Scalar temp = 423.25;
// TODO: no capillary pressure for now
Scalar p0 = 75e5; //CONVERGENCE FAILURE WITH 75
//\Note, for an AD variable, if we multiply it with 2, the derivative will also be scalced with 2,
//\Note, so we should not do it.
if (spatialIdx == inj){
p0 *= 2.0;
}
if (spatialIdx == prod) {
p0 *= 0.5;
}
Evaluation p_init = Evaluation::createVariable(p0, 0);
fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
// It is used here only for calculate the z
fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
fs.setTemperature(temp);
// ParameterCache paramCache;
{
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
}
// ComponentVector zInit(0.); // TODO; zInit needs to be normalized.
// {
// Scalar sumMoles = 0.0;
// for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
// for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
// Scalar tmp = Opm::getValue(fs.molarity(phaseIdx, compIdx) * fs.saturation(phaseIdx));
// zInit[compIdx] += Opm::max(tmp, 1e-8);
// sumMoles += tmp;
// }
// }
// zInit /= sumMoles;
// // initialize the derivatives
// // TODO: the derivative eventually should be from the reservoir flow equations
// Evaluation z_last = 1.;
// for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
// zInit[compIdx] = Evaluation::createVariable(Opm::getValue(zInit[compIdx]), compIdx + 1);
// z_last -= zInit[compIdx];
// }
// zInit[numComponents - 1] = z_last;
// }
// TODO: only, p, z need the derivatives.
const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
//const int flash_verbosity = 1;
const std::string flash_twophase_method = "newton"; // "ssi"
//const std::string flash_twophase_method = "ssi";
// const std::string flash_twophase_method = "ssi+newton";
// Set initial K and L
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation Ktmp = fs.wilsonK_(compIdx);
fs.setKvalue(compIdx, Ktmp);
}
const Evaluation& Ltmp = -1.0;
fs.setLvalue(Ltmp);
}
DimMatrix K_;
Scalar porosity_;
Scalar temperature_;
MaterialLawParams mat_;
DimVector gravity_;
};
} // namespace Opm
#endif

27
opm/models/common/darcyfluxmodule.hh
View File

@ -244,7 +244,9 @@ protected:
Evaluation pStatIn;
if (std::is_same<Scalar, Evaluation>::value ||
interiorDofIdx_ == static_cast<int>(focusDofIdx))
interiorDofIdx_ == static_cast<int>(focusDofIdx) ||
(phaseIdx == 0 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0)) ||
(phaseIdx == 1 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1)))
{
const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
pStatIn = - rhoIn*(gIn*distVecIn);
@ -259,7 +261,9 @@ protected:
Evaluation pStatEx;
if (std::is_same<Scalar, Evaluation>::value ||
exteriorDofIdx_ == static_cast<int>(focusDofIdx))
exteriorDofIdx_ == static_cast<int>(focusDofIdx) ||
(phaseIdx == 0 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1)) ||
(phaseIdx == 1 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0)))
{
const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx);
pStatEx = - rhoEx*(gEx*distVecEx);
@ -274,7 +278,22 @@ protected:
// control volume centers and the length (pStaticExterior -
// pStaticInterior)/distanceInteriorToExterior
Dune::FieldVector<Evaluation, dimWorld> f(distVecTotal);
f *= (pStatEx - pStatIn)/absDistTotalSquared;
if (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0))
f *= -(pStatIn + pStatIn)/absDistTotalSquared;
else if (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1))
f *= (pStatEx + pStatEx)/absDistTotalSquared;
else if (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 0))
f *= 0.0;
else if (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1))
f *= -(pStatIn + pStatIn)/absDistTotalSquared;
else if (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0))
f *= (pStatEx + pStatEx)/absDistTotalSquared;
else if (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 1))
f *= 0.0;
else
f *= (pStatEx - pStatIn)/absDistTotalSquared;
// calculate the final potential gradient
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
@ -559,4 +578,4 @@ protected:
} // namespace Opm
#endif
#endif

12
opm/models/discretization/common/fvbaseadlocallinearizer.hh
View File

@ -273,18 +273,26 @@ protected:
{
const auto& resid = localResidual_.residual();
for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++)
for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++){
residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value();
std::cout << "residual = " << residual_[focusDofIdx][eqIdx] << std::endl;
}
int spatialIdx = elemCtx.globalSpaceIndex(focusDofIdx, /*timeIdx=*/0);
std::cout << "Focus cell = " << spatialIdx << std::endl;
size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
int spatialIdx2 = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
// A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
// the residual function 'eqIdx' for the degree of freedom 'dofIdx'
// with regard to the focus variable 'pvIdx' of the degree of freedom
// 'focusDofIdx'
jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx);
//std::cout << "J[" << spatialIdx2 << "][" << spatialIdx << "][" << eqIdx << "][" << pvIdx << "] = " << jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] << std::endl;
std::cout << "J[" << dofIdx+1 << "," << focusDofIdx+1 << "][" << eqIdx+1 << "," << pvIdx+1 << "] = " << jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] << std::endl;
Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
}
}

8
opm/models/discretization/common/fvbasediscretization.hh
View File

@ -274,12 +274,12 @@ struct EnableConstraints<TypeTag, TTag::FvBaseDiscretization> { static constexpr
// impact because of the intensive quantity cache will cause additional pressure on the
// CPU caches...
template<class TypeTag>
struct EnableIntensiveQuantityCache<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
struct EnableIntensiveQuantityCache<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = true; };
// do not use thermodynamic hints by default. If you enable this, make sure to also
// enable the intensive quantity cache above to avoid getting an exception...
template<class TypeTag>
struct EnableThermodynamicHints<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
struct EnableThermodynamicHints<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = true; };
// if the deflection of the newton method is large, we do not need to solve the linear
// approximation accurately. Assuming that the value for the current solution is quite
@ -309,9 +309,9 @@ struct TimeDiscHistorySize<TypeTag, TTag::FvBaseDiscretization> { static constex
template<class TypeTag>
struct ExtensiveStorageTerm<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
// use volumetric residuals is default
// use volumetric residuals is not default
template<class TypeTag>
struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = true; };
struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
//! eWoms is mainly targeted at research, so experimental features are enabled by
//! default.

4
opm/models/discretization/common/fvbaseelementcontext.hh
View File

@ -554,8 +554,8 @@ protected:
const PrimaryVariables& dofSol = globalSol[globalIdx];
dofVars_[dofIdx].priVars[timeIdx] = &dofSol;
dofVars_[dofIdx].thermodynamicHint[timeIdx] =
model().thermodynamicHint(globalIdx, timeIdx);
dofVars_[dofIdx].thermodynamicHint[timeIdx] = model().thermodynamicHint(globalIdx, timeIdx);
dofVars_[dofIdx].thermodynamicHint[1] = model().thermodynamicHint(globalIdx, 1);
const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
if (cachedIntQuants) {

3
opm/models/discretization/common/fvbasefdlocallinearizer.hh
View File

@ -490,6 +490,8 @@ protected:
unsigned pvIdx)
{
size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
int spatialIdx = elemCtx.globalSpaceIndex(focusDofIdx, /*timeIdx=*/0);
std::cout << "Focus cell = " << spatialIdx << std::endl;
for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
// A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the
@ -497,6 +499,7 @@ protected:
// regard to the primary variable 'pvIdx' of the degree of freedom
// 'focusDofIdx'
jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
std::cout << "J[" << dofIdx << "][" << focusDofIdx << "][" << eqIdx << "][" << pvIdx << "] = " << jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] << std::endl;
Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
}
}

6
opm/models/discretization/common/fvbaselocalresidual.hh
View File

@ -161,6 +161,8 @@ public:
assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
residual = 0.0;
std::cout << "residual from elemtCtx" << std::endl;
// evaluate the flux terms
asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
@ -324,8 +326,10 @@ public:
assert(alpha > 0.0);
assert(isfinite(alpha));
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx){
flux[eqIdx] *= alpha;
std::cout << " flux in eqIdx " << eqIdx << " is: " << flux[eqIdx] << std::endl;
}
// The balance equation for a finite volume is given by
//

46
opm/models/io/vtkcompositionmodule.hh
View File

@ -59,6 +59,10 @@ template<class TypeTag, class MyTypeTag>
struct VtkWriteFugacities { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct VtkWriteFugacityCoeffs { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct VtkWriteLiquidMoleFractions { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct VtkWriteEquilibriumConstants { using type = UndefinedProperty; };
// set default values for what quantities to output
template<class TypeTag>
@ -75,6 +79,10 @@ template<class TypeTag>
struct VtkWriteFugacities<TypeTag, TTag::VtkComposition> { static constexpr bool value = false; };
template<class TypeTag>
struct VtkWriteFugacityCoeffs<TypeTag, TTag::VtkComposition> { static constexpr bool value = false; };
template<class TypeTag>
struct VtkWriteLiquidMoleFractions<TypeTag, TTag::VtkComposition> { static constexpr bool value = false; };
template<class TypeTag>
struct VtkWriteEquilibriumConstants<TypeTag, TTag::VtkComposition> { static constexpr bool value = false; };
} // namespace Opm::Properties
@ -112,6 +120,7 @@ class VtkCompositionModule : public BaseOutputModule<TypeTag>
using ComponentBuffer = typename ParentType::ComponentBuffer;
using PhaseComponentBuffer = typename ParentType::PhaseComponentBuffer;
using ScalarBuffer = typename ParentType::ScalarBuffer;
public:
VtkCompositionModule(const Simulator& simulator)
@ -137,6 +146,10 @@ public:
"Include component fugacities in the VTK output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFugacityCoeffs,
"Include component fugacity coefficients in the VTK output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteLiquidMoleFractions,
"Include liquid mole fractions (L) in the VTK output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteEquilibriumConstants,
"Include equilibrium constants (K) in the VTK output files");
}
/*!
@ -155,7 +168,10 @@ public:
this->resizeComponentBuffer_(totalMoleFrac_);
if (molarityOutput_())
this->resizePhaseComponentBuffer_(molarity_);
if (LOutput_())
this->resizeScalarBuffer_(L_);
if (equilConstOutput_())
this->resizeComponentBuffer_(K_);
if (fugacityOutput_())
this->resizeComponentBuffer_(fugacity_);
if (fugacityCoeffOutput_())
@ -178,6 +194,9 @@ public:
const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
if (LOutput_())
L_[I] = Toolbox::value(fs.L());
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
if (moleFracOutput_())
@ -222,6 +241,10 @@ public:
}
if (fugacityOutput_())
fugacity_[compIdx][I] = Toolbox::value(intQuants.fluidState().fugacity(/*phaseIdx=*/0, compIdx));
if (equilConstOutput_())
K_[compIdx][I] = Toolbox::value(fs.K(compIdx));
}
}
}
@ -236,7 +259,7 @@ public:
return;
}
if (moleFracOutput_())
if (moleFracOutput_())
this->commitPhaseComponentBuffer_(baseWriter, "moleFrac_%s^%s", moleFrac_);
if (massFracOutput_())
this->commitPhaseComponentBuffer_(baseWriter, "massFrac_%s^%s", massFrac_);
@ -246,11 +269,15 @@ public:
this->commitComponentBuffer_(baseWriter, "totalMassFrac^%s", totalMassFrac_);
if (totalMoleFracOutput_())
this->commitComponentBuffer_(baseWriter, "totalMoleFrac^%s", totalMoleFrac_);
if (equilConstOutput_())
this->commitComponentBuffer_(baseWriter, "K^%s", K_);
if (fugacityOutput_())
this->commitComponentBuffer_(baseWriter, "fugacity^%s", fugacity_);
if (fugacityCoeffOutput_())
this->commitPhaseComponentBuffer_(baseWriter, "fugacityCoeff_%s^%s", fugacityCoeff_);
if (LOutput_())
this->commitScalarBuffer_(baseWriter, "L", L_);
}
private:
@ -295,15 +322,30 @@ private:
static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFugacityCoeffs);
return val;
}
static bool LOutput_()
{
static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteLiquidMoleFractions);
return val;
}
static bool equilConstOutput_()
{
static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteEquilibriumConstants);
return val;
}
PhaseComponentBuffer moleFrac_;
PhaseComponentBuffer massFrac_;
PhaseComponentBuffer molarity_;
ComponentBuffer totalMassFrac_;
ComponentBuffer totalMoleFrac_;
ComponentBuffer K_;
ComponentBuffer fugacity_;
PhaseComponentBuffer fugacityCoeff_;
ScalarBuffer L_;
};
} // namespace Opm

4
opm/models/nonlinear/newtonmethod.hh
View File

@ -44,6 +44,7 @@
#include <dune/istl/istlexception.hh>
#include <dune/common/classname.hh>
#include <dune/common/parallel/mpihelper.hh>
//#include <dune/istl/io.hh>
#include <iostream>
#include <sstream>
@ -128,7 +129,7 @@ class NewtonMethod
using Communicator = typename Dune::MPIHelper::MPICommunicator;
using CollectiveCommunication = typename Dune::Communication<typename Dune::MPIHelper::MPICommunicator>;
// using PrintMatrix = typename Dune::printSparseMatrix<typename Dune::>;
public:
NewtonMethod(Simulator& simulator)
: simulator_(simulator)
@ -319,6 +320,7 @@ public:
linearSolver_.getResidual(residual);
solveTimer_.stop();
//Dune::storeMatrixMarket(linearSolver_.overlappingMatrix_, "mymatrix" + ".mm");
// The preSolve_() method usually computes the errors, but it can do
// something else in addition. TODO: should its costs be counted to
// the linearization or to the update?

205
opm/models/ptflash/flashboundaryratevector.hh Normal file
View File

@ -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 <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::FlashBoundaryRateVector
*/
#ifndef EWOMS_FLASH_BOUNDARY_RATE_VECTOR_HH
#define EWOMS_FLASH_BOUNDARY_RATE_VECTOR_HH
#include "flashproperties.hh"
#include <opm/models/common/energymodule.hh>
#include <opm/material/common/Valgrind.hpp>
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 TypeTag>
class FlashBoundaryRateVector : public GetPropType<TypeTag, Properties::RateVector>
{
using ParentType = GetPropType<TypeTag, Properties::RateVector>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
using Toolbox = Opm::MathToolbox<Evaluation>;
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 <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) {
Evaluation& val = this->operator[](eqIdx);
val = Toolbox::min(0.0, val);
}
}
/*!
* \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) {
Evaluation& val = this->operator[](eqIdx);
val = Toolbox::max(0.0, val);
}
}
/*!
* \copydoc ImmiscibleBoundaryRateVector::setNoFlow
*/
void setNoFlow()
{ (*this) = Evaluation(0.0); }
};
} // namespace Opm
#endif

96
opm/models/ptflash/flashextensivequantities.hh Normal file
View File

@ -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 <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::FlashExtensiveQuantities
*/
#ifndef EWOMS_FLASH_EXTENSIVE_QUANTITIES_HH
#define EWOMS_FLASH_EXTENSIVE_QUANTITIES_HH
#include "flashproperties.hh"
#include <opm/models/common/multiphasebaseextensivequantities.hh>
#include <opm/models/common/energymodule.hh>
#include <opm/models/common/diffusionmodule.hh>
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 TypeTag>
class FlashExtensiveQuantities
: public MultiPhaseBaseExtensiveQuantities<TypeTag>
, public EnergyExtensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableEnergy>()>
, public DiffusionExtensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>()>
{
using ParentType = MultiPhaseBaseExtensiveQuantities<TypeTag>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
using DiffusionExtensiveQuantities = Opm::DiffusionExtensiveQuantities<TypeTag, enableDiffusion>;
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using EnergyExtensiveQuantities = Opm::EnergyExtensiveQuantities<TypeTag, enableEnergy>;
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

73
opm/models/ptflash/flashindices.hh Normal file
View File

@ -0,0 +1,73 @@
// -*- 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::FlashIndices
*/
#ifndef EWOMS_FLASH_INDICES_HH
#define EWOMS_FLASH_INDICES_HH
#include "flashproperties.hh"
#include <opm/models/common/energymodule.hh>
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 TypeTag, int PVOffset>
class FlashIndices
: public EnergyIndices<PVOffset + getPropValue<TypeTag, Properties::NumComponents>(),
getPropValue<TypeTag, Properties::EnableEnergy>()>
{
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
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 pressure0Idx = PVOffset;
static const int z0Idx = pressure0Idx+1; // numcomponents-1 variables follow
// equation indices
//! Index of the mass conservation equation for the first
//! component.
static const int conti0EqIdx = PVOffset;
};
} // namespace Opm
#endif

331
opm/models/ptflash/flashintensivequantities.hh Normal file
View File

@ -0,0 +1,331 @@
// -*- 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::FlashIntensiveQuantities
*/
#ifndef EWOMS_FLASH_INTENSIVE_QUANTITIES_HH
#define EWOMS_FLASH_INTENSIVE_QUANTITIES_HH
#include "flashproperties.hh"
#include "flashindices.hh"
#include <opm/models/common/energymodule.hh>
#include <opm/models/common/diffusionmodule.hh>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
namespace Opm {
/*!
* \ingroup FlashModel
* \ingroup IntensiveQuantities
*
* \brief Contains the intensive quantities of the flash-based compositional multi-phase model
*/
template <class TypeTag>
class FlashIntensiveQuantities
: public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
, public DiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
, public EnergyIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableEnergy>() >
, public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
{
using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using FluxModule = GetPropType<TypeTag, Properties::FluxModule>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
// primary variable indices
enum { z0Idx = Indices::z0Idx };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { dimWorld = GridView::dimensionworld };
enum { pressure0Idx = Indices::pressure0Idx };
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using FlashSolver = GetPropType<TypeTag, Properties::FlashSolver>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
using DiffusionIntensiveQuantities = Opm::DiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
using EnergyIntensiveQuantities = Opm::EnergyIntensiveQuantities<TypeTag, enableEnergy>;
public:
//! The type of the object returned by the fluidState() method
using FluidState = Opm::CompositionalFluidState<Evaluation, FluidSystem, enableEnergy>;
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 = 1.e-12;//EWOMS_GET_PARAM(TypeTag, Scalar, FlashTolerance);
int flashVerbosity = 0;//EWOMS_GET_PARAM(TypeTag, int, FlashVerbosity);
std::string flashTwoPhaseMethod = EWOMS_GET_PARAM(TypeTag, std::string, FlashTwoPhaseMethod);
// extract the total molar densities of the components
ComponentVector z(0.);
{
Evaluation lastZ = 1.0;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
z[compIdx] = priVars.makeEvaluation(z0Idx + compIdx, timeIdx);
lastZ -= z[compIdx];
}
z[numComponents - 1] = lastZ;
Evaluation sumz = 0.0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
z[compIdx] = Opm::max(z[compIdx], 1e-8);
sumz +=z[compIdx];
}
z /= sumz;
}
Evaluation p = priVars.makeEvaluation(pressure0Idx, timeIdx);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fluidState_.setPressure(phaseIdx, p);
// Get initial K and L from storage initially (if enabled)
const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx);
const auto *hint2 = elemCtx.thermodynamicHint(dofIdx, 1);
if (hint) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation& Ktmp = hint->fluidState().K(compIdx);
fluidState_.setKvalue(compIdx, Ktmp);
}
const Evaluation& Ltmp = hint->fluidState().L();
fluidState_.setLvalue(Ltmp);
}
else if (hint2) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation& Ktmp = hint2->fluidState().K(compIdx);
fluidState_.setKvalue(compIdx, Ktmp);
}
const Evaluation& Ltmp = hint2->fluidState().L();
fluidState_.setLvalue(Ltmp);
}
else {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation Ktmp = fluidState_.wilsonK_(compIdx);
fluidState_.setKvalue(compIdx, Ktmp);
}
const Evaluation& Ltmp = -1.0;
fluidState_.setLvalue(Ltmp);
}
/////////////
// Compute the phase compositions and densities
/////////////
int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
//FlashSolver::solve(fluidState_, z, spatialIdx, flashVerbosity, flashTwoPhaseMethod, flashTolerance);
//Flash::solve(fluidState_, z, spatialIdx, flashVerbosity, flashTwoPhaseMethod, flashTolerance);
using Flash = Opm::PTFlash<double, FluidSystem>;
FlashSolver::solve(fluidState_, z, spatialIdx, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
//printing of flashresult after solve
// std::cout << " \n After flashsolve for cell " << spatialIdx << std::endl;
// ComponentVector x, y;
// Evaluation L0 = fluidState_.L();
// for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
// x[comp_idx] = fluidState_.moleFraction(FluidSystem::oilPhaseIdx, comp_idx);
// y[comp_idx] = fluidState_.moleFraction(FluidSystem::gasPhaseIdx, comp_idx);
// }
// for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
// std::cout << " x for component: " << comp_idx << " is " << x[comp_idx] << std::endl;
// for (int i = 0; i < 3; ++i) {
// std::cout << " x deriv " << i << " is: " << x[comp_idx].derivative(i) << std::endl;
// }
// std::cout << " y for component: " << comp_idx << "is " << y[comp_idx] << std::endl;
// for (int i = 0; i < 3; ++i) {
// std::cout << " y deriv " << i << " is: " << y[comp_idx].derivative(i) << std::endl;
// }
// }
// std::cout << " L is " << L0 << std::endl;
// for (int i = 0; i < L0.size(); ++i) {
// std::cout << " L deriv " << i << " is: " << L0.derivative(i) << std::endl;
// }
// //end printinting 1
// Update phases
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx);
const Scalar R = Opm::Constants<Scalar>::R;
Evaluation Z_L = (paramCache.molarVolume(FluidSystem::oilPhaseIdx) * fluidState_.pressure(FluidSystem::oilPhaseIdx) )/
(R * fluidState_.temperature(FluidSystem::oilPhaseIdx));
paramCache.updatePhase(fluidState_, FluidSystem::gasPhaseIdx);
Evaluation Z_V = (paramCache.molarVolume(FluidSystem::gasPhaseIdx) * fluidState_.pressure(FluidSystem::gasPhaseIdx) )/
(R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
// Update saturation
Evaluation L = fluidState_.L();
Evaluation So = Opm::max((L*Z_L/(L*Z_L+(1-L)*Z_V)), 0.0);
Evaluation Sg = Opm::max(1-So, 0.0);
Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg);
So /= sumS;
Sg /= sumS;
fluidState_.setSaturation(0, So);
fluidState_.setSaturation(1, Sg);
fluidState_.setCompressFactor(0, Z_L);
fluidState_.setCompressFactor(1, Z_V);
// Print saturation
if (flashVerbosity == 5) {
std::cout << "So = " << So <<std::endl;
std::cout << "Sg = " << Sg <<std::endl;
}
// Print saturation
if (flashVerbosity == 5) {
std::cout << "So = " << So <<std::endl;
std::cout << "Sg = " << Sg <<std::endl;
std::cout << "Z_L = " << Z_L <<std::endl;
std::cout << "Z_V = " << Z_V <<std::endl;
}
/////////////
// Compute rel. perm and viscosities and densities
/////////////
const MaterialLawParams& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
// calculate relative permeabilities
MaterialLaw::relativePermeabilities(relativePermeability_,
materialParams, fluidState_);
Opm::Valgrind::CheckDefined(relativePermeability_);
// set the phase viscosities and density
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
paramCache.updatePhase(fluidState_, phaseIdx);
const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
fluidState_.setViscosity(phaseIdx, mu);
//fluidState_.setViscosity(phaseIdx, newmu);
mobility_[phaseIdx] = relativePermeability_[phaseIdx] / mu;
Opm::Valgrind::CheckDefined(mobility_[phaseIdx]);
const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, phaseIdx);
fluidState_.setDensity(phaseIdx, rho);
}
/////////////
// Compute 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

199
opm/models/ptflash/flashlocalresidual.hh Normal file
View File

@ -0,0 +1,199 @@
// -*- 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::FlashLocalResidual
*/
#ifndef EWOMS_FLASH_LOCAL_RESIDUAL_HH
#define EWOMS_FLASH_LOCAL_RESIDUAL_HH
#include "flashproperties.hh"
#include <opm/models/common/diffusionmodule.hh>
#include <opm/models/common/energymodule.hh>
#include <opm/material/common/Valgrind.hpp>
namespace Opm {
/*!
* \ingroup FlashModel
*
* \brief Calculates the local residual of the compositional multi-phase
* model based on flash calculations.
*/
template <class TypeTag>
class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResidual>
{
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
using DiffusionModule = Opm::DiffusionModule<TypeTag, enableDiffusion>;
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
using Toolbox = Opm::MathToolbox<Evaluation>;
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& 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<LhsEval>(fs.massFraction(phaseIdx, compIdx))
* Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.porosity());
}
EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
}
/*!
* \copydoc FvBaseLocalResidual::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 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<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().density(phaseIdx)
* extQuants.volumeFlux(phaseIdx);
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
flux[conti0EqIdx + compIdx] +=
tmp*up.fluidState().massFraction(phaseIdx, compIdx);
}
}
else {
Evaluation tmp =
Toolbox::value(up.fluidState().density(phaseIdx)
* extQuants.volumeFlux(phaseIdx));
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
flux[conti0EqIdx + compIdx] +=
tmp*Toolbox::value(up.fluidState().massFraction(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

355
opm/models/ptflash/flashmodel.hh Normal file
View File

@ -0,0 +1,355 @@
// -*- 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::FlashModel
*/
#ifndef EWOMS_FLASH_MODEL_HH
#define EWOMS_FLASH_MODEL_HH
#include <opm/material/densead/Math.hpp>
#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 "flashnewtonmethod.hh"
#include <opm/models/common/multiphasebasemodel.hh>
#include <opm/models/common/energymodule.hh>
#include <opm/models/io/vtkcompositionmodule.hh>
#include <opm/models/io/vtkenergymodule.hh>
#include <opm/models/io/vtkdiffusionmodule.hh>
#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
#include <opm/material/constraintsolvers/PTFlash.hpp>
#include <sstream>
#include <string>
namespace Opm {
template <class TypeTag>
class FlashModel;
}
namespace Opm::Properties {
namespace TTag {
//! The type tag for the isothermal single phase problems
struct FlashModel { using InheritsFrom = std::tuple<VtkDiffusion,
VtkEnergy,
VtkComposition,
MultiPhaseBaseModel>; };
} // namespace TTag
//! Use the FlashLocalResidual function for the flash model
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::FlashModel> { using type = Opm::FlashLocalResidual<TypeTag>; };
//! Use the PT flash specific newton method for the flash model
template<class TypeTag>
struct NewtonMethod<TypeTag, TTag::FlashModel> { using type = Opm::FlashNewtonMethod<TypeTag>; };
//! Use the Pt flash solver by default
template<class TypeTag>
struct FlashSolver<TypeTag, TTag::FlashModel>
{ using type = Opm::PTFlash<GetPropType<TypeTag, Properties::Scalar>,
GetPropType<TypeTag, Properties::FluidSystem>>; };
//! Let the flash solver choose its tolerance by default
template<class TypeTag>
struct FlashTolerance<TypeTag, TTag::FlashModel>
{
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = -1.0;
};
// Flash solver verbosity
template<class TypeTag>
struct FlashVerbosity<TypeTag, TTag::FlashModel> { static constexpr int value = 0; };
// Flash two-phase method
template<class TypeTag>
struct FlashTwoPhaseMethod<TypeTag, TTag::FlashModel> { static constexpr auto value = "ssi"; };
//! the Model property
template<class TypeTag>
struct Model<TypeTag, TTag::FlashModel> { using type = Opm::FlashModel<TypeTag>; };
//! the PrimaryVariables property
template<class TypeTag>
struct PrimaryVariables<TypeTag, TTag::FlashModel> { using type = Opm::FlashPrimaryVariables<TypeTag>; };
//! the RateVector property
template<class TypeTag>
struct RateVector<TypeTag, TTag::FlashModel> { using type = Opm::FlashRateVector<TypeTag>; };
//! the BoundaryRateVector property
template<class TypeTag>
struct BoundaryRateVector<TypeTag, TTag::FlashModel> { using type = Opm::FlashBoundaryRateVector<TypeTag>; };
//! the IntensiveQuantities property
template<class TypeTag>
struct IntensiveQuantities<TypeTag, TTag::FlashModel> { using type = Opm::FlashIntensiveQuantities<TypeTag>; };
//! the ExtensiveQuantities property
template<class TypeTag>
struct ExtensiveQuantities<TypeTag, TTag::FlashModel> { using type = Opm::FlashExtensiveQuantities<TypeTag>; };
//! The indices required by the flash-baseed isothermal compositional model
template<class TypeTag>
struct Indices<TypeTag, TTag::FlashModel> { using type = Opm::FlashIndices<TypeTag, /*PVIdx=*/0>; };
// The updates of intensive quantities tend to be _very_ expensive for this
// model, so let's try to minimize the number of required ones
template<class TypeTag>
struct EnableIntensiveQuantityCache<TypeTag, TTag::FlashModel> { static constexpr bool value = 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.
template<class TypeTag>
struct EnableThermodynamicHints<TypeTag, TTag::FlashModel> { static constexpr bool value = true; };
// disable molecular diffusion by default
template<class TypeTag>
struct EnableDiffusion<TypeTag, TTag::FlashModel> { static constexpr bool value = false; };
//! Disable the energy equation by default
template<class TypeTag>
struct EnableEnergy<TypeTag, TTag::FlashModel> { static constexpr bool value = false; };
} // namespace Opm::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
* template<class TypeTag>
struct FluxModule<TypeTag, TTag::MyProblemTypeTag> { using type = 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]
*
* To determine the quanties that occur in the equations above, this
* model uses <i>flash calculations</i>. 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 TypeTag>
class FlashModel
: public MultiPhaseBaseModel<TypeTag>
{
using ParentType = MultiPhaseBaseModel<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
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<TypeTag>::registerParameters();
if (enableDiffusion)
Opm::VtkDiffusionModule<TypeTag>::registerParameters();
if (enableEnergy)
Opm::VtkEnergyModule<TypeTag>::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, Scalar, FlashTolerance,
"The maximum tolerance for the flash solver to "
"consider the solution converged");
EWOMS_REGISTER_PARAM(TypeTag, int, FlashVerbosity,
"Flash solver verbosity level");
EWOMS_REGISTER_PARAM(TypeTag, std::string, FlashTwoPhaseMethod,
"Method for solving vapor-liquid composition");
}
/*!
* \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::z0Idx <= pvIdx && pvIdx < Indices::z0Idx + numComponents - 1)
oss << "z_," << FluidSystem::componentName(/*compIdx=*/pvIdx - Indices::z0Idx);
else if (pvIdx==Indices::pressure0Idx)
oss << "pressure_" << FluidSystem::phaseName(0);
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::conti0EqIdx;
// 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<TypeTag>(this->simulator_));
if (enableDiffusion)
this->addOutputModule(new Opm::VtkDiffusionModule<TypeTag>(this->simulator_));
if (enableEnergy)
this->addOutputModule(new Opm::VtkEnergyModule<TypeTag>(this->simulator_));
}
};
} // namespace Opm
#endif

144
opm/models/ptflash/flashnewtonmethod.hh Normal file
View File

@ -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 <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::FlashNewtonMethod
*/
#ifndef EWOMS_FLASH_NEWTON_METHOD_HH
#define EWOMS_FLASH_NEWTON_METHOD_HH
#include "flashproperties.hh"
#include <opm/models/nonlinear/newtonmethod.hh>
#include <opm/common/Exceptions.hpp>
#include <algorithm>
namespace Opm::Properties {
template <class TypeTag, class MyTypeTag>
struct DiscNewtonMethod;
} // namespace Opm::Properties
namespace Opm {
/*!
* \ingroup FlashModel
*
* \brief A Newton solver specific to the NCP model.
*/
template <class TypeTag>
class FlashNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMethod>
{
using ParentType = GetPropType<TypeTag, Properties::DiscNewtonMethod>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
enum { pressure0Idx = Indices::pressure0Idx };
enum { z0Idx = Indices::z0Idx };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
public:
/*!
* \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& )
*/
FlashNewtonMethod(Simulator& simulator) : ParentType(simulator)
{}
protected:
friend ParentType;
friend NewtonMethod<TypeTag>;
/*!
* \copydoc FvBaseNewtonMethod::updatePrimaryVariables_
*/
void updatePrimaryVariables_(unsigned globalDofIdx,
PrimaryVariables& nextValue,
const PrimaryVariables& currentValue,
const EqVector& update,
const EqVector& currentResidual)
{
// normal Newton-Raphson update
nextValue = currentValue;
nextValue -= update;
////
// Pressure updates
////
// limit pressure reference change to 20% of the total value per iteration
clampValue_(nextValue[pressure0Idx],
currentValue[pressure0Idx]*0.8,
currentValue[pressure0Idx]*1.2);
////
// z updates
////
// restrict update to at most 0.1
Scalar maxDeltaZ = 0.0; // in update vector
Scalar sumDeltaZ = 0.0; // changes in last component (not in update vector)
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
maxDeltaZ = std::max(std::abs(update[z0Idx + compIdx]), maxDeltaZ);
sumDeltaZ += update[z0Idx + compIdx];
}
maxDeltaZ = std::max(std::abs(-sumDeltaZ), maxDeltaZ);
// if max. update is above 0.1, restrict that one to 0.1 and adjust the rest accordingly (s.t. last comp. update is sum of the changes in update vector)
if (maxDeltaZ > 0.1) {
Scalar alpha = 0.2/maxDeltaZ;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
nextValue[z0Idx + compIdx] = currentValue[z0Idx + compIdx] - alpha*update[z0Idx + compIdx];
}
}
// ensure that z-values are less than tol or more than 1-tol
Scalar tol = 1e-8;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol);
}
// normalize s.t. z sums to 1.0
// Scalar lastZ = 1.0;
// Scalar sumZ = 0.0;
// for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
// lastZ -= nextValue[z0Idx + compIdx];
// sumZ += nextValue[z0Idx + compIdx];
// }
// sumZ += lastZ;
// for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
// nextValue[z0Idx + compIdx] /= sumZ;
// }
}
private:
void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const
{ val = std::max(minVal, std::min(val, maxVal)); }
}; // class FlashNewtonMethod
} // namespace Opm
#endif

160
opm/models/ptflash/flashprimaryvariables.hh Normal file
View File

@ -0,0 +1,160 @@
// -*- 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::FlashPrimaryVariables
*/
#ifndef EWOMS_FLASH_PRIMARY_VARIABLES_HH
#define EWOMS_FLASH_PRIMARY_VARIABLES_HH
#include "flashindices.hh"
#include "flashproperties.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/common/Valgrind.hpp>
#include <dune/common/fvector.hh>
#include <iostream>
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 TypeTag>
class FlashPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
{
using ParentType = FvBasePrimaryVariables<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
// primary variable indices
enum { z0Idx = Indices::z0Idx };
enum { pressure0Idx = Indices::pressure0Idx };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
using Toolbox = typename Opm::MathToolbox<Evaluation>;
using EnergyModule = Opm::EnergyModule<TypeTag, getPropValue<TypeTag, Properties::EnableEnergy>()>;
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 <class FluidState>
void assignMassConservative(const FluidState& fluidState,
const MaterialLawParams&,
bool = 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 <class FluidState>
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.
Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar tmp = Opm::getValue(fluidState.molarity(phaseIdx, compIdx) * fluidState.saturation(phaseIdx));
z[compIdx] += Opm::max(tmp, 1e-8);
sumMoles += tmp;
}
}
z /= sumMoles;
for (int i = 0; i < numComponents - 1; ++i)
(*this)[z0Idx + i] = z[i];
(*this)[pressure0Idx] = Opm::getValue(fluidState.pressure(0));
}
/*!
* \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
{
os << "(p_" << FluidSystem::phaseName(0) << " = "
<< this->operator[](pressure0Idx);
for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) {
os << ", z_" << FluidSystem::componentName(compIdx) << " = "
<< this->operator[](z0Idx + compIdx);
}
os << ")" << std::flush;
}
};
} // namespace Opm
#endif

55
opm/models/ptflash/flashproperties.hh Normal file
View File

@ -0,0 +1,55 @@
// -*- 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 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 <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/io/vtkcompositionmodule.hh>
#include <opm/models/io/vtkenergymodule.hh>
#include <opm/models/io/vtkdiffusionmodule.hh>
namespace Opm::Properties {
//! The type of the flash constraint solver
template<class TypeTag, class MyTypeTag>
struct FlashSolver { using type = UndefinedProperty; };
//! The maximum accepted error of the flash solver
template<class TypeTag, class MyTypeTag>
struct FlashTolerance { using type = UndefinedProperty; };
//! The verbosity level of the flash solver
template<class TypeTag, class MyTypeTag>
struct FlashVerbosity { using type = UndefinedProperty; };
//! Two-phase flash method
template<class TypeTag, class MyTypeTag>
struct FlashTwoPhaseMethod { using type = UndefinedProperty; };
} // namespace Opm::Properties
#endif

143
opm/models/ptflash/flashratevector.hh Normal file
View File

@ -0,0 +1,143 @@
// -*- 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::FlashRateVector
*/
#ifndef EWOMS_FLASH_RATE_VECTOR_HH
#define EWOMS_FLASH_RATE_VECTOR_HH
#include <dune/common/fvector.hh>
#include <opm/models/common/energymodule.hh>
#include <opm/material/common/Valgrind.hpp>
#include "flashintensivequantities.hh"
namespace Opm {
/*!
* \ingroup FlashModel
*
* \copydoc Opm::ImmiscibleRateVector
*/
template <class TypeTag>
class FlashRateVector
: public Dune::FieldVector<GetPropType<TypeTag, Properties::Evaluation>,
getPropValue<TypeTag, Properties::NumEq>()>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
using ParentType = Dune::FieldVector<Evaluation, numEq>;
using EnergyModule = Opm::EnergyModule<TypeTag, getPropValue<TypeTag, Properties::EnableEnergy>()>;
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 <class FluidState, class RhsEval>
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 <class RhsEval>
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

6
opm/simulators/linalg/parallelbasebackend.hh
View File

@ -37,6 +37,8 @@
#include <opm/simulators/linalg/overlappingoperator.hh>
#include <opm/simulators/linalg/parallelbasebackend.hh>
#include <opm/simulators/linalg/istlpreconditionerwrappers.hh>
//#include <MatrixMarketSpecializations.hpp>
//#include <dune/istl/matrixmarket.hh>
#include <opm/models/utils/genericguard.hh>
#include <opm/models/utils/propertysystem.hh>
@ -261,6 +263,9 @@ public:
ParallelScalarProduct parScalarProduct(overlappingMatrix_->overlap());
ParallelOperator parOperator(*overlappingMatrix_);
//const auto& matrix = parOperator.getMatrix();
//Dune::storeMatrixMarket(*overlappingMatrix_, std::string("mymatrix.mm"));
// retrieve the linear solver
auto solver = asImp_().prepareSolver_(parOperator,
parScalarProduct,
@ -278,6 +283,7 @@ public:
// copy the result back to the non-overlapping vector
overlappingx_->assignTo(x);
overlappingx_->print();
// return the result of the solver
return result.first;

3
opm/simulators/linalg/parallelbicgstabbackend.hh
View File

@ -156,6 +156,9 @@ protected:
bicgstabSolver->setMaxIterations(EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIterations));
bicgstabSolver->setLinearOperator(&parOperator);
bicgstabSolver->setRhs(this->overlappingb_);
// const auto& matrix = parOperator.getMatrix();
// Dune::storeMatrixMarket(matrix, "mymatrix" + ".mm");
return bicgstabSolver;
}

64
tests/simpletest_ptflash_ecfv.cc Normal file
View File

@ -0,0 +1,64 @@
// -*- 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
*
* \brief Box problem with two phases and multiple components
*/
#include "config.h"
#include <opm/models/utils/start.hh>
//#include <opm/models/immiscible/immisciblemodel.hh>
//#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include "problems/simpletestproblem.hh"
namespace Opm::Properties {
namespace TTag {
struct SimpleTestEcfvProblem
{
using InheritsFrom = std::tuple<SimpleTest, FlashModel>;
};
}
template <class TypeTag>
struct SpatialDiscretizationSplice<TypeTag, TTag::SimpleTestEcfvProblem>
{
using type = TTag::EcfvDiscretization;
};
//TESTAD
template <class TypeTag>
struct LocalLinearizerSplice<TypeTag, TTag::SimpleTestEcfvProblem>
{
using type = TTag::AutoDiffLocalLinearizer;
};
} // namespace Opm::Properties
int main(int argc, char **argv)
{
using EcfvProblemTypeTag = Opm::Properties::TTag::SimpleTestEcfvProblem;
return Opm::start<EcfvProblemTypeTag>(argc, argv);
}