diff --git a/README b/README new file mode 100644 index 00000000..54dbb531 --- /dev/null +++ b/README @@ -0,0 +1,73 @@ +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 deck input and preprosessing +* Fluid properties (basic PVT models and rock properties) +* 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) + + +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. + + +DOWNLOADING +----------- + +git clone git://github.com/OPM/opm-core.git + + +BUILDING +-------- + + cd ../opm-core + autoreconf -i + ./configure + make + sudo make install + + +DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS +------------------------------------------- + +(to be updated) + + +DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS +----------------------------------------- + +blas libblas3 lapack liblapack3 libboost libxml2 gcc automake autoconf +git doxygen umfpack + + +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/configure.ac b/configure.ac index 59d642f0..cc4b92e1 100644 --- a/configure.ac +++ b/configure.ac @@ -8,6 +8,10 @@ AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) +# Needed for automake since version 1.12 because extra-portability +# warnings were then added to -Wall. Ifdef makes it backwards compatible. +m4_ifdef([AM_PROG_AR], [AM_PROG_AR]) + AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([opm/core/grid.h]) AC_CONFIG_HEADERS([config.h]) diff --git a/examples/sim_2p_comp_reorder.cpp b/examples/sim_2p_comp_reorder.cpp index 22e32c37..9f2b516b 100644 --- a/examples/sim_2p_comp_reorder.cpp +++ b/examples/sim_2p_comp_reorder.cpp @@ -277,9 +277,9 @@ main(int argc, char** argv) 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); + 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/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 02f67556..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 @@ -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 67d3a651..954e07fa 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -504,13 +504,13 @@ namespace Opm cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; int was_adjusted = 0; - if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { - was_adjusted = + 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 = + 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], diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index affada84..0b54178e 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -61,7 +61,7 @@ namespace Opm /// and completions does not change during the /// run. However, controls (only) are allowed /// to change. - CompressibleTpfa(const UnstructuredGrid& grid, + CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, @@ -71,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 @@ -111,13 +111,13 @@ namespace Opm void solveIncrement(); double residualNorm() const; double incrementNorm() const; - void computeResults(BlackoilState& state, + void computeResults(BlackoilState& state, WellState& well_state) const; protected: void computeWellPotentials(const BlackoilState& state); // ------ Data that will remain unmodified after construction. ------ - const UnstructuredGrid& grid_; + const UnstructuredGrid& grid_; const BlackoilPropertiesInterface& props_; const RockCompressibility* rock_comp_props_; const LinearSolverInterface& linsolver_; @@ -126,12 +126,12 @@ namespace Opm 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 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_; diff --git a/opm/core/simulator/SimulatorIncompTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp index 9de04f38..1f7e1119 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -77,7 +77,6 @@ namespace Opm private: // Data. - // Parameters for output. bool output_; bool output_vtk_; @@ -133,7 +132,38 @@ namespace Opm return pimpl_->run(timer, state, well_state); } - + 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]) + { + std::cout.precision(5); + const int width = 18; + os << "\nVolume balance report (all numbers relative to total pore volume).\n"; + os << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + os << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + os << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + os << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + os << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + os << " 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; + os << " 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; + os.precision(8); + } static void outputStateVtk(const UnstructuredGrid& grid, const Opm::TwophaseState& state, @@ -145,10 +175,10 @@ namespace Opm vtkfilename << output_dir << "/vtk_files"; boost::filesystem::path fpath(vtkfilename.str()); try { - create_directories(fpath); + create_directories(fpath); } catch (...) { - THROW("Creating directories failed: " << fpath); + THROW("Creating directories failed: " << fpath); } vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); @@ -164,6 +194,27 @@ namespace Opm Opm::writeVtkData(grid, dm, vtkfile); } + 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, @@ -183,10 +234,10 @@ namespace Opm fname << output_dir << "/" << it->first; boost::filesystem::path fpath = fname.str(); try { - create_directories(fpath); + create_directories(fpath); } catch (...) { - THROW("Creating directories failed: " << fpath); + THROW("Creating directories failed: " << fpath); } fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; std::ofstream file(fname.str().c_str()); @@ -377,6 +428,9 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputVectorMatlab(std::string("reorder_it"), + tsolver_.getReorderIterations(), + timer.currentStepNum(), output_dir_); } SimulatorReport sreport; @@ -474,6 +528,14 @@ namespace Opm if (use_segregation_split_) { tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation()); } + 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.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } } transport_timer.stop(); double tt = transport_timer.secsSinceStart(); @@ -486,41 +548,10 @@ namespace Opm 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.saturation(), - timer.currentTime() + timer.currentStepLength(), - well_state.bhp(), well_state.perfRates()); - } + 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); @@ -532,6 +563,9 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputVectorMatlab(std::string("reorder_it"), + tsolver_.getReorderIterations(), + timer.currentStepNum(), output_dir_); outputWaterCut(watercut, output_dir_); if (wells_) { outputWellReport(wellreport, output_dir_); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index c2d50327..697f94a5 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -96,7 +96,7 @@ 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 non-miscibility requirement (only done for first cell). + // 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."); } diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index b3630ea7..2740d065 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -53,6 +53,7 @@ namespace Opm dt_(0.0), saturation_(grid.number_of_cells, -1.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), @@ -101,11 +102,18 @@ namespace Opm &seq[0], &comp[0], &ncomp, &ia_downw_[0], &ja_downw_[0]); #endif - + std::fill(reorder_iterations_.begin(),reorder_iterations_.end(),0); reorderAndTransport(grid_, darcyflux); 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) ) @@ -170,9 +178,11 @@ namespace Opm // if (std::fabs(r0) < tol_) { // return; // } - int iters_used; + int iters_used = 0; // saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); + // add if it is iteration on an out loop + reorder_iterations_[cell] = reorder_iterations_[cell] + iters_used; fractionalflow_[cell] = fracFlow(saturation_[cell], cell); } @@ -544,8 +554,9 @@ namespace Opm const int cell = cells[pos]; GravityResidual res(*this, cells, pos, gravflux); if (std::fabs(res(saturation_[cell])) > tol_) { - int iters_used; + int iters_used = 0; saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); + reorder_iterations_[cell] = reorder_iterations_[cell] + iters_used; } saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]); mobility(saturation_[cell], cell, &mob_[2*cell]); diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 74ed2bfe..97386c8a 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -23,7 +23,7 @@ #include #include #include - +#include struct UnstructuredGrid; namespace Opm @@ -75,6 +75,10 @@ namespace Opm const double dt, std::vector& saturation); + //// 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); @@ -83,7 +87,6 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - private: const UnstructuredGrid& grid_; const IncompPropertiesInterface& props_; @@ -99,6 +102,8 @@ namespace Opm double dt_; std::vector saturation_; // one per cell, only water saturation! std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell + std::vector reorder_iterations_; + //std::vector reorder_fval_; // For gravity segregation. std::vector gravflux_; std::vector mob_; diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 548aaa23..6f4f7556 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -243,9 +243,9 @@ namespace Opm // matrix data. std::fill(surfacevol, surfacevol + n*np, 0.0); for (int i = 0; i < n; ++i) { - for (int row = 0; row < np; ++row) { - for (int col = 0; col < np; ++col) { - surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col]; + 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]; } } }