diff --git a/Makefile.am b/Makefile.am index c1c85dec..712d7c65 100644 --- a/Makefile.am +++ b/Makefile.am @@ -92,6 +92,7 @@ opm/core/pressure/tpfa/compr_source.c \ opm/core/pressure/tpfa/ifs_tpfa.c \ opm/core/pressure/tpfa/trans_tpfa.c \ opm/core/pressure/well.c \ +opm/core/simulator/SimulatorCompressibleTwophase.cpp \ opm/core/simulator/SimulatorIncompTwophase.cpp \ opm/core/simulator/SimulatorReport.cpp \ opm/core/simulator/SimulatorTimer.cpp \ @@ -195,6 +196,7 @@ opm/core/pressure/tpfa/compr_source.h \ opm/core/pressure/tpfa/ifs_tpfa.h \ opm/core/pressure/tpfa/trans_tpfa.h \ opm/core/simulator/BlackoilState.hpp \ +opm/core/simulator/SimulatorCompressibleTwophase.hpp \ opm/core/simulator/SimulatorReport.hpp \ opm/core/simulator/SimulatorIncompTwophase.hpp \ opm/core/simulator/SimulatorTimer.hpp \ diff --git a/README b/README index 9653fb74..8c7c522f 100644 --- a/README +++ b/README @@ -1,28 +1,42 @@ -These are the release notes for opm-core - Open Porous Media Core Library +============================== + +These are release notes for opm-core. + CONTENT +------- + opm-core is the core library within OPM and contains the following -* Eclipse preprosessing +* Eclipse deck input and preprosessing * Fluid properties (basic PVT models and rock properties) -* Grid (generates different grids) -* Linear Algerbra (interface to different linear solvers) -* Pressure solvers (different discretization schemes, different flow models) -* Simulator (some basic examples of simulators based on sequential schemes) -* Transport solvers (different discretization schemes, different flow models) -* Utility (input and output processing, unit conversion) +* Grid handling (cornerpoint grids, unstructured grid interface) +* Linear Algebra (interface to different linear solvers) +* Pressure solvers (various discretization schemes, flow models) +* Simulators (some basic examples of simulators based on sequential splitting schemes) +* Transport solvers (various discretization schemes, flow models) +* Utilities (input and output processing, unit conversion) * Wells (basic well handling) -ON WHAT PLATFORMS DOES IT RUN? -The opm-core module is designed to run on linux platforms. No efforts have -been made to ensure that the code will compile and run on windows platforms. -DOCUMENTATION -Efforts have been made to document the code with Doxygen. +LICENSE +------- + +The library is distributed under the GNU General Public License, +version 3 or later (GPLv3+). + + +PLATFORMS +--------- + +The opm-core module is designed to run on Linux platforms. It is also +regularly run on Mac OS X. No efforts have been made to ensure that +the code will compile and run on windows platforms. + DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS (Debian Squeeze/Ubuntu Precise) +--------------------------------------------------------------------------- # packages necessary for building sudo apt-get install -y build-essential gfortran pkg-config libtool \ @@ -40,7 +54,9 @@ sudo apt-get install -y libboost-all-dev libsuperlu3-dev libsuitesparse-dev # libraries necessary for OPM sudo apt-get install -y libxml0-dev + DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS +----------------------------------------- # libraries sudo zypper install blas libblas3 lapack liblapack3 libboost libxml2 umfpack @@ -48,7 +64,11 @@ sudo zypper install blas libblas3 lapack liblapack3 libboost libxml2 umfpack # tools sudo zypper install gcc automake autoconf git doxygen + RETRIEVING AND BUILDING DUNE PREREQUISITES +------------------------------------------ + +(only necessary if you want to use opm-core as a dune module) # trust DUNE certificate (sic) echo p | svn list https://svn.dune-project.org/svn/dune-common @@ -67,20 +87,39 @@ for module in common istl geometry grid localfunctions; do --make-opts="-j -l 0.8" autogen : configure : make done -DOWNLOAD opm-core +DOWNLOADING +----------- + +For a read-only download: git clone git://github.com/OPM/opm-core.git -BUILDING opm-core +If you want to contribute, fork OPM/opm-core on github. -# note: this is done from the parent directory of opm-core -env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=opm-core \ - --configure-opts="" --make-opts="-j -l 0.8" autogen : configure : make -or, without using dunecontrol: +BUILDING +-------- -pushd opm-core -autoreconf -i -./configure -make -sudo make install +(standalone opm-core:) + + cd ../opm-core + autoreconf -i + ./configure + make + sudo make install + +(using opm-core as a dune module:) + + # note: this is done from the parent directory of opm-core + env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=opm-core \ + --configure-opts="" --make-opts="-j -l 0.8" autogen : configure : make + + + +DOCUMENTATION +------------- + +Efforts have been made to document the code with Doxygen. +In order to build the documentation, enter the command +$ doxygen +in the topmost directory. diff --git a/examples/Makefile.am b/examples/Makefile.am index 6d664bdf..b01f7674 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -17,6 +17,7 @@ $(BOOST_SYSTEM_LIB) noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ +sim_2p_comp_reorder \ sim_2p_incomp_reorder \ sim_wateroil \ wells_example @@ -29,6 +30,7 @@ wells_example # Please maintain sort order from "noinst_PROGRAMS". refine_wells_SOURCES = refine_wells.cpp +sim_2p_comp_reorder_SOURCES = sim_2p_comp_reorder.cpp sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp sim_wateroil_SOURCES = sim_wateroil.cpp wells_example_SOURCES = wells_example.cpp diff --git a/examples/sim_2p_comp_reorder.cpp b/examples/sim_2p_comp_reorder.cpp new file mode 100644 index 00000000..9f2b516b --- /dev/null +++ b/examples/sim_2p_comp_reorder.cpp @@ -0,0 +1,285 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#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 + + +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) +{ + 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; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr rock_comp; + BlackoilState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + deck.reset(new EclipseGridParser(deck_filename)); + // Grid init + grid.reset(new GridManager(*deck)); + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid())); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(*deck)); + // Gravity. + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); + } + initBlackoilSurfvol(*grid->c_grid(), *props, state); + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, 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; + + // Initialising src + int num_cells = grid->c_grid()->number_of_cells; + std::vector src(num_cells, 0.0); + if (use_deck) { + // Do nothing, wells will be the driving force, not source terms. + } else { + // Compute pore volumes, in order to enable specifying injection rate + // terms of total pore volume. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + FlowBCManager bcs; + if (param.getDefault("use_pside", false)) { + int pside = param.get("pside"); + double pside_pressure = param.get("pside_pressure"); + bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure); + } + + // 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 (...) { + THROW("Creating directories failed: " << fpath); + } + std::string filename = output_dir + "/epoch_timing.param"; + epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out); + // open file to clean it. The file is appended to in SimulatorTwophase + filename = output_dir + "/step_timing.param"; + std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out); + step_os.close(); + param.writeParam(output_dir + "/simulation.param"); + } + + + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << " (number of epochs: " + << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + + SimulatorReport rep; + if (!use_deck) { + // Simple simulation without a deck. + WellsManager wells; // no wells. + SimulatorCompressibleTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); + SimulatorTimer simtimer; + simtimer.init(param); + warnIfUnusedParams(param); + WellState well_state; + well_state.init(0, state); + rep = simulator.run(simtimer, state, well_state); + } else { + // With a deck, we may have more epochs etc. + WellState well_state; + int step = 0; + SimulatorTimer simtimer; + // Use timer for last epoch to obtain total time. + deck->setCurrentEpoch(deck->numberOfEpochs() - 1); + simtimer.init(*deck); + const double total_time = simtimer.totalTime(); + for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { + // Set epoch index. + deck->setCurrentEpoch(epoch); + + // Update the timer. + if (deck->hasField("TSTEP")) { + simtimer.init(*deck); + } else { + if (epoch != 0) { + THROW("No TSTEP in deck for epoch " << epoch); + } + simtimer.init(param); + } + simtimer.setCurrentStepNum(step); + simtimer.setTotalTime(total_time); + + // Report on start of epoch. + std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" + << "\n (number of steps: " + << simtimer.numSteps() - step << ")\n\n" << std::flush; + + // Create new wells, well_state + WellsManager wells(*deck, *grid->c_grid(), props->permeability()); + // @@@ HACK: we should really make a new well state and + // properly transfer old well state to it every epoch, + // since number of wells may change etc. + if (epoch == 0) { + well_state.init(wells.c_wells(), state); + } + + // Create and run simulator. + SimulatorCompressibleTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); + if (epoch == 0) { + warnIfUnusedParams(param); + } + SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + if (output) { + epoch_rep.reportParam(epoch_os); + } + // Update total timing report and remember step number. + rep += epoch_rep; + step = simtimer.currentStepNum(); + } + } + + std::cout << "\n\n================ End of simulation ===============\n\n"; + rep.report(std::cout); + + if (output) { + std::string filename = output_dir + "/walltime.param"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + rep.reportParam(tot_os); + } + +} diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 0ef94b05..f7e0172f 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -141,7 +141,6 @@ main(int argc, char** argv) std::cout << "--------------- Reading parameters ---------------" << std::endl; // Reading various control parameters. - const bool use_reorder = param.getDefault("use_reorder", true); const bool output = param.getDefault("output", true); std::string output_dir; int output_interval = 1; @@ -229,28 +228,8 @@ main(int argc, char** argv) } } bool use_segregation_split = false; - bool use_column_solver = false; - bool use_gauss_seidel_gravity = false; - if (use_gravity && use_reorder) { + if (use_gravity) { use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split); - if (use_segregation_split) { - use_column_solver = param.getDefault("use_column_solver", use_column_solver); - if (use_column_solver) { - use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity); - } - } - } - - // Check that rock compressibility is not used with solvers that do not handle it. - // int nl_pressure_maxiter = 0; - // double nl_pressure_tolerance = 0.0; - if (rock_comp->isActive()) { - THROW("No rock compressibility in comp. pressure solver yet."); - if (!use_reorder) { - THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet."); - } - // nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10); - // nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal } // Source-related variables init. @@ -262,11 +241,11 @@ main(int argc, char** argv) // Extra rock init. std::vector porevol; if (rock_comp->isActive()) { - THROW("CompressibleTpfa solver does not handle this."); computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); } else { computePorevolume(*grid->c_grid(), props->porosity(), porevol); } + std::vector initial_porevol = porevol; double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); @@ -293,20 +272,20 @@ main(int argc, char** argv) const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0); const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20); const double *grav = use_gravity ? &gravity[0] : 0; - Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, linsolver, + Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver, nl_press_res_tol, nl_press_change_tol, nl_press_maxiter, grav, wells->c_wells()); // Reordering solver. const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const int nl_maxiter = param.getDefault("nl_maxiter", 30); Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); - if (use_gauss_seidel_gravity) { - reorder_model.initGravity(); + if (use_segregation_split) { + reorder_model.initGravity(grav); } // Column-based gravity segregation solver. std::vector > columns; - if (use_column_solver) { + if (use_segregation_split) { Opm::extractColumn(*grid->c_grid(), columns); } @@ -409,6 +388,11 @@ main(int argc, char** argv) Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, wells->c_wells(), well_state.perfRates(), reorder_src); + // Compute new porevolumes after pressure solve, if necessary. + if (rock_comp->isActive()) { + initial_porevol = porevol; + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } // Solve transport. transport_timer.start(); double stepsize = simtimer.currentStepLength(); @@ -419,11 +403,13 @@ main(int argc, char** argv) for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { // Note that for now we do not handle rock compressibility, // although the transport solver should be able to. - reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0], - &porevol[0], &porevol[0], &reorder_src[0], stepsize, state.saturation()); + reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], + &porevol[0], &initial_porevol[0], &reorder_src[0], stepsize, + state.saturation(), state.surfacevol()); // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); if (use_segregation_split) { - reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation()); + reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0], + stepsize, state.saturation(), state.surfacevol()); } } transport_timer.stop(); diff --git a/generate_doc_figures.py b/generate_doc_figures.py index ae42709d..8d4931ee 100755 --- a/generate_doc_figures.py +++ b/generate_doc_figures.py @@ -1,6 +1,24 @@ + +# This script generate the illustration pictures for the documentation. +# +# To run this script, you have to install paraview, see: +# +# http://www.paraview.org/paraview/resources/software.php +# +# Eventually, set up the paths (figure_path, tutorial_data_path, tutorial_path) according to your own installation. +# (The default values should be ok.) +# +# Make sure that pvpython is in your path of executables. +# +# Run the following command at the root of the directory where +# opm-core is installed (which also corresponds to the directory where +# generate_doc_figures is located): +# +# pvpython generate_doc_figures.py +# + from subprocess import call from paraview.simple import * -# from paraview import servermanager from os import remove, mkdir, curdir from os.path import join, isdir @@ -13,7 +31,6 @@ collected_garbage_file = [] if not isdir(figure_path): mkdir(figure_path) -# connection = servermanager.Connect() # tutorial 1 call(join(tutorial_path, "tutorial1")) diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index f84d1c55..04d8cd04 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -166,6 +166,16 @@ namespace Opm // Construct tensor grid from deck. void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck) { + // Extract logical cartesian size. + std::vector dims; + if (deck.hasField("DIMENS")) { + dims = deck.getIntegerValue("DIMENS"); + } else if (deck.hasField("SPECGRID")) { + dims = deck.getSPECGRID().dimensions; + } else { + THROW("Deck must have either DIMENS or SPECGRID."); + } + // Extract coordinates (or offsets from top, in case of z). const std::vector& dxv = deck.getFloatingPointValue("DXV"); const std::vector& dyv = deck.getFloatingPointValue("DYV"); @@ -174,6 +184,17 @@ namespace Opm std::vector y = coordsFromDeltas(dyv); std::vector z = coordsFromDeltas(dzv); + // Check that number of cells given are consistent with DIMENS/SPECGRID. + if (dims[0] != int(dxv.size())) { + THROW("Number of DXV data points do not match DIMENS or SPECGRID."); + } + if (dims[1] != int(dyv.size())) { + THROW("Number of DYV data points do not match DIMENS or SPECGRID."); + } + if (dims[2] != int(dzv.size())) { + THROW("Number of DZV data points do not match DIMENS or SPECGRID."); + } + // Extract top corner depths, if available. const double* top_depths = 0; std::vector top_depths_vec; diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 6568e8cc..7ba8903f 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -36,8 +36,6 @@ namespace Opm PermeabilityKind fillTensor(const EclipseGridParser& parser, std::vector*>& tensor, std::tr1::array& kmap); - - int numGlobalCells(const EclipseGridParser& parser); } // anonymous namespace @@ -87,7 +85,7 @@ namespace Opm double perm_threshold) { const int dim = 3; - const int num_global_cells = numGlobalCells(parser); + const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2]; const int nc = grid.number_of_cells; ASSERT (num_global_cells > 0); @@ -334,26 +332,6 @@ namespace Opm return kind; } - int numGlobalCells(const EclipseGridParser& parser) - { - int ngc = -1; - - if (parser.hasField("DIMENS")) { - const std::vector& - dims = parser.getIntegerValue("DIMENS"); - - ngc = dims[0] * dims[1] * dims[2]; - } - else if (parser.hasField("SPECGRID")) { - const SPECGRID& sgr = parser.getSPECGRID(); - - ngc = sgr.dimensions[ 0 ]; - ngc *= sgr.dimensions[ 1 ]; - ngc *= sgr.dimensions[ 2 ]; - } - - return ngc; - } } // anonymous namespace } // namespace Opm diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 154f0633..954e07fa 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -58,6 +59,7 @@ namespace Opm /// to change. CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, const double residual_tol, const double change_tol, @@ -66,6 +68,7 @@ namespace Opm const struct Wells* wells) : grid_(grid), props_(props), + rock_comp_props_(rock_comp_props), linsolver_(linsolver), residual_tol_(residual_tol), change_tol_(change_tol), @@ -74,8 +77,8 @@ namespace Opm wells_(wells), htrans_(grid.cell_facepos[ grid.number_of_cells ]), trans_ (grid.number_of_faces), - porevol_(grid.number_of_cells), - allcells_(grid.number_of_cells) + allcells_(grid.number_of_cells), + singular_(false) { if (wells_ && (wells_->number_of_phases != props.numPhases())) { THROW("Inconsistent number of phases specified (wells vs. props): " @@ -86,7 +89,12 @@ namespace Opm UnstructuredGrid* gg = const_cast(&grid_); tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]); tpfa_trans_compute(gg, &htrans_[0], &trans_[0]); - computePorevolume(grid_, props.porosity(), porevol_); + // If we have rock compressibility, pore volumes are updated + // in the compute*() methods, otherwise they are constant and + // hence may be computed here. + if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { + computePorevolume(grid_, props.porosity(), porevol_); + } for (int c = 0; c < grid.number_of_cells; ++c) { allcells_[c] = c; } @@ -182,6 +190,21 @@ namespace Opm + /// @brief After solve(), was the resulting pressure singular. + /// Returns true if the pressure is singular in the following + /// sense: if everything is incompressible and there are no + /// pressure conditions, the absolute values of the pressure + /// solution are arbitrary. (But the differences in pressure + /// are significant.) + bool CompressibleTpfa::singularPressure() const + { + return singular_; + } + + + + + /// Compute well potentials. void CompressibleTpfa::computeWellPotentials(const BlackoilState& state) { @@ -230,6 +253,9 @@ namespace Opm const WellState& /*well_state*/) { computeWellPotentials(state); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_); + } } @@ -252,6 +278,8 @@ namespace Opm // std::vector face_gravcap_; // std::vector wellperf_A_; // std::vector wellperf_phasemob_; + // std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + // std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. computeCellDynamicData(dt, state, well_state); computeFaceDynamicData(dt, state, well_state); computeWellDynamicData(dt, state, well_state); @@ -273,6 +301,8 @@ namespace Opm // std::vector cell_viscosity_; // std::vector cell_phasemob_; // std::vector cell_voldisc_; + // std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + // std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. const int nc = grid_.number_of_cells; const int np = props_.numPhases(); const double* cell_p = &state.pressure()[0]; @@ -296,6 +326,14 @@ namespace Opm // TODO: Check this! cell_voldisc_.clear(); cell_voldisc_.resize(nc, 0.0); + + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_); + rock_comp_.resize(nc); + for (int cell = 0; cell < nc; ++cell) { + rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]); + } + } } @@ -465,9 +503,20 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; - cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], - &face_gravcap_[0], cell_press, well_bhp, - &porevol_[0], h_); + int was_adjusted = 0; + if (! (rock_comp_props_ && rock_comp_props_->isActive())) { + was_adjusted = + cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], h_); + } else { + was_adjusted = + cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], &initial_porevol_[0], + &rock_comp_[0], h_); + } + singular_ = (was_adjusted == 1); } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 8233ee17..054ac2e4 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -33,6 +33,7 @@ namespace Opm class BlackoilState; class BlackoilPropertiesInterface; + class RockCompressibility; class LinearSolverInterface; class WellState; @@ -44,23 +45,25 @@ namespace Opm { public: /// Construct solver. - /// \param[in] grid A 2d or 3d grid. - /// \param[in] props Rock and fluid properties. - /// \param[in] linsolver Linear solver to use. - /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. - /// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller. - /// \param[in] maxiter Maximum acceptable number of iterations. - /// \param[in] gravity Gravity vector. If non-null, the array should - /// have D elements. - /// \param[in] wells The wells argument. Will be used in solution, - /// is ignored if NULL. - /// Note: this class observes the well object, and - /// makes the assumption that the well topology - /// and completions does not change during the - /// run. However, controls (only) are allowed - /// to change. - CompressibleTpfa(const UnstructuredGrid& grid, + /// \param[in] grid A 2d or 3d grid. + /// \param[in] props Rock and fluid properties. + /// \param[in] rock_comp_props Rock compressibility properties. May be null. + /// \param[in] linsolver Linear solver to use. + /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. + /// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller. + /// \param[in] maxiter Maximum acceptable number of iterations. + /// \param[in] gravity Gravity vector. If non-null, the array should + /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL. + /// Note: this class observes the well object, and + /// makes the assumption that the well topology + /// and completions does not change during the + /// run. However, controls (only) are allowed + /// to change. + CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, const double residual_tol, const double change_tol, @@ -68,8 +71,8 @@ namespace Opm const double* gravity, const Wells* wells); - /// Destructor. - ~CompressibleTpfa(); + /// Destructor. + ~CompressibleTpfa(); /// Solve the pressure equation by Newton-Raphson scheme. /// May throw an exception if the number of iterations @@ -78,6 +81,14 @@ namespace Opm BlackoilState& state, WellState& well_state); + /// @brief After solve(), was the resulting pressure singular. + /// Returns true if the pressure is singular in the following + /// sense: if everything is incompressible and there are no + /// pressure conditions, the absolute values of the pressure + /// solution are arbitrary. (But the differences in pressure + /// are significant.) + bool singularPressure() const; + private: void computePerSolveDynamicData(const double dt, const BlackoilState& state, @@ -101,28 +112,29 @@ namespace Opm void solveIncrement(); double residualNorm() const; double incrementNorm() const; - void computeResults(BlackoilState& state, + void computeResults(BlackoilState& state, WellState& well_state) const; // ------ Data that will remain unmodified after construction. ------ - const UnstructuredGrid& grid_; + const UnstructuredGrid& grid_; const BlackoilPropertiesInterface& props_; + const RockCompressibility* rock_comp_props_; const LinearSolverInterface& linsolver_; const double residual_tol_; const double change_tol_; const int maxiter_; const double* gravity_; // May be NULL const Wells* wells_; // May be NULL, outside may modify controls (only) between calls to solve(). - std::vector htrans_; - std::vector trans_ ; - std::vector porevol_; + std::vector htrans_; + std::vector trans_ ; std::vector allcells_; // ------ Internal data for the cfs_tpfa_res solver. ------ - struct cfs_tpfa_res_data* h_; + struct cfs_tpfa_res_data* h_; // ------ Data that will be modified for every solve. ------ std::vector wellperf_gpot_; + std::vector initial_porevol_; // ------ Data that will be modified for every solver iteration. ------ std::vector cell_A_; @@ -135,13 +147,15 @@ namespace Opm std::vector face_gravcap_; std::vector wellperf_A_; std::vector wellperf_phasemob_; + std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. // The update to be applied to the pressures (cell and bhp). std::vector pressure_increment_; - - - - - + // True if the matrix assembled would be singular but for the + // adjustment made in the cfs_*_assemble() calls. This happens + // if everything is incompressible and there are no pressure + // conditions. + bool singular_; }; } // namespace Opm diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 7a52db1e..eb8dc8f3 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , /* ---------------------------------------------------------------------- */ -void +int cfs_tpfa_res_assemble(struct UnstructuredGrid *G , double dt , struct cfs_tpfa_res_forces *forces , @@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , struct cfs_tpfa_res_data *h ) /* ---------------------------------------------------------------------- */ { - int res_is_neumann, well_is_neumann, c, np2; + int res_is_neumann, well_is_neumann, c, np2, singular; csrmatrix_zero( h->J); vector_zero (h->J->m, h->F); @@ -1207,9 +1207,76 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , assemble_sources(dt, forces->src, h); } - if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) { - h->J->sa[0] *= 2; + singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp; + if (singular) { + h->J->sa[0] *= 2.0; } + + return singular; +} + + +/* ---------------------------------------------------------------------- */ +int +cfs_tpfa_res_comprock_assemble( + struct UnstructuredGrid *G , + double dt , + struct cfs_tpfa_res_forces *forces , + const double *zc , + struct compr_quantities_gen *cq , + const double *trans , + const double *gravcap_f, + const double *cpress , + const double *wpress , + const double *porevol , + const double *porevol0 , + const double *rock_comp, + struct cfs_tpfa_res_data *h ) +/* ---------------------------------------------------------------------- */ +{ + /* We want to add this term to the usual residual: + * + * (porevol(pressure)-porevol(initial_pressure))/dt. + * + * Its derivative (for the diagonal term of the Jacobian) is: + * + * porevol(pressure)*rock_comp(pressure)/dt + */ + + int c, rock_is_incomp, singular; + size_t j; + double dpv; + + /* Assemble usual system (without rock compressibility). */ + singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, + cpress, wpress, porevol0, h); + + /* If we made a singularity-removing adjustment in the + regular assembly, we undo it here. */ + if (singular) { + h->J->sa[0] /= 2.0; + } + + /* Add new terms to residual and Jacobian. */ + rock_is_incomp = 1; + for (c = 0; c < G->number_of_cells; c++) { + j = csrmatrix_elm_index(c, c, h->J); + + dpv = (porevol[c] - porevol0[c]); + if (dpv != 0.0 || rock_comp[c] != 0.0) { + rock_is_incomp = 0; + } + + h->J->sa[j] += porevol[c] * rock_comp[c]; + h->F[c] += dpv; + } + + /* Re-do the singularity-removing adjustment if necessary */ + if (rock_is_incomp && singular) { + h->J->sa[0] *= 2.0; + } + + return rock_is_incomp && singular; } diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 4fde3ba3..5eddeed6 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , void cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); -void +/* Return value is 1 if the assembled matrix was adjusted to remove a + singularity. This happens if all fluids are incompressible and + there are no pressure conditions on wells or boundaries. + Otherwise return 0. */ +int cfs_tpfa_res_assemble(struct UnstructuredGrid *G, double dt, struct cfs_tpfa_res_forces *forces, @@ -72,6 +76,27 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G, const double *porevol, struct cfs_tpfa_res_data *h); +/* Return value is 1 if the assembled matrix was adjusted to remove a + singularity. This happens if all fluids are incompressible, the + rock is incompressible, and there are no pressure conditions on + wells or boundaries. + Otherwise return 0. */ +int +cfs_tpfa_res_comprock_assemble( + struct UnstructuredGrid *G, + double dt, + struct cfs_tpfa_res_forces *forces, + const double *zc, + struct compr_quantities_gen *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *wpress, + const double *porevol, + const double *porevol0, + const double *rock_comp, + struct cfs_tpfa_res_data *h); + void cfs_tpfa_res_flux(struct UnstructuredGrid *G , struct cfs_tpfa_res_forces *forces , diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index 94e3d638..0e69bc49 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -771,7 +771,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G , assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok); /* We want to solve a Newton step for the residual - * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible + * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_incompressible * */ diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp new file mode 100644 index 00000000..c5a77eb0 --- /dev/null +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -0,0 +1,533 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#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 + + +namespace Opm +{ + + class SimulatorCompressibleTwophase::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + 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_; + WellsManager& wells_manager_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + const double* gravity_; + // Solvers + CompressibleTpfa psolver_; + TransportModelCompressibleTwophase tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorCompressibleTwophase::SimulatorCompressibleTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp, wells_manager, src, bcs, linsolver, gravity)); + } + + + + + + SimulatorReport SimulatorCompressibleTwophase::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 (...) { + THROW("Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("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(); + 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 (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + 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) { + THROW("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) { + THROW("Failed to open " << fname); + } + wellreport.write(os); + } + + + + // \TODO: make CompressibleTpfa take src and bcs. + SimulatorCompressibleTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_(rock_comp), + wells_manager_(wells_manager), + wells_(wells_manager.c_wells()), + src_(src), + bcs_(bcs), + psolver_(grid, props, rock_comp, linsolver, + param.getDefault("nl_pressure_residual_tolerance", 0.0), + param.getDefault("nl_pressure_change_tolerance", 1.0), + param.getDefault("nl_pressure_maxiter", 10), + gravity, wells_manager.c_wells() /*, src, bcs*/), + 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 (...) { + THROW("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 SimulatorCompressibleTwophase::Impl::run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state) + { + std::vector transport_src; + + // Initialisation. + std::vector porevol; + if (rock_comp_ && rock_comp_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_, 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_satvol[2] = { 0.0 }; + double satvol[2] = { 0.0 }; + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); + std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init + << " " << init_satvol[1]/tot_porevol_init << std::endl; + 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.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); + + // 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; + } + } + + // 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_, + fractional_flows, + well_state.perfRates(), + 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_) { + THROW("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_ && rock_comp_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + } + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), 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; + } + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], + &porevol[0], &initial_porevol[0], &transport_src[0], stepsize, + state.saturation(), state.surfacevol()); + Opm::computeInjectedProduced(props_, + state.pressure(), state.surfacevol(), state.saturation(), + transport_src, stepsize, injected, produced); + if (gravity_ != 0 && use_segregation_split_) { + tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0], + 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.saturation(), satvol); + 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 << "\nVolume balance report (all numbers relative to total pore volume).\n"; + std::cout << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + std::cout << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + std::cout << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + std::cout << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + std::cout << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + std::cout << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; + std::cout << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::endl; + std::cout.precision(8); + + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + 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/core/simulator/SimulatorCompressibleTwophase.hpp b/opm/core/simulator/SimulatorCompressibleTwophase.hpp new file mode 100644 index 00000000..d38d74b2 --- /dev/null +++ b/opm/core/simulator/SimulatorCompressibleTwophase.hpp @@ -0,0 +1,99 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED +#define OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; +struct FlowBoundaryConditions; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class BlackoilPropertiesInterface; + 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 SimulatorCompressibleTwophase + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp if non-null, rock compressibility properties + /// \param[in] well_manager well manager, may manage no (null) wells + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorCompressibleTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + 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. + boost::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED diff --git a/opm/core/simulator/SimulatorIncompTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp index d4e7459f..1f7e1119 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -132,7 +132,7 @@ namespace Opm return pimpl_->run(timer, state, well_state); } - static void reportVolumes(std::ostream &os,double satvol[2],double tot_porevol_init, + static void reportVolumes(std::ostream &os, double satvol[2], double tot_porevol_init, double tot_injected[2], double tot_produced[2], double injected[2], double produced[2], double init_satvol[2]) @@ -164,7 +164,7 @@ namespace Opm << std::endl; os.precision(8); } - + static void outputStateVtk(const UnstructuredGrid& grid, const Opm::TwophaseState& state, const int step, @@ -180,7 +180,7 @@ namespace Opm catch (...) { THROW("Creating directories failed: " << fpath); } - vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu"; + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); @@ -193,7 +193,8 @@ namespace Opm dm["velocity"] = &cell_velocity; Opm::writeVtkData(grid, dm, vtkfile); } - static void outputVectorMatlab(const std::string name, + + static void outputVectorMatlab(const std::string& name, const std::vector& vec, const int step, const std::string& output_dir) @@ -214,29 +215,6 @@ namespace Opm } std::copy(vec.begin(), vec.end(), std::ostream_iterator(file, "\n")); } - - static void outputVectorMatlab(const std::string name, - const std::vector vec, - const int step, - const std::string& output_dir) - { - std::ostringstream fname; - fname << output_dir << "/" << name; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - THROW("Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - THROW("Failed to open " << fname.str()); - } - std::copy(vec.begin(), vec.end(), std::ostream_iterator(file, "\n")); - } - static void outputStateMatlab(const UnstructuredGrid& grid, const Opm::TwophaseState& state, @@ -408,7 +386,7 @@ namespace Opm 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; @@ -528,6 +506,7 @@ namespace Opm // Update pore volumes if rock is compressible. if (rock_comp_ && rock_comp_->isActive()) { + initial_porevol = porevol; computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); } @@ -543,11 +522,11 @@ namespace Opm std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; } for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { - tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], - stepsize, state.saturation()); + tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], + stepsize, state.saturation()); Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced); if (use_segregation_split_) { - tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation()); + tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation()); } watercut.push(timer.currentTime() + timer.currentStepLength(), produced[0]/(produced[0] + produced[1]), @@ -556,7 +535,7 @@ namespace Opm wellreport.push(props_, *wells_, state.saturation(), timer.currentTime() + timer.currentStepLength(), well_state.bhp(), well_state.perfRates()); - } + } } transport_timer.stop(); double tt = transport_timer.secsSinceStart(); @@ -569,11 +548,10 @@ namespace Opm tot_injected[1] += injected[1]; tot_produced[0] += produced[0]; tot_produced[1] += produced[1]; - //reportVolumes(std::cout,satvol[0],&tot_porevol_init[0],&tot_injected[0]); - Opm::reportVolumes(std::cout,satvol, tot_porevol_init, - tot_injected, tot_produced, - injected, produced, - init_satvol); + reportVolumes(std::cout,satvol, tot_porevol_init, + tot_injected, tot_produced, + injected, produced, + init_satvol); sreport.total_time = step_timer.secsSinceStart(); if (output_) { sreport.reportParam(tstep_os); diff --git a/opm/core/simulator/SimulatorIncompTwophase.hpp b/opm/core/simulator/SimulatorIncompTwophase.hpp index ff4fbb68..10c46ff0 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.hpp +++ b/opm/core/simulator/SimulatorIncompTwophase.hpp @@ -96,4 +96,4 @@ namespace Opm } // namespace Opm -#endif // OPM_SIMULATORTWOPHASE_HEADER_INCLUDED +#endif // OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index f55f67ed..697f94a5 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -38,7 +39,8 @@ namespace Opm typedef RegulaFalsi RootFinder; - TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid, + TransportModelCompressibleTwophase::TransportModelCompressibleTwophase( + const UnstructuredGrid& grid, const Opm::BlackoilPropertiesInterface& props, const double tol, const int maxit) @@ -51,6 +53,7 @@ namespace Opm dt_(0.0), saturation_(grid.number_of_cells, -1.0), fractionalflow_(grid.number_of_cells, -1.0), + gravity_(0), mob_(2*grid.number_of_cells, -1.0), ia_upw_(grid.number_of_cells + 1, -1), ja_upw_(grid.number_of_faces, -1), @@ -75,15 +78,15 @@ namespace Opm void TransportModelCompressibleTwophase::solve(const double* darcyflux, const double* pressure, - const double* surfacevol0, const double* porevolume0, const double* porevolume, const double* source, const double dt, - std::vector& saturation) + std::vector& saturation, + std::vector& surfacevol) { darcyflux_ = darcyflux; - surfacevol0_ = surfacevol0; + surfacevol0_ = &surfacevol[0]; porevolume0_ = porevolume0; porevolume_ = porevolume; source_ = source; @@ -93,6 +96,11 @@ namespace Opm props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); + // Check immiscibility requirement (only done for first cell). + if (A_[1] != 0.0 || A_[2] != 0.0) { + THROW("TransportModelCompressibleTwophase requires a property object without miscibility."); + } + std::vector seq(grid_.number_of_cells); std::vector comp(grid_.number_of_cells + 1); int ncomp; @@ -107,6 +115,9 @@ namespace Opm &ia_downw_[0], &ja_downw_[0]); reorderAndTransport(grid_, darcyflux); toBothSat(saturation_, saturation); + + // Compute surface volume as a postprocessing step from saturation and A_ + computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]); } // Residual function r(s) for a single-cell implicit Euler transport @@ -375,20 +386,24 @@ namespace Opm - void TransportModelCompressibleTwophase::initGravity() + void TransportModelCompressibleTwophase::initGravity(const double* grav) { // Set up transmissibilities. std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); const int nf = grid_.number_of_faces; trans_.resize(nf); + gravflux_.resize(nf); tpfa_htrans_compute(const_cast(&grid_), props_.permeability(), &htrans[0]); tpfa_trans_compute(const_cast(&grid_), &htrans[0], &trans_[0]); + + // Remember gravity vector. + gravity_ = grav; } - void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav) + void TransportModelCompressibleTwophase::initGravityDynamic() { // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] @@ -410,7 +425,7 @@ namespace Opm for (int ci = 0; ci < 2; ++ci) { double gdz = 0.0; for (int d = 0; d < dim; ++d) { - gdz += grav[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]); + gdz += gravity_[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]); } gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]); } @@ -494,11 +509,11 @@ namespace Opm const double* pressure, const double* porevolume0, const double dt, - const double* grav, - std::vector& saturation) + std::vector& saturation, + std::vector& surfacevol) { // Assume that solve() has already been called, so that A_ is current. - initGravityDynamic(grav); + initGravityDynamic(); // Initialize mobilities. const int nc = grid_.number_of_cells; @@ -531,6 +546,9 @@ namespace Opm std::cout << "Gauss-Seidel column solver average iterations: " << double(num_iters)/double(columns.size()) << std::endl; toBothSat(saturation_, saturation); + + // Compute surface volume as a postprocessing step from saturation and A_ + computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]); } } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 19c44b22..4ee93e0a 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -30,6 +30,8 @@ namespace Opm class BlackoilPropertiesInterface; + /// Implements a reordering transport solver for compressible, + /// non-miscible two-phase flow. class TransportModelCompressibleTwophase : public TransportModelInterface { public: @@ -52,17 +54,19 @@ namespace Opm /// \param[in] source Transport source term. /// \param[in] dt Time step. /// \param[in, out] saturation Phase saturations. + /// \param[in, out] surfacevol Surface volume densities for each phase. void solve(const double* darcyflux, const double* pressure, - const double* surfacevol0, const double* porevolume0, const double* porevolume, const double* source, const double dt, - std::vector& saturation); + std::vector& saturation, + std::vector& surfacevol); /// Initialise quantities needed by gravity solver. - void initGravity(); + /// \param[in] grav Gravity vector + void initGravity(const double* grav); /// Solve for gravity segregation. /// This uses a column-wise nonlinear Gauss-Seidel approach. @@ -72,14 +76,14 @@ namespace Opm /// \param[in] columns Vector of cell-columns. /// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] dt Time step. - /// \param[in] grav Gravity vector. /// \param[in, out] saturation Phase saturations. + /// \param[in, out] surfacevol Surface volume densities for each phase. void solveGravity(const std::vector >& columns, const double* pressure, const double* porevolume0, const double dt, - const double* grav, - std::vector& saturation); + std::vector& saturation, + std::vector& surfacevol); private: virtual void solveSingleCell(const int cell); @@ -88,7 +92,7 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - void initGravityDynamic(const double* grav); + void initGravityDynamic(); private: const UnstructuredGrid& grid_; @@ -110,6 +114,7 @@ namespace Opm std::vector saturation_; // P (= num. phases) per cell std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell // For gravity segregation. + const double* gravity_; std::vector trans_; std::vector density_; std::vector gravflux_; diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index b6451809..2740d065 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -52,8 +52,8 @@ namespace Opm source_(0), dt_(0.0), saturation_(grid.number_of_cells, -1.0), - reorder_iterations_(grid.number_of_cells, 0), fractionalflow_(grid.number_of_cells, -1.0), + reorder_iterations_(grid.number_of_cells, 0), mob_(2*grid.number_of_cells, -1.0) #ifdef EXPERIMENT_GAUSS_SEIDEL , ia_upw_(grid.number_of_cells + 1, -1), @@ -107,6 +107,13 @@ namespace Opm toBothSat(saturation_, saturation); } + + const std::vector& TransportModelTwophase::getReorderIterations() const + { + return reorder_iterations_; + } + + // Residual function r(s) for a single-cell implicit Euler transport // // r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) @@ -645,10 +652,7 @@ namespace Opm toBothSat(saturation_, saturation); } - void TransportModelTwophase::getReorderIterations() - { - return reorder_iterations_; - }; + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 358baaa3..97386c8a 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -74,10 +74,11 @@ namespace Opm const double* porevolume, const double dt, std::vector& saturation); - //// Return reorder iterations - //// - //// \param[out] vector of iteration per cell - const std::vector& getReorderIterations(){return reorder_iterations_;}; + + //// Return the number of iterations used by the reordering solver. + //// \return vector of iteration per cell + const std::vector& getReorderIterations() const; + private: virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); @@ -86,7 +87,7 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - private: + private: const UnstructuredGrid& grid_; const IncompPropertiesInterface& props_; const double* visc_; diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index af088e05..6f4f7556 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -48,39 +48,39 @@ namespace Opm void computeInjectedProduced(const BlackoilPropertiesInterface& props, const std::vector& press, const std::vector& z, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced) + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); - } - std::fill(injected, injected + np, 0.0); - std::fill(produced, produced + np, 0.0); + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); std::vector visc(np); - std::vector mob(np); - for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; - const double* sat = &s[np*c]; - props.relperm(1, sat, &c, &mob[0], 0); + std::vector mob(np); + for (int c = 0; c < num_cells; ++c) { + if (src[c] > 0.0) { + injected[0] += src[c]*dt; + } else if (src[c] < 0.0) { + const double flux = -src[c]*dt; + const double* sat = &s[np*c]; + props.relperm(1, sat, &c, &mob[0], 0); props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); - double totmob = 0.0; - for (int p = 0; p < np; ++p) { - mob[p] /= visc[p]; - totmob += mob[p]; - } - for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; - } - } - } + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + mob[p] /= visc[p]; + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + } + } } @@ -93,11 +93,11 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& press, const std::vector& z, - const std::vector& s, - std::vector& totmob) + const std::vector& s, + std::vector& totmob) { std::vector pmobc; @@ -126,12 +126,12 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob, - std::vector& omega) + const std::vector& s, + std::vector& totmob, + std::vector& omega) { std::vector pmobc; @@ -185,10 +185,10 @@ namespace Opm props.relperm(nc, &s[0], &cells[0], &pmobc[0], dpmobc); - std::transform(pmobc.begin(), pmobc.end(), - mu.begin(), - pmobc.begin(), - std::divides()); + std::transform(pmobc.begin(), pmobc.end(), + mu.begin(), + pmobc.begin(), + std::divides()); } /// Computes the fractional flow for each cell in the cells argument @@ -220,5 +220,35 @@ namespace Opm } } + /// Computes the surface volume densities from saturations by the formula + /// z = A s + /// for a number of data points, where z is the surface volume density, + /// s is the saturation (both as column vectors) and A is the + /// phase-to-component relation matrix. + /// @param[in] n number of data points + /// @param[in] np number of phases, must be 2 or 3 + /// @param[in] A array containing n square matrices of size num_phases^2, + /// in Fortran ordering, typically the output of a call + /// to the matrix() method of a BlackoilProperties* class. + /// @param[in] saturation concatenated saturation values (for all P phases) + /// @param[out] surfacevol concatenated surface-volume values (for all P phases) + void computeSurfacevol(const int n, + const int np, + const double* A, + const double* saturation, + double* surfacevol) + { + // Note: since this is a simple matrix-vector product, it can + // be done by a BLAS call, but then we have to reorder the A + // matrix data. + std::fill(surfacevol, surfacevol + n*np, 0.0); + for (int i = 0; i < n; ++i) { + for (int col = 0; col < np; ++col) { + for (int row = 0; row < np; ++row) { + surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col]; + } + } + } + } } // namespace Opm diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 929d42ec..565acc6e 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -46,11 +46,11 @@ namespace Opm void computeInjectedProduced(const BlackoilPropertiesInterface& props, const std::vector& p, const std::vector& z, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced); + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced); /// @brief Computes total mobility for a set of saturation values. /// @param[in] props rock and fluid properties @@ -60,11 +60,11 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob); + const std::vector& s, + std::vector& totmob); /// @brief Computes total mobility and omega for a set of saturation values. /// @param[in] props rock and fluid properties @@ -75,12 +75,12 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob, - std::vector& omega); + const std::vector& s, + std::vector& totmob, + std::vector& omega); /// @brief Computes phase mobilities for a set of saturation values. @@ -96,7 +96,7 @@ namespace Opm const std::vector& z, const std::vector& s, std::vector& pmobc); - + /// Computes the fractional flow for each cell in the cells argument /// @param[in] props rock and fluid properties @@ -112,6 +112,25 @@ namespace Opm const std::vector& s, std::vector& fractional_flows); + + /// Computes the surface volume densities from saturations by the formula + /// z = A s + /// for a number of data points, where z is the surface volume density, + /// s is the saturation (both as column vectors) and A is the + /// phase-to-component relation matrix. + /// @param[in] n number of data points + /// @param[in] np number of phases, must be 2 or 3 + /// @param[in] A array containing n square matrices of size num_phases, + /// in Fortran ordering, typically the output of a call + /// to the matrix() method of a BlackoilProperties* class. + /// @param[in] saturation concatenated saturation values (for all P phases) + /// @param[out] surfacevol concatenated surface-volume values (for all P phases) + void computeSurfacevol(const int n, + const int np, + const double* A, + const double* saturation, + double* surfacevol); + } // namespace Opm #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED