diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 2f50c79d7..9a0e78a2e 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -72,7 +72,6 @@ BEGIN_PROPERTIES NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem, FlowNonLinearSolver, FlowIstlSolver, FlowModelParameters, FlowTimeSteppingParameters)); SET_STRING_PROP(EclFlowProblem, OutputDir, ""); -SET_BOOL_PROP(EclFlowProblem, DisableWells, true); SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); // default in flow is to formulate the equations in surface volumes SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume, true); @@ -87,6 +86,8 @@ SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false); SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true); SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false); +SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel); + END_PROPERTIES namespace Opm { @@ -209,8 +210,6 @@ namespace Opm { wasSwitched_.resize(numDof); std::fill(wasSwitched_.begin(), wasSwitched_.end(), false); - wellModel().beginTimeStep(timer.reportStepNum(), timer.simulationTimeElapsed()); - if (param_.update_equations_scaling_) { std::cout << "equation scaling not suported yet" << std::endl; //updateEquationsScaling(); @@ -302,7 +301,7 @@ namespace Opm { // handling well state update before oscillation treatment is a decision based // on observation to avoid some big performance degeneration under some circumstances. // there is no theorectical explanation which way is better for sure. - wellModel().recoverWellSolutionAndUpdateWellState(x); + wellModel().postSolve(x); if (param_.use_update_stabilization_) { // Stabilize the nonlinear update. @@ -343,9 +342,7 @@ namespace Opm { /// \param[in] timer simulation timer void afterStep(const SimulatorTimerInterface& OPM_UNUSED timer) { - wellModel().timeStepSucceeded(timer.simulationTimeElapsed()); ebosSimulator_.problem().endTimeStep(); - } /// Assemble the residual and Jacobian of the nonlinear system. @@ -361,28 +358,12 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - // -------- Current time step length ---------- - const double dt = timer.currentStepLength(); - - // -------- Well equations ---------- - - try - { - // assembles the well equations and applies the wells to - // the reservoir equations as a source term. - wellModel().assemble(iterationIdx, dt); - } - catch ( const Dune::FMatrixError& ) - { - OPM_THROW(Opm::NumericalIssue,"Error encounted when solving well equations"); - } - auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); if (param_.matrix_add_well_contributions_) { wellModel().addWellContributions(ebosJac); } if ( param_.preconditioner_add_well_contributions_ && - ! param_.matrix_add_well_contributions_ ) { + ! param_.matrix_add_well_contributions_ ) { matrix_for_preconditioner_ .reset(new Mat(ebosJac)); wellModel().addWellContributions(*matrix_for_preconditioner_); } @@ -504,9 +485,6 @@ namespace Opm { /// r is the residual. void solveJacobianSystem(BVector& x) const { - const auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator_.model().linearizer().residual(); - // J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well // with the reservoir and D is the wells itself. // The full system is reduced to a number of cells X number of cells system via Schur complement @@ -516,7 +494,9 @@ namespace Opm { // r = [r, r_well], where r is the residual and r_well the well residual. // r -= B^T * D^-1 r_well - // apply well residual to the residual. + auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + wellModel().apply(ebosResid); // set initial guess @@ -1001,9 +981,9 @@ namespace Opm { const BlackoilWellModel& wellModel() const { return well_model_; } - void beginReportStep() + void beginReportStep(bool isRestart) { - ebosSimulator_.problem().beginEpisode(); + ebosSimulator_.problem().beginEpisode(isRestart); } void endReportStep() diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index a58f623ab..ac3335675 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -24,6 +24,7 @@ #ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED #define OPM_BLACKOILWELLMODEL_HEADER_INCLUDED +#include #include #include @@ -60,12 +61,18 @@ #include +BEGIN_PROPERTIES + +NEW_PROP_TAG(EnableTerminalOutput); + +END_PROPERTIES namespace Opm { /// Class for handling the blackoil well model. template - class BlackoilWellModel { + class BlackoilWellModel : public Ewoms::BaseAuxiliaryModule + { public: // --------- Types --------- typedef WellStateFullyImplicitBlackoil WellState; @@ -77,6 +84,11 @@ namespace Opm { typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector; + typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; + + typedef typename Ewoms::BaseAuxiliaryModule::NeighborSet NeighborSet; static const int numEq = Indices::numEq; static const int solventSaturationIdx = Indices::solventSaturationIdx; @@ -101,48 +113,94 @@ namespace Opm { using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; - BlackoilWellModel(Simulator& ebosSimulator, - const ModelParameters& param, - const bool terminal_output); + BlackoilWellModel(Simulator& ebosSimulator); - void initFromRestartFile(const RestartValue& restartValues) + void init(const Opm::EclipseState& eclState, const Opm::Schedule& schedule); + + ///////////// + // + ///////////// + unsigned numDofs() const + // No extra dofs are inserted for wells. (we use a Schur complement.) + { return 0; } + + void addNeighbors(std::vector& neighbors) const; + + void applyInitial() + {} + + void linearize(JacobianMatrix& mat , GlobalEqVector& res); + + void postSolve(GlobalEqVector& deltaX) { - // gives a dummy dynamic_list_econ_limited - DynamicListEconLimited dummyListEconLimited; - const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames(); - WellsManager wellsmanager(eclState(), - schedule(), - // The restart step value is used to identify wells present at the given - // time step. Wells that are added at the same time step as RESTART is initiated - // will not be present in a restart file. Use the previous time step to retrieve - // wells that have information written to the restart file. - std::max(eclState().getInitConfig().getRestartStep() - 1, 0), - Opm::UgGridHelpers::numCells(grid()), - Opm::UgGridHelpers::globalCell(grid()), - Opm::UgGridHelpers::cartDims(grid()), - Opm::UgGridHelpers::dimensions(grid()), - Opm::UgGridHelpers::cell2Faces(grid()), - Opm::UgGridHelpers::beginFaceCentroids(grid()), - dummyListEconLimited, - grid().comm().size() > 1, - defunctWellNames); - - const Wells* wells = wellsmanager.c_wells(); - - const int nw = wells->number_of_wells; - if (nw > 0) { - auto phaseUsage = phaseUsageFromDeck(eclState()); - size_t numCells = Opm::UgGridHelpers::numCells(grid()); - well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step - wellsToState(restartValues.wells, phaseUsage, well_state_); - previous_well_state_ = well_state_; - } + recoverWellSolutionAndUpdateWellState(deltaX); } - // compute the well fluxes and assemble them in to the reservoir equations as source terms - // and in the well equations. - void assemble(const int iterationIdx, - const double dt); + ///////////// + // + ///////////// + + template + void deserialize(Restarter& res) + { + // TODO (?) + } + + /*! + * \brief This method writes the complete state of the well + * to the harddisk. + */ + template + void serialize(Restarter& res) + { + // TODO (?) + } + + void beginEpisode(const Opm::EclipseState& eclState, + const Opm::Schedule& schedule, + bool isRestart) + { + size_t episodeIdx = ebosSimulator_.episodeIndex(); + // beginEpisode in eclProblem advances the episode index + // we don't want this when we are at the beginning of an + // restart. + if (isRestart) + episodeIdx -= 1; + + beginReportStep(episodeIdx); + } + + void beginTimeStep(); + + void beginIteration() + { + assemble(ebosSimulator_.model().newtonMethod().numIterations(), + ebosSimulator_.timeStepSize()); + } + + void endIteration() + { } + + void endTimeStep() + { + timeStepSucceeded(ebosSimulator_.time()); + } + + void endEpisode() + { + endReportStep(); + } + + template + void computeTotalRatesForDof(RateVector& rate, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) const; + + void initFromRestartFile(const RestartValue& restartValues); + + Opm::data::Wells wellData() const + { return well_state_.report(phase_usage_, Opm::UgGridHelpers::globalCell(grid())); } // substract Binv(D)rw from r; void apply( BVector& r) const; @@ -153,10 +211,6 @@ namespace Opm { // apply well model with scaling of alpha void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const; - // using the solution x to recover the solution xw for wells and applying - // xw to update Well State - void recoverWellSolutionAndUpdateWellState(const BVector& x); - // Check if well equations is converged. bool getWellConvergence(const std::vector& B_avg) const; @@ -172,31 +226,20 @@ namespace Opm { // return the internal well state const WellState& wellState() const; - // only use this for restart. - void setRestartWellState(const WellState& well_state); - - // called at the beginning of a time step - void beginTimeStep(const int timeStepIdx,const double simulationTime); - // called at the end of a time step - void timeStepSucceeded(const double& simulationTime); - - // called at the beginning of a report step - void beginReportStep(const int time_step); - - // called at the end of a report step - void endReportStep(); - const SimulatorReport& lastReport() const; - void addWellContributions(Mat& mat) { for ( const auto& well: well_container_ ) { - well->addWellContributions(mat); + well->addWellContributions(mat); } } + // called at the beginning of a report step + void beginReportStep(const int time_step); + protected: + void extractLegacyPressure_(std::vector& cellPressure) const { size_t nc = number_of_cells_; @@ -234,6 +277,9 @@ namespace Opm { // a vector of all the wells. std::vector well_container_; + // map from logically cartesian cell indices to compressed ones + std::vector cartesian_to_compressed_; + // create the well container std::vector createWellContainer(const int time_step); @@ -276,6 +322,21 @@ namespace Opm { const Schedule& schedule() const { return ebosSimulator_.vanguard().schedule(); } + // compute the well fluxes and assemble them in to the reservoir equations as source terms + // and in the well equations. + void assemble(const int iterationIdx, + const double dt); + + // called at the end of a time step + void timeStepSucceeded(const double& simulationTime); + + // called at the end of a report step + void endReportStep(); + + // using the solution x to recover the solution xw for wells and applying + // xw to update Well State + void recoverWellSolutionAndUpdateWellState(const BVector& x); + void updateWellControls(); void updateGroupControls(); @@ -283,7 +344,7 @@ namespace Opm { // setting the well_solutions_ based on well_state. void updatePrimaryVariables(); - void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const; + void setupCartesianToCompressed_(const int* global_cell, int number_of_cells); void computeRepRadiusPerfLength(const Grid& grid); @@ -322,8 +383,7 @@ namespace Opm { void resetWellControlFromState() const; - void assembleWellEq(const double dt, - bool only_wells); + void assembleWellEq(const double dt); // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index d5134ad7e..f7c3b319e 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -4,16 +4,26 @@ namespace Opm { template BlackoilWellModel:: - BlackoilWellModel(Simulator& ebosSimulator, - const ModelParameters& param, - const bool terminal_output) + BlackoilWellModel(Simulator& ebosSimulator) : ebosSimulator_(ebosSimulator) - , param_(param) - , terminal_output_(terminal_output) , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) { - const auto& eclState = ebosSimulator_.vanguard().eclState(); + terminal_output_ = false; + if (ebosSimulator.gridView().comm().rank() == 0) + terminal_output_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput); + } + + template + void + BlackoilWellModel:: + init(const Opm::EclipseState& eclState, const Opm::Schedule& schedule) + { + gravity_ = ebosSimulator_.problem().gravity()[2]; + + extractLegacyCellPvtRegionIndex_(); + extractLegacyDepth_(); + phase_usage_ = phaseUsageFromDeck(eclState); const auto& gridView = ebosSimulator_.gridView(); @@ -28,6 +38,79 @@ namespace Opm { extractLegacyCellPvtRegionIndex_(); extractLegacyDepth_(); initial_step_ = true; + + const auto& grid = ebosSimulator_.vanguard().grid(); + const auto& cartDims = Opm::UgGridHelpers::cartDims(grid); + setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid), + cartDims[0]*cartDims[1]*cartDims[2]); + + // add the eWoms auxiliary module for the wells to the list + ebosSimulator_.model().addAuxiliaryModule(this); + } + + template + void + BlackoilWellModel:: + addNeighbors(std::vector& neighbors) const + { + if (!param_.matrix_add_well_contributions_) { + return; + } + + // Create cartesian to compressed mapping + int last_time_step = schedule().getTimeMap().size() - 1; + const auto& schedule_wells = schedule().getWells(); + const auto& cartesianSize = Opm::UgGridHelpers::cartDims(grid()); + + // initialize the additional cell connections introduced by wells. + for (const auto well : schedule_wells) + { + std::vector wellCells; + // All possible connections of the well + const auto& connectionSet = well->getConnections(last_time_step); + wellCells.reserve(connectionSet.size()); + + for ( size_t c=0; c < connectionSet.size(); c++ ) + { + const auto& connection = connectionSet.get(c); + int i = connection.getI(); + int j = connection.getJ(); + int k = connection.getK(); + int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k); + int compressed_idx = cartesian_to_compressed_.at(cart_grid_idx); + + if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells. + wellCells.push_back(compressed_idx); + } + } + + for (int cellIdx : wellCells) { + neighbors[cellIdx].insert(wellCells.begin(), + wellCells.end()); + } + } + } + + template + void + BlackoilWellModel:: + linearize(JacobianMatrix& mat , GlobalEqVector& res) + { + if (!localWellsActive()) + return; + + // we don't what to add the schur complement + // here since it affects the getConvergence method + /* + for (const auto& well: well_container_) { + if (param_.matrix_add_well_contributions_) + well->addWellContributions(mat); + + // applying the well residual to reservoir residuals + // r = r - duneC_^T * invDuneD_ * resWell_ + well->apply(res); + } + */ } @@ -126,15 +209,17 @@ namespace Opm { template void BlackoilWellModel:: - beginTimeStep(const int timeStepIdx, const double simulationTime) { + beginTimeStep() { well_state_ = previous_well_state_; + const int reportStepIdx = ebosSimulator_.episodeIndex(); + const double simulationTime = ebosSimulator_.time(); // test wells - wellTesting(timeStepIdx, simulationTime); + wellTesting(reportStepIdx, simulationTime); // create the well container - well_container_ = createWellContainer(timeStepIdx); + well_container_ = createWellContainer(reportStepIdx); // do the initialization for all the wells // TODO: to see whether we can postpone of the intialization of the well containers to @@ -205,12 +290,6 @@ namespace Opm { } } - // only use this for restart. - template - void - BlackoilWellModel:: - setRestartWellState(const WellState& well_state) { previous_well_state_ = well_state; } - // called at the end of a report step template void @@ -238,6 +317,60 @@ namespace Opm { previous_well_state_ = well_state_; } + + template + template + void + BlackoilWellModel:: + computeTotalRatesForDof(RateVector& rate, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) const + { + rate = 0; + int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + for (const auto& well : well_container_) + well->addCellRates(rate, elemIdx); + } + + template + void + BlackoilWellModel:: + initFromRestartFile(const RestartValue& restartValues) + { + // gives a dummy dynamic_list_econ_limited + DynamicListEconLimited dummyListEconLimited; + const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames(); + WellsManager wellsmanager(eclState(), + schedule(), + // The restart step value is used to identify wells present at the given + // time step. Wells that are added at the same time step as RESTART is initiated + // will not be present in a restart file. Use the previous time step to retrieve + // wells that have information written to the restart file. + std::max(eclState().getInitConfig().getRestartStep() - 1, 0), + Opm::UgGridHelpers::numCells(grid()), + Opm::UgGridHelpers::globalCell(grid()), + Opm::UgGridHelpers::cartDims(grid()), + Opm::UgGridHelpers::dimensions(grid()), + Opm::UgGridHelpers::cell2Faces(grid()), + Opm::UgGridHelpers::beginFaceCentroids(grid()), + dummyListEconLimited, + grid().comm().size() > 1, + defunctWellNames); + + const Wells* wells = wellsmanager.c_wells(); + + const int nw = wells->number_of_wells; + if (nw > 0) { + auto phaseUsage = phaseUsageFromDeck(eclState()); + size_t numCells = Opm::UgGridHelpers::numCells(grid()); + well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step + wellsToState(restartValues.wells, phaseUsage, well_state_); + previous_well_state_ = well_state_; + } + initial_step_ = false; + } + template std::vector::WellInterfacePtr > BlackoilWellModel:: @@ -401,7 +534,7 @@ namespace Opm { // basically, this is a more updated state from the solveWellEq based on fixed // reservoir state, will tihs be a better place to inialize the explict information? } - assembleWellEq(dt, false); + assembleWellEq(dt); last_report_.converged = true; } @@ -413,16 +546,13 @@ namespace Opm { template void BlackoilWellModel:: - assembleWellEq(const double dt, - bool only_wells) + assembleWellEq(const double dt) { for (auto& well : well_container_) { - well->assembleWellEq(ebosSimulator_, dt, well_state_, only_wells); + well->assembleWellEq(ebosSimulator_, dt, well_state_); } } - // applying the well residual to reservoir residuals - // r = r - duneC_^T * invDuneD_ * resWell_ template void BlackoilWellModel:: @@ -438,9 +568,6 @@ namespace Opm { } - - - // Ax = A x - C D^-1 B x template void @@ -584,7 +711,7 @@ namespace Opm { int it = 0; bool converged; do { - assembleWellEq(dt, true); + assembleWellEq(dt); converged = getWellConvergence(B_avg); @@ -1071,16 +1198,17 @@ namespace Opm { template void BlackoilWellModel:: - setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const + setupCartesianToCompressed_(const int* global_cell, int number_of_cartesian_cells) { + cartesian_to_compressed_.resize(number_of_cartesian_cells, -1); if (global_cell) { - for (int i = 0; i < number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + for (unsigned i = 0; i < number_of_cells_; ++i) { + cartesian_to_compressed_[global_cell[i]] = i; } } else { - for (int i = 0; i < number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(i, i)); + for (unsigned i = 0; i < number_of_cells_; ++i) { + cartesian_to_compressed_[i] = i; } } @@ -1091,16 +1219,8 @@ namespace Opm { BlackoilWellModel:: computeRepRadiusPerfLength(const Grid& grid) { - // TODO, the function does not work for parallel running - // to be fixed later. - const int* global_cell = Opm::UgGridHelpers::globalCell(grid); - - std::map cartesian_to_compressed; - setupCompressedToCartesian(global_cell, number_of_cells_, - cartesian_to_compressed); - for (const auto& well : well_container_) { - well->computeRepRadiusPerfLength(grid, cartesian_to_compressed); + well->computeRepRadiusPerfLength(grid, cartesian_to_compressed_); } } diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 4f22e9c93..f874e0fd6 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -437,6 +437,7 @@ namespace Opm void setupEbosSimulator() { ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false)); + ebosSimulator_->executionTimer().start(); ebosSimulator_->model().applyInitialSolution(); try { diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp index 187a7c716..cb7dd7907 100644 --- a/opm/autodiff/MultisegmentWell.hpp +++ b/opm/autodiff/MultisegmentWell.hpp @@ -33,6 +33,7 @@ namespace Opm { public: typedef WellInterface Base; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; using typename Base::WellState; using typename Base::Simulator; @@ -107,44 +108,43 @@ namespace Opm virtual void init(const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, - const int num_cells); + const int num_cells) override; - virtual void initPrimaryVariablesEvaluation() const; + virtual void initPrimaryVariablesEvaluation() const override; - virtual void assembleWellEq(Simulator& ebosSimulator, + virtual void assembleWellEq(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells); + WellState& well_state) override; /// updating the well state based the control mode specified with current // TODO: later will check wheter we need current - virtual void updateWellStateWithTarget(WellState& well_state) const; + virtual void updateWellStateWithTarget(WellState& well_state) const override; /// check whether the well equations get converged for this well - virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const; + virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const override; /// Ax = Ax - C D^-1 B x - virtual void apply(const BVector& x, BVector& Ax) const; + virtual void apply(const BVector& x, BVector& Ax) const override; /// r = r - C D^-1 Rw - virtual void apply(BVector& r) const; + virtual void apply(BVector& r) const override; /// using the solution x to recover the solution xw for wells and applying /// xw to update Well State virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, - WellState& well_state) const; + WellState& well_state) const override; /// computing the well potentials for group control virtual void computeWellPotentials(const Simulator& ebosSimulator, const WellState& well_state, - std::vector& well_potentials); + std::vector& well_potentials) override; - virtual void updatePrimaryVariables(const WellState& well_state) const; + virtual void updatePrimaryVariables(const WellState& well_state) const override; - virtual void solveEqAndUpdateWellState(WellState& well_state); // const? + virtual void solveEqAndUpdateWellState(WellState& well_state) override; // const? virtual void calculateExplicitQuantities(const Simulator& ebosSimulator, - const WellState& well_state); // should be const? + const WellState& well_state) override; // should be const? /// number of segments for this well /// int number_of_segments_; @@ -152,6 +152,17 @@ namespace Opm int numberOfPerforations() const; + void addCellRates(RateVector& rates, int cellIdx) const override + { + for (int perfIdx = 0; perfIdx < number_of_perforations_; ++perfIdx) { + if (Base::cells()[perfIdx] == cellIdx) { + for (int i = 0; i < RateVector::dimension; ++i) { + rates[i] += connectionRates_[perfIdx][i]; + } + } + } + } + protected: int number_segments_; @@ -253,6 +264,8 @@ namespace Opm std::vector segment_depth_diffs_; + std::vector connectionRates_; + void initMatrixAndVectors(const int num_cells) const; // protected functions @@ -337,14 +350,13 @@ namespace Opm bool accelerationalPressureLossConsidered() const; // TODO: try to make ebosSimulator const, as it should be - void iterateWellEquations(Simulator& ebosSimulator, + void iterateWellEquations(const Simulator& ebosSimulator, const double dt, WellState& well_state); - void assembleWellEqWithoutIteration(Simulator& ebosSimulator, + void assembleWellEqWithoutIteration(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells); + WellState& well_state); }; } diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index 16e0d27a7..74bf8af38 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -52,6 +52,10 @@ namespace Opm if (has_polymer) { OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet"); } + + if (Base::has_energy) { + OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet"); + } // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. // for other facilities needed but not available from parser, we need to process them here @@ -110,6 +114,8 @@ namespace Opm { Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells); + connectionRates_.resize(number_of_perforations_); + // TODO: for StandardWell, we need to update the perf depth here using depth_arg. // for MultisegmentWell, it is much more complicated. // It can be specified directly, it can be calculated from the segment depth, @@ -228,10 +234,9 @@ namespace Opm template void MultisegmentWell:: - assembleWellEq(Simulator& ebosSimulator, + assembleWellEq(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells) + WellState& well_state) { const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_; @@ -239,7 +244,7 @@ namespace Opm iterateWellEquations(ebosSimulator, dt, well_state); } - assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells); + assembleWellEqWithoutIteration(ebosSimulator, dt, well_state); } @@ -1719,7 +1724,7 @@ namespace Opm template void MultisegmentWell:: - iterateWellEquations(Simulator& ebosSimulator, + iterateWellEquations(const Simulator& ebosSimulator, const double dt, WellState& well_state) { @@ -1732,7 +1737,7 @@ namespace Opm int it = 0; for (; it < max_iter_number; ++it) { - assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true); + assembleWellEqWithoutIteration(ebosSimulator, dt, well_state); const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); @@ -1762,19 +1767,16 @@ namespace Opm template void MultisegmentWell:: - assembleWellEqWithoutIteration(Simulator& ebosSimulator, + assembleWellEqWithoutIteration(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells) + WellState& well_state) { // calculate the fluid properties needed. computeSegmentFluidProperties(ebosSimulator); // clear all entries - if (!only_wells) { - duneB_ = 0.0; - duneC_ = 0.0; - } + duneB_ = 0.0; + duneC_ = 0.0; duneD_ = 0.0; resWell_ = 0.0; @@ -1784,9 +1786,6 @@ namespace Opm // // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same. - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); - const bool allow_cf = getAllowCrossFlow(); const int nseg = numberOfSegments(); @@ -1835,31 +1834,24 @@ namespace Opm // the cq_s entering mass balance equations need to consider the efficiency factors. const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_; - if (!only_wells) { - // subtract sum of component fluxes in the reservoir equation. - // need to consider the efficiency factor - ebosResid[cell_idx][comp_idx] -= cq_s_effective.value(); - } + connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective); // subtract sum of phase fluxes in the well equations. resWell_[seg][comp_idx] -= cq_s_effective.value(); // assemble the jacobians for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix - } + + // also need to consider the efficiency factor when manipulating the jacobians. + duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix + // the index name for the D should be eq_idx / pv_idx duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq); } for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx); - duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx); - } + // also need to consider the efficiency factor when manipulating the jacobians. + duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx); } } // TODO: we should save the perforation pressure and preforation rates? diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 854292af7..c1571fbce 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -188,18 +188,13 @@ public: SimulatorReport report; SimulatorReport stepReport; - WellModel wellModel(ebosSimulator_, modelParam_, terminalOutput_); if (isRestart()) { - wellModel.initFromRestartFile(*restartValues); + wellModel_().initFromRestartFile(*restartValues); } - if (modelParam_.matrix_add_well_contributions_ || - modelParam_.preconditioner_add_well_contributions_) - { - ebosSimulator_.model().clearAuxiliaryModules(); - wellAuxMod_.reset(new WellConnectionAuxiliaryModule(schedule(), grid())); - ebosSimulator_.model().addAuxiliaryModule(wellAuxMod_.get()); - } + // beginReportStep(...) wants to know when we are at the + // beginning of a restart + bool firstRestartStep = isRestart(); // Main simulation loop. while (!timer.done()) { @@ -210,30 +205,25 @@ public: OpmLog::debug(ss.str()); } - // Run a multiple steps of the solver depending on the time step control. - solverTimer.start(); - //ebosSimulator_.setEpisodeIndex(timer.reportStepNum()); - wellModel.beginReportStep(timer.currentStepNum()); - - auto solver = createSolver(wellModel); - // write the inital state at the report stage if (timer.initialStep()) { Dune::Timer perfTimer; perfTimer.start(); - // No per cell data is written for initial step, but will be - // for subsequent steps, when we have started simulating - auto localWellData = wellModel.wellState().report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid())); - ebosSimulator_.problem().writeOutput(localWellData, - timer.simulationTimeElapsed(), - /*isSubstep=*/false, - totalTimer.secsSinceStart(), - /*nextStepSize=*/-1.0); + wellModel_().beginReportStep(timer.currentStepNum()); + ebosSimulator_.problem().writeOutput(false); report.output_write_time += perfTimer.stop(); } + // Run a multiple steps of the solver depending on the time step control. + solverTimer.start(); + + auto solver = createSolver(wellModel_()); + + solver->model().beginReportStep(firstRestartStep); + firstRestartStep = false; + if (terminalOutput_) { std::ostringstream stepMsg; boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y"); @@ -246,8 +236,6 @@ public: OpmLog::info(stepMsg.str()); } - solver->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 // @@ -282,7 +270,6 @@ public: } solver->model().endReportStep(); - wellModel.endReportStep(); // take time that was used to solve system for this reportStep solverTimer.stop(); @@ -305,13 +292,8 @@ public: Dune::Timer perfTimer; perfTimer.start(); const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0; - - auto localWellData = wellModel.wellState().report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid())); - ebosSimulator_.problem().writeOutput(localWellData, - timer.simulationTimeElapsed(), - /*isSubstep=*/false, - totalTimer.secsSinceStart(), - nextstep); + ebosSimulator_.problem().setNextTimeStepSize(nextstep); + ebosSimulator_.problem().writeOutput(false); report.output_write_time += perfTimer.stop(); if (terminalOutput_) { @@ -379,6 +361,12 @@ protected: return initconfig.restartRequested(); } + WellModel& wellModel_() + { return ebosSimulator_.problem().wellModel(); } + + const WellModel& wellModel_() const + { return ebosSimulator_.problem().wellModel(); } + // Data. Simulator& ebosSimulator_; std::unique_ptr> wellAuxMod_; diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index e20e5dfd3..da83ef86a 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -37,6 +37,8 @@ namespace Opm public: typedef WellInterface Base; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; + // TODO: some functions working with AD variables handles only with values (double) without // dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value. // And also, it can also be beneficial to make these functions hanle different types of AD variables. @@ -129,52 +131,63 @@ namespace Opm virtual void init(const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, const double gravity_arg, - const int num_cells); + const int num_cells) override; - virtual void initPrimaryVariablesEvaluation() const; + virtual void initPrimaryVariablesEvaluation() const override; - virtual void assembleWellEq(Simulator& ebosSimulator, + virtual void assembleWellEq(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells); + WellState& well_state) override; /// updating the well state based the control mode specified with current // TODO: later will check wheter we need current - virtual void updateWellStateWithTarget(WellState& well_state) const; + virtual void updateWellStateWithTarget(WellState& well_state) const override; /// check whether the well equations get converged for this well - virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const; + virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const override; /// Ax = Ax - C D^-1 B x - virtual void apply(const BVector& x, BVector& Ax) const; + virtual void apply(const BVector& x, BVector& Ax) const override; /// r = r - C D^-1 Rw - virtual void apply(BVector& r) const; + virtual void apply(BVector& r) const override; /// using the solution x to recover the solution xw for wells and applying /// xw to update Well State virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, - WellState& well_state) const; + WellState& well_state) const override; /// computing the well potentials for group control virtual void computeWellPotentials(const Simulator& ebosSimulator, const WellState& well_state, - std::vector& well_potentials) /* const */; + std::vector& well_potentials) /* const */ override; - virtual void updatePrimaryVariables(const WellState& well_state) const; + virtual void updatePrimaryVariables(const WellState& well_state) const override; - virtual void solveEqAndUpdateWellState(WellState& well_state); + virtual void solveEqAndUpdateWellState(WellState& well_state) override; virtual void calculateExplicitQuantities(const Simulator& ebosSimulator, - const WellState& well_state); // should be const? + const WellState& well_state) override; // should be const? - virtual void addWellContributions(Mat& mat) const; + virtual void addWellContributions(Mat& mat) const override; /// \brief Wether the Jacobian will also have well contributions in it. - virtual bool jacobianContainsWellContributions() const + virtual bool jacobianContainsWellContributions() const override { return param_.matrix_add_well_contributions_; } + + void addCellRates(RateVector& rates, int cellIdx) const override + { + for (int perfIdx = 0; perfIdx < number_of_perforations_; ++perfIdx) { + if (Base::cells()[perfIdx] == cellIdx) { + for (int i = 0; i < RateVector::dimension; ++i) { + rates[i] += connectionRates_[perfIdx][i]; + } + } + } + } + protected: // protected functions from the Base class @@ -228,6 +241,8 @@ namespace Opm // diagonal matrix for the well DiagMatWell invDuneD_; + std::vector connectionRates_; + // several vector used in the matrix calculation mutable BVectorWell Bx_; mutable BVectorWell invDrw_; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 950e937fc..a94c56612 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -58,6 +58,8 @@ namespace Opm { Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells); + connectionRates_.resize(number_of_perforations_); + perf_depth_.resize(number_of_perforations_, 0.); for (int perf = 0; perf < number_of_perforations_; ++perf) { const int cell_idx = well_cells_[perf]; @@ -241,6 +243,7 @@ namespace Opm StandardWell:: wellSurfaceVolumeFraction(const int compIdx) const { + EvalWell sum_volume_fraction_scaled = 0.; for (int idx = 0; idx < num_components_; ++idx) { sum_volume_fraction_scaled += wellVolumeFractionScaled(idx); @@ -432,24 +435,18 @@ namespace Opm template void StandardWell:: - assembleWellEq(Simulator& ebosSimulator, + assembleWellEq(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells) + WellState& well_state) { const int np = number_of_phases_; // clear all entries - if (!only_wells) { - duneB_ = 0.0; - duneC_ = 0.0; - } + duneB_ = 0.0; + duneC_ = 0.0; invDuneD_ = 0.0; resWell_ = 0.0; - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); - // TODO: it probably can be static member for StandardWell const double volume = 0.002831684659200; // 0.1 cu ft; @@ -483,34 +480,28 @@ namespace Opm well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; } + if (has_energy) { + connectionRates_[perf][contiEnergyEqIdx] = 0.0; + } + for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { // the cq_s entering mass balance equations need to consider the efficiency factors. const EvalWell cq_s_effective = cq_s[componentIdx] * well_efficiency_factor_; - if (!only_wells) { - // subtract sum of component fluxes in the reservoir equation. - // need to consider the efficiency factor - ebosResid[cell_idx][componentIdx] -= cq_s_effective.value(); - } + connectionRates_[perf][componentIdx] = Base::restrictEval(cq_s_effective); // subtract sum of phase fluxes in the well equations. resWell_[0][componentIdx] -= cq_s_effective.value(); // assemble the jacobians for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix - } + // also need to consider the efficiency factor when manipulating the jacobians. + duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq); } for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx); - duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx); - } + duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx); } // Store the perforation phase flux for later usage. @@ -574,15 +565,7 @@ namespace Opm } // compute the thermal flux cq_r_thermal *= extendEval(fs.enthalpy(phaseIdx)) * extendEval(fs.density(phaseIdx)); - // scale the flux by the scaling factor for the energy equation - cq_r_thermal *= GET_PROP_VALUE(TypeTag, BlackOilEnergyScalingFactor); - - if (!only_wells) { - for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { - ebosJac[cell_idx][cell_idx][contiEnergyEqIdx][pvIdx] -= cq_r_thermal.derivative(pvIdx); - } - ebosResid[cell_idx][contiEnergyEqIdx] -= cq_r_thermal.value(); - } + connectionRates_[perf][contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal); } } @@ -595,12 +578,7 @@ namespace Opm } else { cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); } - if (!only_wells) { - for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { - ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx); - } - ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value(); - } + connectionRates_[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly); } // Store the perforation pressure for later usage. @@ -631,9 +609,16 @@ namespace Opm assembleControlEq(); // do the local inversion of D. - // we do this manually with invertMatrix to always get our - // specializations in for 3x3 and 4x4 matrices. - Dune::ISTLUtility::invertMatrix(invDuneD_[0][0]); + try + { + Dune::ISTLUtility::invertMatrix(invDuneD_[0][0]); + } + catch( ... ) + { + OPM_THROW(Opm::NumericalIssue,"Error when inverting local well equations for well " + name()); + } + + } diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index 4fcf5e4e4..1a7aa3144 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -82,6 +82,7 @@ namespace Opm typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; static const int numEq = Indices::numEq; typedef double Scalar; @@ -121,7 +122,7 @@ namespace Opm const int indexOfWell() const; /// Well cells. - const std::vector& cells() {return well_cells_; } + const std::vector& cells() const {return well_cells_; } /// Well type, INJECTOR or PRODUCER. WellType wellType() const; @@ -142,20 +143,19 @@ namespace Opm virtual void solveEqAndUpdateWellState(WellState& well_state) = 0; - virtual void assembleWellEq(Simulator& ebosSimulator, + virtual void assembleWellEq(const Simulator& ebosSimulator, const double dt, - WellState& well_state, - bool only_wells) = 0; + WellState& well_state) = 0; void updateWellTestState(const WellState& well_state, const double& simulationTime, const bool& writeMessageToOPMLog, - WellTestState& wellTestState) const; - + WellTestState& wellTestState + ) const; void setWellEfficiencyFactor(const double efficiency_factor); - void computeRepRadiusPerfLength(const Grid& grid, const std::map& cartesian_to_compressed); + void computeRepRadiusPerfLength(const Grid& grid, const std::vector& cartesian_to_compressed); /// using the solution x to recover the solution xw for wells and applying /// xw to update Well State @@ -196,6 +196,20 @@ namespace Opm virtual void addWellContributions(Mat&) const {} + virtual void addCellRates(RateVector& rates, int cellIdx) const + {} + + template + Eval restrictEval(const EvalWell& in) const + { + Eval out = 0.0; + out.setValue(in.value()); + for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { + out.setDerivative(eqIdx, in.derivative(eqIdx)); + } + return out; + } + void closeCompletions(WellTestState& wellTestState); const Well* wellEcl() const; diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index e111c63a7..f46f60241 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -875,7 +875,7 @@ namespace Opm void WellInterface:: computeRepRadiusPerfLength(const Grid& grid, - const std::map& cartesian_to_compressed) + const std::vector& cartesian_to_compressed) { const int* cart_dims = Opm::UgGridHelpers::cartDims(grid); auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid); @@ -902,12 +902,12 @@ namespace Opm const int* cpgdim = cart_dims; const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); - const std::map::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); - if (cgit == cartesian_to_compressed.end()) { + const int cell = cartesian_to_compressed[cart_grid_indx]; + + if (cell < 0) { OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' << k << " not found in grid (well = " << name() << ')'); } - const int cell = cgit->second; { double radius = connection.rw(); @@ -1023,7 +1023,7 @@ namespace Opm bool converged; WellState well_state0 = well_state; do { - assembleWellEq(ebosSimulator, dt, well_state, true); + assembleWellEq(ebosSimulator, dt, well_state); auto report = getWellConvergence(B_avg); converged = report.converged(); @@ -1047,6 +1047,7 @@ namespace Opm if ( terminal_output ) { OpmLog::debug("WellTest: Well equation for well" +name() + " solution failed in getting converged with " + std::to_string(it) + " iterations"); } + well_state = well_state0; } } diff --git a/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp b/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp index c982003d3..bd59839ee 100644 --- a/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp @@ -203,7 +203,6 @@ namespace Opm { auto& ebosSimulator = solver.model().ebosSimulator(); auto& ebosProblem = ebosSimulator.problem(); - auto phaseUsage = Opm::phaseUsageFromDeck(ebosSimulator.vanguard().eclState()); // create adaptive step timer with previously used sub step size AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_); @@ -316,13 +315,7 @@ namespace Opm { Opm::time::StopWatch perfTimer; perfTimer.start(); - // The writeOutput expects a local data::solution vector and a local data::well vector. - auto localWellData = solver.model().wellModel().wellState().report(phaseUsage, Opm::UgGridHelpers::globalCell(ebosSimulator.vanguard().grid())); - ebosProblem.writeOutput(localWellData, - substepTimer.simulationTimeElapsed(), - /*isSubstep=*/true, - substepReport.total_time, - /*nextStepSize=*/-1.0); + ebosProblem.writeOutput(/*isSubStep=*/true); report.output_write_time += perfTimer.secsSinceStart(); } @@ -361,6 +354,7 @@ namespace Opm { ++restarts; } + ebosProblem.setNextTimeStepSize(substepTimer.currentStepLength()); }