mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
new simulator for fully implicit compressible twophase
polymer works.
This commit is contained in:
Liu Ming
parent
c37539b3ab
commit
cf28164a5a
272
examples/sim_poly_fi2p_comp_ad.cpp
Normal file
272
examples/sim_poly_fi2p_comp_ad.cpp
Normal file
@ -0,0 +1,272 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
|
||||
#if HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif // HAVE_CONFIG_H
|
||||
|
||||
#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/BlackoilPropertiesBasic.hpp>
|
||||
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
|
||||
#include <opm/core/props/rock/RockCompressibility.hpp>
|
||||
|
||||
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
||||
|
||||
#include <opm/polymer/PolymerBlackoilState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
|
||||
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp>
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp>
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
|
||||
#include <opm/polymer/PolymerProperties.hpp>
|
||||
#include <opm/polymer/PolymerInflow.hpp>
|
||||
#include <opm/polymer/PolymerState.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<BlackoilPropertiesInterface> props;
|
||||
boost::scoped_ptr<BlackoilPropsAdInterface> new_props;
|
||||
boost::scoped_ptr<RockCompressibility> rock_comp;
|
||||
PolymerBlackoilState 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);
|
||||
|
||||
Opm::EclipseWriter outputWriter(param, share_obj(*deck), share_obj(*grid->c_grid()));
|
||||
// Rock and fluid init
|
||||
props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param));
|
||||
new_props.reset(new BlackoilPropsAdFromDeck(*deck, *grid->c_grid()));
|
||||
PolymerProperties polymer_props(*deck);
|
||||
PolymerPropsAd polymer_props_ad(polymer_props);
|
||||
// check_well_controls = param.getDefault("check_well_controls", false);
|
||||
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
|
||||
// Rock compressibility.
|
||||
rock_comp.reset(new RockCompressibility(*deck));
|
||||
// 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);
|
||||
initBlackoilSurfvol(*grid->c_grid(), *props, 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();
|
||||
// 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.");
|
||||
}
|
||||
}
|
||||
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, 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);
|
||||
}
|
||||
|
||||
// Create and run simulator.
|
||||
SimulatorFullyImplicitCompressiblePolymer simulator(param,
|
||||
*grid->c_grid(),
|
||||
*new_props,
|
||||
polymer_props_ad,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
*polymer_inflow,
|
||||
linsolver,
|
||||
grav);
|
||||
if (epoch == 0) {
|
||||
warnIfUnusedParams(param);
|
||||
}
|
||||
SimulatorReport epoch_rep = simulator.run(simtimer, state, 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;
|
||||
}
|
||||
|
||||
@ -39,6 +39,10 @@ namespace Opm
|
||||
concentration_.resize(g.number_of_cells, 0.0);
|
||||
cmax_.resize(g.number_of_cells, 0.0);
|
||||
}
|
||||
int numPhases() const
|
||||
{
|
||||
return state_blackoil_.numPhases();
|
||||
}
|
||||
|
||||
enum ExtremalSat { MinSat = BlackoilState::MinSat, MaxSat = BlackoilState::MaxSat };
|
||||
|
||||
|
||||
588
opm/polymer/fullyimplicit/BlackoilPropsAd.cpp
Normal file
588
opm/polymer/fullyimplicit/BlackoilPropsAd.cpp
Normal file
@ -0,0 +1,588 @@
|
||||
/*
|
||||
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 <config.h>
|
||||
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAd.hpp>
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
|
||||
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
// Making these typedef to make the code more readable.
|
||||
typedef BlackoilPropsAd::ADB ADB;
|
||||
typedef BlackoilPropsAd::V V;
|
||||
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
|
||||
|
||||
/// Constructor wrapping an opm-core black oil interface.
|
||||
BlackoilPropsAd::BlackoilPropsAd(const BlackoilPropertiesInterface& props)
|
||||
: props_(props),
|
||||
pu_(props.phaseUsage())
|
||||
{
|
||||
}
|
||||
|
||||
////////////////////////////
|
||||
// Rock interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return D, the number of spatial dimensions.
|
||||
int BlackoilPropsAd::numDimensions() const
|
||||
{
|
||||
return props_.numDimensions();
|
||||
}
|
||||
|
||||
/// \return N, the number of cells.
|
||||
int BlackoilPropsAd::numCells() const
|
||||
{
|
||||
return props_.numCells();
|
||||
}
|
||||
|
||||
/// \return Array of N porosity values.
|
||||
const double* BlackoilPropsAd::porosity() const
|
||||
{
|
||||
return props_.porosity();
|
||||
}
|
||||
|
||||
/// \return Array of ND^2 permeability values.
|
||||
/// The D^2 permeability values for a cell are organized as a matrix,
|
||||
/// which is symmetric (so ordering does not matter).
|
||||
const double* BlackoilPropsAd::permeability() const
|
||||
{
|
||||
return props_.permeability();
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////
|
||||
// Fluid interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return Number of active phases (also the number of components).
|
||||
int BlackoilPropsAd::numPhases() const
|
||||
{
|
||||
return props_.numPhases();
|
||||
}
|
||||
|
||||
/// \return Object describing the active phases.
|
||||
PhaseUsage BlackoilPropsAd::phaseUsage() const
|
||||
{
|
||||
return props_.phaseUsage();
|
||||
}
|
||||
|
||||
// ------ Density ------
|
||||
|
||||
/// Densities of stock components at surface conditions.
|
||||
/// \return Array of 3 density values.
|
||||
const double* BlackoilPropsAd::surfaceDensity() const
|
||||
{
|
||||
return props_.surfaceDensity();
|
||||
}
|
||||
|
||||
|
||||
// ------ Viscosity ------
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V BlackoilPropsAd::muWat(const V& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block mu(n, np);
|
||||
props_.viscosity(n, pw.data(), z.data(), cells.data(), mu.data(), 0);
|
||||
return mu.col(pu_.phase_pos[Water]);
|
||||
}
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V BlackoilPropsAd::muOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
if (pu_.phase_used[Gas]) {
|
||||
// Faking a z with the right ratio:
|
||||
// rs = zg/zo
|
||||
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
|
||||
z.col(pu_.phase_pos[Gas]) = rs;
|
||||
}
|
||||
Block mu(n, np);
|
||||
props_.viscosity(n, po.data(), z.data(), cells.data(), mu.data(), 0);
|
||||
return mu.col(pu_.phase_pos[Oil]);
|
||||
}
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V BlackoilPropsAd::muGas(const V& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block mu(n, np);
|
||||
props_.viscosity(n, pg.data(), z.data(), cells.data(), mu.data(), 0);
|
||||
return mu.col(pu_.phase_pos[Gas]);
|
||||
}
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB BlackoilPropsAd::muWat(const ADB& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
#if 1
|
||||
return ADB::constant(muWat(pw.value(), cells), pw.blockPattern());
|
||||
#else
|
||||
if (!pu_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.value().size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block mu(n, np);
|
||||
Block dmu(n, np);
|
||||
props_.viscosity(n, pw.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
|
||||
ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Water]));
|
||||
const int num_blocks = pw.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dmu_diag * pw.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu.col(pu_.phase_pos[Water]), jacs);
|
||||
#endif
|
||||
}
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB BlackoilPropsAd::muOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
#if 1
|
||||
return ADB::constant(muOil(po.value(), rs.value(), cells), po.blockPattern());
|
||||
#else
|
||||
if (!pu_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.value().size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
if (pu_.phase_used[Gas]) {
|
||||
// Faking a z with the right ratio:
|
||||
// rs = zg/zo
|
||||
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
|
||||
z.col(pu_.phase_pos[Gas]) = rs.value();
|
||||
}
|
||||
Block mu(n, np);
|
||||
Block dmu(n, np);
|
||||
props_.viscosity(n, po.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
|
||||
ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Oil]));
|
||||
const int num_blocks = po.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
// For now, we deliberately ignore the derivative with respect to rs,
|
||||
// since the BlackoilPropertiesInterface class does not evaluate it.
|
||||
// We would add to the next line: + dmu_drs_diag * rs.derivative()[block]
|
||||
jacs[block] = dmu_diag * po.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu.col(pu_.phase_pos[Oil]), jacs);
|
||||
#endif
|
||||
}
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB BlackoilPropsAd::muGas(const ADB& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
#if 1
|
||||
return ADB::constant(muGas(pg.value(), cells), pg.blockPattern());
|
||||
#else
|
||||
if (!pu_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.value().size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block mu(n, np);
|
||||
Block dmu(n, np);
|
||||
props_.viscosity(n, pg.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
|
||||
ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Gas]));
|
||||
const int num_blocks = pg.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dmu_diag * pg.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
// ------ Formation volume factor (b) ------
|
||||
|
||||
// These methods all call the matrix() method, after which the variable
|
||||
// (also) called 'matrix' contains, in each row, the A = RB^{-1} matrix for
|
||||
// a cell. For three-phase black oil:
|
||||
// A = [ bw 0 0
|
||||
// 0 bo 0
|
||||
// 0 b0*rs bw ]
|
||||
// Where b = B^{-1}.
|
||||
// Therefore, we extract the correct diagonal element, and are done.
|
||||
// When we need the derivatives (w.r.t. p, since we don't do w.r.t. rs),
|
||||
// we also get the following derivative matrix:
|
||||
// A = [ dbw 0 0
|
||||
// 0 dbo 0
|
||||
// 0 db0*rs dbw ]
|
||||
// Again, we just extract a diagonal element.
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V BlackoilPropsAd::bWat(const V& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call bWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block matrix(n, np*np);
|
||||
props_.matrix(n, pw.data(), z.data(), cells.data(), matrix.data(), 0);
|
||||
const int wi = pu_.phase_pos[Water];
|
||||
return matrix.col(wi*np + wi);
|
||||
}
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V BlackoilPropsAd::bOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call bOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
if (pu_.phase_used[Gas]) {
|
||||
// Faking a z with the right ratio:
|
||||
// rs = zg/zo
|
||||
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
|
||||
z.col(pu_.phase_pos[Gas]) = rs;
|
||||
}
|
||||
Block matrix(n, np*np);
|
||||
props_.matrix(n, po.data(), z.data(), cells.data(), matrix.data(), 0);
|
||||
const int oi = pu_.phase_pos[Oil];
|
||||
return matrix.col(oi*np + oi);
|
||||
}
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V BlackoilPropsAd::bGas(const V& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call bGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block matrix(n, np*np);
|
||||
props_.matrix(n, pg.data(), z.data(), cells.data(), matrix.data(), 0);
|
||||
const int gi = pu_.phase_pos[Gas];
|
||||
return matrix.col(gi*np + gi);
|
||||
}
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB BlackoilPropsAd::bWat(const ADB& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.value().size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block matrix(n, np*np);
|
||||
Block dmatrix(n, np*np);
|
||||
props_.matrix(n, pw.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
|
||||
const int phase_ind = pu_.phase_pos[Water];
|
||||
const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
|
||||
ADB::M db_diag = spdiag(dmatrix.col(column));
|
||||
const int num_blocks = pw.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = db_diag * pw.derivative()[block];
|
||||
}
|
||||
return ADB::function(matrix.col(column), jacs);
|
||||
}
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB BlackoilPropsAd::bOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.value().size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
if (pu_.phase_used[Gas]) {
|
||||
// Faking a z with the right ratio:
|
||||
// rs = zg/zo
|
||||
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
|
||||
z.col(pu_.phase_pos[Gas]) = rs.value();
|
||||
}
|
||||
Block matrix(n, np*np);
|
||||
Block dmatrix(n, np*np);
|
||||
props_.matrix(n, po.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
|
||||
const int phase_ind = pu_.phase_pos[Oil];
|
||||
const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
|
||||
ADB::M db_diag = spdiag(dmatrix.col(column));
|
||||
const int num_blocks = po.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
// For now, we deliberately ignore the derivative with respect to rs,
|
||||
// since the BlackoilPropertiesInterface class does not evaluate it.
|
||||
// We would add to the next line: + db_drs_diag * rs.derivative()[block]
|
||||
jacs[block] = db_diag * po.derivative()[block];
|
||||
}
|
||||
return ADB::function(matrix.col(column), jacs);
|
||||
}
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB BlackoilPropsAd::bGas(const ADB& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!pu_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.value().size() == n);
|
||||
const int np = props_.numPhases();
|
||||
Block z = Block::Zero(n, np);
|
||||
Block matrix(n, np*np);
|
||||
Block dmatrix(n, np*np);
|
||||
props_.matrix(n, pg.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
|
||||
const int phase_ind = pu_.phase_pos[Gas];
|
||||
const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
|
||||
ADB::M db_diag = spdiag(dmatrix.col(column));
|
||||
const int num_blocks = pg.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = db_diag * pg.derivative()[block];
|
||||
}
|
||||
return ADB::function(matrix.col(column), jacs);
|
||||
}
|
||||
|
||||
|
||||
// ------ Rs bubble point curve ------
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
V BlackoilPropsAd::rsMax(const V& po,
|
||||
const Cells& cells) const
|
||||
{
|
||||
// Suppress warning about "unused parameters".
|
||||
static_cast<void>(po);
|
||||
static_cast<void>(cells);
|
||||
|
||||
OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
|
||||
}
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
ADB BlackoilPropsAd::rsMax(const ADB& po,
|
||||
const Cells& cells) const
|
||||
{
|
||||
// Suppress warning about "unused parameters".
|
||||
static_cast<void>(po);
|
||||
static_cast<void>(cells);
|
||||
|
||||
OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
|
||||
}
|
||||
|
||||
// ------ Relative permeability ------
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<V> BlackoilPropsAd::relperm(const V& sw,
|
||||
const V& so,
|
||||
const V& sg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int n = cells.size();
|
||||
const int np = props_.numPhases();
|
||||
Block s_all(n, np);
|
||||
if (pu_.phase_used[Water]) {
|
||||
assert(sw.size() == n);
|
||||
s_all.col(pu_.phase_pos[Water]) = sw;
|
||||
}
|
||||
if (pu_.phase_used[Oil]) {
|
||||
assert(so.size() == n);
|
||||
s_all.col(pu_.phase_pos[Oil]) = so;
|
||||
}
|
||||
if (pu_.phase_used[Gas]) {
|
||||
assert(sg.size() == n);
|
||||
s_all.col(pu_.phase_pos[Gas]) = sg;
|
||||
}
|
||||
Block kr(n, np);
|
||||
props_.relperm(n, s_all.data(), cells.data(), kr.data(), 0);
|
||||
std::vector<V> relperms;
|
||||
relperms.reserve(3);
|
||||
for (int phase = 0; phase < 3; ++phase) {
|
||||
if (pu_.phase_used[phase]) {
|
||||
relperms.emplace_back(kr.col(pu_.phase_pos[phase]));
|
||||
} else {
|
||||
relperms.emplace_back();
|
||||
}
|
||||
}
|
||||
return relperms;
|
||||
}
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<ADB> BlackoilPropsAd::relperm(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int n = cells.size();
|
||||
const int np = props_.numPhases();
|
||||
Block s_all(n, np);
|
||||
if (pu_.phase_used[Water]) {
|
||||
assert(sw.value().size() == n);
|
||||
s_all.col(pu_.phase_pos[Water]) = sw.value();
|
||||
}
|
||||
if (pu_.phase_used[Oil]) {
|
||||
assert(so.value().size() == n);
|
||||
s_all.col(pu_.phase_pos[Oil]) = so.value();
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "BlackoilPropsAd::relperm() assumes oil phase is active.");
|
||||
}
|
||||
if (pu_.phase_used[Gas]) {
|
||||
assert(sg.value().size() == n);
|
||||
s_all.col(pu_.phase_pos[Gas]) = sg.value();
|
||||
}
|
||||
Block kr(n, np);
|
||||
Block dkr(n, np*np);
|
||||
props_.relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
|
||||
const int num_blocks = so.numBlocks();
|
||||
std::vector<ADB> relperms;
|
||||
relperms.reserve(3);
|
||||
typedef const ADB* ADBPtr;
|
||||
ADBPtr s[3] = { &sw, &so, &sg };
|
||||
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
||||
if (pu_.phase_used[phase1]) {
|
||||
const int phase1_pos = pu_.phase_pos[phase1];
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
|
||||
}
|
||||
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
||||
if (!pu_.phase_used[phase2]) {
|
||||
continue;
|
||||
}
|
||||
const int phase2_pos = pu_.phase_pos[phase2];
|
||||
// Assemble dkr1/ds2.
|
||||
const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
||||
ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
|
||||
}
|
||||
}
|
||||
relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
|
||||
} else {
|
||||
relperms.emplace_back(ADB::null());
|
||||
}
|
||||
}
|
||||
return relperms;
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
248
opm/polymer/fullyimplicit/BlackoilPropsAd.hpp
Normal file
248
opm/polymer/fullyimplicit/BlackoilPropsAd.hpp
Normal file
@ -0,0 +1,248 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BLACKOILPROPSAD_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILPROPSAD_HEADER_INCLUDED
|
||||
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
class BlackoilPropertiesInterface;
|
||||
|
||||
/// This class implements the AD-adapted fluid interface for
|
||||
/// three-phase black-oil.
|
||||
///
|
||||
/// It is implemented by wrapping a BlackoilPropertiesInterface
|
||||
/// object (the interface class defined in opm-core) and calling
|
||||
/// its methods. This class does not implement rsMax() because the
|
||||
/// required information is not available when wrapping a
|
||||
/// BlackoilPropertiesInterface. Consequently, class
|
||||
/// BlackoilPropsAd cannot be used to simulate problems involving
|
||||
/// miscibility.
|
||||
///
|
||||
/// Most methods are available in two overloaded versions, one
|
||||
/// taking a constant vector and returning the same, and one
|
||||
/// taking an AD type and returning the same. Derivatives are not
|
||||
/// returned separately by any method, only implicitly with the AD
|
||||
/// version of the methods.
|
||||
class BlackoilPropsAd : public BlackoilPropsAdInterface
|
||||
{
|
||||
public:
|
||||
/// Constructor wrapping an opm-core black oil interface.
|
||||
explicit BlackoilPropsAd(const BlackoilPropertiesInterface& props);
|
||||
|
||||
////////////////////////////
|
||||
// Rock interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return D, the number of spatial dimensions.
|
||||
int numDimensions() const;
|
||||
|
||||
/// \return N, the number of cells.
|
||||
int numCells() const;
|
||||
|
||||
/// \return Array of N porosity values.
|
||||
const double* porosity() const;
|
||||
|
||||
/// \return Array of ND^2 permeability values.
|
||||
/// The D^2 permeability values for a cell are organized as a matrix,
|
||||
/// which is symmetric (so ordering does not matter).
|
||||
const double* permeability() const;
|
||||
|
||||
|
||||
////////////////////////////
|
||||
// Fluid interface //
|
||||
////////////////////////////
|
||||
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef std::vector<int> Cells;
|
||||
|
||||
/// \return Number of active phases (also the number of components).
|
||||
virtual int numPhases() const;
|
||||
|
||||
/// \return Object describing the active phases.
|
||||
virtual PhaseUsage phaseUsage() const;
|
||||
|
||||
// ------ Canonical named indices for each phase ------
|
||||
|
||||
/// Canonical named indices for each phase.
|
||||
enum PhaseIndex { Water = 0, Oil = 1, Gas = 2 };
|
||||
|
||||
|
||||
// ------ Density ------
|
||||
|
||||
/// Densities of stock components at surface conditions.
|
||||
/// \return Array of 3 density values.
|
||||
const double* surfaceDensity() const;
|
||||
|
||||
|
||||
// ------ Viscosity ------
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V muWat(const V& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V muOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V muGas(const V& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muWat(const ADB& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muGas(const ADB& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
// ------ Formation volume factor (b) ------
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V bWat(const V& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V bOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V bGas(const V& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bWat(const ADB& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bGas(const ADB& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
// ------ Rs bubble point curve ------
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
V rsMax(const V& po,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
ADB rsMax(const ADB& po,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
// ------ Relative permeability ------
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<V> relperm(const V& sw,
|
||||
const V& so,
|
||||
const V& sg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<ADB> relperm(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const;
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
PhaseUsage pu_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILPROPSAD_HEADER_INCLUDED
|
||||
682
opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.cpp
Normal file
682
opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.cpp
Normal file
@ -0,0 +1,682 @@
|
||||
/*
|
||||
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 <config.h>
|
||||
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp>
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
|
||||
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtDead.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
|
||||
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
// Making these typedef to make the code more readable.
|
||||
typedef BlackoilPropsAdFromDeck::ADB ADB;
|
||||
typedef BlackoilPropsAdFromDeck::V V;
|
||||
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
|
||||
enum { Aqua = BlackoilPhases::Aqua,
|
||||
Liquid = BlackoilPhases::Liquid,
|
||||
Vapour = BlackoilPhases::Vapour };
|
||||
|
||||
/// Constructor wrapping an opm-core black oil interface.
|
||||
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
|
||||
const UnstructuredGrid& grid,
|
||||
const bool init_rock)
|
||||
{
|
||||
if (init_rock){
|
||||
rock_.init(deck, grid);
|
||||
}
|
||||
const int samples = 0;
|
||||
const int region_number = 0;
|
||||
|
||||
phase_usage_ = phaseUsageFromDeck(deck);
|
||||
|
||||
// Surface densities. Accounting for different orders in eclipse and our code.
|
||||
if (deck.hasField("DENSITY")) {
|
||||
const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
|
||||
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
|
||||
if (phase_usage_.phase_used[Aqua]) {
|
||||
densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water];
|
||||
}
|
||||
if (phase_usage_.phase_used[Vapour]) {
|
||||
densities_[phase_usage_.phase_pos[Vapour]] = d[ECL_gas];
|
||||
}
|
||||
if (phase_usage_.phase_used[Liquid]) {
|
||||
densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil];
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
|
||||
}
|
||||
|
||||
// Set the properties.
|
||||
props_.resize(phase_usage_.num_phases);
|
||||
// Water PVT
|
||||
if (phase_usage_.phase_used[Aqua]) {
|
||||
if (deck.hasField("PVTW")) {
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(deck.getPVTW().pvtw_));
|
||||
} else {
|
||||
// Eclipse 100 default.
|
||||
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
|
||||
}
|
||||
}
|
||||
// Oil PVT
|
||||
if (phase_usage_.phase_used[Liquid]) {
|
||||
if (deck.hasField("PVDO")) {
|
||||
if (samples > 0) {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_, samples));
|
||||
} else {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_));
|
||||
}
|
||||
} else if (deck.hasField("PVTO")) {
|
||||
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_));
|
||||
} else if (deck.hasField("PVCDO")) {
|
||||
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(deck.getPVCDO().pvcdo_));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
|
||||
}
|
||||
}
|
||||
// Gas PVT
|
||||
if (phase_usage_.phase_used[Vapour]) {
|
||||
if (deck.hasField("PVDG")) {
|
||||
if (samples > 0) {
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_, samples));
|
||||
} else {
|
||||
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_));
|
||||
}
|
||||
// } else if (deck.hasField("PVTG")) {
|
||||
// props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
|
||||
}
|
||||
}
|
||||
|
||||
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
|
||||
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
|
||||
satprops_.reset(ptr);
|
||||
ptr->init(deck, grid, -1);
|
||||
|
||||
if (phase_usage_.num_phases != satprops_->numPhases()) {
|
||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
|
||||
"Inconsistent number of phases in pvt data (" << phase_usage_.num_phases
|
||||
<< ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////
|
||||
// Rock interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return D, the number of spatial dimensions.
|
||||
int BlackoilPropsAdFromDeck::numDimensions() const
|
||||
{
|
||||
return rock_.numDimensions();
|
||||
}
|
||||
|
||||
/// \return N, the number of cells.
|
||||
int BlackoilPropsAdFromDeck::numCells() const
|
||||
{
|
||||
return rock_.numCells();
|
||||
}
|
||||
|
||||
/// \return Array of N porosity values.
|
||||
const double* BlackoilPropsAdFromDeck::porosity() const
|
||||
{
|
||||
return rock_.porosity();
|
||||
}
|
||||
|
||||
/// \return Array of ND^2 permeability values.
|
||||
/// The D^2 permeability values for a cell are organized as a matrix,
|
||||
/// which is symmetric (so ordering does not matter).
|
||||
const double* BlackoilPropsAdFromDeck::permeability() const
|
||||
{
|
||||
return rock_.permeability();
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////
|
||||
// Fluid interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return Number of active phases (also the number of components).
|
||||
int BlackoilPropsAdFromDeck::numPhases() const
|
||||
{
|
||||
return phase_usage_.num_phases;
|
||||
}
|
||||
|
||||
/// \return Object describing the active phases.
|
||||
PhaseUsage BlackoilPropsAdFromDeck::phaseUsage() const
|
||||
{
|
||||
return phase_usage_;
|
||||
}
|
||||
|
||||
// ------ Density ------
|
||||
|
||||
/// Densities of stock components at surface conditions.
|
||||
/// \return Array of 3 density values.
|
||||
const double* BlackoilPropsAdFromDeck::surfaceDensity() const
|
||||
{
|
||||
return densities_;
|
||||
}
|
||||
|
||||
|
||||
// ------ Viscosity ------
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V BlackoilPropsAdFromDeck::muWat(const V& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.size() == n);
|
||||
V mu(n);
|
||||
V dmudp(n);
|
||||
V dmudr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Water]]->mu(n, pw.data(), rs,
|
||||
mu.data(), dmudp.data(), dmudr.data());
|
||||
return mu;
|
||||
}
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V BlackoilPropsAdFromDeck::muOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
V mu(n);
|
||||
V dmudp(n);
|
||||
V dmudr(n);
|
||||
|
||||
props_[phase_usage_.phase_pos[Oil]]->mu(n, po.data(), rs.data(),
|
||||
mu.data(), dmudp.data(), dmudr.data());
|
||||
return mu;
|
||||
}
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V BlackoilPropsAdFromDeck::muGas(const V& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.size() == n);
|
||||
V mu(n);
|
||||
V dmudp(n);
|
||||
V dmudr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.data(), rs,
|
||||
mu.data(), dmudp.data(), dmudr.data());
|
||||
return mu;
|
||||
}
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB BlackoilPropsAdFromDeck::muWat(const ADB& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.size() == n);
|
||||
V mu(n);
|
||||
V dmudp(n);
|
||||
V dmudr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Water]]->mu(n, pw.value().data(), rs,
|
||||
mu.data(), dmudp.data(), dmudr.data());
|
||||
ADB::M dmudp_diag = spdiag(dmudp);
|
||||
const int num_blocks = pw.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dmudp_diag * pw.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu, jacs);
|
||||
}
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB BlackoilPropsAdFromDeck::muOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
V mu(n);
|
||||
V dmudp(n);
|
||||
V dmudr(n);
|
||||
|
||||
props_[phase_usage_.phase_pos[Oil]]->mu(n, po.value().data(), rs.value().data(),
|
||||
mu.data(), dmudp.data(), dmudr.data());
|
||||
|
||||
ADB::M dmudp_diag = spdiag(dmudp);
|
||||
ADB::M dmudr_diag = spdiag(dmudr);
|
||||
const int num_blocks = po.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dmudp_diag * po.derivative()[block] + dmudr_diag * rs.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu, jacs);
|
||||
}
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB BlackoilPropsAdFromDeck::muGas(const ADB& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.value().size() == n);
|
||||
V mu(n);
|
||||
V dmudp(n);
|
||||
V dmudr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.value().data(), rs,
|
||||
mu.data(), dmudp.data(), dmudr.data());
|
||||
|
||||
ADB::M dmudp_diag = spdiag(dmudp);
|
||||
const int num_blocks = pg.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dmudp_diag * pg.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu, jacs);
|
||||
}
|
||||
|
||||
|
||||
// ------ Formation volume factor (b) ------
|
||||
|
||||
// These methods all call the matrix() method, after which the variable
|
||||
// (also) called 'matrix' contains, in each row, the A = RB^{-1} matrix for
|
||||
// a cell. For three-phase black oil:
|
||||
// A = [ bw 0 0
|
||||
// 0 bo 0
|
||||
// 0 b0*rs bw ]
|
||||
// Where b = B^{-1}.
|
||||
// Therefore, we extract the correct diagonal element, and are done.
|
||||
// When we need the derivatives (w.r.t. p, since we don't do w.r.t. rs),
|
||||
// we also get the following derivative matrix:
|
||||
// A = [ dbw 0 0
|
||||
// 0 dbo 0
|
||||
// 0 db0*rs dbw ]
|
||||
// Again, we just extract a diagonal element.
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V BlackoilPropsAdFromDeck::bWat(const V& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call bWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.size() == n);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
V dbdr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Water]]->b(n, pw.data(), rs,
|
||||
b.data(), dbdp.data(), dbdr.data());
|
||||
|
||||
return b;
|
||||
}
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V BlackoilPropsAdFromDeck::bOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call bOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
V dbdr(n);
|
||||
|
||||
props_[phase_usage_.phase_pos[Oil]]->b(n, po.data(), rs.data(),
|
||||
b.data(), dbdp.data(), dbdr.data());
|
||||
|
||||
return b;
|
||||
}
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V BlackoilPropsAdFromDeck::bGas(const V& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call bGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.size() == n);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
V dbdr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Gas]]->b(n, pg.data(), rs,
|
||||
b.data(), dbdp.data(), dbdr.data());
|
||||
|
||||
return b;
|
||||
}
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB BlackoilPropsAdFromDeck::bWat(const ADB& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Water]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pw.size() == n);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
V dbdr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Water]]->b(n, pw.value().data(), rs,
|
||||
b.data(), dbdp.data(), dbdr.data());
|
||||
|
||||
ADB::M dbdp_diag = spdiag(dbdp);
|
||||
const int num_blocks = pw.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dbdp_diag * pw.derivative()[block];
|
||||
}
|
||||
return ADB::function(b, jacs);
|
||||
}
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB BlackoilPropsAdFromDeck::bOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
V dbdr(n);
|
||||
|
||||
props_[phase_usage_.phase_pos[Oil]]->b(n, po.value().data(), rs.value().data(),
|
||||
b.data(), dbdp.data(), dbdr.data());
|
||||
|
||||
ADB::M dbdp_diag = spdiag(dbdp);
|
||||
ADB::M dbdr_diag = spdiag(dbdr);
|
||||
const int num_blocks = po.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dbdp_diag * po.derivative()[block] + dbdr_diag * rs.derivative()[block];
|
||||
}
|
||||
return ADB::function(b, jacs);
|
||||
}
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB BlackoilPropsAdFromDeck::bGas(const ADB& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Gas]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(pg.size() == n);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
V dbdr(n);
|
||||
const double* rs = 0;
|
||||
|
||||
props_[phase_usage_.phase_pos[Gas]]->b(n, pg.value().data(), rs,
|
||||
b.data(), dbdp.data(), dbdr.data());
|
||||
|
||||
ADB::M dbdp_diag = spdiag(dbdp);
|
||||
const int num_blocks = pg.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dbdp_diag * pg.derivative()[block];
|
||||
}
|
||||
return ADB::function(b, jacs);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// ------ Rs bubble point curve ------
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
V BlackoilPropsAdFromDeck::rsMax(const V& po,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call rsMax(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
V rbub(n);
|
||||
V drbubdp(n);
|
||||
props_[Oil]->rbub(n, po.data(), rbub.data(), drbubdp.data());
|
||||
return rbub;
|
||||
}
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
ADB BlackoilPropsAdFromDeck::rsMax(const ADB& po,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (!phase_usage_.phase_used[Oil]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot call rsMax(): oil phase not present.");
|
||||
}
|
||||
const int n = cells.size();
|
||||
assert(po.size() == n);
|
||||
V rbub(n);
|
||||
V drbubdp(n);
|
||||
props_[Oil]->rbub(n, po.value().data(), rbub.data(), drbubdp.data());
|
||||
ADB::M drbubdp_diag = spdiag(drbubdp);
|
||||
const int num_blocks = po.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = drbubdp_diag * po.derivative()[block];
|
||||
}
|
||||
return ADB::function(rbub, jacs);
|
||||
}
|
||||
|
||||
// ------ Relative permeability ------
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<V> BlackoilPropsAdFromDeck::relperm(const V& sw,
|
||||
const V& so,
|
||||
const V& sg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int n = cells.size();
|
||||
const int np = numPhases();
|
||||
Block s_all(n, np);
|
||||
if (phase_usage_.phase_used[Water]) {
|
||||
assert(sw.size() == n);
|
||||
s_all.col(phase_usage_.phase_pos[Water]) = sw;
|
||||
}
|
||||
if (phase_usage_.phase_used[Oil]) {
|
||||
assert(so.size() == n);
|
||||
s_all.col(phase_usage_.phase_pos[Oil]) = so;
|
||||
}
|
||||
if (phase_usage_.phase_used[Gas]) {
|
||||
assert(sg.size() == n);
|
||||
s_all.col(phase_usage_.phase_pos[Gas]) = sg;
|
||||
}
|
||||
Block kr(n, np);
|
||||
satprops_->relperm(n, s_all.data(), cells.data(), kr.data(), 0);
|
||||
std::vector<V> relperms;
|
||||
relperms.reserve(3);
|
||||
for (int phase = 0; phase < 3; ++phase) {
|
||||
if (phase_usage_.phase_used[phase]) {
|
||||
relperms.emplace_back(kr.col(phase_usage_.phase_pos[phase]));
|
||||
} else {
|
||||
relperms.emplace_back();
|
||||
}
|
||||
}
|
||||
return relperms;
|
||||
}
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<ADB> BlackoilPropsAdFromDeck::relperm(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int n = cells.size();
|
||||
const int np = numPhases();
|
||||
Block s_all(n, np);
|
||||
if (phase_usage_.phase_used[Water]) {
|
||||
assert(sw.value().size() == n);
|
||||
s_all.col(phase_usage_.phase_pos[Water]) = sw.value();
|
||||
}
|
||||
if (phase_usage_.phase_used[Oil]) {
|
||||
assert(so.value().size() == n);
|
||||
s_all.col(phase_usage_.phase_pos[Oil]) = so.value();
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
||||
}
|
||||
if (phase_usage_.phase_used[Gas]) {
|
||||
assert(sg.value().size() == n);
|
||||
s_all.col(phase_usage_.phase_pos[Gas]) = sg.value();
|
||||
}
|
||||
Block kr(n, np);
|
||||
Block dkr(n, np*np);
|
||||
satprops_->relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
|
||||
const int num_blocks = so.numBlocks();
|
||||
std::vector<ADB> relperms;
|
||||
relperms.reserve(3);
|
||||
typedef const ADB* ADBPtr;
|
||||
ADBPtr s[3] = { &sw, &so, &sg };
|
||||
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
||||
if (phase_usage_.phase_used[phase1]) {
|
||||
const int phase1_pos = phase_usage_.phase_pos[phase1];
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
|
||||
}
|
||||
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
||||
if (!phase_usage_.phase_used[phase2]) {
|
||||
continue;
|
||||
}
|
||||
const int phase2_pos = phase_usage_.phase_pos[phase2];
|
||||
// Assemble dkr1/ds2.
|
||||
const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
||||
ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
|
||||
}
|
||||
}
|
||||
relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
|
||||
} else {
|
||||
relperms.emplace_back(ADB::null());
|
||||
}
|
||||
}
|
||||
return relperms;
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
253
opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp
Normal file
253
opm/polymer/fullyimplicit/BlackoilPropsAdFromDeck.hpp
Normal file
@ -0,0 +1,253 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BLACKOILPROPSADFROMDECK_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILPROPSADFROMDECK_HEADER_INCLUDED
|
||||
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
|
||||
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/props/rock/RockFromDeck.hpp>
|
||||
|
||||
#include <boost/shared_ptr.hpp>
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
class SinglePvtInterface;
|
||||
|
||||
/// This class implements the AD-adapted fluid interface for
|
||||
/// three-phase black-oil. It requires an input deck from which it
|
||||
/// reads all relevant property data.
|
||||
///
|
||||
/// Most methods are available in two overloaded versions, one
|
||||
/// taking a constant vector and returning the same, and one
|
||||
/// taking an AD type and returning the same. Derivatives are not
|
||||
/// returned separately by any method, only implicitly with the AD
|
||||
/// version of the methods.
|
||||
class BlackoilPropsAdFromDeck : public BlackoilPropsAdInterface
|
||||
{
|
||||
public:
|
||||
/// Constructor wrapping an opm-core black oil interface.
|
||||
BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
|
||||
const UnstructuredGrid& grid,
|
||||
const bool init_rock = true );
|
||||
|
||||
////////////////////////////
|
||||
// Rock interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return D, the number of spatial dimensions.
|
||||
int numDimensions() const;
|
||||
|
||||
/// \return N, the number of cells.
|
||||
int numCells() const;
|
||||
|
||||
/// \return Array of N porosity values.
|
||||
const double* porosity() const;
|
||||
|
||||
/// \return Array of ND^2 permeability values.
|
||||
/// The D^2 permeability values for a cell are organized as a matrix,
|
||||
/// which is symmetric (so ordering does not matter).
|
||||
const double* permeability() const;
|
||||
|
||||
|
||||
////////////////////////////
|
||||
// Fluid interface //
|
||||
////////////////////////////
|
||||
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef std::vector<int> Cells;
|
||||
|
||||
/// \return Number of active phases (also the number of components).
|
||||
int numPhases() const;
|
||||
|
||||
/// \return Object describing the active phases.
|
||||
PhaseUsage phaseUsage() const;
|
||||
|
||||
// ------ Canonical named indices for each phase ------
|
||||
|
||||
/// Canonical named indices for each phase.
|
||||
enum PhaseIndex { Water = 0, Oil = 1, Gas = 2 };
|
||||
|
||||
|
||||
// ------ Density ------
|
||||
|
||||
/// Densities of stock components at surface conditions.
|
||||
/// \return Array of 3 density values.
|
||||
const double* surfaceDensity() const;
|
||||
|
||||
|
||||
// ------ Viscosity ------
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V muWat(const V& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V muOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
V muGas(const V& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muWat(const ADB& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muGas(const ADB& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
// ------ Formation volume factor (b) ------
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V bWat(const V& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V bOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
V bGas(const V& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bWat(const ADB& pw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bGas(const ADB& pg,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
// ------ Rs bubble point curve ------
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
V rsMax(const V& po,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
ADB rsMax(const ADB& po,
|
||||
const Cells& cells) const;
|
||||
|
||||
|
||||
// ------ Relative permeability ------
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<V> relperm(const V& sw,
|
||||
const V& so,
|
||||
const V& sg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
std::vector<ADB> relperm(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const;
|
||||
|
||||
private:
|
||||
RockFromDeck rock_;
|
||||
boost::scoped_ptr<SaturationPropsInterface> satprops_;
|
||||
PhaseUsage phase_usage_;
|
||||
std::vector<boost::shared_ptr<SinglePvtInterface> > props_;
|
||||
double densities_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILPROPSADFROMDECK_HEADER_INCLUDED
|
||||
24
opm/polymer/fullyimplicit/BlackoilPropsAdInterface.cpp
Normal file
24
opm/polymer/fullyimplicit/BlackoilPropsAdInterface.cpp
Normal file
@ -0,0 +1,24 @@
|
||||
/*
|
||||
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/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
|
||||
Opm::BlackoilPropsAdInterface::~BlackoilPropsAdInterface()
|
||||
{
|
||||
}
|
||||
250
opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp
Normal file
250
opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp
Normal file
@ -0,0 +1,250 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BLACKOILPROPSADINTERFACE_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILPROPSADINTERFACE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/// This class is intended to present a fluid interface for
|
||||
/// three-phase black-oil that is easy to use with the AD-using
|
||||
/// simulators.
|
||||
///
|
||||
/// Most methods are available in two overloaded versions, one
|
||||
/// taking a constant vector and returning the same, and one
|
||||
/// taking an AD type and returning the same. Derivatives are not
|
||||
/// returned separately by any method, only implicitly with the AD
|
||||
/// version of the methods.
|
||||
class BlackoilPropsAdInterface
|
||||
{
|
||||
public:
|
||||
/// Virtual destructor for inheritance.
|
||||
virtual ~BlackoilPropsAdInterface();
|
||||
|
||||
////////////////////////////
|
||||
// Rock interface //
|
||||
////////////////////////////
|
||||
|
||||
/// \return D, the number of spatial dimensions.
|
||||
virtual int numDimensions() const = 0;
|
||||
|
||||
/// \return N, the number of cells.
|
||||
virtual int numCells() const = 0;
|
||||
|
||||
/// \return Array of N porosity values.
|
||||
virtual const double* porosity() const = 0;
|
||||
|
||||
/// \return Array of ND^2 permeability values.
|
||||
/// The D^2 permeability values for a cell are organized as a matrix,
|
||||
/// which is symmetric (so ordering does not matter).
|
||||
virtual const double* permeability() const = 0;
|
||||
|
||||
|
||||
////////////////////////////
|
||||
// Fluid interface //
|
||||
////////////////////////////
|
||||
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef std::vector<int> Cells;
|
||||
|
||||
/// \return Number of active phases (also the number of components).
|
||||
virtual int numPhases() const = 0;
|
||||
|
||||
/// \return Object describing the active phases.
|
||||
virtual PhaseUsage phaseUsage() const = 0;
|
||||
|
||||
// ------ Canonical named indices for each phase ------
|
||||
|
||||
/// Canonical named indices for each phase.
|
||||
enum PhaseIndex { Water = 0, Oil = 1, Gas = 2 };
|
||||
|
||||
// ------ Density ------
|
||||
|
||||
/// Densities of stock components at surface conditions.
|
||||
/// \return Array of 3 density values.
|
||||
virtual const double* surfaceDensity() const = 0;
|
||||
|
||||
|
||||
// ------ Viscosity ------
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
virtual
|
||||
V muWat(const V& pw,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
virtual
|
||||
V muOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
virtual
|
||||
V muGas(const V& pg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Water viscosity.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
virtual
|
||||
ADB muWat(const ADB& pw,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Oil viscosity.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
virtual
|
||||
ADB muOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Gas viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
virtual
|
||||
ADB muGas(const ADB& pg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
|
||||
// ------ Formation volume factor (b) ------
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
virtual
|
||||
V bWat(const V& pw,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
virtual
|
||||
V bOil(const V& po,
|
||||
const V& rs,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
virtual
|
||||
V bGas(const V& pg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Water formation volume factor.
|
||||
/// \param[in] pw Array of n water pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
virtual
|
||||
ADB bWat(const ADB& pw,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Oil formation volume factor.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] rs Array of n gas solution factor values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
virtual
|
||||
ADB bOil(const ADB& po,
|
||||
const ADB& rs,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Gas formation volume factor.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
virtual
|
||||
ADB bGas(const ADB& pg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
|
||||
// ------ Rs bubble point curve ------
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
virtual
|
||||
V rsMax(const V& po,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Bubble point curve for Rs as function of oil pressure.
|
||||
/// \param[in] po Array of n oil pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n bubble point values for Rs.
|
||||
virtual
|
||||
ADB rsMax(const ADB& po,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
// ------ Relative permeability ------
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
virtual
|
||||
std::vector<V> relperm(const V& sw,
|
||||
const V& so,
|
||||
const V& sg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
/// Relative permeabilities for all phases.
|
||||
/// \param[in] sw Array of n water saturation values.
|
||||
/// \param[in] so Array of n oil saturation values.
|
||||
/// \param[in] sg Array of n gas saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
||||
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
||||
virtual
|
||||
std::vector<ADB> relperm(const ADB& sw,
|
||||
const ADB& so,
|
||||
const ADB& sg,
|
||||
const Cells& cells) const = 0;
|
||||
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILPROPSADINTERFACE_HEADER_INCLUDED
|
||||
@ -0,0 +1,964 @@
|
||||
/*
|
||||
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/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
|
||||
|
||||
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
#include <opm/polymer/fullyimplicit/GeoProps.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/linalg/LinearSolverInterface.hpp>
|
||||
#include <opm/core/props/rock/RockCompressibility.hpp>
|
||||
#include <opm/polymer/PolymerBlackoilState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
|
||||
// A debugging utility.
|
||||
#define DUMP(foo) \
|
||||
do { \
|
||||
std::cout << "==========================================\n" \
|
||||
<< #foo ":\n" \
|
||||
<< collapseJacs(foo) << std::endl; \
|
||||
} while (0)
|
||||
|
||||
|
||||
|
||||
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;
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class GeoProps>
|
||||
AutoDiffBlock<double>::M
|
||||
gravityOperator(const UnstructuredGrid& grid,
|
||||
const HelperOps& ops ,
|
||||
const GeoProps& geo )
|
||||
{
|
||||
const int nc = grid.number_of_cells;
|
||||
|
||||
std::vector<int> f2hf(2 * grid.number_of_faces, -1);
|
||||
for (int c = 0, i = 0; c < nc; ++c) {
|
||||
for (; i < grid.cell_facepos[c + 1]; ++i) {
|
||||
const int f = grid.cell_faces[ i ];
|
||||
const int p = 0 + (grid.face_cells[2*f + 0] != c);
|
||||
|
||||
f2hf[2*f + p] = i;
|
||||
}
|
||||
}
|
||||
|
||||
typedef AutoDiffBlock<double>::V V;
|
||||
typedef AutoDiffBlock<double>::M M;
|
||||
|
||||
const V& gpot = geo.gravityPotential();
|
||||
const V& trans = geo.transmissibility();
|
||||
|
||||
const HelperOps::IFaces::Index ni = ops.internal_faces.size();
|
||||
|
||||
typedef Eigen::Triplet<double> Tri;
|
||||
std::vector<Tri> grav; grav.reserve(2 * ni);
|
||||
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
|
||||
const int f = ops.internal_faces[ i ];
|
||||
const int c1 = grid.face_cells[2*f + 0];
|
||||
const int c2 = grid.face_cells[2*f + 1];
|
||||
|
||||
assert ((c1 >= 0) && (c2 >= 0));
|
||||
|
||||
const double dG1 = gpot[ f2hf[2*f + 0] ];
|
||||
const double dG2 = gpot[ f2hf[2*f + 1] ];
|
||||
const double t = trans[ f ];
|
||||
|
||||
grav.push_back(Tri(i, c1, t * dG1));
|
||||
grav.push_back(Tri(i, c2, - t * dG2));
|
||||
}
|
||||
|
||||
M G(ni, nc); G.setFromTriplets(grav.begin(), grav.end());
|
||||
|
||||
return G;
|
||||
}
|
||||
|
||||
|
||||
|
||||
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
|
||||
|
||||
|
||||
|
||||
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const DerivedGeology& geo ,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
const PolymerPropsAd& polymer_props_ad,
|
||||
const Wells& wells,
|
||||
const LinearSolverInterface& linsolver)
|
||||
: grid_ (grid)
|
||||
, fluid_ (fluid)
|
||||
, geo_ (geo)
|
||||
, rock_comp_props_(rock_comp_props)
|
||||
, polymer_props_ad_(polymer_props_ad)
|
||||
, wells_ (wells)
|
||||
, linsolver_ (linsolver)
|
||||
, cells_ (buildAllCells(grid.number_of_cells))
|
||||
, ops_ (grid)
|
||||
, wops_ (wells)
|
||||
, grav_ (gravityOperator(grid_, ops_, geo_))
|
||||
, rq_ (fluid.numPhases() + 1)
|
||||
, residual_ ( { std::vector<ADB>(fluid.numPhases() + 1, ADB::null()),
|
||||
ADB::null(),
|
||||
ADB::null() } )
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
step(const double dt,
|
||||
PolymerBlackoilState& x ,
|
||||
WellState& xw,
|
||||
const std::vector<double>& polymer_inflow)
|
||||
{
|
||||
const V pvdt = geo_.poreVolume() / dt;
|
||||
|
||||
{
|
||||
const SolutionState state = constantState(x, xw);
|
||||
computeAccum(state, 0);
|
||||
}
|
||||
|
||||
const double atol = 1.0e-12;
|
||||
const double rtol = 5.0e-8;
|
||||
const int maxit = 15;
|
||||
|
||||
assemble(pvdt, x, xw, polymer_inflow);
|
||||
|
||||
const double r0 = residualNorm();
|
||||
int it = 0;
|
||||
std::cout << "\nIteration Residual\n"
|
||||
<< std::setw(9) << it << std::setprecision(9)
|
||||
<< std::setw(18) << r0 << std::endl;
|
||||
bool resTooLarge = r0 > atol;
|
||||
while (resTooLarge && (it < maxit)) {
|
||||
const V dx = solveJacobianSystem();
|
||||
|
||||
updateState(dx, x, xw);
|
||||
|
||||
assemble(pvdt, x, xw, polymer_inflow);
|
||||
|
||||
const double r = residualNorm();
|
||||
|
||||
resTooLarge = (r > atol) && (r > rtol*r0);
|
||||
|
||||
it += 1;
|
||||
std::cout << std::setw(9) << it << std::setprecision(9)
|
||||
<< std::setw(18) << r << std::endl;
|
||||
}
|
||||
|
||||
if (resTooLarge) {
|
||||
std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
|
||||
// OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FullyImplicitCompressiblePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant()
|
||||
: accum(2, ADB::null())
|
||||
, mflux( ADB::null())
|
||||
, b ( ADB::null())
|
||||
, head ( ADB::null())
|
||||
, mob ( ADB::null())
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FullyImplicitCompressiblePolymerSolver::SolutionState::SolutionState(const int np)
|
||||
: pressure ( ADB::null())
|
||||
, saturation(np, ADB::null())
|
||||
, concentration( ADB::null())
|
||||
, qs ( ADB::null())
|
||||
, bhp ( ADB::null())
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
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());
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FullyImplicitCompressiblePolymerSolver::SolutionState
|
||||
FullyImplicitCompressiblePolymerSolver::constantState(const PolymerBlackoilState& 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 (if water present)
|
||||
// polymer concentration
|
||||
// well rates per active phase and well
|
||||
// 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);
|
||||
|
||||
// Pressure.
|
||||
assert (not x.pressure().empty());
|
||||
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
|
||||
state.pressure = ADB::constant(p, bpat);
|
||||
|
||||
// Saturation.
|
||||
assert (not x.saturation().empty());
|
||||
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
V so = V::Ones(nc, 1);
|
||||
const V sw = s.col(0);
|
||||
so -= sw;
|
||||
state.saturation[0] = ADB::constant(sw, bpat);
|
||||
state.saturation[1] = ADB::constant(so, bpat);
|
||||
|
||||
// Concentration
|
||||
assert(not x.concentration().empty());
|
||||
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);
|
||||
|
||||
// Well 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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FullyImplicitCompressiblePolymerSolver::SolutionState
|
||||
FullyImplicitCompressiblePolymerSolver::variableState(const PolymerBlackoilState& x,
|
||||
const WellState& xw)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const int np = x.numPhases();
|
||||
|
||||
std::vector<V> vars0;
|
||||
|
||||
// Initial pressure.
|
||||
assert (not x.pressure().empty());
|
||||
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
|
||||
vars0.push_back(p);
|
||||
|
||||
// Initial saturation.
|
||||
assert (not x.saturation().empty());
|
||||
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
const V sw = s.col(0);
|
||||
vars0.push_back(sw);
|
||||
|
||||
// Initial concentration.
|
||||
assert (not x.concentration().empty());
|
||||
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().empty());
|
||||
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);
|
||||
|
||||
// Pressure.
|
||||
int nextvar = 0;
|
||||
state.pressure = vars[ nextvar++ ];
|
||||
|
||||
// Saturation.
|
||||
const std::vector<int>& bpat = vars[0].blockPattern();
|
||||
{
|
||||
ADB so = ADB::constant(V::Ones(nc, 1), bpat);
|
||||
ADB& sw = vars[ nextvar++ ];
|
||||
state.saturation[0] = sw;
|
||||
so = so - sw;
|
||||
state.saturation[1] = so;
|
||||
}
|
||||
|
||||
// Concentration.
|
||||
state.concentration = vars[nextvar++];
|
||||
|
||||
// Qs.
|
||||
state.qs = vars[ nextvar++ ];
|
||||
|
||||
// Bhp.
|
||||
state.bhp = vars[ nextvar++ ];
|
||||
|
||||
assert(nextvar == int(vars.size()));
|
||||
|
||||
return state;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void
|
||||
FullyImplicitCompressiblePolymerSolver::computeAccum(const SolutionState& state,
|
||||
const int aix )
|
||||
{
|
||||
|
||||
const ADB& press = state.pressure;
|
||||
const std::vector<ADB>& sat = state.saturation;
|
||||
const ADB& c = state.concentration;
|
||||
const ADB pv_mult = poroMult(press);
|
||||
|
||||
for (int phase = 0; phase < 2; ++phase) {
|
||||
rq_[phase].b = fluidReciprocFVF(phase, press, cells_);
|
||||
}
|
||||
rq_[0].accum[aix] = pv_mult * rq_[0].b * sat[0];
|
||||
rq_[1].accum[aix] = pv_mult * rq_[1].b * sat[1];
|
||||
rq_[2].accum[aix] = pv_mult * rq_[0].b * sat[0] * c;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
ADB
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
computeCmax(const ADB& c) const
|
||||
{
|
||||
const int nc = c.value().size();
|
||||
V cmax(nc);
|
||||
|
||||
for (int i = 0; i < nc; ++i) {
|
||||
cmax(i) = (cmax(i) > c.value()(i)) ? cmax(i) : c.value()(i);
|
||||
}
|
||||
|
||||
return ADB::constant(cmax, c.blockPattern());
|
||||
}
|
||||
|
||||
void
|
||||
FullyImplicitCompressiblePolymerSolver::
|
||||
assemble(const V& pvdt,
|
||||
const PolymerBlackoilState& x ,
|
||||
const WellState& xw,
|
||||
const std::vector<double>& polymer_inflow)
|
||||
{
|
||||
// Create the primary variables.
|
||||
const SolutionState state = variableState(x, xw);
|
||||
|
||||
// -------- Mass balance equations --------
|
||||
|
||||
// Compute b_p and the accumulation term b_p*s_p for each phase,
|
||||
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
|
||||
// These quantities are stored in rq_[phase].accum[1].
|
||||
// The corresponding accumulation terms from the start of
|
||||
// the timestep (b^0_p*s^0_p etc.) were already computed
|
||||
// in step() and stored in rq_[phase].accum[0].
|
||||
computeAccum(state, 1);
|
||||
|
||||
// Set up the common parts of the mass balance equations
|
||||
// for each active phase.
|
||||
const V trans = subset(geo_.transmissibility(), ops_.internal_faces);
|
||||
const std::vector<ADB> kr = computeRelPerm(state);
|
||||
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
|
||||
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]);
|
||||
const ADB mc = computeMc(state);
|
||||
computeMassFlux(trans, mc, kr[1], krw_eff, state);
|
||||
const double rho_rock = polymer_props_ad_.rockDensity();
|
||||
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], grid_.number_of_cells, 1);
|
||||
residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
|
||||
+ ops_.div*rq_[0].mflux;
|
||||
residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0])
|
||||
+ ops_.div*rq_[1].mflux;
|
||||
residual_.mass_balance[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) * (1. - dead_pore_vol)
|
||||
+ pvdt * rho_rock * (1. - phi) / phi * ads
|
||||
+ ops_.div*rq_[2].mflux;
|
||||
|
||||
|
||||
// -------- Extra (optional) sg or rs equation, and rs contributions to the mass balance equations --------
|
||||
|
||||
// Add the extra (flux) terms to the gas mass balance equations
|
||||
// from gas dissolved in the oil phase.
|
||||
// The extra terms in the accumulation part of the equation are already handled.
|
||||
|
||||
// -------- 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;
|
||||
const double* g = geo_.gravity();
|
||||
if (g) {
|
||||
// Guard against gravity in anything but last dimension.
|
||||
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) {
|
||||
const ADB cell_rho = fluidDensity(phase, state.pressure, 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, state.pressure, cells_);
|
||||
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, g ? g[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(rq_[0].mob, well_cells) + subset(rq_[1].mob, 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 = rq_[phase].mob;
|
||||
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*perf_b);
|
||||
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*perf_b, 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 and SURFACE_RATE 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);
|
||||
// DUMP(residual_.well_eq);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
V FullyImplicitCompressiblePolymerSolver::solveJacobianSystem() const
|
||||
{
|
||||
ADB mass_res = vertcat(residual_.mass_balance[0], residual_.mass_balance[1]);
|
||||
mass_res = vertcat(mass_res, residual_.mass_balance[2]);
|
||||
const ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
|
||||
const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
|
||||
// DUMP(total_residual);
|
||||
|
||||
const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
|
||||
|
||||
V dx(V::Zero(total_residual.size()));
|
||||
Opm::LinearSolverInterface::LinearSolverReport rep
|
||||
= linsolver_.solve(matr.rows(), matr.nonZeros(),
|
||||
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
|
||||
total_residual.value().data(), dx.data());
|
||||
if (!rep.converged) {
|
||||
OPM_THROW(std::runtime_error,
|
||||
"FullyImplicitCompressiblePolymerSolver::solveJacobianSystem(): "
|
||||
"Linear solver convergence failure.");
|
||||
}
|
||||
return dx;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
namespace {
|
||||
struct Chop01 {
|
||||
double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); }
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void FullyImplicitCompressiblePolymerSolver::updateState(const V& dx,
|
||||
PolymerBlackoilState& state,
|
||||
WellState& well_state) const
|
||||
{
|
||||
const int np = fluid_.numPhases();
|
||||
const int nc = grid_.number_of_cells;
|
||||
const int nw = wells_.number_of_wells;
|
||||
const V one = V::Constant(nc, 1.0);
|
||||
const V zero = V::Zero(nc);
|
||||
// Extract parts of dx corresponding to each part.
|
||||
const V dp = subset(dx, Span(nc));
|
||||
int varstart = nc;
|
||||
const V dsw =subset(dx, Span(nc, 1, varstart));
|
||||
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());
|
||||
|
||||
// Pressure update.
|
||||
const double dpmaxrel = 0.8;
|
||||
const V p_old = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
|
||||
const V absdpmax = dpmaxrel*p_old.abs();
|
||||
const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
|
||||
const V p = (p_old - dp_limited).max(zero);
|
||||
std::copy(&p[0], &p[0] + nc, state.pressure().begin());
|
||||
|
||||
// Saturation updates.
|
||||
const double dsmax = 0.3;
|
||||
const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
|
||||
V so = one;
|
||||
const V sw_old = s_old.col(0);
|
||||
const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax);
|
||||
const V sw = (sw_old - dsw_limited).unaryExpr(Chop01());
|
||||
so -= sw;
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
state.saturation()[c*np] = sw[c];
|
||||
state.saturation()[c*np + 1] = so[c];
|
||||
}
|
||||
|
||||
// Concentration updates.
|
||||
const V c_old = Eigen::Map<const V>(&state.concentration()[0], nc);
|
||||
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());
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
std::vector<ADB>
|
||||
FullyImplicitCompressiblePolymerSolver::computeRelPerm(const SolutionState& state) const
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const std::vector<int>& bpat = state.pressure.blockPattern();
|
||||
|
||||
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
||||
|
||||
const ADB sw = state.saturation[0];
|
||||
const ADB so = state.saturation[1];
|
||||
const ADB sg = null;
|
||||
return fluid_.relperm(sw, so, sg, cells_);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
std::vector<ADB>
|
||||
FullyImplicitCompressiblePolymerSolver::computeRelPermWells(const SolutionState& state,
|
||||
const DataBlock& well_s,
|
||||
const std::vector<int>& well_cells) const
|
||||
{
|
||||
const int nw = wells_.number_of_wells;
|
||||
const int nperf = wells_.well_connpos[nw];
|
||||
const std::vector<int>& bpat = state.pressure.blockPattern();
|
||||
|
||||
const ADB null = ADB::constant(V::Zero(nperf), bpat);
|
||||
|
||||
const ADB sw = state.saturation[0];
|
||||
const ADB so = state.saturation[1];
|
||||
const ADB sg = null;
|
||||
return fluid_.relperm(sw, so, sg, well_cells);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void
|
||||
FullyImplicitCompressiblePolymerSolver::computeMassFlux(
|
||||
const V& transi,
|
||||
const ADB& mc,
|
||||
const ADB& kro,
|
||||
const ADB& krw_eff,
|
||||
const SolutionState& state )
|
||||
{
|
||||
const ADB tr_mult = transMult(state.pressure);
|
||||
|
||||
const ADB mu_w = fluidViscosity(0, state.pressure, 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, state.pressure, cells_);
|
||||
rq_[1].mob = tr_mult * kro / mu_o;
|
||||
for (int phase = 0; phase < 2; ++phase) {
|
||||
const ADB rho = fluidDensity(phase, state.pressure, cells_);
|
||||
ADB& head = rq_[ phase ].head;
|
||||
// compute gravity potensial using the face average as in eclipse and MRST
|
||||
const ADB rhoavg = ops_.caver * rho;
|
||||
const ADB dp = ops_.ngrad * state.pressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
|
||||
head = transi*dp;
|
||||
UpwindSelector<double> upwind(grid_, ops_, head.value());
|
||||
const ADB& b = rq_[ phase ].b;
|
||||
const ADB& mob = rq_[ phase ].mob;
|
||||
rq_[ phase ].mflux = upwind.select(b * mob) * head;
|
||||
}
|
||||
rq_[2].b = rq_[0].b;
|
||||
rq_[2].head = rq_[0].head;
|
||||
UpwindSelector<double> upwind(grid_, ops_, rq_[2].head.value());
|
||||
rq_[2].mflux = upwind.select(rq_[2].b * rq_[2].mob) * rq_[2].head;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
double
|
||||
FullyImplicitCompressiblePolymerSolver::residualNorm() const
|
||||
{
|
||||
double r = 0;
|
||||
for (std::vector<ADB>::const_iterator
|
||||
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
|
||||
FullyImplicitCompressiblePolymerSolver::fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
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);
|
||||
case Oil: {
|
||||
return fluid_.muOil(p, null, cells);
|
||||
}
|
||||
default:
|
||||
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
ADB
|
||||
FullyImplicitCompressiblePolymerSolver::fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
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);
|
||||
case Oil: {
|
||||
return fluid_.bOil(p, null, cells);
|
||||
}
|
||||
default:
|
||||
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
ADB
|
||||
FullyImplicitCompressiblePolymerSolver::fluidDensity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const
|
||||
{
|
||||
const double* rhos = fluid_.surfaceDensity();
|
||||
ADB b = fluidReciprocFVF(phase, p, cells);
|
||||
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
|
||||
return rho;
|
||||
}
|
||||
|
||||
|
||||
// here mc means m(c) * c.
|
||||
ADB
|
||||
FullyImplicitCompressiblePolymerSolver::computeMc(const SolutionState& state) const
|
||||
{
|
||||
ADB c = state.concentration;
|
||||
return polymer_props_ad_.polymerWaterVelocityRatio(c);
|
||||
}
|
||||
|
||||
|
||||
|
||||
ADB
|
||||
FullyImplicitCompressiblePolymerSolver::poroMult(const ADB& p) const
|
||||
{
|
||||
const int n = p.size();
|
||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||
V pm(n);
|
||||
V dpm(n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
pm[i] = rock_comp_props_->poroMult(p.value()[i]);
|
||||
dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
|
||||
}
|
||||
ADB::M dpm_diag = spdiag(dpm);
|
||||
const int num_blocks = p.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dpm_diag * p.derivative()[block];
|
||||
}
|
||||
return ADB::function(pm, jacs);
|
||||
} else {
|
||||
return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
ADB
|
||||
FullyImplicitCompressiblePolymerSolver::transMult(const ADB& p) const
|
||||
{
|
||||
const int n = p.size();
|
||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||
V tm(n);
|
||||
V dtm(n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
tm[i] = rock_comp_props_->transMult(p.value()[i]);
|
||||
dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
|
||||
}
|
||||
ADB::M dtm_diag = spdiag(dtm);
|
||||
const int num_blocks = p.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
jacs[block] = dtm_diag * p.derivative()[block];
|
||||
}
|
||||
return ADB::function(tm, jacs);
|
||||
} else {
|
||||
return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@ -0,0 +1,229 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
|
||||
#define OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
|
||||
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
|
||||
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
#include <opm/polymer/PolymerProperties.hpp>
|
||||
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
|
||||
|
||||
struct UnstructuredGrid;
|
||||
struct Wells;
|
||||
|
||||
namespace Opm {
|
||||
|
||||
class DerivedGeology;
|
||||
class RockCompressibility;
|
||||
class LinearSolverInterface;
|
||||
class PolymerBlackoilState;
|
||||
class WellState;
|
||||
|
||||
/// A fully implicit solver for the black-oil problem.
|
||||
///
|
||||
/// The simulator is capable of handling three-phase problems
|
||||
/// where gas can be dissolved in oil (but not vice versa). It
|
||||
/// uses an industry-standard TPFA discretization with per-phase
|
||||
/// upwind weighting of mobilities.
|
||||
///
|
||||
/// It uses automatic differentiation via the class AutoDiffBlock
|
||||
/// to simplify assembly of the jacobian matrix.
|
||||
class FullyImplicitCompressiblePolymerSolver
|
||||
{
|
||||
public:
|
||||
/// Construct a solver. It will retain references to the
|
||||
/// arguments of this functions, and they are expected to
|
||||
/// remain in scope for the lifetime of the solver.
|
||||
/// \param[in] grid grid data structure
|
||||
/// \param[in] fluid fluid properties
|
||||
/// \param[in] geo rock properties
|
||||
/// \param[in] rock_comp_props if non-null, rock compressibility properties
|
||||
/// \param[in] wells well structure
|
||||
/// \param[in] linsolver linear solver
|
||||
FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const DerivedGeology& geo ,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
const PolymerPropsAd& polymer_props_ad,
|
||||
const Wells& wells,
|
||||
const LinearSolverInterface& linsolver);
|
||||
|
||||
/// Take a single forward step, modifiying
|
||||
/// state.pressure()
|
||||
/// state.faceflux()
|
||||
/// state.saturation()
|
||||
/// state.gasoilratio()
|
||||
/// wstate.bhp()
|
||||
/// \param[in] dt time step size
|
||||
/// \param[in] state reservoir state
|
||||
/// \param[in] wstate well state
|
||||
void
|
||||
step(const double dt ,
|
||||
PolymerBlackoilState& state ,
|
||||
WellState& wstate,
|
||||
const std::vector<double>& polymer_inflow);
|
||||
|
||||
private:
|
||||
// Types and enums
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef Eigen::Array<double,
|
||||
Eigen::Dynamic,
|
||||
Eigen::Dynamic,
|
||||
Eigen::RowMajor> DataBlock;
|
||||
|
||||
struct ReservoirResidualQuant {
|
||||
ReservoirResidualQuant();
|
||||
std::vector<ADB> accum; // Accumulations
|
||||
ADB mflux; // Mass flux (surface conditions)
|
||||
ADB b; // Reciprocal FVF
|
||||
ADB head; // Pressure drop across int. interfaces
|
||||
ADB mob; // Phase mobility (per cell)
|
||||
};
|
||||
|
||||
struct SolutionState {
|
||||
SolutionState(const int np);
|
||||
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)
|
||||
};
|
||||
|
||||
enum { Water = BlackoilPropsAdInterface::Water,
|
||||
Oil = BlackoilPropsAdInterface::Oil };
|
||||
|
||||
// Member data
|
||||
const UnstructuredGrid& grid_;
|
||||
const BlackoilPropsAdInterface& fluid_;
|
||||
const DerivedGeology& geo_;
|
||||
const RockCompressibility* rock_comp_props_;
|
||||
const PolymerPropsAd& polymer_props_ad_;
|
||||
const Wells& wells_;
|
||||
const LinearSolverInterface& linsolver_;
|
||||
const std::vector<int> cells_; // All grid cells
|
||||
HelperOps ops_;
|
||||
const WellOps wops_;
|
||||
const M grav_;
|
||||
|
||||
std::vector<ReservoirResidualQuant> rq_;
|
||||
|
||||
// The mass_balance vector has one element for each active phase,
|
||||
// each of which has size equal to the number of cells.
|
||||
// The well_eq has size equal to the number of wells.
|
||||
struct {
|
||||
std::vector<ADB> mass_balance;
|
||||
ADB well_flux_eq;
|
||||
ADB well_eq;
|
||||
} residual_;
|
||||
|
||||
// Private methods.
|
||||
SolutionState
|
||||
constantState(const PolymerBlackoilState& x,
|
||||
const WellState& xw);
|
||||
|
||||
SolutionState
|
||||
variableState(const PolymerBlackoilState& x,
|
||||
const WellState& xw);
|
||||
|
||||
void
|
||||
computeAccum(const SolutionState& state,
|
||||
const int aix );
|
||||
|
||||
void
|
||||
assemble(const V& pvdt,
|
||||
const PolymerBlackoilState& x,
|
||||
const WellState& xw,
|
||||
const std::vector<double>& polymer_inflow);
|
||||
|
||||
V solveJacobianSystem() const;
|
||||
|
||||
void updateState(const V& dx,
|
||||
PolymerBlackoilState& state,
|
||||
WellState& well_state) const;
|
||||
|
||||
std::vector<ADB>
|
||||
computeRelPerm(const SolutionState& state) const;
|
||||
|
||||
std::vector<ADB>
|
||||
computeRelPermWells(const SolutionState& state,
|
||||
const DataBlock& well_s,
|
||||
const std::vector<int>& well_cells) const;
|
||||
|
||||
void
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const std::vector<ADB>& kr ,
|
||||
const SolutionState& state );
|
||||
void
|
||||
computeMassFlux(const V& trans,
|
||||
const ADB& mc,
|
||||
const ADB& kro,
|
||||
const ADB& krw_eff,
|
||||
const SolutionState& state);
|
||||
|
||||
std::vector<ADB>
|
||||
computeFracFlow(const ADB& kro,
|
||||
const ADB& krw_eff,
|
||||
const ADB& c) const;
|
||||
ADB
|
||||
computeCmax(const ADB& c) const;
|
||||
ADB
|
||||
computeMc(const SolutionState& state) const;
|
||||
ADB
|
||||
rockPorosity(const ADB& p) const;
|
||||
ADB
|
||||
rockPermeability(const ADB& p) const;
|
||||
double
|
||||
residualNorm() const;
|
||||
|
||||
ADB
|
||||
fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const;
|
||||
|
||||
ADB
|
||||
fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const;
|
||||
|
||||
ADB
|
||||
fluidDensity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const;
|
||||
|
||||
ADB
|
||||
poroMult(const ADB& p) const;
|
||||
|
||||
ADB
|
||||
transMult(const ADB& p) const;
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
|
||||
@ -698,8 +698,6 @@ namespace {
|
||||
so -= sw;
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
state.saturation()[c*np] = sw[c];
|
||||
}
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
state.saturation()[c*np + 1] = so[c];
|
||||
}
|
||||
|
||||
@ -793,7 +791,7 @@ namespace {
|
||||
V trans(grid_.cell_facepos[nc]);
|
||||
UnstructuredGrid* ug = const_cast<UnstructuredGrid*>(& grid_);
|
||||
tpfa_htrans_compute(ug, fluid_.permeability(), htrans.data());
|
||||
tpfa_trans_compute (ug, htrans.data() , trans.data());
|
||||
tpfa_trans_compute (ug, htrans.data(), trans.data());
|
||||
|
||||
return trans;
|
||||
}
|
||||
|
||||
110
opm/polymer/fullyimplicit/GeoProps.hpp
Normal file
110
opm/polymer/fullyimplicit/GeoProps.hpp
Normal file
@ -0,0 +1,110 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_GEOPROPS_HEADER_INCLUDED
|
||||
#define OPM_GEOPROPS_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||
#include <Eigen/Eigen>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/// Class containing static geological properties that are
|
||||
/// derived from grid and petrophysical properties:
|
||||
/// - pore volume
|
||||
/// - transmissibilities
|
||||
/// - gravity potentials
|
||||
class DerivedGeology
|
||||
{
|
||||
public:
|
||||
typedef Eigen::ArrayXd Vector;
|
||||
|
||||
/// Construct contained derived geological properties
|
||||
/// from grid and property information.
|
||||
template <class Props>
|
||||
DerivedGeology(const UnstructuredGrid& grid,
|
||||
const Props& props ,
|
||||
const double* grav = 0)
|
||||
: pvol_ (grid.number_of_cells)
|
||||
, trans_(grid.number_of_faces)
|
||||
, gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
|
||||
, z_(grid.number_of_cells)
|
||||
{
|
||||
// Pore volume
|
||||
const typename Vector::Index nc = grid.number_of_cells;
|
||||
std::transform(grid.cell_volumes, grid.cell_volumes + nc,
|
||||
props.porosity(), pvol_.data(),
|
||||
std::multiplies<double>());
|
||||
|
||||
// Transmissibility
|
||||
Vector htrans(grid.cell_facepos[nc]);
|
||||
UnstructuredGrid* ug = const_cast<UnstructuredGrid*>(& grid);
|
||||
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
|
||||
tpfa_trans_compute (ug, htrans.data() , trans_.data());
|
||||
|
||||
// Compute z coordinates
|
||||
for (int c = 0; c<nc; ++c){
|
||||
z_[c] = grid.cell_centroids[c*3 + 2];
|
||||
}
|
||||
|
||||
|
||||
// Gravity potential
|
||||
std::fill(gravity_, gravity_ + 3, 0.0);
|
||||
if (grav != 0) {
|
||||
const typename Vector::Index nd = grid.dimensions;
|
||||
|
||||
for (typename Vector::Index c = 0; c < nc; ++c) {
|
||||
const double* const cc = & grid.cell_centroids[c*nd + 0];
|
||||
|
||||
const int* const p = grid.cell_facepos;
|
||||
for (int i = p[c]; i < p[c + 1]; ++i) {
|
||||
const int f = grid.cell_faces[i];
|
||||
|
||||
const double* const fc = & grid.face_centroids[f*nd + 0];
|
||||
|
||||
for (typename Vector::Index d = 0; d < nd; ++d) {
|
||||
gpot_[i] += grav[d] * (fc[d] - cc[d]);
|
||||
}
|
||||
}
|
||||
}
|
||||
std::copy(grav, grav + nd, gravity_);
|
||||
}
|
||||
}
|
||||
|
||||
const Vector& poreVolume() const { return pvol_ ;}
|
||||
const Vector& transmissibility() const { return trans_ ;}
|
||||
const Vector& gravityPotential() const { return gpot_ ;}
|
||||
const Vector& z() const { return z_ ;}
|
||||
const double* gravity() const { return gravity_;}
|
||||
|
||||
private:
|
||||
Vector pvol_ ;
|
||||
Vector trans_;
|
||||
Vector gpot_ ;
|
||||
Vector z_;
|
||||
double gravity_[3]; // Size 3 even if grid is 2-dim.
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
#endif // OPM_GEOPROPS_HEADER_INCLUDED
|
||||
@ -0,0 +1,522 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
|
||||
#if HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif // HAVE_CONFIG_H
|
||||
|
||||
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <opm/polymer/fullyimplicit/GeoProps.hpp>
|
||||
#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
|
||||
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/wells.h>
|
||||
#include <opm/core/pressure/flow_bc.h>
|
||||
|
||||
#include <opm/core/simulator/SimulatorReport.hpp>
|
||||
#include <opm/core/simulator/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/StopWatch.hpp>
|
||||
#include <opm/core/io/eclipse/EclipseWriter.hpp>
|
||||
#include <opm/core/io/vtk/writeVtkData.hpp>
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
|
||||
|
||||
#include <opm/core/wells/WellsManager.hpp>
|
||||
|
||||
#include <opm/core/props/rock/RockCompressibility.hpp>
|
||||
|
||||
#include <opm/core/grid/ColumnExtract.hpp>
|
||||
#include <opm/polymer/PolymerBlackoilState.hpp>
|
||||
#include <opm/polymer/PolymerInflow.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
|
||||
|
||||
#include <boost/filesystem.hpp>
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
#include <boost/lexical_cast.hpp>
|
||||
|
||||
#include <numeric>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
class SimulatorFullyImplicitCompressiblePolymer::Impl
|
||||
{
|
||||
public:
|
||||
Impl(const parameter::ParameterGroup& param,
|
||||
const UnstructuredGrid& grid,
|
||||
const BlackoilPropsAdInterface& props,
|
||||
const PolymerPropsAd& polymer_props,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
WellsManager& wells_manager,
|
||||
PolymerInflowInterface& polymer_inflow,
|
||||
LinearSolverInterface& linsolver,
|
||||
const double* gravity);
|
||||
|
||||
SimulatorReport run(SimulatorTimer& timer,
|
||||
PolymerBlackoilState& state,
|
||||
WellState& well_state);
|
||||
|
||||
private:
|
||||
// Data.
|
||||
|
||||
// Parameters for output.
|
||||
bool output_;
|
||||
bool output_vtk_;
|
||||
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 BlackoilPropsAdInterface& props_;
|
||||
const PolymerPropsAd& polymer_props_;
|
||||
const RockCompressibility* rock_comp_props_;
|
||||
WellsManager& wells_manager_;
|
||||
const Wells* wells_;
|
||||
PolymerInflowInterface& polymer_inflow_;
|
||||
const double* gravity_;
|
||||
// Solvers
|
||||
DerivedGeology geo_;
|
||||
FullyImplicitCompressiblePolymerSolver solver_;
|
||||
// Misc. data
|
||||
std::vector<int> allcells_;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
SimulatorFullyImplicitCompressiblePolymer::
|
||||
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
|
||||
const UnstructuredGrid& grid,
|
||||
const BlackoilPropsAdInterface& props,
|
||||
const PolymerPropsAd& polymer_props,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
WellsManager& wells_manager,
|
||||
PolymerInflowInterface& polymer_inflow,
|
||||
LinearSolverInterface& linsolver,
|
||||
const double* gravity)
|
||||
|
||||
{
|
||||
pimpl_.reset(new Impl(param, grid, props, polymer_props, rock_comp_props, wells_manager, polymer_inflow, linsolver, gravity));
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer,
|
||||
PolymerBlackoilState& state,
|
||||
WellState& well_state)
|
||||
{
|
||||
return pimpl_->run(timer, state, well_state);
|
||||
}
|
||||
|
||||
|
||||
|
||||
static void outputStateVtk(const UnstructuredGrid& grid,
|
||||
const Opm::PolymerBlackoilState& state,
|
||||
const int step,
|
||||
const std::string& output_dir)
|
||||
{
|
||||
// Write data in VTK format.
|
||||
std::ostringstream vtkfilename;
|
||||
vtkfilename << output_dir << "/vtk_files";
|
||||
boost::filesystem::path fpath(vtkfilename.str());
|
||||
try {
|
||||
create_directories(fpath);
|
||||
}
|
||||
catch (...) {
|
||||
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
|
||||
}
|
||||
vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
|
||||
std::ofstream vtkfile(vtkfilename.str().c_str());
|
||||
if (!vtkfile) {
|
||||
OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
|
||||
}
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
std::vector<double> cell_velocity;
|
||||
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
|
||||
dm["velocity"] = &cell_velocity;
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
}
|
||||
|
||||
|
||||
static void outputStateMatlab(const UnstructuredGrid& grid,
|
||||
const Opm::PolymerBlackoilState& state,
|
||||
const int step,
|
||||
const std::string& output_dir)
|
||||
{
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
dm["surfvolume"] = &state.surfacevol();
|
||||
std::vector<double> cell_velocity;
|
||||
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
|
||||
dm["velocity"] = &cell_velocity;
|
||||
|
||||
// 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"));
|
||||
}
|
||||
}
|
||||
static void outputWellStateMatlab(const Opm::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"));
|
||||
}
|
||||
}
|
||||
|
||||
#if 0
|
||||
static void outputWaterCut(const Opm::Watercut& watercut,
|
||||
const std::string& output_dir)
|
||||
{
|
||||
// Write water cut curve.
|
||||
std::string fname = output_dir + "/watercut.txt";
|
||||
std::ofstream os(fname.c_str());
|
||||
if (!os) {
|
||||
OPM_THROW(std::runtime_error, "Failed to open " << fname);
|
||||
}
|
||||
watercut.write(os);
|
||||
}
|
||||
|
||||
static void outputWellReport(const Opm::WellReport& wellreport,
|
||||
const std::string& output_dir)
|
||||
{
|
||||
// Write well report.
|
||||
std::string fname = output_dir + "/wellreport.txt";
|
||||
std::ofstream os(fname.c_str());
|
||||
if (!os) {
|
||||
OPM_THROW(std::runtime_error, "Failed to open " << fname);
|
||||
}
|
||||
wellreport.write(os);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
// \TODO: Treat bcs.
|
||||
SimulatorFullyImplicitCompressiblePolymer::Impl::Impl(const parameter::ParameterGroup& param,
|
||||
const UnstructuredGrid& grid,
|
||||
const BlackoilPropsAdInterface& props,
|
||||
const PolymerPropsAd& polymer_props,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
WellsManager& wells_manager,
|
||||
PolymerInflowInterface& polymer_inflow,
|
||||
LinearSolverInterface& linsolver,
|
||||
const double* gravity)
|
||||
: grid_(grid),
|
||||
props_(props),
|
||||
polymer_props_(polymer_props),
|
||||
rock_comp_props_(rock_comp_props),
|
||||
wells_manager_(wells_manager),
|
||||
wells_(wells_manager.c_wells()),
|
||||
polymer_inflow_(polymer_inflow),
|
||||
gravity_(gravity),
|
||||
geo_(grid_, props_, gravity_),
|
||||
solver_(grid_, props_, geo_, rock_comp_props, polymer_props, *wells_manager.c_wells(), linsolver)
|
||||
|
||||
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
|
||||
param.getDefault("nl_pressure_change_tolerance", 1.0),
|
||||
param.getDefault("nl_pressure_maxiter", 10),
|
||||
gravity, */
|
||||
{
|
||||
// For output.
|
||||
output_ = param.getDefault("output", true);
|
||||
if (output_) {
|
||||
output_vtk_ = param.getDefault("output_vtk", true);
|
||||
output_dir_ = param.getDefault("output_dir", std::string("output"));
|
||||
// Ensure that output dir exists
|
||||
boost::filesystem::path fpath(output_dir_);
|
||||
try {
|
||||
create_directories(fpath);
|
||||
}
|
||||
catch (...) {
|
||||
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
|
||||
}
|
||||
output_interval_ = param.getDefault("output_interval", 1);
|
||||
}
|
||||
|
||||
// Well control related init.
|
||||
check_well_controls_ = param.getDefault("check_well_controls", false);
|
||||
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
|
||||
|
||||
// Misc init.
|
||||
const int num_cells = grid.number_of_cells;
|
||||
allcells_.resize(num_cells);
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
allcells_[cell] = cell;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer,
|
||||
PolymerBlackoilState& state,
|
||||
WellState& well_state)
|
||||
{
|
||||
|
||||
// Initialisation.
|
||||
std::vector<double> porevol;
|
||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
||||
} 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);
|
||||
// Main simulation loop.
|
||||
Opm::time::StopWatch solver_timer;
|
||||
double stime = 0.0;
|
||||
Opm::time::StopWatch step_timer;
|
||||
Opm::time::StopWatch total_timer;
|
||||
total_timer.start();
|
||||
#if 0
|
||||
// These must be changed for three-phase.
|
||||
double init_surfvol[2] = { 0.0 };
|
||||
double inplace_surfvol[2] = { 0.0 };
|
||||
double tot_injected[2] = { 0.0 };
|
||||
double tot_produced[2] = { 0.0 };
|
||||
Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
|
||||
Opm::Watercut watercut;
|
||||
watercut.push(0.0, 0.0, 0.0);
|
||||
Opm::WellReport wellreport;
|
||||
#endif
|
||||
std::vector<double> fractional_flows;
|
||||
std::vector<double> well_resflows_phase;
|
||||
if (wells_) {
|
||||
well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
|
||||
#if 0
|
||||
wellreport.push(props_, *wells_,
|
||||
state.pressure(), state.surfacevol(), state.saturation(),
|
||||
0.0, well_state.bhp(), well_state.perfRates());
|
||||
#endif
|
||||
}
|
||||
std::fstream tstep_os;
|
||||
if (output_) {
|
||||
std::string filename = output_dir_ + "/step_timing.param";
|
||||
tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
|
||||
}
|
||||
while (!timer.done()) {
|
||||
// Report timestep and (optionally) write state to disk.
|
||||
step_timer.start();
|
||||
timer.report(std::cout);
|
||||
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
||||
if (output_vtk_) {
|
||||
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
||||
}
|
||||
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
||||
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
|
||||
|
||||
}
|
||||
|
||||
SimulatorReport sreport;
|
||||
|
||||
// Solve pressure equation.
|
||||
// if (check_well_controls_) {
|
||||
// computeFractionalFlow(props_, allcells_,
|
||||
// state.pressure(), state.surfacevol(), state.saturation(),
|
||||
// fractional_flows);
|
||||
// wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
|
||||
// }
|
||||
bool well_control_passed = !check_well_controls_;
|
||||
int well_control_iteration = 0;
|
||||
do {
|
||||
// Run solver.
|
||||
const double current_time = timer.currentTime();
|
||||
double stepsize = timer.currentStepLength();
|
||||
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
||||
solver_timer.start();
|
||||
std::vector<double> initial_pressure = state.pressure();
|
||||
solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
|
||||
|
||||
// Stop timer and report.
|
||||
solver_timer.stop();
|
||||
const double st = solver_timer.secsSinceStart();
|
||||
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
|
||||
|
||||
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.
|
||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||
initial_porevol = porevol;
|
||||
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
||||
}
|
||||
|
||||
// The reports below are geared towards two phases only.
|
||||
#if 0
|
||||
// Report mass balances.
|
||||
double injected[2] = { 0.0 };
|
||||
double produced[2] = { 0.0 };
|
||||
Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
|
||||
injected, produced);
|
||||
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
|
||||
tot_injected[0] += injected[0];
|
||||
tot_injected[1] += injected[1];
|
||||
tot_produced[0] += produced[0];
|
||||
tot_produced[1] += produced[1];
|
||||
std::cout.precision(5);
|
||||
const int width = 18;
|
||||
std::cout << "\nMass balance report.\n";
|
||||
std::cout << " Injected surface volumes: "
|
||||
<< std::setw(width) << injected[0]
|
||||
<< std::setw(width) << injected[1] << std::endl;
|
||||
std::cout << " Produced surface volumes: "
|
||||
<< std::setw(width) << produced[0]
|
||||
<< std::setw(width) << produced[1] << std::endl;
|
||||
std::cout << " Total inj surface volumes: "
|
||||
<< std::setw(width) << tot_injected[0]
|
||||
<< std::setw(width) << tot_injected[1] << std::endl;
|
||||
std::cout << " Total prod surface volumes: "
|
||||
<< std::setw(width) << tot_produced[0]
|
||||
<< std::setw(width) << tot_produced[1] << std::endl;
|
||||
const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
|
||||
init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] };
|
||||
std::cout << " Initial - inplace + inj - prod: "
|
||||
<< std::setw(width) << balance[0]
|
||||
<< std::setw(width) << balance[1]
|
||||
<< std::endl;
|
||||
std::cout << " Relative mass error: "
|
||||
<< std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
|
||||
<< std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
|
||||
<< std::endl;
|
||||
std::cout.precision(8);
|
||||
|
||||
// Make well reports.
|
||||
watercut.push(timer.currentTime() + timer.currentStepLength(),
|
||||
produced[0]/(produced[0] + produced[1]),
|
||||
tot_produced[0]/tot_porevol_init);
|
||||
if (wells_) {
|
||||
wellreport.push(props_, *wells_,
|
||||
state.pressure(), state.surfacevol(), state.saturation(),
|
||||
timer.currentTime() + timer.currentStepLength(),
|
||||
well_state.bhp(), well_state.perfRates());
|
||||
}
|
||||
#endif
|
||||
sreport.total_time = step_timer.secsSinceStart();
|
||||
if (output_) {
|
||||
sreport.reportParam(tstep_os);
|
||||
|
||||
if (output_vtk_) {
|
||||
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_) {
|
||||
outputWellReport(wellreport, output_dir_);
|
||||
}
|
||||
#endif
|
||||
tstep_os.close();
|
||||
}
|
||||
|
||||
// advance to next timestep before reporting at this location
|
||||
++timer;
|
||||
|
||||
// write an output file for later inspection
|
||||
}
|
||||
|
||||
total_timer.stop();
|
||||
|
||||
SimulatorReport report;
|
||||
report.pressure_time = stime;
|
||||
report.transport_time = 0.0;
|
||||
report.total_time = total_timer.secsSinceStart();
|
||||
return report;
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@ -0,0 +1,98 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
||||
#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
||||
|
||||
#include <boost/shared_ptr.hpp>
|
||||
#include <vector>
|
||||
|
||||
struct UnstructuredGrid;
|
||||
struct Wells;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace parameter { class ParameterGroup; }
|
||||
class BlackoilPropsAdInterface;
|
||||
class RockCompressibility;
|
||||
class WellsManager;
|
||||
class LinearSolverInterface;
|
||||
class SimulatorTimer;
|
||||
class PolymerBlackoilState;
|
||||
class WellState;
|
||||
class PolymerPropsAd;
|
||||
class PolymerInflowInterface;
|
||||
struct SimulatorReport;
|
||||
|
||||
/// Class collecting all necessary components for a two-phase simulation.
|
||||
class SimulatorFullyImplicitCompressiblePolymer
|
||||
{
|
||||
public:
|
||||
/// Initialise from parameters and objects to observe.
|
||||
/// \param[in] param parameters, this class accepts the following:
|
||||
/// parameter (default) effect
|
||||
/// -----------------------------------------------------------
|
||||
/// output (true) write output to files?
|
||||
/// output_dir ("output") output directoty
|
||||
/// output_interval (1) output every nth step
|
||||
/// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
|
||||
/// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal)
|
||||
/// nl_pressure_maxiter (10) max nonlinear iterations in pressure
|
||||
/// nl_maxiter (30) max nonlinear iterations in transport
|
||||
/// nl_tolerance (1e-9) transport solver absolute residual tolerance
|
||||
/// num_transport_substeps (1) number of transport steps per pressure step
|
||||
/// use_segregation_split (false) solve for gravity segregation (if false,
|
||||
/// segregation is ignored).
|
||||
///
|
||||
/// \param[in] grid grid data structure
|
||||
/// \param[in] props fluid and rock properties
|
||||
/// \param[in] rock_comp_props if non-null, rock compressibility properties
|
||||
/// \param[in] well_manager well manager, may manage no (null) wells
|
||||
/// \param[in] linsolver linear solver
|
||||
/// \param[in] gravity if non-null, gravity vector
|
||||
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
|
||||
const UnstructuredGrid& grid,
|
||||
const BlackoilPropsAdInterface& props,
|
||||
const PolymerPropsAd& polymer_props,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
WellsManager& wells_manager,
|
||||
PolymerInflowInterface& polymer_inflow,
|
||||
LinearSolverInterface& linsolver,
|
||||
const double* gravity);
|
||||
|
||||
/// Run the simulation.
|
||||
/// This will run succesive timesteps until timer.done() is true. It will
|
||||
/// modify the reservoir and well states.
|
||||
/// \param[in,out] timer governs the requested reporting timesteps
|
||||
/// \param[in,out] state state of reservoir: pressure, fluxes
|
||||
/// \param[in,out] well_state state of wells: bhp, perforation rates
|
||||
/// \return simulation report, with timing data
|
||||
SimulatorReport run(SimulatorTimer& timer,
|
||||
PolymerBlackoilState& state,
|
||||
WellState& well_state);
|
||||
|
||||
private:
|
||||
class Impl;
|
||||
// Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
|
||||
boost::shared_ptr<Impl> pimpl_;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
||||
@ -146,6 +146,7 @@ namespace Opm
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
std::vector<double> cell_velocity;
|
||||
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
|
||||
dm["velocity"] = &cell_velocity;
|
||||
@ -161,6 +162,7 @@ namespace Opm
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
std::vector<double> cell_velocity;
|
||||
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
|
||||
dm["velocity"] = &cell_velocity;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user