diff --git a/opm/polymer/fullyimplicit/utilities.cpp b/opm/polymer/fullyimplicit/utilities.cpp index cec4ed7df..bb03ff1cd 100644 --- a/opm/polymer/fullyimplicit/utilities.cpp +++ b/opm/polymer/fullyimplicit/utilities.cpp @@ -3,7 +3,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -80,6 +82,94 @@ namespace Opm } } + /// @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] state state variables (pressure, fluxes etc.) + /// @param[in] transport_src if < 0: total reservoir volume outflow, + /// if > 0: first phase *surface volume* inflow. + /// @param[in] inj_c injected concentration by cell + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/transport_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 IncompPropsAdInterface& props, + const Opm::PolymerPropsAd& polymer_props, + const PolymerState& state, + const std::vector& transport_src, + const std::vector& inj_c, + const double dt, + double* injected, + double* produced, + double& polyinj, + double& polyprod) + { + const int num_cells = transport_src.size(); + if (props.numCells() != num_cells) { + OPM_THROW(std::runtime_error, "Size of transport_src vector does not match number of cells in props."); + } + const int np = props.numPhases(); + if (int(state.saturation().size()) != num_cells*np) { + OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells."); + } + std::vector cells(num_cells); + const V p = Eigen::Map(&state.pressure()[0], num_cells, 1); + const DataBlock s = Eigen::Map(&state.saturation()[0], num_cells, np); + const V sw = s.col(0); + const V so = s.col(1); + const V c = Eigen::Map(&state.concentration()[0], num_cells, 1); + const V cmax = Eigen::Map(&state.maxconcentration()[0], num_cells, 1); + const V trans_src = Eigen::Map(&transport_src[0], num_cells, 1); + V src = V::Constant(num_cells, -1.0); // negative is injec, positive is producer. + for (int cell = 0; cell < num_cells; ++cell) { + cells[cell] = cell; + if(transport_src[cell] > 0.0) { + src[cell] = 1.0; + } + } + const Selector src_selector(src); + const V one = V::Constant(num_cells, 1.0); + const V zero = V::Zero(num_cells); + const std::vector kr = props.relperm(sw, so, cells); + + const V krw_eff = polymer_props.effectiveRelPerm(c, cmax, kr[0]); + const double* mus = props.viscosity(); + 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]; + + const V watmob_c = src_selector.select(mob[0], one); + const V oilmob_c = src_selector.select(mob[1], zero); + const V flux = trans_src * dt; + const V totmob_c = watmob_c + oilmob_c; + const V wat_src = flux * (watmob_c / totmob_c); + const V oil_src = flux * (oilmob_c / totmob_c); + const V mc = polymer_props.polymerWaterVelocityRatio(c); + + polyinj = 0.0; + polyprod = 0.0; + std::fill(injected, injected + np , 0.0); + std::fill(produced, produced + np , 0.0); + for (int cell = 0; cell < num_cells; ++cell) { + if (wat_src[cell] < 0) { + injected[0] += wat_src[cell]; + polyinj += injected[0] * inj_c[cell]; + } else { + produced[0] += wat_src[cell]; + produced[1] += oil_src[cell]; + polyprod += produced[0] * mc[cell]; + } + } + } + /// @brief Computes injected and produced volumes of all phases, /// and injected and produced polymer mass - in the compressible case. diff --git a/opm/polymer/fullyimplicit/utilities.hpp b/opm/polymer/fullyimplicit/utilities.hpp index 685396d7b..1b74b83df 100644 --- a/opm/polymer/fullyimplicit/utilities.hpp +++ b/opm/polymer/fullyimplicit/utilities.hpp @@ -5,7 +5,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -46,6 +48,35 @@ namespace Opm const WellState& well_state, std::vector& transport_src); + /// @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] state state variables (pressure, fluxes etc.) + /// @param[in] transport_src if < 0: total reservoir volume outflow, + /// if > 0: first phase *surface volume* inflow. + /// @param[in] inj_c injected concentration by cell + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/transport_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 IncompPropsAdInterface& props, + const Opm::PolymerPropsAd& polymer_props, + const PolymerState& state, + const std::vector& transport_src, + const std::vector& inj_c, + const double dt, + double* injected, + double* produced, + double& polyinj, + double& polyprod); + /// @brief Computes injected and produced volumes of all phases, /// and injected and produced polymer mass - in the compressible case.