diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index 565082e6f..62eb090a9 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -339,11 +339,11 @@ namespace Opm &transport_src[0], stepsize, inflow_c, state.saturation(), state.surfacevol(), state.concentration(), state.maxconcentration()); - - // Computeinjectedproduced function does not take into account polymer. - Opm::computeInjectedProduced(props_, + Opm::computeInjectedProduced(props_, poly_props_, state.pressure(), state.surfacevol(), state.saturation(), - transport_src, stepsize, injected, produced); + state.concentration(), state.maxconcentration(), + transport_src, stepsize, inflow_c, injected, produced, + polyinj, polyprod); if (gravity_ != 0 && use_segregation_split_) { tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol(), diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp index ad7f9ba5c..72b42f578 100644 --- a/opm/polymer/polymerUtilities.cpp +++ b/opm/polymer/polymerUtilities.cpp @@ -92,7 +92,7 @@ namespace Opm /// @brief Computes injected and produced volumes of all phases, - /// and injeced and produced polymer mass. + /// and injected and produced polymer mass. /// Note 1: assumes that only the first phase is injected. /// Note 2: assumes that transport has been done with an /// implicit method, i.e. that the current state @@ -101,6 +101,7 @@ namespace Opm /// @param[in] polyprops polymer properties /// @param[in] s saturation values (for all P phases) /// @param[in] c polymer concentration + /// @param[in] cmax polymer maximum concentration /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. /// @param[in] dt timestep used /// @param[in] inj_c injected concentration @@ -134,6 +135,7 @@ namespace Opm const double* visc = props.viscosity(); std::vector kr_cell(np); double mob[2]; + double mc; for (int cell = 0; cell < num_cells; ++cell) { if (src[cell] > 0.0) { injected[0] += src[cell]*dt; @@ -148,11 +150,89 @@ namespace Opm for (int p = 0; p < np; ++p) { produced[p] += (mob[p]/totmob)*flux; } - polyprod += (mob[0]/totmob)*flux*c[cell]; // TODO check this term. + polyprops.computeMc(c[cell], mc); + polyprod += (mob[0]/totmob)*flux*mc; } } } + /// @brief Computes injected and produced volumes of all phases, + /// and injected and produced polymer mass - in the compressible case. + /// Note 1: assumes that only the first phase is injected. + /// Note 2: assumes that transport has been done with an + /// implicit method, i.e. that the current state + /// gives the mobilities used for the preceding timestep. + /// @param[in] props fluid and rock properties. + /// @param[in] polyprops polymer properties + /// @param[in] press pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all P phases) + /// @param[in] c polymer concentration + /// @param[in] cmax polymer maximum concentration + /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. + /// @param[in] dt timestep used + /// @param[in] inj_c injected concentration + /// + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. + /// @param[out] polyinj injected mass of polymer + /// @param[out] polyprod produced mass of polymer + + void computeInjectedProduced(const BlackoilPropertiesInterface& props, + const Opm::PolymerProperties& polyprops, + const std::vector& press, + const std::vector& z, + const std::vector& s, + const std::vector& c, + const std::vector& cmax, + const std::vector& src, + const double dt, + const double inj_c, + double* injected, + double* produced, + double& polyinj, + double& polyprod) + { + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); + polyinj = 0.0; + polyprod = 0.0; + std::vector visc(np); + std::vector kr_cell(np); + std::vector mob(np); + double mc; + for (int cell = 0; cell < num_cells; ++cell) { + if (src[cell] > 0.0) { + injected[0] += src[cell]*dt; + polyinj += src[cell]*dt*inj_c; + } else if (src[cell] < 0.0) { + const double flux = -src[cell]*dt; + const double* sat = &s[np*cell]; + props.relperm(1, sat, &cell, &kr_cell[0], 0); + props.viscosity(1, &press[cell], &z[np*cell], &cell, &visc[0], 0); + polyprops.effectiveMobilities(c[cell], cmax[cell], &visc[0], + &kr_cell[0], &mob[0]); + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + polyprops.computeMc(c[cell], mc); + polyprod += (mob[0]/totmob)*flux*mc; + } + } + } + + + /// @brief Computes total polymer mass over all grid cells. /// @param[in] pv the pore volume by cell. diff --git a/opm/polymer/polymerUtilities.hpp b/opm/polymer/polymerUtilities.hpp index 0ae16a3fc..064b4808a 100644 --- a/opm/polymer/polymerUtilities.hpp +++ b/opm/polymer/polymerUtilities.hpp @@ -97,6 +97,44 @@ namespace Opm double& polyinj, double& polyprod); + /// @brief Computes injected and produced volumes of all phases, + /// and injected and produced polymer mass - in the compressible case. + /// Note 1: assumes that only the first phase is injected. + /// Note 2: assumes that transport has been done with an + /// implicit method, i.e. that the current state + /// gives the mobilities used for the preceding timestep. + /// @param[in] props fluid and rock properties. + /// @param[in] polyprops polymer properties + /// @param[in] press pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all P phases) + /// @param[in] c polymer concentration + /// @param[in] cmax polymer maximum concentration + /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. + /// @param[in] dt timestep used + /// @param[in] inj_c injected concentration + /// + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. + /// @param[out] polyinj injected mass of polymer + /// @param[out] polyprod produced mass of polymer + + void computeInjectedProduced(const BlackoilPropertiesInterface& props, + const Opm::PolymerProperties& polyprops, + const std::vector& press, + const std::vector& z, + const std::vector& s, + const std::vector& c, + const std::vector& cmax, + const std::vector& src, + const double dt, + const double inj_c, + double* injected, + double* produced, + double& polyinj, + double& polyprod); + /// @brief Computes total (free) polymer mass over all grid cells. /// @param[in] pv the pore volume by cell. /// @param[in] s saturation values (for all P phases)