OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

since two phase polymer simulator support EclipseWriter, need to remove historical

Browse Source 
files for outputing water cut.
This commit is contained in:
Liu Ming 2014-10-29 13:38:15 +08:00
parent ac764bd8a8
commit 58bdc701e3
5 changed files with 17 additions and 457 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -197,9 +197,8 @@ namespace {
FullyImplicitCompressiblePolymerSolver::
step(const double dt,
PolymerBlackoilState& x ,
WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow,
std::vector<double>& src)
WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow)
{
const SolutionState state = constantState(x, xw);
@ -209,7 +208,7 @@ namespace {
const double atol = 1.0e-12;
const double rtol = 5.0e-8;
const int maxit = 15;
assemble(dt, x, xw, polymer_inflow, src);
assemble(dt, x, xw, polymer_inflow);
const double r0 = residualNorm();
const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm<Eigen::Infinity>();
@ -223,7 +222,7 @@ namespace {
const V dx = solveJacobianSystem();
updateState(dx, x, xw);
assemble(dt, x, xw, polymer_inflow, src);
assemble(dt, x, xw, polymer_inflow);
const double r = residualNorm();
@ -494,10 +493,9 @@ namespace {
void
FullyImplicitCompressiblePolymerSolver::
assemble(const double dt,
const PolymerBlackoilState& x ,
const WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow,
std::vector<double>& src)
const PolymerBlackoilState& x,
const WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow)
{
// Create the primary variables.
//
@ -541,9 +539,6 @@ namespace {
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
for (int i = 0; i < nc; ++i) {
src[i] = 0.0;
}
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
const V transw = Eigen::Map<const V>(wells_.WI, nperf);
@ -621,9 +616,6 @@ namespace {
well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc);
// DUMP(well_contribs[phase]);
residual_.material_balance_eq[phase] += well_contribs[phase];
for (int cell = 0; cell < nc; ++cell) {
src[cell] += well_contribs[phase].value()[cell];
}
}
// well rates contribs to polymer mass balance eqn.

11
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@ -79,13 +79,11 @@ namespace Opm {
/// \param[in] state reservoir state
/// \param[in] wstate well state
/// \param[in] polymer_inflow polymer influx
/// \param[in] src to caculate wc
void
step(const double dt,
PolymerBlackoilState& state ,
WellStateFullyImplicitBlackoil& wstate,
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
WellStateFullyImplicitBlackoil& wstate,
const std::vector<double>& polymer_inflow);
private:
typedef AutoDiffBlock<double> ADB;
@ -159,9 +157,8 @@ namespace Opm {
void
assemble(const double dt,
const PolymerBlackoilState& x,
const WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
const WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow);
V solveJacobianSystem() const;

15
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp
View File

@ -32,7 +32,6 @@
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
#include <opm/polymer/fullyimplicit/utilities.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
@ -232,11 +231,9 @@ namespace Opm
} else {
computePorevolume(grid_, props_.porosity(), porevol);
}
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
std::vector<double> initial_porevol = porevol;
std::vector<double> polymer_inflow_c(grid_.number_of_cells);
std::vector<double> transport_src(grid_.number_of_cells);
// Main simulation loop.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
@ -248,11 +245,11 @@ namespace Opm
//Main simulation loop.
while (!timer.done()) {
#if 0
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
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
if (wells_) {
@ -318,7 +315,7 @@ namespace Opm
// Run solver.
solver_timer.start();
FullyImplicitCompressiblePolymerSolver solver(grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells_manager.c_wells(), linsolver_);
solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src);
solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
// Stop timer and report.
solver_timer.stop();
const double st = solver_timer.secsSinceStart();
@ -330,12 +327,11 @@ namespace Opm
initial_porevol = porevol;
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
}
/*
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(),
@ -363,12 +359,12 @@ namespace Opm
std::cout << " Total prod reservoir volumes: "
<< std::setw(width) << tot_produced[0]
<< std::setw(width) << tot_produced[1] << std::endl;
*/
if (output_) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
outputWaterCut(watercut, output_dir_);
}
++timer;
prev_well_state = well_state;
@ -464,7 +460,7 @@ namespace Opm
}
}
#if 0
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
@ -476,6 +472,7 @@ namespace Opm
}
watercut.write(os);
}
#endif
}
} // namespace Opm

293
opm/polymer/fullyimplicit/utilities.cpp
View File

@ -1,293 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/linalg/blas_lapack.h>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
//#include <opm/autodiff/IncompPropsAdInterface.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/polymer/fullyimplicit/utilities.hpp>
#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <iostream>
#include <iterator>
namespace Opm
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (reservoir volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropsAdInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src)
{
int nc = props.numCells();
transport_src.clear();
transport_src.resize(nc, 0.0);
// Well contributions.
if (wells) {
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
if (np != 2) {
OPM_THROW(std::runtime_error, "computeTransportSource() requires a 2 phase case.");
}
std::vector<double> A(np*np);
for (int w = 0; w < nw; ++w) {
const double* comp_frac = wells->comp_frac + np*w;
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
const int perf_cell = wells->well_cells[perf];
double perf_rate = well_state.perfRates()[perf];
if (perf_rate > 0.0) {
// perf_rate is a total inflow reservoir rate, we want a surface water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Reservoir rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
assert(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
perf_rate *= comp_frac[0]; // Water reservoir volume rate.
}
}
transport_src[perf_cell] += perf_rate;
}
}
}
}
/// @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
// This function need a incompProps based on Ad.
/*
void computeInjectedProduced(const IncompPropsAdInterface& props,
const Opm::PolymerPropsAd& polymer_props,
const PolymerState& state,
const std::vector<double>& transport_src,
const std::vector<double>& 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<int> cells(num_cells);
const V p = Eigen::Map<const V>(&state.pressure()[0], num_cells, 1);
const DataBlock s = Eigen::Map<const DataBlock>(&state.saturation()[0], num_cells, np);
const V sw = s.col(0);
const V so = s.col(1);
const V c = Eigen::Map<const V>(&state.concentration()[0], num_cells, 1);
const V cmax = Eigen::Map<const V>(&state.maxconcentration()[0], num_cells, 1);
const V trans_src = Eigen::Map<const V>(&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<double> src_selector(src);
const V one = V::Constant(num_cells, 1.0);
const V zero = V::Zero(num_cells);
const std::vector<V> 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<V> mob(np);
mob[0] = krw_eff * inv_muw_eff;
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);
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.
/// 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 BlackoilPropsAdInterface& props,
const Opm::PolymerPropsAd& polymer_props,
const PolymerBlackoilState& state,
const std::vector<double>& transport_src,
const std::vector<double>& 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<int> cells(num_cells);
const V p = Eigen::Map<const V>(&state.pressure()[0], num_cells, 1);
const DataBlock s = Eigen::Map<const DataBlock>(&state.saturation()[0], num_cells, np);
const V sw = s.col(0);
const V so = s.col(1);
const V c = Eigen::Map<const V>(&state.concentration()[0], num_cells, 1);
const V cmax = Eigen::Map<const V>(&state.maxconcentration()[0], num_cells, 1);
const V trans_src = Eigen::Map<const V>(&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;
}
}
//Add PhasePresence make muOil() happy.
std::vector<PhasePresence> phaseCondition(num_cells);
for (int c = 0; c < num_cells; ++c) {
phaseCondition[c] = PhasePresence();
phaseCondition[c].setFreeWater();
phaseCondition[c].setFreeOil();
}
const Selector<double> src_selector(src);
const V one = V::Constant(num_cells, 1.0);
const V zero = V::Zero(num_cells);
const std::vector<V> kr = props.relperm(sw, so, zero, cells);
const V muw = props.muWat(p, cells);
const V muo = props.muOil(p, zero, phaseCondition, cells);
const V krw_eff = polymer_props.effectiveRelPerm(c, cmax, kr[0]);
const V inv_muw_eff = polymer_props.effectiveInvWaterVisc(c, muw.data());
std::vector<V> mob(np);
mob[0] = krw_eff * inv_muw_eff;
mob[1] = kr[1] / muo;
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];
}
}
}
} //namespace Opm

133
opm/polymer/fullyimplicit/utilities.hpp
View File

@ -1,133 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_UTILITIES_HEADER_INCLUDED
#define OPM_UTILITIES_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
//#include <opm/autodiff/IncompPropsAdInterface.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <iostream>
#include <iterator>
namespace Opm
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (reservoir volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropsAdInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& 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
// This function need a incompProps based on Ad.
/*
void computeInjectedProduced(const IncompPropsAdInterface& props,
const Opm::PolymerPropsAd& polymer_props,
const PolymerState& state,
const std::vector<double>& transport_src,
const std::vector<double>& 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.
/// 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 BlackoilPropsAdInterface& props,
const Opm::PolymerPropsAd& polymer_props,
const PolymerBlackoilState& state,
const std::vector<double>& transport_src,
const std::vector<double>& inj_c,
const double dt,
double* injected,
double* produced,
double& polyinj,
double& polyprod);
} //namespace Opm
#endif //OPM_UTILITIES_HEADER_INCLUDED