diff --git a/opm/autodiff/BlackoilDetails.hpp b/opm/autodiff/BlackoilDetails.hpp index d0afc2e8a..8a95b64fe 100644 --- a/opm/autodiff/BlackoilDetails.hpp +++ b/opm/autodiff/BlackoilDetails.hpp @@ -291,6 +291,91 @@ namespace detail { } } + + 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 +#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] = 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. diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index a0f5b5eaa..1a83b7046 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -51,7 +51,6 @@ #include #include #include -#include #include #include #include @@ -104,10 +103,6 @@ namespace Opm { typedef BlackoilModelEbos ThisType; public: // --------- Types and enums --------- - typedef AutoDiffBlock ADB; - typedef ADB::V V; - typedef ADB::M M; - typedef BlackoilState ReservoirState; typedef WellStateFullyImplicitBlackoil WellState; typedef BlackoilModelParameters ModelParameters; @@ -124,8 +119,8 @@ namespace Opm { typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; typedef double Scalar; - typedef Dune::FieldVector VectorBlockType; - typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector BVector; //typedef typename SolutionVector :: value_type PrimaryVariables ; @@ -635,6 +630,8 @@ namespace Opm { /// \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_; @@ -643,50 +640,56 @@ namespace Opm { const int nc = Opm::AutoDiffGrid::numCells(grid_); const int np = numPhases(); - const V& pv = geo_.poreVolume(); + const auto& pv = geo_.poreVolume(); - std::vector R_sum(np); - std::vector B_avg(np); - std::vector maxCoeff(np); - std::vector maxNormWell(np); - Eigen::Array B(nc, np); - Eigen::Array R(nc, np); - Eigen::Array R2(nc, np); - Eigen::Array tempV(nc, np); + 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(); - auto ebosResid = ebosSimulator_.model().linearizer().residual(); for ( int idx = 0; idx < np; ++idx ) { - V b(nc); - V r(nc); + 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(); - int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); - int ebosCompIdx = flowPhaseToEbosCompIdx(idx); - - b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; - r[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; - + B_idx [cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; + R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; } - R2.col(idx) = r; - B.col(idx) = b; } for ( int idx = 0; idx < np; ++idx ) { - tempV.col(idx) = R2.col(idx).abs()/pv; + //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 ]; + } } - std::vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + Vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); + Vector wellResidual = wellModel().residual(); const double pvSum = detail::convergenceReduction(B, tempV, R2, R_sum, maxCoeff, B_avg, maxNormWell, - nc, np, pv_vector, wellModel().residual()); + nc, np, pv_vector, wellResidual ); - std::vector CNV(np); - std::vector mass_balance_residual(np); - std::vector well_flux_residual(np); + Vector CNV(np); + Vector mass_balance_residual(np); + Vector well_flux_residual(np); bool converged_MB = true; bool converged_CNV = true; @@ -787,13 +790,6 @@ namespace Opm { protected: - // --------- Types and enums --------- - - typedef Eigen::Array DataBlock; - // --------- Data members --------- Simulator& ebosSimulator_; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 945c6b8f4..2a24c3364 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -215,7 +215,8 @@ public: props_.permeability(), dynamic_list_econ_limited, is_parallel_run_, - well_potentials); + well_potentials ); + const Wells* wells = wells_manager.c_wells(); WellState well_state; well_state.init(wells, state, prev_well_state); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp index 0415ab26e..01c2e686b 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp @@ -30,7 +30,6 @@ #include #include #include -#include #include #include diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index fb916c901..271982afe 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -452,7 +452,7 @@ namespace Opm { /// Diff to bhp for each well perforation. std::vector wellPerforationPressureDiffs() const - { + { return well_perforation_pressure_diffs_; } @@ -683,9 +683,10 @@ namespace Opm { const int nw = wells().number_of_wells; std::vector res(np*nw); for( int p=0; p