From 1cf1e3383bd9368b328cc1f53bb44a7c5ad541d4 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 27 May 2013 14:55:32 +0200 Subject: [PATCH] Add new interface to BlackoilPropsAdFromDeck -- not finished -- --- CMakeLists_files.cmake | 2 + opm/autodiff/BlackoilPropsAdFromDeck.cpp | 636 +++++++++++++++++++++++ opm/autodiff/BlackoilPropsAdFromDeck.hpp | 246 +++++++++ 3 files changed, 884 insertions(+) create mode 100644 opm/autodiff/BlackoilPropsAdFromDeck.cpp create mode 100644 opm/autodiff/BlackoilPropsAdFromDeck.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index c4bc41b60..d4bc5a097 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -32,6 +32,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/SimulatorCompressibleAd.cpp opm/autodiff/SimulatorIncompTwophaseAdfi.cpp opm/autodiff/TransportSolverTwophaseAd.cpp + opm/autodiff/BlackoilPropsAdFromDeck.cpp ) # originally generated with the command: @@ -69,6 +70,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/AutoDiffHelpers.hpp opm/autodiff/AutoDiff.hpp opm/autodiff/BlackoilPropsAd.hpp + opm/autodiff/BlackoilPropsAdFromDeck.hpp opm/autodiff/BlackoilPropsAdInterface.hpp opm/autodiff/GeoProps.hpp opm/autodiff/ImpesTPFAAD.hpp diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp new file mode 100644 index 000000000..805aaeff8 --- /dev/null +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -0,0 +1,636 @@ +/* + 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 . +*/ + +#include + +#include +#include +#include +#include +#include + +namespace Opm +{ + + // Making these typedef to make the code more readable. + typedef BlackoilPropsAdFromDeck::ADB ADB; + typedef BlackoilPropsAdFromDeck::V V; + typedef Eigen::Array Block; + + /// Constructor wrapping an opm-core black oil interface. + BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropertiesInterface& props) + { + if (init_rock){ + rock_.init(deck, grid); + } + + + + phase_usage_ = phaseUsageFromDeck(deck); + // 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 { + THROW("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 { + THROW("Input is missing PVDG or PVTG\n"); + } + } + + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, 200); + + if (pvt_.numPhases() != satprops_->numPhases()) { + THROW("BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") 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 // + //////////////////////////// + + + // ------ 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]) { + THROW("Cannot call muWat(): water phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + double mu[n]; + double dmudp[n]; + double dmudr[n]; + double rs[n]; + + props_[phase_usage_.phase_pos[Water]]->mu(n, pw.data(), rs, mu,dmudp,dmudr); + + 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]) { + THROW("Cannot call muOil(): oil phase not present."); + } + const int n = cells.size(); + ASSERT(po.size() == n); + double mu[n]; + double dmudp[n]; + double dmudr[n]; + + props_[phase_usage_.phase_pos[Oil]]->mu(n, po.data(), rs.data(), mu,dmudp,dmudr); + + 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]) { + THROW("Cannot call muGas(): gas phase not present."); + } + const int n = cells.size(); + ASSERT(po.size() == n); + double mu[n]; + double dmudp[n]; + double dmudr[n]; + double rs[n]; + + props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.data(), rs.data(), mu,dmudp,dmudr); + + 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]) { + THROW("Cannot call muWat(): water phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + double mu[n]; + double dmudp[n]; + double dmudr[n]; + double rs[n]; + + props_[phase_usage_.phase_pos[Water]]->mu(n, pw.data(), rs, mu,dmudp,dmudr); + + ADB::M dmu_diag = spdiag(dmudp); + const int num_blocks = pw.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dmu_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]) { + THROW("Cannot call muOil(): oil phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + double mu[n]; + double dmudp[n]; + double dmudr[n]; + + props_[phase_usage_.phase_pos[Oil]]->mu(n, po.data(), rs, mu,dmudp,dmudr); + + ADB::M dmu_diag = spdiag(dmudp); + ADB::M dmu_drs_diag = spdiag(dmudr); + const int num_blocks = po.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dmu_diag * po.derivative()[block] + dmu_drs_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]) { + THROW("Cannot call muGas(): gas phase not present."); + } + const int n = cells.size(); + ASSERT(pg.value().size() == n); + double mu[n]; + double dmudp[n]; + double dmudr[n]; + + props_[phase_usage_.phase_pos[Gas]]->mu(n, pg.data(), rs, mu,dmudp,dmudr); + + ADB::M dmu_diag = spdiag(dmudp); + ADB::M dmu_drs_diag = spdiag(dmudr); + const int num_blocks = pg.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dmu_diag * pg.derivative()[block] + dmu_drs_diag * rs.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]) { + THROW("Cannot call bWat(): water phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + + double b[n]; + double dbdr[n]; + double dbdp[n]; + double rs[n]; + + props_[phase_usage_.phase_pos[Water]]->b(n, pw, rs, b,dbdp,dbdr); + + 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]) { + THROW("Cannot call bOil(): oil phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + + double b[n]; + double dbdr[n]; + double dbdp[n]; + + props_[phase_usage_.phase_pos[Oil]]->b(n, po, rs, b,dbdp,dbdr); + + 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]) { + THROW("Cannot call bGas(): gas phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + + double b[n]; + double dbdr[n]; + double dbdp[n]; + + props_[phase_usage_.phase_pos[Gas]]->b(n, pg, rs, b,dbdp,dbdr); + + 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]) { + THROW("Cannot call muWat(): water phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + + double b[n]; + double dbdr[n]; + double dbdp[n]; + double rs[n]; + + props_[phase_usage_.phase_pos[Water]]->b(n, pw, rs, b,dbdp,dbdr); + + ADB::M db_diag = spdiag(dbdp); + const int num_blocks = pw.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = db_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]) { + THROW("Cannot call muOil(): oil phase not present."); + } + const int n = cells.size(); + ASSERT(po.size() == n); + + double b[n]; + double dbdr[n]; + double dbdp[n]; + + props_[phase_usage_.phase_pos[Oil]]->b(n, po, rs, b,dbdp,dbdr); + + ADB::M db_diag = spdiag(dbdp); + ADB::M db_dr_diag = spdiag(dbdr); + const int num_blocks = po.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = db_diag * po.derivative()[block] + db_dr_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]) { + THROW("Cannot call muGas(): gas phase not present."); + } + const int n = cells.size(); + ASSERT(pg.size() == n); + + double b[n]; + double dbdr[n]; + double dbdp[n]; + + props_[phase_usage_.phase_pos[Gas]]->b(n, pg, rs, b,dbdp,dbdr); + + ADB::M db_diag = spdiag(dbdp); + ADB::M db_dr_diag = spdiag(dbdr); + const int num_blocks = pg.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = db_diag * pg.derivative()[block] + db_dr_diag * rs.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]) { + THROW("Cannot call muOil(): oil phase not present."); + } + const int n = cells.size(); + ASSERT(pw.size() == n); + double rbub[n]; + double drbubdp[n]; + props_[phase] ->rbub(n,po,rbub,drbubdp); + 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]) { + THROW("Cannot call muOil(): oil phase not present."); + } + const int n = cells.size(); + ASSERT(po.size() == n); + double rbub[n]; + double drbubdp[n]; + props_[phase] ->rbub(n,po,rbub,drbubdp); + ADB::M drbub_diag = spdiag(drbubdp); + const int num_blocks = po.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = drbub_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 BlackoilPropsAdFromDeck::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 (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); + props_.relperm(n, s_all.data(), cells.data(), kr.data(), 0); + std::vector 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 BlackoilPropsAdFromDeck::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 (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 { + THROW("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); + props_.relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data()); + const int num_blocks = so.numBlocks(); + std::vector 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 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 + diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp new file mode 100644 index 000000000..69d0647ae --- /dev/null +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -0,0 +1,246 @@ +/* + 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 . +*/ + +#ifndef OPM_BLACKOILPROPSAD_HEADER_INCLUDED +#define OPM_BLACKOILPROPSAD_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include + + +namespace Opm +{ + + /// class BlackoilPropertiesInterface; + + /// 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 BlackoilPropsAdFromDeck : public BlackoilPropsAdInterface + { + public: + /// Constructor wrapping an opm-core black oil interface. + BlackoilPropsAdFromDeck(const EclipseGridParser& deck, + const UnstructuredGrid& grid, 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 AutoDiff::ForwardBlock ADB; + typedef ADB::V V; + typedef std::vector Cells; + + + // ------ 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 ------ +#if 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. + 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; +#endif + + // ------ 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 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 relperm(const ADB& sw, + const ADB& so, + const ADB& sg, + const Cells& cells) const; + + private: + RockFromDeck rock_; + boost::scoped_ptr satprops_; + PhaseUsage phase_usage_; + std::vector > props_; + double densities_[MaxNumPhases]; + + + }; + +} // namespace Opm + +#endif // OPM_BLACKOILPROPSAD_HEADER_INCLUDED