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

moved the utility classes to opm-core.

Browse Source
This commit is contained in:
Robert K 2014-10-06 13:59:21 +02:00
parent faf191b9f1
commit d4802121d3
2 changed files with 14 additions and 264 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

87
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -28,7 +28,6 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/TimeStepControl.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
@ -47,6 +46,7 @@
#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>
@ -54,6 +54,7 @@
#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>
@ -294,12 +295,14 @@ namespace Opm
std::string tstep_filename = output_dir_ + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
double lastSubStep = timer.currentStepLength();
typename FullyImplicitBlackoilSolver<T>::SolverParameter solverParam( param_ );
// sub stepping flag, if false just the normal time steps will be used
const bool subStepping = param_.getDefault("substepping", bool(false) );
// sub stepping
std::unique_ptr< AdaptiveTimeStepping > subStepping;
if( param_.getDefault("timestep.adaptive", bool(false) ) )
{
subStepping = std::unique_ptr< AdaptiveTimeStepping > (new AdaptiveTimeStepping( param_ ));
}
// create time step control object, TODO introduce parameter
std::unique_ptr< TimeStepControlInterface >
@ -362,74 +365,12 @@ namespace Opm
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are to large for the solver to converge
// \Note: The report steps are met in any case
if( subStepping )
{
// create sub step simulator timer with previously used sub step size
const double start_time = timer.simulationTimeElapsed();
const double end_time = start_time + timer.currentStepLength();
AdaptiveSimulatorTimer subStepper( start_time, end_time, lastSubStep );
// copy states in case solver has to be restarted (to be revised)
BlackoilState last_state( state );
WellStateFullyImplicitBlackoil last_well_state( well_state );
// sub step time loop
while( ! subStepper.done() )
{
// initialize time step control in case current state is needed later
timeStepControl->initialize( state );
int linearIterations = -1;
try {
// (linearIterations < 0 means on convergence in solver)
linearIterations = solver.step(subStepper.currentStepLength(), state, well_state);
// report number of linear iterations
std::cout << "Overall linear iterations used: " << linearIterations << std::endl;
}
catch (Opm::NumericalProblem)
{
// since linearIterations is < 0 this will restart the solver
}
// (linearIterations < 0 means on convergence in solver)
if( linearIterations >= 0 )
{
// advance by current dt
subStepper.advance();
// compute new time step estimate
const double dtEstimate =
timeStepControl->computeTimeStepSize( subStepper.currentStepLength(), linearIterations, state );
std::cout << "Suggested time step size = " << dtEstimate/86400.0 << " (days)" << std::endl;
// set new time step length
subStepper.provideTimeStepEstimate( dtEstimate );
// update states
last_state = state ;
last_well_state = well_state;
}
else // in case of no convergence
{
// we need to revise this
subStepper.provideTimeStepEstimate( 0.1 * subStepper.currentStepLength() );
std::cerr << "Solver convergence failed, restarting solver with half time step ("<< subStepper.currentStepLength()<<" days)." << std::endl;
// reset states
state = last_state;
well_state = last_well_state;
}
}
subStepper.report( std::cout );
// store last small time step for next reportStep
lastSubStep = subStepper.suggestedAverage();
std::cout << "Last suggested step size = " << lastSubStep/86400.0 << " (days)" << std::endl;
if( ! std::isfinite( lastSubStep ) ) // check for NaN
lastSubStep = timer.currentStepLength();
//
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if( subStepping ) {
subStepping->step( solver, state, well_state,
timer.simulationTimeElapsed(), timer.currentStepLength() );
}
else {
// solve for complete report step

191
opm/autodiff/TimeStepControl.hpp
View File

@ -1,191 +0,0 @@
/*
Copyright 2014 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TIMESTEPCONTROL_HEADER_INCLUDED
#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED
namespace Opm
{
///////////////////////////////////////////////////////////////////
///
/// TimeStepControlInterface
///
///////////////////////////////////////////////////////////////////
class TimeStepControlInterface
{
protected:
TimeStepControlInterface() {}
public:
/// \param state simulation state before computing update in the solver (default is empty)
virtual void initialize( const SimulatorState& state ) {}
/// compute new time step size suggestions based on the PID controller
/// \param dt time step size used in the current step
/// \param iterations number of iterations used (linear/nonlinear)
/// \param state new solution state
///
/// \return suggested time step size for the next step
virtual double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const = 0;
/// virtual destructor (empty)
virtual ~TimeStepControlInterface () {}
};
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
///
/// PID controller based adaptive time step control as suggested in:
/// Turek and Kuzmin. Algebraic Flux Correction III. Incompressible Flow Problems. Uni Dortmund.
///
/// See also:
/// D. Kuzmin and S.Turek. Numerical simulation of turbulent bubbly flows. Techreport Uni Dortmund. 2004
///
/// and the original article:
/// Valli, Coutinho, and Carey. Adaptive Control for Time Step Selection in Finite Element
/// Simulation of Coupled Viscous Flow and Heat Transfer. Proc of the 10th
/// International Conference on Numerical Methods in Fluids. 1998.
///
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
class PIDTimeStepControl : public TimeStepControlInterface
{
protected:
mutable std::vector<double> p0_;
mutable std::vector<double> sat0_;
const double tol_;
mutable std::vector< double > errors_;
const bool verbose_;
public:
/// \brief constructor
/// \param tol tolerance for the relative changes of the numerical solution to be accepted
/// in one time step (default is 1e-3)
PIDTimeStepControl( const double tol = 1e-3, const bool verbose = false )
: p0_()
, sat0_()
, tol_( tol )
, errors_( 3, tol_ )
, verbose_( verbose )
{}
/// \brief \copydoc TimeStepControlInterface::initialize
void initialize( const SimulatorState& state )
{
// store current state for later time step computation
p0_ = state.pressure();
sat0_ = state.saturation();
}
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
double computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const
{
const size_t size = p0_.size();
assert( state.pressure().size() == size );
assert( state.saturation().size() == size );
assert( sat0_.size() == size );
// compute u^n - u^n+1
for( size_t i=0; i<size; ++i )
{
p0_[ i ] -= state.pressure()[ i ];
sat0_[ i ] -= state.saturation()[ i ];
}
// compute || u^n - u^n+1 ||
const double stateOld = inner_product( p0_.begin(), p0_.end() ) +
inner_product( sat0_.begin(), sat0_.end() );
// compute || u^n+1 ||
const double stateNew = inner_product( state.pressure().begin(), state.pressure().end() ) +
inner_product( state.saturation().begin(), state.saturation().end() );
// shift errors
for( int i=0; i<2; ++i )
errors_[ i ] = errors_[i+1];
// store new error
const double error = stateOld / stateNew;
errors_[ 2 ] = error ;
if( error > tol_ )
{
// adjust dt by given tolerance
if( verbose_ )
std::cout << "Computed step size (tol): " << (dt * tol_ / error )/86400.0 << " (days)" << std::endl;
return (dt * tol_ / error );
}
else
{
// values taking from turek time stepping paper
const double kP = 0.075 ;
const double kI = 0.175 ;
const double kD = 0.01 ;
double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
std::pow( tol_ / errors_[ 2 ], kI ) *
std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD ));
if( verbose_ )
std::cout << "Computed step size (pow): " << newDt/86400.0 << " (days)" << std::endl;
return newDt;
}
}
protected:
// return inner product for given container, here std::vector
template <class Iterator>
double inner_product( Iterator it, const Iterator end ) const
{
double product = 0.0 ;
for( ; it != end; ++it )
product += ( *it * *it );
return product;
}
};
class PIDAndIterationCountTimeStepControl : public PIDTimeStepControl
{
typedef PIDTimeStepControl BaseType;
protected:
const int targetIterationCount_;
public:
explicit PIDAndIterationCountTimeStepControl( const int target_iterations = 20,
const double tol = 1e-3,
const bool verbose = false)
: BaseType( tol, verbose )
, targetIterationCount_( target_iterations )
{}
double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const
{
double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, state );
// further reduce step size if to many iterations were used
if( iterations > targetIterationCount_ )
{
// if iterations was the same or dts were the same, do some magic
dtEstimate *= double( targetIterationCount_ ) / double(iterations);
}
return dtEstimate;
}
};
} // end namespace OPM
#endif