From 5b2845bf9a0e4f662a751a555e331c8933d62a82 Mon Sep 17 00:00:00 2001 From: Trine Mykkeltvedt Date: Mon, 20 Jun 2022 09:04:17 +0200 Subject: [PATCH] made additional test with co2 and brine working, similar as juliacode --- opm/material/constraintsolvers/ChiFlash.hpp | 23 +-- .../chifluid/co2brinefluidsystem.hh | 20 +-- .../juliathreecomponentfluidsystem.hh | 3 + tests/test_chiflash.cpp | 1 + tests/test_co2brine_flash.cpp | 144 ++++++++++++++++++ 5 files changed, 162 insertions(+), 29 deletions(-) create mode 100644 tests/test_co2brine_flash.cpp diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp index 6e18e4bd0..e8ba707c5 100644 --- a/opm/material/constraintsolvers/ChiFlash.hpp +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -61,24 +61,17 @@ class ChiFlash //using Problem = GetPropType; 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 /*! diff --git a/opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh b/opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh index 3a5bf06f2..e06d8feb7 100644 --- a/opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh +++ b/opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh @@ -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; using Comp1 = Opm::ChiwomsBrine; - using Comp2 = Opm::JuliaC10; + //using Comp2 = Opm::JuliaC10; template using ParameterCache = Opm::ChiParameterCache>; @@ -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); diff --git a/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh b/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh index 80a5d0840..a08101d9b 100644 --- a/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh +++ b/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh @@ -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); diff --git a/tests/test_chiflash.cpp b/tests/test_chiflash.cpp index d8c79fd66..68cce489a 100644 --- a/tests/test_chiflash.cpp +++ b/tests/test_chiflash.cpp @@ -30,6 +30,7 @@ #include #include + #include #include #include diff --git a/tests/test_co2brine_flash.cpp b/tests/test_co2brine_flash.cpp new file mode 100644 index 000000000..63b12776a --- /dev/null +++ b/tests/test_co2brine_flash.cpp @@ -0,0 +1,144 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief This is test for the ChiFlash flash solver. + */ +#include "config.h" + +#include +#include +#include +#include +#include +#include + +#include +// the following include should be removed later +// #include + +void testCo2BrineFlash() +{ + using Scalar = double; + using FluidSystem = Opm::Co2BrineFluidSystem; + + constexpr auto numComponents = FluidSystem::numComponents; + using Evaluation = Opm::DenseAd::Evaluation; + typedef Dune::FieldVector ComponentVector; + typedef Opm::CompositionalFluidState 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 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 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; + // 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; +}