Loading...
Searching...
No Matches
ImplicitTransport.hpp
/*===========================================================================
//
// File: ImplicitTransport.hpp
//
// Created: 2011-09-29 10:38:42+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
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_IMPLICITTRANSPORT_HPP_HEADER
#define OPM_IMPLICITTRANSPORT_HPP_HEADER
#include <opm/core/transport/implicit/ImplicitAssembly.hpp>
#include <iostream>
namespace Opm {
namespace ImplicitTransportDetails {
struct NRControl {
NRControl()
: max_it(1),
atol(1.0e-6),
rtol(5.0e-7),
dxtol(1.0e-8),
max_it_ls(5),
verbosity(0)
{}
int max_it;
double atol ;
double rtol ;
double dxtol ;
int max_it_ls;
int verbosity;
};
struct NRReport {
int nit;
int flag;
double norm_res;
double norm_dx;
};
template <class Ostream>
Ostream& operator<<(Ostream& os, const NRReport& rpt)
{
os << "Number of linear solves: " << rpt.nit << '\n'
<< "Process converged: " << (rpt.flag > 0) << '\n'
<< "Convergence flag: " << rpt.flag;
switch (rpt.flag) {
case 1:
os << " (absolute tolerance satisfied)\n";
break;
case 2:
os << " (relative tolerance satisfied)\n";
break;
case 3:
os << " (solution change tolerance satisfied)\n";
break;
case -1:
os << " (failed to converge)\n";
break;
default:
os << " (unknown convergence flag)\n";
}
os << "Final residual norm: " << rpt.norm_res << '\n'
<< "Final increment norm: " << rpt.norm_dx << '\n';
return os;
}
}
template <class Model ,
class JacobianSystem ,
template <class> class VNorm ,
template <class> class VNeg ,
template <class> class VZero ,
template <class> class MZero ,
template <class> class VAsgn >
class ImplicitTransport {
public:
explicit ImplicitTransport(Model& model)
: model_(model),
asm_ (model)
{}
template <class Grid ,
class SourceTerms ,
class ReservoirState,
class LinearSolver >
void solve(const Grid& g ,
const SourceTerms* src ,
const double dt ,
const ImplicitTransportDetails::NRControl& ctrl ,
ReservoirState& state ,
LinearSolver& linsolve,
ImplicitTransportDetails::NRReport& rpt )
{
typedef typename JacobianSystem::vector_type vector_type;
typedef typename JacobianSystem::matrix_type matrix_type;
asm_.createSystem(g, sys_);
model_.initStep(state, g, sys_);
model_.initIteration(state, g, sys_);
MZero<matrix_type>::zero(sys_.writableMatrix());
VZero<vector_type>::zero(sys_.vector().writableResidual());
asm_.assemble(state, g, src, dt, sys_);
const double nrm_res0 =
VNorm<vector_type>::norm(sys_.vector().residual());
rpt.norm_res = nrm_res0;
rpt.norm_dx = -1.0;
rpt.nit = 0;
bool done = false; //rpt.norm_res < ctrl.atol;
while (! done) {
VZero<vector_type>::zero(sys_.vector().writableIncrement());
linsolve.solve(sys_.matrix(),
sys_.vector().residual(),
sys_.vector().writableIncrement());
VNeg<vector_type>::negate(sys_.vector().writableIncrement());
//model_.finishIteration(state, g, sys_.vector());
rpt.norm_dx =
VNorm<vector_type>::norm(sys_.vector().increment());
// Begin line search
double residual;
int lin_it=0;
bool finished=rpt.norm_res<ctrl.atol;
double alpha=2.0;
// store old solution and increment before line search
vector_type dx_old(sys_.vector().increment());
vector_type x_old(sys_.vector().solution());
while(! finished){
alpha/=2.0;
VAsgn<vector_type>::assign(alpha, dx_old,
sys_.vector().writableIncrement());
VAsgn<vector_type>::assign(x_old,
sys_.vector().writableSolution());
sys_.vector().addIncrement();
bool init = model_.initIteration(state, g, sys_);
if (init) {
MZero<matrix_type>::zero(sys_.writableMatrix());
VZero<vector_type>::zero(sys_.vector().writableResidual());
asm_.assemble(state, g, src, dt, sys_);
residual = VNorm<vector_type>::norm(sys_.vector().residual());
if (ctrl.verbosity > 1){
std::cout << "Line search iteration " << std::scientific << lin_it
<< " norm :" << residual << " alpha " << alpha << '\n';
}
}else{
if (ctrl.verbosity > 1){
std::cout << "Line search iteration " << std::scientific << lin_it
<< " Value out of range, continue search. alpha " << alpha << '\n';
}
residual = 1e99;
}
finished=(residual < rpt.norm_res) || (lin_it> ctrl.max_it_ls);
lin_it +=1;
}// End line search
rpt.norm_res =
VNorm<vector_type>::norm(sys_.vector().residual());
rpt.nit++;
if (ctrl.verbosity > 0){
std::cout << "Iteration " << std::scientific << rpt.nit
<< " norm :" << rpt.norm_res << " alpha " << alpha << std::endl;
}
done = (rpt.norm_res < ctrl.atol) ||
(rpt.norm_res < ctrl.rtol * nrm_res0) ||
(rpt.norm_dx < ctrl.dxtol) ||
(lin_it > ctrl.max_it_ls) ||
(rpt.nit == ctrl.max_it);
}
model_.finishStep(g, sys_.vector().solution(), state);
if (rpt.norm_res < ctrl.atol) { rpt.flag = 1; }
else if (rpt.norm_res < ctrl.rtol * nrm_res0) { rpt.flag = 2; }
else if (rpt.norm_dx < ctrl.dxtol) { rpt.flag = 3; }
else { rpt.flag = -1; }
}
private:
ImplicitTransport (const ImplicitTransport&);
ImplicitTransport& operator=(const ImplicitTransport&);
#if 0
using Model::initStep;
using Model::initIteration;
using Model::finishIteration;
using Model::finishStep;
#endif
Model& model_;
ImplicitAssembly<Model> asm_;
JacobianSystem sys_;
};
}
#endif /* OPM_IMPLICITTRANSPORT_HPP_HEADER */
Opm::ImplicitAssembly
Definition ImplicitAssembly.hpp:45
Opm::ReservoirState
Definition ImplicitTransportDefs.hpp:52
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::ImplicitTransportDetails::NRControl
Definition ImplicitTransport.hpp:45
Opm::ImplicitTransportDetails::NRReport
Definition ImplicitTransport.hpp:63