OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

checking after introducing linesearch, do maybe work, but did not help our objektives

Browse Source
This commit is contained in:
Halvor M. Nilsen 2011-10-18 15:08:30 +02:00
parent 6ce89ee3c3
commit 4750e17699
3 changed files with 45 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
src/ImplicitTransport.hpp
View File

@ -37,7 +37,7 @@
#define OPM_IMPLICITTRANSPORT_HPP_HEADER
#include "ImplicitAssembly.hpp"
#include <boost/lambda/lambda.hpp>
namespace Opm {
namespace ImplicitTransportDetails {
struct NRControl {
@ -45,13 +45,15 @@ namespace Opm {
: max_it(1),
atol(1.0e-6),
rtol(5.0e-7),
dxtol(1.0e-8)
dxtol(1.0e-8),
max_it_ls(5)
{}
int max_it;
double atol ;
double rtol ;
double dxtol ;
int max_it_ls;
};
struct NRReport {
@ -89,7 +91,7 @@ namespace Opm {
typedef typename JacobianSystem::vector_type vector_type;
typedef typename JacobianSystem::matrix_type matrix_type;
typedef typename JacobianSystem::vector_collection_type vector_collection_type;
asm_.createSystem(g, sys_);
model_.initStep(state, g, sys_);
@ -107,7 +109,7 @@ namespace Opm {
rpt.norm_dx = -1.0;
rpt.nit = 0;
bool done = rpt.norm_res < ctrl.atol;
bool done = false;//rpt.norm_res < ctrl.atol;
while (! done) {
linsolve.solve(sys_.matrix(),
@ -121,25 +123,46 @@ namespace Opm {
rpt.norm_dx =
VNorm<vector_type>::norm(sys_.vector().increment());
//std::cerr << rpt.norm_dx << '\n';
sys_.vector().addIncrement();
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_);
int lin_it=0;
double residual=VNorm<vector_type>::norm(sys_.vector().residual());
bool finnished=false;//residual < rpt.norm_res;
double alpha=2.0;
// store old solution and increasement before line search
vector_type dx_old(sys_.vector().increment());
vector_type x_old(sys_.vector().solution());
while(! finnished){
alpha/=2.0;
sys_.vector().writableIncrement()=dx_old;
sys_.vector().writableIncrement()*=alpha;
sys_.vector().writableSolution()=x_old;
/*
* should be used if vector_type is std::vector<double>
std::vector<double> operator*=(std::vector<double>& vec,const double a){
for_each(vec.begin(),vec.end(), boost::lambda::_1*=a );
return vec;
}
*/
sys_.vector().addIncrement();
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_);
residual=VNorm<vector_type>::norm(sys_.vector().residual());
lin_it +=1;
finnished=(residual < rpt.norm_res) || (lin_it> ctrl.max_it_ls);
std::cerr << "Line search iteration " << std::scientific << lin_it << " norm :" << residual << " alpha " << alpha << '\n';
}
rpt.norm_res =
VNorm<vector_type>::norm(sys_.vector().residual());
rpt.nit++;
std::cerr << "Iteration " << std::scientific << rpt.nit << " norm :" << rpt.norm_res << '\n';
std::cerr << "Iteration " << std::scientific << rpt.nit << " norm :" << rpt.norm_res << " alpha " << alpha << '\n';
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);
}

7
src/JacobianSystem.hpp
View File

@ -54,7 +54,6 @@ namespace Opm {
static void
add(const BaseVec& x, BaseVec& y) {
typedef typename BaseVec::value_type VT;
::std::transform(x.begin(), x.end(),
y.begin(),
y.begin(),
@ -69,10 +68,9 @@ namespace Opm {
static void
negate(BaseVec& x) {
typedef typename BaseVec::value_type VT;
::std::transform(x.begin(), x.end(),
x.begin(),
::std::negate<VT>());
x.begin(),
::std::negate<VT>());
}
};
@ -208,6 +206,7 @@ namespace Opm {
typedef Matrix matrix_type;
typedef MatrixBlockAssembler<Matrix> assembler_type;
typedef NVecCollection vector_collection_type;
typedef typename NVecCollection::vector_type vector_type;
assembler_type& matasm() { return mba_ ; }

7
src/SinglePointUpwindTwoPhase.hpp
View File

@ -291,7 +291,7 @@ namespace Opm {
// Impose s=0.5 at next time level as an NR initial value.
const ::std::vector<double>& s = state.saturation();
//const ::std::vector<double>& s = state.saturation();
typename JacobianSystem::vector_type& x =
sys.vector().writableSolution();
@ -343,6 +343,7 @@ namespace Opm {
NewtonIterate& it ) {
// Nothing to do at end of iteration in this model.
(void) state; (void) g; (void) it;
typedef typename NewtonIterate::vector_type vector_t;
}
template <class Grid ,
@ -356,8 +357,10 @@ namespace Opm {
double *s = &state.saturation()[0*2 + 0];
for (int c = 0; c < g.number_of_cells; ++c, s += 2) {
s[0] += x[c] ;
s[0] += x[c] ;
s[1] = 1 - s[0];
assert(s[0]<=1);
assert(s[0]<=1);
}
}