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

Remove unused classes.

Browse Source 
After this, the two affected tests fail due to bugs in PVT region
support in BlackoilPropsAdFromDeck.
This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-23 09:40:00 +01:00
parent eb89236552
commit d9ce8625cf
9 changed files with 38 additions and 2227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists_files.cmake
View File

@ -26,14 +26,12 @@
# originally generated with the command:
# find opm -name '*.c*' -printf '\t%p\n' | sort
list (APPEND MAIN_SOURCE_FILES
opm/autodiff/BlackoilPropsAd.cpp
opm/autodiff/BlackoilPropsAdInterface.cpp
opm/autodiff/ExtractParallelGridInformationToISTL.cpp
opm/autodiff/NewtonIterationBlackoilCPR.cpp
opm/autodiff/NewtonIterationBlackoilSimple.cpp
opm/autodiff/GridHelpers.cpp
opm/autodiff/ImpesTPFAAD.cpp
opm/autodiff/SimulatorCompressibleAd.cpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
opm/autodiff/SimulatorIncompTwophaseAd.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp
@ -75,7 +73,6 @@ endif()
list (APPEND EXAMPLE_SOURCE_FILES
examples/find_zero.cpp
examples/sim_fibo_ad.cpp
examples/sim_2p_comp_ad.cpp
examples/sim_2p_incomp_ad.cpp
examples/sim_simple.cpp
examples/opm_init_check.cpp
@ -96,7 +93,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/AutoDiffHelpers.hpp
opm/autodiff/AutoDiff.hpp
opm/autodiff/BackupRestore.hpp
opm/autodiff/BlackoilPropsAd.hpp
opm/autodiff/BlackoilPropsAdFromDeck.hpp
opm/autodiff/BlackoilPropsAdInterface.hpp
opm/autodiff/CPRPreconditioner.hpp
@ -113,8 +109,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/NewtonIterationBlackoilSimple.hpp
opm/autodiff/LinearisedBlackoilResidual.hpp
opm/autodiff/RateConverter.hpp
opm/autodiff/RedistributeDataHandles.hpp
opm/autodiff/SimulatorCompressibleAd.hpp
opm/autodiff/RedistributeDataHandles.hpp
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
opm/autodiff/SimulatorIncompTwophaseAd.hpp

211
examples/sim_2p_comp_ad.cpp
View File

@ -1,211 +0,0 @@
/*
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/well_controls.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/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/SimulatorCompressibleAd.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <memory>
#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 weakly compressible two-phase flow ===============\n\n";
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
Opm::DeckConstPtr deck;
std::unique_ptr<GridManager> grid;
std::unique_ptr<BlackoilPropertiesInterface> props;
std::unique_ptr<RockCompressibility> rock_comp;
EclipseStateConstPtr eclipseState;
BlackoilState 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");
Opm::ParserPtr parser(new Opm::Parser());
deck = parser->parseFile( deck_filename );
eclipseState.reset(new EclipseState(deck));
// Grid init
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
// 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, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
initBlackoilSurfvol(*grid->c_grid(), *props, 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";
SimulatorReport rep;
// With a deck, we may have more epochs etc.
WellState well_state;
Opm::TimeMapConstPtr timeMap = eclipseState->getSchedule()->getTimeMap();
Opm::DerivedGeology geology(*grid->c_grid(), *props, eclipseState,false);
SimulatorTimer simtimer;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of report step.
std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------"
<< "\n (number of steps left: "
<< timeMap->numTimesteps() - reportStepIdx << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(eclipseState, reportStepIdx, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every report step,
// since number of wells may change etc.
if (reportStepIdx == 0) {
well_state.init(wells.c_wells(), state);
}
simtimer.setCurrentStepNum(reportStepIdx);
// Create and run simulator.
SimulatorCompressibleAd simulator(param,
*grid->c_grid(),
geology,
*props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
linsolver,
grav);
if (reportStepIdx == 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;
}
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;
}

940
opm/autodiff/BlackoilPropsAd.cpp
View File

@ -1,940 +0,0 @@
/*
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/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/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(int regionIdx) const
{
// this class only supports a single PVT region for now...
static_cast<void>(regionIdx);
assert(regionIdx == 0);
return props_.surfaceDensity();
}
// ------ Viscosity ------
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
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(), T.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] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& /*cond*/,
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(), T.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] T Array of n temperature 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 V& T,
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(), T.data(), z.data(), cells.data(), mu.data(), 0);
return mu.col(pu_.phase_pos[Gas]);
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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::muGas(const V& pg,
const V& T,
const V& rv,
const std::vector<PhasePresence>& /*cond*/,
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);
if (pu_.phase_used[Oil]) {
// Faking a z with the right ratio:
// rv = zo/zg
z.col(pu_.phase_pos[Oil]) = rv;
z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
}
Block mu(n, np);
props_.viscosity(n, pg.data(), T.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] T Array of n temperature 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 ADB& T,
const Cells& cells) const
{
#if 1
return ADB::constant(muWat(pw.value(), T.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(), T.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] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
{
#if 1
return ADB::constant(muOil(po.value(), T.value(), rs.value(), cond, 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] T Array of n temperature 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 ADB& T,
const Cells& cells) const
{
#if 1
return ADB::constant(muGas(pg.value(), T.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
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
{
#if 1
return ADB::constant(muGas(pg.value(), T.value(), rv.value(),cond,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);
if (pu_.phase_used[Oil]) {
// Faking a z with the right ratio:
// rv = zo/zg
z.col(pu_.phase_pos[Oil]) = rv;
z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
}
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] T Array of n temperature 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 V& T,
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(), T.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] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& /*cond*/,
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(), T.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] T Array of n temperature 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 V& /* T */,
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(), pg.data(), z.data(), cells.data(), matrix.data(), 0);
const int gi = pu_.phase_pos[Gas];
return matrix.col(gi*np + gi);
}
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& /*cond*/,
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);
if (pu_.phase_used[Oil]) {
// Faking a z with the right ratio:
// rv = zo/zg
z.col(pu_.phase_pos[Oil]) = rv;
z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
}
Block matrix(n, np*np);
props_.matrix(n, pg.data(), T.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] T Array of n temperature 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 ADB& T,
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(), T.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] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& /*cond*/,
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(), T.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] T Array of n temperature 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 ADB& T,
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(), T.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);
}
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& /*cond*/,
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);
if (pu_.phase_used[Oil]) {
// Faking a z with the right ratio:
// rv = zo/zg
z.col(pu_.phase_pos[Oil]) = rv.value();
z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
}
Block matrix(n, np*np);
Block dmatrix(n, np*np);
props_.matrix(n, pg.value().data(), T.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::rsSat(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] so Array of n oil saturation 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::rsSat(const V& po,
const V& so,
const Cells& cells) const
{
// Suppress warning about "unused parameters".
static_cast<void>(po);
static_cast<void>(so);
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::rsSat(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.");
}
/// Bubble point curve for Rs as function of oil pressure.
/// \param[in] po Array of n oil pressure values.
/// \param[in] so Array of n oil saturation 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::rsSat(const ADB& po,
const ADB& so,
const Cells& cells) const
{
// Suppress warning about "unused parameters".
static_cast<void>(po);
static_cast<void>(so);
static_cast<void>(cells);
OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
}
// ------ 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::rvSat(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] so Array of n oil saturation 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::rvSat(const V& po,
const V& so,
const Cells& cells) const
{
// Suppress warning about "unused parameters".
static_cast<void>(po);
static_cast<void>(so);
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::rvSat(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.");
}
/// Bubble point curve for Rs as function of oil pressure.
/// \param[in] po Array of n oil pressure values.
/// \param[in] so Array of n oil saturation 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::rvSat(const ADB& po,
const ADB& so,
const Cells& cells) const
{
// Suppress warning about "unused parameters".
static_cast<void>(po);
static_cast<void>(so);
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;
}
std::vector<ADB> BlackoilPropsAd::capPress(const ADB& sw,
const ADB& so,
const ADB& sg,
const Cells& cells) const
{
const int numCells = cells.size();
const int numActivePhases = numPhases();
const int numBlocks = so.numBlocks();
Block activeSat(numCells, numActivePhases);
if (pu_.phase_used[Water]) {
assert(sw.value().size() == numCells);
activeSat.col(pu_.phase_pos[Water]) = sw.value();
}
if (pu_.phase_used[Oil]) {
assert(so.value().size() == numCells);
activeSat.col(pu_.phase_pos[Oil]) = so.value();
} else {
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
}
if (pu_.phase_used[Gas]) {
assert(sg.value().size() == numCells);
activeSat.col(pu_.phase_pos[Gas]) = sg.value();
}
Block pc(numCells, numActivePhases);
Block dpc(numCells, numActivePhases*numActivePhases);
props_.capPress(numCells, activeSat.data(), cells.data(), pc.data(), dpc.data());
std::vector<ADB> adbCapPressures;
adbCapPressures.reserve(3);
const ADB* 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(numBlocks);
for (int block = 0; block < numBlocks; ++block) {
jacs[block] = ADB::M(numCells, 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 dpc1/ds2.
const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
for (int block = 0; block < numBlocks; ++block) {
jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
}
}
adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
} else {
adbCapPressures.emplace_back(ADB::null());
}
}
return adbCapPressures;
}
/// Saturation update for hysteresis behavior.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
void BlackoilPropsAd::updateSatHyst(const std::vector<double>& /* saturation */,
const std::vector<int>& /* cells */)
{
OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support hysteresis.");
}
/// Update for max oil saturation.
void BlackoilPropsAd::updateSatOilMax(const std::vector<double>& /*saturation*/)
{
OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support this functionality.");
}
} // namespace Opm

397
opm/autodiff/BlackoilPropsAd.hpp
View File

@ -1,397 +0,0 @@
/*
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/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/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;
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 3 density values.
const double* surfaceDensity(int regionIdx = 0) const;
// ------ Viscosity ------
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
// ------ Formation volume factor (b) ------
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
// ------ Rs bubble point curve ------
/// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
/// \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 rsSat(const V& po,
const Cells& cells) const;
/// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
/// \param[in] po Array of n oil pressure values.
/// \param[in] so Array of n oil saturation 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 rsSat(const V& po,
const V& so,
const Cells& cells) const;
/// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
/// \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 rsSat(const ADB& po,
const Cells& cells) const;
/// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
/// \param[in] po Array of n oil pressure values.
/// \param[in] so Array of n oil saturation 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 rsSat(const ADB& po,
const ADB& so,
const Cells& cells) const;
// ------ Rv condensation curve ------
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
/// \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 rvSat(const V& po,
const Cells& cells) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
/// \param[in] po Array of n oil pressure values.
/// \param[in] so Array of n oil saturation 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 rvSat(const V& po,
const V& so,
const Cells& cells) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
/// \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 rvSat(const ADB& po,
const Cells& cells) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
/// \param[in] po Array of n oil pressure values.
/// \param[in] so Array of n oil saturation 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 rvSat(const ADB& po,
const ADB& so,
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;
/// Capillary pressure 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 capillary pressure values,
/// containing the offsets for each p_g, p_o, p_w. The capillary pressure between
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
std::vector<ADB> capPress(const ADB& sw,
const ADB& so,
const ADB& sg,
const Cells& cells) const;
/// Saturation update for hysteresis behavior.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
void updateSatHyst(const std::vector<double>& saturation,
const std::vector<int>& cells);
/// Update for max oil saturation.
void updateSatOilMax(const std::vector<double>& saturation);
private:
const BlackoilPropertiesInterface& props_;
PhaseUsage pu_;
};
} // namespace Opm
#endif // OPM_BLACKOILPROPSAD_HEADER_INCLUDED

539
opm/autodiff/SimulatorCompressibleAd.cpp
View File

@ -1,539 +0,0 @@
/*
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/autodiff/SimulatorCompressibleAd.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/ImpesTPFAAD.hpp>
#include <opm/autodiff/BlackoilPropsAd.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/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/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/grid/ColumnExtract.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
namespace Opm
{
class SimulatorCompressibleAd::Impl
{
public:
Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity);
SimulatorReport run(SimulatorTimer& timer,
BlackoilState& 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_;
// Parameters for transport solver.
int num_transport_substeps_;
bool use_segregation_split_;
// Observed objects.
const UnstructuredGrid& grid_;
const BlackoilPropertiesInterface& props_;
const RockCompressibility* rock_comp_props_;
WellsManager& wells_manager_;
const Wells* wells_;
const double* gravity_;
// Solvers
BlackoilPropsAd fluid_;
const DerivedGeology& geo_;
ImpesTPFAAD psolver_;
TransportSolverCompressibleTwophaseReorder tsolver_;
// Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_;
// Misc. data
std::vector<int> allcells_;
};
SimulatorCompressibleAd::SimulatorCompressibleAd(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, wells_manager, linsolver, gravity));
}
SimulatorReport SimulatorCompressibleAd::run(SimulatorTimer& timer,
BlackoilState& state,
WellState& well_state)
{
return pimpl_->run(timer, state, well_state);
}
static void outputStateVtk(const UnstructuredGrid& grid,
const Opm::BlackoilState& 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();
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::BlackoilState& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
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 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);
}
// \TODO: make CompressibleTpfa take bcs.
SimulatorCompressibleAd::Impl::Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity)
: grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
wells_manager_(wells_manager),
wells_(wells_manager.c_wells()),
gravity_(gravity),
fluid_(props_),
geo_(geo),
psolver_(grid_, fluid_, geo_, *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, */
tsolver_(grid, props,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30))
{
// 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);
// Transport related init.
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity != 0 && use_segregation_split_){
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}
// 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 SimulatorCompressibleAd::Impl::run(SimulatorTimer& timer,
BlackoilState& state,
WellState& well_state)
{
std::vector<double> transport_src;
// 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;
// Main simulation loop.
Opm::time::StopWatch pressure_timer;
double ptime = 0.0;
Opm::time::StopWatch transport_timer;
double ttime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
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;
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);
wellreport.push(props_, *wells_,
state.pressure(), state.surfacevol(), state.saturation(),
0.0, well_state.bhp(), well_state.perfRates());
}
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);
}
for (; !timer.done(); ++timer) {
// 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_);
}
SimulatorReport sreport;
// Solve pressure equation.
if (check_well_controls_) {
computeFractionalFlow(props_, allcells_,
state.pressure(), state.temperature(), 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.
pressure_timer.start();
std::vector<double> initial_pressure = state.pressure();
psolver_.solve(timer.currentStepLength(), state, well_state);
#if 0
// Renormalize pressure if both fluids and rock are
// incompressible, and there are no pressure
// conditions (bcs or wells). It is deemed sufficient
// for now to renormalize using geometric volume
// instead of pore volume.
if (psolver_.singularPressure()) {
// Compute average pressures of previous and last
// step, and total volume.
double av_prev_press = 0.0;
double av_press = 0.0;
double tot_vol = 0.0;
const int num_cells = grid_.number_of_cells;
for (int cell = 0; cell < num_cells; ++cell) {
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
tot_vol += grid_.cell_volumes[cell];
}
// Renormalization constant
const double ren_const = (av_prev_press - av_press)/tot_vol;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += ren_const;
}
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
for (int well = 0; well < num_wells; ++well) {
well_state.bhp()[well] += ren_const;
}
}
#endif
// Stop timer and report.
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
sreport.pressure_time = pt;
// 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);
}
// Process transport sources from well flows.
Opm::computeTransportSource(props_, wells_, well_state, transport_src);
// Solve transport.
transport_timer.start();
double stepsize = timer.currentStepLength();
if (num_transport_substeps_ != 1) {
stepsize /= double(num_transport_substeps_);
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
}
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], &state.temperature()[0],
&initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
state.saturation(), state.surfacevol());
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
substep_injected, substep_produced);
injected[0] += substep_injected[0];
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol());
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
sreport.transport_time = tt;
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Report volume balances.
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);
watercut.push(timer.simulationTimeElapsed() + 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.simulationTimeElapsed() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
sreport.total_time = step_timer.secsSinceStart();
if (output_) {
sreport.reportParam(tstep_os);
}
}
if (output_) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWaterCut(watercut, output_dir_);
if (wells_) {
outputWellReport(wellreport, output_dir_);
}
tstep_os.close();
}
total_timer.stop();
SimulatorReport report;
report.pressure_time = ptime;
report.transport_time = ttime;
report.total_time = total_timer.secsSinceStart();
return report;
}
} // namespace Opm

98
opm/autodiff/SimulatorCompressibleAd.hpp
View File

@ -1,98 +0,0 @@
/*
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_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED
#define OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED
#include <memory>
#include <vector>
struct UnstructuredGrid;
struct Wells;
struct FlowBoundaryConditions;
namespace Opm
{
namespace parameter { class ParameterGroup; }
class BlackoilPropertiesInterface;
class DerivedGeology;
class RockCompressibility;
class WellsManager;
class LinearSolverInterface;
class SimulatorTimer;
class BlackoilState;
class WellState;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation.
class SimulatorCompressibleAd
{
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] geo the "ready to use" geological properties of the reservoir
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] well_manager well manager
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
SimulatorCompressibleAd(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
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,
BlackoilState& state,
WellState& well_state);
private:
class Impl;
// Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
std::shared_ptr<Impl> pimpl_;
};
} // namespace Opm
#endif // OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED

11
tests/fluid.data
View File

@ -3,6 +3,7 @@ RUNSPEC
OIL
WATER
GAS
METRIC
@ -47,11 +48,21 @@ PVCDO
1 1 0 1000 0
/
PVDG
-- Pg Bg(Pg) mug
1 1 1
/
SWOF
0 0 1 0
1 1 0 0
/
SGOF
0 0 1 0
1 1 0 0
/
DENSITY
800 1000 1
/

47
tests/test_boprops_ad.cpp
View File

@ -26,13 +26,11 @@
#define BOOST_TEST_MODULE FluidPropertiesTest
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <boost/test/unit_test.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
@ -70,8 +68,7 @@ struct TestFixture : public Setup
TestFixture()
: Setup()
, grid (deck)
, props(deck, eclState, *grid.c_grid(), param,
param.getDefault("init_rock", false))
, boprops_ad(deck, eclState, *grid.c_grid(), param.getDefault("init_rock", false))
{
}
@ -79,8 +76,8 @@ struct TestFixture : public Setup
using Setup::deck;
using Setup::eclState;
Opm::GridManager grid;
Opm::BlackoilPropertiesFromDeck props;
Opm::GridManager grid;
Opm::BlackoilPropsAdFromDeck boprops_ad;
};
template <class Setup>
@ -105,7 +102,6 @@ struct TestFixtureAd : public Setup
BOOST_FIXTURE_TEST_CASE(Construction, TestFixture<SetupSimple>)
{
Opm::BlackoilPropsAd boprops_ad(props);
}
BOOST_FIXTURE_TEST_CASE(SubgridConstruction, TestFixtureAd<SetupSimple>)
@ -115,29 +111,24 @@ BOOST_FIXTURE_TEST_CASE(SubgridConstruction, TestFixtureAd<SetupSimple>)
BOOST_FIXTURE_TEST_CASE(SurfaceDensity, TestFixture<SetupSimple>)
{
Opm::BlackoilPropsAd boprops_ad(props);
const double* rho0 = props .surfaceDensity();
const double* rho0AD = boprops_ad.surfaceDensity();
enum { Water = Opm::BlackoilPropsAd::Water };
BOOST_CHECK_EQUAL(rho0AD[ Water ], rho0[ Water ]);
enum { Water = Opm::BlackoilPropsAdFromDeck::Water };
BOOST_CHECK_EQUAL(rho0AD[ Water ], 1000.0);
enum { Oil = Opm::BlackoilPropsAd::Oil };
BOOST_CHECK_EQUAL(rho0AD[ Oil ], rho0[ Oil ]);
enum { Oil = Opm::BlackoilPropsAdFromDeck::Oil };
BOOST_CHECK_EQUAL(rho0AD[ Oil ], 800.0);
enum { Gas = Opm::BlackoilPropsAd::Gas };
BOOST_CHECK_EQUAL(rho0AD[ Gas ], rho0[ Gas ]);
enum { Gas = Opm::BlackoilPropsAdFromDeck::Gas };
BOOST_CHECK_EQUAL(rho0AD[ Gas ], 1.0);
}
BOOST_FIXTURE_TEST_CASE(ViscosityValue, TestFixture<SetupSimple>)
{
Opm::BlackoilPropsAd boprops_ad(props);
const Opm::BlackoilPropsAdFromDeck::Cells cells(5, 0);
const Opm::BlackoilPropsAd::Cells cells(5, 0);
typedef Opm::BlackoilPropsAd::V V;
typedef Opm::BlackoilPropsAdFromDeck::V V;
V Vpw;
Vpw.resize(cells.size());
@ -150,7 +141,11 @@ BOOST_FIXTURE_TEST_CASE(ViscosityValue, TestFixture<SetupSimple>)
// standard temperature
V T = V::Constant(cells.size(), 273.15+20);
const Opm::BlackoilPropsAd::V VmuWat = boprops_ad.muWat(Vpw, T, cells);
BOOST_REQUIRE_EQUAL(Vpw.size(), cells.size());
const Opm::BlackoilPropsAdFromDeck::V VmuWat = boprops_ad.muWat(Vpw, T, cells);
BOOST_REQUIRE_EQUAL(Vpw.size(), cells.size());
// Zero pressure dependence in water viscosity
for (V::Index i = 0, n = VmuWat.size(); i < n; ++i) {
@ -161,11 +156,9 @@ BOOST_FIXTURE_TEST_CASE(ViscosityValue, TestFixture<SetupSimple>)
BOOST_FIXTURE_TEST_CASE(ViscosityAD, TestFixture<SetupSimple>)
{
Opm::BlackoilPropsAd boprops_ad(props);
const Opm::BlackoilPropsAdFromDeck::Cells cells(5, 0);
const Opm::BlackoilPropsAd::Cells cells(5, 0);
typedef Opm::BlackoilPropsAd::V V;
typedef Opm::BlackoilPropsAdFromDeck::V V;
V Vpw;
Vpw.resize(cells.size());
@ -178,13 +171,13 @@ BOOST_FIXTURE_TEST_CASE(ViscosityAD, TestFixture<SetupSimple>)
// standard temperature
V T = V::Constant(cells.size(), 273.15+20);
typedef Opm::BlackoilPropsAd::ADB ADB;
typedef Opm::BlackoilPropsAdFromDeck::ADB ADB;
const V VmuWat = boprops_ad.muWat(Vpw, T, cells);
for (V::Index i = 0, n = Vpw.size(); i < n; ++i) {
const std::vector<int> bp(1, grid.c_grid()->number_of_cells);
const Opm::BlackoilPropsAd::Cells c(1, 0);
const Opm::BlackoilPropsAdFromDeck::Cells c(1, 0);
const V pw = V(1, 1) * Vpw[i];
const ADB Apw = ADB::variable(0, pw, bp);
const ADB AT = ADB::constant(T);

15
tests/test_rateconverter.cpp
View File

@ -28,7 +28,7 @@
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <boost/test/unit_test.hpp>
@ -69,8 +69,7 @@ struct TestFixture : public Setup
TestFixture()
: Setup()
, grid (deck)
, props(deck, eclState, *grid.c_grid(), param,
param.getDefault("init_rock", false))
, ad_props(deck, eclState, *grid.c_grid(), param.getDefault("init_rock", false))
{
}
@ -78,20 +77,19 @@ struct TestFixture : public Setup
using Setup::deck;
using Setup::eclState;
Opm::GridManager grid;
Opm::BlackoilPropertiesFromDeck props;
Opm::GridManager grid;
Opm::BlackoilPropsAdFromDeck ad_props;
};
BOOST_FIXTURE_TEST_CASE(Construction, TestFixture<SetupSimple>)
{
typedef std::vector<int> Region;
typedef Opm::BlackoilPropsAd Props;
typedef Opm::BlackoilPropsAdFromDeck Props;
typedef Opm::RateConverter::
SurfaceToReservoirVoidage<Props, Region> RCvrt;
Region reg{ 0 };
Props ad_props(props);
RCvrt cvrt(ad_props, reg);
}
@ -100,12 +98,11 @@ BOOST_FIXTURE_TEST_CASE(TwoPhaseII, TestFixture<SetupSimple>)
{
// Immiscible and incompressible two-phase fluid
typedef std::vector<int> Region;
typedef Opm::BlackoilPropsAd Props;
typedef Opm::BlackoilPropsAdFromDeck Props;
typedef Opm::RateConverter::
SurfaceToReservoirVoidage<Props, Region> RCvrt;
Region reg{ 0 };
Props ad_props(props);
RCvrt cvrt(ad_props, reg);
Opm::BlackoilState x;