OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #389 from andlaus/simplify_simulator

Browse Source 
Simplify simulator
This commit is contained in:
Atgeirr Flø Rasmussen 2015-05-29 15:46:38 +02:00
commit af9a5992a3
5 changed files with 331 additions and 289 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
opm/autodiff/NewtonSolver.hpp
View File

@ -23,6 +23,7 @@
#include <opm/autodiff/AutoDiffBlock.hpp> #include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <memory>
namespace Opm { namespace Opm {
@ -63,10 +64,14 @@ namespace Opm {
// --------- Public methods --------- // --------- Public methods ---------
/// Construct solver for a given model. /// Construct solver for a given model.
///
/// The model is a std::unique_ptr because the object to which model points to is
/// not allowed to be deleted as long as the NewtonSolver object exists.
///
/// \param[in] param parameters controlling nonlinear Newton process /// \param[in] param parameters controlling nonlinear Newton process
/// \param[in, out] model physical simulation model /// \param[in, out] model physical simulation model.
explicit NewtonSolver(const SolverParameters& param, explicit NewtonSolver(const SolverParameters& param,
PhysicalModel& model); std::unique_ptr<PhysicalModel> model);
/// Take a single forward step, after which the states will be modified /// Take a single forward step, after which the states will be modified
/// according to the physical model. /// according to the physical model.
@ -94,7 +99,7 @@ namespace Opm {
private: private:
// --------- Data members --------- // --------- Data members ---------
SolverParameters param_; SolverParameters param_;
PhysicalModel& model_; std::unique_ptr<PhysicalModel> model_;
unsigned int newtonIterations_; unsigned int newtonIterations_;
unsigned int linearIterations_; unsigned int linearIterations_;
unsigned int newtonIterationsLast_; unsigned int newtonIterationsLast_;

34
opm/autodiff/NewtonSolver_impl.hpp
View File

@ -29,9 +29,9 @@ namespace Opm
{ {
template <class PhysicalModel> template <class PhysicalModel>
NewtonSolver<PhysicalModel>::NewtonSolver(const SolverParameters& param, NewtonSolver<PhysicalModel>::NewtonSolver(const SolverParameters& param,
PhysicalModel& model) std::unique_ptr<PhysicalModel> model)
: param_(param), : param_(param),
model_(model), model_(std::move(model)),
newtonIterations_(0), newtonIterations_(0),
linearIterations_(0) linearIterations_(0)
{ {
@ -58,21 +58,21 @@ namespace Opm
WellState& well_state) WellState& well_state)
{ {
// Do model-specific once-per-step calculations. // Do model-specific once-per-step calculations.
model_.prepareStep(dt, reservoir_state, well_state); model_->prepareStep(dt, reservoir_state, well_state);
// For each iteration we store in a vector the norms of the residual of // 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. // the mass balance for each active phase, the well flux and the well equations.
std::vector<std::vector<double>> residual_norms_history; std::vector<std::vector<double>> residual_norms_history;
// Assemble residual and Jacobian, store residual norms. // Assemble residual and Jacobian, store residual norms.
model_.assemble(reservoir_state, well_state, true); model_->assemble(reservoir_state, well_state, true);
residual_norms_history.push_back(model_.computeResidualNorms()); residual_norms_history.push_back(model_->computeResidualNorms());
// Set up for main Newton loop. // Set up for main Newton loop.
double omega = 1.0; double omega = 1.0;
int iteration = 0; int iteration = 0;
bool converged = model_.getConvergence(dt, iteration); bool converged = model_->getConvergence(dt, iteration);
const int sizeNonLinear = model_.sizeNonLinear(); const int sizeNonLinear = model_->sizeNonLinear();
V dxOld = V::Zero(sizeNonLinear); V dxOld = V::Zero(sizeNonLinear);
bool isOscillate = false; bool isOscillate = false;
bool isStagnate = false; bool isStagnate = false;
@ -82,17 +82,17 @@ namespace Opm
// ---------- Main Newton loop ---------- // ---------- Main Newton loop ----------
while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) { while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
// Compute the Newton update to the primary variables. // Compute the Newton update to the primary variables.
V dx = model_.solveJacobianSystem(); V dx = model_->solveJacobianSystem();
// Store number of linear iterations used. // Store number of linear iterations used.
linearIterations += model_.linearIterationsLastSolve(); linearIterations += model_->linearIterationsLastSolve();
// Stabilize the Newton update. // Stabilize the Newton update.
detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate); detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
if (isOscillate) { if (isOscillate) {
omega -= relaxIncrement(); omega -= relaxIncrement();
omega = std::max(omega, relaxMax()); omega = std::max(omega, relaxMax());
if (model_.terminalOutputEnabled()) { if (model_->terminalOutputEnabled()) {
std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl; std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
} }
} }
@ -100,20 +100,20 @@ namespace Opm
// Apply the update, the model may apply model-dependent // Apply the update, the model may apply model-dependent
// limitations and chopping of the update. // limitations and chopping of the update.
model_.updateState(dx, reservoir_state, well_state); model_->updateState(dx, reservoir_state, well_state);
// Assemble residual and Jacobian, store residual norms. // Assemble residual and Jacobian, store residual norms.
model_.assemble(reservoir_state, well_state, false); model_->assemble(reservoir_state, well_state, false);
residual_norms_history.push_back(model_.computeResidualNorms()); residual_norms_history.push_back(model_->computeResidualNorms());
// increase iteration counter // increase iteration counter
++iteration; ++iteration;
converged = model_.getConvergence(dt, iteration); converged = model_->getConvergence(dt, iteration);
} }
if (!converged) { if (!converged) {
if (model_.terminalOutputEnabled()) { if (model_->terminalOutputEnabled()) {
std::cerr << "WARNING: Failed to compute converged solution in " << iteration << " iterations." << std::endl; std::cerr << "WARNING: Failed to compute converged solution in " << iteration << " iterations." << std::endl;
} }
return -1; // -1 indicates that the solver has to be restarted return -1; // -1 indicates that the solver has to be restarted
@ -125,7 +125,7 @@ namespace Opm
newtonIterationsLast_ = iteration; newtonIterationsLast_ = iteration;
// Do model-specific post-step actions. // Do model-specific post-step actions.
model_.afterStep(dt, reservoir_state, well_state); model_->afterStep(dt, reservoir_state, well_state);
return linearIterations; return linearIterations;
} }
@ -197,7 +197,7 @@ namespace Opm
const std::vector<double>& F0 = residual_history[it]; const std::vector<double>& F0 = residual_history[it];
const std::vector<double>& F1 = residual_history[it - 1]; const std::vector<double>& F1 = residual_history[it - 1];
const std::vector<double>& F2 = residual_history[it - 2]; const std::vector<double>& F2 = residual_history[it - 2];
for (int p= 0; p < model_.numPhases(); ++p){ for (int p= 0; p < model_->numPhases(); ++p){
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]); const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]); const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);

196
opm/autodiff/SimulatorBase.hpp Normal file
View File

@ -0,0 +1,196 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORBASE_HEADER_INCLUDED
#define OPM_SIMULATORBASE_HEADER_INCLUDED
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/NewtonSolver.hpp>
#include <opm/autodiff/BlackoilModel.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm
{
template <class Simulator>
struct SimulatorTraits;
/// Class collecting all necessary components for a two-phase simulation.
template <class Implementation>
class SimulatorBase
{
typedef SimulatorTraits<Implementation> Traits;
public:
typedef typename Traits::ReservoirState ReservoirState;
typedef typename Traits::WellState WellState;
typedef typename Traits::OutputWriter OutputWriter;
typedef typename Traits::Grid Grid;
typedef typename Traits::Solver Solver;
/// Initialise from parameters and objects to observe.
/// \param[in] param parameters, this class accepts the following:
/// parameter (default) effect
/// -----------------------------------------------------------
/// output (true) write output to files?
/// output_dir ("output") output directoty
/// output_interval (1) output every nth step
/// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
/// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal)
/// nl_pressure_maxiter (10) max nonlinear iterations in pressure
/// nl_maxiter (30) max nonlinear iterations in transport
/// nl_tolerance (1e-9) transport solver absolute residual tolerance
/// num_transport_substeps (1) number of transport steps per pressure step
/// use_segregation_split (false) solve for gravity segregation (if false,
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] geo derived geological properties
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] disgas true for dissolved gas option
/// \param[in] vapoil true for vaporized oil option
/// \param[in] eclipse_state
/// \param[in] output_writer
/// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow
SimulatorBase(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool disgas,
const bool vapoil,
std::shared_ptr<EclipseState> eclipse_state,
OutputWriter& output_writer,
const std::vector<double>& threshold_pressures_by_face);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will
/// modify the reservoir and well states.
/// \param[in,out] timer governs the requested reporting timesteps
/// \param[in,out] state state of reservoir: pressure, fluxes
/// \param[in,out] well_state state of wells: bhp, perforation rates
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
ReservoirState& state);
protected:
Implementation& asImpl() { return *static_cast<Implementation*>(this); }
const Implementation& asImpl() const { return *static_cast<const Implementation*>(this); }
void handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& wells_manager,
WellState& well_state,
const Wells* wells);
std::unique_ptr<Solver> createSolver(const Wells* wells);
void
computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellState& xw);
// Data.
typedef RateConverter::
SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
std::vector<int> > RateConverterType;
const parameter::ParameterGroup param_;
// Observed objects.
const Grid& grid_;
BlackoilPropsAdInterface& props_;
const RockCompressibility* rock_comp_props_;
const double* gravity_;
// Solvers
const DerivedGeology& geo_;
NewtonIterationBlackoilInterface& solver_;
// Misc. data
std::vector<int> allcells_;
const bool has_disgas_;
const bool has_vapoil_;
bool terminal_output_;
// eclipse_state
std::shared_ptr<EclipseState> eclipse_state_;
// output_writer
OutputWriter& output_writer_;
RateConverterType rateConverter_;
// Threshold pressures.
std::vector<double> threshold_pressures_by_face_;
// Whether this a parallel simulation or not
bool is_parallel_run_;
};
} // namespace Opm
#include "SimulatorBase_impl.hpp"
#endif // OPM_SIMULATORBASE_HEADER_INCLUDED

258
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp → opm/autodiff/SimulatorBase_impl.hpp
View File

@ -1,6 +1,7 @@
/* /*
Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 IRIS AS Copyright 2014 IRIS AS
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -18,167 +19,24 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/NewtonSolver.hpp>
#include <opm/autodiff/BlackoilModel.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm> #include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm namespace Opm
{ {
template<class T>
class SimulatorFullyImplicitBlackoil<T>::Impl
{
public:
Impl(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
bool has_disgas,
bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
const std::vector<double>& threshold_pressures_by_face);
SimulatorReport run(SimulatorTimer& timer, template <class Implementation>
BlackoilState& state); SimulatorBase<Implementation>::SimulatorBase(const parameter::ParameterGroup& param,
const Grid& grid,
private: const DerivedGeology& geo,
// Data. BlackoilPropsAdInterface& props,
typedef RateConverter:: const RockCompressibility* rock_comp_props,
SurfaceToReservoirVoidage< BlackoilPropsAdInterface, NewtonIterationBlackoilInterface& linsolver,
std::vector<int> > RateConverterType; const double* gravity,
const bool has_disgas,
const parameter::ParameterGroup param_; const bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
// Observed objects. OutputWriter& output_writer,
const Grid& grid_; const std::vector<double>& threshold_pressures_by_face)
BlackoilPropsAdInterface& props_;
const RockCompressibility* rock_comp_props_;
const double* gravity_;
// Solvers
const DerivedGeology& geo_;
NewtonIterationBlackoilInterface& solver_;
// Misc. data
std::vector<int> allcells_;
const bool has_disgas_;
const bool has_vapoil_;
bool terminal_output_;
// eclipse_state
std::shared_ptr<EclipseState> eclipse_state_;
// output_writer
BlackoilOutputWriter& output_writer_;
RateConverterType rateConverter_;
// Threshold pressures.
std::vector<double> threshold_pressures_by_face_;
// Whether this a parallel simulation or not
bool is_parallel_run_;
void
computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw);
};
template<class T>
SimulatorFullyImplicitBlackoil<T>::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
const std::vector<double>& threshold_pressures_by_face)
{
pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil,
eclipse_state, output_writer, threshold_pressures_by_face));
}
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoil<T>::run(SimulatorTimer& timer,
BlackoilState& state)
{
return pimpl_->run(timer, state);
}
// \TODO: Treat bcs.
template<class T>
SimulatorFullyImplicitBlackoil<T>::Impl::Impl(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
const std::vector<double>& threshold_pressures_by_face)
: param_(param), : param_(param),
grid_(grid), grid_(grid),
props_(props), props_(props),
@ -214,14 +72,11 @@ namespace Opm
#endif #endif
} }
template <class Implementation>
SimulatorReport SimulatorBase<Implementation>::run(SimulatorTimer& timer,
ReservoirState& state)
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
BlackoilState& state)
{ {
WellStateFullyImplicitBlackoil prev_well_state; WellState prev_well_state;
// Create timers and file for writing timing info. // Create timers and file for writing timing info.
Opm::time::StopWatch solver_timer; Opm::time::StopWatch solver_timer;
@ -232,14 +87,6 @@ namespace Opm
std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt"; std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str()); std::ofstream tstep_os(tstep_filename.c_str());
typedef T Grid;
typedef BlackoilModel<Grid> Model;
typedef typename Model::ModelParameters ModelParams;
ModelParams modelParams( param_ );
typedef NewtonSolver<Model> Solver;
typedef typename Solver::SolverParameters SolverParams;
SolverParams solverParams( param_ );
// adaptive time stepping // adaptive time stepping
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
if( param_.getDefault("timestep.adaptive", true ) ) if( param_.getDefault("timestep.adaptive", true ) )
@ -282,9 +129,12 @@ namespace Opm
props_.permeability(), props_.permeability(),
is_parallel_run_); is_parallel_run_);
const Wells* wells = wells_manager.c_wells(); const Wells* wells = wells_manager.c_wells();
WellStateFullyImplicitBlackoil well_state; WellState well_state;
well_state.init(wells, state, prev_well_state); well_state.init(wells, state, prev_well_state);
// give the polymer and surfactant simulators the chance to do their stuff
asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
// write simulation state at the report stage // write simulation state at the report stage
output_writer_.writeTimeStep( timer, state, well_state ); output_writer_.writeTimeStep( timer, state, well_state );
@ -293,16 +143,12 @@ namespace Opm
props_.updateSatHyst(state.saturation(), allcells_); props_.updateSatHyst(state.saturation(), allcells_);
// Compute reservoir volumes for RESV controls. // Compute reservoir volumes for RESV controls.
computeRESV(timer.currentStepNum(), wells, state, well_state); asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);
// Run a multiple steps of the solver depending on the time step control. // Run a multiple steps of the solver depending on the time step control.
solver_timer.start(); solver_timer.start();
Model model(modelParams, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_, terminal_output_); auto solver = asImpl().createSolver(wells);
if (!threshold_pressures_by_face_.empty()) {
model.setThresholdPressures(threshold_pressures_by_face_);
}
Solver solver(solverParams, model);
// If sub stepping is enabled allow the solver to sub cycle // If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are to large for the solver to converge // in case the report steps are to large for the solver to converge
@ -310,19 +156,19 @@ namespace Opm
// \Note: The report steps are met in any case // \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables // \Note: The sub stepping will require a copy of the state variables
if( adaptiveTimeStepping ) { if( adaptiveTimeStepping ) {
adaptiveTimeStepping->step( timer, solver, state, well_state, output_writer_ ); adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
} }
else { else {
// solve for complete report step // solve for complete report step
solver.step(timer.currentStepLength(), state, well_state); solver->step(timer.currentStepLength(), state, well_state);
} }
// take time that was used to solve system for this reportStep // take time that was used to solve system for this reportStep
solver_timer.stop(); solver_timer.stop();
// accumulate the number of Newton and Linear Iterations // accumulate the number of Newton and Linear Iterations
totalNewtonIterations += solver.newtonIterations(); totalNewtonIterations += solver->newtonIterations();
totalLinearIterations += solver.linearIterations(); totalLinearIterations += solver->linearIterations();
// Report timing. // Report timing.
const double st = solver_timer.secsSinceStart(); const double st = solver_timer.secsSinceStart();
@ -475,13 +321,47 @@ namespace Opm
} }
} // namespace SimFIBODetails } // namespace SimFIBODetails
template <class T> template <class Implementation>
void void SimulatorBase<Implementation>::handleAdditionalWellInflow(SimulatorTimer& /* timer */,
SimulatorFullyImplicitBlackoil<T>:: WellsManager& /* wells_manager */,
Impl::computeRESV(const std::size_t step, WellState& /* well_state */,
const Wells* wells, const Wells* /* wells */)
const BlackoilState& x, { }
WellStateFullyImplicitBlackoil& xw)
template <class Implementation>
auto SimulatorBase<Implementation>::createSolver(const Wells* wells)
-> std::unique_ptr<Solver>
{
typedef typename Traits::Model Model;
typedef typename Model::ModelParameters ModelParams;
ModelParams modelParams( param_ );
typedef NewtonSolver<Model> Solver;
auto model = std::unique_ptr<Model>(new Model(modelParams,
grid_,
props_,
geo_,
rock_comp_props_,
wells,
solver_,
has_disgas_,
has_vapoil_,
terminal_output_));
if (!threshold_pressures_by_face_.empty()) {
model->setThresholdPressures(threshold_pressures_by_face_);
}
typedef typename Solver::SolverParameters SolverParams;
SolverParams solverParams( param_ );
return std::unique_ptr<Solver>(new Solver(solverParams, std::move(model)));
}
template <class Implementation>
void SimulatorBase<Implementation>::computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellState& xw)
{ {
typedef SimFIBODetails::WellMap WellMap; typedef SimFIBODetails::WellMap WellMap;

121
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -1,5 +1,6 @@
/* /*
Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2015 Andreas Lauser
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -20,91 +21,51 @@
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#include <memory> #include "SimulatorBase.hpp"
#include <vector>
struct UnstructuredGrid; #include "NewtonSolver.hpp"
struct Wells;
struct FlowBoundaryConditions;
namespace Opm namespace Opm {
template <class GridT>
class SimulatorFullyImplicitBlackoil;
template <class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoil<GridT> >
{ {
namespace parameter { class ParameterGroup; } typedef WellStateFullyImplicitBlackoil WellState;
class BlackoilPropsAdInterface; typedef BlackoilState ReservoirState;
class RockCompressibility; typedef BlackoilOutputWriter OutputWriter;
class DerivedGeology; typedef GridT Grid;
class NewtonIterationBlackoilInterface; typedef BlackoilModel<Grid> Model;
class SimulatorTimer; typedef NewtonSolver<Model> Solver;
class BlackoilState; };
class WellStateFullyImplicitBlackoil;
class EclipseState;
class BlackoilOutputWriter;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation. /// a simulator for the blackoil model
template<class T> template <class GridT>
class SimulatorFullyImplicitBlackoil class SimulatorFullyImplicitBlackoil
{ : public SimulatorBase<SimulatorFullyImplicitBlackoil<GridT> >
public: {
/// \brief The type of the grid that we use. typedef SimulatorBase<SimulatorFullyImplicitBlackoil<GridT> > Base;
typedef T Grid; public:
/// Initialise from parameters and objects to observe. // forward the constructor to the base class
/// \param[in] param parameters, this class accepts the following: SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
/// parameter (default) effect const typename Base::Grid& grid,
/// ----------------------------------------------------------- const DerivedGeology& geo,
/// output (true) write output to files? BlackoilPropsAdInterface& props,
/// output_dir ("output") output directoty const RockCompressibility* rock_comp_props,
/// output_interval (1) output every nth step NewtonIterationBlackoilInterface& linsolver,
/// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) const double* gravity,
/// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) const bool disgas,
/// nl_pressure_maxiter (10) max nonlinear iterations in pressure const bool vapoil,
/// nl_maxiter (30) max nonlinear iterations in transport std::shared_ptr<EclipseState> eclipse_state,
/// nl_tolerance (1e-9) transport solver absolute residual tolerance BlackoilOutputWriter& output_writer,
/// num_transport_substeps (1) number of transport steps per pressure step const std::vector<double>& threshold_pressures_by_face)
/// use_segregation_split (false) solve for gravity segregation (if false, : Base(param, grid, geo, props, rock_comp_props, linsolver, gravity, disgas, vapoil,
/// segregation is ignored). eclipse_state, output_writer, threshold_pressures_by_face)
/// {}
/// \param[in] grid grid data structure };
/// \param[in] geo derived geological properties
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] disgas true for dissolved gas option
/// \param[in] vapoil true for vaporized oil option
/// \param[in] eclipse_state
/// \param[in] output_writer
/// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow
SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool disgas,
const bool vapoil,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
const std::vector<double>& threshold_pressures_by_face);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will
/// modify the reservoir and well states.
/// \param[in,out] timer governs the requested reporting timesteps
/// \param[in,out] state state of reservoir: pressure, fluxes
/// \param[in,out] well_state state of wells: bhp, perforation rates
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state);
private:
class Impl;
// Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
std::shared_ptr<Impl> pimpl_;
};
} // namespace Opm } // namespace Opm
#include "SimulatorFullyImplicitBlackoil_impl.hpp"
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED #endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED