From 076164d3f3a9c2fcf35ee6bbcbddba717182f08c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 17 Sep 2012 07:54:50 +0200 Subject: [PATCH 01/10] Add "output_vtk" parameter, split output function in two. Vtk and Matlab output now happens in two different functions. --- opm/polymer/SimulatorPolymer.cpp | 77 +++++++++++++++++++++++++++----- 1 file changed, 65 insertions(+), 12 deletions(-) diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index d48067b52..f5e7a7d0e 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -64,15 +64,18 @@ namespace Opm namespace { - void outputState(const UnstructuredGrid& grid, - const Opm::PolymerState& state, - const int step, - const std::string& output_dir); + void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir); + void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir); void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir); void outputWellReport(const Opm::WellReport& wellreport, const std::string& output_dir); - } // anonymous namespace @@ -100,6 +103,7 @@ namespace Opm // Parameters for output. bool output_; + bool output_vtk_; std::string output_dir_; int output_interval_; // Parameters for transport solver. @@ -191,6 +195,7 @@ namespace Opm // For output. output_ = param.getDefault("output", true); if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); output_dir_ = param.getDefault("output_dir", std::string("output")); // Ensure that output dir exists boost::filesystem::path fpath(output_dir_); @@ -283,7 +288,10 @@ namespace Opm // Report timestep and (optionally) write state to disk. timer.report(std::cout); if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - outputState(grid_, state, timer.currentStepNum(), output_dir_); + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); } // Solve pressure. @@ -398,7 +406,10 @@ namespace Opm } if (output_) { - outputState(grid_, state, timer.currentStepNum(), output_dir_); + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); outputWaterCut(watercut, output_dir_); if (wells_) { outputWellReport(wellreport, output_dir_); @@ -423,14 +434,22 @@ namespace Opm namespace { - void outputState(const UnstructuredGrid& grid, - const Opm::PolymerState& state, - const int step, - const std::string& output_dir) + void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir) { // Write data in VTK format. std::ostringstream vtkfilename; - vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); @@ -458,6 +477,40 @@ namespace Opm } } + void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["concentration"] = &state.concentration(); + dm["cmax"] = &state.maxconcentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir) From 678611166d8288a526ec9a30db863a614b1b7446 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 14:25:20 +0200 Subject: [PATCH 02/10] Improve reporting of volume balance. Should now give more accurate results when taking multiple substeps. --- opm/polymer/SimulatorPolymer.cpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index f5e7a7d0e..88013e73e 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -329,13 +329,24 @@ namespace Opm stepsize /= double(num_transport_substeps_); std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; } + double substep_injected[2] = { 0.0 }; + double substep_produced[2] = { 0.0 }; + double substep_polyinj = 0.0; + double substep_polyprod = 0.0; + injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0; for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], stepsize, inflow_c, state.saturation(), state.concentration(), state.maxconcentration()); Opm::computeInjectedProduced(props_, poly_props_, state.saturation(), state.concentration(), state.maxconcentration(), - transport_src, timer.currentStepLength(), 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 (use_segregation_split_) { tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation(), state.concentration(), state.maxconcentration()); From 349c81fb1e57ad356b82bab0777f779876d98e76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 15:59:18 +0200 Subject: [PATCH 03/10] Stop writing ascii output twice. --- opm/polymer/SimulatorPolymer.cpp | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index 88013e73e..83e759d56 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -474,18 +474,6 @@ namespace Opm Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; Opm::writeVtkData(grid, dm, vtkfile); - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; - std::ofstream file(fname.str().c_str()); - if (!file) { - THROW("Failed to open " << fname.str()); - } - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } } void outputStateMatlab(const UnstructuredGrid& grid, From 2990b0a3e8dea32f47203a5dea348b374f404f3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 16:00:13 +0200 Subject: [PATCH 04/10] Improve function comment. --- opm/polymer/SimulatorPolymer.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/polymer/SimulatorPolymer.hpp b/opm/polymer/SimulatorPolymer.hpp index 448ae7319..091d386bf 100644 --- a/opm/polymer/SimulatorPolymer.hpp +++ b/opm/polymer/SimulatorPolymer.hpp @@ -83,7 +83,8 @@ namespace Opm /// This will run succesive timesteps until timer.done() is true. It will /// modify the reservoir and well states. /// \param[in,out] timer governs the requested reporting timesteps - /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] state state of reservoir: pressure, fluxes, polymer concentration, + /// saturations. /// \param[in,out] well_state state of wells: bhp, perforation rates /// \return simulation report, with timing data SimulatorReport run(SimulatorTimer& timer, From 753a0c69498944a581c17863d0072e634e706abe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 16:00:51 +0200 Subject: [PATCH 05/10] Reinstate 10% ad-hoc c-interval overestimation. --- opm/polymer/TransportModelPolymer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index 84361766a..5227f5411 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -196,7 +196,7 @@ namespace Opm fractionalflow_(grid.number_of_cells, -1.0), mc_(grid.number_of_cells, -1.0), method_(method), - adhoc_safety_(1.0) + adhoc_safety_(1.1) { if (props.numPhases() != 2) { THROW("Property object must have 2 phases"); From dc3b0dd6f471724a23f2ed90d01eeaf8bcb11ec9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 16:01:54 +0200 Subject: [PATCH 06/10] Improve function comment (similar to in SimulatorPolymer). --- opm/polymer/SimulatorCompressiblePolymer.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/polymer/SimulatorCompressiblePolymer.hpp b/opm/polymer/SimulatorCompressiblePolymer.hpp index cd94e46f4..cff17c48c 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.hpp +++ b/opm/polymer/SimulatorCompressiblePolymer.hpp @@ -83,7 +83,8 @@ namespace Opm /// This will run succesive timesteps until timer.done() is true. It will /// modify the reservoir and well states. /// \param[in,out] timer governs the requested reporting timesteps - /// \param[in,out] state state of reservoir: pressure, fluxes, polymer concentration + /// \param[in,out] state state of reservoir: pressure, fluxes, polymer concentration, + /// saturations, surface volumes. /// \param[in,out] well_state state of wells: bhp, perforation rates /// \return simulation report, with timing data SimulatorReport run(SimulatorTimer& timer, From 25e528df0dcf48119896dc6b8da65b9b0d7ecd4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 16:02:27 +0200 Subject: [PATCH 07/10] Make output work the same as in other sims. I.e. putting each field in its own directory, adding the output_vtk parameter etc. --- opm/polymer/SimulatorCompressiblePolymer.cpp | 75 +++++++++++++++----- 1 file changed, 58 insertions(+), 17 deletions(-) diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index 62eb090a9..8702d76f3 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -66,10 +66,14 @@ namespace Opm namespace { - void outputState(const UnstructuredGrid& grid, - const Opm::PolymerBlackoilState& state, - const int step, - const std::string& output_dir); + void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir); + void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir); void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir); void outputWellReport(const Opm::WellReport& wellreport, @@ -102,6 +106,7 @@ namespace Opm // Parameters for output. bool output_; + bool output_vtk_; std::string output_dir_; int output_interval_; // Parameters for transport solver. @@ -194,6 +199,7 @@ namespace Opm // For output. output_ = param.getDefault("output", true); if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); output_dir_ = param.getDefault("output_dir", std::string("output")); // Ensure that output dir exists boost::filesystem::path fpath(output_dir_); @@ -289,13 +295,13 @@ namespace Opm // Report timestep and (optionally) write state to disk. timer.report(std::cout); if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - outputState(grid_, state, timer.currentStepNum(), output_dir_); - } - - if (rock_comp_props_ && rock_comp_props_->isActive()) { - initial_pressure = state.pressure(); + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); } + initial_pressure = state.pressure(); // Solve pressure. do { @@ -416,7 +422,10 @@ namespace Opm } if (output_) { - outputState(grid_, state, timer.currentStepNum(), output_dir_); + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); outputWaterCut(watercut, output_dir_); if (wells_) { outputWellReport(wellreport, output_dir_); @@ -441,14 +450,22 @@ namespace Opm namespace { - void outputState(const UnstructuredGrid& grid, - const Opm::PolymerBlackoilState& state, - const int step, - const std::string& output_dir) + void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir) { // Write data in VTK format. std::ostringstream vtkfilename; - vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); @@ -458,15 +475,40 @@ namespace Opm dm["pressure"] = &state.pressure(); dm["concentration"] = &state.concentration(); dm["cmax"] = &state.maxconcentration(); + dm["surfvol"] = &state.surfacevol(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; Opm::writeVtkData(grid, dm, vtkfile); + } + + void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["concentration"] = &state.concentration(); + dm["cmax"] = &state.maxconcentration(); + dm["surfvol"] = &state.surfacevol(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; // Write data (not grid) in Matlab format for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { std::ostringstream fname; - fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; std::ofstream file(fname.str().c_str()); if (!file) { THROW("Failed to open " << fname.str()); @@ -476,7 +518,6 @@ namespace Opm } } - void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir) { From dbd4c9e7b5bc79cde21db3e5779abdc60296935b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 16:16:08 +0200 Subject: [PATCH 08/10] Reinstate 10% ad-hoc c-interval overestimation. Similar to in TransportModelPolymer. --- opm/polymer/TransportModelCompressiblePolymer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/polymer/TransportModelCompressiblePolymer.cpp b/opm/polymer/TransportModelCompressiblePolymer.cpp index 74cd8f631..755e357e9 100644 --- a/opm/polymer/TransportModelCompressiblePolymer.cpp +++ b/opm/polymer/TransportModelCompressiblePolymer.cpp @@ -173,7 +173,7 @@ namespace Opm tol_(tol), maxit_(maxit), method_(method), - adhoc_safety_(1.0), + adhoc_safety_(1.1), concentration_(0), cmax_(0), fractionalflow_(grid.number_of_cells, -1.0), From c6b76e715d129ea64eaa6497447dd085854560d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Oct 2012 16:17:54 +0200 Subject: [PATCH 09/10] Fix for case with incompressible rock. Made rock comp argument optional in computePolymerAdsorbed(). Check inside func for active rock comp. --- opm/polymer/SimulatorCompressiblePolymer.cpp | 2 +- opm/polymer/polymerUtilities.cpp | 8 ++++++-- opm/polymer/polymerUtilities.hpp | 4 ++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index 8702d76f3..efb09ea61 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -365,7 +365,7 @@ namespace Opm Opm::computeSaturatedVol(porevol, state.saturation(), satvol); polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_, - state, *rock_comp_props_); + state, rock_comp_props_); tot_injected[0] += injected[0]; tot_injected[1] += injected[1]; tot_produced[0] += produced[0]; diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp index 72b42f578..910556a7c 100644 --- a/opm/polymer/polymerUtilities.cpp +++ b/opm/polymer/polymerUtilities.cpp @@ -294,13 +294,17 @@ namespace Opm const BlackoilPropertiesInterface& props, const Opm::PolymerProperties& polyprops, const PolymerBlackoilState& state, - const RockCompressibility& rock_comp + const RockCompressibility* rock_comp ) { const int num_cells = props.numCells(); const double rhor = polyprops.rockDensity(); std::vector porosity; - computePorosity(grid, props.porosity(), rock_comp, state.pressure(), porosity); + if (rock_comp && rock_comp->isActive()) { + computePorosity(grid, props.porosity(), *rock_comp, state.pressure(), porosity); + } else { + porosity.assign(props.porosity(), props.porosity() + num_cells); + } double abs_mass = 0.0; const std::vector& cmax = state.maxconcentration(); for (int cell = 0; cell < num_cells; ++cell) { diff --git a/opm/polymer/polymerUtilities.hpp b/opm/polymer/polymerUtilities.hpp index 662139fe8..66729af4c 100644 --- a/opm/polymer/polymerUtilities.hpp +++ b/opm/polymer/polymerUtilities.hpp @@ -163,13 +163,13 @@ namespace Opm /// @param[in] props fluid and rock properties. /// @param[in] polyprops polymer properties /// @param[in] state State variables - /// @param[in] rock_comp Rock compressibility + /// @param[in] rock_comp Rock compressibility (optional) /// @return total absorbed polymer mass. double computePolymerAdsorbed(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const Opm::PolymerProperties& polyprops, const PolymerBlackoilState& state, - const RockCompressibility& rock_comp); + const RockCompressibility* rock_comp); /// @brief Functor giving the injected amount of polymer as a function of time. From b32c212e5615f472f606bc54948a0cd40cd64b67 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 4 Oct 2012 11:25:42 +0200 Subject: [PATCH 10/10] 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. --- opm/polymer/SimulatorCompressiblePolymer.cpp | 89 +++++++++++--------- 1 file changed, 51 insertions(+), 38 deletions(-) diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index efb09ea61..8e07c55a3 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -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);