diff --git a/src/ImplicitTransport.hpp b/src/ImplicitTransport.hpp index 08d60204..110cb0d4 100644 --- a/src/ImplicitTransport.hpp +++ b/src/ImplicitTransport.hpp @@ -46,7 +46,8 @@ namespace Opm { atol(1.0e-6), rtol(5.0e-7), dxtol(1.0e-8), - max_it_ls(5) + max_it_ls(5), + verbose(0) {} int max_it; @@ -54,6 +55,7 @@ namespace Opm { double rtol ; double dxtol ; int max_it_ls; + bool verbose; }; struct NRReport { @@ -88,7 +90,7 @@ namespace Opm { ReservoirState& state , LinearSolver& linsolve, ImplicitTransportDetails::NRReport& rpt ) { - bool s_range; + bool s_range; // for use in line search typedef typename JacobianSystem::vector_type vector_type; typedef typename JacobianSystem::matrix_type matrix_type; typedef typename JacobianSystem::vector_collection_type vector_collection_type; @@ -108,7 +110,7 @@ namespace Opm { rpt.norm_dx = -1.0; rpt.nit = 0; - bool done = false;//rpt.norm_res < ctrl.atol; + bool done = false; //rpt.norm_res < ctrl.atol; while (! done) { VZero::zero(sys_.vector().writableIncrement()); @@ -124,14 +126,15 @@ namespace Opm { rpt.norm_dx = VNorm::norm(sys_.vector().increment()); - int lin_it=0; + // Begin line search. The line search should be moved to model double residual=VNorm::norm(sys_.vector().residual()); - bool finnished=rpt.norm_res::zero(sys_.writableMatrix()); VZero::zero(sys_.vector().writableResidual()); asm_.assemble(state, g, src, dt, sys_); - residual=VNorm::norm(sys_.vector().residual()); - - std::cout << "Line search iteration " << std::scientific << lin_it << " norm :" << residual << " alpha " << alpha << '\n'; - - + residual = VNorm::norm(sys_.vector().residual()); + if (ctrl.verbose){ + std::cout << "Line search iteration " << std::scientific << lin_it + << " norm :" << residual << " alpha " << alpha << '\n'; + } }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; } - 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; - - } + }// End line search rpt.norm_res = VNorm::norm(sys_.vector().residual()); rpt.nit++; - - std::cout << "Iteration " << std::scientific << rpt.nit - << " norm :" << rpt.norm_res << " alpha " << alpha << std::endl; - + if (ctrl.verbose){ + 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) || diff --git a/src/SinglePointUpwindTwoPhase.hpp b/src/SinglePointUpwindTwoPhase.hpp index 11b13525..ce41355a 100644 --- a/src/SinglePointUpwindTwoPhase.hpp +++ b/src/SinglePointUpwindTwoPhase.hpp @@ -97,8 +97,8 @@ namespace Opm { double& dpc(int c) { return dpc_[c] ; } double dpc(int c) const { return dpc_[c] ; } - double& trans(int f) { return trans_[f] ; } - double trans(int f) const { return trans_[f] ; } + double& trans(int f) { return trans_[f] ; } + double trans(int f) const { return trans_[f] ; } private: double drho_ ; @@ -124,20 +124,13 @@ namespace Opm { const Grid& g , const std::vector& porevol , const double* grav = 0) - // const double* htrans = 0) : fluid_ (fluid) , gravity_(grav) , f2hf_ (2 * g.number_of_faces, -1) , store_ (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_) { store_.drho() = fluid_.density(0) - fluid_.density(1); //this->computeStaticGravity(g, gravity_); @@ -151,7 +144,6 @@ namespace Opm { } } - std::copy(porevol.begin(), porevol.end(), store_.porevol()); } void makefhfQPeriodic( const std::vector& p_faces,const std::vector& hf_faces, @@ -239,7 +231,7 @@ namespace Opm { 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 *F += dt * f1 * v1; @@ -264,16 +256,6 @@ namespace Opm { /* dF/dm_2 \cdot dm_2/ds */ *J[ pix[1] ] -= dt * f1 / mt * v1 * 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 @@ -364,7 +346,6 @@ namespace Opm { JacobianSystem& sys ) { // Impose s=0.5 at next time level as an NR initial value. - //const ::std::vector& s = state.saturation(); typename JacobianSystem::vector_type& x = sys.vector().writableSolution(); @@ -401,26 +382,14 @@ namespace Opm { double s_min = fluid_.s_min(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){ - 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; - */ + 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; + } + s_range = false; //line search fails } s[0] = std::max(s_min, s[0]);