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

No functional changes.

Browse Source
This commit is contained in:
Ingeborg Ligaarden 2011-11-25 11:22:09 +01:00
parent d999a8c259
commit 844c5b2b29
2 changed files with 35 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/ImplicitTransport.hpp
View File

@ -46,7 +46,8 @@ namespace Opm {
atol(1.0e-6), atol(1.0e-6),
rtol(5.0e-7), rtol(5.0e-7),
dxtol(1.0e-8), dxtol(1.0e-8),
max_it_ls(5) max_it_ls(5),
verbose(0)
{} {}
int max_it; int max_it;
@ -54,6 +55,7 @@ namespace Opm {
double rtol ; double rtol ;
double dxtol ; double dxtol ;
int max_it_ls; int max_it_ls;
bool verbose;
}; };
struct NRReport { struct NRReport {
@ -88,7 +90,7 @@ namespace Opm {
ReservoirState& state , ReservoirState& state ,
LinearSolver& linsolve, LinearSolver& linsolve,
ImplicitTransportDetails::NRReport& rpt ) { ImplicitTransportDetails::NRReport& rpt ) {
bool s_range; bool s_range; // for use in line search
typedef typename JacobianSystem::vector_type vector_type; typedef typename JacobianSystem::vector_type vector_type;
typedef typename JacobianSystem::matrix_type matrix_type; typedef typename JacobianSystem::matrix_type matrix_type;
typedef typename JacobianSystem::vector_collection_type vector_collection_type; typedef typename JacobianSystem::vector_collection_type vector_collection_type;
@ -108,7 +110,7 @@ namespace Opm {
rpt.norm_dx = -1.0; rpt.norm_dx = -1.0;
rpt.nit = 0; rpt.nit = 0;
bool done = false;//rpt.norm_res < ctrl.atol; bool done = false; //rpt.norm_res < ctrl.atol;
while (! done) { while (! done) {
VZero<vector_type>::zero(sys_.vector().writableIncrement()); VZero<vector_type>::zero(sys_.vector().writableIncrement());
@ -124,14 +126,15 @@ namespace Opm {
rpt.norm_dx = rpt.norm_dx =
VNorm<vector_type>::norm(sys_.vector().increment()); VNorm<vector_type>::norm(sys_.vector().increment());
int lin_it=0; // Begin line search. The line search should be moved to model
double residual=VNorm<vector_type>::norm(sys_.vector().residual()); double residual=VNorm<vector_type>::norm(sys_.vector().residual());
bool finnished=rpt.norm_res<ctrl.atol;//residual < rpt.norm_res; int lin_it=0;
bool finished=rpt.norm_res<ctrl.atol;//residual < rpt.norm_res;
double alpha=2.0; double alpha=2.0;
// store old solution and increasement before line search // store old solution and increment before line search
vector_type dx_old(sys_.vector().increment()); vector_type dx_old(sys_.vector().increment());
vector_type x_old(sys_.vector().solution()); vector_type x_old(sys_.vector().solution());
while(! finnished){ while(! finished){
alpha/=2.0; alpha/=2.0;
s_range = true; s_range = true;
sys_.vector().writableIncrement()=dx_old; sys_.vector().writableIncrement()=dx_old;
@ -149,27 +152,29 @@ namespace Opm {
MZero<matrix_type>::zero(sys_.writableMatrix()); MZero<matrix_type>::zero(sys_.writableMatrix());
VZero<vector_type>::zero(sys_.vector().writableResidual()); VZero<vector_type>::zero(sys_.vector().writableResidual());
asm_.assemble(state, g, src, dt, sys_); asm_.assemble(state, g, src, dt, sys_);
residual=VNorm<vector_type>::norm(sys_.vector().residual()); residual = VNorm<vector_type>::norm(sys_.vector().residual());
if (ctrl.verbose){
std::cout << "Line search iteration " << std::scientific << lin_it << " norm :" << residual << " alpha " << alpha << '\n'; std::cout << "Line search iteration " << std::scientific << lin_it
<< " norm :" << residual << " alpha " << alpha << '\n';
}
}else{ }else{
std::cout << "Line search iteration " << std::scientific << lin_it << " Saturation out of range, continue. Alpha " << alpha << '\n'; if (ctrl.verbose){
std::cout << "Line search iteration " << std::scientific << lin_it
<< " Saturation out of range, continue. Alpha " << alpha << '\n';
}
residual = 1e99; residual = 1e99;
} }
finnished=(residual < rpt.norm_res) || (lin_it> ctrl.max_it_ls); finished=(residual < rpt.norm_res) || (lin_it> ctrl.max_it_ls);
lin_it +=1; lin_it +=1;
}// End line search
}
rpt.norm_res = rpt.norm_res =
VNorm<vector_type>::norm(sys_.vector().residual()); VNorm<vector_type>::norm(sys_.vector().residual());
rpt.nit++; rpt.nit++;
if (ctrl.verbose){
std::cout << "Iteration " << std::scientific << rpt.nit std::cout << "Iteration " << std::scientific << rpt.nit
<< " norm :" << rpt.norm_res << " alpha " << alpha << std::endl; << " norm :" << rpt.norm_res << " alpha " << alpha << std::endl;
}
done = (rpt.norm_res < ctrl.atol) || done = (rpt.norm_res < ctrl.atol) ||
(rpt.norm_res < ctrl.rtol * nrm_res0) || (rpt.norm_res < ctrl.rtol * nrm_res0) ||
(rpt.norm_dx < ctrl.dxtol) || (rpt.norm_dx < ctrl.dxtol) ||

51
src/SinglePointUpwindTwoPhase.hpp
View File

@ -97,8 +97,8 @@ namespace Opm {
double& dpc(int c) { return dpc_[c] ; } double& dpc(int c) { return dpc_[c] ; }
double dpc(int c) const { return dpc_[c] ; } double dpc(int c) const { return dpc_[c] ; }
double& trans(int f) { return trans_[f] ; } double& trans(int f) { return trans_[f] ; }
double trans(int f) const { return trans_[f] ; } double trans(int f) const { return trans_[f] ; }
private: private:
double drho_ ; double drho_ ;
@ -124,20 +124,13 @@ namespace Opm {
const Grid& g , const Grid& g ,
const std::vector<double>& porevol , const std::vector<double>& porevol ,
const double* grav = 0) const double* grav = 0)
// const double* htrans = 0)
: fluid_ (fluid) , : fluid_ (fluid) ,
gravity_(grav) , gravity_(grav) ,
f2hf_ (2 * g.number_of_faces, -1) , f2hf_ (2 * g.number_of_faces, -1) ,
store_ (g.number_of_cells, store_ (g.number_of_cells,
g.cell_facepos[ g.number_of_cells ]) g.cell_facepos[ g.number_of_cells ])
{/* {
int n_hf=g.cell_facepos[ g.number_of_cells ];
if(htrans){
for (int hf = 0; hf < n_hf; ++hf) {
store_.htrans(hf)=htrans[hf];
}
}
*/
if (gravity_) { if (gravity_) {
store_.drho() = fluid_.density(0) - fluid_.density(1); store_.drho() = fluid_.density(0) - fluid_.density(1);
//this->computeStaticGravity(g, gravity_); //this->computeStaticGravity(g, gravity_);
@ -151,7 +144,6 @@ namespace Opm {
} }
} }
std::copy(porevol.begin(), porevol.end(), store_.porevol()); std::copy(porevol.begin(), porevol.end(), store_.porevol());
} }
void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces, void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces,
@ -239,7 +231,7 @@ namespace Opm {
double f1 = m[0] / mt; double f1 = m[0] / mt;
const double v1 = dflux + m[1]*gflux;// gflux+ sgn*trans (pc(n[0]) - pc(n[1])) const double v1 = dflux + m[1]*gflux;
// Assemble residual contributions // Assemble residual contributions
*F += dt * f1 * v1; *F += dt * f1 * v1;
@ -264,16 +256,6 @@ namespace Opm {
/* dF/dm_2 \cdot dm_2/ds */ /* dF/dm_2 \cdot dm_2/ds */
*J[ pix[1] ] -= dt * f1 / mt * v1 * dm[1]; *J[ pix[1] ] -= dt * f1 / mt * v1 * dm[1];
*J[ pix[1] ] += dt * f1 * gflux * dm[1]; *J[ pix[1] ] += dt * f1 * gflux * dm[1];
/* contribution from dpcflux */
//if(sgn>0){
// J1[0*2 + 0] += sgn*dt * f1 * dpcflux[0] * m[1];
// J2[0*2 + 0] -= sgn*dt * f1 * dpcflux[1] * m[1];
//}else{
// J1[0*2 + 0] += sgn*dt * f1 * dpcflux[1] * m[1];
// J2[0*2 + 0] -= sgn*dt * f1 * dpcflux[0] * m[1];
//}
} }
template <class Grid> template <class Grid>
@ -364,7 +346,6 @@ namespace Opm {
JacobianSystem& sys ) { JacobianSystem& sys ) {
// Impose s=0.5 at next time level as an NR initial value. // 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 = typename JacobianSystem::vector_type& x =
sys.vector().writableSolution(); sys.vector().writableSolution();
@ -401,26 +382,14 @@ namespace Opm {
double s_min = fluid_.s_min(c); double s_min = fluid_.s_min(c);
double s_max = fluid_.s_max(c); double s_max = fluid_.s_max(c);
if ( s[0] < (s_min - 1e-5) || s[0] > (s_max + 1e-5) ) { if ( s[0] < (s_min - 1e-5) || s[0] > (s_max + 1e-5) ) {
if (s[0] < s_min){ if (s[0] < s_min){
std::cout << "Warning: s out of range:" << s[0]-s_min << std::endl;} std::cout << "Warning: s out of range:" << s[0]-s_min << std::endl;
if (s[0] > s_max){
std::cout << "Warning: s out of range:" << s[0]-s_max << std::endl;}
//std::cout << "Warning: s0<0." << std::endl;
s_range = false;
if (std::abs(x[c]) > 0){
double m1;
/* if (x[c] < 0){
m1 = -(sat[c*2 + 0]-s_min)/x[c];
} else {
m1 = (s_max-sat[c*2 + 0])/x[c];
}
max_alpha = std::min(max_alpha, m1);
std::cout << "max alpha :" << max_alpha << std::endl;
*/
} }
if (s[0] > s_max){
std::cout << "Warning: s out of range:" << s[0]-s_max << std::endl;
}
s_range = false; //line search fails
} }
s[0] = std::max(s_min, s[0]); s[0] = std::max(s_min, s[0]);