diff --git a/opm/polymer/fullyimplicit/.utilities.cpp.swp b/opm/polymer/fullyimplicit/.utilities.cpp.swp deleted file mode 100644 index b7460437f..000000000 Binary files a/opm/polymer/fullyimplicit/.utilities.cpp.swp and /dev/null differ diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index c3e99d1b2..f24ac775a 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -150,7 +150,8 @@ namespace { step(const double dt, PolymerState& x, WellState& xw, - const std::vector& polymer_inflow) + const std::vector& polymer_inflow, + std::vector& src) { V pvol(grid_.number_of_cells); @@ -169,7 +170,7 @@ namespace { const double atol = 1.0e-12; const double rtol = 5.0e-8; const int maxit = 40; - assemble(pvdt, old_state, x, xw, polymer_inflow); + assemble(pvdt, x, xw, polymer_inflow, src); const double r0 = residualNorm(); int it = 0; @@ -181,7 +182,7 @@ namespace { const V dx = solveJacobianSystem(); updateState(dx, x, xw); - assemble(pvdt, old_state, x, xw, polymer_inflow); + assemble(pvdt, x, xw, polymer_inflow, src); const double r = residualNorm(); @@ -393,10 +394,10 @@ namespace { void FullyImplicitTwophasePolymerSolver:: assemble(const V& pvdt, - const SolutionState& old_state, const PolymerState& x, const WellState& xw, - const std::vector& polymer_inflow) + const std::vector& polymer_inflow, + std::vector& src) { // Create the primary variables. const SolutionState state = variableState(x, xw); @@ -404,27 +405,10 @@ namespace { // -------- Mass balance equations for water and oil -------- const V trans = subset(transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); - const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); - // const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); const ADB mc = computeMc(state); computeMassFlux(trans, mc, kr[1], krw_eff, state); - //const std::vector source = accumSource(kr[1], krw_eff, state.concentration, src, polymer_inflow); - // const std::vector source = polymerSource(); - // const double rho_r = polymer_props_ad_.rockDensity(); -// const V phi = V::Constant(pvdt.size(), 1, *fluid_.porosity()); - -// const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); -// residual_.mass_balance[0] = pvdt * (state.saturation[0] - old_state.saturation[0]) -// + ops_.div * mflux[0]; -// residual_.mass_balance[1] = pvdt * (state.saturation[1] - old_state.saturation[1]) -// + ops_.div * mflux[1]; - // Mass balance equation for polymer -// residual_.mass_balance[2] = pvdt * (state.saturation[0] * state.concentration -// - old_state.saturation[0] * old_state.concentration) * (1. - dead_pore_vol) -// + pvdt * rho_r * (1. - phi) / phi * ads - // + ops_.div * mflux[2]; residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0]) + ops_.div*rq_[0].mflux; residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0]) @@ -516,14 +500,17 @@ namespace { well_contribs[phase] = superset(perf_flux, well_cells, nc); // DUMP(well_contribs[phase]); residual_.mass_balance[phase] += well_contribs[phase]; + for (int i = 0; i < nc; ++i) { + src[i] += well_contribs[phase].value()[i]; + } } // well rates contribs to polymer mass balance eqn. // for injection wells. const V polyin = Eigen::Map(& polymer_inflow[0], nc); const V poly_in_perf = subset(polyin, well_cells); - const V poly_c_cell = subset(state.concentration, well_cells).value(); - const V poly_c = producer.select(poly_c_cell, poly_in_perf); + const V poly_mc_cell = subset(mc, well_cells).value(); + const V poly_c = producer.select(poly_mc_cell, poly_in_perf); residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_c, well_cells, nc); // Set the well flux equation diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp index d96b19750..525b80a7d 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -29,7 +29,8 @@ namespace Opm { void step(const double dt, PolymerState& state, WellState& well_state, - const std::vector& polymer_inflow); + const std::vector& polymer_inflow, + std::vector& src); private: typedef AutoDiffBlock ADB; typedef ADB::V V; @@ -88,10 +89,10 @@ namespace Opm { const WellState& xw); void assemble(const V& pvdt, - const SolutionState& old_state, const PolymerState& x, const WellState& xw, - const std::vector& polymer_inflow); + const std::vector& polymer_inflow, + std::vector& src); V solveJacobianSystem() const; void updateState(const V& dx, PolymerState& x, diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp index 34d19a8ea..291207b08 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -146,6 +147,7 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); dm["concentration"] = &state.concentration(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); @@ -162,6 +164,7 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); dm["concentration"] = &state.concentration(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); @@ -191,7 +194,6 @@ namespace Opm -/* static void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir) { @@ -203,6 +205,7 @@ namespace Opm } watercut.write(os); } +/* static void outputWellReport(const Opm::WellReport& wellreport, const std::string& output_dir) { @@ -298,8 +301,9 @@ namespace Opm std::vector porevol; Opm::computePorevolume(grid_, props_.porosity(), porevol); - // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); std::vector polymer_inflow_c(grid_.number_of_cells); + std::vector transport_src(grid_.number_of_cells); // Main simulation loop. Opm::time::StopWatch solver_timer; @@ -307,6 +311,10 @@ namespace Opm Opm::time::StopWatch step_timer; Opm::time::StopWatch total_timer; total_timer.start(); + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); #if 0 // These must be changed for three-phase. double init_surfvol[2] = { 0.0 }; @@ -333,7 +341,7 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); + // outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); } @@ -348,7 +356,7 @@ namespace Opm polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); solver_timer.start(); std::vector initial_pressure = state.pressure(); - solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); + solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src); // Stop timer and report. solver_timer.stop(); @@ -378,6 +386,38 @@ namespace Opm } } } while (!well_control_passed); + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double polyinj = 0; + double polyprod = 0; + + Opm::computeInjectedProduced(props_, polymer_props_, + state, + transport_src, polymer_inflow_c, timer.currentStepLength(), + injected, produced, + polyinj, polyprod); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance report.\n"; + std::cout << " Injected reservoir volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced reservoir volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj reservoir volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod reservoir volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; // Update pore volumes if rock is compressible. @@ -440,9 +480,9 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); -#if 0 outputWaterCut(watercut, output_dir_); +#if 0 + outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); if (wells_) { outputWellReport(wellreport, output_dir_); } diff --git a/opm/polymer/fullyimplicit/utilities.cpp b/opm/polymer/fullyimplicit/utilities.cpp index bb03ff1cd..2008795a0 100644 --- a/opm/polymer/fullyimplicit/utilities.cpp +++ b/opm/polymer/fullyimplicit/utilities.cpp @@ -144,7 +144,7 @@ namespace Opm const V inv_muw_eff = polymer_props.effectiveInvWaterVisc(c, mus); std::vector mob(np); mob[0] = krw_eff * inv_muw_eff; - mob[1] = kr[1] / mu[1]; + mob[1] = kr[1] / mus[1]; const V watmob_c = src_selector.select(mob[0], one); const V oilmob_c = src_selector.select(mob[1], zero);