diff --git a/Makefile.am b/Makefile.am index 1ffe69a8..c1c85dec 100644 --- a/Makefile.am +++ b/Makefile.am @@ -92,9 +92,9 @@ 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/SimulatorIncompTwophase.cpp \ opm/core/simulator/SimulatorReport.cpp \ opm/core/simulator/SimulatorTimer.cpp \ -opm/core/simulator/SimulatorTwophase.cpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \ opm/core/transport/reorder/TransportModelInterface.cpp \ opm/core/transport/reorder/TransportModelTwophase.cpp \ @@ -196,8 +196,8 @@ opm/core/pressure/tpfa/ifs_tpfa.h \ opm/core/pressure/tpfa/trans_tpfa.h \ opm/core/simulator/BlackoilState.hpp \ opm/core/simulator/SimulatorReport.hpp \ +opm/core/simulator/SimulatorIncompTwophase.hpp \ opm/core/simulator/SimulatorTimer.hpp \ -opm/core/simulator/SimulatorTwophase.hpp \ opm/core/simulator/TwophaseState.hpp \ opm/core/simulator/WellState.hpp \ opm/core/transport/CSRMatrixBlockAssembler.hpp \ diff --git a/README b/README new file mode 100644 index 00000000..21eb98a2 --- /dev/null +++ b/README @@ -0,0 +1,41 @@ +These are the release notes for opm-core + +Open Porous Media Core Library + +CONTENT +opm-core is the core library within OPM and contains the following + +* Eclipse 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) +* 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. + +DOWNLOAD opm-core +git clone git://github.com/OPM/opm-core.git + +BUILDING opm-core +cd ../opm-core + autoreconf -i + ./configure + make + sudo make install + + +DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS + + +DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS +blas libblas3 lapack liblapack3 libboost libxml2 gcc automake autoconf git +doxygen umfpack 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_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp index c3f7136a..0a7c3a71 100644 --- a/examples/sim_2p_incomp_reorder.cpp +++ b/examples/sim_2p_incomp_reorder.cpp @@ -43,7 +43,7 @@ #include #include -#include +#include #include #include @@ -173,23 +173,22 @@ main(int argc, char** argv) 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"); + 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"); } @@ -200,15 +199,16 @@ main(int argc, char** argv) SimulatorReport rep; if (!use_deck) { // Simple simulation without a deck. - SimulatorTwophase simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - 0, // wells - src, - bcs.c_bcs(), - linsolver, - grav); + WellsManager wells; // no wells. + SimulatorIncompTwophase 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); @@ -255,37 +255,35 @@ main(int argc, char** argv) } // Create and run simulator. - SimulatorTwophase simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells.c_wells(), - src, - bcs.c_bcs(), - linsolver, - grav); + SimulatorIncompTwophase 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); + if (output) { + epoch_rep.reportParam(epoch_os); } // Update total timing report and remember step number. rep += epoch_rep; step = simtimer.currentStepNum(); } } - - epoch_os.close(); + 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); - tot_os.close(); } - + } diff --git a/opm/core/grid/cpgpreprocess/preprocess.c b/opm/core/grid/cpgpreprocess/preprocess.c index 899ccfe7..84e9e8fe 100644 --- a/opm/core/grid/cpgpreprocess/preprocess.c +++ b/opm/core/grid/cpgpreprocess/preprocess.c @@ -587,8 +587,8 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum, z1 = sign*zcorn[i+2*nx*(j+2*ny*(k))]; z2 = sign*zcorn[i+2*nx*(j+2*ny*(k+1))]; - c1 = i/2 + nx*(j/2 + ny*k/2); - c2 = i/2 + nx*(j/2 + ny*(k+1)/2); + c1 = i/2 + nx*(j/2 + ny*(k/2)); + c2 = i/2 + nx*(j/2 + ny*((k+1)/2)); if (((actnum == NULL) || (actnum[c1] && actnum[c2])) diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp similarity index 72% rename from opm/core/simulator/SimulatorTwophase.cpp rename to opm/core/simulator/SimulatorIncompTwophase.cpp index a2a0ad70..35b3e84b 100644 --- a/opm/core/simulator/SimulatorTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -21,7 +21,7 @@ #include "config.h" #endif // HAVE_CONFIG_H -#include +#include #include #include @@ -37,6 +37,8 @@ #include #include +#include + #include #include @@ -56,14 +58,14 @@ namespace Opm { - class SimulatorTwophase::Impl + class SimulatorIncompTwophase::Impl { public: Impl(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropertiesInterface& props, const RockCompressibility* rock_comp, - const Wells* wells, + WellsManager& wells_manager, const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, @@ -80,6 +82,9 @@ namespace Opm 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_; @@ -87,11 +92,10 @@ namespace Opm const UnstructuredGrid& grid_; const IncompPropertiesInterface& props_; const RockCompressibility* rock_comp_; + WellsManager& wells_manager_; const Wells* wells_; const std::vector& src_; - //const FlowBoundaryConditions* bcs_; - //const LinearSolverInterface& linsolver_; - //const double* gravity_; + const FlowBoundaryConditions* bcs_; // Solvers IncompTpfa psolver_; TransportModelTwophase tsolver_; @@ -104,26 +108,26 @@ namespace Opm - SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp, - const Wells* wells, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, - const double* gravity) + SimulatorIncompTwophase::SimulatorIncompTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& 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, src, bcs, linsolver, gravity)); + pimpl_.reset(new Impl(param, grid, props, rock_comp, wells_manager, src, bcs, linsolver, gravity)); } - SimulatorReport SimulatorTwophase::run(SimulatorTimer& timer, - TwophaseState& state, - WellState& well_state) + SimulatorReport SimulatorIncompTwophase::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) { return pimpl_->run(timer, state, well_state); } @@ -293,28 +297,59 @@ namespace Opm } - SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp, - const Wells* wells, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, - const double* gravity) + static bool allNeumannBCs(const FlowBoundaryConditions* bcs) + { + if (bcs == NULL) { + return true; + } else { + return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE) + == bcs->type + bcs->nbc; + } + } + + + static bool allRateWells(const Wells* wells) + { + if (wells == NULL) { + return true; + } + const int nw = wells->number_of_wells; + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells->ctrls[w]; + if (wc->current >= 0) { + if (wc->type[wc->current] == BHP) { + return false; + } + } + } + return true; + } + + + + + + SimulatorIncompTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& 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_(wells), + wells_manager_(wells_manager), + wells_(wells_manager.c_wells()), src_(src), - //bcs_(bcs), - //linsolver_(linsolver), - //gravity_(gravity), + 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, src, bcs), + gravity, wells_manager.c_wells(), src, bcs), tsolver_(grid, props, param.getDefault("nl_tolerance", 1e-9), param.getDefault("nl_maxiter", 30)) @@ -335,6 +370,10 @@ namespace Opm 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); @@ -354,9 +393,9 @@ namespace Opm - SimulatorReport SimulatorTwophase::Impl::run(SimulatorTimer& timer, - TwophaseState& state, - WellState& well_state) + SimulatorReport SimulatorIncompTwophase::Impl::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) { std::vector transport_src; @@ -397,9 +436,9 @@ namespace Opm wellreport.push(props_, *wells_, 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); + 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. @@ -414,18 +453,77 @@ namespace Opm tsolver_.getReorderIterations(), timer.currentStepNum(), output_dir_); } + SimulatorReport sreport; - - // Solve pressure. + + // Solve pressure equation. + if (check_well_controls_) { + computeFractionalFlow(props_, allcells_, 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 rock is 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 ((rock_comp_ == NULL || !rock_comp_->isActive()) + && allNeumannBCs(bcs_) && allRateWells(wells_)) { + // 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; - } while (false); + + // 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()) { @@ -476,11 +574,9 @@ namespace Opm injected, produced, init_satvol); sreport.total_time = step_timer.secsSinceStart(); - if(output_){ - sreport.reportParam(tstep_os); + if (output_) { + sreport.reportParam(tstep_os); } - - } if (output_) { @@ -503,7 +599,7 @@ namespace Opm SimulatorReport report; report.pressure_time = ptime; report.transport_time = ttime; - report.total_time = total_timer.secsSinceStart(); + report.total_time = total_timer.secsSinceStart(); return report; } diff --git a/opm/core/simulator/SimulatorTwophase.hpp b/opm/core/simulator/SimulatorIncompTwophase.hpp similarity index 72% rename from opm/core/simulator/SimulatorTwophase.hpp rename to opm/core/simulator/SimulatorIncompTwophase.hpp index 818d40db..ff4fbb68 100644 --- a/opm/core/simulator/SimulatorTwophase.hpp +++ b/opm/core/simulator/SimulatorIncompTwophase.hpp @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED -#define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED +#ifndef OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED +#define OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED #include #include @@ -32,6 +32,7 @@ namespace Opm namespace parameter { class ParameterGroup; } class IncompPropertiesInterface; class RockCompressibility; + class WellsManager; class LinearSolverInterface; class SimulatorTimer; class TwophaseState; @@ -39,7 +40,7 @@ namespace Opm struct SimulatorReport; /// Class collecting all necessary components for a two-phase simulation. - class SimulatorTwophase + class SimulatorIncompTwophase { public: /// Initialise from parameters and objects to observe. @@ -58,23 +59,23 @@ namespace Opm /// 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] wells if non-null, wells data structure - /// \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 - SimulatorTwophase(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp, - const Wells* wells, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, - const double* gravity); + /// \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 + SimulatorIncompTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& 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