making test_chiflash_scalar to be CO2, C10, C1 test

This commit is contained in:
Kai Bao 2021-12-02 15:21:09 +01:00 committed by Trine Mykkeltvedt
parent 7c2694bba4
commit 96c3462571
3 changed files with 173 additions and 14 deletions

View File

@ -1,8 +1,6 @@
#ifndef COMPONENTS_HH
#define COMPONENTS_HH
#include "chiwoms.h"
#include <opm/material/IdealGas.hpp>
#include <opm/material/components/Component.hpp>
#include <opm/material/components/SimpleCO2.hpp>

View File

@ -2,6 +2,11 @@
#define OPM_JULIATHREECOMPONENTFLUIDSYSTEM_HH
#include <opm/material/fluidsystems/BaseFluidSystem.hpp>
#include <opm/material/fluidsystems/chifluid/components.hh>
// TODO: this is something else need to check
#include "ChiParameterCache.hpp"
#include "LBCviscosity.hpp"
namespace Opm {
/*!
@ -10,9 +15,163 @@ namespace Opm {
* \brief A two phase three component fluid system from the Julia test
*/
template <class Scalar>
class JuliaThreeComponentFluidSystem
: public Opm::BaseFluidSystem<Scalar, JuliaThreeComponentFluidSystem<Scalar> > {
template<class Scalar>
class JuliaThreeComponentFluidSystem
: public Opm::BaseFluidSystem<Scalar, JuliaThreeComponentFluidSystem<Scalar> > {
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;
// 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;
// TODO: needs to be more general
using Comp0 = Opm::JuliaCO2<Scalar>;
using Comp1 = Opm::JuliaC1<Scalar>;
using Comp2 = Opm::JuliaC10<Scalar>;
template <class ValueType>
using ParameterCache = Opm::ChiParameterCache<ValueType, JuliaThreeComponentFluidSystem<Scalar>>;
using LBCviscosity = typename Opm::LBCviscosity<Scalar, JuliaThreeComponentFluidSystem<Scalar>>;
using PengRobinsonMixture = typename Opm::PengRobinsonMixture<Scalar, JuliaThreeComponentFluidSystem<Scalar>>;
/*!
* \brief The acentric factor of a component [].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar acentricFactor(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::acentricFactor();
case Comp1Idx: return Comp1::acentricFactor();
case Comp2Idx: return Comp2::acentricFactor();
default: throw std::runtime_error("Illegal component index for acentricFactor");
}
}
/*!
* \brief Critical temperature of a component [K].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar criticalTemperature(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::criticalTemperature();
case Comp1Idx: return Comp1::criticalTemperature();
case Comp2Idx: return Comp2::criticalTemperature();
default: throw std::runtime_error("Illegal component index for criticalTemperature");
}
}
/*!
* \brief Critical pressure of a component [Pa].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar criticalPressure(unsigned compIdx) {
switch (compIdx) {
case Comp0Idx: return Comp0::criticalPressure();
case Comp1Idx: return Comp1::criticalPressure();
case Comp2Idx: return Comp2::criticalPressure();
default: throw std::runtime_error("Illegal component index for criticalPressure");
}
}
/*!
* \brief Critical volume of a component [m3].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar criticalVolume(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::criticalVolume();
case Comp1Idx: return Comp1::criticalVolume();
case Comp2Idx: return Comp2::criticalVolume();
default: throw std::runtime_error("Illegal component index for criticalVolume");
}
}
//! \copydoc BaseFluidSystem::molarMass
static Scalar molarMass(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::molarMass();
case Comp1Idx: return Comp1::molarMass();
case Comp2Idx: return Comp2::molarMass();
default: throw std::runtime_error("Illegal component index for molarMass");
}
}
/*!
* \brief Returns the interaction coefficient for two components.
*.
*/
static Scalar interactionCoefficient(unsigned /*comp1Idx*/, unsigned /*comp2Idx*/)
{
return 0.0;
}
/*!
* \copydoc BaseFluidSystem::density
*/
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval density(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx)
{
LhsEval dens;
if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) {
// paramCache.updatePhase(fluidState, phaseIdx);
dens = fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx);
}
return dens;
}
//! \copydoc BaseFluidSystem::viscosity
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval viscosity(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx)
{
// Use LBC method to calculate viscosity
LhsEval mu;
// if (phaseIdx == gasPhaseIdx) {
mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx);
// }
// else {
// const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
// const auto& p = Opm::decay<LhsEval>(fluidState.pressure(0));
// mu = Brine::liquidViscosity(T, p);
// }
return mu;
}
//! \copydoc BaseFluidSystem::fugacityCoefficient
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval fugacityCoefficient(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx,
unsigned compIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
Scalar phi = Opm::getValue(
PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
return phi;
throw std::invalid_argument("crap!");
}
};
}
#endif //OPM_JULIATHREECOMPONENTFLUIDSYSTEM_HH

View File

@ -28,24 +28,23 @@
#include "config.h"
#include <opm/material/constraintsolvers/ChiFlash.hpp>
#include <opm/material/fluidsystems/chifluid/twophasefluidsystem.hh>
#include <opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
// #include <opm/material/constraintsolvers/NcpFlash.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
// #include <opm/material/fluidsystems/Spe5FluidSystem.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
// #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
#include <dune/common/parallel/mpihelper.hh>
// the following include should be removed later
// #include <opm/material/fluidsystems/chifluid/chiwoms.h>
void testChiFlash()
{
using Scalar = double;
using FluidSystem = Opm::TwoPhaseTwoComponentFluidSystem<Scalar>;
using FluidSystem = Opm::JuliaThreeComponentFluidSystem<Scalar>;
constexpr auto numComponents = FluidSystem::numComponents;
//using Evaluation = Opm::DenseAd::Evaluation<double, numComponents>;
@ -53,13 +52,14 @@ void testChiFlash()
typedef Opm::CompositionalFluidState<Scalar, FluidSystem> FluidState;
// input
Scalar p_init = 100.*1.e5; // 100 bar
Scalar p_init = 10e5; // 10 bar
ComponentVector comp;
comp[0] = MFCOMP0;
comp[1] = MFCOMP1;
comp[0] = 0.5;
comp[1] = 0.3;
comp[2] = 0.2;
ComponentVector sat;
sat[0] = 1.0; sat[1] = 1.0-sat[0];
Scalar temp = 303.0;
Scalar temp = 300.0;
// From co2-compositional branch, it uses
// typedef typename FluidSystem::template ParameterCache<Scalar> ParameterCache;
@ -71,9 +71,11 @@ void testChiFlash()
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]);
fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);