From 4ecd6ca64a7364df62a001f5778ccb81f929af29 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 12 Sep 2016 23:28:44 +0200 Subject: [PATCH] fix some serious screw-ups almost all of them were caused by recent changes in the master branch: - there were methods added which depend on the types `V` and `DataBlock`. these do not make much sense in the context of the frankenstein simulator. Also, these types are defined globally for the whole Opm namespace in `BlackoilModelBase_impl.hpp` (which should be prosecuted as a fellony IMO)! Besides this, their names are useless; 'V' is the letter which comes after `U` in the alphabet and when it comes to computers basically everything can be seen as a chunk of data (i.e., a `DataBlock`). - it seems like the new and shiny dense-AD based well model was never compiled with assertations enabled, at least some asserts referenced non-existing variables. - the recent output-related API changes were pretty unfortunate because they had the effect of tying the (sub-optimal, IMO) internal structure of the model even closer to the output code: as far as I can see, `rq` does only make sense if the model works *exactly* like BlackoilModelBase and friends. (for flow_ebos, this could be replicated, but first it would be another unnecessary conversion step and second, most of the quantities in `rq` are of type `ADB` and much of the "frankenstein" excercise is devoted to getting rid of these.) I thus reverted back to an old version of the output code and created a `frankenstein` branch in my personal `opm-output` github fork. --- CMakeLists_files.cmake | 1 + opm/autodiff/BlackoilModelBase.hpp | 4 +- opm/autodiff/BlackoilModelEbos.hpp | 24 +- opm/autodiff/BlackoilPressureModel.hpp | 1 - opm/autodiff/FlowMain.hpp | 13 +- .../NewtonIterationBlackoilInterleaved.cpp | 13 +- opm/autodiff/NewtonIterationUtilities.cpp | 12 - opm/autodiff/NonlinearSolver.hpp | 6 +- opm/autodiff/NonlinearSolver_impl.hpp | 13 +- opm/autodiff/SimFIBODetails.hpp | 151 ++++++ opm/autodiff/SimulatorBase.hpp | 4 +- opm/autodiff/SimulatorBase_impl.hpp | 1 + .../SimulatorFullyImplicitBlackoilEbos.hpp | 37 +- ...mulatorFullyImplicitBlackoilOutputEbos.cpp | 282 ++++++++++ ...mulatorFullyImplicitBlackoilOutputEbos.hpp | 500 ++++++++++++++++++ opm/autodiff/StandardWellsDense.hpp | 6 +- .../WellStateFullyImplicitBlackoil.hpp | 1 - ...FullyImplicitCompressiblePolymerSolver.cpp | 2 +- 18 files changed, 998 insertions(+), 73 deletions(-) create mode 100644 opm/autodiff/SimFIBODetails.hpp create mode 100644 opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp create mode 100644 opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index d8025c6b0..8b1edaa84 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 diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 1f990212f..37c82f2b2 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -27,13 +27,13 @@ #include #include -#include -#include #include #include #include #include #include +#include +#include #include #include diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index edf51702d..88dcd3da7 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -42,6 +42,7 @@ #include #include #include +#include #include #include @@ -52,6 +53,8 @@ #include #include #include +#include +#include #include #include #include @@ -221,7 +224,6 @@ namespace Opm { ReservoirState& reservoir_state, WellState& well_state) { - const double dt = timer.currentStepLength(); 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. @@ -302,13 +304,6 @@ namespace Opm { OPM_THROW(Opm::NumericalProblem,"no convergence"); } - typedef double Scalar; - typedef Dune::FieldVector VectorBlockType; - typedef Dune::FieldMatrix MatrixBlockType; - typedef Dune::BCRSMatrix Mat; - typedef Dune::BlockVector BVector; - - typedef Dune::MatrixAdapter Operator; auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); @@ -927,6 +922,16 @@ namespace Opm { } } + std::vector > + computeFluidInPlace(const ReservoirState& x, + const std::vector& fipnum) const + { + OPM_THROW(std::logic_error, + "computeFluidInPlace() not implemented by BlackoilModelEbos!"); + } + + const Simulator& ebosSimulator() const + { return ebosSimulator_; } protected: @@ -1244,8 +1249,7 @@ namespace Opm { { const int numPhases = wells().number_of_phases; const int numCells = ebosJac.N(); - const int cols = ebosJac.M(); - assert( numCells == cols ); + assert( numCells == static_cast(ebosJac.M()) ); // write the right-hand-side values from the ebosJac into the objects // allocated above. diff --git a/opm/autodiff/BlackoilPressureModel.hpp b/opm/autodiff/BlackoilPressureModel.hpp index 05847ee8d..716705fb6 100644 --- a/opm/autodiff/BlackoilPressureModel.hpp +++ b/opm/autodiff/BlackoilPressureModel.hpp @@ -26,7 +26,6 @@ #include #include #include -#include #include #include diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 61fbc984d..58d6f6dae 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -173,6 +173,7 @@ namespace Opm typedef BlackoilPropsAdFromDeck FluidProps; typedef FluidProps::MaterialLawManager MaterialLawManager; typedef typename Simulator::ReservoirState ReservoirState; + typedef typename Simulator::OutputWriter OutputWriter; // ------------ Data members ------------ @@ -207,7 +208,7 @@ namespace Opm // distributeData() boost::any parallel_information_; // setupOutputWriter() - std::unique_ptr output_writer_; + std::unique_ptr output_writer_; // setupLinearSolver std::unique_ptr fis_solver_; // createSimulator() @@ -673,11 +674,11 @@ 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_, - Opm::phaseUsageFromDeck(deck_), - fluidprops_->permeability())); + output_writer_.reset(new OutputWriter(grid_init_->grid(), + param_, + eclipse_state_, + Opm::phaseUsageFromDeck(deck_), + fluidprops_->permeability())); } diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 1d96d6ea3..1ab7f4906 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -390,7 +390,7 @@ namespace Opm //const int np = residual.material_balance_eq.size(); assert( np == int(residual.material_balance_eq.size()) ); std::vector eqs; - eqs.reserve(np + 1); + eqs.reserve(np + 2); for (int phase = 0; phase < np; ++phase) { eqs.push_back(residual.material_balance_eq[phase]); } @@ -401,14 +401,14 @@ namespace Opm if( hasWells ) { eqs.push_back(residual.well_flux_eq); - //eqs.push_back(residual.well_eq); + eqs.push_back(residual.well_eq); // Eliminate the well-related unknowns, and corresponding equations. - elim_eqs.reserve(1); + elim_eqs.reserve(2); elim_eqs.push_back(eqs[np]); eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns. - //elim_eqs.push_back(eqs[np]); - //eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns. + elim_eqs.push_back(eqs[np]); + eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns. assert(int(eqs.size()) == np); } @@ -500,7 +500,7 @@ namespace Opm if ( hasWells ) { // Compute full solution using the eliminated equations. // Recovery in inverse order of elimination. - //dx = recoverVariable(elim_eqs[1], dx, np); + dx = recoverVariable(elim_eqs[1], dx, np); dx = recoverVariable(elim_eqs[0], dx, np); } return dx; @@ -575,7 +575,6 @@ namespace Opm computePressureIncrement(const LinearisedBlackoilResidual& residual) { typedef LinearisedBlackoilResidual::ADB ADB; - typedef ADB::V V; // Build the vector of equations (should be just a single material balance equation // in which the pressure equation is stored). diff --git a/opm/autodiff/NewtonIterationUtilities.cpp b/opm/autodiff/NewtonIterationUtilities.cpp index 8f68b5b48..05bf0dc34 100644 --- a/opm/autodiff/NewtonIterationUtilities.cpp +++ b/opm/autodiff/NewtonIterationUtilities.cpp @@ -85,16 +85,9 @@ namespace Opm for (int eq = 0; eq < num_eq; ++eq) { jacs[eq].reserve(num_eq - 1); const std::vector& Je = eqs[eq].derivative(); - Sp Bb; const M& B = Je[n]; - B.toSparse(Bb); - //std::cout << "B eigen" << std::endl; - //std::cout << Bb << std::endl; - // Update right hand side. vals[eq] = eqs[eq].value().matrix() - B * Dibn; - //std::cout << "vals " << eq << std::endl; - //std::cout << vals[eq][0] << std::endl; } for (int var = 0; var < num_eq; ++var) { if (var == n) { @@ -156,11 +149,6 @@ namespace Opm V equation_value = equation.value(); ADB eq_coll = collapseJacs(ADB::function(std::move(equation_value), std::move(C_jacs))); const M& C = eq_coll.derivative()[0]; - typedef Eigen::SparseMatrix Sp; - Sp Cc; - C.toSparse(Cc); - //std::cout << "C eigen" << std::endl; - //std::cout << Cc < Sp; diff --git a/opm/autodiff/NonlinearSolver.hpp b/opm/autodiff/NonlinearSolver.hpp index 09af664b9..78a8e7a33 100644 --- a/opm/autodiff/NonlinearSolver.hpp +++ b/opm/autodiff/NonlinearSolver.hpp @@ -131,7 +131,11 @@ 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); + } + /// Reference to physical model. const PhysicalModel& model() const; diff --git a/opm/autodiff/NonlinearSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp index d1dd65ab8..925152700 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:: @@ -119,6 +112,7 @@ namespace Opm } + template int NonlinearSolver:: @@ -177,6 +171,7 @@ namespace Opm } + template void NonlinearSolver::SolverParameters:: reset() 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.hpp b/opm/autodiff/SimulatorBase.hpp index 493d6769c..35f315e99 100644 --- a/opm/autodiff/SimulatorBase.hpp +++ b/opm/autodiff/SimulatorBase.hpp @@ -21,14 +21,12 @@ #ifndef OPM_SIMULATORBASE_HEADER_INCLUDED #define OPM_SIMULATORBASE_HEADER_INCLUDED -#include #include -#include -#include #include #include #include +#include #include #include #include diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index 96c16fcc5..e3fe40a31 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -26,6 +26,7 @@ #include #include #include +#include namespace Opm { diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 98cfc7fb4..674a6a199 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -18,17 +18,26 @@ along with OPM. If not, see . */ -#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED -#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED +#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 +#include namespace Opm { @@ -47,7 +56,7 @@ public: typedef WellStateFullyImplicitBlackoil WellState; typedef BlackoilState ReservoirState; - typedef BlackoilOutputWriter OutputWriter; + typedef BlackoilOutputWriterEbos OutputWriter; typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; @@ -90,7 +99,7 @@ public: const bool has_disgas, const bool has_vapoil, std::shared_ptr eclipse_state, - BlackoilOutputWriter& output_writer, + BlackoilOutputWriterEbos& output_writer, const std::vector& threshold_pressures_by_face) : ebosSimulator_(ebosSimulator), param_(param), @@ -166,15 +175,11 @@ public: adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) ); } - - - output_writer_.writeInit( geo_.simProps(grid()) , geo_.nonCartesianConnections( ) ); - 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) ); +// const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); // output_writer_.restore( timer, // state, @@ -221,11 +226,6 @@ public: // give the polymer and surfactant simulators the chance to do their stuff handleAdditionalWellInflow(timer, wells_manager, well_state, wells); - // write the inital state at the report stage - if (timer.initialStep()) { - output_writer_.writeTimeStep( timer, state, well_state ); - } - // Compute reservoir volumes for RESV controls. computeRESV(timer.currentStepNum(), wells, state, well_state); @@ -237,6 +237,11 @@ public: auto solver = createSolver(well_model); + // write the inital state at the report stage + if (timer.initialStep()) { + output_writer_.writeTimeStep( timer, state, well_state, solver->model() ); + } + if( terminal_output_ ) { std::ostringstream step_msg; @@ -326,7 +331,7 @@ public: ++timer; // write simulation state at the report stage - output_writer_.writeTimeStep( timer, state, well_state ); + output_writer_.writeTimeStep( timer, state, well_state, solver->model() ); prev_well_state = well_state; // The well potentials are only computed if they are needed diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp new file mode 100644 index 000000000..b4a38e187 --- /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 + +#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_; + std::vector simProps_; + const bool substep_; + + explicit WriterCallEbos( BlackoilOutputWriterEbos& writer, + const SimulatorTimerInterface& timer, + const SimulationDataContainer& state, + const WellState& wellState, + const std::vector& 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) + { + std::vector noCellProperties; + writeTimeStepWithCellProperties(timer, localState, localWellState, noCellProperties, substep); + } + + + + + + void + BlackoilOutputWriterEbos:: + writeTimeStepWithCellProperties( + const SimulatorTimerInterface& timer, + const SimulationDataContainer& localState, + const WellState& localWellState, + const std::vector& cellData, + 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, 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, cellData, substep ) ); + } + else { + // just write the data to disk + writeTimeStepSerial( timer, state, wellState, cellData, substep ); + } + } + } + + + + void + BlackoilOutputWriterEbos:: + writeTimeStepSerial(const SimulatorTimerInterface& timer, + const SimulationDataContainer& state, + const WellState& wellState, + const std::vector& simProps, + 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 { + eclWriter_->writeTimeStep(timer.reportStepNum(), + substep, + timer.simulationTimeElapsed(), + simToSolution( state, phaseUsage_ ), + wellState.report(), + simProps); + } + } + + // 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, + WellStateFullyImplicitBlackoil& wellState, + const std::string& filename, + const int desiredResportStep ) + { + std::ifstream restorefile( filename.c_str() ); + if( restorefile ) + { + std::cout << "============================================================================"<getInitConfig(); + return initconfig.restartRequested(); + } +} diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp new file mode 100644 index 000000000..11053751e --- /dev/null +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp @@ -0,0 +1,500 @@ +/* + Copyright (c) 2014 SINTEF ICT, Applied Mathematics. + Copyright (c) 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_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, + Opm::EclipseStateConstPtr eclipseState, + 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 std::vector& simProps, + 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 std::vector& 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, + WellStateFullyImplicitBlackoil& wellState, + const std::string& filename, + const int desiredReportStep); + + + template + void initFromRestartFile(const PhaseUsage& phaseusage, + const double* permeability, + const Grid& grid, + SimulationDataContainer& simulatorstate, + WellStateFullyImplicitBlackoil& 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< EclipseWriter > eclWriter_; + EclipseStateConstPtr eclipseState_; + + std::unique_ptr< ThreadHandle > asyncOutput_; + }; + + + ////////////////////////////////////////////////////////////// + // + // Implementation + // + ////////////////////////////////////////////////////////////// + template + inline + BlackoilOutputWriterEbos:: + BlackoilOutputWriterEbos(const Grid& grid, + const parameter::ParameterGroup& param, + Opm::EclipseStateConstPtr eclipseState, + const Opm::PhaseUsage &phaseUsage, + const double* permeability ) + : output_( param.getDefault("output", true) ), + parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, permeability ) : 0 ), + outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ), + output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ), + lastBackupReportStep_( -1 ), + phaseUsage_( phaseUsage ), + eclWriter_( output_ && parallelOutput_->isIORank() && + param.getDefault("output_ecl", true) ? + new EclipseWriter(eclipseState,UgGridHelpers::createEclipseGrid( grid , *eclipseState->getInputGrid())) + : 0 ), + eclipseState_(eclipseState), + asyncOutput_() + { + // For output. + if (output_ && parallelOutput_->isIORank() ) { + // 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, + WellStateFullyImplicitBlackoil& 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); + + 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, wellstate ); + } + + + + + + namespace detail { + template + std::vector getCellDataEbos( + const Opm::PhaseUsage& phaseUsage, + const Model& model, + const RestartConfig& restartConfig, + const int reportStepNum) + { + typedef typename Model::FluidSystem FluidSystem; + + std::vector simProps; + + //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); + } + + //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]; + + const auto& ebosModel = model.ebosSimulator().model(); + + // extract everything which can possibly be written to disk + int numCells = ebosModel.numGridDof(); + 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); + + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + const auto& intQuants = *ebosModel.cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value; + bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value; + bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value; + + rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value; + rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value; + rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value; + + muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value; + muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value; + muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value; + + krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value; + krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value; + krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value; + + Rs[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs, + FluidSystem::oilPhaseIdx, + intQuants.pvtRegionIndex(), + /*maxOilSaturation=*/1.0).value; + Rv[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs, + FluidSystem::gasPhaseIdx, + intQuants.pvtRegionIndex(), + /*maxOilSaturation=*/1.0).value; + } + + /** + * Formation volume factors for water, oil, gas + */ + if (aqua_active && outKeywords["BW"] > 0) { + outKeywords["BW"] = 0; + simProps.emplace_back(data::CellData{ + "1OVERBW", + Opm::UnitSystem::measure::water_inverse_formation_volume_factor, + std::move(bWater)}); + } + if (liquid_active && outKeywords["BO"] > 0) { + outKeywords["BO"] = 0; + simProps.emplace_back(data::CellData{ + "1OVERBO", + Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, + std::move(bOil)}); + } + if (vapour_active && outKeywords["BG"] > 0) { + outKeywords["BG"] = 0; + simProps.emplace_back(data::CellData{ + "1OVERBG", + Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, + std::move(bGas)}); + } + + /** + * Densities for water, oil gas + */ + if (outKeywords["DEN"] > 0) { + outKeywords["DEN"] = 0; + if (aqua_active) { + simProps.emplace_back(data::CellData{ + "WAT_DEN", + Opm::UnitSystem::measure::density, + std::move(rhoWater)}); + } + if (liquid_active) { + simProps.emplace_back(data::CellData{ + "OIL_DEN", + Opm::UnitSystem::measure::density, + std::move(rhoOil)}); + } + if (vapour_active) { + simProps.emplace_back(data::CellData{ + "GAS_DEN", + Opm::UnitSystem::measure::density, + std::move(rhoGas)}); + } + } + + /** + * Viscosities for water, oil gas + */ + if (outKeywords["VISC"] > 0) { + outKeywords["VISC"] = 0; + if (aqua_active) { + simProps.emplace_back(data::CellData{ + "WAT_VISC", + Opm::UnitSystem::measure::viscosity, + std::move(muWater)}); + } + if (liquid_active) { + simProps.emplace_back(data::CellData{ + "OIL_VISC", + Opm::UnitSystem::measure::viscosity, + std::move(muOil)}); + } + if (vapour_active) { + simProps.emplace_back(data::CellData{ + "GAS_VISC", + Opm::UnitSystem::measure::viscosity, + std::move(muGas)}); + } + } + + /** + * Relative permeabilities for water, oil, gas + */ + if (aqua_active && outKeywords["KRW"] > 0) { + outKeywords["KRW"] = 0; + simProps.emplace_back(data::CellData{ + "WATKR", + Opm::UnitSystem::measure::permeability, + std::move(krWater)}); + } + if (liquid_active && outKeywords["KRO"] > 0) { + outKeywords["KRO"] = 0; + simProps.emplace_back(data::CellData{ + "OILKR", + Opm::UnitSystem::measure::permeability, + std::move(krOil)}); + } + if (vapour_active && outKeywords["KRG"] > 0) { + outKeywords["KRG"] = 0; + simProps.emplace_back(data::CellData{ + "GASKR", + Opm::UnitSystem::measure::permeability, + std::move(krGas)}); + } + + /** + * Vaporized and dissolved gas/oil ratio + */ + if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) { + outKeywords["RSSAT"] = 0; + simProps.emplace_back(data::CellData{ + "RSSAT", + Opm::UnitSystem::measure::gas_oil_ratio, + std::move(Rs)}); + } + if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) { + outKeywords["RVSAT"] = 0; + simProps.emplace_back(data::CellData{ + "RVSAT", + Opm::UnitSystem::measure::oil_gas_ratio, + std::move(Rv)}); + } + + + /** + * Bubble point and dew point pressures + */ + if (vapour_active && liquid_active && outKeywords["PBPD"] > 0) { + 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 + 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); + } + } + + return simProps; + } + } + + + + + template + inline void + BlackoilOutputWriterEbos:: + writeTimeStep(const SimulatorTimerInterface& timer, + const SimulationDataContainer& localState, + const WellState& localWellState, + const Model& physicalModel, + bool substep) + { + const RestartConfig& restartConfig = eclipseState_->getRestartConfig(); + const int reportStepNum = timer.reportStepNum(); + std::vector cellData = detail::getCellDataEbos( phaseUsage_, physicalModel, restartConfig, reportStepNum ); + writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep); + } +} +#endif diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index f0b3a075f..687e197e1 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -47,7 +47,7 @@ #include #include #include - +#include #include #include @@ -305,7 +305,7 @@ namespace Opm { } void addRhs(BVector& x, Mat& jac) const { - assert(x.size() == rhs.size()); + assert(x.size() == rhs_.size()); x += rhs_; // jac = A + duneA jac = matAdd( jac, duneA_ ); @@ -955,7 +955,6 @@ namespace Opm { dx_new_eigen(idx) = dx_new[w][flowPhaseToEbosCompIdx(p)]; } } - assert(dx.size() == total_residual_v.size()); updateWellState(dx_new_eigen.array(), well_state); updateWellControls(well_state); setWellVariables(well_state); @@ -1995,7 +1994,6 @@ namespace Opm { } EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const { - assert(fluid_.numPhases() == 3); const int nw = wells().number_of_wells; if (phaseIdx == Water) { return wellVariables_[nw + wellIdx]; diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 6de120ba8..f9e83cd5c 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -225,7 +225,6 @@ namespace Opm } } } - // perfPressures if( num_perf_old_well == num_perf_this_well ) { diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index ca3c1fa74..68147b5b0 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -205,8 +205,8 @@ namespace { PolymerBlackoilState& x , WellStateFullyImplicitBlackoilPolymer& xw) { - const double dt = timer.currentStepLength(); const std::vector& polymer_inflow = xw.polymerInflow(); + const double dt = timer.currentStepLength(); // Initial max concentration of this time step from PolymerBlackoilState. cmax_ = Eigen::Map(&x.getCellData( x.CMAX )[0], Opm::AutoDiffGrid::numCells(grid_));