OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

made additional test with co2 and brine working, similar as juliacode

Browse Source
This commit is contained in:
Trine Mykkeltvedt 2022-06-20 09:04:17 +02:00
parent 4fd4fc7029
commit 5b2845bf9a
5 changed files with 162 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
opm/material/constraintsolvers/ChiFlash.hpp
View File

@ -61,24 +61,17 @@ class ChiFlash
//using Problem = GetPropType<TypeTag, Properties::Problem>;
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
// enum { Comp0Idx = FluidSystem::Comp0Idx }; //rename for generic ?
// enum { Comp1Idx = FluidSystem::Comp1Idx }; //rename for generic ?
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx};
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx};
enum { numMiscibleComponents = 3}; //octane, co2 // should be brine instead of brine here.
enum { numMisciblePhases = 2}; //oil, gas
enum { numMiscibleComponents = FluidSystem::numMiscibleComponents}; //octane, co2 // should be brine instead of brine here.
enum { numMisciblePhases = FluidSystem::numMisciblePhases}; //oil, gas
enum {
numEq =
numMisciblePhases+
numMisciblePhases*numMiscibleComponents
};//pressure, saturation, composition
/* enum {
// p0PvIdx = 0, // pressure first phase primary variable index
// S0PvIdx = 1, // saturation first phase primary variable index
// x00PvIdx = S0PvIdx + 1, // molefraction first phase first component primary variable index
//numMiscibleComponennets*numMisciblePhases-1 molefractions/primvar follow
}; */
public:
/*!
@ -203,16 +196,6 @@ public:
updateDerivatives_(fluid_state_scalar, z, fluid_state, single);
//print summary after flash
if (verbosity >= 1) {
std::cout << " ------ SUMMARY AFTER FLASH ------ " << std::endl;
std::cout << " L " << fluid_state.L() << std::endl;
std::cout << " K " << fluid_state.K(0) << ", " << fluid_state.K(1) << ", " << fluid_state.K(2) << std::endl;
std::cout << " x " << fluid_state.moleFraction(oilPhaseIdx, 0) << ", " << fluid_state.moleFraction(oilPhaseIdx, 1) << ", " << fluid_state.moleFraction(oilPhaseIdx, 2) << std::endl;
std::cout << " y " << fluid_state.moleFraction(gasPhaseIdx, 0) << ", " << fluid_state.moleFraction(gasPhaseIdx, 1) << ", " << fluid_state.moleFraction(gasPhaseIdx, 2) << std::endl;
}
}//end solve
/*!

20
opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh
View File

@ -20,19 +20,21 @@ namespace Opm {
public:
// TODO: I do not think these should be constant in fluidsystem, will try to make it non-constant later
static constexpr int numPhases=2;
static constexpr int numComponents = 3;
static constexpr int numComponents = 2;
static constexpr int numMisciblePhases=2;
static constexpr int numMiscibleComponents = 2;
// TODO: phase location should be more general
static constexpr int oilPhaseIdx = 0;
static constexpr int gasPhaseIdx = 1;
static constexpr int Comp0Idx = 0;
static constexpr int Comp1Idx = 1;
static constexpr int Comp2Idx = 2;
//static constexpr int Comp2Idx = 2;
// TODO: needs to be more general
using Comp0 = Opm::JuliaCO2<Scalar>;
using Comp1 = Opm::ChiwomsBrine<Scalar>;
using Comp2 = Opm::JuliaC10<Scalar>;
//using Comp2 = Opm::JuliaC10<Scalar>;
template <class ValueType>
using ParameterCache = Opm::ChiParameterCache<ValueType, Co2BrineFluidSystem<Scalar>>;
@ -49,7 +51,7 @@ namespace Opm {
switch (compIdx) {
case Comp0Idx: return Comp0::acentricFactor();
case Comp1Idx: return Comp1::acentricFactor();
case Comp2Idx: return Comp2::acentricFactor();
// case Comp2Idx: return Comp2::acentricFactor();
default: throw std::runtime_error("Illegal component index for acentricFactor");
}
}
@ -63,7 +65,7 @@ namespace Opm {
switch (compIdx) {
case Comp0Idx: return Comp0::criticalTemperature();
case Comp1Idx: return Comp1::criticalTemperature();
case Comp2Idx: return Comp2::criticalTemperature();
// case Comp2Idx: return Comp2::criticalTemperature();
default: throw std::runtime_error("Illegal component index for criticalTemperature");
}
}
@ -76,7 +78,7 @@ namespace Opm {
switch (compIdx) {
case Comp0Idx: return Comp0::criticalPressure();
case Comp1Idx: return Comp1::criticalPressure();
case Comp2Idx: return Comp2::criticalPressure();
// case Comp2Idx: return Comp2::criticalPressure();
default: throw std::runtime_error("Illegal component index for criticalPressure");
}
}
@ -90,7 +92,7 @@ namespace Opm {
switch (compIdx) {
case Comp0Idx: return Comp0::criticalVolume();
case Comp1Idx: return Comp1::criticalVolume();
case Comp2Idx: return Comp2::criticalVolume();
// case Comp2Idx: return Comp2::criticalVolume();
default: throw std::runtime_error("Illegal component index for criticalVolume");
}
}
@ -101,7 +103,7 @@ namespace Opm {
switch (compIdx) {
case Comp0Idx: return Comp0::molarMass();
case Comp1Idx: return Comp1::molarMass();
case Comp2Idx: return Comp2::molarMass();
// case Comp2Idx: return Comp2::molarMass();
default: throw std::runtime_error("Illegal component index for molarMass");
}
}
@ -131,7 +133,7 @@ namespace Opm {
static const char* name[] = {
Comp0::name(),
Comp1::name(),
Comp2::name(),
// Comp2::name(),
};
assert(0 <= compIdx && compIdx < 3);

3
opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh
View File

@ -23,6 +23,8 @@ namespace Opm {
// TODO: I do not think these should be constant in fluidsystem, will try to make it non-constant later
static constexpr int numPhases=2;
static constexpr int numComponents = 3;
static constexpr int numMisciblePhases=2;
static constexpr int numMiscibleComponents = 3;
// TODO: phase location should be more general
static constexpr int oilPhaseIdx = 0;
static constexpr int gasPhaseIdx = 1;
@ -133,6 +135,7 @@ namespace Opm {
static const char* name[] = {
Comp0::name(),
Comp1::name(),
Comp2::name(),
};
assert(0 <= compIdx && compIdx < 3);

1
tests/test_chiflash.cpp
View File

@ -30,6 +30,7 @@
#include <opm/material/constraintsolvers/ChiFlash.hpp>
#include <opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>

144
tests/test_co2brine_flash.cpp 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
*
* \brief This is test for the ChiFlash flash solver.
*/
#include "config.h"
#include <opm/material/constraintsolvers/ChiFlash.hpp>
#include <opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <dune/common/parallel/mpihelper.hh>
// the following include should be removed later
// #include <opm/material/fluidsystems/chifluid/chiwoms.h>
void testCo2BrineFlash()
{
using Scalar = double;
using FluidSystem = Opm::Co2BrineFluidSystem<Scalar>;
constexpr auto numComponents = FluidSystem::numComponents;
using Evaluation = Opm::DenseAd::Evaluation<double, numComponents>;
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
typedef Opm::CompositionalFluidState<Evaluation, FluidSystem> FluidState;
// input
Evaluation p_init = Evaluation::createVariable(10e5, 0); // 10 bar
ComponentVector comp;
comp[0] = Evaluation::createVariable(0.5, 1);
comp[1] = 1. - comp[0];//Evaluation::createVariable(0.1, 1);
//comp[2] = 0;//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?
Scalar temp = 300.0;
// From co2-compositional branch, it uses
// typedef typename FluidSystem::template ParameterCache<Scalar> ParameterCache;
FluidState fs;
// TODO: no capillary pressure for now
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));
}
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 = "ssi"; // "ssi"
//const std::string flash_twophase_method = "newton";
const std::string flash_twophase_method = "newton";
// TODO: should we set these?
// 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.;
fs.setLvalue(Ltmp);
const int spatialIdx = 0;
using Flash = Opm::ChiFlash<double, FluidSystem>;
// TODO: here the zInit does not have the proper derivatives
Flash::solve(fs, zInit, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance);
}
int main(int argc, char **argv)
{
Dune::MPIHelper::instance(argc, argv);
testCo2BrineFlash();
return 0;
}