Merge remote-tracking branch 'opm/master' into wpolymer_schedule

This commit is contained in:
Liu Ming 2014-12-03 16:56:38 +08:00
commit 242ec598c8
15 changed files with 86 additions and 34 deletions

View File

@ -17,6 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/CompressibleTpfaPolymer.hpp>
@ -130,12 +131,13 @@ namespace Opm
const int nc = grid_.number_of_cells;
const int np = props_.numPhases();
const double* cell_p = &state.pressure()[0];
const double* cell_T = &state.temperature()[0];
const double* cell_z = &state.surfacevol()[0];
cell_A_.resize(nc*np*np);
cell_dA_.resize(nc*np*np);
props_.matrix(nc, cell_p, cell_z, &allcells_[0], &cell_A_[0], &cell_dA_[0]);
props_.matrix(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_A_[0], &cell_dA_[0]);
cell_viscosity_.resize(nc*np);
props_.viscosity(nc, cell_p, cell_z, &allcells_[0], &cell_viscosity_[0], 0);
props_.viscosity(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_viscosity_[0], 0);
cell_phasemob_.resize(nc*np);
for (int cell = 0; cell < nc; ++cell) {
poly_props_.effectiveVisc((*c_)[cell], &cell_viscosity_[np*cell + 0], cell_eff_viscosity_[np*cell + 0]);

View File

@ -18,6 +18,8 @@
*/
#include <config.h>
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>

View File

@ -61,6 +61,7 @@ namespace Opm
}
std::vector<double>& pressure () { return state_blackoil_.pressure(); }
std::vector<double>& temperature () { return state_blackoil_.temperature(); }
std::vector<double>& surfacevol () { return state_blackoil_.surfacevol(); }
std::vector<double>& facepressure() { return state_blackoil_.facepressure(); }
std::vector<double>& faceflux () { return state_blackoil_.faceflux(); }
@ -71,6 +72,7 @@ namespace Opm
std::vector<double>& maxconcentration() { return cmax_; }
const std::vector<double>& pressure () const { return state_blackoil_.pressure(); }
const std::vector<double>& temperature () const { return state_blackoil_.temperature(); }
const std::vector<double>& surfacevol () const { return state_blackoil_.surfacevol(); }
const std::vector<double>& facepressure() const { return state_blackoil_.facepressure(); }
const std::vector<double>& faceflux () const { return state_blackoil_.faceflux(); }

View File

@ -17,6 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/wells.h>
#include <opm/parser/eclipse/Deck/Deck.hpp>

View File

@ -17,6 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/polymer/PolymerProperties.hpp>
#include <cmath>
#include <vector>

View File

@ -306,7 +306,7 @@ namespace Opm
// Solve pressure equation.
if (check_well_controls_) {
computeFractionalFlow(props_, poly_props_, allcells_,
state.pressure(), state.surfacevol(), state.saturation(),
state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
state.concentration(), state.maxconcentration(),
fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
@ -400,7 +400,7 @@ namespace Opm
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],
state.pressure(), state.temperature(), &initial_porevol[0], &porevol[0],
&transport_src[0], &polymer_inflow_c[0], stepsize,
state.saturation(), state.surfacevol(),
state.concentration(), state.maxconcentration());

View File

@ -17,6 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
@ -214,6 +215,7 @@ namespace Opm
void TransportSolverTwophaseCompressiblePolymer::solve(const double* darcyflux,
const std::vector<double>& initial_pressure,
const std::vector<double>& pressure,
const std::vector<double>& temperature,
const double* porevolume0,
const double* porevolume,
const double* source,
@ -238,9 +240,9 @@ namespace Opm
res_counts.clear();
#endif
props_.viscosity(grid_.number_of_cells, &pressure[0], NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(grid_.number_of_cells, &initial_pressure[0], NULL, &allcells_[0], &A0_[0], NULL);
props_.matrix(grid_.number_of_cells, &pressure[0], NULL, &allcells_[0], &A_[0], NULL);
props_.viscosity(grid_.number_of_cells, &pressure[0], &temperature[0], NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(grid_.number_of_cells, &initial_pressure[0], &temperature[0], NULL, &allcells_[0], &A0_[0], NULL);
props_.matrix(grid_.number_of_cells, &pressure[0], &temperature[0], NULL, &allcells_[0], &A_[0], NULL);
// Check immiscibility requirement (only done for first cell).
if (A_[1] != 0.0 || A_[2] != 0.0) {

View File

@ -74,6 +74,7 @@ namespace Opm
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] initial_pressure Array with pressure at start of timestep.
/// \param[in] pressure Array with pressure.
/// \param[in] temperature Array with temperature.
/// \param[in] porevolume0 Array with pore volume at start of timestep.
/// \param[in] porevolume Array with pore volume.
/// \param[in] source Transport source term, to be interpreted by sign:
@ -91,6 +92,7 @@ namespace Opm
void solve(const double* darcyflux,
const std::vector<double>& initial_pressure,
const std::vector<double>& pressure,
const std::vector<double>& temperature,
const double* porevolume0,
const double* porevolume,
const double* source,

View File

@ -18,6 +18,8 @@
*/
#include <config.h>
#include <opm/polymer/TransportSolverTwophasePolymer.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/grid.h>

View File

@ -127,6 +127,7 @@ namespace Opm {
struct SolutionState {
SolutionState(const int np);
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB rs;
ADB rv;
@ -277,6 +278,7 @@ namespace Opm {
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -285,6 +287,7 @@ namespace Opm {
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -293,6 +296,7 @@ namespace Opm {
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,

View File

@ -365,6 +365,7 @@ namespace {
template<class T>
FullyImplicitBlackoilPolymerSolver<T>::SolutionState::SolutionState(const int np)
: pressure ( ADB::null())
, temperature( ADB::null())
, saturation(np, ADB::null())
, rs ( ADB::null())
, rv ( ADB::null())
@ -422,6 +423,7 @@ namespace {
// automatically consistent with variableState() (and doing
// things automatically is all the rage in this module ;)
state.pressure = ADB::constant(state.pressure.value());
state.temperature = ADB::constant(state.temperature.value());
state.rs = ADB::constant(state.rs.value());
state.rv = ADB::constant(state.rv.value());
state.concentration = ADB::constant(state.concentration.value());
@ -527,6 +529,10 @@ namespace {
int nextvar = 0;
state.pressure = vars[ nextvar++ ];
// Temperature.
const V temp = Eigen::Map<const V>(& x.temperature()[0], nc, 1);
state.pressure = ADB::constant(temp);
// Saturations
const std::vector<int>& bpat = vars[0].blockPattern();
{
@ -592,6 +598,7 @@ namespace {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const std::vector<ADB>& sat = state.saturation;
const ADB& rs = state.rs;
const ADB& rv = state.rv;
@ -606,7 +613,7 @@ namespace {
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
const int pos = pu.phase_pos[ phase ];
rq_[pos].b = fluidReciprocFVF(phase, pressure[phase], rs, rv, cond, cells_);
rq_[pos].b = fluidReciprocFVF(phase, pressure[phase], temp, rs, rv, cond, cells_);
rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
// DUMP(rq_[pos].b);
// DUMP(rq_[pos].accum[aix]);
@ -666,6 +673,7 @@ namespace {
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
// Compute b, rsmax, rvmax values for perforations.
const ADB perf_press = subset(state.pressure, well_cells);
const ADB perf_temp = subset(state.temperature, well_cells);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
@ -676,21 +684,21 @@ namespace {
std::vector<double> rssat_perf(nperf, 0.0);
std::vector<double> rvsat_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const ADB bw = fluid_.bWat(perf_press, well_cells);
const ADB bw = fluid_.bWat(perf_press, perf_temp, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
}
assert(active_[Oil]);
const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_temp, perf_rs, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells);
rssat_perf.assign(rssat.data(), rssat.data() + nperf);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_temp, perf_rv, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells);
rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
@ -1688,10 +1696,10 @@ namespace {
const std::vector<PhasePresence> cond = phaseCondition();
const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond, cells_);
rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond, cells_);
ADB& head = rq_[phase].head;
@ -1984,6 +1992,7 @@ namespace {
ADB
FullyImplicitBlackoilPolymerSolver<T>::fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -1991,12 +2000,12 @@ namespace {
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, temp, cells);
case Oil: {
return fluid_.muOil(p, rs, cond, cells);
return fluid_.muOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.muGas(p, rv, cond, cells);
return fluid_.muGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -2010,6 +2019,7 @@ namespace {
ADB
FullyImplicitBlackoilPolymerSolver<T>::fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -2017,12 +2027,12 @@ namespace {
{
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, temp, cells);
case Oil: {
return fluid_.bOil(p, rs, cond, cells);
return fluid_.bOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.bGas(p, rv, cond, cells);
return fluid_.bGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -2036,13 +2046,14 @@ namespace {
ADB
FullyImplicitBlackoilPolymerSolver<T>::fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidReciprocFVF(phase, p, rs, rv, cond, cells);
ADB b = fluidReciprocFVF(phase, p, temp, rs, rv, cond, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
if (phase == Oil && active_[Gas]) {
// It is correct to index into rhos with canonical phase indices.

View File

@ -263,6 +263,7 @@ namespace {
FullyImplicitCompressiblePolymerSolver::SolutionState::SolutionState(const int np)
: pressure ( ADB::null())
, temperature( ADB::null())
, saturation(np, ADB::null())
, concentration( ADB::null())
, qs ( ADB::null())
@ -331,6 +332,10 @@ namespace {
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
state.pressure = ADB::constant(p, bpat);
// Temperature.
const V T = Eigen::Map<const V>(& x.temperature()[0], nc, 1);
state.temperature = ADB::constant(T);
// Saturation.
assert (not x.saturation().empty());
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
@ -449,6 +454,7 @@ namespace {
{
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const std::vector<ADB>& sat = state.saturation;
const ADB& c = state.concentration;
@ -458,7 +464,7 @@ namespace {
const ADB pv_mult = poroMult(press);
for (int phase = 0; phase < 2; ++phase) {
rq_[phase].b = fluidReciprocFVF(phase, pressure[phase], cond, cells_);
rq_[phase].b = fluidReciprocFVF(phase, pressure[phase], temp, cond, cells_);
}
rq_[0].accum[aix] = pv_mult * rq_[0].b * sat[0];
rq_[1].accum[aix] = pv_mult * rq_[1].b * sat[1];
@ -566,16 +572,17 @@ namespace {
}
ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern());
std::vector<ADB> press = computePressures(state);
const ADB& temp = state.temperature;
const std::vector<PhasePresence> cond = phaseCondition();
for (int phase = 0; phase < 2; ++phase) {
const ADB cell_rho = fluidDensity(phase, press[phase], cond, cells_);
const ADB cell_rho = fluidDensity(phase, press[phase], temp, cond, cells_);
cell_rho_total += state.saturation[phase] * cell_rho;
}
ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern());
assert(np == wells_.number_of_phases);
const DataBlock compi = Eigen::Map<const DataBlock>(wells_.comp_frac, nw, np);
for (int phase = 0; phase < 2; ++phase) {
const ADB cell_rho = fluidDensity(phase, press[phase], cond, cells_);
const ADB cell_rho = fluidDensity(phase, press[phase], temp, cond, cells_);
const V fraction = compi.col(phase);
inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells);
}
@ -840,15 +847,16 @@ namespace {
const ADB tr_mult = transMult(state.pressure);
const std::vector<PhasePresence> cond = phaseCondition();
std::vector<ADB> press = computePressures(state);
const ADB& temp = state.temperature;
const ADB mu_w = fluidViscosity(0, press[0], cond, cells_);
const ADB mu_w = fluidViscosity(0, press[0], temp, cond, cells_);
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu_w.value().data());
rq_[0].mob = tr_mult * krw_eff * inv_wat_eff_vis;
rq_[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis;
const ADB mu_o = fluidViscosity(1, press[1], cond, cells_);
const ADB mu_o = fluidViscosity(1, press[1], temp, cond, cells_);
rq_[1].mob = tr_mult * kro / mu_o;
for (int phase = 0; phase < 2; ++phase) {
const ADB rho = fluidDensity(phase, press[phase], cond, cells_);
const ADB rho = fluidDensity(phase, press[phase], temp, cond, cells_);
ADB& head = rq_[ phase ].head;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rho;
@ -893,15 +901,16 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::fluidViscosity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern());
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, T, cells);
case Oil: {
return fluid_.muOil(p, null, cond, cells);
return fluid_.muOil(p, T, null, cond, cells);
}
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
@ -915,15 +924,16 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern());
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, T, cells);
case Oil: {
return fluid_.bOil(p, null, cond, cells);
return fluid_.bOil(p, T, null, cond, cells);
}
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
@ -937,11 +947,12 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::fluidDensity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidReciprocFVF(phase, p, cond, cells);
ADB b = fluidReciprocFVF(phase, p, T, cond, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
return rho;
}

View File

@ -107,6 +107,7 @@ namespace Opm {
struct SolutionState {
SolutionState(const int np);
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB concentration;
ADB qs;
@ -213,18 +214,21 @@ namespace Opm {
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;

View File

@ -17,6 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/polymer/polymerUtilities.hpp>
#include <opm/core/utility/miscUtilities.hpp>
@ -141,6 +143,7 @@ namespace Opm
const Opm::PolymerProperties& polyprops,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& c,
@ -157,7 +160,7 @@ namespace Opm
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
std::vector<double> mu(num_cells*num_phases);
props.viscosity(num_cells, &p[0], &z[0], &cells[0], &mu[0], 0);
props.viscosity(num_cells, &p[0], &T[0], &z[0], &cells[0], &mu[0], 0);
double mob[2]; // here we assume num_phases=2
for (int cell = 0; cell < num_cells; ++cell) {
double* kr_cell = &kr[2*cell];
@ -276,6 +279,7 @@ namespace Opm
OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
}
const std::vector<double>& press = state.pressure();
const std::vector<double>& temp = state.temperature();
const std::vector<double>& s = state.saturation();
const std::vector<double>& z = state.surfacevol();
const std::vector<double>& c = state.concentration();
@ -303,8 +307,8 @@ namespace Opm
const double flux = -transport_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);
props.matrix(1, &press[cell], &z[np*cell], &cell, &A[0], 0);
props.viscosity(1, &press[cell], &temp[cell], &z[np*cell], &cell, &visc[0], 0);
props.matrix(1, &press[cell], &temp[cell], &z[np*cell], &cell, &A[0], 0);
polyprops.effectiveMobilities(c[cell], cmax[cell], &visc[0],
&kr_cell[0], &mob[0]);
double totmob = 0.0;

View File

@ -90,6 +90,7 @@ namespace Opm
/// @param[in] polyprops polymer properties
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell)
/// @param[in] T temperature (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases)
/// @param[in] c concentration values
@ -99,6 +100,7 @@ namespace Opm
const Opm::PolymerProperties& polyprops,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& T,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& c,