OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-28 16:54:24 -06:00
Code Issues Packages Projects Releases Wiki Activity

revert to upstream/master.

Browse Source
This commit is contained in:
Robert K 2014-10-01 12:54:30 +02:00
parent 23a12421a8
commit 739be5a873
5 changed files with 51 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@ -156,11 +156,57 @@ namespace {
} // Anonymous namespace
template<class T>
void FullyImplicitBlackoilSolver<T>::SolverParameter::
reset()
{
// default values for the solver parameters
dp_max_rel_ = 1.0e9;
ds_max_ = 0.2;
drs_max_rel_ = 1.0e9;
relax_type_ = DAMPEN;
relax_max_ = 0.5;
relax_increment_ = 0.1;
relax_rel_tol_ = 0.2;
max_iter_ = 15;
}
template<class T>
FullyImplicitBlackoilSolver<T>::SolverParameter::
SolverParameter()
{
// set default values
reset();
}
template<class T>
FullyImplicitBlackoilSolver<T>::SolverParameter::
SolverParameter( const parameter::ParameterGroup& param )
{
// set default values
reset();
// overload with given parameters
dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
ds_max_ = param.getDefault("ds_max", ds_max_);
drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_);
relax_max_ = param.getDefault("relax_max", relax_max_);
max_iter_ = param.getDefault("max_iter", max_iter_);
std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
if (relaxation_type == "dampen") {
relax_type_ = DAMPEN;
} else if (relaxation_type == "sor") {
relax_type_ = SOR;
} else {
OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
}
}
template<class T>
FullyImplicitBlackoilSolver<T>::
FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param,
FullyImplicitBlackoilSolver(const SolverParameter& param,
const Grid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
@ -182,36 +228,13 @@ namespace {
, wops_ (wells)
, has_disgas_(has_disgas)
, has_vapoil_(has_vapoil)
, dp_max_rel_ (1.0e9)
, ds_max_ (0.2)
, drs_max_rel_ (1.0e9)
, relax_type_ (DAMPEN)
, relax_max_ (0.5)
, relax_increment_ (0.1)
, relax_rel_tol_ (0.2)
, max_iter_ (15)
, use_threshold_pressure_(false)
, rq_ (fluid.numPhases())
, phaseCondition_(AutoDiffGrid::numCells(grid))
, residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
ADB::null(),
ADB::null() } )
, timeStepControl_()
{
dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
ds_max_ = param.getDefault("ds_max", ds_max_);
drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_);
relax_max_ = param.getDefault("relax_max", relax_max_);
max_iter_ = param.getDefault("max_iter", max_iter_);
std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
if (relaxation_type == "dampen") {
relax_type_ = DAMPEN;
} else if (relaxation_type == "sor") {
relax_type_ = SOR;
} else {
OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
}
}
@ -235,8 +258,10 @@ namespace {
}
template<class T>
double
void
FullyImplicitBlackoilSolver<T>::
step(const double dt,
BlackoilState& x ,
@ -278,14 +303,10 @@ namespace {
bool isOscillate = false;
bool isStagnate = false;
const enum RelaxType relaxtype = relaxType();
int linearIterations = 0 ;
while ((!converged) && (it < maxIter())) {
V dx = solveJacobianSystem();
// store number of linear iterations used
linearIterations += linsolver_.iterationCount();
detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate);
if (isOscillate) {
@ -315,9 +336,6 @@ namespace {
std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
// OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
}
std::cout << "Iterations count " << linearIterations << std::endl;
return timeStepControl_.computeTimeStepSize( dt, linearIterations );
}

6
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@ -108,7 +108,6 @@ namespace Opm
/// Construct a system solver.
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
: iteration_count_( 0 )
{
use_amg_ = param.getDefault("cpr_use_amg", false);
use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
@ -192,7 +191,7 @@ namespace Opm
// Construct linear solver.
const double tolerance = 1e-3;
const int maxit = 5000;
const int verbosity = 0;
const int verbosity = 1;
const int restart = 40;
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
@ -200,9 +199,6 @@ namespace Opm
Dune::InverseOperatorResult result;
linsolve.apply(x, istlb, result);
// store number of iterations used
iteration_count_ = result.iterations;
// Check for failure of linear solver.
if (!result.converged) {
OPM_THROW(std::runtime_error, "Convergence failure for linear solver.");

2
opm/autodiff/NewtonIterationBlackoilCPR.hpp
View File

@ -54,9 +54,7 @@ namespace Opm
/// \return the solution x
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
virtual int iterationCount () const { return iteration_count_; }
private:
mutable int iteration_count_;
bool use_amg_;
bool use_bicgstab_;
};

2
opm/autodiff/NewtonIterationBlackoilInterface.hpp
View File

@ -38,8 +38,6 @@ namespace Opm
/// \param[in] residual residual object containing A and b.
/// \return the solution x
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const = 0;
virtual int iterationCount () const { return 0; }
};
} // namespace Opm

103
opm/autodiff/TimeStepControl.hpp
View File

@ -1,103 +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
{
class TimeStepControlInterface
{
protected:
TimeStepControlInterface() {}
public:
virtual double computeTimeStepSize( const double dt, const int iterations ) const = 0;
virtual ~TimeStepControlInterface () {}
};
class IterationCountTimeStepControl : public TimeStepControlInterface
{
protected:
mutable double prevDt_;
mutable int prevIterations_;
const int targetIterationCount_;
const double adjustmentFactor_;
const int upperTargetIterationCount_;
const int lowerTargetIterationCount_;
public:
IterationCountTimeStepControl()
: prevDt_( 0.0 ), prevIterations_( 0 ),
targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ),
upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 )
{}
double computeTimeStepSize( const double dt, const int iterations ) const
{
// make sure dt is somewhat reliable
assert( dt > 0 && dt == dt );
double newDt = dt;
double derivation = double(std::abs( iterations - targetIterationCount_ )) / double(targetIterationCount_);
if( derivation > 0.1 )
{
if( iterations < targetIterationCount_ )
{
newDt = dt * adjustmentFactor_;
}
else
newDt = dt / adjustmentFactor_;
}
/*
if( prevDt_ > 0 && std::abs( dt - prevDt_ ) > 1e-12 ) {
const double dFdt = double(iterations - prevIterations_) / ( dt - prevDt_ );
if( std::abs( dFdt ) > 1e-12 )
newDt = dt + (targetIterationCount_ - iterations) / dFdt;
else
// if iterations was the same or dts were the same, do some magic
newDt = dt * double( targetIterationCount_ ) / double(targetIterationCount_ - iterations);
}
if( newDt < 0 || ! (prevDt_ > 0) || ( iterations == prevIterations_) )
{
if( iterations > upperTargetIterationCount_ )
newDt = dt / adjustmentFactor_;
else if( iterations < lowerTargetIterationCount_ )
newDt = dt * adjustmentFactor_;
else
newDt = dt;
}
*/
assert( newDt == newDt );
//std::cout << "dt = " << dt << " " << prevDt_ << std::endl;
prevDt_ = dt;
prevIterations_ = iterations;
return newDt;
}
};
} // end namespace OPM
#endif