OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-22 15:33:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

various cleaning up and refactoring for co2 ptflash simulation

Browse Source
This commit is contained in:
Kai Bao 2023-10-03 10:00:43 +02:00
parent df5ee732fa
commit 9f0fb32713
24 changed files with 427 additions and 929 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
tests/simpletest_ptflash_ecfv.cc → examples/co2_ptflash_ecfv.cpp
View File

@ -23,33 +23,31 @@
/*!
* \file
*
* \brief Box problem with two phases and multiple components
* \brief Box problem with two phases and multiple components.
* Solved with a PTFlash two phase solver.
*/
#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"
#include "problems/co2ptflashproblem.hh"
namespace Opm::Properties {
namespace TTag {
struct SimpleTestEcfvProblem
{
using InheritsFrom = std::tuple<SimpleTest, FlashModel>;
struct CO2PTEcfvProblem {
using InheritsFrom = std::tuple<CO2PTBaseProblem, FlashModel>;
};
}
template <class TypeTag>
struct SpatialDiscretizationSplice<TypeTag, TTag::SimpleTestEcfvProblem>
struct SpatialDiscretizationSplice<TypeTag, TTag::CO2PTEcfvProblem>
{
using type = TTag::EcfvDiscretization;
};
//TESTAD
template <class TypeTag>
struct LocalLinearizerSplice<TypeTag, TTag::SimpleTestEcfvProblem>
struct LocalLinearizerSplice<TypeTag, TTag::CO2PTEcfvProblem>
{
using type = TTag::AutoDiffLocalLinearizer;
};
@ -59,6 +57,6 @@ struct LocalLinearizerSplice<TypeTag, TTag::SimpleTestEcfvProblem>
int main(int argc, char **argv)
{
using EcfvProblemTypeTag = Opm::Properties::TTag::SimpleTestEcfvProblem;
using EcfvProblemTypeTag = Opm::Properties::TTag::CO2PTEcfvProblem;
return Opm::start<EcfvProblemTypeTag>(argc, argv);
}

66
examples/problems/co2injectiontestflash.hh
View File

@ -1,66 +0,0 @@
// -*- 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

284
examples/problems/simpletestproblem.hh → examples/problems/co2ptflashproblem.hh
View File

@ -23,10 +23,10 @@
/*!
* \file
*
* \copydoc Opm::simpletestproblem
* \copydoc Opm::co2ptflashproblem
*/
#ifndef EWOMS_SIMPLETEST_PROBLEM_HH
#define EWOMS_SIMPLETEST_PROBLEM_HH
#ifndef OPM_CO2PTFLASH_PROBLEM_HH
#define OPM_CO2PTFLASH_PROBLEM_HH
#include <opm/common/Exceptions.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
@ -54,16 +54,15 @@
namespace Opm {
template <class TypeTag>
class SimpleTest;
} // namespace Opm
class CO2PTProblem;
} // namespace Opm */
namespace Opm::Properties {
namespace TTag {
struct SimpleTest {};
struct CO2PTBaseProblem {};
} // 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>
@ -74,18 +73,18 @@ struct EpisodeLength { using type = UndefinedProperty;};
template <class TypeTag, class MyTypeTag>
struct Initialpressure { using type = UndefinedProperty;};
// Set the grid type: --->1D
// Set the grid type: --->2D
template <class TypeTag>
struct Grid<TypeTag, TTag::SimpleTest> { using type = Dune::YaspGrid</*dim=*/2>; };
struct Grid<TypeTag, TTag::CO2PTBaseProblem> { using type = Dune::YaspGrid</*dim=*/2>; };
// Set the problem property
template <class TypeTag>
struct Problem<TypeTag, TTag::SimpleTest>
{ using type = Opm::SimpleTest<TypeTag>; };
struct Problem<TypeTag, TTag::CO2PTBaseProblem>
{ using type = Opm::CO2PTProblem<TypeTag>; };
// Set flash solver
template <class TypeTag>
struct FlashSolver<TypeTag, TTag::SimpleTest> {
struct FlashSolver<TypeTag, TTag::CO2PTBaseProblem> {
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@ -97,7 +96,7 @@ public:
// Set fluid configuration
template <class TypeTag>
struct FluidSystem<TypeTag, TTag::SimpleTest>
struct FluidSystem<TypeTag, TTag::CO2PTBaseProblem>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@ -108,8 +107,7 @@ public:
// Set the material Law
template <class TypeTag>
struct MaterialLaw<TypeTag, TTag::SimpleTest>
{
struct MaterialLaw<TypeTag, TTag::CO2PTBaseProblem> {
private:
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
@ -117,12 +115,11 @@ private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
// /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, TODO
// /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, // TODO
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
// define the material law which is parameterized by effective saturations
// define the material law which is parameterized by effective saturation
using EffMaterialLaw = Opm::NullMaterial<Traits>;
//using EffMaterialLaw = Opm::BrooksCorey<Traits>;
@ -132,129 +129,122 @@ public:
// Write the Newton convergence behavior to disk?
template <class TypeTag>
struct NewtonWriteConvergence<TypeTag, TTag::SimpleTest> {
struct NewtonWriteConvergence<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = false; };
// Enable gravity false
template <class TypeTag>
struct EnableGravity<TypeTag, TTag::SimpleTest> { static constexpr bool value = false;
struct EnableGravity<TypeTag, TTag::CO2PTBaseProblem> { static constexpr bool value = false;
};
// set the defaults for the problem specific properties
template <class TypeTag>
struct Temperature<TypeTag, TTag::SimpleTest> {
struct Temperature<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 423.25;//TODO
};
template <class TypeTag>
struct Initialpressure<TypeTag, TTag::SimpleTest> {
struct Initialpressure<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e5;
};
template <class TypeTag>
struct SimulationName<TypeTag, TTag::SimpleTest> {
static constexpr auto value = "simpletest";
struct SimulationName<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr auto value = "co2_ptflash";
};
// The default for the end time of the simulation
template <class TypeTag>
struct EndTime<TypeTag, TTag::SimpleTest> {
struct EndTime<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 3 * 24. * 60. * 60.;//3 days
static constexpr type value = 60. * 60.;
};
// convergence control
template <class TypeTag>
struct InitialTimeStepSize<TypeTag, TTag::SimpleTest> {
struct InitialTimeStepSize<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
//static constexpr type value = 30;
static constexpr type value = 1 * 24. * 60. * 60.;
static constexpr type value = 0.1 * 60. * 60.;
};
template <class TypeTag>
struct LinearSolverTolerance<TypeTag, TTag::SimpleTest> {
struct LinearSolverTolerance<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e-3;
};
template <class TypeTag>
struct LinearSolverAbsTolerance<TypeTag, TTag::SimpleTest> {
struct LinearSolverAbsTolerance<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 0.;
};
template <class TypeTag>
struct NewtonTolerance<TypeTag, TTag::SimpleTest> {
struct NewtonTolerance<TypeTag, TTag::CO2PTBaseProblem> {
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> {
struct NewtonMaxIterations<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 30;
};
template <class TypeTag>
struct NewtonTargetIterations<TypeTag, TTag::SimpleTest> {
struct NewtonTargetIterations<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 6;
};
// output
template <class TypeTag>
struct VtkWriteFilterVelocities<TypeTag, TTag::SimpleTest> {
struct VtkWriteFilterVelocities<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWritePotentialGradients<TypeTag, TTag::SimpleTest> {
struct VtkWritePotentialGradients<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteTotalMassFractions<TypeTag, TTag::SimpleTest> {
struct VtkWriteTotalMassFractions<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteTotalMoleFractions<TypeTag, TTag::SimpleTest> {
struct VtkWriteTotalMoleFractions<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteFugacityCoeffs<TypeTag, TTag::SimpleTest> {
struct VtkWriteFugacityCoeffs<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteLiquidMoleFractions<TypeTag, TTag::SimpleTest> {
struct VtkWriteLiquidMoleFractions<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
template <class TypeTag>
struct VtkWriteEquilibriumConstants<TypeTag, TTag::SimpleTest> {
struct VtkWriteEquilibriumConstants<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
// write restart for every hour
// this is kinds of telling the report step length
template <class TypeTag>
struct EpisodeLength<TypeTag, TTag::SimpleTest> {
struct EpisodeLength<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 60. * 60.;
static constexpr type value = 0.1 * 60. * 60.;
};
// mesh grid
template <class TypeTag>
struct Vanguard<TypeTag, TTag::SimpleTest> {
struct Vanguard<TypeTag, TTag::CO2PTBaseProblem> {
using type = Opm::StructuredGridVanguard<TypeTag>;
};
@ -262,34 +252,36 @@ struct Vanguard<TypeTag, TTag::SimpleTest> {
//\Note: DomainSizeX is 3.0 meters
//\Note: DomainSizeY is 1.0 meters
template <class TypeTag>
struct DomainSizeX<TypeTag, TTag::SimpleTest> {
struct DomainSizeX<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 3; // meter
static constexpr type value = 300; // meter
};
template <class TypeTag>
struct DomainSizeY<TypeTag, TTag::SimpleTest> {
struct DomainSizeY<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0;
};
// DomainSizeZ is not needed, while to keep structuredgridvanguard.hh compile
template <class TypeTag>
struct DomainSizeZ<TypeTag, TTag::SimpleTest> {
struct DomainSizeZ<TypeTag, TTag::CO2PTBaseProblem> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0;
};
template<class TypeTag>
struct CellsX<TypeTag, TTag::SimpleTest> { static constexpr int value = 3; };
struct CellsX<TypeTag, TTag::CO2PTBaseProblem> { static constexpr int value = 30; };
template<class TypeTag>
struct CellsY<TypeTag, TTag::SimpleTest> { static constexpr int value = 1; };
struct CellsY<TypeTag, TTag::CO2PTBaseProblem> { static constexpr int value = 1; };
// CellsZ is not needed, while to keep structuredgridvanguard.hh compile
template<class TypeTag>
struct CellsZ<TypeTag, TTag::SimpleTest> { static constexpr int value = 1; };
struct CellsZ<TypeTag, TTag::CO2PTBaseProblem> { static constexpr int value = 1; };
// compositional, with diffusion
template <class TypeTag>
struct EnableEnergy<TypeTag, TTag::SimpleTest> {
struct EnableEnergy<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = false;
};
@ -306,7 +298,7 @@ namespace Opm {
*
*/
template <class TypeTag>
class SimpleTest : public GetPropType<TypeTag, Properties::BaseProblem>
class CO2PTProblem : public GetPropType<TypeTag, Properties::BaseProblem>
{
using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
@ -352,10 +344,10 @@ public:
/*!
* \copydoc Doxygen::defaultProblemConstructor
*/
SimpleTest(Simulator& simulator)
explicit CO2PTProblem(Simulator& simulator)
: ParentType(simulator)
{
Scalar epi_len = EWOMS_GET_PARAM(TypeTag, Scalar, EpisodeLength);
const Scalar epi_len = EWOMS_GET_PARAM(TypeTag, Scalar, EpisodeLength);
simulator.setEpisodeLength(epi_len);
}
@ -391,7 +383,7 @@ public:
}
/*!
* \copydoc FvBaseMultiPhaseProblem::registerParameters
* \copydoc co2ptflashproblem::registerParameters
*/
static void registerParameters()
{
@ -476,7 +468,6 @@ public:
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
@ -502,11 +493,12 @@ public:
{
int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
int inj = 0;
int prod = 2;
if (spatialIdx == inj || spatialIdx == prod)
int prod = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX) - 1;
if (spatialIdx == inj || spatialIdx == prod) {
return 1.0;
else
} else {
return porosity_;
}
}
/*!
@ -547,113 +539,75 @@ private:
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>;
// 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 = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX) - 1;
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];
const Scalar temp = 423.25;
constexpr auto numComponents = FluidSystem::numComponents;
using Evaluation = Opm::DenseAd::Evaluation<double, numComponents>;
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
Scalar p0 = 75e5; // CONVERGENCE FAILURE WITH 75
// 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;
//\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);
// TODO: no capillary pressure for now
Scalar p0 = 75e5; //CONVERGENCE FAILURE WITH 75
fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
//\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.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, 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.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
fs.setTemperature(temp);
// It is used here only for calculate the z
fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
// 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));
}
fs.setTemperature(temp);
// Set initial K and L
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation Ktmp = fs.wilsonK_(compIdx);
fs.setKvalue(compIdx, Ktmp);
}
// 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);
const Evaluation& Ltmp = -1.0;
fs.setLvalue(Ltmp);
}
DimMatrix K_;
@ -664,4 +618,4 @@ private:
};
} // namespace Opm
#endif
#endif

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

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

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

@ -273,26 +273,18 @@ 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 = true; };
struct EnableIntensiveQuantityCache<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
// 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 = true; };
struct EnableThermodynamicHints<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
// 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 not default
// use volumetric residuals is default
template<class TypeTag>
struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = true; };
//! 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[1] = model().thermodynamicHint(globalIdx, 1);
dofVars_[dofIdx].thermodynamicHint[timeIdx] =
model().thermodynamicHint(globalIdx, timeIdx);
const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
if (cachedIntQuants) {

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

@ -490,8 +490,6 @@ 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
@ -499,7 +497,6 @@ 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,8 +161,6 @@ 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);
@ -326,10 +324,8 @@ 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,10 +59,6 @@ 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>
@ -79,10 +75,6 @@ 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
@ -120,7 +112,6 @@ 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)
@ -146,10 +137,6 @@ 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");
}
/*!
@ -168,10 +155,7 @@ 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_())
@ -194,9 +178,6 @@ 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_())
@ -241,10 +222,6 @@ public:
}
if (fugacityOutput_())
fugacity_[compIdx][I] = Toolbox::value(intQuants.fluidState().fugacity(/*phaseIdx=*/0, compIdx));
if (equilConstOutput_())
K_[compIdx][I] = Toolbox::value(fs.K(compIdx));
}
}
}
@ -259,7 +236,7 @@ public:
return;
}
if (moleFracOutput_())
if (moleFracOutput_())
this->commitPhaseComponentBuffer_(baseWriter, "moleFrac_%s^%s", moleFrac_);
if (massFracOutput_())
this->commitPhaseComponentBuffer_(baseWriter, "massFrac_%s^%s", massFrac_);
@ -269,15 +246,11 @@ 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:
@ -322,30 +295,15 @@ 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

181
opm/models/io/vtkptflashmodule.hh Normal file
View File

@ -0,0 +1,181 @@
// -*- 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::VtkPTFlashModule
*/
#ifndef OPM_VTK_PTFLASH_MODULE_HH
#define OPM_VTK_PTFLASH_MODULE_HH
#include "vtkmultiwriter.hh"
#include "baseoutputmodule.hh"
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/parametersystem.hh>
#include <opm/material/common/MathToolbox.hpp>
namespace Opm::Properties {
namespace TTag {
// create new type tag for the VTK PTFlash output
struct VtkPTFlash {};
} // namespace TTag
// create the property tags needed for the composition module
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>
struct VtkWriteLiquidMoleFractions<TypeTag, TTag::VtkPTFlash> { static constexpr bool value = false; };
template<class TypeTag>
struct VtkWriteEquilibriumConstants<TypeTag, TTag::VtkPTFlash> { static constexpr bool value = false; };
} // namespace Opm::Properties
namespace Opm {
/*!
* \ingroup Vtk
*
* \brief VTK output module for the PT Flash calculation
* This module deals with the following quantities:
* K, equilibrium ratio for all the components
* L, liquid fraction in the two-phase system
*/
template <class TypeTag>
class VtkPTFlashModule: public BaseOutputModule<TypeTag>
{
using ParentType = BaseOutputModule<TypeTag>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
static const int vtkFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
using VtkMultiWriter = ::Opm::VtkMultiWriter<GridView, vtkFormat>;
using ComponentBuffer = typename ParentType::ComponentBuffer;
using ScalarBuffer = typename ParentType::ScalarBuffer;
public:
explicit VtkPTFlashModule(const Simulator& simulator)
: ParentType(simulator)
{ }
/*!
* \brief Register all run-time parameters for the Vtk output module.
*/
static void registerParameters()
{
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");
}
/*!
* \brief Allocate memory for the scalar fields we would like to
* write to the VTK file.
*/
void allocBuffers()
{
if (LOutput_())
this->resizeScalarBuffer_(L_);
if (equilConstOutput_())
this->resizeComponentBuffer_(K_);
}
/*!
* \brief Modify the internal buffers according to the intensive quantities relevant
* for an element
*/
void processElement(const ElementContext& elemCtx)
{
using Toolbox = MathToolbox<Evaluation>;
if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput))
return;
for (unsigned i = 0; i < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++i) {
unsigned I = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
if (LOutput_())
L_[I] = Toolbox::value(fs.L());
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
if (equilConstOutput_())
K_[compIdx][I] = Toolbox::value(fs.K(compIdx));
}
}
}
/*!
* \brief Add all buffers to the VTK output writer.
*/
void commitBuffers(BaseOutputWriter& baseWriter)
{
auto *vtkWriter = dynamic_cast<VtkMultiWriter*>(&baseWriter);
if (!vtkWriter) {
return;
}
if (equilConstOutput_())
this->commitComponentBuffer_(baseWriter, "K^%s", K_);
if (LOutput_())
this->commitScalarBuffer_(baseWriter, "L", L_);
}
private:
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;
}
ComponentBuffer K_;
ScalarBuffer L_;
};
} // namespace Opm
#endif

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

@ -44,7 +44,6 @@
#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>
@ -129,7 +128,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)
@ -320,7 +319,6 @@ 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
View File

@ -1,205 +0,0 @@
// -*- 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
View File

@ -1,96 +0,0 @@
// -*- 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

16
opm/models/ptflash/flashindices.hh
View File

@ -25,10 +25,9 @@
*
* \copydoc Opm::FlashIndices
*/
#ifndef EWOMS_FLASH_INDICES_HH
#define EWOMS_FLASH_INDICES_HH
#ifndef OPM_PTFLASH_INDICES_HH
#define OPM_PTFLASH_INDICES_HH
#include "flashproperties.hh"
#include <opm/models/common/energymodule.hh>
namespace Opm {
@ -37,7 +36,7 @@ namespace Opm {
* \ingroup FlashModel
*
* \brief Defines the primary variable and equation indices for the
* compositional multi-phase model based on flash calculations.
* compositional multi-phase model based on PT flash calculations.
*
* \tparam PVOffset The first index in a primary variable vector.
*/
@ -56,10 +55,11 @@ public:
// 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
//! Index of the pressure
static constexpr int pressure0Idx = PVOffset;
//! Index of the molefraction of the first component
static constexpr int z0Idx = pressure0Idx + 1;
// equation indices

106
opm/models/ptflash/flashintensivequantities.hh
View File

@ -25,8 +25,8 @@
*
* \copydoc Opm::FlashIntensiveQuantities
*/
#ifndef EWOMS_FLASH_INTENSIVE_QUANTITIES_HH
#define EWOMS_FLASH_INTENSIVE_QUANTITIES_HH
#ifndef OPM_FLASH_INTENSIVE_QUANTITIES_HH
#define OPM_FLASH_INTENSIVE_QUANTITIES_HH
#include "flashproperties.hh"
#include "flashindices.hh"
@ -46,7 +46,7 @@ namespace Opm {
* \ingroup FlashModel
* \ingroup IntensiveQuantities
*
* \brief Contains the intensive quantities of the flash-based compositional multi-phase model
* \brief Contains the intensive quantities of the ptflash-based compositional multi-phase model
*/
template <class TypeTag>
class FlashIntensiveQuantities
@ -78,7 +78,7 @@ class FlashIntensiveQuantities
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>;
@ -90,8 +90,7 @@ public:
//! The type of the object returned by the fluidState() method
using FluidState = Opm::CompositionalFluidState<Evaluation, FluidSystem, enableEnergy>;
FlashIntensiveQuantities()
{ }
FlashIntensiveQuantities() = default;
FlashIntensiveQuantities(const FlashIntensiveQuantities& other) = default;
@ -107,11 +106,11 @@ public:
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);
constexpr Scalar flashTolerance = getPropValue<TypeTag, Properties::FlashTolerance>();
constexpr int flashVerbosity = getPropValue<TypeTag, Properties::FlashVerbosity>();
const std::string flashTwoPhaseMethod = getPropValue<TypeTag, Properties::FlashTwoPhaseMethod>();
// extract the total molar densities of the components
ComponentVector z(0.);
{
@ -128,7 +127,6 @@ public:
sumz +=z[compIdx];
}
z /= sumz;
}
Evaluation p = priVars.makeEvaluation(pressure0Idx, timeIdx);
@ -164,45 +162,40 @@ public:
}
/////////////
// Compute the phase compositions and densities
// 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);
if (flashVerbosity >= 1) {
const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
std::cout << " updating the intensive quantities for Cell " << spatialIdx << std::endl;
}
FlashSolver::solve(fluidState_, z, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
using Flash = Opm::PTFlash<double, FluidSystem>;
FlashSolver::solve(fluidState_, z, spatialIdx, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
if (flashVerbosity >= 5) {
// printing of flash result after solve
const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
std::cout << " \n After flash solve for cell " << spatialIdx << std::endl;
ComponentVector x, y;
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:" << std::endl;
std::cout << x[comp_idx] << std::endl;
//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:" << std::endl;
std::cout << y[comp_idx] << std::endl;
}
const Evaluation& L = fluidState_.L();
std::cout << " L is:" << std::endl;
std::cout << L << 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
// 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));
@ -213,12 +206,12 @@ public:
// 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);
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);
@ -226,37 +219,36 @@ public:
fluidState_.setCompressFactor(1, Z_V);
// Print saturation
if (flashVerbosity == 5) {
if (flashVerbosity >= 5) {
std::cout << "So = " << So <<std::endl;
std::cout << "Sg = " << Sg <<std::endl;
}
// Print saturation
if (flashVerbosity == 5) {
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
// Compute rel. perm and viscosity and densities
/////////////
const MaterialLawParams& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
// calculate relative permeabilities
// calculate relative permeability
MaterialLaw::relativePermeabilities(relativePermeability_,
materialParams, fluidState_);
Opm::Valgrind::CheckDefined(relativePermeability_);
// set the phase viscosities and density
// set the phase viscosity 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]);
@ -322,8 +314,8 @@ private:
DimMatrix intrinsicPerm_;
FluidState fluidState_;
Evaluation porosity_;
Evaluation relativePermeability_[numPhases];
Evaluation mobility_[numPhases];
std::array<Evaluation,numPhases> relativePermeability_;
std::array<Evaluation,numPhases> mobility_;
};
} // namespace Opm

12
opm/models/ptflash/flashlocalresidual.hh
View File

@ -25,10 +25,8 @@
*
* \copydoc Opm::FlashLocalResidual
*/
#ifndef EWOMS_FLASH_LOCAL_RESIDUAL_HH
#define EWOMS_FLASH_LOCAL_RESIDUAL_HH
#include "flashproperties.hh"
#ifndef OPM_PTFLASH_LOCAL_RESIDUAL_HH
#define OPM_PTFLASH_LOCAL_RESIDUAL_HH
#include <opm/models/common/diffusionmodule.hh>
#include <opm/models/common/energymodule.hh>
@ -40,7 +38,7 @@ namespace Opm {
* \ingroup FlashModel
*
* \brief Calculates the local residual of the compositional multi-phase
* model based on flash calculations.
* model based on PTFlash calculations.
*/
template <class TypeTag>
class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResidual>
@ -88,7 +86,7 @@ public:
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.porosity());
}
EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
}
@ -137,7 +135,7 @@ public:
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));
auto 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

54
opm/models/ptflash/flashmodel.hh
View File

@ -25,18 +25,18 @@
*
* \copydoc Opm::FlashModel
*/
#ifndef EWOMS_FLASH_MODEL_HH
#define EWOMS_FLASH_MODEL_HH
#ifndef OPM_PTFLASH_MODEL_HH
#define OPM_PTFLASH_MODEL_HH
#include <opm/material/densead/Math.hpp>
#include "flashproperties.hh"
#include "flashprimaryvariables.hh"
#include "flashlocalresidual.hh"
#include "flashratevector.hh"
#include "flashboundaryratevector.hh"
#include <opm/models/flash/flashratevector.hh>
#include <opm/models/flash/flashboundaryratevector.hh>
#include "flashintensivequantities.hh"
#include "flashextensivequantities.hh"
#include <opm/models/flash/flashextensivequantities.hh>
#include "flashindices.hh"
#include "flashnewtonmethod.hh"
@ -45,6 +45,7 @@
#include <opm/models/io/vtkcompositionmodule.hh>
#include <opm/models/io/vtkenergymodule.hh>
#include <opm/models/io/vtkdiffusionmodule.hh>
#include <opm/models/io/vtkptflashmodule.hh>
#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
#include <opm/material/constraintsolvers/PTFlash.hpp>
@ -64,6 +65,7 @@ namespace TTag {
struct FlashModel { using InheritsFrom = std::tuple<VtkDiffusion,
VtkEnergy,
VtkComposition,
VtkPTFlash,
MultiPhaseBaseModel>; };
} // namespace TTag
@ -86,7 +88,7 @@ template<class TypeTag>
struct FlashTolerance<TypeTag, TTag::FlashModel>
{
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = -1.0;
static constexpr type value = 1.e-12;
};
// Flash solver verbosity
@ -184,28 +186,11 @@ struct FluxModule<TypeTag, TTag::MyProblemTypeTag> { using type = Opm::Forchheim
* 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
* calculates the compositions, saturation 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.)
* some model assumptions internally. Here a constant pressure, constant temperature,
* two-phase flash calculation method is used.
*
* 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
@ -227,7 +212,7 @@ class FlashModel
using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
public:
FlashModel(Simulator& simulator)
explicit FlashModel(Simulator& simulator)
: ParentType(simulator)
{}
@ -240,6 +225,7 @@ public:
// register runtime parameters of the VTK output modules
Opm::VtkCompositionModule<TypeTag>::registerParameters();
Opm::VtkPTFlashModule<TypeTag>::registerParameters();
if (enableDiffusion)
Opm::VtkDiffusionModule<TypeTag>::registerParameters();
@ -253,22 +239,17 @@ public:
EWOMS_REGISTER_PARAM(TypeTag, int, FlashVerbosity,
"Flash solver verbosity level");
EWOMS_REGISTER_PARAM(TypeTag, std::string, FlashTwoPhaseMethod,
"Method for solving vapor-liquid composition");
"Method for solving vapor-liquid composition. Available options include:"
"ssi, newton, ssi+newton");
}
/*!
* \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 != "")
if (!tmp.empty())
return tmp;
std::ostringstream oss;
@ -288,7 +269,7 @@ public:
std::string eqName(unsigned eqIdx) const
{
const std::string& tmp = EnergyModule::eqName(eqIdx);
if (tmp != "")
if (!tmp.empty())
return tmp;
std::ostringstream oss;
@ -343,6 +324,7 @@ public:
// add the VTK output modules which are meaningful for the model
this->addOutputModule(new Opm::VtkCompositionModule<TypeTag>(this->simulator_));
this->addOutputModule(new Opm::VtkPTFlashModule<TypeTag>(this->simulator_));
if (enableDiffusion)
this->addOutputModule(new Opm::VtkDiffusionModule<TypeTag>(this->simulator_));
if (enableEnergy)

42
opm/models/ptflash/flashnewtonmethod.hh
View File

@ -25,10 +25,8 @@
*
* \copydoc Opm::FlashNewtonMethod
*/
#ifndef EWOMS_FLASH_NEWTON_METHOD_HH
#define EWOMS_FLASH_NEWTON_METHOD_HH
#include "flashproperties.hh"
#ifndef OPM_FLASH_NEWTON_METHOD_HH
#define OPM_FLASH_NEWTON_METHOD_HH
#include <opm/models/nonlinear/newtonmethod.hh>
@ -47,7 +45,7 @@ namespace Opm {
/*!
* \ingroup FlashModel
*
* \brief A Newton solver specific to the NCP model.
* \brief A Newton solver specific to the PTFlash model.
*/
template <class TypeTag>
@ -79,11 +77,11 @@ protected:
/*!
* \copydoc FvBaseNewtonMethod::updatePrimaryVariables_
*/
void updatePrimaryVariables_(unsigned globalDofIdx,
void updatePrimaryVariables_(unsigned /* globalDofIdx */,
PrimaryVariables& nextValue,
const PrimaryVariables& currentValue,
const EqVector& update,
const EqVector& currentResidual)
const EqVector& /* currentResidual */)
{
// normal Newton-Raphson update
nextValue = currentValue;
@ -107,13 +105,15 @@ protected:
maxDeltaZ = std::max(std::abs(update[z0Idx + compIdx]), maxDeltaZ);
sumDeltaZ += update[z0Idx + compIdx];
}
maxDeltaZ = std::max(std::abs(-sumDeltaZ), maxDeltaZ);
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;
// if max. update is above 0.2, restrict that one to 0.2 and adjust the rest accordingly (s.t. last comp. update is sum of the changes in update vector)
// \Note: original code uses 0.1, while 0.1 looks like having problem make it converged. So there is some more to investigate here
constexpr Scalar deltaz_limit = 0.2;
if (maxDeltaZ > deltaz_limit) {
Scalar alpha = deltaz_limit / maxDeltaZ;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
nextValue[z0Idx + compIdx] = currentValue[z0Idx + compIdx] - alpha*update[z0Idx + compIdx];
nextValue[z0Idx + compIdx] = currentValue[z0Idx + compIdx] - alpha * update[z0Idx + compIdx];
}
}
@ -122,23 +122,13 @@ protected:
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)); }
{
val = std::clamp(val, minVal, maxVal);
}
}; // class FlashNewtonMethod
} // namespace Opm
#endif
#endif

8
opm/models/ptflash/flashprimaryvariables.hh
View File

@ -25,11 +25,10 @@
*
* \copydoc Opm::FlashPrimaryVariables
*/
#ifndef EWOMS_FLASH_PRIMARY_VARIABLES_HH
#define EWOMS_FLASH_PRIMARY_VARIABLES_HH
#ifndef OPM_PTFLASH_PRIMARY_VARIABLES_HH
#define OPM_PTFLASH_PRIMARY_VARIABLES_HH
#include "flashindices.hh"
#include "flashproperties.hh"
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
#include <opm/models/common/energymodule.hh>
@ -121,7 +120,7 @@ public:
// the energy module
EnergyModule::setPriVarTemperatures(*this, fluidState);
// determine the phase presence.
// determine the component fractions
Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
@ -135,6 +134,7 @@ public:
for (int i = 0; i < numComponents - 1; ++i)
(*this)[z0Idx + i] = z[i];
(*this)[pressure0Idx] = Opm::getValue(fluidState.pressure(0));
}

4
opm/models/ptflash/flashproperties.hh
View File

@ -27,8 +27,8 @@
* \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
#ifndef OPM_PTFLASH_PROPERTIES_HH
#define OPM_PTFLASH_PROPERTIES_HH
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/io/vtkcompositionmodule.hh>

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

@ -1,143 +0,0 @@
// -*- 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,8 +37,6 @@
#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>
@ -263,9 +261,6 @@ 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,
@ -283,7 +278,6 @@ 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,9 +156,6 @@ 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;
}