diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 5198da706..32a7c1df3 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/examples/sim_2p_comp_ad.cpp b/examples/sim_2p_comp_ad.cpp deleted file mode 100644 index 558656b2e..000000000 --- a/examples/sim_2p_comp_ad.cpp +++ /dev/null @@ -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 . -*/ - - -#if HAVE_CONFIG_H -#include "config.h" -#endif // HAVE_CONFIG_H - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include - -#include -#include -#include -#include - -#include -#include -#include - - -#include - -#include -#include -#include -#include -#include - - -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 grid; - std::unique_ptr props; - std::unique_ptr 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("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; -} - diff --git a/opm/autodiff/BlackoilPropsAd.cpp b/opm/autodiff/BlackoilPropsAd.cpp deleted file mode 100644 index 436d7ca85..000000000 --- a/opm/autodiff/BlackoilPropsAd.cpp +++ /dev/null @@ -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 . -*/ - -#include - -#include -#include -#include -#include -#include - -namespace Opm -{ - - // Making these typedef to make the code more readable. - typedef BlackoilPropsAd::ADB ADB; - typedef BlackoilPropsAd::V V; - typedef Eigen::Array 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(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& /*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& /*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 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& 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 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 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& 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 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& /*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& /*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 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& /*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 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 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& /*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 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(po); - static_cast(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(po); - static_cast(so); - static_cast(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(po); - static_cast(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(po); - static_cast(so); - static_cast(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(po); - static_cast(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(po); - static_cast(so); - static_cast(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(po); - static_cast(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(po); - static_cast(so); - static_cast(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 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 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 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 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 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 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 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 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& /* saturation */, - const std::vector& /* cells */) - { - OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support hysteresis."); - } - - /// Update for max oil saturation. - void BlackoilPropsAd::updateSatOilMax(const std::vector& /*saturation*/) - { - OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support this functionality."); - } - -} // namespace Opm - diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp deleted file mode 100644 index 0511355ee..000000000 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_BLACKOILPROPSAD_HEADER_INCLUDED -#define OPM_BLACKOILPROPSAD_HEADER_INCLUDED - -#include -#include -#include - -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 ADB; - typedef ADB::V V; - typedef std::vector 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& 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& 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& 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& 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& 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& 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& 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& 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 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; - - /// 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 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& saturation, - const std::vector& cells); - - - /// Update for max oil saturation. - void updateSatOilMax(const std::vector& saturation); - - private: - const BlackoilPropertiesInterface& props_; - PhaseUsage pu_; - }; - -} // namespace Opm - -#endif // OPM_BLACKOILPROPSAD_HEADER_INCLUDED diff --git a/opm/autodiff/SimulatorCompressibleAd.cpp b/opm/autodiff/SimulatorCompressibleAd.cpp deleted file mode 100644 index a97d2118a..000000000 --- a/opm/autodiff/SimulatorCompressibleAd.cpp +++ /dev/null @@ -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 . -*/ - -#if HAVE_CONFIG_H -#include "config.h" -#endif // HAVE_CONFIG_H - -#include -#include -#include - -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include - -#include -#include - -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include - - -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 > columns_; - // Misc. data - std::vector 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 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 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& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(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 transport_src; - - // Initialisation. - std::vector 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 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 fractional_flows; - std::vector 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 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 diff --git a/opm/autodiff/SimulatorCompressibleAd.hpp b/opm/autodiff/SimulatorCompressibleAd.hpp deleted file mode 100644 index c7895f796..000000000 --- a/opm/autodiff/SimulatorCompressibleAd.hpp +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED -#define OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED - -#include -#include - -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 pimpl_; - }; - -} // namespace Opm - -#endif // OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED diff --git a/tests/fluid.data b/tests/fluid.data index 1843d8b10..81f49299c 100644 --- a/tests/fluid.data +++ b/tests/fluid.data @@ -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 / diff --git a/tests/test_boprops_ad.cpp b/tests/test_boprops_ad.cpp index 1066bb1ab..e83324f43 100644 --- a/tests/test_boprops_ad.cpp +++ b/tests/test_boprops_ad.cpp @@ -26,13 +26,11 @@ #define BOOST_TEST_MODULE FluidPropertiesTest -#include #include #include #include -#include #include #include @@ -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 @@ -105,7 +102,6 @@ struct TestFixtureAd : public Setup BOOST_FIXTURE_TEST_CASE(Construction, TestFixture) { - Opm::BlackoilPropsAd boprops_ad(props); } BOOST_FIXTURE_TEST_CASE(SubgridConstruction, TestFixtureAd) @@ -115,29 +111,24 @@ BOOST_FIXTURE_TEST_CASE(SubgridConstruction, TestFixtureAd) BOOST_FIXTURE_TEST_CASE(SurfaceDensity, TestFixture) { - 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) { - 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) // 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) BOOST_FIXTURE_TEST_CASE(ViscosityAD, TestFixture) { - 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) // 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 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); diff --git a/tests/test_rateconverter.cpp b/tests/test_rateconverter.cpp index d4c70346a..08edb2c73 100644 --- a/tests/test_rateconverter.cpp +++ b/tests/test_rateconverter.cpp @@ -28,7 +28,7 @@ #include -#include +#include #include @@ -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) { typedef std::vector Region; - typedef Opm::BlackoilPropsAd Props; + typedef Opm::BlackoilPropsAdFromDeck Props; typedef Opm::RateConverter:: SurfaceToReservoirVoidage RCvrt; Region reg{ 0 }; - Props ad_props(props); RCvrt cvrt(ad_props, reg); } @@ -100,12 +98,11 @@ BOOST_FIXTURE_TEST_CASE(TwoPhaseII, TestFixture) { // Immiscible and incompressible two-phase fluid typedef std::vector Region; - typedef Opm::BlackoilPropsAd Props; + typedef Opm::BlackoilPropsAdFromDeck Props; typedef Opm::RateConverter:: SurfaceToReservoirVoidage RCvrt; Region reg{ 0 }; - Props ad_props(props); RCvrt cvrt(ad_props, reg); Opm::BlackoilState x;