diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index fe995939d..45b8d3414 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -26,6 +26,7 @@ # originally generated with the command: # find opm -name '*.c*' -printf '\t%p\n' | sort list (APPEND MAIN_SOURCE_FILES + opm/autodiff/Compat.cpp opm/autodiff/ExtractParallelGridInformationToISTL.cpp opm/autodiff/NewtonIterationBlackoilCPR.cpp opm/autodiff/NewtonIterationBlackoilInterleaved.cpp diff --git a/opm/autodiff/Compat.cpp b/opm/autodiff/Compat.cpp new file mode 100644 index 000000000..2b9d8c00b --- /dev/null +++ b/opm/autodiff/Compat.cpp @@ -0,0 +1,257 @@ +/* + Copyright 2016 Statoil ASA + Copyright 2016 IRIS + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + + 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 3 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 . +*/ + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + +std::vector< double > destripe( const std::vector< double >& v, + size_t stride, + size_t offset ) { + + std::vector< double > dst( v.size() / stride ); + + size_t di = 0; + for( size_t i = offset; i < v.size(); i += stride ) { + dst[ di++ ] = v[ i ]; + } + + return dst; +} + + + + + + +std::vector< double >& stripe( const std::vector< double >& v, + size_t stride, + size_t offset, + std::vector< double >& dst ) { + + /* does little range checking etc; for future revisions */ + size_t vi = 0; + for( size_t i = offset; i < dst.size(); i += stride ) { + dst[ i ] = v[ vi++ ]; + } + + return dst; +} + + + + + + +data::Solution simToSolution( const SimulationDataContainer& reservoir, + PhaseUsage phases ) { + data::Solution sol; + sol.insert( "PRESSURE", UnitSystem::measure::pressure, reservoir.pressure() , data::TargetType::RESTART_SOLUTION); + sol.insert( "TEMP" , UnitSystem::measure::temperature, reservoir.temperature() , data::TargetType::RESTART_SOLUTION ); + + const auto ph = reservoir.numPhases(); + const auto& sat = reservoir.saturation(); + + const auto aqua = BlackoilPhases::Aqua; + const auto vapour = BlackoilPhases::Vapour; + + if( phases.phase_used[ aqua ] ) { + sol.insert( "SWAT", UnitSystem::measure::identity, destripe( sat, ph, phases.phase_pos[ aqua ] ) , data::TargetType::RESTART_SOLUTION ); + } + + if( phases.phase_used[ vapour ] ) { + sol.insert( "SGAS", UnitSystem::measure::identity, destripe( sat, ph, phases.phase_pos[ vapour ] ) , data::TargetType::RESTART_SOLUTION ); + } + + if( reservoir.hasCellData( BlackoilState::GASOILRATIO ) ) { + sol.insert( "RS", UnitSystem::measure::gas_oil_ratio, reservoir.getCellData( BlackoilState::GASOILRATIO ) , data::TargetType::RESTART_SOLUTION ); + } + + if( reservoir.hasCellData( BlackoilState::RV ) ) { + sol.insert( "RV", UnitSystem::measure::oil_gas_ratio, reservoir.getCellData( BlackoilState::RV ) , data::TargetType::RESTART_SOLUTION ); + } + + if (reservoir.hasCellData( BlackoilSolventState::SSOL)) { + sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilSolventState::SSOL ) , data::TargetType::RESTART_SOLUTION ); + } + + return sol; +} + + + + + + +void solutionToSim( const data::Solution& sol, + PhaseUsage phases, + SimulationDataContainer& state ) { + + const auto stride = phases.num_phases; + + if( sol.has( "SWAT" ) ) { + stripe( sol.data( "SWAT" ), + stride, + phases.phase_pos[ BlackoilPhases::Aqua ], + state.saturation() ); + } + + if( sol.has( "SGAS" ) ) { + stripe( sol.data( "SGAS" ), + stride, + phases.phase_pos[ BlackoilPhases::Vapour ], + state.saturation() ); + } + + for (size_t c = 0; c < state.numCells(); ++c) { + double& so = state.saturation()[phases.num_phases*c + phases.phase_pos[ BlackoilPhases::Liquid ]]; + so = 1.0; + if (phases.phase_used[ BlackoilPhases::Aqua]) { + so -= state.saturation()[phases.num_phases*c + phases.phase_pos[ BlackoilPhases::Aqua ]]; + } + if (phases.phase_used[ BlackoilPhases::Vapour]) { + so -= state.saturation()[phases.num_phases*c + phases.phase_pos[ BlackoilPhases::Vapour ]]; + } + } + + + if( sol.has( "PRESSURE" ) ) { + state.pressure() = sol.data( "PRESSURE" ); + } + + if( sol.has( "TEMP" ) ) { + state.temperature() = sol.data( "TEMP" ); + } + + if( sol.has( "RS" ) ) { + state.getCellData( "GASOILRATIO" ) = sol.data( "RS" ); + } + + if( sol.has( "RV" ) ) { + state.getCellData( "RV" ) = sol.data( "RV" ); + } + + if (sol.has( "SSOL" ) ) { + state.getCellData("SSOL") = sol.data("SSOL"); + } + +} + + + + + + +void wellsToState( const data::Wells& wells, + PhaseUsage phases, + WellStateFullyImplicitBlackoil& state ) { + + using rt = data::Rates::opt; + + const auto np = phases.num_phases; + + std::vector< rt > phs( np ); + if( phases.phase_used[BlackoilPhases::Aqua] ) { + phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat; + } + + if( phases.phase_used[BlackoilPhases::Liquid] ) { + phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil; + } + + if( phases.phase_used[BlackoilPhases::Vapour] ) { + phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; + } + + for( const auto& wm : state.wellMap() ) { + const auto well_index = wm.second[ 0 ]; + const auto& well = wells.at( wm.first ); + + state.bhp()[ well_index ] = well.bhp; + state.temperature()[ well_index ] = well.temperature; + state.currentControls()[ well_index ] = well.control; + + const auto wellrate_index = well_index * np; + for( size_t i = 0; i < phs.size(); ++i ) { + assert( well.rates.has( phs[ i ] ) ); + state.wellRates()[ wellrate_index + i ] = well.rates.get( phs[ i ] ); + } + + const auto perforation_pressure = []( const data::Completion& comp ) { + return comp.pressure; + }; + + const auto perforation_reservoir_rate = []( const data::Completion& comp ) { + return comp.reservoir_rate; + }; + + std::transform( well.completions.begin(), + well.completions.end(), + state.perfPress().begin() + wm.second[ 1 ], + perforation_pressure ); + + std::transform( well.completions.begin(), + well.completions.end(), + state.perfRates().begin() + wm.second[ 1 ], + perforation_reservoir_rate ); + + int local_comp_index = 0; + for (const data::Completion& comp : well.completions) { + const int global_comp_index = wm.second[1] + local_comp_index; + for (int phase_index = 0; phase_index < np; ++phase_index) { + state.perfPhaseRates()[global_comp_index*np + phase_index] = comp.rates.get(phs[phase_index]); + } + ++local_comp_index; + } + } +} + + + + + + +void wellsToState( const data::Wells& wells, + PhaseUsage phases, + WellStateFullyImplicitBlackoilDense& state ) +{ + // Set base class variables. + wellsToState(wells, phases, static_cast(state)); + + // Set wellSolution() variable. + state.setWellSolutions(phases); +} + + +} // namespace Opm diff --git a/opm/autodiff/Compat.hpp b/opm/autodiff/Compat.hpp index 46c97c3a2..76ee2d033 100644 --- a/opm/autodiff/Compat.hpp +++ b/opm/autodiff/Compat.hpp @@ -21,229 +21,56 @@ #ifndef OPM_SIMULATORS_COMPAT_HPP #define OPM_SIMULATORS_COMPAT_HPP -#include -#include - -#include -#include -#include -#include -#include -#include #include #include +#include +#include namespace Opm { -inline std::vector< double > destripe( const std::vector< double >& v, - size_t stride, - size_t offset ) { - - std::vector< double > dst( v.size() / stride ); - - size_t di = 0; - for( size_t i = offset; i < v.size(); i += stride ) { - dst[ di++ ] = v[ i ]; - } - - return dst; -} - - - - - -inline std::vector< double >& stripe( const std::vector< double >& v, - size_t stride, - size_t offset, - std::vector< double >& dst ) { - - /* does little range checking etc; for future revisions */ - size_t vi = 0; - for( size_t i = offset; i < dst.size(); i += stride ) { - dst[ i ] = v[ vi++ ]; - } - - return dst; -} - - - - - - - - inline data::Solution simToSolution( const SimulationDataContainer& reservoir, - PhaseUsage phases ) { - data::Solution sol; - sol.insert( "PRESSURE", UnitSystem::measure::pressure, reservoir.pressure() , data::TargetType::RESTART_SOLUTION); - sol.insert( "TEMP" , UnitSystem::measure::temperature, reservoir.temperature() , data::TargetType::RESTART_SOLUTION ); - - const auto ph = reservoir.numPhases(); - const auto& sat = reservoir.saturation(); - - const auto aqua = BlackoilPhases::Aqua; - const auto vapour = BlackoilPhases::Vapour; - - if( phases.phase_used[ aqua ] ) { - sol.insert( "SWAT", UnitSystem::measure::identity, destripe( sat, ph, phases.phase_pos[ aqua ] ) , data::TargetType::RESTART_SOLUTION ); - } - - if( phases.phase_used[ vapour ] ) { - sol.insert( "SGAS", UnitSystem::measure::identity, destripe( sat, ph, phases.phase_pos[ vapour ] ) , data::TargetType::RESTART_SOLUTION ); - } - - if( reservoir.hasCellData( BlackoilState::GASOILRATIO ) ) { - sol.insert( "RS", UnitSystem::measure::gas_oil_ratio, reservoir.getCellData( BlackoilState::GASOILRATIO ) , data::TargetType::RESTART_SOLUTION ); - } - - if( reservoir.hasCellData( BlackoilState::RV ) ) { - sol.insert( "RV", UnitSystem::measure::oil_gas_ratio, reservoir.getCellData( BlackoilState::RV ) , data::TargetType::RESTART_SOLUTION ); - } - - if (reservoir.hasCellData( BlackoilSolventState::SSOL)) { - sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilSolventState::SSOL ) , data::TargetType::RESTART_SOLUTION ); - } - - return sol; -} - - - - - - - - -inline void solutionToSim( const data::Solution& sol, - PhaseUsage phases, - SimulationDataContainer& state ) { - - const auto stride = phases.num_phases; - - if( sol.has( "SWAT" ) ) { - stripe( sol.data( "SWAT" ), - stride, - phases.phase_pos[ BlackoilPhases::Aqua ], - state.saturation() ); - } - - if( sol.has( "SGAS" ) ) { - stripe( sol.data( "SGAS" ), - stride, - phases.phase_pos[ BlackoilPhases::Vapour ], - state.saturation() ); - } - - for (size_t c = 0; c < state.numCells(); ++c) { - double& so = state.saturation()[phases.num_phases*c + phases.phase_pos[ BlackoilPhases::Liquid ]]; - so = 1.0; - if (phases.phase_used[ BlackoilPhases::Aqua]) { - so -= state.saturation()[phases.num_phases*c + phases.phase_pos[ BlackoilPhases::Aqua ]]; - } - if (phases.phase_used[ BlackoilPhases::Vapour]) { - so -= state.saturation()[phases.num_phases*c + phases.phase_pos[ BlackoilPhases::Vapour ]]; - } - } - - - if( sol.has( "PRESSURE" ) ) { - state.pressure() = sol.data( "PRESSURE" ); - } - - if( sol.has( "TEMP" ) ) { - state.temperature() = sol.data( "TEMP" ); - } - - if( sol.has( "RS" ) ) { - state.getCellData( "GASOILRATIO" ) = sol.data( "RS" ); - } - - if( sol.has( "RV" ) ) { - state.getCellData( "RV" ) = sol.data( "RV" ); - } - - if (sol.has( "SSOL" ) ) { - state.getCellData("SSOL") = sol.data("SSOL"); - } - -} - - - - - - - - - - - - - -inline void wellsToState( const data::Wells& wells, - PhaseUsage phases, - WellStateFullyImplicitBlackoil& state ) { - - using rt = data::Rates::opt; - - const auto np = phases.num_phases; - - std::vector< rt > phs( np ); - if( phases.phase_used[BlackoilPhases::Aqua] ) { - phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat; - } - - if( phases.phase_used[BlackoilPhases::Liquid] ) { - phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil; - } - - if( phases.phase_used[BlackoilPhases::Vapour] ) { - phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; - } - - for( const auto& wm : state.wellMap() ) { - const auto well_index = wm.second[ 0 ]; - const auto& well = wells.at( wm.first ); - - state.bhp()[ well_index ] = well.bhp; - state.temperature()[ well_index ] = well.temperature; - state.currentControls()[ well_index ] = well.control; - - const auto wellrate_index = well_index * np; - for( size_t i = 0; i < phs.size(); ++i ) { - assert( well.rates.has( phs[ i ] ) ); - state.wellRates()[ wellrate_index + i ] = well.rates.get( phs[ i ] ); - } - - const auto perforation_pressure = []( const data::Completion& comp ) { - return comp.pressure; - }; - - const auto perforation_reservoir_rate = []( const data::Completion& comp ) { - return comp.reservoir_rate; - }; - - std::transform( well.completions.begin(), - well.completions.end(), - state.perfPress().begin() + wm.second[ 1 ], - perforation_pressure ); - - std::transform( well.completions.begin(), - well.completions.end(), - state.perfRates().begin() + wm.second[ 1 ], - perforation_reservoir_rate ); - - int local_comp_index = 0; - for (const data::Completion& comp : well.completions) { - const int global_comp_index = wm.second[1] + local_comp_index; - for (int phase_index = 0; phase_index < np; ++phase_index) { - state.perfPhaseRates()[global_comp_index*np + phase_index] = comp.rates.get(phs[phase_index]); - } - ++local_comp_index; - } - } -} + // Forward declarations + class SimulationDataContainer; + class WellStateFullyImplicitBlackoil; + class WellStateFullyImplicitBlackoilDense; + + /// Extract single data vector from striped data. + /// \return u such that u[i] = v[offset + i*stride]. + std::vector< double > destripe( const std::vector< double >& v, + size_t stride, + size_t offset ); + + /// Inject single data vector into striped data. + /// \return reference to dst input, that is changed so that + /// dst[offset + i*stride] = v[i]. This is done for + /// i = 0..(dst.size()/stride). + std::vector< double >& stripe( const std::vector< double >& v, + size_t stride, + size_t offset, + std::vector< double >& dst ); + + /// Returns Solution with the following fields: + /// PRESSURE, TEMP (unconditionally) + /// SWAT, SGAS, RS, RV, SSOL (if appropriate fields present in input) + data::Solution simToSolution( const SimulationDataContainer& reservoir, + PhaseUsage phases ); + + /// Copies the following fields from sol into state (all conditionally): + /// PRESSURE, TEMP, SWAT, SGAS, RS, RV, SSOL + void solutionToSim( const data::Solution& sol, + PhaseUsage phases, + SimulationDataContainer& state ); + + /// Copies the following fields from wells into state. + /// bhp, temperature, currentControls, wellRates, perfPress, perfRates, perfPhaseRates + void wellsToState( const data::Wells& wells, + PhaseUsage phases, + WellStateFullyImplicitBlackoil& state ); + + /// As the WellStateFullyImplicitBlackoil overload, but also sets + /// the wellSolution field from the values of the other fields. + void wellsToState( const data::Wells& wells, + PhaseUsage phases, + WellStateFullyImplicitBlackoilDense& state ); } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 23d1331e7..95a624c9c 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -395,13 +395,13 @@ namespace Opm } - template + template inline void BlackoilOutputWriter:: initFromRestartFile( const PhaseUsage& phaseUsage, const Grid& grid, SimulationDataContainer& simulatorstate, - WellStateFullyImplicitBlackOel& wellstate) + WellState& wellstate) { std::map solution_keys {{"PRESSURE" , UnitSystem::measure::pressure}, {"SWAT" , UnitSystem::measure::identity}, diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index 6f041f6b6..2f6710d61 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -35,6 +35,7 @@ #include #include #include +#include namespace Opm { @@ -67,52 +68,94 @@ namespace Opm // call init on base class BaseType :: init(wells, state, prevState); - // if there are no well, do nothing in init - if (wells == 0) { + setWellSolutions(pu); + setWellSolutionsFromPrevState(prevState); + } + + + template + void setWellSolutionsFromPrevState(const PrevState& prevState) + { + // Set nw and np, or return if no wells. + if (wells_.get() == nullptr) { return; } + const int nw = wells_->number_of_wells; + if (nw == 0) { + return; + } + const int np = wells_->number_of_phases; - const int nw = wells->number_of_wells; - if( nw == 0 ) return ; + // intialize wells that have been there before + // order may change so the mapping is based on the well name + // if there are no well, do nothing in init + if( ! prevState.wellMap().empty() ) + { + typedef typename WellMapType :: const_iterator const_iterator; + const_iterator end = prevState.wellMap().end(); + int nw_old = prevState.bhp().size(); + for (int w = 0; w < nw; ++w) { + std::string name( wells_->name[ w ] ); + const_iterator it = prevState.wellMap().find( name ); + if( it != end ) + { + const int oldIndex = (*it).second[ 0 ]; + const int newIndex = w; - const int np = wells->number_of_phases; + // wellSolutions + for( int i = 0; i < np; ++i) + { + wellSolutions()[ i*nw + newIndex ] = prevState.wellSolutions()[i * nw_old + oldIndex ]; + } + + } + } + } + } + + + /// Set wellSolutions() based on the base class members. + void setWellSolutions(const PhaseUsage& pu) + { + // Set nw and np, or return if no wells. + if (wells_.get() == nullptr) { + return; + } + const int nw = wells_->number_of_wells; + if (nw == 0) { + return; + } + const int np = wells_->number_of_phases; well_solutions_.clear(); well_solutions_.resize(nw * np, 0.0); std::vector g = {1.0,1.0,0.01}; for (int w = 0; w < nw; ++w) { - WellControls* wc = wells->ctrls[w]; + WellControls* wc = wells_->ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. const int current = currentControls()[w]; well_controls_set_current( wc, current); - const WellType& well_type = wells->type[w]; + const WellType& well_type = wells_->type[w]; switch (well_controls_iget_type(wc, current)) { case THP: // Intentional fall-through case BHP: - { if (well_type == INJECTOR) { for (int p = 0; p < np; ++p) { - well_solutions_[w] += wellRates()[np*w + p] * wells->comp_frac[np*w + p]; + well_solutions_[w] += wellRates()[np*w + p] * wells_->comp_frac[np*w + p]; } } else { for (int p = 0; p < np; ++p) { well_solutions_[w] += g[p] * wellRates()[np*w + p]; } } - } break; - - case RESERVOIR_RATE: // Intentional fall-through case SURFACE_RATE: - { wellSolutions()[w] = bhp()[w]; - - } break; } @@ -135,41 +178,16 @@ namespace Opm } } else { if( pu.phase_used[Water] ) { - wellSolutions()[nw + w] = wells->comp_frac[np*w + waterpos]; + wellSolutions()[nw + w] = wells_->comp_frac[np*w + waterpos]; } if( pu.phase_used[Gas] ) { - wellSolutions()[2*nw + w] = wells->comp_frac[np*w + gaspos]; - } - } - } - - - // intialize wells that have been there before - // order may change so the mapping is based on the well name - if( ! prevState.wellMap().empty() ) - { - typedef typename WellMapType :: const_iterator const_iterator; - const_iterator end = prevState.wellMap().end(); - int nw_old = prevState.bhp().size(); - for (int w = 0; w < nw; ++w) { - std::string name( wells->name[ w ] ); - const_iterator it = prevState.wellMap().find( name ); - if( it != end ) - { - const int oldIndex = (*it).second[ 0 ]; - const int newIndex = w; - - // wellSolutions - for( int i = 0; i < np; ++i) - { - wellSolutions()[ i*nw + newIndex ] = prevState.wellSolutions()[i * nw_old + oldIndex ]; - } - + wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos]; } } } } + template void resize(const Wells* wells, const State& state, const PhaseUsage& pu ) { const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes @@ -186,6 +204,7 @@ namespace Opm return res; } + private: std::vector well_solutions_; };