OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-02 12:36:54 -06:00
Code Issues Packages Projects Releases Wiki Activity

add well controls for polymer, but rate control just for

Browse Source 
water phase and oil phase.
This commit is contained in:
Liu Ming 2013-12-24 17:31:11 +08:00
parent 7a874427af
commit c37539b3ab
7 changed files with 900 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

242
examples/sim_fi2p_incomp_ad.cpp Normal file
View File

@ -0,0 +1,242 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
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/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp>
#include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp>
#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
#include <opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp>
#include <opm/core/utility/share_obj.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
std::cout << "\n================ Test program for fully implicit three-phase black-oil flow ===============\n\n";
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
if (!use_deck) {
OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
"Specify the deck with deck_filename=deckname.data (for example).");
}
boost::scoped_ptr<EclipseGridParser> deck;
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<IncompPropertiesInterface> props;
boost::scoped_ptr<IncompPropsAdInterface> new_props;
TwophaseState state;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
std::string deck_filename = param.get<std::string>("deck_filename");
deck.reset(new EclipseGridParser(deck_filename));
// Grid init
grid.reset(new GridManager(*deck));
// use the capitalized part of the deck's filename between the
// last '/' and the last '.' character as base name.
std::string baseName = deck_filename;
auto charPos = baseName.rfind('/');
if (charPos != std::string::npos)
baseName = baseName.substr(charPos + 1);
charPos = baseName.rfind('.');
if (charPos != std::string::npos)
baseName = baseName.substr(0, charPos);
baseName = boost::to_upper_copy(baseName);
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
new_props.reset(new IncompPropsAdFromDeck(*deck, *grid->c_grid()));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
// Gravity.
gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
}
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
// Linear solver.
LinearSolverFactory linsolver(param);
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::ofstream epoch_os;
std::string output_dir;
if (output) {
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
std::string filename = output_dir + "/epoch_timing.param";
epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
// open file to clean it. The file is appended to in SimulatorTwophase
filename = output_dir + "/step_timing.param";
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
step_os.close();
param.writeParam(output_dir + "/simulation.param");
}
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< " (number of epochs: "
<< (deck->numberOfEpochs()) << ")\n\n" << std::flush;
SimulatorReport rep;
// With a deck, we may have more epochs etc.
WellState well_state;
int step = 0;
SimulatorTimer simtimer;
// Use timer for last epoch to obtain total time.
deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
simtimer.init(*deck);
const double total_time = simtimer.totalTime();
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
// Set epoch index.
deck->setCurrentEpoch(epoch);
// Update the timer.
if (deck->hasField("TSTEP")) {
simtimer.init(*deck);
} else {
if (epoch != 0) {
OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
}
simtimer.init(param);
}
simtimer.setCurrentStepNum(step);
simtimer.setTotalTime(total_time);
// Report on start of epoch.
std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------"
<< "\n (number of steps: "
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.
if (epoch == 0) {
well_state.init(wells.c_wells(), state);
}
// Create and run simulator.
std::vector<double> src(grid->c_grid()->number_of_cells, 0.0);
src[0] = 10. / Opm::unit::day;
src[grid->c_grid()->number_of_cells-1] = -10. / Opm::unit::day;
SimulatorFullyImplicitTwophase simulator(param,
*grid->c_grid(),
*new_props,
wells,
linsolver,
src,
grav);
if (epoch == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, src, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}
// Update total timing report and remember step number.
rep += epoch_rep;
step = simtimer.currentStepNum();
}
std::cout << "\n\n================ End of simulation ===============\n\n";
rep.report(std::cout);
if (output) {
std::string filename = output_dir + "/walltime.param";
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
rep.reportParam(tot_os);
}
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

100
examples/sim_poly_fi2p_incomp_ad.cpp
View File

@ -98,8 +98,8 @@ try
boost::scoped_ptr<IncompPropsAdInterface> new_props;
// boost::scoped_ptr<PolymerPropsAd> polymer_props;
PolymerState state;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
std::string deck_filename = param.get<std::string>("deck_filename");
deck.reset(new EclipseGridParser(deck_filename));
@ -123,9 +123,8 @@ try
PolymerProperties polymer_props(*deck);
PolymerPropsAd polymer_props_ad(polymer_props);
// polymer_props.reset(new PolymerPropsAd(*deck, *grid->c_grid()));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Gravity.
gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
@ -137,7 +136,7 @@ try
}
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
const double* grav = use_gravity ? &gravity[0] : 0;
// Linear solver.
LinearSolverFactory linsolver(param);
@ -171,62 +170,85 @@ try
<< (deck->numberOfEpochs()) << ")\n\n" << std::flush;
SimulatorReport rep;
// With a deck, we may have more epochs etc.
// WellState well_state;
int step = 0;
SimulatorTimer simtimer;
// Use timer for last epoch to obtain total time.
deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
simtimer.init(*deck);
const double total_time = simtimer.totalTime();
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
// Set epoch index.
deck->setCurrentEpoch(epoch);
// Update the timer.
if (deck->hasField("TSTEP")) {
simtimer.init(*deck);
} else {
if (epoch != 0) {
OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
// With a deck, we may have more epochs etc.
WellState well_state;
int step = 0;
SimulatorTimer simtimer;
// Use timer for last epoch to obtain total time.
deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
simtimer.init(*deck);
const double total_time = simtimer.totalTime();
// Check for WPOLYMER presence in last epoch to decide
// polymer injection control type.
const bool use_wpolymer = deck->hasField("WPOLYMER");
if (use_wpolymer) {
if (param.has("poly_start_days")) {
OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
"You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
}
simtimer.init(param);
}
simtimer.setCurrentStepNum(step);
simtimer.setTotalTime(total_time);
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
// Set epoch index.
deck->setCurrentEpoch(epoch);
// Report on start of epoch.
std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------"
<< "\n (number of steps: "
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Update the timer.
if (deck->hasField("TSTEP")) {
simtimer.init(*deck);
} else {
if (epoch != 0) {
OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
}
simtimer.init(param);
}
simtimer.setCurrentStepNum(step);
simtimer.setTotalTime(total_time);
// Create new wells, well_state
// WellsManager wells(*deck, *grid->c_grid(), props->permeability());
// Report on start of epoch.
std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------"
<< "\n (number of steps: "
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Create new wells, polymer inflow controls.
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
boost::scoped_ptr<PolymerInflowInterface> polymer_inflow;
if (use_wpolymer) {
if (wells.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells()));
} else {
polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", polymer_props.cMax())));
}
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.
// if (epoch == 0) {
// well_state.init(wells.c_wells(), state);
// }
if (epoch == 0) {
well_state.init(wells.c_wells(), state);
}
// Create and run simulator.
#if 0
std::vector<double> src(grid->c_grid()->number_of_cells, 0.0);
src[0] = 10. / Opm::unit::day;
src[grid->c_grid()->number_of_cells-1] = -10. / Opm::unit::day;
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", polymer_props.cMax()));
#endif
SimulatorFullyImplicitTwophasePolymer simulator(param,
*grid->c_grid(),
*new_props,
polymer_props_ad,
linsolver,
polymer_inflow,
src);
wells,
*polymer_inflow,
grav);
if (epoch == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state);
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}

431
opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp
View File

@ -25,23 +25,6 @@
#include <algorithm>
namespace Opm {
namespace {
std::vector<int>
buildAllCells(const int nc)
{
std::vector<int> all_cells(nc);
for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
return all_cells;
}
struct Chop01 {
double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); }
};
}//anonymous namespace
@ -56,19 +39,74 @@ typedef Eigen::Array<double,
namespace {
std::vector<int>
buildAllCells(const int nc)
{
std::vector<int> all_cells(nc);
for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
return all_cells;
}
struct Chop01 {
double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); }
};
V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
{
const int nw = wells.number_of_wells;
const int nperf = wells.well_connpos[nw];
const int dim = grid.dimensions;
V wdp = V::Zero(nperf,1);
assert(wdp.size() == rho.size());
// Main loop, iterate over all perforations,
// using the following formula:
// wdp(perf) = g*(perf_z - well_ref_z)*rho(perf)
// where the total density rho(perf) is taken to be
// sum_p (rho_p*saturation_p) in the perforation cell.
// [although this is computed on the outside of this function].
for (int w = 0; w < nw; ++w) {
const double ref_depth = wells.depth_ref[w];
for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) {
const int cell = wells.well_cells[j];
const double cell_depth = grid.cell_centroids[dim * cell + dim - 1];
wdp[j] = rho[j]*grav*(cell_depth - ref_depth);
}
}
return wdp;
}
}//anonymous namespace
FullyImplicitTwophasePolymerSolver::
FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid,
const IncompPropsAdInterface& fluid,
const PolymerPropsAd& polymer_props_ad,
const LinearSolverInterface& linsolver)
const IncompPropsAdInterface& fluid,
const PolymerPropsAd& polymer_props_ad,
const LinearSolverInterface& linsolver,
const Wells& wells,
const double* gravity)
: grid_ (grid)
, fluid_(fluid)
, polymer_props_ad_ (polymer_props_ad)
, linsolver_(linsolver)
, wells_(wells)
, gravity_(gravity)
, cells_ (buildAllCells(grid.number_of_cells))
, ops_(grid)
, residual_(std::vector<ADB>(3, ADB::null()))
, wops_(wells)
, mob_(std::vector<ADB>(fluid.numPhases() + 1, ADB::null()))
, residual_( { std::vector<ADB>(fluid.numPhases() + 1, ADB::null()), ADB::null(), ADB::null()})
{
}
@ -76,11 +114,42 @@ typedef Eigen::Array<double,
FullyImplicitTwophasePolymerSolver::
WellOps::WellOps(const Wells& wells)
: w2p(wells.well_connpos[ wells.number_of_wells ],
wells.number_of_wells)
, p2w(wells.number_of_wells,
wells.well_connpos[ wells.number_of_wells ])
{
const int nw = wells.number_of_wells;
const int* const wpos = wells.well_connpos;
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> scatter, gather;
scatter.reserve(wpos[nw]);
gather .reserve(wpos[nw]);
for (int w = 0, i = 0; w < nw; ++w) {
for (; i < wpos[ w + 1 ]; ++i) {
scatter.push_back(Tri(i, w, 1.0));
gather .push_back(Tri(w, i, 1.0));
}
}
w2p.setFromTriplets(scatter.begin(), scatter.end());
p2w.setFromTriplets(gather .begin(), gather .end());
}
void
FullyImplicitTwophasePolymerSolver::
step(const double dt,
PolymerState& x,
const std::vector<double>& src,
PolymerState& x,
WellState& xw,
const std::vector<double>& polymer_inflow)
{
@ -94,11 +163,11 @@ typedef Eigen::Array<double,
const V pvdt = pvol / dt;
const SolutionState old_state = constantState(x);
const SolutionState old_state = constantState(x, xw);
const double atol = 1.0e-12;
const double rtol = 5.0e-8;
const int maxit = 40;
assemble(pvdt, old_state, x, src, polymer_inflow);
assemble(pvdt, old_state, x, xw, polymer_inflow);
const double r0 = residualNorm();
int it = 0;
@ -108,9 +177,9 @@ typedef Eigen::Array<double,
bool resTooLarge = r0 > atol;
while (resTooLarge && (it < maxit)) {
const V dx = solveJacobianSystem();
updateState(dx, x);
updateState(dx, x, xw);
assemble(pvdt, old_state, x, src, polymer_inflow);
assemble(pvdt, old_state, x, xw, polymer_inflow);
const double r = residualNorm();
@ -135,6 +204,8 @@ typedef Eigen::Array<double,
: pressure ( ADB::null())
, saturation (np, ADB::null())
, concentration ( ADB::null())
, qs ( ADB::null())
, bhp ( ADB::null())
{
}
@ -143,10 +214,23 @@ typedef Eigen::Array<double,
FullyImplicitTwophasePolymerSolver::SolutionState
FullyImplicitTwophasePolymerSolver::constantState(const PolymerState& x)
FullyImplicitTwophasePolymerSolver::constantState(const PolymerState& x,
const WellState& xw)
{
const int nc = grid_.number_of_cells;
const int np = x.numPhases();
// The block pattern assumes the following primary variables:
// pressure
// water saturation
// polymer concentration
// well surface rates
// well bottom-hole pressure
// Note that oil is assumed to always be present, but is never
// a primary variable.
std::vector<int> bpat(np + 1, nc);
bpat.push_back(xw.bhp().size() * np);
bpat.push_back(xw.bhp().size());
SolutionState state(np);
@ -167,6 +251,20 @@ typedef Eigen::Array<double,
const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
state.concentration = ADB::constant(c);
// Well rates.
assert (not xw.wellRates().empty());
// Need to reshuffle well rates, from ordered by wells, then phase,
// to ordered by phase, then wells.
const int nw = wells_.number_of_wells;
// The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
const V qs = Eigen::Map<const V>(wrates.data(), nw * np);
state.qs = ADB::constant(qs, bpat);
// Bottom hole pressure.
assert (not xw.bhp().empty());
const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
state.bhp = ADB::constant(bhp, bpat);
return state;
}
@ -175,13 +273,14 @@ typedef Eigen::Array<double,
FullyImplicitTwophasePolymerSolver::SolutionState
FullyImplicitTwophasePolymerSolver::variableState(const PolymerState& x)
FullyImplicitTwophasePolymerSolver::variableState(const PolymerState& x,
const WellState& xw)
{
const int nc = grid_.number_of_cells;
const int np = x.numPhases();
std::vector<V> vars0;
vars0.reserve(np);
vars0.reserve(np + 3);
// Initial pressure.
assert (not x.pressure().empty());
@ -199,6 +298,21 @@ typedef Eigen::Array<double,
const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
vars0.push_back(c);
// Initial well rates.
assert (not xw.wellRates().empty());
// Need to reshuffle well rates, from ordered by wells, then phase,
// to ordered by phase, then wells.
const int nw = wells_.number_of_wells;
// The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
const V qs = Eigen::Map<const V>(wrates.data(), nw*np);
vars0.push_back(qs);
// Initial well bottom hole pressure.
assert (not xw.bhp().size());
const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
vars0.push_back(bhp);
std::vector<ADB> vars = ADB::variables(vars0);
SolutionState state(np);
@ -220,6 +334,10 @@ typedef Eigen::Array<double,
// Concentration.
state.concentration = vars[nextvar++];
// Qs.
state.qs = vars[ nextvar++ ];
// BHP.
state.bhp = vars[ nextvar++ ];
assert(nextvar == int(vars.size()));
return state;
@ -237,7 +355,7 @@ typedef Eigen::Array<double,
cmax(i) = (cmax(i) > c.value()(i)) ? cmax(i) : c.value()(i);
}
return ADB::constant(cmax);
return ADB::constant(cmax, c.blockPattern());
}
@ -247,11 +365,11 @@ typedef Eigen::Array<double,
assemble(const V& pvdt,
const SolutionState& old_state,
const PolymerState& x,
const std::vector<double>& src,
const WellState& xw,
const std::vector<double>& polymer_inflow)
{
// Create the primary variables.
const SolutionState state = variableState(x);
const SolutionState state = variableState(x, xw);
// -------- Mass balance equations for water and oil --------
const V trans = subset(transmissibility(), ops_.internal_faces);
@ -260,27 +378,147 @@ typedef Eigen::Array<double,
const ADB cmax = computeCmax(state.concentration);
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]);
std::cout << "krw krw_eff\n";
for (int i = 0; i < grid_.number_of_cells; ++i) {
std::cout <<kr[0].value()(i)<<" "<<krw_eff.value()(i)<< std::endl;
}
const ADB mc = computeMc(state);
// const std::vector<ADB> mflux = computeMassFlux(trans, mc, kr[0], krw_eff, state);
const std::vector<ADB> mflux = computeMassFlux(trans, mc, kr[1], krw_eff, state);
const std::vector<ADB> source = accumSource(kr[1], krw_eff, state.concentration, src, polymer_inflow);
//const std::vector<ADB> source = accumSource(kr[1], krw_eff, state.concentration, src, polymer_inflow);
// const std::vector<ADB> 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_[0] = pvdt * (state.saturation[0] - old_state.saturation[0])
+ ops_.div * mflux[0] - source[0];
residual_[1] = pvdt * (state.saturation[1] - old_state.saturation[1])
+ ops_.div * mflux[1] - source[1];
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_[2] = pvdt * (state.saturation[0] * state.concentration
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] - source[2];
+ ops_.div * mflux[2];
// -------- Well equation, and well contributions to the mass balance equations --------
// Contribution to mass balance will have to wait.
const int nc = grid_.number_of_cells;
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
const V transw = Eigen::Map<const V>(wells_.WI, nperf);
const ADB& bhp = state.bhp;
const DataBlock well_s = wops_.w2p * Eigen::Map<const DataBlock>(wells_.comp_frac, nw, np).matrix();
// Extract variables for perforation cell pressures
// and corresponding perforation well pressures.
const ADB p_perfcell = subset(state.pressure, well_cells);
// Finally construct well perforation pressures and well flows.
// Compute well pressure differentials.
// Construct pressure difference vector for wells.
const int dim = grid_.dimensions;
if (gravity_) {
for (int dd = 0; dd < dim -1; ++dd) {
assert(g[dd] == 0.0);
}
}
ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern());
for (int phase = 0; phase < 2; ++phase) {
// For incompressible flow cell rho is the same.
const ADB cell_rho = fluidDensity(phase, state.pressure);
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, state.pressure);
const V fraction = compi.col(phase);
inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells);
}
const V rho_perf_cell = subset(cell_rho_total, well_cells).value();
const V rho_perf_well = inj_rho_total.value();
V prodperfs = V::Constant(nperf, -1.0);
for (int w = 0; w < nw; ++w) {
if (wells_.type[w] == PRODUCER) {
std::fill(prodperfs.data() + wells_.well_connpos[w],
prodperfs.data() + wells_.well_connpos[w+1], 1.0);
}
}
const Selector<double> producer(prodperfs);
const V rho_perf = producer.select(rho_perf_cell, rho_perf_well);
const V well_perf_dp = computePerfPress(grid_, wells_, rho_perf, gravity_ ? gravity_[dim - 1] : 0.0);
const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp;
const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
// DUMP(nkgradp_well);
const Selector<double> cell_to_well_selector(nkgradp_well.value());
ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern());
ADB perf_total_mob = subset(mob_[0], well_cells) + subset(mob_[1], well_cells);
std::vector<ADB> well_contribs(np, ADB::null());
std::vector<ADB> well_perf_rates(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
// const ADB& cell_b = rq_[phase].b;
// const ADB perf_b = subset(cell_b, well_cells);
const ADB& cell_mob = mob_[phase];
const V well_fraction = compi.col(phase);
// Using total mobilities for all phases for injection.
const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob;
const ADB perf_mob = producer.select(subset(cell_mob, well_cells),
perf_mob_injector);
const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations.
well_perf_rates[phase] = perf_flux;
const ADB well_rates = wops_.p2w * well_perf_rates[phase];
well_rates_all += superset(well_rates, Span(nw, 1, phase*nw), nw*np);
// const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
well_contribs[phase] = superset(perf_flux, well_cells, nc);
// DUMP(well_contribs[phase]);
residual_.mass_balance[phase] += well_contribs[phase];
}
// well rates contribs to polymer mass balance eqn.
// for injection wells.
const V polyin = Eigen::Map<const V>(& 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);
residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_c, well_cells, nc);
// Set the well flux equation
residual_.well_flux_eq = state.qs + well_rates_all;
// DUMP(residual_.well_flux_eq);
// Handling BHP and SURFACE_RATE wells.
V bhp_targets(nw);
V rate_targets(nw);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
if (wc->type[wc->current] == BHP) {
bhp_targets[w] = wc->target[wc->current];
rate_targets[w] = -1e100;
} else if (wc->type[wc->current] == SURFACE_RATE) {
bhp_targets[w] = -1e100;
rate_targets[w] = wc->target[wc->current];
for (int phase = 0; phase < np; ++phase) {
rate_distr.insert(w, phase*nw + w) = wc->distr[phase];
}
} else {
OPM_THROW(std::runtime_error, "Can only handle BHP type controls.");
}
}
const ADB bhp_residual = bhp - bhp_targets;
const ADB rate_residual = rate_distr * state.qs - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// residual_.well_eq = bhp_residual;
}
@ -289,24 +527,39 @@ typedef Eigen::Array<double,
const ADB& mc,
const ADB& kro,
const ADB& krw_eff,
const SolutionState& state ) const
const SolutionState& state )
{
const double* mus = fluid_.viscosity();
std::vector<ADB> mflux;
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus);
ADB wat_mob = krw_eff * inv_wat_eff_vis;
ADB oil_mob = kro / V::Constant(kro.size(), 1, mus[1]);
ADB poly_mob = mc * krw_eff * inv_wat_eff_vis;
mob_[0] = krw_eff * inv_wat_eff_vis;
mob_[1] = kro / V::Constant(kro.size(), 1, mus[1]);
mob_[2] = mc * krw_eff * inv_wat_eff_vis;
const ADB dp = ops_.ngrad * state.pressure;
const int nc = grid_.number_of_cells;
V z(nc);
// Compute z coordinates
for (int c = 0; c < nc; ++c){
z[c] = grid_.cell_centroids[c * 3 + 2];
}
for (int phase = 0; phase < 2; ++phase) {
const ADB rho = fluidDensity(phase, state.pressure);
const ADB rhoavg = ops_.caver * rho;
const ADB dp = ops_.ngrad * state.pressure
- gravity_[2] * (rhoavg * (ops_.ngrad * z.matrix()));
const ADB head = trans * dp;
UpwindSelector<double> upwind(grid_, ops_, head.value());
mflux.push_back(upwind.select(mob_[phase])*head);
}
// polymer mass flux.
const ADB rho = fluidDensity(0, state.pressure);
const ADB rhoavg = ops_.caver * rho;
const ADB dp = ops_.ngrad * state.pressure
- gravity_[2] * (rhoavg * (ops_.ngrad * z.matrix()));
const ADB head = trans * dp;
UpwindSelector<double> upwind(grid_, ops_, head.value());
mflux.push_back(upwind.select(wat_mob)*head);
mflux.push_back(upwind.select(oil_mob)*head);
mflux.push_back(upwind.select(poly_mob)*head);
mflux.push_back(upwind.select(mob_[2])*head);
return mflux;
}
@ -363,16 +616,12 @@ typedef Eigen::Array<double,
const ADB& krw_eff,
const ADB& c) const
{
const double* mus = fluid_.viscosity();
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(c, mus);
ADB wat_mob = krw_eff * inv_wat_eff_vis;
ADB oil_mob = kro / V::Constant(kro.size(), 1, mus[1]);
ADB total_mob = wat_mob + oil_mob;
ADB total_mob = mob_[0] + mob_[1];
std::vector<ADB> fracflow;
fracflow.push_back(wat_mob / total_mob);
fracflow.push_back(oil_mob / total_mob);
fracflow.push_back(mob_[0] / total_mob);
fracflow.push_back(mob_[1] / total_mob);
return fracflow;
}
@ -388,15 +637,17 @@ typedef Eigen::Array<double,
if (np != 2) {
OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwophasePolymerSolver.");
}
ADB mass_phase_res = vertcat(residual_[0], residual_[1]);
ADB mass_res = collapseJacs(vertcat(mass_phase_res, residual_[2]));
ADB mass_res = vertcat(residual_.mass_balance[0], residual_.mass_balance[1]);
mass_res = vertcat(mass_res, residual_.mass_balance[2]);
ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
ADB total_res = collapseJacs(vertcat(mass_res, well_res));
const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = mass_res.derivative()[0];
V dx(V::Zero(mass_res.size()));
const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_res.derivative()[0];
V dx(V::Zero(total_res.size()));
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
mass_res.value().data(), dx.data());
total_res.value().data(), dx.data());
if (!rep.converged) {
OPM_THROW(std::runtime_error,
"FullyImplicitBlackoilSolver::solveJacobianSystem(): "
@ -410,13 +661,12 @@ typedef Eigen::Array<double,
void FullyImplicitTwophasePolymerSolver::updateState(const V& dx,
PolymerState& state) const
PolymerState& state,
WellState& well_state) const
{
const int np = fluid_.numPhases();
const int nc = grid_.number_of_cells;
const V null;
assert(null.size() == 0);
const V zero = V::Zero(nc);
const int nw = wells_.number_of_wells;
const V one = V::Constant(nc, 1.0);
// Extract parts of dx corresponding to each part.
@ -426,6 +676,10 @@ typedef Eigen::Array<double,
varstart += dsw.size();
const V dc = subset(dx, Span(nc, 1, varstart));
varstart += dc.size();
const V dqs = subset(dx, Span(np*nw, 1, varstart));
varstart += dqs.size();
const V dbhp = subset(dx, Span(nw, 1, varstart));
varstart += dbhp.size();
assert(varstart == dx.size());
@ -454,6 +708,24 @@ typedef Eigen::Array<double,
const V c = c_old - dc;
std::copy(&c[0], &c[0] + nc, state.concentration().begin());
// Qs update.
// Since we need to update the wellrates, that are ordered by wells,
// from dqs which are ordered by phase, the simplest is to compute
// dwr, which is the data from dqs but ordered by wells.
const DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
const V dwr = Eigen::Map<const V>(wwr.data(), nw*np);
const V wr_old = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np);
const V wr = wr_old - dwr;
std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
// Bhp update.
const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
const V bhp = bhp_old - dbhp;
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
}
@ -486,18 +758,31 @@ typedef Eigen::Array<double,
{
double r = 0;
for (std::vector<ADB>::const_iterator
b = residual_.begin(),
e = residual_.end();
b = residual_.mass_balance.begin(),
e = residual_.mass_balance.end();
b != e; ++b)
{
r = std::max(r, (*b).value().matrix().norm());
}
r = std::max(r, residual_.well_flux_eq.value().matrix().norm());
r = std::max(r, residual_.well_eq.value().matrix().norm());
return r;
}
ADB
FullyImplicitTwophasePolymerSolver::fluidDensity(const int phase,
const ADB p) const
{
const double* rhos = fluid_.surfaceDensity();
ADB rho = ADB::constant(V::Constant(grid_.number_of_cells, 1, rhos[phase]),
p.blockPattern());
return rho;
}
V

47
opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp
View File

@ -12,10 +12,11 @@
struct UnstructuredGrid;
struct Wells;
namespace Opm {
class LinearSolverInterface;
class PolymerState;
class WellState;
class FullyImplicitTwophasePolymerSolver
{
@ -23,13 +24,14 @@ namespace Opm {
FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid,
const IncompPropsAdInterface& fluid,
const PolymerPropsAd& polymer_props_ad,
const LinearSolverInterface& linsolver);
const LinearSolverInterface& linsolver,
const Wells& wells,
const double* gravity);
void step(const double dt,
PolymerState& state,
const std::vector<double>& src,
const std::vector<double>& polymer_inflow
);
WellState& well_state,
const std::vector<double>& polymer_inflow);
private:
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
@ -43,29 +45,47 @@ namespace Opm {
ADB pressure;
std::vector<ADB> saturation;
ADB concentration;
ADB qs;
ADB bhp;
};
struct WellOps {
WellOps(const Wells& wells);
M w2p; // well -> perf (scatter)
M p2w; // perf -> well (gather)
};
const UnstructuredGrid& grid_;
const IncompPropsAdInterface& fluid_;
const PolymerPropsAd& polymer_props_ad_;
const LinearSolverInterface& linsolver_;
const Wells& wells_;
const double* gravity_;
const std::vector<int> cells_;
HelperOps ops_;
std::vector<ADB> residual_;
const WellOps wops_;
std::vector<ADB> mob_;
struct {
std::vector<ADB> mass_balance;
ADB well_eq;
ADB well_flux_eq;
} residual_;
SolutionState
constantState(const PolymerState& x);
constantState(const PolymerState& x,
const WellState& xw);
SolutionState
variableState(const PolymerState& x);
variableState(const PolymerState& x,
const WellState& xw);
void
assemble(const V& pvdt,
const SolutionState& old_state,
const PolymerState& x,
const std::vector<double>& src,
const PolymerState& x,
const WellState& xw,
const std::vector<double>& polymer_inflow);
V solveJacobianSystem() const;
void updateState(const V& dx,
PolymerState& x) const;
PolymerState& x,
WellState& xw) const;
std::vector<ADB>
computeRelPerm(const SolutionState& state) const;
V
@ -76,7 +96,7 @@ namespace Opm {
const ADB& mc,
const ADB& kro,
const ADB& krw_eff,
const SolutionState& state ) const;
const SolutionState& state );
std::vector<ADB>
accumSource(const ADB& kro,
@ -109,6 +129,9 @@ namespace Opm {
const double
fluidDensity(const int phase) const;
ADB
fluidDensity(const int phase,
const ADB p) const;
ADB
transMult(const ADB& p) const;
};

104
opm/polymer/fullyimplicit/PolymerWellState.hpp Normal file
View File

@ -0,0 +1,104 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
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_WELLSTATE_HEADER_INCLUDED
#define OPM_WELLSTATE_HEADER_INCLUDED
#include <opm/core/wells.h>
#include <vector>
namespace Opm
{
/// The state of a set of wells.
class WellState
{
public:
/// Allocate and initialize if wells is non-null.
/// Also tries to give useful initial values to the bhp() and
/// wellRates() fields, depending on controls. The
/// perfRates() field is filled with zero, and perfPress()
/// with -1e100.
template <class State>
void init(const Wells* wells, const State& state)
{
if (wells) {
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases + 1;
bhp_.resize(nw);
wellrates_.resize(nw * np, 0.0);
for (int w = 0; w < nw; ++w) {
const WellControls* ctrl = wells->ctrls[w];
// Initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a little
// above or below (depending on if the well is an
// injector or producer) pressure in first perforation
// cell.
if ((ctrl->current < 0) || // SHUT
(ctrl->type[ctrl->current] != BHP)) {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*state.pressure()[first_cell];
} else {
bhp_[w] = ctrl->target[ctrl->current];
}
// Initialize well rates to match controls if type is SURFACE_RATE
if ((ctrl->current >= 0) && // open well
(ctrl->type[ctrl->current] == SURFACE_RATE)) {
const double rate_target = ctrl->target[ctrl->current];
for (int p = 0; p < np; ++p) {
const double phase_distr = ctrl->distr[np * ctrl->current + p];
wellrates_[np*w + p] = rate_target * phase_distr;
}
}
}
// The perforation rates and perforation pressures are
// not expected to be consistent with bhp_ and wellrates_
// after init().
perfrates_.resize(wells->well_connpos[nw], 0.0);
perfpress_.resize(wells->well_connpos[nw], -1e100);
}
}
/// One bhp pressure per well.
std::vector<double>& bhp() { return bhp_; }
const std::vector<double>& bhp() const { return bhp_; }
/// One rate per well and phase.
std::vector<double>& wellRates() { return wellrates_; }
const std::vector<double>& wellRates() const { return wellrates_; }
/// One rate per well connection.
std::vector<double>& perfRates() { return perfrates_; }
const std::vector<double>& perfRates() const { return perfrates_; }
/// One pressure per well connection.
std::vector<double>& perfPress() { return perfpress_; }
const std::vector<double>& perfPress() const { return perfpress_; }
private:
std::vector<double> bhp_;
std::vector<double> wellrates_;
std::vector<double> perfrates_;
std::vector<double> perfpress_;
};
} // namespace Opm
#endif // OPM_WELLSTATE_HEADER_INCLUDED

114
opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp
View File

@ -26,6 +26,7 @@
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
@ -38,6 +39,7 @@
#include <opm/core/grid/ColumnExtract.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
#include <boost/filesystem.hpp>
@ -61,11 +63,13 @@ namespace Opm
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
const PolymerInflowInterface& polymer_inflow,
std::vector<double>& src);
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
const double* gravity);
SimulatorReport run(SimulatorTimer& timer,
PolymerState& state);
PolymerState& state,
WellState& well_state);
private:
@ -75,12 +79,15 @@ namespace Opm
std::string output_dir_;
int output_interval_;
// Parameters for well control
bool check_well_controls_;
int max_well_control_iterations_;
// Observed objects.
const UnstructuredGrid& grid_;
const IncompPropsAdInterface& props_;
const PolymerPropsAd& polymer_props_;
const PolymerInflowInterface& polymer_inflow_;
const std::vector<double>& src_;
WellsManager& wells_manager_;
const Wells* wells_;
PolymerInflowInterface& polymer_inflow_;
// Solvers
FullyImplicitTwophasePolymerSolver solver_;
// Misc. data
@ -90,15 +97,17 @@ namespace Opm
SimulatorFullyImplicitTwophasePolymer::SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
const PolymerInflowInterface& polymer_inflow,
std::vector<double>& src)
SimulatorFullyImplicitTwophasePolymer::
SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, props, polymer_props, linsolver, polymer_inflow, src));
pimpl_.reset(new Impl(param, grid, props, polymer_props, linsolver, wells_manager, polymer_inflow, gravity));
}
@ -106,9 +115,10 @@ namespace Opm
SimulatorReport SimulatorFullyImplicitTwophasePolymer::run(SimulatorTimer& timer,
PolymerState& state)
PolymerState& state,
WellState& well_state)
{
return pimpl_->run(timer, state);
return pimpl_->run(timer, state, well_state);
}
@ -204,6 +214,37 @@ namespace Opm
}
*/
static void outputWellStateMatlab(WellState& well_state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["bhp"] = &well_state.bhp();
dm["wellrates"] = &well_state.wellRates();
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath);
}
fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
std::ofstream file(fname.str().c_str());
if (!file) {
OPM_THROW(std::runtime_error,"Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
SimulatorFullyImplicitTwophasePolymer::Impl::Impl(const parameter::ParameterGroup& param,
@ -211,14 +252,16 @@ namespace Opm
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
const PolymerInflowInterface& polymer_inflow,
std::vector<double>& src)
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
const double* gravity)
: grid_(grid),
props_(props),
polymer_props_(polymer_props),
wells_manager_(wells_manager),
wells_(wells_manager.c_wells()),
polymer_inflow_(polymer_inflow),
src_ (src),
solver_(grid_, props_, polymer_props_, linsolver)
solver_(grid_, props_, polymer_props_, linsolver, *wells_manager.c_wells(), gravity)
{
// For output.
@ -245,7 +288,8 @@ namespace Opm
}
SimulatorReport SimulatorFullyImplicitTwophasePolymer::Impl::run(SimulatorTimer& timer,
PolymerState& state)
PolymerState& state,
WellState& well_state)
{
// Initialisation.
@ -287,18 +331,22 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
}
SimulatorReport sreport;
bool well_control_passed = !check_well_controls_;
int well_control_iteration = 0;
do {
// Run solver.
// Find inflow rate.
const double current_time = timer.currentTime();
double stepsize = timer.currentStepLength();
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
solver_timer.start();
solver_.step(timer.currentStepLength(), state, src_, polymer_inflow_c);
std::vector<double> initial_pressure = state.pressure();
solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
// Stop timer and report.
solver_timer.stop();
@ -308,6 +356,27 @@ namespace Opm
stime += st;
sreport.pressure_time = st;
// Optionally, check if well controls are satisfied.
if (check_well_controls_) {
Opm::computePhaseFlowRatesPerWell(*wells_,
well_state.perfRates(),
fractional_flows,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Update pore volumes if rock is compressible.
@ -369,6 +438,7 @@ 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 (wells_) {

12
opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp
View File

@ -24,7 +24,7 @@
#include <vector>
struct UnstructuredGrid;
struct Wells;
namespace Opm
{
namespace parameter { class ParameterGroup; }
@ -34,6 +34,8 @@ namespace Opm
class PolymerState;
class PolymerPropsAd;
class PolymerInflowInterface;
class WellsManager;
class WellState;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation.
@ -64,8 +66,9 @@ namespace Opm
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
const PolymerInflowInterface& polymer_inflow,
std::vector<double>& src);
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
const double* gravity);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will
@ -75,7 +78,8 @@ namespace Opm
/// \param[in,out] well_state state of wells: bhp, perforation rates
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
PolymerState& state);
PolymerState& state,
WellState& well_state);
private:
class Impl;