diff --git a/CMakeLists.txt b/CMakeLists.txt index 88c2acc1c..d30161752 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,6 +81,7 @@ include (CMakeLists_files.cmake) macro (config_hook) opm_need_version_of ("dune-common") opm_need_version_of ("dune-istl") + opm_need_version_of ("ewoms") endmacro (config_hook) macro (prereqs_hook) @@ -148,10 +149,7 @@ if (HAVE_OPM_DATA) endif() # create a symbolic link from flow to flow_legacy +message("create symlink") ADD_CUSTOM_TARGET(flow ALL - COMMAND ${CMAKE_COMMAND} -E create_symlink "flow_legacy" "${CMAKE_BINARY_DIR}/bin/flow") -install( - FILES "${CMAKE_BINARY_DIR}/bin/flow" - DESTINATION "${CMAKE_INSTALL_PREFIX}/bin" - ) + COMMAND ${CMAKE_COMMAND} -E create_symlink "flow_legacy" "${CMAKE_BINARY_DIR}/bin/flow") diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 9d719ea0b..74fba857f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -36,6 +36,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/moduleVersion.cpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp + opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp opm/autodiff/SimulatorIncompTwophaseAd.cpp opm/autodiff/TransportSolverTwophaseAd.cpp opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -49,7 +50,6 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/WellMultiSegment.cpp opm/autodiff/MultisegmentWells.cpp opm/autodiff/BlackoilSolventState.cpp - opm/autodiff/ThreadHandle.hpp opm/autodiff/MissingFeatures.cpp opm/polymer/PolymerState.cpp opm/polymer/PolymerBlackoilState.cpp @@ -105,6 +105,7 @@ list (APPEND EXAMPLE_SOURCE_FILES examples/find_zero.cpp examples/flow_legacy.cpp examples/flow_sequential.cpp + examples/flow_ebos.cpp examples/flow_multisegment.cpp examples/flow_solvent.cpp examples/sim_2p_incomp.cpp @@ -124,6 +125,7 @@ list (APPEND PROGRAM_SOURCE_FILES examples/sim_2p_incomp.cpp examples/sim_2p_incomp_ad.cpp examples/sim_2p_comp_reorder.cpp + examples/flow_ebos.cpp examples/flow_legacy.cpp examples/flow_sequential.cpp examples/flow_solvent.cpp @@ -143,6 +145,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/AutoDiffMatrix.hpp opm/autodiff/AutoDiff.hpp opm/autodiff/BackupRestore.hpp + opm/autodiff/BlackoilDetails.hpp opm/autodiff/BlackoilModel.hpp opm/autodiff/BlackoilModelBase.hpp opm/autodiff/BlackoilModelBase_impl.hpp @@ -155,6 +158,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/Compat.hpp opm/autodiff/CPRPreconditioner.hpp opm/autodiff/createGlobalCellArray.hpp + opm/autodiff/DefaultBlackoilSolutionState.hpp opm/autodiff/BlackoilSequentialModel.hpp opm/autodiff/BlackoilSolventModel.hpp opm/autodiff/BlackoilSolventModel_impl.hpp @@ -166,6 +170,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/DuneMatrix.hpp opm/autodiff/ExtractParallelGridInformationToISTL.hpp opm/autodiff/FlowMain.hpp + opm/autodiff/FlowMainEbos.hpp opm/autodiff/FlowMainPolymer.hpp opm/autodiff/FlowMainSequential.hpp opm/autodiff/FlowMainSolvent.hpp @@ -173,6 +178,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/GridHelpers.hpp opm/autodiff/GridInit.hpp opm/autodiff/ImpesTPFAAD.hpp + opm/autodiff/ISTLSolver.hpp + opm/autodiff/IterationReport.hpp opm/autodiff/moduleVersion.hpp opm/autodiff/NewtonIterationBlackoilCPR.hpp opm/autodiff/NewtonIterationBlackoilInterface.hpp @@ -187,18 +194,22 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp opm/autodiff/RateConverter.hpp opm/autodiff/RedistributeDataHandles.hpp + opm/autodiff/SimFIBODetails.hpp opm/autodiff/SimulatorBase.hpp opm/autodiff/SimulatorBase_impl.hpp + opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp opm/autodiff/SimulatorFullyImplicitBlackoil.hpp opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp + opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp opm/autodiff/SimulatorIncompTwophaseAd.hpp opm/autodiff/SimulatorSequentialBlackoil.hpp opm/autodiff/TransportSolverTwophaseAd.hpp opm/autodiff/WellDensitySegmented.hpp opm/autodiff/WellStateFullyImplicitBlackoil.hpp + opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp opm/autodiff/VFPProperties.hpp @@ -212,9 +223,11 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellHelpers.hpp opm/autodiff/StandardWells.hpp opm/autodiff/StandardWells_impl.hpp + opm/autodiff/StandardWellsDense.hpp opm/autodiff/StandardWellsSolvent.hpp opm/autodiff/StandardWellsSolvent_impl.hpp opm/autodiff/MissingFeatures.hpp + opm/autodiff/ThreadHandle.hpp opm/polymer/CompressibleTpfaPolymer.hpp opm/polymer/GravityColumnSolverPolymer.hpp opm/polymer/GravityColumnSolverPolymer_impl.hpp @@ -243,6 +256,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/ParallelFileMerger.hpp opm/simulators/SimulatorCompressibleTwophase.hpp opm/simulators/SimulatorIncompTwophase.hpp + opm/simulators/thresholdPressures.hpp opm/simulators/WellSwitchingLogger.hpp ) diff --git a/dune.module b/dune.module index 9b6ca05d6..026125417 100644 --- a/dune.module +++ b/dune.module @@ -10,4 +10,4 @@ Label: 2017.04-pre Maintainer: atgeirr@sintef.no MaintainerName: Atgeirr F. Rasmussen Url: http://opm-project.org -Depends: opm-common opm-parser opm-output opm-material opm-core opm-grid dune-istl (>=2.2) +Depends: opm-common opm-parser opm-output opm-material opm-core opm-grid dune-istl (>=2.2) ewoms diff --git a/examples/flow_ebos.cpp b/examples/flow_ebos.cpp new file mode 100644 index 000000000..9ff12a482 --- /dev/null +++ b/examples/flow_ebos.cpp @@ -0,0 +1,39 @@ +/* + Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 IRIS AS + + 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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include +#include +#include + + +// ----------------- Main program ----------------- +int main(int argc, char** argv) +{ + Opm::FlowMainEbos mainfunc; + return mainfunc.execute(argc, argv); +} diff --git a/jenkins/build.sh b/jenkins/build.sh index f785f1c8d..04158ebc5 100755 --- a/jenkins/build.sh +++ b/jenkins/build.sh @@ -7,7 +7,8 @@ upstreams=(opm-common opm-output opm-material opm-core - opm-grid) + opm-grid + ewoms) declare -A upstreamRev upstreamRev[opm-common]=master @@ -17,6 +18,7 @@ upstreamRev[opm-material]=master upstreamRev[opm-core]=master upstreamRev[opm-grid]=master upstreamRev[opm-output]=master +upstreamRev[ewoms]=master if grep -q "opm-common=" <<< $ghprbCommentBody then diff --git a/jenkins/run-norne.sh b/jenkins/run-norne.sh index cac96c625..ce12723a7 100755 --- a/jenkins/run-norne.sh +++ b/jenkins/run-norne.sh @@ -2,10 +2,11 @@ pushd . cd deps/opm-data +test -z $SIM && SIM=flow # Run the norne case cd norne -$WORKSPACE/$configuration/build-opm-simulators/bin/flow deck_filename=NORNE_ATW2013.DATA output_dir=OPM +$WORKSPACE/$configuration/build-opm-simulators/bin/$SIM deck_filename=NORNE_ATW2013.DATA output_dir=OPM test $? -eq 0 || exit 1 PATH=$WORKSPACE/$configuration/install/bin:$PATH ./plotwells.sh diff --git a/opm/autodiff/BlackoilDetails.hpp b/opm/autodiff/BlackoilDetails.hpp new file mode 100644 index 000000000..7422c2b73 --- /dev/null +++ b/opm/autodiff/BlackoilDetails.hpp @@ -0,0 +1,363 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2014, 2015 Statoil ASA. + Copyright 2015 NTNU + Copyright 2015 IRIS AS + + 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 . +*/ + +#ifndef OPM_BLACKOILDETAILS_HEADER_INCLUDED +#define OPM_BLACKOILDETAILS_HEADER_INCLUDED + +#include + +#include +#include + +namespace Opm { +namespace detail { + + + inline + std::vector + buildAllCells(const int nc) + { + std::vector all_cells(nc); + + for (int c = 0; c < nc; ++c) { all_cells[c] = c; } + + return all_cells; + } + + + + template + std::vector + activePhases(const PU& pu) + { + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + std::vector active(maxnp, false); + + for (int p = 0; p < pu.MaxNumPhases; ++p) { + active[ p ] = pu.phase_used[ p ] != 0; + } + + return active; + } + + + + template + std::vector + active2Canonical(const PU& pu) + { + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + std::vector act2can(maxnp, -1); + + for (int phase = 0; phase < maxnp; ++phase) { + if (pu.phase_used[ phase ]) { + act2can[ pu.phase_pos[ phase ] ] = phase; + } + } + + return act2can; + } + + + + inline + double getGravity(const double* g, const int dim) { + double grav = 0.0; + if (g) { + // Guard against gravity in anything but last dimension. + for (int dd = 0; dd < dim - 1; ++dd) { + assert(g[dd] == 0.0); + } + grav = g[dim - 1]; + } + return grav; + } + + /// \brief Compute the L-infinity norm of a vector + /// \warn This function is not suitable to compute on the well equations. + /// \param a The container to compute the infinity norm on. + /// It has to have one entry for each cell. + /// \param info In a parallel this holds the information about the data distribution. + template + inline + double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() ) + { + static_cast(pinfo); // Suppress warning in non-MPI case. +#if HAVE_MPI + if ( pinfo.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& real_info = + boost::any_cast(pinfo); + double result=0; + real_info.computeReduction(a.value(), Reduction::makeLInfinityNormFunctor(), result); + return result; + } + else +#endif + { + if( a.value().size() > 0 ) { + return a.value().matrix().template lpNorm (); + } + else { // this situation can occur when no wells are present + return 0.0; + } + } + } + + /// \brief Compute the Euclidian norm of a vector + /// \warning In the case that num_components is greater than 1 + /// an interleaved ordering is assumed. E.g. for each cell + /// all phases of that cell are stored consecutively. First + /// the ones for cell 0, then the ones for cell 1, ... . + /// \param it begin iterator for the given vector + /// \param end end iterator for the given vector + /// \param num_components number of components (i.e. phases) in the vector + /// \param pinfo In a parallel this holds the information about the data distribution. + template + inline + double euclidianNormSquared( Iterator it, const Iterator end, int num_components, const boost::any& pinfo = boost::any() ) + { + static_cast(num_components); // Suppress warning in the serial case. + static_cast(pinfo); // Suppress warning in non-MPI case. +#if HAVE_MPI + if ( pinfo.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(pinfo); + typedef typename Iterator::value_type Scalar; + Scalar product = 0.0; + int size_per_component = (end - it); + size_per_component /= num_components; // two lines to supresse unused warning. + assert((end - it) == num_components * size_per_component); + + if( num_components == 1 ) + { + auto component_container = + boost::make_iterator_range(it, end); + info.computeReduction(component_container, + Opm::Reduction::makeInnerProductFunctor(), + product); + } + else + { + auto& maskContainer = info.getOwnerMask(); + auto mask = maskContainer.begin(); + assert(static_cast(maskContainer.size()) == size_per_component); + + for(int cell = 0; cell < size_per_component; ++cell, ++mask) + { + Scalar cell_product = (*it) * (*it); + ++it; + for(int component=1; component < num_components; + ++component, ++it) + { + cell_product += (*it) * (*it); + } + product += cell_product * (*mask); + } + } + return info.communicator().sum(product); + } + else +#endif + { + double product = 0.0 ; + for( ; it != end; ++it ) { + product += ( *it * *it ); + } + return product; + } + } + /// \brief Compute the reduction within the convergence check. + /// \param[in] B A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. B.col(i) contains the values + /// for phase i. + /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. tempV.col(i) contains the + /// values + /// for phase i. + /// \param[in] R A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. B.col(i) contains the values + /// for phase i. + /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum + /// of R for the phase i. + /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the + /// maximum of tempV for the phase i. + /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average + /// of B for the phase i. + /// \param[out] maxNormWell The maximum of the well flux equations for each phase. + /// \param[in] nc The number of cells of the local grid. + /// \return The total pore volume over all cells. + inline + double + convergenceReduction(const Eigen::Array& B, + const Eigen::Array& tempV, + const Eigen::Array& R, + std::vector& R_sum, + std::vector& maxCoeff, + std::vector& B_avg, + std::vector& maxNormWell, + int nc, + int np, + const std::vector pv, + std::vector residual_well) + { + const int nw = residual_well.size() / np; + assert(nw * np == int(residual_well.size())); + + // Do the global reductions +#if 0 // HAVE_MPI + if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(linsolver_.parallelInformation()); + + // Compute the global number of cells and porevolume + std::vector v(nc, 1); + auto nc_and_pv = std::tuple(0, 0.0); + auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), + Opm::Reduction::makeGlobalSumFunctor()); + auto nc_and_pv_containers = std::make_tuple(v, pv); + info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv); + + for ( int idx = 0; idx < np; ++idx ) + { + auto values = std::tuple(0.0 ,0.0 ,0.0); + auto containers = std::make_tuple(B.col(idx), + tempV.col(idx), + R.col(idx)); + auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), + Opm::Reduction::makeGlobalMaxFunctor(), + Opm::Reduction::makeGlobalSumFunctor()); + info.computeReduction(containers, operators, values); + B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv); + maxCoeff[idx] = std::get<1>(values); + R_sum[idx] = std::get<2>(values); + assert(np >= np); + if (idx < np) { + maxNormWell[idx] = 0.0; + for ( int w = 0; w < nw; ++w ) { + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w])); + } + } + } + info.communicator().max(maxNormWell.data(), np); + // Compute pore volume + return std::get<1>(nc_and_pv); + } + else +#endif + { + B_avg.resize(np); + maxCoeff.resize(np); + R_sum.resize(np); + maxNormWell.resize(np); + for ( int idx = 0; idx < np; ++idx ) + { + B_avg[idx] = B.col(idx).sum()/nc; + maxCoeff[idx] = tempV.col(idx).maxCoeff(); + R_sum[idx] = R.col(idx).sum(); + + assert(np >= np); + if (idx < np) { + maxNormWell[idx] = 0.0; + for ( int w = 0; w < nw; ++w ) { + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w])); + } + } + } + // Compute total pore volume + return std::accumulate(pv.begin(), pv.end(), 0.0); + } + } + + + template + inline + double + convergenceReduction(const std::vector< std::vector< Scalar > >& B, + const std::vector< std::vector< Scalar > >& tempV, + const std::vector< std::vector< Scalar > >& R, + std::vector< Scalar >& R_sum, + std::vector< Scalar >& maxCoeff, + std::vector< Scalar >& B_avg, + std::vector< Scalar >& maxNormWell, + const int nc, + const int np, + const std::vector< Scalar >& pv, + const std::vector< Scalar >& residual_well) + { + const int nw = residual_well.size() / np; + assert(nw * np == int(residual_well.size())); + + // Do the global reductions + { + B_avg.resize(np); + maxCoeff.resize(np); + R_sum.resize(np); + maxNormWell.resize(np); + for ( int idx = 0; idx < np; ++idx ) + { + B_avg[idx] = std::accumulate( B[ idx ].begin(), B[ idx ].end(), 0.0 ) / nc; + R_sum[idx] = std::accumulate( R[ idx ].begin(), R[ idx ].end(), 0.0 ); + maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() )); + + assert(np >= np); + if (idx < np) { + maxNormWell[idx] = 0.0; + for ( int w = 0; w < nw; ++w ) { + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w])); + } + } + } + // Compute total pore volume + return std::accumulate(pv.begin(), pv.end(), 0.0); + } + } + + /// \brief Compute the L-infinity norm of a vector representing a well equation. + /// \param a The container to compute the infinity norm on. + /// \param info In a parallel this holds the information about the data distribution. + template + inline + double infinityNormWell( const ADB& a, const boost::any& pinfo ) + { + static_cast(pinfo); // Suppress warning in non-MPI case. + double result=0; + if( a.value().size() > 0 ) { + result = a.value().matrix().template lpNorm (); + } +#if HAVE_MPI + if ( pinfo.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& real_info = + boost::any_cast(pinfo); + result = real_info.communicator().max(result); + } +#endif + return result; + } + } // namespace detail +} // namespace Opm + +#endif // OPM_BLACKOILDETAILS_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 87a96e696..3f4e01e66 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -33,6 +33,8 @@ #include #include #include +#include +#include #include #include @@ -49,47 +51,6 @@ namespace Opm { class VFPProperties; class SimulationDataContainer; - /// Struct for containing iteration variables. - struct DefaultBlackoilSolutionState - { - typedef AutoDiffBlock ADB; - explicit DefaultBlackoilSolutionState(const int np) - : pressure ( ADB::null()) - , temperature( ADB::null()) - , saturation(np, ADB::null()) - , rs ( ADB::null()) - , rv ( ADB::null()) - , qs ( ADB::null()) - , bhp ( ADB::null()) - , canonical_phase_pressures(3, ADB::null()) - { - } - ADB pressure; - ADB temperature; - std::vector saturation; - ADB rs; - ADB rv; - ADB qs; - ADB bhp; - // Below are quantities stored in the state for optimization purposes. - std::vector canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active. - }; - - - - - /// Class used for reporting the outcome of a nonlinearIteration() call. - struct IterationReport - { - bool failed; - bool converged; - int linear_iterations; - int well_iterations; - }; - - - - /// Traits to encapsulate the types used by classes using or /// extending this model. Forward declared here, must be /// specialised for each concrete model class. diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index dbbcc7322..cd989f39b 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -25,6 +25,7 @@ #define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED #include +#include #include #include @@ -93,73 +94,6 @@ typedef Eigen::Array DataBlock; - -namespace detail { - - - inline - std::vector - buildAllCells(const int nc) - { - std::vector all_cells(nc); - - for (int c = 0; c < nc; ++c) { all_cells[c] = c; } - - return all_cells; - } - - - - template - std::vector - activePhases(const PU& pu) - { - const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - std::vector active(maxnp, false); - - for (int p = 0; p < pu.MaxNumPhases; ++p) { - active[ p ] = pu.phase_used[ p ] != 0; - } - - return active; - } - - - - template - std::vector - active2Canonical(const PU& pu) - { - const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - std::vector act2can(maxnp, -1); - - for (int phase = 0; phase < maxnp; ++phase) { - if (pu.phase_used[ phase ]) { - act2can[ pu.phase_pos[ phase ] ] = phase; - } - } - - return act2can; - } - - - - inline - double getGravity(const double* g, const int dim) { - double grav = 0.0; - if (g) { - // Guard against gravity in anything but last dimension. - for (int dd = 0; dd < dim - 1; ++dd) { - assert(g[dd] == 0.0); - } - grav = g[dim - 1]; - } - return grav; - } - -} // namespace detail - - template BlackoilModelBase:: BlackoilModelBase(const ModelParameters& param, @@ -1155,135 +1089,6 @@ namespace detail { return linsolver_.computeNewtonIncrement(residual_); } - - - - - namespace detail - { - /// \brief Compute the L-infinity norm of a vector - /// \warn This function is not suitable to compute on the well equations. - /// \param a The container to compute the infinity norm on. - /// It has to have one entry for each cell. - /// \param info In a parallel this holds the information about the data distribution. - inline - double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() ) - { - static_cast(pinfo); // Suppress warning in non-MPI case. -#if HAVE_MPI - if ( pinfo.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& real_info = - boost::any_cast(pinfo); - double result=0; - real_info.computeReduction(a.value(), Reduction::makeLInfinityNormFunctor(), result); - return result; - } - else -#endif - { - if( a.value().size() > 0 ) { - return a.value().matrix().lpNorm (); - } - else { // this situation can occur when no wells are present - return 0.0; - } - } - } - - /// \brief Compute the Euclidian norm of a vector - /// \warning In the case that num_components is greater than 1 - /// an interleaved ordering is assumed. E.g. for each cell - /// all phases of that cell are stored consecutively. First - /// the ones for cell 0, then the ones for cell 1, ... . - /// \param it begin iterator for the given vector - /// \param end end iterator for the given vector - /// \param num_components number of components (i.e. phases) in the vector - /// \param pinfo In a parallel this holds the information about the data distribution. - template - inline - double euclidianNormSquared( Iterator it, const Iterator end, int num_components, const boost::any& pinfo = boost::any() ) - { - static_cast(num_components); // Suppress warning in the serial case. - static_cast(pinfo); // Suppress warning in non-MPI case. -#if HAVE_MPI - if ( pinfo.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& info = - boost::any_cast(pinfo); - typedef typename Iterator::value_type Scalar; - Scalar product = 0.0; - int size_per_component = (end - it); - size_per_component /= num_components; // two lines to supresse unused warning. - assert((end - it) == num_components * size_per_component); - - if( num_components == 1 ) - { - auto component_container = - boost::make_iterator_range(it, end); - info.computeReduction(component_container, - Opm::Reduction::makeInnerProductFunctor(), - product); - } - else - { - auto& maskContainer = info.getOwnerMask(); - auto mask = maskContainer.begin(); - assert(static_cast(maskContainer.size()) == size_per_component); - - for(int cell = 0; cell < size_per_component; ++cell, ++mask) - { - Scalar cell_product = (*it) * (*it); - ++it; - for(int component=1; component < num_components; - ++component, ++it) - { - cell_product += (*it) * (*it); - } - product += cell_product * (*mask); - } - } - return info.communicator().sum(product); - } - else -#endif - { - double product = 0.0 ; - for( ; it != end; ++it ) { - product += ( *it * *it ); - } - return product; - } - } - - /// \brief Compute the L-infinity norm of a vector representing a well equation. - /// \param a The container to compute the infinity norm on. - /// \param info In a parallel this holds the information about the data distribution. - inline - double infinityNormWell( const ADB& a, const boost::any& pinfo ) - { - static_cast(pinfo); // Suppress warning in non-MPI case. - double result=0; - if( a.value().size() > 0 ) { - result = a.value().matrix().lpNorm (); - } -#if HAVE_MPI - if ( pinfo.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& real_info = - boost::any_cast(pinfo); - result = real_info.communicator().max(result); - } -#endif - return result; - } - - } // namespace detail - - - - - template void BlackoilModelBase:: diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp new file mode 100644 index 000000000..b3b34f6b3 --- /dev/null +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -0,0 +1,1459 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2014, 2015 Statoil ASA. + Copyright 2015 NTNU + Copyright 2015 IRIS AS + + 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 . +*/ + +#ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED +#define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +//#include + + + +namespace Ewoms { +namespace Properties { +NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); +SET_BOOL_PROP(EclFlowProblem, DisableWells, true); +SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); +}} + +namespace Opm { + + + namespace parameter { class ParameterGroup; } + class DerivedGeology; + class RockCompressibility; + class NewtonIterationBlackoilInterface; + class VFPProperties; + class SimulationDataContainer; + + + + + /// A model implementation for three-phase black oil. + /// + /// The simulator is capable of handling three-phase problems + /// where gas can be dissolved in oil and vice versa. It + /// uses an industry-standard TPFA discretization with per-phase + /// upwind weighting of mobilities. + class BlackoilModelEbos + { + typedef BlackoilModelEbos ThisType; + public: + // --------- Types and enums --------- + typedef BlackoilState ReservoirState; + typedef WellStateFullyImplicitBlackoilDense WellState; + typedef BlackoilModelParameters ModelParameters; + typedef DefaultBlackoilSolutionState SolutionState; + + typedef typename TTAG(EclFlowProblem) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator ; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables ; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; + + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; + + typedef ISTLSolver< MatrixBlockType, VectorBlockType > ISTLSolverType; + //typedef typename SolutionVector :: value_type PrimaryVariables ; + + struct FIPData { + enum FipId { + FIP_AQUA = Opm::Water, + FIP_LIQUID = Opm::Oil, + FIP_VAPOUR = Opm::Gas, + FIP_DISSOLVED_GAS = 3, + FIP_VAPORIZED_OIL = 4, + FIP_PV = 5, //< Pore volume + FIP_WEIGHTED_PRESSURE = 6 + }; + std::array, 7> fip; + }; + + // --------- Public methods --------- + + /// Construct the model. It will retain references to the + /// arguments of this functions, and they are expected to + /// remain in scope for the lifetime of the solver. + /// \param[in] param parameters + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] geo rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] wells well structure + /// \param[in] vfp_properties Vertical flow performance tables + /// \param[in] linsolver linear solver + /// \param[in] eclState eclipse state + /// \param[in] terminal_output request output to cout/cerr + BlackoilModelEbos(Simulator& ebosSimulator, + const ModelParameters& param, + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const StandardWellsDense& well_model, + const NewtonIterationBlackoilInterface& linsolver, + const bool terminal_output) + : ebosSimulator_(ebosSimulator) + , grid_(ebosSimulator_.gridManager().grid()) + , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) ) + , fluid_ (fluid) + , geo_ (geo) + , vfp_properties_( + eclState().getTableManager().getVFPInjTables(), + eclState().getTableManager().getVFPProdTables()) + , active_(detail::activePhases(fluid.phaseUsage())) + , has_disgas_(FluidSystem::enableDissolvedGas()) + , has_vapoil_(FluidSystem::enableVaporizedOil()) + , param_( param ) + , well_model_ (well_model) + , terminal_output_ (terminal_output) + , current_relaxation_(1.0) + , dx_old_(AutoDiffGrid::numCells(grid_)) + , isBeginReportStep_(false) + , invalidateIntensiveQuantitiesCache_(true) + { + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const std::vector pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + const std::vector depth(geo_.z().data(), geo_.z().data() + geo_.z().size()); + well_model_.init(&fluid_, &active_, &vfp_properties_, gravity, depth, pv); + wellModel().setWellsActive( localWellsActive() ); + global_nc_ = Opm::AutoDiffGrid::numCells(grid_); + // compute global sum of number of cells + global_nc_ = grid_.comm().sum( global_nc_ ); + + if (!istlSolver_) + { + OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed"); + } + } + + bool + isParallel() const + { + #if HAVE_MPI + if ( istlSolver().parallelInformation().type() != + typeid(ParallelISTLInformation) ) + { + return false; + } + else + { + const auto& comm =boost::any_cast(istlSolver().parallelInformation()).communicator(); + return comm.size() > 1; + } + #else + return false; + #endif + } + + + const EclipseState& eclState() const + { return *ebosSimulator_.gridManager().eclState(); } + + /// Called once before each time step. + /// \param[in] timer simulation timer + /// \param[in, out] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + void prepareStep(const SimulatorTimerInterface& /*timer*/, + const ReservoirState& /*reservoir_state*/, + const WellState& /* well_state */) + { + } + + + /// Called once per nonlinear iteration. + /// This model will perform a Newton-Raphson update, changing reservoir_state + /// and well_state. It will also use the nonlinear_solver to do relaxation of + /// updates if necessary. + /// \param[in] iteration should be 0 for the first call of a new timestep + /// \param[in] timer simulation timer + /// \param[in] nonlinear_solver nonlinear solver used (for oscillation/relaxation control) + /// \param[in, out] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + template + IterationReport nonlinearIteration(const int iteration, + const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver, + ReservoirState& reservoir_state, + WellState& well_state) + { + if (iteration == 0) { + // For each iteration we store in a vector the norms of the residual of + // the mass balance for each active phase, the well flux and the well equations. + residual_norms_history_.clear(); + current_relaxation_ = 1.0; + dx_old_ = 0.0; + } + + // reset intensive quantities cache useless other options are set + // further down + invalidateIntensiveQuantitiesCache_ = true; + + IterationReport iter_report = assemble(timer, iteration, reservoir_state, well_state); + std::vector residual_norms; + const bool converged = getConvergence(timer, iteration,residual_norms); + residual_norms_history_.push_back(residual_norms); + bool must_solve = (iteration < nonlinear_solver.minIter()) || (!converged); + + // don't solve if we have reached the maximum number of iteration. + must_solve = must_solve && (iteration < nonlinear_solver.maxIter()); + if (must_solve) { + + // enable single precision for solvers when dt is smaller then 20 days + //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ; + + // Compute the nonlinear update. + const int nc = AutoDiffGrid::numCells(grid_); + const int nw = wellModel().wells().number_of_wells; + BVector x(nc); + BVector xw(nw); + + solveJacobianSystem(x, xw); + + // Stabilize the nonlinear update. + bool isOscillate = false; + bool isStagnate = false; + nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate); + if (isOscillate) { + current_relaxation_ -= nonlinear_solver.relaxIncrement(); + current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); + if (terminalOutputEnabled()) { + std::string msg = " Oscillating behavior detected: Relaxation set to " + + std::to_string(current_relaxation_); + OpmLog::info(msg); + } + } + nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_); + + // Apply the update, applying model-dependent + // limitations and chopping of the update. + updateState(x,reservoir_state); + wellModel().updateWellState(xw, well_state); + + // since the solution was changed, the cache for the intensive quantities are invalid + // ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); + } + + if( converged && (iteration >= nonlinear_solver.minIter()) ) + { + // in case of convergence we do not need to reset intensive quantities + invalidateIntensiveQuantitiesCache_ = false ; + } + + const bool failed = false; // Not needed in this model. + const int linear_iters = must_solve ? linearIterationsLastSolve() : 0; + return IterationReport{ failed, converged, linear_iters, iter_report.well_iterations }; + } + + void printIf(int c, double x, double y, double eps, std::string type) { + if (std::abs(x-y) > eps) { + std::cout << type << " " < p0 ( previous.pressure() ); + std::vector< double > sat0( previous.saturation() ); + + const std::size_t pSize = p0.size(); + const std::size_t satSize = sat0.size(); + + // compute u^n - u^n+1 + for( std::size_t i=0; i 0.0 ) { + return stateOld / stateNew ; + } + else { + return 0.0; + } + } + + + /// The size (number of unknowns) of the nonlinear system of equations. + int sizeNonLinear() const + { + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int nw = wellModel().wells().number_of_wells; + return numPhases() * (nc + nw); + } + + /// Number of linear iterations used in last call to solveJacobianSystem(). + int linearIterationsLastSolve() const + { + return istlSolver().iterations(); + } + + template + void applyWellModelAdd(const X& x, Y& y ) + { + wellModel().apply(x, y); + } + + template + void applyWellModelScaleAdd(const Scalar alpha, const X& x, Y& y ) + { + wellModel().applyScaleAdd(alpha, x, y); + } + + + /// Solve the Jacobian system Jx = r where J is the Jacobian and + /// r is the residual. + void solveJacobianSystem(BVector& x, BVector& xw) const + { + const auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + + // apply well residual to the residual. + wellModel().apply(ebosResid); + + // set initial guess + x = 0.0; + + // Solve system. + if( isParallel() ) + { + typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, true > Operator; + Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() ); + assert( opA.comm() ); + istlSolver().solve( opA, x, ebosResid, *(opA.comm()) ); + } + else + { + typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, false > Operator; + Operator opA(ebosJac, const_cast< ThisType& > (*this) ); + istlSolver().solve( opA, x, ebosResid ); + } + + // recover wells. + xw = 0.0; + wellModel().recoverVariable(x, xw); + } + + //===================================================================== + // Implementation for ISTL-matrix based operator + //===================================================================== + + /*! + \brief Adapter to turn a matrix into a linear operator. + + Adapts a matrix to the assembled linear operator interface + */ + template + class WellModelMatrixAdapter : public Dune::AssembledLinearOperator + { + typedef Dune::AssembledLinearOperator BaseType; + + public: + typedef M matrix_type; + typedef X domain_type; + typedef Y range_type; + typedef typename X::field_type field_type; + +#if HAVE_MPI + typedef Dune::OwnerOverlapCopyCommunication communication_type; +#else + typedef Dune::CollectiveCommunication< Grid > communication_type; +#endif + + enum { + //! \brief The solver category. + category = overlapping ? + Dune::SolverCategory::overlapping : + Dune::SolverCategory::sequential + }; + + //! constructor: just store a reference to a matrix + WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation = boost::any() ) + : A_( A ), wellMod_( wellMod ), comm_() + { +#if HAVE_MPI + if( parallelInformation.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation); + comm_.reset( new communication_type( info.communicator() ) ); + } +#endif + } + + virtual void apply( const X& x, Y& y ) const + { + A_.mv( x, y ); + // add well model modification to y + wellMod_.applyWellModelAdd(x, y ); + +#if HAVE_MPI + if( comm_ ) + comm_->project( y ); +#endif + } + + // y += \alpha * A * x + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const + { + A_.usmv(alpha,x,y); + // add scaled well model modification to y + wellMod_.applyWellModelScaleAdd( alpha, x, y ); + +#if HAVE_MPI + if( comm_ ) + comm_->project( y ); +#endif + } + + virtual const matrix_type& getmat() const { return A_; } + + communication_type* comm() + { + return comm_.operator->(); + } + + protected: + const matrix_type& A_ ; + WellModel& wellMod_; + std::unique_ptr< communication_type > comm_; + }; + + /// Apply an update to the primary variables, chopped if appropriate. + /// \param[in] dx updates to apply to primary variables + /// \param[in, out] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + void updateState(const BVector& dx, + ReservoirState& reservoir_state) + { + using namespace Opm::AutoDiffGrid; + const int np = fluid_.numPhases(); + const int nc = numCells(grid_); + + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const double& dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)]; + //reservoir_state.pressure()[cell_idx] -= dp; + double& p = reservoir_state.pressure()[cell_idx]; + const double& dp_rel_max = dpMaxRel(); + const int sign_dp = dp > 0 ? 1: -1; + p -= sign_dp * std::min(std::abs(dp), std::abs(p)*dp_rel_max); + p = std::max(p, 0.0); + + // Saturation updates. + const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0; + const int xvar_ind = active_[Water] ? 2 : 1; + const double dxvar = active_[Gas] ? dx[cell_idx][flowPhaseToEbosCompIdx(xvar_ind)] : 0.0; + + double dso = 0.0; + double dsg = 0.0; + double drs = 0.0; + double drv = 0.0; + + double maxVal = 0.0; + // water phase + maxVal = std::max(std::abs(dsw),maxVal); + dso -= dsw; + // gas phase + switch (reservoir_state.hydroCarbonState()[cell_idx]) { + case HydroCarbonState::GasAndOil: + dsg = dxvar; + break; + case HydroCarbonState::OilOnly: + drs = dxvar; + break; + case HydroCarbonState::GasOnly: + dsg -= dsw; + drv = dxvar; + break; + default: + OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << reservoir_state.hydroCarbonState()[cell_idx]); + } + dso -= dsg; + + // Appleyard chop process. + maxVal = std::max(std::abs(dsg),maxVal); + double step = dsMax()/maxVal; + step = std::min(step, 1.0); + + + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + if (active_[Water]) { + double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; + sw -= step * dsw; + } + if (active_[Gas]) { + double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; + sg -= step * dsg; + } + double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]]; + so -= step * dso; + + // const double drmaxrel = drMaxRel(); + // Update rs and rv + if (has_disgas_) { + double& rs = reservoir_state.gasoilratio()[cell_idx]; + rs -= drs; + rs = std::max(rs, 0.0); + + } + if (has_vapoil_) { + double& rv = reservoir_state.rv()[cell_idx]; + rv -= drv; + rv = std::max(rv, 0.0); + } + + // Sg is used as primal variable for water only cells. + const double epsilon = 1e-4; //std::sqrt(std::numeric_limits::epsilon()); + double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; + double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; + double& rs = reservoir_state.gasoilratio()[cell_idx]; + double& rv = reservoir_state.rv()[cell_idx]; + + + // phase translation sg <-> rs + const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx]; + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + switch (hydroCarbonState) { + case HydroCarbonState::GasAndOil: { + + if (sw > (1.0 - epsilon)) // water only i.e. do nothing + break; + + if (sg <= 0.0 && has_disgas_) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs + sg = 0; + so = 1.0 - sw - sg; + const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + double& rs = reservoir_state.gasoilratio()[cell_idx]; + rs = rsSat*(1-epsilon); + } else if (so <= 0.0 && has_vapoil_) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv + so = 0; + sg = 1.0 - sw - so; + double& rv = reservoir_state.rv()[cell_idx]; + // use gas pressure? + const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + rv = rvSat*(1-epsilon); + } + break; + } + case HydroCarbonState::OilOnly: { + if (sw > (1.0 - epsilon)) { + // water only change to Sg + rs = 0; + rv = 0; + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + break; + } + + + + + const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + if (rs > ( rsSat * (1+epsilon) ) ) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + sg = epsilon; + so -= epsilon; + rs = rsSat; + } + break; + } + case HydroCarbonState::GasOnly: { + if (sw > (1.0 - epsilon)) { + // water only change to Sg + rs = 0; + rv = 0; + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + break; + } + const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + if (rv > rvSat * (1+epsilon) ) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + so = epsilon; + rv = rvSat; + sg -= epsilon; + } + break; + } + + default: + OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState); + } + } + + } + + + /// Return true if output to cout is wanted. + bool terminalOutputEnabled() const + { + return terminal_output_; + } + + template + double convergenceReduction(const CollectiveCommunication& comm, + const long int ncGlobal, + const int np, + const std::vector< std::vector< Scalar > >& B, + const std::vector< std::vector< Scalar > >& tempV, + const std::vector< std::vector< Scalar > >& R, + const std::vector< Scalar >& pv, + const std::vector< Scalar >& residual_well, + std::vector< Scalar >& R_sum, + std::vector< Scalar >& maxCoeff, + std::vector< Scalar >& B_avg, + std::vector< Scalar >& maxNormWell ) + { + const int nw = residual_well.size() / np; + assert(nw * np == int(residual_well.size())); + + // Do the global reductions + B_avg.resize(np); + maxCoeff.resize(np); + R_sum.resize(np); + maxNormWell.resize(np); + + // computation + for ( int idx = 0; idx < np; ++idx ) + { + B_avg[idx] = std::accumulate( B[ idx ].begin(), B[ idx ].end(), 0.0 ) / double(ncGlobal); + R_sum[idx] = std::accumulate( R[ idx ].begin(), R[ idx ].end(), 0.0 ); + maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() )); + + assert(np >= np); + if (idx < np) { + maxNormWell[idx] = 0.0; + for ( int w = 0; w < nw; ++w ) { + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w])); + } + } + } + + // Compute total pore volume + double pvSum = std::accumulate(pv.begin(), pv.end(), 0.0); + + if( comm.size() > 1 ) + { + // global reduction + std::vector< Scalar > sumBuffer; + std::vector< Scalar > maxBuffer; + sumBuffer.reserve( B_avg.size() + R_sum.size() + 1 ); + maxBuffer.reserve( maxCoeff.size() + maxNormWell.size() ); + for( int idx = 0; idx < np; ++idx ) + { + sumBuffer.push_back( B_avg[ idx ] ); + sumBuffer.push_back( R_sum[ idx ] ); + maxBuffer.push_back( maxCoeff[ idx ] ); + maxBuffer.push_back( maxNormWell[ idx ] ); + } + + // Compute total pore volume + sumBuffer.push_back( pvSum ); + + // compute global sum + comm.sum( sumBuffer.data(), sumBuffer.size() ); + + // compute global max + comm.max( maxBuffer.data(), maxBuffer.size() ); + + // restore values to local variables + for( int idx = 0, buffIdx = 0; idx < np; ++idx, ++buffIdx ) + { + B_avg[ idx ] = sumBuffer[ buffIdx ]; + maxCoeff[ idx ] = maxBuffer[ buffIdx ]; + ++buffIdx; + + R_sum[ idx ] = sumBuffer[ buffIdx ]; + maxNormWell[ idx ] = maxBuffer[ buffIdx ]; + } + + // restore global pore volume + pvSum = sumBuffer.back(); + } + + // return global pore volume + return pvSum; + } + + /// Compute convergence based on total mass balance (tol_mb) and maximum + /// residual mass balance (tol_cnv). + /// \param[in] timer simulation timer + /// \param[in] dt timestep length + /// \param[in] iteration current iteration number + bool getConvergence(const SimulatorTimerInterface& timer, const int iteration, std::vector& residual_norms) + { + typedef std::vector< double > Vector; + + const double dt = timer.currentStepLength(); + const double tol_mb = param_.tolerance_mb_; + const double tol_cnv = param_.tolerance_cnv_; + const double tol_wells = param_.tolerance_wells_; + + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int np = numPhases(); + + const auto& pv = geo_.poreVolume(); + + Vector R_sum(np); + Vector B_avg(np); + Vector maxCoeff(np); + Vector maxNormWell(np); + + std::vector< Vector > B( np, Vector( nc ) ); + std::vector< Vector > R( np, Vector( nc ) ); + std::vector< Vector > R2( np, Vector( nc ) ); + std::vector< Vector > tempV( np, Vector( nc ) ); + + const auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + + for ( int idx = 0; idx < np; ++idx ) + { + Vector& R2_idx = R2[ idx ]; + Vector& B_idx = B[ idx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + const int ebosCompIdx = flowPhaseToEbosCompIdx(idx); + + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + B_idx [cell_idx] = 1 / fs.invB(ebosPhaseIdx).value(); + R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; + } + } + + for ( int idx = 0; idx < np; ++idx ) + { + //tempV.col(idx) = R2.col(idx).abs()/pv; + Vector& tempV_idx = tempV[ idx ]; + Vector& R2_idx = R2[ idx ]; + for( int cell_idx = 0; cell_idx < nc; ++cell_idx ) + { + tempV_idx[ cell_idx ] = std::abs( R2_idx[ cell_idx ] ) / pv[ cell_idx ]; + } + } + + Vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + Vector wellResidual = wellModel().residual(); + + const double pvSum = convergenceReduction(grid_.comm(), global_nc_, np, + B, tempV, R2, pv_vector, wellResidual, + R_sum, maxCoeff, B_avg, maxNormWell ); + + Vector CNV(np); + Vector mass_balance_residual(np); + Vector well_flux_residual(np); + + bool converged_MB = true; + bool converged_CNV = true; + bool converged_Well = true; + // Finish computation + for ( int idx = 0; idx < np; ++idx ) + { + CNV[idx] = B_avg[idx] * dt * maxCoeff[idx]; + mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum; + converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb); + converged_CNV = converged_CNV && (CNV[idx] < tol_cnv); + // Well flux convergence is only for fluid phases, not other materials + // in our current implementation. + assert(np >= np); + if (idx < np) { + well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx]; + converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); + } + residual_norms.push_back(CNV[idx]); + } + + const bool converged = converged_MB && converged_CNV && converged_Well; + + if ( terminal_output_ ) + { + // Only rank 0 does print to std::cout + if (iteration == 0) { + std::string msg = "Iter"; + + std::vector< std::string > key( np ); + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + key[ phaseIdx ] = std::toupper( phaseName.front() ); + } + + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + msg += " MB(" + key[ phaseIdx ] + ") "; + } + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + msg += " CNV(" + key[ phaseIdx ] + ") "; + } + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + msg += " W-FLUX(" + key[ phaseIdx ] + ")"; + } + OpmLog::note(msg); + } + std::ostringstream ss; + const std::streamsize oprec = ss.precision(3); + const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); + ss << std::setw(4) << iteration; + for (int idx = 0; idx < np; ++idx) { + ss << std::setw(11) << mass_balance_residual[idx]; + } + for (int idx = 0; idx < np; ++idx) { + ss << std::setw(11) << CNV[idx]; + } + for (int idx = 0; idx < np; ++idx) { + ss << std::setw(11) << well_flux_residual[idx]; + } + ss.precision(oprec); + ss.flags(oflags); + OpmLog::note(ss.str()); + } + + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + + if (std::isnan(mass_balance_residual[phaseIdx]) + || std::isnan(CNV[phaseIdx]) + || (phaseIdx < np && std::isnan(well_flux_residual[phaseIdx]))) { + OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); + } + if (mass_balance_residual[phaseIdx] > maxResidualAllowed() + || CNV[phaseIdx] > maxResidualAllowed() + || (phaseIdx < np && well_flux_residual[phaseIdx] > maxResidualAllowed())) { + OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); + } + } + + return converged; + } + + + /// The number of active fluid phases in the model. + int numPhases() const + { + return fluid_.numPhases(); + } + + std::vector > + computeFluidInPlace(const std::vector& fipnum) const + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + //const ADB pv_mult = poroMult(pressure); + const auto& pv = geo_.poreVolume(); + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + + for (int i = 0; i<7; i++) { + fip_.fip[i].resize(nc,0.0); + } + + for (int c = 0; c < nc; ++c) { + const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + for (int phase = 0; phase < maxnp; ++phase) { + const double& b = fs.invB(flowPhaseToEbosPhaseIdx(phase)).value(); + const double& s = fs.saturation(flowPhaseToEbosPhaseIdx(phase)).value(); + const double pv_mult = 1.0; //todo + fip_.fip[phase][c] = pv_mult * b * s * pv[c]; + } + + if (active_[ Oil ] && active_[ Gas ]) { + // Account for gas dissolved in oil and vaporized oil + fip_.fip[FIPData::FIP_DISSOLVED_GAS][c] = fs.Rs().value() * fip_.fip[FIPData::FIP_LIQUID][c]; + fip_.fip[FIPData::FIP_VAPORIZED_OIL][c] = fs.Rv().value() * fip_.fip[FIPData::FIP_VAPOUR][c]; + } + } + + // For a parallel run this is just a local maximum and needs to be updated later + int dims = *std::max_element(fipnum.begin(), fipnum.end()); + std::vector> values(dims, std::vector(7,0.0)); + + std::vector hcpv(dims, 0.0); + std::vector pres(dims, 0.0); + + if ( !isParallel() ) + { + //Accumulate phases for each region + for (int phase = 0; phase < maxnp; ++phase) { + if (active_[ phase ]) { + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1) { + values[region][phase] += fip_.fip[phase][c]; + } + } + } + } + + //Accumulate RS and RV-volumes for each region + if (active_[ Oil ] && active_[ Gas ]) { + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1) { + values[region][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][c]; + values[region][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][c]; + } + } + } + + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1) { + const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); + hcpv[region] += pv[c] * hydrocarbon; + pres[region] += pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value(); + } + } + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1) { + + fip_.fip[FIPData::FIP_PV][c] = pv[c]; + const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); + + //Compute hydrocarbon pore volume weighted average pressure. + //If we have no hydrocarbon in region, use pore volume weighted average pressure instead + if (hcpv[region] != 0) { + fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region]; + } else { + fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c]; + } + + values[region][FIPData::FIP_PV] += fip_.fip[FIPData::FIP_PV][c]; + values[region][FIPData::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c]; + } + } + } + else + { +#if HAVE_MPI + // mask[c] is 1 if we need to compute something in parallel + const auto & pinfo = + boost::any_cast(istlSolver().parallelInformation()); + const auto& mask = pinfo.getOwnerMask(); + auto comm = pinfo.communicator(); + // Compute the global dims value and resize values accordingly. + dims = comm.max(dims); + values.resize(dims, std::vector(7,0.0)); + + //Accumulate phases for each region + for (int phase = 0; phase < maxnp; ++phase) { + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1 && mask[c]) { + values[region][phase] += fip_.fip[phase][c]; + } + } + } + + //Accumulate RS and RV-volumes for each region + if (active_[ Oil ] && active_[ Gas ]) { + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1 && mask[c]) { + values[region][FIPData::FIP_DISSOLVED_GAS] += fip_.fip[FIPData::FIP_DISSOLVED_GAS][c]; + values[region][FIPData::FIP_VAPORIZED_OIL] += fip_.fip[FIPData::FIP_VAPORIZED_OIL][c]; + } + } + } + + hcpv = std::vector(dims, 0.0); + pres = std::vector(dims, 0.0); + + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1 && mask[c]) { + const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); + hcpv[region] += pv[c] * hydrocarbon; + pres[region] += pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value(); + } + } + + comm.sum(hcpv.data(), hcpv.size()); + comm.sum(pres.data(), pres.size()); + + for (int c = 0; c < nc; ++c) { + const int region = fipnum[c] - 1; + if (region != -1 && mask[c]) { + fip_.fip[FIPData::FIP_PV][c] = pv[c]; + const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(c, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); + + if (hcpv[region] != 0) { + fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbon / hcpv[region]; + } else { + fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c]; + } + + values[region][FIPData::FIP_PV] += fip_.fip[FIPData::FIP_PV][c]; + values[region][FIPData::FIP_WEIGHTED_PRESSURE] += fip_.fip[FIPData::FIP_WEIGHTED_PRESSURE][c]; + } + } + + // For the frankenstein branch we hopefully can turn values into a vanilla + // std::vector, use some index magic above, use one communication + // to sum up the vector entries instead of looping over the regions. + for(int reg=0; reg < dims; ++reg) + { + comm.sum(values[reg].data(), values[reg].size()); + } +#else + // This should never happen! + OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel"); +#endif + } + + return values; + } + + const FIPData& getFIPData() const { + return fip_; + } + + + + const Simulator& ebosSimulator() const + { return ebosSimulator_; } + + protected: + const ISTLSolverType& istlSolver() const + { + assert( istlSolver_ ); + return *istlSolver_; + } + + + // --------- Data members --------- + + Simulator& ebosSimulator_; + const Grid& grid_; + const ISTLSolverType* istlSolver_; + const BlackoilPropsAdInterface& fluid_; + const DerivedGeology& geo_; + VFPProperties vfp_properties_; + // For each canonical phase -> true if active + const std::vector active_; + // Size = # active phases. Maps active -> canonical phase indices. + const std::vector cells_; // All grid cells + const bool has_disgas_; + const bool has_vapoil_; + + ModelParameters param_; + + // Well Model + StandardWellsDense well_model_; + + /// \brief Whether we print something to std::cout + bool terminal_output_; + /// \brief The number of cells of the global grid. + long int global_nc_; + + std::vector> residual_norms_history_; + double current_relaxation_; + BVector dx_old_; + mutable FIPData fip_; + + + + // --------- Protected methods --------- + + public: + + /// return the StandardWells object + StandardWellsDense& wellModel() { return well_model_; } + const StandardWellsDense& wellModel() const { return well_model_; } + + /// return the Well struct in the StandardWells + const Wells& wells() const { return well_model_.wells(); } + + /// return true if wells are available in the reservoir + bool wellsActive() const { return well_model_.wellsActive(); } + + /// return true if wells are available on this process + bool localWellsActive() const { return well_model_.localWellsActive(); } + + + void convertInput( const int iterationIdx, + const ReservoirState& reservoirState, + Simulator& simulator ) const + { + SolutionVector& solution = simulator.model().solution( 0 /* timeIdx */ ); + const Opm::PhaseUsage pu = fluid_.phaseUsage(); + + const int numCells = reservoirState.numCells(); + const int numPhases = fluid_.numPhases(); + const auto& oilPressure = reservoirState.pressure(); + const auto& saturations = reservoirState.saturation(); + const auto& rs = reservoirState.gasoilratio(); + const auto& rv = reservoirState.rv(); + for( int cellIdx = 0; cellIdx gas only with vaporized oil in the gas) is + // relatively expensive as it requires to compute the capillary + // pressure in order to get the gas phase pressure. (the reason why + // ebos uses the gas pressure here is that it makes the common case + // of the primary variable switching code fast because to determine + // whether the oil phase appears one needs to compute the Rv value + // for the saturated gas phase and if this is not available as a + // primary variable, it needs to be computed.) luckily for here, the + // gas-only case is not too common, so the performance impact of this + // is limited. + typedef Opm::SimpleModularFluidState SatOnlyFluidState; + SatOnlyFluidState fluidState; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]); + + double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; + const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx); + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]); + + cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx]; + cellPv[BlackoilIndices::pressureSwitchIdx] = pg; + cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv ); + } + else + { + assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil); + cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]]; + cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ]; + cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg ); + } + } + + if( iterationIdx == 0 ) + { + simulator.model().solution( 1 /* timeIdx */ ) = solution; + } + } + + public: + int ebosCompToFlowPhaseIdx( const int compIdx ) const + { + const int compToPhase[ 3 ] = { Oil, Water, Gas }; + return compToPhase[ compIdx ]; + } + + int flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + return flowToEbos[ flowPv ]; + } + + int flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + return phaseToComp[ phaseIdx ]; + } + + + + + private: + + void convertResults(BVector& ebosResid, Mat& ebosJac) const + { + const int numPhases = wells().number_of_phases; + const int numCells = ebosJac.N(); + assert( numCells == static_cast(ebosJac.M()) ); + + // write the right-hand-side values from the ebosJac into the objects + // allocated above. + const auto endrow = ebosJac.end(); + for( int cellIdx = 0; cellIdx < numCells; ++cellIdx ) + { + const double cellVolume = ebosSimulator_.model().dofTotalVolume(cellIdx); + auto& cellRes = ebosResid[ cellIdx ]; + + for( int flowPhaseIdx = 0; flowPhaseIdx < numPhases; ++flowPhaseIdx ) + { + const double refDens = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( flowPhaseIdx ), 0 ); + cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] /= refDens; + cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] *= cellVolume; + } + } + + for( auto row = ebosJac.begin(); row != endrow; ++row ) + { + const int rowIdx = row.index(); + const double cellVolume = ebosSimulator_.model().dofTotalVolume(rowIdx); + + + // translate the Jacobian of the residual from the format used by ebos to + // the one expected by flow + const auto endcol = row->end(); + for( auto col = row->begin(); col != endcol; ++col ) + { + for( int flowPhaseIdx = 0; flowPhaseIdx < numPhases; ++flowPhaseIdx ) + { + const double refDens = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( flowPhaseIdx ), 0 ); + for( int pvIdx=0; pvIdx("deck_filename"); } diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index 2f880cc16..558dae91f 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -20,6 +20,8 @@ #ifndef OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED #define OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED +#include + namespace Opm { @@ -59,6 +61,9 @@ namespace Opm /// Try to detect oscillation or stagnation. bool use_update_stabilization_; + // The file name of the deck + std::string deck_file_name_; + /// Construct from user parameters or defaults. explicit BlackoilModelParameters( const parameter::ParameterGroup& param ); diff --git a/opm/autodiff/DefaultBlackoilSolutionState.hpp b/opm/autodiff/DefaultBlackoilSolutionState.hpp new file mode 100644 index 000000000..06d7a4aef --- /dev/null +++ b/opm/autodiff/DefaultBlackoilSolutionState.hpp @@ -0,0 +1,57 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Statoil ASA. + Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services + Copyright 2015 NTNU + + 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 . +*/ +#ifndef OPM_DEFAULTBLACKOILSOLUTIONSTATE_HEADER_INCLUDED +#define OPM_DEFAULTBLACKOILSOLUTIONSTATE_HEADER_INCLUDED + +#include + +namespace Opm { + /// Struct for containing iteration variables. + struct DefaultBlackoilSolutionState + { + typedef AutoDiffBlock ADB; + explicit DefaultBlackoilSolutionState(const int np) + : pressure ( ADB::null()) + , temperature( ADB::null()) + , saturation(np, ADB::null()) + , rs ( ADB::null()) + , rv ( ADB::null()) + , qs ( ADB::null()) + , bhp ( ADB::null()) + , wellVariables ( ADB::null()) + , canonical_phase_pressures(3, ADB::null()) + { + } + ADB pressure; + ADB temperature; + std::vector saturation; + ADB rs; + ADB rv; + ADB qs; + ADB bhp; + ADB wellVariables; + // Below are quantities stored in the state for optimization purposes. + std::vector canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active. + }; +} // namespace Opm + +#endif // OPM_DEFAULTBLACKOILSOLUTIONSTATE_HEADER_INCLUDED diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index a254fd033..59cc99d3d 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -193,6 +193,7 @@ namespace Opm typedef BlackoilPropsAdFromDeck FluidProps; typedef FluidProps::MaterialLawManager MaterialLawManager; typedef typename Simulator::ReservoirState ReservoirState; + typedef typename Simulator::OutputWriter OutputWriter; // ------------ Data members ------------ @@ -229,7 +230,7 @@ namespace Opm boost::any parallel_information_; // setupOutputWriter() std::unique_ptr eclipse_writer_; - std::unique_ptr output_writer_; + std::unique_ptr output_writer_; // setupLinearSolver std::unique_ptr fis_solver_; // createSimulator() @@ -770,12 +771,12 @@ namespace Opm // create output writer after grid is distributed, otherwise the parallel output // won't work correctly since we need to create a mapping from the distributed to // the global view - output_writer_.reset(new BlackoilOutputWriter(grid_init_->grid(), - param_, - *eclipse_state_, - std::move(eclipse_writer_), - Opm::phaseUsageFromDeck(*deck_), - fluidprops_->permeability())); + output_writer_.reset(new OutputWriter(grid_init_->grid(), + param_, + *eclipse_state_, + std::move(eclipse_writer_), + Opm::phaseUsageFromDeck(*deck_), + fluidprops_->permeability())); } diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp new file mode 100644 index 000000000..f41fd098a --- /dev/null +++ b/opm/autodiff/FlowMainEbos.hpp @@ -0,0 +1,141 @@ +/* + Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 IRIS AS + Copyright 2014 STATOIL ASA. + + 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 . +*/ + +#ifndef OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED +#define OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED + +#include +#include + +#include + +namespace Opm +{ + // The FlowMain class is the ebos based black-oil simulator. + class FlowMainEbos : public FlowMainBase + { + protected: + typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator; + typedef FlowMainBase Base; + friend Base; + + typedef typename TTAG(EclFlowProblem) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator; + + // Print startup message if on output rank. + void printStartupMessage() + { + if (output_cout_) { + const int lineLen = 70; + const std::string version = moduleVersionName(); + const std::string banner = "This is flow_ebos (version "+version+")"; + const std::string ewomsVersion = "(eWoms version: " + Ewoms::versionString() + ")"; + const int bannerPreLen = (lineLen - 2 - banner.size())/2; + const int bannerPostLen = bannerPreLen + (lineLen - 2 - banner.size())%2; + const int eVPreLen = (lineLen - 2 - ewomsVersion.size())/2; + const int eVPostLen = eVPreLen + (lineLen - 2 - ewomsVersion.size())%2; + std::cout << "**********************************************************************\n"; + std::cout << "* *\n"; + std::cout << "*" << std::string(bannerPreLen, ' ') << banner << std::string(bannerPostLen, ' ') << "*\n"; + std::cout << "*" << std::string(eVPreLen, ' ') << ewomsVersion << std::string(eVPostLen, ' ') << "*\n"; + std::cout << "* *\n"; + std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n"; + std::cout << "* and is part of OPM. For more information see: *\n"; + std::cout << "* http://opm-project.org *\n"; + std::cout << "* *\n"; + std::cout << "**********************************************************************\n\n"; + } + } + + // Parser the input and creates the Deck and EclipseState objects. + // Writes to: + // deck_ + // eclipse_state_ + // May throw if errors are encountered, here configured to be somewhat tolerant. + void readDeckInput() + { + std::string progName("flow_ebos"); + std::string deckFile("--ecl-deck-file-name="); + deckFile += param_.get("deck_filename"); + char* ptr[2]; + ptr[ 0 ] = const_cast< char * > (progName.c_str()); + ptr[ 1 ] = const_cast< char * > (deckFile.c_str()); + EbosSimulator::registerParameters(); + Ewoms::setupParameters_< TypeTag > ( 2, ptr ); + ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false)); + ebosSimulator_->model().applyInitialSolution(); + + Base::deck_ = ebosSimulator_->gridManager().deck(); + Base::eclipse_state_ = ebosSimulator_->gridManager().eclState(); + IOConfig& ioConfig = Base::eclipse_state_->getIOConfig(); + ioConfig.setOutputDir(Base::output_dir_); + + // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step) + if (Base::param_.has("output_interval")) { + const int output_interval = Base::param_.get("output_interval"); + eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) ); + + } + + // Possible to force initialization only behavior (NOSIM). + if (Base::param_.has("nosim")) { + const bool nosim = Base::param_.get("nosim"); + ioConfig.overrideNOSIM( nosim ); + } + } + + // Setup linear solver. + // Writes to: + // fis_solver_ + void setupLinearSolver() + { + typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType; + Base::fis_solver_.reset( new ISTLSolverType( Base::param_, Base::parallel_information_ ) ); + } + + /// This is the main function of Flow. + // Create simulator instance. + // Writes to: + // simulator_ + void createSimulator() + { + // Create the simulator instance. + Base::simulator_.reset(new Simulator(*ebosSimulator_, + Base::param_, + *Base::geoprops_, + *Base::fluidprops_, + Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr, + *Base::fis_solver_, + Base::gravity_.data(), + Base::deck_->hasKeyword("DISGAS"), + Base::deck_->hasKeyword("VAPOIL"), + Base::eclipse_state_, + *Base::output_writer_, + Base::threshold_pressures_)); + } + + private: + std::unique_ptr ebosSimulator_; + }; +} // namespace Opm + +#endif // OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp new file mode 100644 index 000000000..a5132c4a1 --- /dev/null +++ b/opm/autodiff/ISTLSolver.hpp @@ -0,0 +1,413 @@ +/* + Copyright 2016 IRIS AS + + 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 . +*/ + +#ifndef OPM_ISTLSOLVER_HEADER_INCLUDED +#define OPM_ISTLSOLVER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +namespace Dune +{ + +namespace ISTLUtility { + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + +//! invert matrix by calling matrix.invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + matrix.invert(); +} + +} // end ISTLUtility + +template +class MatrixBlock : public Dune::FieldMatrix +{ +public: + typedef Dune::FieldMatrix BaseType; + + using BaseType :: operator= ; + using BaseType :: rows; + using BaseType :: cols; + explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {} + void invert() + { + ISTLUtility::invertMatrix( *this ); + } + const BaseType& asBase() const { return static_cast< const BaseType& > (*this); } + BaseType& asBase() { return static_cast< BaseType& > (*this); } +}; + +template +void +print_row (std::ostream& s, const MatrixBlock& A, + typename FieldMatrix::size_type I, + typename FieldMatrix::size_type J, + typename FieldMatrix::size_type therow, int width, + int precision) +{ + print_row(s, A.asBase(), I, J, therow, width, precision); +} + +template +K& firstmatrixelement (MatrixBlock& A) +{ + return firstmatrixelement( A.asBase() ); +} + + + +template +struct MatrixDimension< MatrixBlock< Scalar, n, m > > +: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType > +{ +}; + +} // end namespace Dune + +namespace Opm +{ + /// This class solves the fully implicit black-oil system by + /// solving the reduced system (after eliminating well variables) + /// as a block-structured matrix (one block for all cell variables) for a fixed + /// number of cell variables np . + template < class MatrixBlockType, class VectorBlockType > + class ISTLSolver : public NewtonIterationBlackoilInterface + { + typedef typename MatrixBlockType :: field_type Scalar; + + typedef Dune::BCRSMatrix Matrix; + typedef Dune::BlockVector Vector; + + public: + typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; + + typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector; + /// Construct a system solver. + /// \param[in] param parameters controlling the behaviour of the linear solvers + /// \param[in] parallelInformation In the case of a parallel run + /// with dune-istl the information about the parallelization. + ISTLSolver(const NewtonIterationBlackoilInterleavedParameters& param, + const boost::any& parallelInformation_arg=boost::any()) + : iterations_( 0 ), + parallelInformation_(parallelInformation_arg), + isIORank_(isIORank(parallelInformation_arg)), + parameters_( param ) + { + } + + /// Construct a system solver. + /// \param[in] param ParameterGroup controlling the behaviour of the linear solvers + /// \param[in] parallelInformation In the case of a parallel run + /// with dune-istl the information about the parallelization. + ISTLSolver(const parameter::ParameterGroup& param, + const boost::any& parallelInformation_arg=boost::any()) + : iterations_( 0 ), + parallelInformation_(parallelInformation_arg), + isIORank_(isIORank(parallelInformation_arg)), + parameters_( param ) + { + } + + // dummy method that is not implemented for this class + SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const + { + OPM_THROW(std::logic_error,"This method is not implemented"); + return SolutionVector(); + } + + /// Solve the system of linear equations Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] residual residual object containing A and b. + /// \return the solution x + + /// \copydoc NewtonIterationBlackoilInterface::iterations + int iterations () const { return iterations_; } + + /// \copydoc NewtonIterationBlackoilInterface::parallelInformation + const boost::any& parallelInformation() const { return parallelInformation_; } + + public: + /// \brief construct the CPR preconditioner and the solver. + /// \tparam P The type of the parallel information. + /// \param parallelInformation the information about the parallelization. + template + void constructPreconditionerAndSolve(LinearOperator& linearOperator, + Vector& x, Vector& istlb, + const POrComm& parallelInformation_arg, + Dune::InverseOperatorResult& result) const + { + // Construct scalar product. + typedef Dune::ScalarProductChooser ScalarProductChooser; + typedef std::unique_ptr SPPointer; + SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg)); + + // Communicate if parallel. + parallelInformation_arg.copyOwnerToAll(istlb, istlb); + +#if ! HAVE_UMFPACK + const bool useAmg = false ; + if( useAmg ) + { + typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; + typedef typename CPRSelectorType::AMG AMG; + typedef typename CPRSelectorType::Operator MatrixOperator; + + std::unique_ptr< AMG > amg; + std::unique_ptr< MatrixOperator > opA; + + if( ! std::is_same< LinearOperator, MatrixOperator > :: value ) + { + // create new operator in case linear operator and matrix operator differ + opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) ); + } + + const double relax = 1.0; + + // Construct preconditioner. + constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax ); + + // Solve. + solve(linearOperator, x, istlb, *sp, *amg, result); + } + else +#endif + { + // Construct preconditioner. + auto precond = constructPrecond(linearOperator, parallelInformation_arg); + + // Solve. + solve(linearOperator, x, istlb, *sp, *precond, result); + } + } + + typedef Dune::SeqILU0 SeqPreconditioner; + + template + std::unique_ptr constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const + { + const double relax = 0.9; + std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), relax)); + return precond; + } + +#if HAVE_MPI + typedef Dune::OwnerOverlapCopyCommunication Comm; + typedef ParallelOverlappingILU0 ParPreconditioner; + template + std::unique_ptr + constructPrecond(Operator& opA, const Comm& comm) const + { + typedef std::unique_ptr Pointer; + const double relax = 0.9; + return Pointer(new ParPreconditioner(opA.getmat(), comm, relax)); + } +#endif + + template + void + constructAMGPrecond(LinearOperator& linearOperator, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const + { + ISTLUtility::createAMGPreconditionerPointer( *opA, relax, comm, amg ); + } + + + template + void + constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const + { + ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg ); + } + + /// \brief Solve the system using the given preconditioner and scalar product. + template + void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const + { + // TODO: Revise when linear solvers interface opm-core is done + // Construct linear solver. + // GMRes solver + if ( parameters_.newton_use_gmres_ ) { + Dune::RestartedGMResSolver linsolve(opA, sp, precond, + parameters_.linear_solver_reduction_, + parameters_.linear_solver_restart_, + parameters_.linear_solver_maxiter_, + parameters_.linear_solver_verbosity_); + // Solve system. + linsolve.apply(x, istlb, result); + } + else { // BiCGstab solver + Dune::BiCGSTABSolver linsolve(opA, sp, precond, + parameters_.linear_solver_reduction_, + parameters_.linear_solver_maxiter_, + parameters_.linear_solver_verbosity_); + // Solve system. + linsolve.apply(x, istlb, result); + } + } + + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] A matrix A + /// \param[inout] x solution to be computed x + /// \param[in] b right hand side b + void solve(Matrix& A, Vector& x, Vector& b ) const + { + // Parallel version is deactivated until we figure out how to do it properly. +#if HAVE_MPI + if (parallelInformation_.type() == typeid(ParallelISTLInformation)) + { + typedef Dune::OwnerOverlapCopyCommunication Comm; + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation_); + Comm istlComm(info.communicator()); + + // Construct operator, scalar product and vectors needed. + typedef Dune::OverlappingSchwarzOperator Operator; + Operator opA(A, istlComm); + solve( opA, x, b, istlComm ); + } + else +#endif + { + // Construct operator, scalar product and vectors needed. + Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A ); + solve( opA, x, b ); + } + } + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] A matrix A + /// \param[inout] x solution to be computed x + /// \param[in] b right hand side b + template + void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const + { + Dune::InverseOperatorResult result; + // Parallel version is deactivated until we figure out how to do it properly. +#if HAVE_MPI + if (parallelInformation_.type() == typeid(ParallelISTLInformation)) + { + const size_t size = opA.getmat().N(); + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation_); + + // As we use a dune-istl with block size np the number of components + // per parallel is only one. + info.copyValuesTo(comm.indexSet(), comm.remoteIndices(), + size, 1); + // Construct operator, scalar product and vectors needed. + constructPreconditionerAndSolve(opA, x, b, comm, result); + } + else +#endif + { + OPM_THROW(std::logic_error,"this method if for parallel solve only"); + } + + checkConvergence( result ); + } + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] A matrix A + /// \param[inout] x solution to be computed x + /// \param[in] b right hand side b + template + void solve(Operator& opA, Vector& x, Vector& b ) const + { + Dune::InverseOperatorResult result; + // Construct operator, scalar product and vectors needed. + Dune::Amg::SequentialInformation info; + constructPreconditionerAndSolve(opA, x, b, info, result); + checkConvergence( result ); + } + + void checkConvergence( const Dune::InverseOperatorResult& result ) const + { + // store number of iterations + iterations_ = result.iterations; + + // Check for failure of linear solver. + if (!parameters_.ignoreConvergenceFailure_ && !result.converged) { + const std::string msg("Convergence failure for linear solver."); + if (isIORank_) { + OpmLog::problem(msg); + } + OPM_THROW_NOLOG(LinearSolverProblem, msg); + } + } + protected: + mutable int iterations_; + boost::any parallelInformation_; + bool isIORank_; + + NewtonIterationBlackoilInterleavedParameters parameters_; + }; // end ISTLSolver + +} // namespace Opm +#endif diff --git a/opm/autodiff/IterationReport.hpp b/opm/autodiff/IterationReport.hpp new file mode 100644 index 000000000..dc1ea55d0 --- /dev/null +++ b/opm/autodiff/IterationReport.hpp @@ -0,0 +1,37 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Statoil ASA. + Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services + Copyright 2015 NTNU + + 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 . +*/ + +#ifndef OPM_ITERATION_REPORT_HEADER_INCLUDED +#define OPM_ITERATION_REPORT_HEADER_INCLUDED + +namespace Opm { + /// Class used for reporting the outcome of a nonlinearIteration() call. + struct IterationReport + { + bool failed; + bool converged; + int linear_iterations; + int well_iterations; + }; +} // namespace Opm + +#endif // OPM_ITERATION_REPORT_HEADER_INCLUDED diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 83f855833..7568c203c 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -34,13 +34,9 @@ #include #include +#include + #include -#include -#include -#include -#include -#include -#include #if HAVE_UMFPACK #include @@ -49,90 +45,6 @@ #endif #include - -namespace Dune -{ - -namespace ISTLUtility { - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling FMatrixHelp::invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - FieldMatrix A ( matrix ); - FMatrixHelp::invertMatrix(A, matrix ); -} - -//! invert matrix by calling matrix.invert -template -static inline void invertMatrix (FieldMatrix &matrix) -{ - matrix.invert(); -} - -} // end ISTLUtility - -template -class MatrixBlock : public Dune::FieldMatrix -{ -public: - typedef Dune::FieldMatrix BaseType; - - using BaseType :: operator= ; - using BaseType :: rows; - using BaseType :: cols; - explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {} - void invert() - { - ISTLUtility::invertMatrix( *this ); - } - const BaseType& asBase() const { return static_cast< const BaseType& > (*this); } - BaseType& asBase() { return static_cast< BaseType& > (*this); } -}; - -template -void -print_row (std::ostream& s, const MatrixBlock& A, - typename FieldMatrix::size_type I, - typename FieldMatrix::size_type J, - typename FieldMatrix::size_type therow, int width, - int precision) -{ - print_row(s, A.asBase(), I, J, therow, width, precision); -} - -template -K& firstmatrixelement (MatrixBlock& A) -{ - return firstmatrixelement( A.asBase() ); -} - - - -template -struct MatrixDimension< MatrixBlock< Scalar, n, m > > -: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType > -{ -}; - -} // end namespace Dune - namespace Opm { @@ -161,9 +73,12 @@ namespace Opm typedef Dune::FieldVector VectorBlockType; typedef Dune::MatrixBlock MatrixBlockType; + typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector Vector; + typedef Opm::ISTLSolver< MatrixBlockType, VectorBlockType > ISTLSolverType; + public: typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector; /// Construct a system solver. @@ -172,9 +87,7 @@ namespace Opm /// with dune-istl the information about the parallelization. NewtonIterationBlackoilInterleavedImpl(const NewtonIterationBlackoilInterleavedParameters& param, const boost::any& parallelInformation_arg=boost::any()) - : iterations_( 0 ), - parallelInformation_(parallelInformation_arg), - isIORank_(isIORank(parallelInformation_arg)), + : istlSolver_( param, parallelInformation_arg ), parameters_( param ) { } @@ -186,110 +99,12 @@ namespace Opm /// \return the solution x /// \copydoc NewtonIterationBlackoilInterface::iterations - int iterations () const { return iterations_; } + int iterations () const { return istlSolver_.iterations(); } /// \copydoc NewtonIterationBlackoilInterface::parallelInformation - const boost::any& parallelInformation() const { return parallelInformation_; } + const boost::any& parallelInformation() const { return istlSolver_.parallelInformation(); } public: - /// \brief construct the CPR preconditioner and the solver. - /// \tparam P The type of the parallel information. - /// \param parallelInformation the information about the parallelization. - template - void constructPreconditionerAndSolve(O& opA, - Vector& x, Vector& istlb, - const POrComm& parallelInformation_arg, - Dune::InverseOperatorResult& result) const - { - // Construct scalar product. - typedef Dune::ScalarProductChooser ScalarProductChooser; - typedef std::unique_ptr SPPointer; - SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg)); - - // Communicate if parallel. - parallelInformation_arg.copyOwnerToAll(istlb, istlb); - -#if ! HAVE_UMFPACK - const bool useAmg = false ; - if( useAmg ) - { - typedef ISTLUtility::CPRSelector< Mat, Vector, Vector, POrComm> CPRSelectorType; - typedef typename CPRSelectorType::AMG AMG; - std::unique_ptr< AMG > amg; - // Construct preconditioner. - constructAMGPrecond(opA, parallelInformation_arg, amg); - - // Solve. - solve(opA, x, istlb, *sp, *amg, result); - } - else -#endif - { - // Construct preconditioner. - auto precond = constructPrecond(opA, parallelInformation_arg); - - // Solve. - solve(opA, x, istlb, *sp, *precond, result); - } - } - - typedef Dune::SeqILU0 SeqPreconditioner; - - template - std::unique_ptr constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const - { - const double relax = 0.9; - std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), relax)); - return precond; - } - -#if HAVE_MPI - typedef Dune::OwnerOverlapCopyCommunication Comm; - typedef ParallelOverlappingILU0 ParPreconditioner; - template - std::unique_ptr - constructPrecond(Operator& opA, const Comm& comm) const - { - typedef std::unique_ptr Pointer; - const double relax = 0.9; - return Pointer(new ParPreconditioner(opA.getmat(), comm, relax)); - } -#endif - - template - void - constructAMGPrecond(Operator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg ) const - { - const double relax = 1.0; - ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg ); - } - - /// \brief Solve the system using the given preconditioner and scalar product. - template - void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const - { - // TODO: Revise when linear solvers interface opm-core is done - // Construct linear solver. - // GMRes solver - if ( parameters_.newton_use_gmres_ ) { - Dune::RestartedGMResSolver linsolve(opA, sp, precond, - parameters_.linear_solver_reduction_, - parameters_.linear_solver_restart_, - parameters_.linear_solver_maxiter_, - parameters_.linear_solver_verbosity_); - // Solve system. - linsolve.apply(x, istlb, result); - } - else { // BiCGstab solver - Dune::BiCGSTABSolver linsolve(opA, sp, precond, - parameters_.linear_solver_reduction_, - parameters_.linear_solver_maxiter_, - parameters_.linear_solver_verbosity_); - // Solve system. - linsolve.apply(x, istlb, result); - } - } - void formInterleavedSystem(const std::vector& eqs, Mat& istlA) const { @@ -455,45 +270,8 @@ namespace Opm Vector x(istlA.M()); x = 0.0; - Dune::InverseOperatorResult result; - // Parallel version is deactivated until we figure out how to do it properly. -#if HAVE_MPI - if (parallelInformation_.type() == typeid(ParallelISTLInformation)) - { - typedef Dune::OwnerOverlapCopyCommunication Comm; - const ParallelISTLInformation& info = - boost::any_cast( parallelInformation_); - Comm istlComm(info.communicator()); - // As we use a dune-istl with block size np the number of components - // per parallel is only one. - info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(), - size, 1); - // Construct operator, scalar product and vectors needed. - typedef Dune::OverlappingSchwarzOperator Operator; - Operator opA(istlA, istlComm); - constructPreconditionerAndSolve(opA, x, istlb, istlComm, result); - } - else -#endif - { - // Construct operator, scalar product and vectors needed. - typedef Dune::MatrixAdapter Operator; - Operator opA(istlA); - Dune::Amg::SequentialInformation info; - constructPreconditionerAndSolve(opA, x, istlb, info, result); - } - - // store number of iterations - iterations_ = result.iterations; - - // Check for failure of linear solver. - if (!parameters_.ignoreConvergenceFailure_ && !result.converged) { - const std::string msg("Convergence failure for linear solver."); - if (isIORank_) { - OpmLog::problem(msg); - } - OPM_THROW_NOLOG(LinearSolverProblem, msg); - } + // solve linear system using ISTL methods + istlSolver_.solve( istlA, x, istlb ); // Copy solver output to dx. for (int i = 0; i < size; ++i) { @@ -512,10 +290,7 @@ namespace Opm } protected: - mutable int iterations_; - boost::any parallelInformation_; - bool isIORank_; - + ISTLSolverType istlSolver_; NewtonIterationBlackoilInterleavedParameters parameters_; }; // end NewtonIterationBlackoilInterleavedImpl diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index f59ddfb2a..dcd34b6d5 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -66,7 +66,7 @@ namespace Opm { newton_use_gmres_ = false; linear_solver_reduction_ = 1e-2; - linear_solver_maxiter_ = 75; + linear_solver_maxiter_ = 150; linear_solver_restart_ = 40; linear_solver_verbosity_ = 0; require_full_sparsity_pattern_ = false; diff --git a/opm/autodiff/NonlinearSolver.hpp b/opm/autodiff/NonlinearSolver.hpp index 09af664b9..fb414017b 100644 --- a/opm/autodiff/NonlinearSolver.hpp +++ b/opm/autodiff/NonlinearSolver.hpp @@ -24,6 +24,9 @@ #include #include #include +#include +#include +#include #include namespace Opm { @@ -40,6 +43,10 @@ namespace Opm { typedef ADB::V V; typedef ADB::M M; + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::BlockVector BVector; + // Available relaxation scheme types. enum RelaxType { DAMPEN, SOR }; @@ -131,7 +138,17 @@ namespace Opm { /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas. std::vector computeFluidInPlace(const ReservoirState& x, - const std::vector& fipnum) const; + const std::vector& fipnum) const + { + return model_->computeFluidInPlace(x, fipnum); + } + + std::vector> + computeFluidInPlace(const std::vector& fipnum) const + { + return model_->computeFluidInPlace(fipnum); + } + /// Reference to physical model. const PhysicalModel& model() const; @@ -146,6 +163,8 @@ namespace Opm { /// Apply a stabilization to dx, depending on dxOld and relaxation parameters. void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const; + void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const; + /// The greatest relaxation factor (i.e. smallest factor) allowed. double relaxMax() const { return param_.relax_max_; } diff --git a/opm/autodiff/NonlinearSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp index 8c2cffb9b..41d999591 100644 --- a/opm/autodiff/NonlinearSolver_impl.hpp +++ b/opm/autodiff/NonlinearSolver_impl.hpp @@ -24,6 +24,8 @@ #define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED #include +#include +#include namespace Opm { @@ -99,15 +101,6 @@ namespace Opm return wellIterationsLast_; } - template - std::vector - NonlinearSolver::computeFluidInPlace(const ReservoirState& x, - const std::vector& fipnum) const - { - return model_->computeFluidInPlace(x, fipnum); - } - - template int NonlinearSolver:: @@ -289,6 +282,38 @@ namespace Opm return; } + template + void + NonlinearSolver::stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const + { + // The dxOld is updated with dx. + // If omega is equal to 1., no relaxtion will be appiled. + + BVector tempDxOld = dxOld; + dxOld = dx; + + switch (relaxType()) { + case DAMPEN: + if (omega == 1.) { + return; + } + dx *= omega; + return; + case SOR: + if (omega == 1.) { + return; + } + dx *= omega; + tempDxOld *= (1.-omega); + dx += tempDxOld; + return; + default: + OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type."); + } + + return; + } + } // namespace Opm diff --git a/opm/autodiff/SimFIBODetails.hpp b/opm/autodiff/SimFIBODetails.hpp new file mode 100644 index 000000000..27e985714 --- /dev/null +++ b/opm/autodiff/SimFIBODetails.hpp @@ -0,0 +1,151 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2014-2016 IRIS AS + Copyright 2015 Andreas Lauser + + 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 . +*/ +#ifndef OPM_SIM_FIBO_DETAILS_HPP +#define OPM_SIM_FIBO_DETAILS_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + namespace SimFIBODetails { + typedef std::unordered_map WellMap; + + inline WellMap + mapWells(const std::vector< const Well* >& wells) + { + WellMap wmap; + + for (std::vector< const Well* >::const_iterator + w = wells.begin(), e = wells.end(); + w != e; ++w) + { + wmap.insert(std::make_pair((*w)->name(), *w)); + } + + return wmap; + } + + inline int + resv_control(const WellControls* ctrl) + { + int i, n = well_controls_get_num(ctrl); + + bool match = false; + for (i = 0; (! match) && (i < n); ++i) { + match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE; + } + + if (! match) { i = 0; } + + return i - 1; // -1 if no match, undo final "++" otherwise + } + + inline bool + is_resv(const Wells& wells, + const int w) + { + return (0 <= resv_control(wells.ctrls[w])); + } + + inline bool + is_resv(const WellMap& wmap, + const std::string& name, + const std::size_t step) + { + bool match = false; + + WellMap::const_iterator i = wmap.find(name); + + if (i != wmap.end()) { + const Well* wp = i->second; + + match = (wp->isProducer(step) && + wp->getProductionProperties(step) + .hasProductionControl(WellProducer::RESV)) + || (wp->isInjector(step) && + wp->getInjectionProperties(step) + .hasInjectionControl(WellInjector::RESV)); + } + + return match; + } + + inline std::vector + resvWells(const Wells* wells, + const std::size_t step, + const WellMap& wmap) + { + std::vector resv_wells; + if( wells ) + { + for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) { + if (is_resv(*wells, w) || + ((wells->name[w] != 0) && + is_resv(wmap, wells->name[w], step))) + { + resv_wells.push_back(w); + } + } + } + + return resv_wells; + } + + inline void + historyRates(const PhaseUsage& pu, + const WellProductionProperties& p, + std::vector& rates) + { + assert (! p.predictionMode); + assert (rates.size() == + std::vector::size_type(pu.num_phases)); + + if (pu.phase_used[ BlackoilPhases::Aqua ]) { + const std::vector::size_type + i = pu.phase_pos[ BlackoilPhases::Aqua ]; + + rates[i] = p.WaterRate; + } + + if (pu.phase_used[ BlackoilPhases::Liquid ]) { + const std::vector::size_type + i = pu.phase_pos[ BlackoilPhases::Liquid ]; + + rates[i] = p.OilRate; + } + + if (pu.phase_used[ BlackoilPhases::Vapour ]) { + const std::vector::size_type + i = pu.phase_pos[ BlackoilPhases::Vapour ]; + + rates[i] = p.GasRate; + } + } + } // namespace SimFIBODetails +} // namespace Opm + +#endif // OPM_SIM_FIBO_DETAILS_HPP diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index 7e70ac78f..fa1673e9e 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -27,6 +27,7 @@ #include #include #include +#include namespace Opm { diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index 9832fcda4..ebb82da14 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp @@ -22,6 +22,7 @@ #define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED #include +#include #include namespace Opm { diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp new file mode 100644 index 000000000..353b7196c --- /dev/null +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -0,0 +1,755 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2015 Andreas Lauser + + 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 . +*/ + +#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED + +//#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +namespace Opm { + +class SimulatorFullyImplicitBlackoilEbos; +//class StandardWellsDense; + +/// a simulator for the blackoil model +class SimulatorFullyImplicitBlackoilEbos +{ +public: + typedef typename TTAG(EclFlowProblem) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + + typedef WellStateFullyImplicitBlackoilDense WellState; + typedef BlackoilState ReservoirState; + typedef BlackoilOutputWriterEbos OutputWriter; + typedef BlackoilModelEbos Model; + typedef BlackoilModelParameters ModelParameters; + typedef NonlinearSolver Solver; + typedef StandardWellsDense WellModel; + + + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] geo derived geological properties + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + /// \param[in] has_disgas true for dissolved gas option + /// \param[in] has_vapoil true for vaporized oil option + /// \param[in] eclipse_state the object which represents an internalized ECL deck + /// \param[in] output_writer + /// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow + SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator, + const parameter::ParameterGroup& param, + DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + const bool has_disgas, + const bool has_vapoil, + std::shared_ptr eclipse_state, + BlackoilOutputWriterEbos& output_writer, + const std::vector& threshold_pressures_by_face) + : ebosSimulator_(ebosSimulator), + param_(param), + model_param_(param), + solver_param_(param), + props_(props), + rock_comp_props_(rock_comp_props), + gravity_(gravity), + geo_(geo), + solver_(linsolver), + has_disgas_(has_disgas), + has_vapoil_(has_vapoil), + terminal_output_(param.getDefault("output_terminal", true)), + output_writer_(output_writer), + threshold_pressures_by_face_(threshold_pressures_by_face), + is_parallel_run_( false ) + { + + // Misc init. + const int num_cells = AutoDiffGrid::numCells(grid()); + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + + rateConverter_.reset(new RateConverterType(props_, std::vector(AutoDiffGrid::numCells(grid()), 0))); + +#if HAVE_MPI + if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(solver_.parallelInformation()); + // Only rank 0 does print to std::cout + terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 ); + is_parallel_run_ = ( info.communicator().size() > 1 ); + } +#endif + } + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + ReservoirState& state) + { + WellState prev_well_state; + + if (output_writer_.isRestart()) { + // This is a restart, populate WellState and ReservoirState state objects from restart file + output_writer_.initFromRestartFile(props_.phaseUsage(), props_.permeability(), grid(), state, prev_well_state); + initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_); + } + + // Create timers and file for writing timing info. + Opm::time::StopWatch solver_timer; + double stime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt"; + std::ofstream tstep_os(tstep_filename.c_str()); + + // adaptive time stepping + std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; + if( param_.getDefault("timestep.adaptive", true ) ) + { + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) ); + } + + std::string restorefilename = param_.getDefault("restorefile", std::string("") ); + if( ! restorefilename.empty() ) + { + // -1 means that we'll take the last report step that was written + //const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); + + // output_writer_.restore( timer, + // state, + // prev_well_state, + // restorefilename, + // desiredRestoreStep ); + } + + unsigned int totalLinearizations = 0; + unsigned int totalNonlinearIterations = 0; + unsigned int totalLinearIterations = 0; + bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false ); + std::vector well_potentials; + DynamicListEconLimited dynamic_list_econ_limited; + + bool ooip_computed = false; + std::vector fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData(); + //Get compressed cell fipnum. + std::vector fipnum(Opm::UgGridHelpers::numCells(grid())); + if (fipnum_global.empty()) { + std::fill(fipnum.begin(), fipnum.end(), 0); + } else { + for (size_t c = 0; c < fipnum.size(); ++c) { + fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]]; + } + } + std::vector> OOIP; + + // Main simulation loop. + while (!timer.done()) { + // Report timestep. + step_timer.start(); + if ( terminal_output_ ) + { + std::ostringstream ss; + timer.report(ss); + OpmLog::note(ss.str()); + } + + // Create wells and well state. + WellsManager wells_manager(eclState(), + timer.currentStepNum(), + Opm::UgGridHelpers::numCells(grid()), + Opm::UgGridHelpers::globalCell(grid()), + Opm::UgGridHelpers::cartDims(grid()), + Opm::UgGridHelpers::dimensions(grid()), + Opm::UgGridHelpers::cell2Faces(grid()), + Opm::UgGridHelpers::beginFaceCentroids(grid()), + props_.permeability(), + dynamic_list_econ_limited, + is_parallel_run_, + well_potentials ); + + const Wells* wells = wells_manager.c_wells(); + WellState well_state; + well_state.init(wells, state, prev_well_state); + + // give the polymer and surfactant simulators the chance to do their stuff + handleAdditionalWellInflow(timer, wells_manager, well_state, wells); + + // Compute reservoir volumes for RESV controls. + computeRESV(timer.currentStepNum(), wells, state, well_state); + + // Run a multiple steps of the solver depending on the time step control. + solver_timer.start(); + + const std::vector pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + const WellModel well_model(wells, model_param_, terminal_output_, pv); + + auto solver = createSolver(well_model); + + // write the inital state at the report stage + if (timer.initialStep()) { + + // calculate Intensive Quantities + const auto& gridManager = ebosSimulator_.gridManager(); + const auto& gridView = gridManager.gridView(); + auto elemIt = gridView.template begin<0>(); + auto elemEndIt = gridView.template end<0>(); + ElementContext elemCtx(ebosSimulator_); + for (; elemIt != elemEndIt; ++ elemIt) { + // this is convenient, but slightly inefficient: one only needs to update + // the primary intensive quantities + elemCtx.updateAll(*elemIt); + } + + // No per cell data is written for initial step, but will be + // for subsequent steps, when we have started simulating + output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state ); + } + + // Compute orignal FIP; + if (!ooip_computed) { + OOIP = solver->computeFluidInPlace(fipnum); + FIPUnitConvert(eclState().getUnits(), OOIP); + ooip_computed = true; + } + + if( terminal_output_ ) + { + std::ostringstream step_msg; + boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y"); + step_msg.imbue(std::locale(std::locale::classic(), facet)); + step_msg << "\nTime step " << std::setw(4) <model().beginReportStep(); + + // If sub stepping is enabled allow the solver to sub cycle + // in case the report steps are too large for the solver to converge + // + // \Note: The report steps are met in any case + // \Note: The sub stepping will require a copy of the state variables + if( adaptiveTimeStepping ) { + adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ ); + } + else { + // solve for complete report step + solver->step(timer, state, well_state); + + if( terminal_output_ ) + { + std::ostringstream iter_msg; + iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day); + if (solver->wellIterations() != 0) { + iter_msg << " days well iterations = " << solver->wellIterations() << ", "; + } + iter_msg << "non-linear iterations = " << solver->nonlinearIterations() + << ", total linear iterations = " << solver->linearIterations() + << "\n"; + OpmLog::info(iter_msg.str()); + } + } + + solver->model().endReportStep(); + + // take time that was used to solve system for this reportStep + solver_timer.stop(); + + // accumulate the number of nonlinear and linear Iterations + totalLinearizations += solver->linearizations(); + totalNonlinearIterations += solver->nonlinearIterations(); + totalLinearIterations += solver->linearIterations(); + + // Report timing. + const double st = solver_timer.secsSinceStart(); + + // Compute current FIP. + std::vector> COIP; + COIP = solver->computeFluidInPlace(fipnum); + FIPUnitConvert(eclState().getUnits(), COIP); + std::vector OOIP_totals = FIPTotals(OOIP, state); + std::vector COIP_totals = FIPTotals(COIP, state); + + if ( terminal_output_ ) + { + outputFluidInPlace(OOIP_totals, COIP_totals,eclState().getUnits(), 0); + for (size_t reg = 0; reg < OOIP.size(); ++reg) { + outputFluidInPlace(OOIP[reg], COIP[reg], eclState().getUnits(), reg+1); + } + } + + // accumulate total time + stime += st; + + if ( terminal_output_ ) + { + std::string msg; + msg = "Fully implicit solver took: " + std::to_string(st) + " seconds. Total solver time taken: " + std::to_string(stime) + " seconds."; + OpmLog::note(msg); + } + + if ( output_writer_.output() ) { + SimulatorReport step_report; + step_report.pressure_time = st; + step_report.total_time = step_timer.secsSinceStart(); + step_report.reportParam(tstep_os); + } + + // Increment timer, remember well state. + ++timer; + + // write simulation state at the report stage + output_writer_.writeTimeStep( timer, state, well_state, solver->model() ); + + prev_well_state = well_state; + // The well potentials are only computed if they are needed + // For now thay are only used to determine default guide rates for group controlled wells + if ( is_well_potentials_computed ) { + computeWellPotentials(wells, well_state, well_potentials); + } + + updateListEconLimited(solver, eclState().getSchedule(), timer.currentStepNum(), wells, + well_state, dynamic_list_econ_limited); + } + + // Stop timer and create timing report + total_timer.stop(); + SimulatorReport report; + report.pressure_time = stime; + report.transport_time = 0.0; + report.total_time = total_timer.secsSinceStart(); + report.total_linearizations = totalLinearizations; + report.total_newton_iterations = totalNonlinearIterations; + report.total_linear_iterations = totalLinearIterations; + return report; + } + + const Grid& grid() const + { return ebosSimulator_.gridManager().grid(); } + +protected: + void handleAdditionalWellInflow(SimulatorTimer& timer, + WellsManager& wells_manager, + WellState& well_state, + const Wells* wells) + { } + + std::unique_ptr createSolver(const WellModel& well_model) + { + auto model = std::unique_ptr(new Model(ebosSimulator_, + model_param_, + props_, + geo_, + rock_comp_props_, + well_model, + solver_, + terminal_output_)); + + return std::unique_ptr(new Solver(solver_param_, std::move(model))); + } + + void computeRESV(const std::size_t step, + const Wells* wells, + const BlackoilState& x, + WellState& xw) + { + typedef SimFIBODetails::WellMap WellMap; + + const auto w_ecl = eclState().getSchedule().getWells(step); + const WellMap& wmap = SimFIBODetails::mapWells(w_ecl); + + const std::vector& resv_wells = SimFIBODetails::resvWells(wells, step, wmap); + + const std::size_t number_resv_wells = resv_wells.size(); + std::size_t global_number_resv_wells = number_resv_wells; +#if HAVE_MPI + if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const auto& info = + boost::any_cast(solver_.parallelInformation()); + global_number_resv_wells = info.communicator().sum(global_number_resv_wells); + if ( global_number_resv_wells ) + { + // At least one process has resv wells. Therefore rate converter needs + // to calculate averages over regions that might cross process + // borders. This needs to be done by all processes and therefore + // outside of the next if statement. + rateConverter_->defineState(x, boost::any_cast(solver_.parallelInformation())); + } + } + else +#endif + { + if ( global_number_resv_wells ) + { + rateConverter_->defineState(x); + } + } + + if (! resv_wells.empty()) { + const PhaseUsage& pu = props_.phaseUsage(); + const std::vector::size_type np = props_.numPhases(); + + std::vector distr (np); + std::vector hrates(np); + std::vector prates(np); + + for (std::vector::const_iterator + rp = resv_wells.begin(), e = resv_wells.end(); + rp != e; ++rp) + { + WellControls* ctrl = wells->ctrls[*rp]; + const bool is_producer = wells->type[*rp] == PRODUCER; + + // RESV control mode, all wells + { + const int rctrl = SimFIBODetails::resv_control(ctrl); + + if (0 <= rctrl) { + const std::vector::size_type off = (*rp) * np; + + if (is_producer) { + // Convert to positive rates to avoid issues + // in coefficient calculations. + std::transform(xw.wellRates().begin() + (off + 0*np), + xw.wellRates().begin() + (off + 1*np), + prates.begin(), std::negate()); + } else { + std::copy(xw.wellRates().begin() + (off + 0*np), + xw.wellRates().begin() + (off + 1*np), + prates.begin()); + } + + const int fipreg = 0; // Hack. Ignore FIP regions. + rateConverter_->calcCoeff(prates, fipreg, distr); + + well_controls_iset_distr(ctrl, rctrl, & distr[0]); + } + } + + // RESV control, WCONHIST wells. A bit of duplicate + // work, regrettably. + if (is_producer && wells->name[*rp] != 0) { + WellMap::const_iterator i = wmap.find(wells->name[*rp]); + + if (i != wmap.end()) { + const auto* wp = i->second; + + const WellProductionProperties& p = + wp->getProductionProperties(step); + + if (! p.predictionMode) { + // History matching (WCONHIST/RESV) + SimFIBODetails::historyRates(pu, p, hrates); + + const int fipreg = 0; // Hack. Ignore FIP regions. + rateConverter_->calcCoeff(hrates, fipreg, distr); + + // WCONHIST/RESV target is sum of all + // observed phase rates translated to + // reservoir conditions. Recall sign + // convention: Negative for producers. + const double target = + - std::inner_product(distr.begin(), distr.end(), + hrates.begin(), 0.0); + + well_controls_clear(ctrl); + well_controls_assert_number_of_phases(ctrl, int(np)); + + static const double invalid_alq = -std::numeric_limits::max(); + static const int invalid_vfp = -std::numeric_limits::max(); + + const int ok_resv = + well_controls_add_new(RESERVOIR_RATE, target, + invalid_alq, invalid_vfp, + & distr[0], ctrl); + + // For WCONHIST the BHP limit is set to 1 atm. + // or a value specified using WELTARG + double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm); + const int ok_bhp = + well_controls_add_new(BHP, bhp_limit, + invalid_alq, invalid_vfp, + NULL, ctrl); + + if (ok_resv != 0 && ok_bhp != 0) { + xw.currentControls()[*rp] = 0; + well_controls_set_current(ctrl, 0); + } + } + } + } + } + } + + if( wells ) + { + for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) { + WellControls* ctrl = wells->ctrls[w]; + const bool is_producer = wells->type[w] == PRODUCER; + if (!is_producer && wells->name[w] != 0) { + WellMap::const_iterator i = wmap.find(wells->name[w]); + if (i != wmap.end()) { + const auto* wp = i->second; + const WellInjectionProperties& injector = wp->getInjectionProperties(step); + if (!injector.predictionMode) { + //History matching WCONINJEH + static const double invalid_alq = -std::numeric_limits::max(); + static const int invalid_vfp = -std::numeric_limits::max(); + // For WCONINJEH the BHP limit is set to a large number + // or a value specified using WELTARG + double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits::max(); + const int ok_bhp = + well_controls_add_new(BHP, bhp_limit, + invalid_alq, invalid_vfp, + NULL, ctrl); + if (!ok_bhp) { + OPM_THROW(std::runtime_error, "Failed to add well control."); + } + } + } + } + } + } + } + + + void computeWellPotentials(const Wells* wells, + const WellState& xw, + std::vector& well_potentials) + { + const int nw = wells->number_of_wells; + const int np = wells->number_of_phases; + well_potentials.clear(); + well_potentials.resize(nw*np,0.0); + for (int w = 0; w < nw; ++w) { + for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { + for (int phase = 0; phase < np; ++phase) { + well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase]; + } + } + } + } + + + void updateListEconLimited(const std::unique_ptr& solver, + const Schedule& schedule, + const int current_step, + const Wells* wells, + const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const + { + solver->model().wellModel().updateListEconLimited(schedule, current_step, wells, + well_state, list_econ_limited); + } + + void FIPUnitConvert(const UnitSystem& units, + std::vector>& fip) + { + if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { + for (size_t i = 0; i < fip.size(); ++i) { + fip[i][0] = unit::convert::to(fip[i][0], unit::stb); + fip[i][1] = unit::convert::to(fip[i][1], unit::stb); + fip[i][2] = unit::convert::to(fip[i][2], 1000*unit::cubic(unit::feet)); + fip[i][3] = unit::convert::to(fip[i][3], 1000*unit::cubic(unit::feet)); + fip[i][4] = unit::convert::to(fip[i][4], unit::stb); + fip[i][5] = unit::convert::to(fip[i][5], unit::stb); + fip[i][6] = unit::convert::to(fip[i][6], unit::psia); + } + } + if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { + for (size_t i = 0; i < fip.size(); ++i) { + fip[i][6] = unit::convert::to(fip[i][6], unit::barsa); + } + } + } + + + std::vector FIPTotals(const std::vector>& fip, const ReservoirState& state) + { + std::vector totals(7,0.0); + for (int i = 0; i < 5; ++i) { + for (size_t reg = 0; reg < fip.size(); ++reg) { + totals[i] += fip[reg][i]; + } + } + const int numCells = Opm::AutoDiffGrid::numCells(grid()); + const auto& pv = geo_.poreVolume(); + double pv_hydrocarbon_sum = 0.0; + double p_pv_hydrocarbon_sum = 0.0; + + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + const double& p = fs.pressure(FluidSystem::oilPhaseIdx).value(); + const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value(); + if ( ! is_parallel_run_ ) + { + totals[5] += pv[cellIdx]; + pv_hydrocarbon_sum += pv[cellIdx] * hydrocarbon; + p_pv_hydrocarbon_sum += p * pv[cellIdx] * hydrocarbon; + } + else { + OPM_THROW(std::logic_error, "FIP not yet implemented for MPI"); + } + } + totals[6] = unit::convert::to( (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum), unit::barsa); + return totals; + } + + + + void outputFluidInPlace(const std::vector& oip, const std::vector& cip, const UnitSystem& units, const int reg) + { + std::ostringstream ss; + if (!reg) { + ss << " ===================================================\n" + << " : Field Totals :\n"; + } else { + ss << " ===================================================\n" + << " : FIPNUM report region " + << std::setw(2) << reg << " :\n"; + } + if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) { + ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n" + << std::fixed << std::setprecision(0) + << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n"; + if (!reg) { + ss << " : Pressure is weighted by hydrocarbon pore volume :\n" + << " : Porv volumes are taken at reference conditions :\n"; + } + ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n"; + } + if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { + ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n" + << std::fixed << std::setprecision(0) + << " : PORV =" << std::setw(14) << cip[5] << " RB :\n"; + if (!reg) { + ss << " : Pressure is weighted by hydrocarbon pore voulme :\n" + << " : Pore volumes are taken at reference conditions :\n"; + } + ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n"; + } + ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n" + << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n" + << ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":" + << std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n" + << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n" + << ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":" + << std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n" + << ":========================:==========================================:================:==========================================:\n"; + OpmLog::note(ss.str()); + } + + + const EclipseState& eclState() const + { return *ebosSimulator_.gridManager().eclState(); } + + EclipseState& eclState() + { return *ebosSimulator_.gridManager().eclState(); } + + // Data. + Simulator& ebosSimulator_; + + typedef RateConverter:: + SurfaceToReservoirVoidage< BlackoilPropsAdInterface, + std::vector > RateConverterType; + typedef typename Solver::SolverParameters SolverParameters; + + const parameter::ParameterGroup param_; + ModelParameters model_param_; + SolverParameters solver_param_; + + // Observed objects. + BlackoilPropsAdInterface& props_; + const RockCompressibility* rock_comp_props_; + const double* gravity_; + // Solvers + DerivedGeology& geo_; + NewtonIterationBlackoilInterface& solver_; + // Misc. data + std::vector allcells_; + const bool has_disgas_; + const bool has_vapoil_; + bool terminal_output_; + // output_writer + OutputWriter& output_writer_; + std::unique_ptr rateConverter_; + // Threshold pressures. + std::vector threshold_pressures_by_face_; + // Whether this a parallel simulation or not + bool is_parallel_run_; + +}; + +} // namespace Opm + +#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp new file mode 100644 index 000000000..0c9c32f76 --- /dev/null +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp @@ -0,0 +1,282 @@ +/* + Copyright (c) 2014 SINTEF ICT, Applied Mathematics. + Copyright (c) 2015-2016 IRIS AS + + 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 "config.h" + +#include "SimulatorFullyImplicitBlackoilOutputEbos.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +//For OutputWriterHelper +#include +#include + + +#ifdef HAVE_OPM_GRID +#include +#include +#include +#include +#endif +namespace Opm +{ + namespace detail { + + struct WriterCallEbos : public ThreadHandle :: ObjectInterface + { + BlackoilOutputWriterEbos& writer_; + std::unique_ptr< SimulatorTimerInterface > timer_; + const SimulationDataContainer state_; + const WellState wellState_; + data::Solution simProps_; + const bool substep_; + + explicit WriterCallEbos( BlackoilOutputWriterEbos& writer, + const SimulatorTimerInterface& timer, + const SimulationDataContainer& state, + const WellState& wellState, + const data::Solution& simProps, + bool substep ) + : writer_( writer ), + timer_( timer.clone() ), + state_( state ), + wellState_( wellState ), + simProps_( simProps ), + substep_( substep ) + { + } + + // callback to writer's serial writeTimeStep method + void run () + { + // write data + writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, substep_ ); + } + }; + } + + void + BlackoilOutputWriterEbos:: + writeTimeStepWithoutCellProperties( + const SimulatorTimerInterface& timer, + const SimulationDataContainer& localState, + const WellState& localWellState, + bool substep) + { + data::Solution noCellProperties; + writeTimeStepWithCellProperties(timer, localState, localWellState, noCellProperties, substep); + } + + + + + + void + BlackoilOutputWriterEbos:: + writeTimeStepWithCellProperties( + const SimulatorTimerInterface& timer, + const SimulationDataContainer& localState, + const WellState& localWellState, + const data::Solution& sol, + bool substep) + { + bool isIORank = output_ ; + if( parallelOutput_ && parallelOutput_->isParallel() ) + { + // If this is not the initial write and no substep, then the well + // state used in the computation is actually the one of the last + // step. We need that well state for the gathering. Otherwise + // It an exception with a message like "global state does not + // contain well ..." might be thrown. + int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ? + (timer.reportStepNum() - 1) : timer.reportStepNum(); + // collect all solutions to I/O rank + isIORank = parallelOutput_->collectToIORank( localState, localWellState, sol, wellStateStepNumber ); + } + + const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState; + const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; + + // serial output is only done on I/O rank + if( isIORank ) + { + if( asyncOutput_ ) { + // dispatch the write call to the extra thread + asyncOutput_->dispatch( detail::WriterCallEbos( *this, timer, state, wellState, sol, substep ) ); + } + else { + // just write the data to disk + writeTimeStepSerial( timer, state, wellState, sol, substep ); + } + } + } + + + + void + BlackoilOutputWriterEbos:: + writeTimeStepSerial(const SimulatorTimerInterface& timer, + const SimulationDataContainer& state, + const WellState& wellState, + const data::Solution& sol, + bool substep) + { + // ECL output + if ( eclWriter_ ) + { + const auto& initConfig = eclipseState_.getInitConfig(); + if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) { + std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; + } else { + data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...) + combined_sol.insert(sol.begin(), sol.end()); // ... insert "extra" data (KR, VISC, ...) + eclWriter_->writeTimeStep(timer.reportStepNum(), + substep, + timer.simulationTimeElapsed(), + combined_sol, + wellState.report(phaseUsage_)); + } + } + + // write backup file + if( backupfile_.is_open() ) + { + int reportStep = timer.reportStepNum(); + int currentTimeStep = timer.currentStepNum(); + if( (reportStep == currentTimeStep || // true for SimulatorTimer + currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep + timer.done() ) // true for AdaptiveSimulatorTimer at reportStep + && lastBackupReportStep_ != reportStep ) // only backup report step once + { + // store report step + lastBackupReportStep_ = reportStep; + // write resport step number + backupfile_.write( (const char *) &reportStep, sizeof(int) ); + + try { + backupfile_ << state; + + const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); + backupfile_ << boWellState; + } + catch ( const std::bad_cast& e ) + { + } + + backupfile_ << std::flush; + } + } // end backup + } + + void + BlackoilOutputWriterEbos:: + restore(SimulatorTimerInterface& timer, + BlackoilState& state, + WellStateFullyImplicitBlackoilDense& wellState, + const std::string& filename, + const int desiredResportStep ) + { + std::ifstream restorefile( filename.c_str() ); + if( restorefile ) + { + std::cout << "============================================================================"<. +*/ +#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUTEBOS_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUTEBOS_HEADER_INCLUDED +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + + +#include +#include +#include +#include +#include +#include + +#include + +#ifdef HAVE_OPM_GRID +#include +#endif +namespace Opm +{ + class BlackoilState; + + /** \brief Wrapper class for VTK, Matlab, and ECL output. */ + class BlackoilOutputWriterEbos + { + + public: + // constructor creating different sub writers + template + BlackoilOutputWriterEbos(const Grid& grid, + const parameter::ParameterGroup& param, + const EclipseState& eclipseState, + std::unique_ptr&& eclWriter, + const Opm::PhaseUsage &phaseUsage, + const double* permeability ); + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection with + * visualization tools like ResInsight. This function will extract the + * requested output cell properties specified by the RPTRST keyword + * and write these to file. + */ + template + void writeTimeStep(const SimulatorTimerInterface& timer, + const SimulationDataContainer& reservoirState, + const Opm::WellState& wellState, + const Model& physicalModel, + bool substep = false); + + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection with + * visualization tools like ResInsight. This function will write all + * CellData in simProps to the file as well. + */ + void writeTimeStepWithCellProperties( + const SimulatorTimerInterface& timer, + const SimulationDataContainer& reservoirState, + const Opm::WellState& wellState, + const data::Solution& sol, + bool substep = false); + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection with + * visualization tools like ResInsight. This function will not write + * any cell properties (e.g., those requested by RPTRST keyword) + */ + void writeTimeStepWithoutCellProperties( + const SimulatorTimerInterface& timer, + const SimulationDataContainer& reservoirState, + const Opm::WellState& wellState, + bool substep = false); + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection withS + * visualization tools like ResInsight. This is the function which does + * the actual write to file. + */ + void writeTimeStepSerial(const SimulatorTimerInterface& timer, + const SimulationDataContainer& reservoirState, + const Opm::WellState& wellState, + const data::Solution& simProps, + bool substep); + + /** \brief return output directory */ + const std::string& outputDirectory() const { return outputDir_; } + + /** \brief return true if output is enabled */ + bool output () const { return output_; } + + void restore(SimulatorTimerInterface& timer, + BlackoilState& state, + WellStateFullyImplicitBlackoilDense& wellState, + const std::string& filename, + const int desiredReportStep); + + + template + void initFromRestartFile(const PhaseUsage& phaseusage, + const double* permeability, + const Grid& grid, + SimulationDataContainer& simulatorstate, + WellStateFullyImplicitBlackoilDense& wellstate); + + bool isRestart() const; + + protected: + const bool output_; + std::unique_ptr< ParallelDebugOutputInterface > parallelOutput_; + + // Parameters for output. + const std::string outputDir_; + const int output_interval_; + + int lastBackupReportStep_; + + std::ofstream backupfile_; + Opm::PhaseUsage phaseUsage_; + std::unique_ptr eclWriter_; + const EclipseState& eclipseState_; + + std::unique_ptr< ThreadHandle > asyncOutput_; + }; + + + ////////////////////////////////////////////////////////////// + // + // Implementation + // + ////////////////////////////////////////////////////////////// + template + inline + BlackoilOutputWriterEbos:: + BlackoilOutputWriterEbos(const Grid& grid, + const parameter::ParameterGroup& param, + const Opm::EclipseState& eclipseState, + std::unique_ptr&& eclWriter, + const Opm::PhaseUsage &phaseUsage, + const double* permeability ) + : output_( param.getDefault("output", true) ), + parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, permeability, phaseUsage ) : 0 ), + outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ), + output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ), + lastBackupReportStep_( -1 ), + phaseUsage_( phaseUsage ), + eclipseState_(eclipseState), + asyncOutput_() + { + // For output. + if (output_ && parallelOutput_->isIORank() ) { + eclWriter_ = std::move(eclWriter); + + // Ensure that output dir exists + boost::filesystem::path fpath(outputDir_); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + + // create output thread if enabled and rank is I/O rank + // async output is enabled by default if pthread are enabled +#if HAVE_PTHREAD + const bool asyncOutputDefault = false; +#else + const bool asyncOutputDefault = false; +#endif + if( param.getDefault("async_output", asyncOutputDefault ) ) + { +#if HAVE_PTHREAD + asyncOutput_.reset( new ThreadHandle() ); +#else + OPM_THROW(std::runtime_error,"Pthreads were not found, cannot enable async_output"); +#endif + } + + std::string backupfilename = param.getDefault("backupfile", std::string("") ); + if( ! backupfilename.empty() ) + { + backupfile_.open( backupfilename.c_str() ); + } + } + } + + + template + inline void + BlackoilOutputWriterEbos:: + initFromRestartFile( const PhaseUsage& phaseusage, + const double* permeability, + const Grid& grid, + SimulationDataContainer& simulatorstate, + WellStateFullyImplicitBlackoilDense& wellstate) + { + // gives a dummy dynamic_list_econ_limited + DynamicListEconLimited dummy_list_econ_limited; + WellsManager wellsmanager(eclipseState_, + eclipseState_.getInitConfig().getRestartStep(), + Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::globalCell(grid), + Opm::UgGridHelpers::cartDims(grid), + Opm::UgGridHelpers::dimensions(grid), + Opm::UgGridHelpers::cell2Faces(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + permeability, + dummy_list_econ_limited + // We need to pass the optionaly arguments + // as we get the following error otherwise + // with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11 + // converting to ‘const std::unordered_set >’ from initializer list would use explicit constructo + , false, + std::vector(), + std::unordered_set()); + + const Wells* wells = wellsmanager.c_wells(); + wellstate.resize(wells, simulatorstate); //Resize for restart step + auto restarted = Opm::init_from_restart_file( + eclipseState_, + Opm::UgGridHelpers::numCells(grid) ); + + solutionToSim( restarted.first, phaseusage, simulatorstate ); + wellsToState( restarted.second, phaseusage, wellstate ); + } + + + + + + namespace detail { + template + void getOutputDataEbos(data::Solution& output, + const Opm::PhaseUsage& phaseUsage, + const Model& model, + const RestartConfig& restartConfig, + const int reportStepNum, + const bool log) + { + typedef typename Model::FluidSystem FluidSystem; + + //Get the value of each of the keys + std::map outKeywords = restartConfig.getRestartKeywords(reportStepNum); + for (auto& keyValue : outKeywords) { + keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum); + } + + const auto& ebosModel = model.ebosSimulator().model(); + + //Get shorthands for water, oil, gas + const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua]; + const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid]; + const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour]; + + // extract everything which can possibly be written to disk + int numCells = ebosModel.numGridDof(); + + std::vector pressureOil(numCells); + std::vector temperature(numCells); + + std::vector satWater(numCells); + std::vector satGas(numCells); + + std::vector bWater(numCells); + std::vector bOil(numCells); + std::vector bGas(numCells); + + std::vector rhoWater(numCells); + std::vector rhoOil(numCells); + std::vector rhoGas(numCells); + + std::vector muWater(numCells); + std::vector muOil(numCells); + std::vector muGas(numCells); + + std::vector krWater(numCells); + std::vector krOil(numCells); + std::vector krGas(numCells); + + std::vector Rs(numCells); + std::vector Rv(numCells); + std::vector RsSat(numCells); + std::vector RvSat(numCells); + + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + const auto& intQuants = *ebosModel.cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + pressureOil[cellIdx] = fs.pressure(FluidSystem::oilPhaseIdx).value(); + + temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value(); + + if (aqua_active) { + satWater[cellIdx] = fs.saturation(FluidSystem::waterPhaseIdx).value(); + bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value(); + rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value(); + muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value(); + krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value(); + } + if (vapour_active) { + satGas[cellIdx] = fs.saturation(FluidSystem::gasPhaseIdx).value(); + bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value(); + rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value(); + muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value(); + krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value(); + Rs[cellIdx] = fs.Rs().value(); + Rv[cellIdx] = fs.Rv().value(); + RsSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs, + FluidSystem::oilPhaseIdx, + intQuants.pvtRegionIndex(), + /*maxOilSaturation=*/1.0).value(); + RvSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs, + FluidSystem::gasPhaseIdx, + intQuants.pvtRegionIndex(), + /*maxOilSaturation=*/1.0).value(); + } + + // oil is always active + bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value(); + rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value(); + muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value(); + krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value(); + } + + /** + * Oil Pressures + */ + outKeywords["PRESSURE"] = 0; + output.insert("PRESSURE", + UnitSystem::measure::pressure, + std::move(pressureOil), + data::TargetType::RESTART_SOLUTION); + + /** + * Temperatures + */ + outKeywords["TEMP"] = 0; + output.insert("TEMP", + UnitSystem::measure::temperature, + std::move(temperature), + data::TargetType::RESTART_SOLUTION); + + /** + * Water and gas saturation. + */ + outKeywords["SWAT"] = 0; + outKeywords["SGAS"] = 0; + output.insert("SWAT", + UnitSystem::measure::identity, + std::move(satWater), + data::TargetType::RESTART_SOLUTION); + output.insert("SGAS", + UnitSystem::measure::identity, + std::move(satGas), + data::TargetType::RESTART_SOLUTION); + + /** + * the dissolution factors + */ + outKeywords["RS"] = 0; + outKeywords["RV"] = 0; + output.insert("RS", + UnitSystem::measure::gas_oil_ratio, + std::move(Rs), + data::TargetType::RESTART_SOLUTION); + output.insert("RV", + UnitSystem::measure::oil_gas_ratio, + std::move(Rv), + data::TargetType::RESTART_SOLUTION); + + /** + * Formation volume factors for water, oil, gas + */ + if (aqua_active && outKeywords["BW"] > 0) { + outKeywords["BW"] = 0; + output.insert("BW", + Opm::UnitSystem::measure::water_inverse_formation_volume_factor, + std::move(bWater), + data::TargetType::RESTART_AUXILLARY); + + } + if (liquid_active && outKeywords["BO"] > 0) { + outKeywords["BO"] = 0; + output.insert("BO", + Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, + std::move(bOil), + data::TargetType::RESTART_AUXILLARY); + } + if (vapour_active && outKeywords["BG"] > 0) { + outKeywords["BG"] = 0; + output.insert("BG", + Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, + std::move(bGas), + data::TargetType::RESTART_AUXILLARY); + } + + /** + * Densities for water, oil gas + */ + if (outKeywords["DEN"] > 0) { + outKeywords["DEN"] = 0; + if (aqua_active) { + output.insert("WAT_DEN", + Opm::UnitSystem::measure::density, + std::move(rhoWater), + data::TargetType::RESTART_AUXILLARY); + } + if (liquid_active) { + output.insert("OIL_DEN", + Opm::UnitSystem::measure::density, + std::move(rhoOil), + data::TargetType::RESTART_AUXILLARY); + } + if (vapour_active) { + output.insert("GAS_DEN", + Opm::UnitSystem::measure::density, + std::move(rhoGas), + data::TargetType::RESTART_AUXILLARY); + } + } + + /** + * Viscosities for water, oil gas + */ + if (outKeywords["VISC"] > 0) { + outKeywords["VISC"] = 0; + if (aqua_active) { + output.insert("WAT_VISC", + Opm::UnitSystem::measure::viscosity, + std::move(muWater), + data::TargetType::RESTART_AUXILLARY); + } + if (liquid_active) { + output.insert("OIL_VISC", + Opm::UnitSystem::measure::viscosity, + std::move(muOil), + data::TargetType::RESTART_AUXILLARY); + } + if (vapour_active) { + output.insert("GAS_VISC", + Opm::UnitSystem::measure::viscosity, + std::move(muGas), + data::TargetType::RESTART_AUXILLARY); + } + } + + /** + * Relative permeabilities for water, oil, gas + */ + if (aqua_active && outKeywords["KRW"] > 0) { + outKeywords["KRW"] = 0; + output.insert("WATKR", + Opm::UnitSystem::measure::permeability, + std::move(krWater), + data::TargetType::RESTART_AUXILLARY); + } + if (liquid_active && outKeywords["KRO"] > 0) { + outKeywords["KRO"] = 0; + output.insert("OILKR", + Opm::UnitSystem::measure::permeability, + std::move(krOil), + data::TargetType::RESTART_AUXILLARY); + } + if (vapour_active && outKeywords["KRG"] > 0) { + outKeywords["KRG"] = 0; + output.insert("GASKR", + Opm::UnitSystem::measure::permeability, + std::move(krGas), + data::TargetType::RESTART_AUXILLARY); + } + + /** + * Vaporized and dissolved gas/oil ratio + */ + if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) { + outKeywords["RSSAT"] = 0; + output.insert("RSSAT", + Opm::UnitSystem::measure::gas_oil_ratio, + std::move(RsSat), + data::TargetType::RESTART_AUXILLARY); + } + if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) { + outKeywords["RVSAT"] = 0; + output.insert("RVSAT", + Opm::UnitSystem::measure::oil_gas_ratio, + std::move(RvSat), + data::TargetType::RESTART_AUXILLARY); + } + + + /** + * Bubble point and dew point pressures + */ + if (log && vapour_active && + liquid_active && outKeywords["PBPD"] > 0) { + Opm::OpmLog::warning("Bubble/dew point pressure output unsupported", + "Writing bubble points and dew points (PBPD) to file is unsupported, " + "as the simulator does not use these internally."); + } + + //Warn for any unhandled keyword + if (log) { + for (auto& keyValue : outKeywords) { + if (keyValue.second > 0) { + std::string logstring = "Keyword '"; + logstring.append(keyValue.first); + logstring.append("' is unhandled for output to file."); + Opm::OpmLog::warning("Unhandled output keyword", logstring); + } + } + } + + } + + /** + * Checks if the summaryConfig has a keyword with the standardized field, region, or block prefixes. + */ + inline bool hasFRBKeyword(const SummaryConfig& summaryConfig, const std::string keyword) { + std::string field_kw = "F" + keyword; + std::string region_kw = "R" + keyword; + std::string block_kw = "B" + keyword; + return summaryConfig.hasKeyword(field_kw) + || summaryConfig.hasKeyword(region_kw) + || summaryConfig.hasKeyword(block_kw); + } + + + /** + * Returns the data as asked for in the summaryConfig + */ + template + void getSummaryData(data::Solution& output, + const Opm::PhaseUsage& phaseUsage, + const Model& physicalModel, + const SummaryConfig& summaryConfig) { + + const typename Model::FIPData& fip = physicalModel.getFIPData(); + + //Get shorthands for water, oil, gas + const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua]; + const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid]; + const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour]; + + /** + * Now process all of the summary config files + */ + // Water in place + if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) { + output.insert("WIP", + Opm::UnitSystem::measure::volume, + fip.fip[Model::FIPData::FIP_AQUA], + data::TargetType::SUMMARY ); + } + if (liquid_active) { + const std::vector& oipl = fip.fip[Model::FIPData::FIP_LIQUID]; + const int size = oipl.size(); + + const std::vector& oipg = vapour_active ? fip.fip[Model::FIPData::FIP_VAPORIZED_OIL] : std::vector(size,0.0); + std::vector oip = oipl; + if (vapour_active) { + oip.insert(oip.end(), oipg.begin(), oipg.end()); + } + + //Oil in place (liquid phase only) + if (hasFRBKeyword(summaryConfig, "OIPL")) { + output.insert("OIPL", + Opm::UnitSystem::measure::volume, + oipl, + data::TargetType::SUMMARY ); + } + + //Oil in place (gas phase only) + if (hasFRBKeyword(summaryConfig, "OIPG")) { + output.insert("OIPG", + Opm::UnitSystem::measure::volume, + oipg, + data::TargetType::SUMMARY ); + } + // Oil in place (in liquid and gas phases) + if (hasFRBKeyword(summaryConfig, "OIP")) { + output.insert("OIP", + Opm::UnitSystem::measure::volume, + oip, + data::TargetType::SUMMARY ); + } + + + + } + if (vapour_active) { + const std::vector& gipg = fip.fip[Model::FIPData::FIP_VAPOUR]; + const int size = gipg.size(); + + const std::vector& gipl= liquid_active ? fip.fip[Model::FIPData::FIP_DISSOLVED_GAS] : std::vector(size,0.0); + std::vector gip = gipg; + if (liquid_active) { + gip.insert(gip.end(), gipl.begin(), gipl.end()); + } + + // Gas in place (gas phase only) + if (hasFRBKeyword(summaryConfig, "GIPG")) { + output.insert("GIPG", + Opm::UnitSystem::measure::volume, + gipg, + data::TargetType::SUMMARY ); + } + // Gas in place (liquid phase only) + if (hasFRBKeyword(summaryConfig, "GIPL")) { + output.insert("GIPL", + Opm::UnitSystem::measure::volume, + gipl, + data::TargetType::SUMMARY ); + } + // Gas in place (in both liquid and gas phases) + if (hasFRBKeyword(summaryConfig, "GIP")) { + output.insert("GIP", + Opm::UnitSystem::measure::volume, + gip, + data::TargetType::SUMMARY ); + } + } + // Cell pore volume in reservoir conditions + if (hasFRBKeyword(summaryConfig, "RPV")) { + output.insert("RPV", + Opm::UnitSystem::measure::volume, + fip.fip[Model::FIPData::FIP_PV], + data::TargetType::SUMMARY ); + } + // Pressure averaged value (hydrocarbon pore volume weighted) + if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) { + output.insert("PRH", + Opm::UnitSystem::measure::pressure, + fip.fip[Model::FIPData::FIP_WEIGHTED_PRESSURE], + data::TargetType::SUMMARY ); + } + } + + } + + + + + template + inline void + BlackoilOutputWriterEbos:: + writeTimeStep(const SimulatorTimerInterface& timer, + const SimulationDataContainer& localState, + const WellState& localWellState, + const Model& physicalModel, + bool substep) + { + data::Solution cellData{}; + const RestartConfig& restartConfig = eclipseState_.getRestartConfig(); + const SummaryConfig& summaryConfig = eclipseState_.getSummaryConfig(); + const int reportStepNum = timer.reportStepNum(); + bool logMessages = output_ && parallelOutput_->isIORank(); + + if( output_ && !parallelOutput_->isParallel() ) + { + + detail::getOutputDataEbos( cellData, phaseUsage_, physicalModel, + restartConfig, reportStepNum, logMessages ); + detail::getSummaryData( cellData, phaseUsage_, physicalModel, summaryConfig ); + } + else + { + if ( logMessages ) + { + std::map rstKeywords = restartConfig.getRestartKeywords(reportStepNum); + std::vector keywords = + { "WIP", "OIPL", "OIPG", "OIP", "GIPG", "GIPL", "GIP", + "RPV", "FRPH", "RPRH"}; + + std::ostringstream str; + str << "Output of restart/summary config not supported in parallel. Requested keywords were "; + std::size_t no_kw = 0; + + auto func = [&] (const char* kw) + { + if ( detail::hasFRBKeyword(summaryConfig, kw) ) + { + str << kw << " "; + ++ no_kw; + } + }; + + std::for_each(keywords.begin(), keywords.end(), func); + + for (auto& keyValue : rstKeywords) + { + str << keyValue.first << " "; + ++ no_kw; + } + + if ( no_kw ) + { + Opm::OpmLog::warning("Unhandled ouput request", str.str()); + } + } + } + + writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep); + } +} +#endif diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp new file mode 100644 index 000000000..ea5161d90 --- /dev/null +++ b/opm/autodiff/StandardWellsDense.hpp @@ -0,0 +1,1752 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + Copyright 2016 IRIS AS + + 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 . +*/ + + +#ifndef OPM_STANDARDWELLSDENSE_HEADER_INCLUDED +#define OPM_STANDARDWELLSDENSE_HEADER_INCLUDED + +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Opm { + + + /// Class for handling the standard well model. + template + class StandardWellsDense { + public: + + // --------- Types --------- + typedef DenseAd::Evaluation EvalWell; + typedef WellStateFullyImplicitBlackoilDense WellState; + typedef BlackoilModelParameters ModelParameters; + + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; + + // --------- Public methods --------- + StandardWellsDense(const Wells* wells_arg, + const ModelParameters& param, + const bool terminal_output, + const std::vector& pv) + : wells_active_(wells_arg!=nullptr) + , wells_(wells_arg) + , param_(param) + , terminal_output_(terminal_output) + , fluid_(nullptr) + , active_(nullptr) + , vfp_properties_(nullptr) + , well_perforation_densities_(wells_arg->well_connpos[wells_arg->number_of_wells]) + , well_perforation_pressure_diffs_(wells_arg->well_connpos[wells_arg->number_of_wells]) + , wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases) + , F0_(wells_arg->number_of_wells * wells_arg->number_of_phases) + { + invDuneD_.setBuildMode( Mat::row_wise ); + duneC_.setBuildMode( Mat::row_wise ); + duneB_.setBuildMode( Mat::row_wise ); + } + + void init(const BlackoilPropsAdInterface* fluid_arg, + const std::vector* active_arg, + const VFPProperties* vfp_properties_arg, + const double gravity_arg, + const std::vector& depth_arg, + const std::vector& pv_arg) + { + + if ( ! localWellsActive() ) { + return; + } + + fluid_ = fluid_arg; + active_ = active_arg; + vfp_properties_ = vfp_properties_arg; + gravity_ = gravity_arg; + cell_depths_ = extractPerfData(depth_arg); + pv_ = pv_arg; + + // setup sparsity pattern for the matrices + //[A B^T [x = [ res + // C D] x_well] res_well] + + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const int nc = numCells(); + + // set invDuneD + invDuneD_.setSize( nw, nw, nw ); + + // set duneC + duneC_.setSize( nw, nc, nperf ); + + // set duneB + duneB_.setSize( nw, nc, nperf ); + + for(auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + row.insert(row.index()); + } + + for(auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) { + // Add nonzeros for diagonal + for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + row.insert(cell_idx); + } + } + + // make the B^T matrix + for(auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) { + for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + row.insert(cell_idx); + } + } + + resWell_.resize( nw ); + + // resize temporary class variables + Cx_.resize( duneC_.N() ); + invDrw_.resize( invDuneD_.N() ); + } + + + template + IterationReport assemble(Simulator& ebosSimulator, + const int iterationIdx, + const double dt, + WellState& well_state) { + + IterationReport iter_report = {false, false, 0, 0}; + if ( ! localWellsActive() ) { + return iter_report; + } + + resetWellControlFromState(well_state); + updateWellControls(well_state); + // Set the primary variables for the wells + setWellVariables(well_state); + + if (iterationIdx == 0) { + computeWellConnectionPressures(ebosSimulator, well_state); + computeAccumWells(); + } + + if (param_.solve_welleq_initially_ && iterationIdx == 0) { + // solve the well equations as a pre-processing step + iter_report = solveWellEq(ebosSimulator, dt, well_state); + } + assembleWellEq(ebosSimulator, dt, well_state, false); + + + if (param_.compute_well_potentials_) { + //wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state); + } + return iter_report; + } + + template + void assembleWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + + // clear all entries + duneB_ = 0.0; + duneC_ = 0.0; + invDuneD_ = 0.0; + resWell_ = 0.0; + + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + const double volume = 0.002831684659200; // 0.1 cu ft; + for (int w = 0; w < nw; ++w) { + bool allow_cf = allow_cross_flow(w, ebosSimulator); + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + + const int cell_idx = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + std::vector cq_s(np,0.0); + computeWellFlux(w, wells().WI[perf], intQuants, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + + for (int p1 = 0; p1 < np; ++p1) { + + if (!only_wells) { + // subtract sum of phase fluxes in the reservoir equation. + ebosResid[cell_idx][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value(); + } + + // subtract sum of phase fluxes in the well equations. + resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value(); + + // assemble the jacobians + for (int p2 = 0; p2 < np; ++p2) { + if (!only_wells) { + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2); + duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].derivative(p2+3); // intput in transformed matrix + duneC_[w][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2); + } + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+3); + } + + // Store the perforation phase flux for later usage. + well_state.perfPhaseRates()[perf*np + p1] = cq_s[p1].value(); + } + // Store the perforation pressure for later usage. + well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; + } + + // add vol * dF/dt + Q to the well equations; + for (int p1 = 0; p1 < np; ++p1) { + EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; + resWell_loc += getQs(w, p1); + for (int p2 = 0; p2 < np; ++p2) { + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+3); + } + resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value(); + } + } + + // do the local inversion of D. + localInvert( invDuneD_ ); + + } + + template + bool allow_cross_flow(const int w, Simulator& ebosSimulator) const { + + if (wells().allow_cf[w]) { + return true; + } + // check for special case where all perforations have cross flow + // then the wells must allow for cross flow + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + EvalWell bhp = getBhp(w); + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + wellPerforationPressureDiffs()[perf]; + EvalWell drawdown = pressure - well_pressure; + + if ( drawdown.value() < 0 && wells().type[w] == INJECTOR) { + return false; + } + if ( drawdown.value() > 0 && wells().type[w] == PRODUCER) { + return false; + } + } + return true; + } + + void localInvert(Mat& istlA) const { + for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { + for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { + //std::cout << (*col) << std::endl; + (*col).invert(); + } + } + } + + void print(Mat& istlA) const { + for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { + for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { + std::cout << row.index() << " " << col.index() << "/n \n"<<(*col) << std::endl; + } + } + } + + // substract Binv(D)rw from r; + void apply( BVector& r) const { + if ( ! localWellsActive() ) { + return; + } + + assert( invDrw_.size() == invDuneD_.N() ); + + invDuneD_.mv(resWell_,invDrw_); + duneB_.mmtv(invDrw_, r); + } + + // subtract B*inv(D)*C * x from A*x + void apply(const BVector& x, BVector& Ax) { + if ( ! localWellsActive() ) { + return; + } + assert( Cx_.size() == duneC_.N() ); + + BVector& invDCx = invDrw_; + assert( invDCx.size() == invDuneD_.N()); + + duneC_.mv(x, Cx_); + invDuneD_.mv(Cx_, invDCx); + duneB_.mmtv(invDCx,Ax); + } + + // apply well model with scaling of alpha + void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) + { + if ( ! localWellsActive() ) { + return; + } + + if( scaleAddRes_.size() != Ax.size() ) { + scaleAddRes_.resize( Ax.size() ); + } + + scaleAddRes_ = 0.0; + apply( x, scaleAddRes_ ); + Ax.axpy( alpha, scaleAddRes_ ); + } + + // xw = inv(D)*(rw - C*x) + void recoverVariable(const BVector& x, BVector& xw) const { + if ( ! localWellsActive() ) { + return; + } + BVector resWell = resWell_; + duneC_.mmv(x, resWell); + invDuneD_.mv(resWell, xw); + } + + + int flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + return phaseToComp[ phaseIdx ]; + } + + int flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + return flowToEbos[ flowPv ]; + } + + int ebosCompToFlowPhaseIdx( const int compIdx ) const + { + const int compToPhase[ 3 ] = { Oil, Water, Gas }; + return compToPhase[ compIdx ]; + } + + + + std::vector + extractPerfData(const std::vector& in) const { + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + std::vector out(nperf); + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + const int well_idx = wells().well_cells[perf]; + out[perf] = in[well_idx]; + } + } + return out; + } + + int numPhases() const { return wells().number_of_phases; } + + int numCells() const { return pv_.size();} + + template + void resetWellControlFromState(WellState xw) { + const int nw = wells_->number_of_wells; + for (int w = 0; w < nw; ++w) { + WellControls* wc = wells_->ctrls[w]; + well_controls_set_current( wc, xw.currentControls()[w]); + } + } + + const Wells& wells() const + { + assert(wells_ != 0); + return *(wells_); + } + + const Wells* wellsPointer() const + { + return wells_; + } + + /// return true if wells are available in the reservoir + bool wellsActive() const + { + return wells_active_; + } + void setWellsActive(const bool wells_active) + { + wells_active_ = wells_active; + } + /// return true if wells are available on this process + bool localWellsActive() const + { + return wells_ ? (wells_->number_of_wells > 0 ) : false; + } + + int numWellVars() const + { + if ( !localWellsActive() ) + { + return 0; + } + + // For each well, we have a bhp variable, and one flux per phase. + const int nw = wells().number_of_wells; + return (numPhases() + 1) * nw; + } + + /// Density of each well perforation + std::vector wellPerforationDensities() const { + return well_perforation_densities_; + } + + /// Diff to bhp for each well perforation. + std::vector wellPerforationPressureDiffs() const + { + return well_perforation_pressure_diffs_; + } + + + typedef DenseAd::Evaluation Eval; + EvalWell extendEval(Eval in) const { + EvalWell out = 0.0; + out.setValue(in.value()); + for(int i = 0;i<3;++i) { + out.setDerivative(i, in.derivative(flowToEbosPvIdx(i))); + } + return out; + } + + void + setWellVariables(const WellState& xw) { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int w = 0; w < nw; ++w) { + wellVariables_[w + nw*phaseIdx] = 0.0; + wellVariables_[w + nw*phaseIdx].setValue(xw.wellSolutions()[w + nw* phaseIdx]); + wellVariables_[w + nw*phaseIdx].setDerivative(np + phaseIdx, 1.0); + } + } + } + + void print(EvalWell in) const { + std::cout << in.value() << std::endl; + for (int i = 0; i < in.size; ++i) { + std::cout << in.derivative(i) << std::endl; + } + } + + void + computeAccumWells() { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + for (int w = 0; w < nw; ++w) { + F0_[w + nw * phaseIdx] = wellVolumeFraction(w,phaseIdx).value(); + } + } + } + + + + template + void + computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const double& cdp, const bool& allow_cf, std::vector& cq_s) const + { + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + EvalWell bhp = getBhp(w); + const int np = wells().number_of_phases; + std::vector cmix_s(np,0.0); + for (int phase = 0; phase < np; ++phase) { + //int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + cmix_s[phase] = wellVolumeFraction(w,phase); + } + const auto& fs = intQuants.fluidState(); + EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + EvalWell rs = extendEval(fs.Rs()); + EvalWell rv = extendEval(fs.Rv()); + std::vector b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); + mob_perfcells_dense[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + } + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + cdp; + EvalWell drawdown = pressure - well_pressure; + + // injection perforations + if ( drawdown.value() > 0 ) { + + //Do nothing if crossflow is not allowed + if (!allow_cf && wells().type[w] == INJECTOR) + return; + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + const EvalWell cq_p = - Tw * (mob_perfcells_dense[phase] * drawdown); + cq_ps[phase] = b_perfcells_dense[phase] * cq_p; + } + + if ((*active_)[Oil] && (*active_)[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const EvalWell cq_psOil = cq_ps[oilpos]; + const EvalWell cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs * cq_psOil; + cq_ps[oilpos] += rv * cq_psGas; + } + + // map to ADB + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase]; + } + + } else { + //Do nothing if crossflow is not allowed + if (!allow_cf && wells().type[w] == PRODUCER) + return; + + // Using total mobilities + EvalWell total_mob_dense = mob_perfcells_dense[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob_dense += mob_perfcells_dense[phase]; + } + // injection perforations total volume rates + const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown); + + // compute volume ratio between connection at standard conditions + EvalWell volumeRatio = 0.0; + if ((*active_)[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + } + + if ((*active_)[Oil] && (*active_)[Gas]) { + EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx)); + EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + EvalWell rvPerf = 0.0; + if (cmix_s[gaspos] > 0) + rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; + + if (rvPerf.value() > rvSatEval.value()) { + rvPerf = rvSatEval; + //rvPerf.setValue(rvSatEval.value()); + } + + EvalWell rsPerf = 0.0; + if (cmix_s[oilpos] > 0) + rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; + + if (rsPerf.value() > rsSatEval.value()) { + //rsPerf = 0.0; + rsPerf= rsSatEval; + } + + // Incorporate RS/RV factors if both oil and gas active + const EvalWell d = 1.0 - rvPerf * rsPerf; + + const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + //std::cout << "tmp_oil " < + IterationReport solveWellEq(Simulator& ebosSimulator, + const double dt, + WellState& well_state) + { + const int nw = wells().number_of_wells; + WellState well_state0 = well_state; + + int it = 0; + bool converged; + do { + assembleWellEq(ebosSimulator, dt, well_state, true); + converged = getWellConvergence(ebosSimulator, it); + + if (converged) { + break; + } + + ++it; + if( localWellsActive() ) + { + BVector dx_well (nw); + invDuneD_.mv(resWell_, dx_well); + + updateWellState(dx_well, well_state); + updateWellControls(well_state); + setWellVariables(well_state); + } + } while (it < 15); + + if (!converged) { + well_state = well_state0; + } + + const bool failed = false; // Not needed in this method. + const int linear_iters = 0; // Not needed in this method + return IterationReport{failed, converged, linear_iters, it}; + } + + void printIf(int c, double x, double y, double eps, std::string type) { + if (std::abs(x-y) > eps) { + std::cout << type << " " << c << ": "< residual() { + + const int np = numPhases(); + const int nw = wells().number_of_wells; + std::vector res(np*nw); + for( int p=0; p + bool getWellConvergence(Simulator& ebosSimulator, + const int iteration) + { + typedef std::vector< double > Vector; + const int np = numPhases(); + const int nc = numCells(); + const double tol_wells = param_.tolerance_wells_; + const double maxResidualAllowed = param_.max_residual_allowed_; + + Vector R_sum(np); + Vector B_avg(np); + Vector maxCoeff(np); + Vector maxNormWell(np); + + std::vector< Vector > B( np, Vector( nc ) ); + std::vector< Vector > R2( np, Vector( nc ) ); + std::vector< Vector > tempV( np, Vector( nc ) ); + + for ( int idx = 0; idx < np; ++idx ) + { + Vector& B_idx = B[ idx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + B_idx [cell_idx] = 1 / fs.invB(ebosPhaseIdx).value(); + } + } + + detail::convergenceReduction(B, tempV, R2, + R_sum, maxCoeff, B_avg, maxNormWell, + nc, np, pv_, residual() ); + + + Vector well_flux_residual(np); + + bool converged_Well = true; + // Finish computation + for ( int idx = 0; idx < np; ++idx ) + { + well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx]; + converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); + } + + // if one of the residuals is NaN, throw exception, so that the solver can be restarted + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + + if (std::isnan(well_flux_residual[phaseIdx])) { + OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName); + } + if (well_flux_residual[phaseIdx] > maxResidualAllowed) { + OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName); + } + } + + if ( terminal_output_ ) + { + // Only rank 0 does print to std::cout + if (iteration == 0) { + std::string msg; + msg = "Iter"; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); + msg += " W-FLUX(" + phaseName + ")"; + } + OpmLog::note(msg); + } + std::ostringstream ss; + const std::streamsize oprec = ss.precision(3); + const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); + ss << std::setw(4) << iteration; + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + ss << std::setw(11) << well_flux_residual[phaseIdx]; + } + ss.precision(oprec); + ss.flags(oflags); + OpmLog::note(ss.str()); + } + return converged_Well; + } + + + + template + void + computeWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw) + { + if( ! localWellsActive() ) return ; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(ebosSimulator, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, cell_depths_, gravity_); + + } + + template + void + computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + const PhaseUsage& pu = fluid_->phaseUsage(); + const int np = fluid_->numPhases(); + b_perf.resize(nperf*np); + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + surf_dens_perf.resize(nperf*np); + + // Compute the average pressure in each well block + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + + const int cell_idx = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + const double p_above = perf == wells().well_connpos[w] ? xw.bhp()[w] : xw.perfPress()[perf - 1]; + const double p_avg = (xw.perfPress()[perf] + p_above)/2; + double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); + + if (pu.phase_used[BlackoilPhases::Aqua]) { + b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * pu.num_phases] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; + int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + if (oilrate > 0) { + const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + double rv = 0.0; + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; + int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + if (pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + if (gasrate > 0) { + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); + double rs = 0.0; + if (oilrate > 0) { + rs = gasrate / oilrate; + } + rs = std::min(rs, rsmax_perf[perf]); + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + // Surface density. + for (int p = 0; p < pu.num_phases; ++p) { + surf_dens_perf[np*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + } + } + } + } + + template + void updateWellState(const BVector& dwells, + WellState& well_state) + { + if( localWellsActive() ) + { + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + double dFLimit = 0.2; + double dBHPLimit = 2.0; + std::vector xvar_well_old = well_state.wellSolutions(); + + for (int w = 0; w < nw; ++w) { + + // update the second and third well variable (The flux fractions) + std::vector F(np,0.0); + const int sign2 = dwells[w][flowPhaseToEbosCompIdx(1)] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(1)]),dFLimit); + well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited; + const int sign3 = dwells[w][flowPhaseToEbosCompIdx(2)] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(2)]),dFLimit); + well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited; + F[Water] = well_state.wellSolutions()[nw + w]; + F[Gas] = well_state.wellSolutions()[2*nw + w]; + F[Oil] = 1.0 - F[Water] - F[Gas]; + + if (F[Water] < 0.0) { + F[Gas] /= (1.0 - F[Water]); + F[Oil] /= (1.0 - F[Water]); + F[Water] = 0.0; + } + if (F[Gas] < 0.0) { + F[Water] /= (1.0 - F[Gas]); + F[Oil] /= (1.0 - F[Gas]); + F[Gas] = 0.0; + } + if (F[Oil] < 0.0) { + F[Water] /= (1.0 - F[Oil]); + F[Gas] /= (1.0 - F[Oil]); + F[Oil] = 0.0; + } + well_state.wellSolutions()[nw + w] = F[Water]; + well_state.wellSolutions()[2*nw + w] = F[Gas]; + //std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + + // The interpretation of the first well variable depends on the well control + const 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 = well_state.currentControls()[w]; + const double target_rate = well_controls_iget_target(wc, current); + + std::vector g = {1,1,0.01}; + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + const double* distr = well_controls_iget_distr(wc, current); + for (int p = 0; p < np; ++p) { + F[p] /= distr[p]; + } + } else { + for (int p = 0; p < np; ++p) { + F[p] /= g[p]; + } + } + + switch (well_controls_iget_type(wc, current)) { + case THP: // The BHP and THP both uses the total rate as first well variable. + case BHP: + { + well_state.wellSolutions()[w] = xvar_well_old[w] - dwells[w][flowPhaseToEbosCompIdx(0)]; + + switch (wells().type[w]) { + case INJECTOR: + for (int p = 0; p < np; ++p) { + const double comp_frac = wells().comp_frac[np*w + p]; + well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w]; + } + break; + case PRODUCER: + for (int p = 0; p < np; ++p) { + well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p]; + } + break; + } + + if (well_controls_iget_type(wc, current) == THP) { + + // Calculate bhp from thp control and well rates + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + // pick the density in the top layer + const int perf = wells().well_connpos[w]; + const double rho = well_perforation_densities_[perf]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + well_state.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + well_state.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + } + + } + break; + case SURFACE_RATE: // Both rate controls use bhp as first well variable + case RESERVOIR_RATE: + { + const int sign1 = dwells[w][flowPhaseToEbosCompIdx(0)] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(0)]),std::abs(xvar_well_old[w])*dBHPLimit); + well_state.wellSolutions()[w] = std::max(xvar_well_old[w] - dx1_limited,1e5); + well_state.bhp()[w] = well_state.wellSolutions()[w]; + + if (well_controls_iget_type(wc, current) == SURFACE_RATE) { + if (wells().type[w]==PRODUCER) { + + double F_target = 0.0; + for (int p = 0; p < np; ++p) { + F_target += wells().comp_frac[np*w + p] * F[p]; + } + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np*w + p] = F[p] * target_rate /F_target; + } + } else { + + for (int p = 0; p < np; ++p) { + well_state.wellRates()[w*np + p] = wells().comp_frac[np*w + p] * target_rate; + } + } + } else { // RESERVOIR_RATE + for (int p = 0; p < np; ++p) { + well_state.wellRates()[np*w + p] = F[p] * target_rate; + } + } + } + break; + } + } + } + } + + + + template + void updateWellControls(WellState& xw) + { + if( !localWellsActive() ) return ; + + std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + #pragma omp parallel for schedule(dynamic) + for (int w = 0; w < nw; ++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. + int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (wellhelpers::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + w, np, wells().type[w], wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + // We disregard terminal_ouput here as with it only messages + // for wells on one process will be printed. + std::ostringstream ss; + ss << "Switching control mode for well " << wells().name[w] + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + OpmLog::info(ss.str()); + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + well_controls_set_current( wc, current); + + + + // Updating well state and primary variables if constraint is broken + + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + break; + + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if ((*active_)[ Oil ]) { + liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if ((*active_)[ Gas ]) { + vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + // pick the density in the top layer + const int perf = wells().well_connpos[w]; + const double rho = well_perforation_densities_[perf]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + rho, gravity_); + + xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + break; + } + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; + + case SURFACE_RATE: + // assign target value as initial guess for injectors and + // single phase producers (orat, grat, wrat) + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * w + phase]; + //if (compi > 0.0) { + xw.wellRates()[np*w + phase] = target * compi; + //} + } + } else if (well_type == PRODUCER) { + + // only set target as initial rates for single phase + // producers. (orat, grat and wrat, and not lrat) + // lrat will result in numPhasesWithTargetsUnderThisControl == 2 + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + + break; + } + std::vector g = {1,1,0.01}; + if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { + const double* distr = well_controls_iget_distr(wc, current); + for (int phase = 0; phase < np; ++phase) { + g[phase] = distr[phase]; + } + } + switch (well_controls_iget_type(wc, current)) { + case THP: + case BHP: + { + const WellType& well_type = wells().type[w]; + xw.wellSolutions()[w] = 0.0; + if (well_type == INJECTOR) { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; + } + } else { + for (int p = 0; p < np; ++p) { + xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; + } + } + } + break; + + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + xw.wellSolutions()[w] = xw.bhp()[w]; + } + break; + } + + double tot_well_rate = 0.0; + for (int p = 0; p < np; ++p) { + tot_well_rate += g[p] * xw.wellRates()[np*w + p]; + } + if(std::abs(tot_well_rate) > 0) { + xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; + xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; + } else { + xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; + xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; + } + } + } + } + + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + return flowToEbos[ phaseIdx ]; + } + + /// upate the dynamic lists related to economic limits + template + void + updateListEconLimited(const Schedule& schedule, + const int current_step, + const Wells* wells_struct, + const WellState& well_state, + DynamicListEconLimited& list_econ_limited) const + { + const int nw = wells_struct->number_of_wells; + + for (int w = 0; w < nw; ++w) { + // flag to check if the mim oil/gas rate limit is violated + bool rate_limit_violated = false; + const std::string& well_name = wells_struct->name[w]; + const Well* well_ecl = schedule.getWell(well_name); + const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step); + + // economic limits only apply for production wells. + if (wells_struct->type[w] != PRODUCER) { + continue; + } + + // if no limit is effective here, then continue to the next well + if ( !econ_production_limits.onAnyEffectiveLimit() ) { + continue; + } + // for the moment, we only handle rate limits, not handling potential limits + // the potential limits should not be difficult to add + const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit(); + if (quantity_limit == WellEcon::POTN) { + const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n") + + std::string("All the limits will be evaluated based on RATE. "); + OpmLog::warning("NOT_SUPPORTING_POTN", msg); + } + + const WellMapType& well_map = well_state.wellMap(); + const typename WellMapType::const_iterator i_well = well_map.find(well_name); + assert(i_well != well_map.end()); // should always be found? + const WellMapEntryType& map_entry = i_well->second; + const int well_number = map_entry[0]; + + if (econ_production_limits.onAnyRateLimit()) { + rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number); + } + + if (rate_limit_violated) { + if (econ_production_limits.endRun()) { + const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n") + + std::string("the program will keep running after ") + well_name + std::string(" is closed"); + OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message); + } + + if (econ_production_limits.validFollowonWell()) { + OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet"); + } + + if (well_ecl->getAutomaticShutIn()) { + list_econ_limited.addShutWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit"); + OpmLog::info(msg); + } else { + list_econ_limited.addStoppedWell(well_name); + const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit"); + OpmLog::info(msg); + } + // the well is closed, not need to check other limits + continue; + } + + // checking for ratio related limits, mostly all kinds of ratio. + bool ratio_limits_violated = false; + RatioCheckTuple ratio_check_return; + + if (econ_production_limits.onAnyRatioLimit()) { + ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry); + ratio_limits_violated = std::get<0>(ratio_check_return); + } + + if (ratio_limits_violated) { + const bool last_connection = std::get<1>(ratio_check_return); + const int worst_offending_connection = std::get<2>(ratio_check_return); + + const int perf_start = map_entry[1]; + + assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2])); + + const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection]; + list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection); + const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ") + + well_name + std::string(" will be closed due to economic limit"); + OpmLog::info(msg); + + if (last_connection) { + list_econ_limited.addShutWell(well_name); + const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed"); + OpmLog::info(msg2); + } + } + + } + } + + template + void computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav) { + // Compute densities + well_perforation_densities_ = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid_->phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // Compute pressure deltas + well_perforation_pressure_diffs_ = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), depth_perf, well_perforation_densities_, grav); + + + + } + + protected: + bool wells_active_; + const Wells* wells_; + ModelParameters param_; + bool terminal_output_; + + const BlackoilPropsAdInterface* fluid_; + const std::vector* active_; + const VFPProperties* vfp_properties_; + double gravity_; + // the depth of the all the cell centers + // for standard Wells, it the same with the perforation depth + std::vector cell_depths_; + std::vector pv_; + + std::vector well_perforation_densities_; + std::vector well_perforation_pressure_diffs_; + + std::vector wellVariables_; + std::vector F0_; + + Mat duneB_; + Mat duneC_; + Mat invDuneD_; + + BVector resWell_; + + mutable BVector Cx_; + mutable BVector invDrw_; + mutable BVector scaleAddRes_; + + // protected methods + EvalWell getBhp(const int wellIdx) const { + const WellControls* wc = wells().ctrls[wellIdx]; + if (well_controls_get_current_type(wc) == BHP) { + EvalWell bhp = 0.0; + const double target_rate = well_controls_get_current_target(wc); + bhp.setValue(target_rate); + return bhp; + } else if (well_controls_get_current_type(wc) == THP) { + const int control = well_controls_get_current(wc); + const double thp = well_controls_get_current_target(wc); + const double alq = well_controls_iget_alq(wc, control); + const int table_id = well_controls_iget_vfp(wc, control); + EvalWell aqua = 0.0; + EvalWell liquid = 0.0; + EvalWell vapour = 0.0; + EvalWell bhp = 0.0; + double vfp_ref_depth = 0.0; + + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + + if ((*active_)[ Water ]) { + aqua = getQs(wellIdx, pu.phase_pos[ Water]); + } + if ((*active_)[ Oil ]) { + liquid = getQs(wellIdx, pu.phase_pos[ Oil ]); + } + if ((*active_)[ Gas ]) { + vapour = getQs(wellIdx, pu.phase_pos[ Gas ]); + } + if (wells().type[wellIdx] == INJECTOR) { + bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp); + vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + } else { + bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); + vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + } + + // pick the density in the top layer + const int perf = wells().well_connpos[wellIdx]; + const double rho = well_perforation_densities_[perf]; + const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_); + bhp -= dp; + return bhp; + + } + return wellVariables_[wellIdx]; + } + + EvalWell getQs(const int wellIdx, const int phaseIdx) const { + EvalWell qs = 0.0; + const WellControls* wc = wells().ctrls[wellIdx]; + const int np = wells().number_of_phases; + const double target_rate = well_controls_get_current_target(wc); + + if (wells().type[wellIdx] == INJECTOR) { + const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; + if (comp_frac == 0.0) + return qs; + + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { + return wellVariables_[wellIdx]; + } + qs.setValue(target_rate); + return qs; + } + + // Producers + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { + return wellVariables_[wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); + } + if (well_controls_get_current_type(wc) == SURFACE_RATE) { + const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; + + if (comp_frac == 1.0) { + qs.setValue(target_rate); + return qs; + } + int currentControlIdx = 0; + for (int i = 0; i < np; ++i) { + currentControlIdx += wells().comp_frac[np*wellIdx + i] * i; + } + + const double eps = 1e-6; + if (wellVolumeFractionScaled(wellIdx,currentControlIdx) < eps) { + return qs; + } + return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx)); + } + // ReservoirRate + return target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx); + } + + EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const { + assert(wells().number_of_phases == 3); + const int nw = wells().number_of_wells; + if (phaseIdx == Water) { + return wellVariables_[nw + wellIdx]; + } + + if (phaseIdx == Gas) { + return wellVariables_[2*nw + wellIdx]; + } + + // Oil + return 1.0 - wellVariables_[nw + wellIdx] - wellVariables_[2 * nw + wellIdx]; + } + + EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const { + const WellControls* wc = wells().ctrls[wellIdx]; + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + const double* distr = well_controls_get_current_distr(wc); + return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx]; + } + std::vector g = {1,1,0.01}; + return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); + } + + + + template + bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const int well_number) const + { + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int np = well_state.numPhases(); + + if (econ_production_limits.onMinOilRate()) { + assert((*active_)[Oil]); + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double min_oil_rate = econ_production_limits.minOilRate(); + if (std::abs(oil_rate) < min_oil_rate) { + return true; + } + } + + if (econ_production_limits.onMinGasRate() ) { + assert((*active_)[Gas]); + const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ]; + const double min_gas_rate = econ_production_limits.minGasRate(); + if (std::abs(gas_rate) < min_gas_rate) { + return true; + } + } + + if (econ_production_limits.onMinLiquidRate() ) { + assert((*active_)[Oil]); + assert((*active_)[Water]); + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + const double min_liquid_rate = econ_production_limits.minLiquidRate(); + if (std::abs(liquid_rate) < min_liquid_rate) { + return true; + } + } + + if (econ_production_limits.onMinReservoirFluidRate()) { + OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet"); + } + + return false; + } + + + using WellMapType = typename WellState::WellMapType; + using WellMapEntryType = typename WellState::mapentry_t; + + // a tuple type for ratio limit check. + // first value indicates whether ratio limit is violated, when the ratio limit is not violated, the following three + // values should not be used. + // second value indicates whehter there is only one connection left. + // third value indicates the indx of the worst-offending connection. + // the last value indicates the extent of the violation for the worst-offending connection, which is defined by + // the ratio of the actual value to the value of the violated limit. + using RatioCheckTuple = std::tuple; + + enum ConnectionIndex { + INVALIDCONNECTION = -10000 + }; + + + template + RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const WellMapEntryType& map_entry) const + { + // TODO: not sure how to define the worst-offending connection when more than one + // ratio related limit is violated. + // The defintion used here is that we define the violation extent based on the + // ratio between the value and the corresponding limit. + // For each violated limit, we decide the worst-offending connection separately. + // Among the worst-offending connections, we use the one has the biggest violation + // extent. + + bool any_limit_violated = false; + bool last_connection = false; + int worst_offending_connection = INVALIDCONNECTION; + double violation_extent = -1.0; + + if (econ_production_limits.onMaxWaterCut()) { + const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry); + bool water_cut_violated = std::get<0>(water_cut_return); + if (water_cut_violated) { + any_limit_violated = true; + const double violation_extent_water_cut = std::get<3>(water_cut_return); + if (violation_extent_water_cut > violation_extent) { + violation_extent = violation_extent_water_cut; + worst_offending_connection = std::get<2>(water_cut_return); + last_connection = std::get<1>(water_cut_return); + } + } + } + + if (econ_production_limits.onMaxGasOilRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxWaterGasRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!"); + } + + if (econ_production_limits.onMaxGasLiquidRatio()) { + OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!"); + } + + if (any_limit_violated) { + assert(worst_offending_connection >=0); + assert(violation_extent > 1.); + } + + return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + template + RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, + const WellState& well_state, + const WellMapEntryType& map_entry) const + { + bool water_cut_limit_violated = false; + int worst_offending_connection = INVALIDCONNECTION; + bool last_connection = false; + double violation_extent = -1.0; + + const int np = well_state.numPhases(); + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int well_number = map_entry[0]; + + assert((*active_)[Oil]); + assert((*active_)[Water]); + + const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; + const double liquid_rate = oil_rate + water_rate; + double water_cut; + if (std::abs(liquid_rate) != 0.) { + water_cut = water_rate / liquid_rate; + } else { + water_cut = 0.0; + } + + const double max_water_cut_limit = econ_production_limits.maxWaterCut(); + if (water_cut > max_water_cut_limit) { + water_cut_limit_violated = true; + } + + if (water_cut_limit_violated) { + // need to handle the worst_offending_connection + const int perf_start = map_entry[1]; + const int perf_number = map_entry[2]; + + std::vector water_cut_perf(perf_number); + for (int perf = 0; perf < perf_number; ++perf) { + const int i_perf = perf_start + perf; + const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ]; + const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ]; + const double liquid_perf_rate = oil_perf_rate + water_perf_rate; + if (std::abs(liquid_perf_rate) != 0.) { + water_cut_perf[perf] = water_perf_rate / liquid_perf_rate; + } else { + water_cut_perf[perf] = 0.; + } + } + + last_connection = (perf_number == 1); + if (last_connection) { + worst_offending_connection = 0; + violation_extent = water_cut_perf[0] / max_water_cut_limit; + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + double max_water_cut_perf = 0.; + for (int perf = 0; perf < perf_number; ++perf) { + if (water_cut_perf[perf] > max_water_cut_perf) { + worst_offending_connection = perf; + max_water_cut_perf = water_cut_perf[perf]; + } + } + + assert(max_water_cut_perf != 0.); + assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number)); + + violation_extent = max_water_cut_perf / max_water_cut_limit; + } + + return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent); + } + + }; + + +} // namespace Opm +#endif diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 983ecce7d..10c65e7de 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -675,17 +675,18 @@ namespace Opm int table_id = well_controls_iget_vfp(wc, ctrl_index); const WellType& well_type = wells().type[w]; + const int perf = wells().well_connpos[w]; //first perforation. if (well_type == INJECTOR) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities(), gravity_); + wellPerforationDensities()[perf], gravity_); well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); } else if (well_type == PRODUCER) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(), - wellPerforationDensities(), gravity_); + wellPerforationDensities()[perf], gravity_); well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); } @@ -790,17 +791,18 @@ namespace Opm //Set *BHP* target by calculating bhp from THP const WellType& well_type = wells().type[w]; + const int perf = wells().well_connpos[w]; // first perforation. if (well_type == INJECTOR) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); + wellPerforationDensities()[perf], gravity_); xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; } else if (well_type == PRODUCER) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); + wellPerforationDensities()[perf], gravity_); xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; } @@ -1116,11 +1118,12 @@ namespace Opm //Set *BHP* target by calculating bhp from THP const WellType& well_type = wells().type[w]; + const int perf = wells().well_connpos[w]; //first perforation if (well_type == INJECTOR) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); + wellPerforationDensities()[perf], gravity_); const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; // apply the strictest of the bhp controlls i.e. smallest bhp for injectors if ( bhp < bhps[w]) { @@ -1130,7 +1133,7 @@ namespace Opm else if (well_type == PRODUCER) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity_); + wellPerforationDensities()[perf], gravity_); const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; // apply the strictest of the bhp controlls i.e. largest bhp for producers diff --git a/opm/autodiff/VFPHelpers.hpp b/opm/autodiff/VFPHelpers.hpp index 6220b71de..6449dc805 100644 --- a/opm/autodiff/VFPHelpers.hpp +++ b/opm/autodiff/VFPHelpers.hpp @@ -25,7 +25,8 @@ #include #include #include - +#include +#include /** * This file contains a set of helper functions used by VFPProd / VFPInj. @@ -35,14 +36,7 @@ namespace detail { typedef AutoDiffBlock ADB; - - - - - - - - +typedef DenseAd::Evaluation EvalWell; /** @@ -52,7 +46,12 @@ inline double zeroIfNan(const double& value) { return (std::isnan(value)) ? 0.0 : value; } - +/** + * Returns zero if input value is NaN + */ +inline double zeroIfNan(const EvalWell& value) { + return (std::isnan(value.value())) ? 0.0 : value.value(); +} diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index 7027c9957..68e758857 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -89,7 +89,36 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, +VFPInjProperties::EvalWell VFPInjProperties::bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp) const { + //Get the table + const VFPInjTable* table = detail::getTable(m_tables, table_id); + EvalWell bhp = 0.0; + + //Find interpolation variables + EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + + //Compute the BHP for each well independently + if (table != nullptr) { + //First, find the values to interpolate between + //Value of FLO is negative in OPM for producers, but positive in VFP table + auto flo_i = detail::findInterpData(flo.value(), table->getFloAxis()); + auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant + + detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i); + + bhp = bhp_val.dflo * flo; + bhp.setValue(bhp_val.value); // thp is assumed constant i.e. + } + else { + bhp.setValue(-1e100); //Signal that this value has not been calculated properly, due to "missing" table + } + return bhp; +} diff --git a/opm/autodiff/VFPInjProperties.hpp b/opm/autodiff/VFPInjProperties.hpp index 7a6e4c71e..20c79af19 100644 --- a/opm/autodiff/VFPInjProperties.hpp +++ b/opm/autodiff/VFPInjProperties.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include @@ -91,6 +93,13 @@ public: const ADB& vapour, const ADB& thp) const; + typedef DenseAd::Evaluation EvalWell; + EvalWell bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp) const; + /** * Linear interpolation of bhp as a function of the input parameters * @param table_id Table number to use diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index 8da911c72..c9b8aebc4 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -24,7 +24,8 @@ #include #include #include - +#include +#include #include @@ -74,7 +75,42 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, } +VFPProdProperties::EvalWell VFPProdProperties::bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp, + const double& alq) const { + //Get the table + const VFPProdTable* table = detail::getTable(m_tables, table_id); + EvalWell bhp = 0.0; + + //Find interpolation variables + EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + EvalWell wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); + EvalWell gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); + + //Compute the BHP for each well independently + if (table != nullptr) { + //First, find the values to interpolate between + //Value of FLO is negative in OPM for producers, but positive in VFP table + auto flo_i = detail::findInterpData(-flo.value(), table->getFloAxis()); + auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant + auto wfr_i = detail::findInterpData( wfr.value(), table->getWFRAxis()); + auto gfr_i = detail::findInterpData( gfr.value(), table->getGFRAxis()); + auto alq_i = detail::findInterpData( alq, table->getALQAxis()); //assume constant + + detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + + bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo * flo); + bhp.setValue(bhp_val.value); + } + else { + bhp.setValue(-1e100); //Signal that this value has not been calculated properly, due to "missing" table + } + return bhp; +} diff --git a/opm/autodiff/VFPProdProperties.hpp b/opm/autodiff/VFPProdProperties.hpp index 76580ecdf..33a5bbba7 100644 --- a/opm/autodiff/VFPProdProperties.hpp +++ b/opm/autodiff/VFPProdProperties.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include @@ -99,6 +101,15 @@ public: const ADB& thp, const ADB& alq) const; + + typedef DenseAd::Evaluation EvalWell; + EvalWell bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const double& thp, + const double& alq) const; + /** * Linear interpolation of bhp as a function of the input parameters * @param table_id Table number to use diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp index 23f423915..6e0728c88 100644 --- a/opm/autodiff/WellHelpers.hpp +++ b/opm/autodiff/WellHelpers.hpp @@ -133,7 +133,7 @@ namespace Opm { */ inline double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth, - const Vector& well_perforation_densities, const double gravity) { + const double& rho, const double gravity) { if ( wells.well_connpos[w] == wells.well_connpos[w+1] ) { // This is a well with no perforations. @@ -144,13 +144,12 @@ namespace Opm { } const double well_ref_depth = wells.depth_ref[w]; const double dh = vfp_ref_depth - well_ref_depth; - const int perf = wells.well_connpos[w]; - const double rho = well_perforation_densities[perf]; const double dp = rho*gravity*dh; return dp; } + inline Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth, const Vector& well_perforation_densities, const double gravity) { @@ -159,7 +158,23 @@ namespace Opm { #pragma omp parallel for schedule(static) for (int i=0; i computeHydrostaticCorrection(const Wells& wells, const std::vector& vfp_ref_depth, + const std::vector& well_perforation_densities, const double gravity) { + const int nw = wells.number_of_wells; + std::vector retval(nw,0.0); + +#pragma omp parallel for schedule(static) + for (int i=0; i. +*/ + +#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED +#define OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + /// The state of a set of wells, tailored for use by the fully + /// implicit blackoil simulator. + class WellStateFullyImplicitBlackoilDense + : public WellStateFullyImplicitBlackoil + { + typedef WellStateFullyImplicitBlackoil BaseType; + public: + typedef BaseType :: WellMapType WellMapType; + + using BaseType :: wellRates; + using BaseType :: bhp; + using BaseType :: perfPress; + using BaseType :: wellMap; + using BaseType :: numWells; + using BaseType :: numPhases; + using BaseType :: perfPhaseRates; + using BaseType :: currentControls; + using BaseType :: wellPotentials; + + /// Allocate and initialize if wells is non-null. Also tries + /// to give useful initial values to the bhp(), wellRates() + /// and perfPhaseRates() fields, depending on controls + template + void init(const Wells* wells, const State& state, const PrevState& prevState) + { + // call init on base class + BaseType :: init(wells, state, prevState); + + // if there are no well, do nothing in init + if (wells == 0) { + 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]; + + // 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]; + + 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]; + } + } 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; + } + assert(np == 3); + double total_rates = 0.0; + for (int p = 0; p < np; ++p) { + total_rates += g[p] * wellRates()[np*w + p]; + } + + if(std::abs(total_rates) > 0) { + wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water; + wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas + } else { + wellSolutions()[nw + w] = wells->comp_frac[np*w + Water]; + wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas]; + } + } + + + // 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 ]; + } + + } + } + } + } + + template + void resize(const Wells* wells, const State& state) { + const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes + init(wells, state, dummy_state) ; + } + + + /// One rate per phase and well connection. + std::vector& wellSolutions() { return well_solutions_; } + const std::vector& wellSolutions() const { return well_solutions_; } + + data::Wells report(const PhaseUsage& pu) const override { + data::Wells res = BaseType::report(pu); + return res; + } + + private: + std::vector well_solutions_; + }; + +} // namespace Opm + + +#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED