From 31bfb2227375ad13fff3101e69f7dc8d3930f382 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 17 Nov 2017 14:49:01 +0100 Subject: [PATCH] Adapt to the refactoring of the Equil initialization code Most noteably the blackoilState is not returned by the initialization code anymore, instead the initialization class is returned. --- opm/autodiff/FlowMain.hpp | 45 ++++++++++++++++++++++++++++++-- opm/autodiff/FlowMainEbos.hpp | 49 ++++++++++++++++++++++++++++++++++- 2 files changed, 91 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 651007713..7d78c30f4 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -637,8 +637,23 @@ namespace Opm Opm::UgGridHelpers::numFaces(grid), props.numPhases())); - initStateEquil(grid, props, *deck_, *eclipse_state_, gravity_[2], *state_); - //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); + + typedef Opm::FluidSystems::BlackOil FluidSystem; + FluidSystem::initFromDeck(*deck_ , *eclipse_state_); + typedef EQUIL::DeckDependent::InitialStateComputer ISC; + + ISC isc(material_law_manager_, *eclipse_state_, grid, gravity_[2]); + + const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const int oilpos = FluidSystem::oilPhaseIdx; + const int waterpos = FluidSystem::waterPhaseIdx; + const int ref_phase = oil ? oilpos : waterpos; + + state_->pressure() = isc.press()[ref_phase]; + convertSats(state_->saturation(), isc.saturation(), pu); + state_->gasoilratio() = isc.rs(); + state_->rv() = isc.rv(); + } else { state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::numFaces(grid), @@ -674,6 +689,31 @@ namespace Opm } + template + void convertSats(std::vector& sat_interleaved, const std::vector< std::vector >& sat, const PhaseUsage& pu) + { + assert(sat.size() == 3); + const auto nc = sat[0].size(); + const auto np = sat_interleaved.size() / nc; + for (size_t c = 0; c < nc; ++c) { + if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sat_p = sat[ FluidSystem::oilPhaseIdx]; + sat_interleaved[np*c + opos] = sat_p[c]; + } + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const std::vector& sat_p = sat[ FluidSystem::waterPhaseIdx]; + sat_interleaved[np*c + wpos] = sat_p[c]; + } + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const std::vector& sat_p = sat[ FluidSystem::gasPhaseIdx]; + sat_interleaved[np*c + gpos] = sat_p[c]; + } + } + } + @@ -925,6 +965,7 @@ namespace Opm Base::threshold_pressures_, Base::defunct_well_names_)); } + }; diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 8b4141a78..6eccaa91c 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -42,6 +42,10 @@ #include +#include +#include +#include + #include #include #include @@ -533,7 +537,21 @@ namespace Opm Opm::UgGridHelpers::numFaces(grid), props.numPhases())); - initStateEquil(grid, props, deck(), eclState(), gravity(), *state_); + typedef EQUIL::DeckDependent::InitialStateComputer ISC; + + ISC isc(materialLawManager(), eclState(), grid, gravity()); + + const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const int oilpos = FluidSystem::oilPhaseIdx; + const int waterpos = FluidSystem::waterPhaseIdx; + const int ref_phase = oil ? oilpos : waterpos; + + state_->pressure() = isc.press()[ref_phase]; + convertSats(state_->saturation(), isc.saturation(), pu); + state_->gasoilratio() = isc.rs(); + state_->rv() = isc.rv(); + + //initStateEquil(grid, materialLawManager(), eclState(), gravity(), pu, *state_); //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); } else { state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), @@ -946,6 +964,35 @@ namespace Opm } } + /// Convert saturations from a vector of individual phase saturation vectors + /// to an interleaved format where all values for a given cell come before all + /// values for the next cell, all in a single vector. + template + void convertSats(std::vector& sat_interleaved, const std::vector< std::vector >& sat, const PhaseUsage& pu) + { + assert(sat.size() == 3); + const auto nc = sat[0].size(); + const auto np = sat_interleaved.size() / nc; + for (size_t c = 0; c < nc; ++c) { + if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sat_p = sat[ FluidSystem::oilPhaseIdx]; + sat_interleaved[np*c + opos] = sat_p[c]; + } + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const std::vector& sat_p = sat[ FluidSystem::waterPhaseIdx]; + sat_interleaved[np*c + wpos] = sat_p[c]; + } + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const std::vector& sat_p = sat[ FluidSystem::gasPhaseIdx]; + sat_interleaved[np*c + gpos] = sat_p[c]; + } + } + } + + std::unique_ptr ebosSimulator_; int mpi_rank_ = 0; bool output_cout_ = false;