00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
00024 #define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
00025
00026 #include <opm/autodiff/NonlinearSolver.hpp>
00027 #include <opm/common/Exceptions.hpp>
00028 #include <opm/common/ErrorMacros.hpp>
00029
00030 namespace Opm
00031 {
00032 template <class PhysicalModel>
00033 NonlinearSolver<PhysicalModel>::NonlinearSolver(const SolverParameters& param,
00034 std::unique_ptr<PhysicalModel> model_arg)
00035 : param_(param),
00036 model_(std::move(model_arg)),
00037 linearizations_(0),
00038 nonlinearIterations_(0),
00039 linearIterations_(0),
00040 wellIterations_(0),
00041 nonlinearIterationsLast_(0),
00042 linearIterationsLast_(0),
00043 wellIterationsLast_(0)
00044 {
00045 if (!model_) {
00046 OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
00047 }
00048 }
00049
00050 template <class PhysicalModel>
00051 int NonlinearSolver<PhysicalModel>::linearizations() const
00052 {
00053 return linearizations_;
00054 }
00055
00056 template <class PhysicalModel>
00057 int NonlinearSolver<PhysicalModel>::nonlinearIterations() const
00058 {
00059 return nonlinearIterations_;
00060 }
00061
00062 template <class PhysicalModel>
00063 int NonlinearSolver<PhysicalModel>::linearIterations() const
00064 {
00065 return linearIterations_;
00066 }
00067
00068 template <class PhysicalModel>
00069 int NonlinearSolver<PhysicalModel>::wellIterations() const
00070 {
00071 return wellIterations_;
00072 }
00073
00074 template <class PhysicalModel>
00075 const PhysicalModel& NonlinearSolver<PhysicalModel>::model() const
00076 {
00077 return *model_;
00078 }
00079
00080 template <class PhysicalModel>
00081 PhysicalModel& NonlinearSolver<PhysicalModel>::model()
00082 {
00083 return *model_;
00084 }
00085
00086 template <class PhysicalModel>
00087 int NonlinearSolver<PhysicalModel>::nonlinearIterationsLastStep() const
00088 {
00089 return nonlinearIterationsLast_;
00090 }
00091
00092 template <class PhysicalModel>
00093 int NonlinearSolver<PhysicalModel>::linearIterationsLastStep() const
00094 {
00095 return linearIterationsLast_;
00096 }
00097
00098 template <class PhysicalModel>
00099 int NonlinearSolver<PhysicalModel>::wellIterationsLastStep() const
00100 {
00101 return wellIterationsLast_;
00102 }
00103
00104 template <class PhysicalModel>
00105 SimulatorReport
00106 NonlinearSolver<PhysicalModel>::
00107 step(const SimulatorTimerInterface& timer,
00108 ReservoirState& reservoir_state,
00109 WellState& well_state)
00110 {
00111 return step(timer, reservoir_state, well_state, reservoir_state, well_state);
00112 }
00113
00114
00115
00116 template <class PhysicalModel>
00117 SimulatorReport
00118 NonlinearSolver<PhysicalModel>::
00119 step(const SimulatorTimerInterface& timer,
00120 const ReservoirState& initial_reservoir_state,
00121 const WellState& initial_well_state,
00122 ReservoirState& reservoir_state,
00123 WellState& well_state)
00124 {
00125 SimulatorReport iterReport;
00126 SimulatorReport report;
00127 failureReport_ = SimulatorReport();
00128
00129
00130 model_->prepareStep(timer, initial_reservoir_state, initial_well_state);
00131
00132 int iteration = 0;
00133
00134
00135
00136
00137 bool converged = false;
00138
00139
00140 do {
00141 try {
00142
00143
00144
00145 iterReport = model_->nonlinearIteration(iteration, timer, *this, reservoir_state, well_state);
00146
00147 report += iterReport;
00148 report.converged = iterReport.converged;
00149
00150 converged = report.converged;
00151 iteration += 1;
00152 }
00153 catch (...) {
00154
00155
00156 failureReport_ += report;
00157 failureReport_ += model_->failureReport();
00158 throw;
00159 }
00160 } while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
00161
00162 if (!converged) {
00163 failureReport_ += report;
00164
00165 std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
00166 OPM_THROW_NOLOG(Opm::TooManyIterations, msg);
00167 }
00168
00169
00170 model_->afterStep(timer, reservoir_state, well_state);
00171 report.converged = true;
00172
00173 return report;
00174 }
00175
00176
00177
00178 template <class PhysicalModel>
00179 void NonlinearSolver<PhysicalModel>::SolverParameters::
00180 reset()
00181 {
00182
00183 relax_type_ = DAMPEN;
00184 relax_max_ = 0.5;
00185 relax_increment_ = 0.1;
00186 relax_rel_tol_ = 0.2;
00187 max_iter_ = 10;
00188 min_iter_ = 1;
00189 }
00190
00191 template <class PhysicalModel>
00192 NonlinearSolver<PhysicalModel>::SolverParameters::
00193 SolverParameters()
00194 {
00195
00196 reset();
00197 }
00198
00199 template <class PhysicalModel>
00200 NonlinearSolver<PhysicalModel>::SolverParameters::
00201 SolverParameters( const ParameterGroup& param )
00202 {
00203
00204 reset();
00205
00206
00207 relax_max_ = param.getDefault("relax_max", relax_max_);
00208 max_iter_ = param.getDefault("max_iter", max_iter_);
00209 min_iter_ = param.getDefault("min_iter", min_iter_);
00210
00211 std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
00212 if (relaxation_type == "dampen") {
00213 relax_type_ = DAMPEN;
00214 } else if (relaxation_type == "sor") {
00215 relax_type_ = SOR;
00216 } else {
00217 OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
00218 }
00219 }
00220
00221 template <class PhysicalModel>
00222 void
00223 NonlinearSolver<PhysicalModel>::detectOscillations(const std::vector<std::vector<double>>& residual_history,
00224 const int it,
00225 bool& oscillate, bool& stagnate) const
00226 {
00227
00228
00229
00230
00231
00232 if ( it < 2 ) {
00233 oscillate = false;
00234 stagnate = false;
00235 return;
00236 }
00237
00238 stagnate = true;
00239 int oscillatePhase = 0;
00240 const std::vector<double>& F0 = residual_history[it];
00241 const std::vector<double>& F1 = residual_history[it - 1];
00242 const std::vector<double>& F2 = residual_history[it - 2];
00243 for (int p= 0; p < model_->numPhases(); ++p){
00244 const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
00245 const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
00246
00247 oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
00248
00249
00250
00251 stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
00252 }
00253
00254 oscillate = (oscillatePhase > 1);
00255 }
00256
00257
00258 template <class PhysicalModel>
00259 template <class BVector>
00260 void
00261 NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
00262 {
00263
00264
00265
00266 BVector tempDxOld = dxOld;
00267 dxOld = dx;
00268
00269 switch (relaxType()) {
00270 case DAMPEN: {
00271 if (omega == 1.) {
00272 return;
00273 }
00274 auto i = dx.size();
00275 for (i = 0; i < dx.size(); ++i) {
00276 dx[i] *= omega;
00277 }
00278 return;
00279 }
00280 case SOR: {
00281 if (omega == 1.) {
00282 return;
00283 }
00284 auto i = dx.size();
00285 for (i = 0; i < dx.size(); ++i) {
00286 dx[i] *= omega;
00287 tempDxOld[i] *= (1.-omega);
00288 dx[i] += tempDxOld[i];
00289 }
00290 return;
00291 }
00292 default:
00293 OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
00294 }
00295
00296 return;
00297 }
00298 }
00299
00300
00301 #endif // OPM_FULLYIMPLICITSOLVER_IMPL_HEADER_INCLUDED