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

Merged changes.

Browse Source
This commit is contained in:
Ingeborg Ligaarden 2011-11-25 14:43:03 +01:00
commit 9dbd4d7773
4 changed files with 21 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/cart_grid.c
View File

@ -100,6 +100,7 @@ create_cart_grid(int nx, int ny, int nz)
cfacepos = G->cell_facepos;
ccentroids = G->cell_centroids;
cvolumes = G->cell_volumes;
cfacepos[0] = 0;
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {

25
src/ImplicitTransport.hpp
View File

@ -47,7 +47,7 @@ namespace Opm {
rtol(5.0e-7),
dxtol(1.0e-8),
max_it_ls(5),
verbose(0)
verbosity(0)
{}
int max_it;
@ -55,7 +55,7 @@ namespace Opm {
double rtol ;
double dxtol ;
int max_it_ls;
bool verbose;
int verbosity;
};
struct NRReport {
@ -90,13 +90,13 @@ namespace Opm {
ReservoirState& state ,
LinearSolver& linsolve,
ImplicitTransportDetails::NRReport& rpt ) {
bool s_range; // for use in line search
bool init;
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_);
model_.initIteration(state, g, sys_, s_range);
init = model_.initIteration(state, g, sys_);
MZero<matrix_type>::zero(sys_.writableMatrix());
VZero<vector_type>::zero(sys_.vector().writableResidual());
@ -126,17 +126,16 @@ namespace Opm {
rpt.norm_dx =
VNorm<vector_type>::norm(sys_.vector().increment());
// Begin line search. The line search should be moved to model
// Begin line search
double residual=VNorm<vector_type>::norm(sys_.vector().residual());
int lin_it=0;
bool finished=rpt.norm_res<ctrl.atol;//residual < rpt.norm_res;
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;
s_range = true;
sys_.vector().writableIncrement()=dx_old;
sys_.vector().writableIncrement()*=alpha;
sys_.vector().writableSolution()=x_old;
@ -147,20 +146,20 @@ namespace Opm {
return vec;
*/
sys_.vector().addIncrement();
model_.initIteration(state, g, sys_, s_range);
if (s_range) {
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.verbose){
if (ctrl.verbosity > 1){
std::cout << "Line search iteration " << std::scientific << lin_it
<< " norm :" << residual << " alpha " << alpha << '\n';
}
}else{
if (ctrl.verbose){
if (ctrl.verbosity > 1){
std::cout << "Line search iteration " << std::scientific << lin_it
<< " Saturation out of range, continue. Alpha " << alpha << '\n';
<< " Value out of range, continue search. alpha " << alpha << '\n';
}
residual = 1e99;
}
@ -171,7 +170,7 @@ namespace Opm {
VNorm<vector_type>::norm(sys_.vector().residual());
rpt.nit++;
if (ctrl.verbose){
if (ctrl.verbosity > 0){
std::cout << "Iteration " << std::scientific << rpt.nit
<< " norm :" << rpt.norm_res << " alpha " << alpha << std::endl;
}

3
src/SimpleFluid2pWrapper.hpp
View File

@ -36,7 +36,6 @@
#ifndef OPM_SimpleFluid2pWrapper_HPP_HEADER
#define OPM_SimpleFluid2pWrapper_HPP_HEADER
#include <array>
#include <cmath>
namespace Opm {
@ -67,8 +66,6 @@ namespace Opm {
private:
ReservoirProperties& resprop_;
std::array<double, 2> mu_ ;
std::array<double, 2> rho_;
};
}

17
src/SinglePointUpwindTwoPhase.hpp
View File

@ -361,11 +361,10 @@ namespace Opm {
template <class ReservoirState,
class Grid ,
class JacobianSystem>
void
bool
initIteration(const ReservoirState& state,
const Grid& g ,
JacobianSystem& sys,
bool& s_range ) {
JacobianSystem& sys) {
double s[2], mob[2], dmob[2 * 2], pc, dpc;
@ -373,7 +372,7 @@ namespace Opm {
sys.vector().solution();
const ::std::vector<double>& sat = state.saturation();
double max_alpha = 1;
bool in_range = true;
for (int c = 0; c < g.number_of_cells; ++c) {
store_.ds(c) = x[c]; // Store sat-change for accumulation().
@ -384,14 +383,13 @@ namespace Opm {
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;
std::cout << "Warning: s out of range, s-s_min = " << s_min-s[0] << std::endl;
}
if (s[0] > s_max){
std::cout << "Warning: s out of range:" << s[0]-s_max << std::endl;
std::cout << "Warning: s out of range, s-s_max = " << s[0]-s_max << std::endl;
}
s_range = false; //line search fails
in_range = false; //line search fails
}
s[0] = std::max(s_min, s[0]);
s[0] = std::min(s_max, s[0]);
s[1] = 1 - s[0];
@ -406,6 +404,7 @@ namespace Opm {
store_.pc(c) = pc;
store_.dpc(c) = dpc;
}
return in_range;
}
template <class ReservoirState,
@ -532,7 +531,7 @@ namespace Opm {
}
TwophaseFluid fluid_ ;
const double* gravity_;
const double* gravity_;
std::vector<int> f2hf_ ;
spu_2p::ModelParameterStorage store_ ;
};