Fix mass balance report.

Now reports surface volumes for water and oil, instead of useless
(for mass balance) reservoir volumes. Also makes the case with
multiple transport steps more accurate.
This commit is contained in:
Atgeirr Flø Rasmussen 2012-10-04 11:25:42 +02:00
parent c6b76e715d
commit b32c212e56

View File

@ -265,22 +265,16 @@ namespace Opm
double ttime = 0.0;
Opm::time::StopWatch total_timer;
total_timer.start();
double init_satvol[2] = { 0.0 };
double init_surfvol[2] = { 0.0 };
double init_polymass = 0.0;
double satvol[2] = { 0.0 };
double inplace_surfvol[2] = { 0.0 };
double polymass = 0.0;
double polymass_adsorbed = 0.0;
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double polyinj = 0.0;
double polyprod = 0.0;
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
double tot_polyinj = 0.0;
double tot_polyprod = 0.0;
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init
<< " " << init_satvol[1]/tot_porevol_init << std::endl;
Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
@ -339,17 +333,32 @@ namespace Opm
stepsize /= double(num_transport_substeps_);
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
}
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double polyinj = 0.0;
double polyprod = 0.0;
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], initial_pressure,
state.pressure(), &initial_porevol[0], &porevol[0],
&transport_src[0], stepsize, inflow_c,
state.saturation(), state.surfacevol(),
state.concentration(), state.maxconcentration());
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
double substep_polyinj = 0.0;
double substep_polyprod = 0.0;
Opm::computeInjectedProduced(props_, poly_props_,
state.pressure(), state.surfacevol(), state.saturation(),
state.concentration(), state.maxconcentration(),
transport_src, stepsize, inflow_c, injected, produced,
polyinj, polyprod);
transport_src, stepsize, inflow_c,
substep_injected, substep_produced,
substep_polyinj, substep_polyprod);
injected[0] += substep_injected[0];
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
polyinj += substep_polyinj;
polyprod += substep_polyprod;
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, stepsize,
state.saturation(), state.surfacevol(),
@ -362,7 +371,7 @@ namespace Opm
ttime += tt;
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
state, rock_comp_props_);
@ -374,40 +383,44 @@ namespace Opm
tot_polyprod += polyprod;
std::cout.precision(5);
const int width = 18;
std::cout << "\nVolume and polymer mass balance: "
" water(pv) oil(pv) polymer(kg)\n";
std::cout << " Saturated volumes: "
<< std::setw(width) << satvol[0]/tot_porevol_init
<< std::setw(width) << satvol[1]/tot_porevol_init
std::cout << "\nMass balance: "
" water(surfvol) oil(surfvol) polymer(kg)\n";
std::cout << " In-place: "
<< std::setw(width) << inplace_surfvol[0]
<< std::setw(width) << inplace_surfvol[1]
<< std::setw(width) << polymass << std::endl;
std::cout << " Adsorbed volumes: "
std::cout << " Adsorbed: "
<< std::setw(width) << 0.0
<< std::setw(width) << 0.0
<< std::setw(width) << polymass_adsorbed << std::endl;
std::cout << " Injected volumes: "
<< std::setw(width) << injected[0]/tot_porevol_init
<< std::setw(width) << injected[1]/tot_porevol_init
std::cout << " Injected: "
<< std::setw(width) << injected[0]
<< std::setw(width) << injected[1]
<< std::setw(width) << polyinj << std::endl;
std::cout << " Produced volumes: "
<< std::setw(width) << produced[0]/tot_porevol_init
<< std::setw(width) << produced[1]/tot_porevol_init
std::cout << " Produced: "
<< std::setw(width) << produced[0]
<< std::setw(width) << produced[1]
<< std::setw(width) << polyprod << std::endl;
std::cout << " Total inj volumes: "
<< std::setw(width) << tot_injected[0]/tot_porevol_init
<< std::setw(width) << tot_injected[1]/tot_porevol_init
std::cout << " Total inj: "
<< std::setw(width) << tot_injected[0]
<< std::setw(width) << tot_injected[1]
<< std::setw(width) << tot_polyinj << std::endl;
std::cout << " Total prod volumes: "
<< std::setw(width) << tot_produced[0]/tot_porevol_init
<< std::setw(width) << tot_produced[1]/tot_porevol_init
std::cout << " Total prod: "
<< std::setw(width) << tot_produced[0]
<< std::setw(width) << tot_produced[1]
<< std::setw(width) << tot_polyprod << std::endl;
std::cout << " In-place + prod - inj: "
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init
<< std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl;
std::cout << " Init - now - pr + inj: "
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
<< std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed)
const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1],
init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed };
std::cout << " Initial - inplace + inj - prod: "
<< std::setw(width) << balance[0]
<< std::setw(width) << balance[1]
<< std::setw(width) << balance[2]
<< std::endl;
std::cout << " Relative mass error: "
<< std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
<< std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
<< std::setw(width) << balance[2]/(init_polymass + tot_polyinj)
<< std::endl;
std::cout.precision(8);