From 742303534a4488128fac9ba626f9bfbe0269b76c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 May 2012 10:53:50 +0200 Subject: [PATCH 01/31] Added preliminary sketch of compressible tpfa solver. This solver is: - using the residual based assembler, - aiming to include the nonlinear Newton iterations (therefore also the (re-)evaluation of fluid data). --- Makefile.am | 2 + opm/core/pressure/CompressibleTpfa.cpp | 189 +++++++++++++++++++++++++ opm/core/pressure/CompressibleTpfa.hpp | 102 +++++++++++++ 3 files changed, 293 insertions(+) create mode 100644 opm/core/pressure/CompressibleTpfa.cpp create mode 100644 opm/core/pressure/CompressibleTpfa.hpp diff --git a/Makefile.am b/Makefile.am index ccf5503e..e1a5e3c1 100644 --- a/Makefile.am +++ b/Makefile.am @@ -89,6 +89,7 @@ opm/core/pressure/well.c \ opm/core/pressure/fsh_common_impl.c \ opm/core/pressure/fsh.c \ opm/core/pressure/IncompTpfa.cpp \ +opm/core/pressure/CompressibleTpfa.cpp \ opm/core/pressure/tpfa/ifs_tpfa.c \ opm/core/pressure/tpfa/compr_source.c \ opm/core/pressure/tpfa/cfs_tpfa.c \ @@ -201,6 +202,7 @@ opm/core/WellsManager.hpp \ opm/core/pressure/fsh.h \ opm/core/pressure/HybridPressureSolver.hpp \ opm/core/pressure/IncompTpfa.hpp \ +opm/core/pressure/CompressibleTpfa.hpp \ opm/core/pressure/TPFAPressureSolver.hpp \ opm/core/pressure/tpfa/compr_quant.h \ opm/core/pressure/tpfa/compr_source.h \ diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp new file mode 100644 index 00000000..cdc4db54 --- /dev/null +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -0,0 +1,189 @@ +/* + 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 . +*/ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace Opm +{ + + + /// Construct solver. + /// \param[in] g A 2d or 3d grid. + /// \param[in] permeability Array of permeability tensors, the array + /// should have size N*D^2, if D == g.dimensions + /// and N == g.number_of_cells. + /// \param[in] gravity Gravity vector. If nonzero, the array should + /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL + CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& g, + const double* permeability, + const double* /* gravity */, // ??? + const LinearSolverInterface& linsolver, + const struct Wells* wells, + const int num_phases) + : grid_(g), + linsolver_(linsolver), + htrans_(g.cell_facepos[ g.number_of_cells ]), + trans_ (g.number_of_faces), + wells_(wells) + { + if (wells_ && (wells_->number_of_phases != num_phases)) { + THROW("Inconsistent number of phases specified: " + << wells_->number_of_phases << " != " << num_phases); + } + const int num_dofs = g.number_of_cells + wells ? wells->number_of_wells : 0; + pressure_increment_.resize(num_dofs); + UnstructuredGrid* gg = const_cast(&grid_); + tpfa_htrans_compute(gg, permeability, &htrans_[0]); + tpfa_trans_compute(gg, &htrans_[0], &trans_[0]); + cfs_tpfa_res_wells w; + w.W = const_cast(wells_); + w.data = NULL; + h_ = cfs_tpfa_res_construct(gg, &w, num_phases); + } + + + + + /// Destructor. + CompressibleTpfa::~CompressibleTpfa() + { + cfs_tpfa_res_destroy(h_); + } + + + + + /// Solve pressure equation, by Newton iterations. + void CompressibleTpfa::solve() + { + // Set up dynamic data. + computeDynamicData(); + + // Assemble J and F. + assemble(); + + bool residual_ok = false; // Replace with tolerance check. + while (!residual_ok) { + // Solve for increment in Newton method: + // incr = x_{n+1} - x_{n} = -J^{-1}F + // (J is Jacobian matrix, F is residual) + solveIncrement(); + + // Update pressure vars with increment. + + // Set up dynamic data. + computeDynamicData(); + + // Assemble J and F. + assemble(); + + // Check for convergence. + // Include both tolerance check for residual + // and solution change. + } + + // Write to output parameters. + // computeResults(...); + } + + + + + /// Solve pressure equation, by Newton iterations. + void CompressibleTpfa::computeDynamicData() + { + } + + + + + /// Compute the residual and Jacobian. + void CompressibleTpfa::assemble() + { + // Arguments or members? + const double dt = 0.0; + const double* z = NULL; + const double* cell_press = NULL; + const double* well_bhp = NULL; + const double* porevol = NULL; + + UnstructuredGrid* gg = const_cast(&grid_); + CompletionData completion_data; + completion_data.gpot = 0; // TODO + completion_data.A = 0; // TODO + completion_data.phasemob = 0; // TODO + cfs_tpfa_res_wells wells_tmp; + wells_tmp.W = const_cast(wells_); + wells_tmp.data = &completion_data; + cfs_tpfa_res_forces forces; + forces.wells = &wells_tmp; + forces.src = NULL; // Check if it is legal to leave it as NULL. + compr_quantities_gen cq; + cq.Ac = 0; // TODO + cq.dAc = 0; // TODO + cq.Af = 0; // TODO + cq.phasemobf = 0; // TODO + cq.voldiscr = 0; // TODO + // TODO: gravcapf_ must be set already. + cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &gravcapf_[0], &cell_press[0], &well_bhp[0], + &porevol[0], h_); + } + + + + + /// Computes pressure_increment_. + void CompressibleTpfa::solveIncrement() + { + // Increment is equal to -J^{-1}F + linsolver_.solve(h_->J, h_->F, &pressure_increment_[0]); + std::transform(pressure_increment_.begin(), pressure_increment_.end(), + pressure_increment_.begin(), std::negate()); + } + + + + + /// Compute the output. + void CompressibleTpfa::computeResults(std::vector& // pressure + , + std::vector& // faceflux + , + std::vector& // well_bhp + , + std::vector& // well_rate + ) + { + } + +} // namespace Opm diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp new file mode 100644 index 00000000..fdacf426 --- /dev/null +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -0,0 +1,102 @@ +/* + 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_COMPRESSIBLETPFA_HEADER_INCLUDED +#define OPM_COMPRESSIBLETPFA_HEADER_INCLUDED + + +#include + +struct UnstructuredGrid; +struct cfs_tpfa_res_data; +struct Wells; +struct FlowBoundaryConditions; + +namespace Opm +{ + + class LinearSolverInterface; + + /// Encapsulating a tpfa pressure solver for the compressible-fluid case. + /// Supports gravity, wells and simple sources as driving forces. + /// Below we use the shortcuts D for the number of dimensions, N + /// for the number of cells and F for the number of faces. + class CompressibleTpfa + { + public: + /// Construct solver. + /// \param[in] g A 2d or 3d grid. + /// \param[in] permeability Array of permeability tensors, the array + /// should have size N*D^2, if D == g.dimensions + /// and N == g.number_of_cells. + /// \param[in] gravity Gravity vector. If nonzero, the array should + /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL + /// \param[in] num_phases Must be 2 or 3. + CompressibleTpfa(const UnstructuredGrid& g, + const double* permeability, + const double* gravity, + const LinearSolverInterface& linsolver, + const struct Wells* wells, + const int num_phases); + + /// Destructor. + ~CompressibleTpfa(); + + void solve(); + + private: + void computeDynamicData(); + void assemble(); + void solveIncrement(); + + void computeResults(std::vector& pressure, + std::vector& faceflux, + std::vector& well_bhp, + std::vector& well_rate); + + // ------ Data that will remain unmodified after construction. ------ + const UnstructuredGrid& grid_; + const LinearSolverInterface& linsolver_; + std::vector htrans_; + std::vector trans_ ; + const Wells* wells_; // Outside may modify controls (only) between calls to solve(). + + // ------ Internal data for the cfs_tpfa_res solver. ------ + struct cfs_tpfa_res_data* h_; + + // ------ Data that will be modified for every solve. ------ + + // ------ Data that will be modified for every solver iteration. ------ + // The update to be applied to the pressures (cell and bhp). + std::vector pressure_increment_; + + + + // Gravity and capillary contributions (per face). + std::vector gravcapf_; + + + }; + +} // namespace Opm + + +#endif // OPM_COMPRESSIBLETPFA_HEADER_INCLUDED From 3c3ce528505060b1cec40f91cbf48f0c3babc3ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 May 2012 10:55:09 +0200 Subject: [PATCH 02/31] Added (non-compiling) test program for compressible fluid case. --- examples/Makefile.am | 7 +- examples/sim_wateroil.cpp | 552 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 557 insertions(+), 2 deletions(-) create mode 100644 examples/sim_wateroil.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index ae2650ee..703e4bd3 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -7,10 +7,13 @@ LDADD = $(top_builddir)/libopmcore.la noinst_PROGRAMS = \ scaneclipsedeck \ spu_2p \ -wells_example - +wells_example \ +sim_wateroil spu_2p_SOURCES = spu_2p.cpp spu_2p_LDADD = $(LDADD) $(LIBS) $(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) $(LAPACK_LIBS) $(LIBS) $(LIBS) +sim_wateroil_SOURCES = sim_wateroil.cpp +sim_wateroil_LDADD = $(LDADD) $(LIBS) $(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) $(LAPACK_LIBS) $(LIBS) $(LIBS) + wells_example_SOURCES = wells_example.cpp diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp new file mode 100644 index 00000000..96fa70c1 --- /dev/null +++ b/examples/sim_wateroil.cpp @@ -0,0 +1,552 @@ +/* + 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 +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define COMPR_INIT_FIXED 0 +#define PRESSURE_SOLVER_FIXED 0 +#define TRANSPORT_SOLVER_FIXED 0 + + + +static void outputState(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) +{ + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/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); + + // 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 << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + 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); +} + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + std::cout << "\n================ Test program for weakly compressible two-phase flow ===============\n\n"; + Opm::parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // Reading various control parameters. + const bool guess_old_solution = param.getDefault("guess_old_solution", false); + const bool use_reorder = param.getDefault("use_reorder", true); + const bool output = param.getDefault("output", true); + std::string output_dir; + int output_interval = 1; + if (output) { + 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", output_interval); + } + const int num_transport_substeps = param.getDefault("num_transport_substeps", 1); + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr wells; + boost::scoped_ptr rock_comp; + Opm::SimulatorTimer simtimer; + Opm::TwophaseState 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"); + Opm::EclipseGridParser deck(deck_filename); + // Grid init + grid.reset(new Opm::GridManager(deck)); + // Rock and fluid init + const int* gc = grid->c_grid()->global_cell; + std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); + props.reset(new Opm::BlackoilPropertiesFromDeck(deck, global_cell)); + // Wells init. + wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); + check_well_controls = param.getDefault("check_well_controls", false); + max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Timer init. + if (deck.hasField("TSTEP")) { + simtimer.init(deck); + } else { + simtimer.init(param); + } + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(deck)); + // Gravity. + gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; + // Init state variables (saturation and pressure). +#if COMPR_INIT_FIXED + if (param.has("init_saturation")) { + initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + } +#endif // COMPR_INIT_FIXED + } 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 Opm::GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new Opm::BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Wells init. + wells.reset(new Opm::WellsManager()); + // Timer init. + simtimer.init(param); + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). +#if COMPR_INIT_FIXED + initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); +#endif // COMPR_INIT_FIXED + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + if (use_gravity) { + if (props->density()[0] == props->density()[1]) { + std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; + } + } + bool use_segregation_split = false; + bool use_column_solver = false; + bool use_gauss_seidel_gravity = false; + if (use_gravity && use_reorder) { + 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. + int num_cells = grid->c_grid()->number_of_cells; + std::vector totmob; + std::vector omega; // Will remain empty if no gravity. + std::vector rc; // Will remain empty if no rock compressibility. + + // Extra rock init. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), *props, porevol); + } + double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + + // We need a separate reorder_sat, because the reorder + // code expects a scalar sw, not both sw and so. + std::vector reorder_sat(num_cells); + std::vector src(num_cells, 0.0); + + // Initialising src + if (wells->c_wells()) { + // Do nothing, wells will be the driving force, not source terms. + // Opm::wellsToSrc(*wells->c_wells(), num_cells, src); + } else { + 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/Opm::unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + std::vector reorder_src = src; + + // Solvers init. + // Linear solver. + Opm::LinearSolverFactory linsolver(param); + // Pressure solver. + const double *grav = use_gravity ? &gravity[0] : 0; + Opm::CompressibleTpfa psolver(*grid->c_grid(), props->permeability(), grav, + linsolver, wells->c_wells(), props->numPhases()); + // Reordering solver. + const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); + const int nl_maxiter = param.getDefault("nl_maxiter", 30); +#if TRANSPORT_SOLVER_FIXED + Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); + if (use_gauss_seidel_gravity) { + reorder_model.initGravity(grav); + } +#endif // TRANSPORT_SOLVER_FIXED + // Column-based gravity segregation solver. + std::vector > columns; + if (use_column_solver) { + Opm::extractColumn(*grid->c_grid(), columns); + } + + // The allcells vector is used in calls to computeTotalMobility() + // and computeTotalMobilityOmega(). + std::vector allcells(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells[cell] = cell; + } + + // Warn if any parameters are unused. + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + + // Write parameters used for later reference. + if (output) { + param.writeParam(output_dir + "/spu_2p.param"); + } + + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl; + 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 well_bhp; + std::vector well_perfrates; + std::vector fractional_flows; + std::vector well_resflows_phase; + int num_wells = 0; + if (wells->c_wells()) { + num_wells = wells->c_wells()->number_of_wells; + well_bhp.resize(num_wells, 0.0); + well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0); + well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); + wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates); + } + for (; !simtimer.done(); ++simtimer) { + // Report timestep and (optionally) write state to disk. + simtimer.report(std::cout); + if (output && (simtimer.currentStepNum() % output_interval == 0)) { + outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); + } + + // Solve pressure. + if (use_gravity) { + computeTotalMobilityOmega(*props, allcells, state.saturation(), totmob, omega); + } else { + computeTotalMobility(*props, allcells, state.saturation(), totmob); + } + std::vector wdp; + if (wells->c_wells()) { + Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp); + } + if (check_well_controls) { + computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows); + } + if (check_well_controls) { + wells->applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } + bool well_control_passed = !check_well_controls; + int well_control_iteration = 0; + do { // Well control outer loop. + pressure_timer.start(); +#if PRESSURE_SOLVER_FIXED + if (rock_comp->isActive()) { + rc.resize(num_cells); + std::vector initial_pressure = state.pressure(); + std::vector initial_porevolume(num_cells); + computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume); + std::vector pressure_increment(num_cells + num_wells); + std::vector prev_pressure(num_cells + num_wells); + for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { + + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp->rockComp(state.pressure()[cell]); + } + computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); + std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells); + // prev_pressure = state.pressure(); + + // compute pressure increment + psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, + prev_pressure, initial_porevolume, simtimer.currentStepLength(), + pressure_increment); + + double max_change = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += pressure_increment[cell]; + max_change = std::max(max_change, std::fabs(pressure_increment[cell])); + } + for (int well = 0; well < num_wells; ++well) { + well_bhp[well] += pressure_increment[num_cells + well]; + max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well])); + } + + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance) { + break; + } + } + psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), + well_bhp, well_perfrates); + } else { + psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), + well_bhp, well_perfrates); + } +#endif // PRESSURE_SOLVER_FIXED + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + + + if (check_well_controls) { + Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), + fractional_flows, + well_perfrates, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells->conditionsMet(well_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); + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, + wells->c_wells(), well_perfrates, reorder_src); + + // Solve transport. + transport_timer.start(); +#if TRANSPORT_SOLVER_FIXED + double stepsize = simtimer.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) { + Opm::toWaterSat(state.saturation(), reorder_sat); + reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0], + stepsize, &reorder_sat[0]); + Opm::toBothSat(reorder_sat, state.saturation()); + Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); + if (use_segregation_split) { + reorder_model.solveGravity(columns, &porevol[0], stepsize, reorder_sat); + Opm::toBothSat(reorder_sat, state.saturation()); + } + } +#endif // TRANSPORT_SOLVER_FIXED + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + 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(simtimer.currentTime() + simtimer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells->c_wells()) { + wellreport.push(*props, *wells->c_wells(), state.saturation(), + simtimer.currentTime() + simtimer.currentStepLength(), + well_bhp, well_perfrates); + } + } + total_timer.stop(); + + std::cout << "\n\n================ End of simulation ===============\n" + << "Total time taken: " << total_timer.secsSinceStart() + << "\n Pressure time: " << ptime + << "\n Transport time: " << ttime << std::endl; + + if (output) { + outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); + outputWaterCut(watercut, output_dir); + if (wells->c_wells()) { + outputWellReport(wellreport, output_dir); + } + } +} From 26fc6a335a9b439bdd404cca0d70692254e2bd0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 May 2012 21:53:36 +0200 Subject: [PATCH 03/31] Adapt to changed computePorevolume() interface. --- examples/sim_wateroil.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 96fa70c1..f37d56a4 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -136,6 +136,8 @@ static void outputWellReport(const Opm::WellReport& wellreport, int main(int argc, char** argv) { + using namespace Opm; + std::cout << "\n================ Test program for weakly compressible two-phase flow ===============\n\n"; Opm::parameter::ParameterGroup param(argc, argv, false); std::cout << "--------------- Reading parameters ---------------" << std::endl; @@ -230,7 +232,7 @@ main(int argc, char** argv) // Warn if gravity but no density difference. bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); if (use_gravity) { - if (props->density()[0] == props->density()[1]) { + if (props->surfaceDensity()[0] == props->surfaceDensity()[1]) { std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; } } @@ -268,9 +270,9 @@ main(int argc, char** argv) // Extra rock init. std::vector porevol; if (rock_comp->isActive()) { - computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); } else { - computePorevolume(*grid->c_grid(), *props, porevol); + computePorevolume(*grid->c_grid(), props->porosity(), porevol); } double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); @@ -398,7 +400,7 @@ main(int argc, char** argv) rc.resize(num_cells); std::vector initial_pressure = state.pressure(); std::vector initial_porevolume(num_cells); - computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume); + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, initial_pressure, initial_porevolume); std::vector pressure_increment(num_cells + num_wells); std::vector prev_pressure(num_cells + num_wells); for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { @@ -406,7 +408,7 @@ main(int argc, char** argv) for (int cell = 0; cell < num_cells; ++cell) { rc[cell] = rock_comp->rockComp(state.pressure()[cell]); } - computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells); // prev_pressure = state.pressure(); From c5ef2eb975317d1ca2ad3522941bd0b1c0e7aff0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 11:54:59 +0200 Subject: [PATCH 04/31] WellsManager(): Support grids that do not define "global_cell" Assume that the "global_cell" in that case corresponds to the identity mapping. --- opm/core/WellsManager.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index 72354796..6fe66a1e 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -250,8 +250,16 @@ namespace Opm const int* global_cell = grid.global_cell; const int* cpgdim = grid.cartdims; std::map cartesian_to_compressed; - for (int i = 0; i < grid.number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + + if (global_cell) { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } } // Get COMPDAT data From df5dbb35577f7a8a229ab2d4b700dec20ec82396 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 12:07:04 +0200 Subject: [PATCH 05/31] assignPermeability(): Don't build an Inspector to count global cells. Thus, assingPermeability() is applicable to simulation decks that aren't based on corner-point descriptions (e.g., the (DXV,DYV,DZV) grid type). --- opm/core/fluid/RockFromDeck.cpp | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 29c7a913..4e45b0d1 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -35,6 +35,8 @@ namespace Opm PermeabilityKind fillTensor(const EclipseGridParser& parser, std::vector*>& tensor, std::tr1::array& kmap); + + int numGlobalCells(const EclipseGridParser& parser); } // anonymous namespace @@ -82,10 +84,8 @@ namespace Opm const std::vector& global_cell, double perm_threshold) { - const int dim = 3; - EclipseGridInspector insp(parser); - std::tr1::array dims = insp.gridSize(); - int num_global_cells = dims[0]*dims[1]*dims[2]; + const int dim = 3; + const int num_global_cells = numGlobalCells(parser); ASSERT (num_global_cells > 0); permeability_.assign(dim * dim * global_cell.size(), 0.0); @@ -330,6 +330,26 @@ 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 From 2acc750d93b28bcf65995a1f016488036a7e26cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 12:08:15 +0200 Subject: [PATCH 06/31] Add input support for ECLIPSE tensor grids (DXV &c). --- opm/core/eclipse/EclipseGridParser.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index ed1d5e10..29d8bb05 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -83,7 +83,9 @@ namespace EclipseKeywords string("BULKMOD"), string("YOUNGMOD"), string("LAMEMOD"), string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"), string("MULTPV"), string("PRESSURE"), string("SGAS"), - string("SWAT"), string("SOIL"), string("RS") + string("SWAT"), string("SOIL"), string("RS"), + string("DXV"), string("DYV"), string("DZV"), + string("DEPTHZ") }; const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); @@ -436,7 +438,9 @@ void EclipseGridParser::convertToSI() // Find the right unit. double unit = 1e100; bool do_convert = true; - if (key == "COORD" || key == "ZCORN") { + if (key == "COORD" || key == "ZCORN" || + key == "DXV" || key == "DYV" || key == "DZV" || + key == "DEPTHZ") { unit = units_.length; } else if (key == "PERMX" || key == "PERMY" || key == "PERMZ" || key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" || From 8fd9826467019ff1d7f223f63d8491d6d1abdc97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 May 2012 12:48:16 +0200 Subject: [PATCH 07/31] Added satRange() method to BlackoilPropertiesInterface. --- opm/core/fluid/BlackoilPropertiesBasic.cpp | 15 +++++++++++++++ opm/core/fluid/BlackoilPropertiesBasic.hpp | 13 +++++++++++++ opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 15 +++++++++++++++ opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 13 +++++++++++++ opm/core/fluid/BlackoilPropertiesInterface.hpp | 13 +++++++++++++ 5 files changed, 69 insertions(+) diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp index 25b520d6..be8211c1 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.cpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.cpp @@ -214,6 +214,21 @@ namespace Opm } + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// \param[in] n Number of data points. + /// \param[in] cells Array of n cell indices. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + void BlackoilPropertiesBasic::satRange(const int n, + const int* /*cells*/, + double* smin, + double* smax) const + { + satprops_.satRange(n, smin, smax); + } + } // namespace Opm diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp index 76dc4bd4..d879036e 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.hpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -150,6 +150,19 @@ namespace Opm double* pc, double* dpcds) const; + + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// \param[in] n Number of data points. + /// \param[in] cells Array of n cell indices. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + virtual void satRange(const int n, + const int* cells, + double* smin, + double* smax) const; + private: RockBasic rock_; PvtPropertiesBasic pvt_; diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 5aeb58fc..0eb23cb0 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -258,6 +258,21 @@ namespace Opm } + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// \param[in] n Number of data points. + /// \param[in] cells Array of n cell indices. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + void BlackoilPropertiesFromDeck::satRange(const int n, + const int* cells, + double* smin, + double* smax) const + { + satprops_.satRange(n, cells, smin, smax); + } + } // namespace Opm diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index a6408e06..9c20b568 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -146,6 +146,19 @@ namespace Opm double* pc, double* dpcds) const; + + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// \param[in] n Number of data points. + /// \param[in] cells Array of n cell indices. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + virtual void satRange(const int n, + const int* cells, + double* smin, + double* smax) const; + private: RockFromDeck rock_; BlackoilPvtProperties pvt_; diff --git a/opm/core/fluid/BlackoilPropertiesInterface.hpp b/opm/core/fluid/BlackoilPropertiesInterface.hpp index 393036a2..dd8a857f 100644 --- a/opm/core/fluid/BlackoilPropertiesInterface.hpp +++ b/opm/core/fluid/BlackoilPropertiesInterface.hpp @@ -136,6 +136,19 @@ namespace Opm const int* cells, double* pc, double* dpcds) const = 0; + + + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// \param[in] n Number of data points. + /// \param[in] cells Array of n cell indices. + /// \param[out] smin Array of nP minimum s values, array must be valid before calling. + /// \param[out] smax Array of nP maximum s values, array must be valid before calling. + virtual void satRange(const int n, + const int* cells, + double* smin, + double* smax) const = 0; }; From 762a6083a2005ebe78bc03e77952bf88dca7b28a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 May 2012 12:49:15 +0200 Subject: [PATCH 08/31] Added utility functions for compressible fluid case. --- Makefile.am | 2 + opm/core/utility/miscUtilitiesBlackoil.cpp | 224 +++++++++++++++++++++ opm/core/utility/miscUtilitiesBlackoil.hpp | 117 +++++++++++ 3 files changed, 343 insertions(+) create mode 100644 opm/core/utility/miscUtilitiesBlackoil.cpp create mode 100644 opm/core/utility/miscUtilitiesBlackoil.hpp diff --git a/Makefile.am b/Makefile.am index e1a5e3c1..b481d607 100644 --- a/Makefile.am +++ b/Makefile.am @@ -69,6 +69,7 @@ opm/core/utility/parameters/tinyxml/tinyxml.cpp \ opm/core/utility/parameters/tinyxml/tinyxmlerror.cpp \ opm/core/utility/parameters/tinyxml/tinyxmlparser.cpp \ opm/core/utility/miscUtilities.cpp \ +opm/core/utility/miscUtilitiesBlackoil.cpp \ opm/core/utility/SimulatorTimer.cpp \ opm/core/utility/StopWatch.cpp \ opm/core/utility/newwells.c \ @@ -182,6 +183,7 @@ opm/core/utility/initState_impl.hpp \ opm/core/utility/linInt.hpp \ opm/core/utility/linearInterpolation.hpp \ opm/core/utility/miscUtilities.hpp \ +opm/core/utility/miscUtilitiesBlackoil.hpp \ opm/core/utility/writeVtkData.hpp \ opm/core/linalg/sparse_sys.h \ opm/core/linalg/blas_lapack.h \ diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp new file mode 100644 index 00000000..070cba05 --- /dev/null +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -0,0 +1,224 @@ +/* + 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 . +*/ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + /// @brief Computes injected and produced volumes of all phases. + /// Note 1: assumes that only the first phase is injected. + /// Note 2: assumes that transport has been done with an + /// implicit method, i.e. that the current state + /// gives the mobilities used for the preceding timestep. + /// @param[in] props fluid and rock properties. + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all P phases) + /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. + 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 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); + props.viscosity(1, &p[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; + } + } + } + } + + + + /// @brief Computes total mobility for a set of saturation values. + /// @param[in] props rock and fluid properties + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @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& p, + const std::vector& z, + const std::vector& s, + std::vector& totmob) + { + std::vector pmobc; + + computePhaseMobilities(props, cells, p, z, s, pmobc); + + const std::size_t np = props.numPhases(); + const std::vector::size_type nc = cells.size(); + + totmob.clear(); + totmob.resize(nc, 0.0); + + for (std::vector::size_type c = 0; c < nc; ++c) { + for (std::size_t p = 0; p < np; ++p) { + totmob[ c ] += pmobc[c*np + p]; + } + } + } + + /* + /// @brief Computes total mobility and omega for a set of saturation values. + /// @param[in] props rock and fluid properties + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all phases) + /// @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& p, + const std::vector& z, + const std::vector& s, + std::vector& totmob, + std::vector& omega) + { + std::vector pmobc; + + computePhaseMobilities(props, cells, p, z, s, pmobc); + + const std::size_t np = props.numPhases(); + const std::vector::size_type nc = cells.size(); + + totmob.clear(); + totmob.resize(nc, 0.0); + omega.clear(); + omega.resize(nc, 0.0); + + const double* rho = props.density(); + for (std::vector::size_type c = 0; c < nc; ++c) { + for (std::size_t p = 0; p < np; ++p) { + totmob[ c ] += pmobc[c*np + p]; + omega [ c ] += pmobc[c*np + p] * rho[ p ]; + } + + omega[ c ] /= totmob[ c ]; + } + } + */ + + /// @brief Computes phase mobilities for a set of saturation values. + /// @param[in] props rock and fluid properties + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all phases) + /// @param[out] pmobc phase mobilities (for all phases). + void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props, + const std::vector& cells, + const std::vector& p, + const std::vector& z, + const std::vector& s, + std::vector& pmobc) + { + const int nc = props.numCells(); + const int np = props.numPhases(); + + ASSERT (int(s.size()) == nc * np); + + std::vector mu(nc*np); + props.viscosity(nc, &p[0], &z[0], &cells[0], &mu[0], 0); + + pmobc.clear(); + pmobc.resize(nc*np, 0.0); + double* dpmobc = 0; + props.relperm(nc, &s[0], &cells[0], + &pmobc[0], dpmobc); + + std::transform(pmobc.begin(), pmobc.end(), + mu.begin(), + pmobc.begin(), + std::divides()); + } + + /// Computes the fractional flow for each cell in the cells argument + /// @param[in] props rock and fluid properties + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all phases) + /// @param[out] fractional_flow the fractional flow for each phase for each cell. + void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props, + const std::vector& cells, + const std::vector& p, + const std::vector& z, + const std::vector& s, + std::vector& fractional_flows) + { + const int num_phases = props.numPhases(); + std::vector pc_mobs(cells.size() * num_phases); + computePhaseMobilities(props, cells, p, z, s, pc_mobs); + fractional_flows.resize(cells.size() * num_phases); + for (size_t i = 0; i < cells.size(); ++i) { + double phase_sum = 0.0; + for (int phase = 0; phase < num_phases; ++phase) { + phase_sum += pc_mobs[i * num_phases + phase]; + } + for (int phase = 0; phase < num_phases; ++phase) { + fractional_flows[i * num_phases + phase] = pc_mobs[i * num_phases + phase] / phase_sum; + } + } + } + + +} // namespace Opm diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp new file mode 100644 index 00000000..929d42ec --- /dev/null +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -0,0 +1,117 @@ +/* + 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_MISCUTILITIESBLACKOIL_HEADER_INCLUDED +#define OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED + +#include + +struct UnstructuredGrid; + +namespace Opm +{ + + class BlackoilPropertiesInterface; + + /// @brief Computes injected and produced volumes of all phases. + /// Note 1: assumes that only the first phase is injected. + /// Note 2: assumes that transport has been done with an + /// implicit method, i.e. that the current state + /// gives the mobilities used for the preceding timestep. + /// @param[in] props fluid and rock properties. + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all P phases) + /// @param[in] src if < 0: total outflow, if > 0: first phase inflow. + /// @param[in] dt timestep used + /// @param[out] injected must point to a valid array with P elements, + /// where P = s.size()/src.size(). + /// @param[out] produced must also point to a valid array with P elements. + 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); + + /// @brief Computes total mobility for a set of saturation values. + /// @param[in] props rock and fluid properties + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @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& p, + const std::vector& z, + 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 + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all phases) + /// @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& p, + const std::vector& z, + const std::vector& s, + std::vector& totmob, + std::vector& omega); + + + /// @brief Computes phase mobilities for a set of saturation values. + /// @param[in] props rock and fluid properties + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all phases) + /// @param[out] pmobc phase mobilities (for all phases). + void computePhaseMobilities(const Opm::BlackoilPropertiesInterface& props, + const std::vector& cells, + const std::vector& p, + 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 + /// @param[in] cells cells with which the saturation values are associated + /// @param[in] p pressure (one value per cell) + /// @param[in] z surface-volume values (for all P phases) + /// @param[in] s saturation values (for all phases) + /// @param[out] fractional_flow the fractional flow for each phase for each cell. + void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props, + const std::vector& cells, + const std::vector& p, + const std::vector& z, + const std::vector& s, + std::vector& fractional_flows); + +} // namespace Opm + +#endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED From 89b4f8ceb0f60a3fe52a239ac8e9fc164b3f65b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 May 2012 12:50:02 +0200 Subject: [PATCH 09/31] Added WellReport::push() overload taking BlackoilPropertiesInterface. --- opm/core/utility/miscUtilities.cpp | 55 ++++++++++++++++++++++++++++++ opm/core/utility/miscUtilities.hpp | 10 ++++++ 2 files changed, 65 insertions(+) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index b84f94e6..289b68ba 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -614,6 +615,60 @@ namespace Opm + void WellReport::push(const BlackoilPropertiesInterface& props, + const Wells& wells, + const std::vector& p, + const std::vector& z, + const std::vector& s, + const double time, + const std::vector& well_bhp, + const std::vector& well_perfrates) + { + // TODO: refactor, since this is almost identical to the other push(). + int nw = well_bhp.size(); + ASSERT(nw == wells.number_of_wells); + if (props.numPhases() != 2) { + THROW("WellReport for now assumes two phase flow."); + } + std::vector data_now; + data_now.reserve(1 + 3*nw); + data_now.push_back(time/unit::day); + for (int w = 0; w < nw; ++w) { + data_now.push_back(well_bhp[w]/(unit::barsa)); + double well_rate_total = 0.0; + double well_rate_water = 0.0; + for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { + const double perf_rate = well_perfrates[perf]*(unit::day/unit::second); + well_rate_total += perf_rate; + if (perf_rate > 0.0) { + // Injection. + well_rate_water += perf_rate*wells.comp_frac[0]; + } else { + // Production. + const int cell = wells.well_cells[perf]; + double mob[2]; + props.relperm(1, &s[2*cell], &cell, mob, 0); + double visc[2]; + props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0); + mob[0] /= visc[0]; + mob[1] /= visc[1]; + const double fracflow = mob[0]/(mob[0] + mob[1]); + well_rate_water += perf_rate*fracflow; + } + } + data_now.push_back(well_rate_total); + if (well_rate_total == 0.0) { + data_now.push_back(0.0); + } else { + data_now.push_back(well_rate_water/well_rate_total); + } + } + data_.push_back(data_now); + } + + + + void WellReport::write(std::ostream& os) const { const int sz = data_.size(); diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp index 40fe0059..37b252bc 100644 --- a/opm/core/utility/miscUtilities.hpp +++ b/opm/core/utility/miscUtilities.hpp @@ -30,6 +30,7 @@ namespace Opm { class IncompPropertiesInterface; + class BlackoilPropertiesInterface; class RockCompressibility; /// @brief Computes pore volume of all cells in a grid. @@ -247,6 +248,15 @@ namespace Opm const double time, const std::vector& well_bhp, const std::vector& well_perfrates); + /// Add a report point (compressible fluids). + void push(const BlackoilPropertiesInterface& props, + const Wells& wells, + const std::vector& p, + const std::vector& z, + const std::vector& s, + const double time, + const std::vector& well_bhp, + const std::vector& well_perfrates); /// Write report to a stream. void write(std::ostream& os) const; private: From c2e7928758a1d4becfcbaf5498bb7cc14b215bc7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 May 2012 12:52:33 +0200 Subject: [PATCH 10/31] Added BlackoilState class. --- Makefile.am | 1 + opm/core/BlackoilState.hpp | 103 +++++++++++++++++++++++++++++++++++++ 2 files changed, 104 insertions(+) create mode 100644 opm/core/BlackoilState.hpp diff --git a/Makefile.am b/Makefile.am index b481d607..7d365179 100644 --- a/Makefile.am +++ b/Makefile.am @@ -196,6 +196,7 @@ opm/core/WellsGroup.hpp \ opm/core/WellCollection.hpp \ opm/core/InjectionSpecification.hpp \ opm/core/ProductionSpecification.hpp \ +opm/core/BlackoilState.hpp \ opm/core/ColumnExtract.hpp \ opm/core/GridAdapter.hpp \ opm/core/GridManager.hpp \ diff --git a/opm/core/BlackoilState.hpp b/opm/core/BlackoilState.hpp new file mode 100644 index 00000000..6a9e6118 --- /dev/null +++ b/opm/core/BlackoilState.hpp @@ -0,0 +1,103 @@ +/* + 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_BLACKOILSTATE_HEADER_INCLUDED +#define OPM_BLACKOILSTATE_HEADER_INCLUDED + +#include +#include +#include + +namespace Opm +{ + + /// Simulator state for a blackoil simulator. + class BlackoilState + { + public: + + void init(const UnstructuredGrid& g, const int num_phases) + { + num_phases_ = num_phases; + press_.resize(g.number_of_cells, 0.0); + fpress_.resize(g.number_of_faces, 0.0); + flux_.resize(g.number_of_faces, 0.0); + sat_.resize(num_phases * g.number_of_cells, 0.0); + for (int cell = 0; cell < g.number_of_cells; ++cell) { + // Defaulting the second saturation to 1.0. + // This will usually be oil in a water-oil case, + // gas in an oil-gas case. + // For proper initialization, one should not rely on this, + // but use available phase information instead. + sat_[num_phases*cell + 1] = 1.0; + } + } + + enum ExtremalSat { MinSat, MaxSat }; + + /// Set the first saturation to either its min or max value in + /// the indicated cells. The second saturation value s2 is set + /// to (1.0 - s1) for each cell. Any further saturation values + /// are unchanged. + void setFirstSat(const std::vector& cells, + const Opm::BlackoilPropertiesInterface& props, + ExtremalSat es) + { + const int n = cells.size(); + std::vector smin(num_phases_*n); + std::vector smax(num_phases_*n); + props.satRange(n, &cells[0], &smin[0], &smax[0]); + const double* svals = (es == MinSat) ? &smin[0] : &smax[0]; + for (int ci = 0; ci < n; ++ci) { + const int cell = cells[ci]; + sat_[num_phases_*cell] = svals[num_phases_*ci]; + sat_[num_phases_*cell + 1] = 1.0 - sat_[num_phases_*cell]; + } + } + + int numPhases() const + { + return num_phases_; + } + + std::vector& pressure () { return press_ ; } + std::vector& facepressure() { return fpress_; } + std::vector& faceflux () { return flux_ ; } + std::vector& surfacevol () { return surfvol_; } + std::vector& saturation () { return sat_ ; } + + const std::vector& pressure () const { return press_ ; } + const std::vector& facepressure() const { return fpress_; } + const std::vector& faceflux () const { return flux_ ; } + const std::vector& surfacevol () const { return surfvol_; } + const std::vector& saturation () const { return sat_ ; } + + private: + int num_phases_; + std::vector press_ ; + std::vector fpress_; + std::vector flux_ ; + std::vector surfvol_; + std::vector sat_ ; + }; + +} // namespace Opm + + +#endif // OPM_BLACKOILSTATE_HEADER_INCLUDED From 0cba414782b7f26d89039dbf14f2427c7f1951aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 May 2012 12:53:04 +0200 Subject: [PATCH 11/31] First compiling version, large sections are disabled. --- examples/sim_wateroil.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index f37d56a4..9b63475e 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -44,7 +44,7 @@ #include #include -#include +#include #include #include @@ -71,9 +71,9 @@ #define TRANSPORT_SOLVER_FIXED 0 - +template static void outputState(const UnstructuredGrid& grid, - const Opm::TwophaseState& state, + const State& state, const int step, const std::string& output_dir) { @@ -169,7 +169,7 @@ main(int argc, char** argv) boost::scoped_ptr wells; boost::scoped_ptr rock_comp; Opm::SimulatorTimer simtimer; - Opm::TwophaseState state; + Opm::BlackoilState state; bool check_well_controls = false; int max_well_control_iterations = 0; double gravity[3] = { 0.0 }; @@ -366,7 +366,9 @@ main(int argc, char** argv) well_bhp.resize(num_wells, 0.0); well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0); well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); - wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates); + wellreport.push(*props, *wells->c_wells(), + state.pressure(), state.surfacevol(), state.saturation(), + 0.0, well_bhp, well_perfrates); } for (; !simtimer.done(); ++simtimer) { // Report timestep and (optionally) write state to disk. @@ -376,6 +378,7 @@ main(int argc, char** argv) } // Solve pressure. +#if PRESSURE_SOLVER_FIXED if (use_gravity) { computeTotalMobilityOmega(*props, allcells, state.saturation(), totmob, omega); } else { @@ -395,7 +398,6 @@ main(int argc, char** argv) int well_control_iteration = 0; do { // Well control outer loop. pressure_timer.start(); -#if PRESSURE_SOLVER_FIXED if (rock_comp->isActive()) { rc.resize(num_cells); std::vector initial_pressure = state.pressure(); @@ -439,13 +441,11 @@ main(int argc, char** argv) psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), well_bhp, well_perfrates); } -#endif // PRESSURE_SOLVER_FIXED pressure_timer.stop(); double pt = pressure_timer.secsSinceStart(); std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; ptime += pt; - if (check_well_controls) { Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), fractional_flows, @@ -465,6 +465,7 @@ main(int argc, char** argv) } } } while (!well_control_passed); +#endif // PRESSURE_SOLVER_FIXED // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, @@ -532,7 +533,8 @@ main(int argc, char** argv) produced[0]/(produced[0] + produced[1]), tot_produced[0]/tot_porevol_init); if (wells->c_wells()) { - wellreport.push(*props, *wells->c_wells(), state.saturation(), + wellreport.push(*props, *wells->c_wells(), + state.pressure(), state.surfacevol(), state.saturation(), simtimer.currentTime() + simtimer.currentStepLength(), well_bhp, well_perfrates); } From 20cf0eda368929ebb980b74c446ce7409923bcb5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 14:34:56 +0200 Subject: [PATCH 12/31] compute_cell_contrib(): Advance derivative pointer for each connection. Failing to do this operation resulted in incorrect matrices in cases with anisotropic tensors and/or non-cube cells. The error has been present since the inception of this implementation. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 6bf0c7a5..e1e42c66 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -606,6 +606,8 @@ compute_cell_contrib(struct UnstructuredGrid *G , pimpl->ratio->mat_row[ 0 ] += s * dt * dF1; pimpl->ratio->mat_row[ off ] += s * dt * dF2; + + dv += 2 * np; /* '2' == number of one-sided derivatives. */ } } } From e66aa8b06a6154125a74d7cf144f99c64b8b5db1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 17:19:35 +0200 Subject: [PATCH 13/31] toInjectorType(): Distinguish types based on first characters only. The manual states that injection types in the "WCONINJE" keyword need only be specified using a single character. --- opm/core/WellsGroup.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index f892c83e..d221e91e 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -964,13 +964,13 @@ namespace Opm InjectionSpecification::InjectorType toInjectorType(std::string type) { - if (type == "OIL") { + if (type[0] == 'O') { return InjectionSpecification::OIL; } - if (type == "WATER") { + if (type[0] == 'W') { return InjectionSpecification::WATER; } - if (type == "GAS") { + if (type[0] == 'G') { return InjectionSpecification::GAS; } THROW("Unknown type " << type << ", could not convert to SurfaceComponent"); From 1b6e9d40e2e11f4efb289141e54cfaaa6c251cda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 17:22:05 +0200 Subject: [PATCH 14/31] toInjectorType(): Pass parameter as "reference-to-const" rather than copied object. --- opm/core/WellsGroup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index d221e91e..bfedccaf 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -962,7 +962,7 @@ namespace Opm namespace { - InjectionSpecification::InjectorType toInjectorType(std::string type) + InjectionSpecification::InjectorType toInjectorType(const std::string& type) { if (type[0] == 'O') { return InjectionSpecification::OIL; From aaeedd867933fed8e583e52c856e48d67d1ca186 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 16 May 2012 00:50:23 +0200 Subject: [PATCH 15/31] WCONINJE: Distinguish injectors based on first character only. This is the completion of change-set e6015b19c4e8 from WellsGroup.cpp . --- opm/core/WellsManager.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index 6fe66a1e..8253d36a 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -467,17 +467,17 @@ namespace Opm // Set well component fraction. double cf[3] = { 0.0, 0.0, 0.0 }; - if (wci_line.injector_type_ == "WATER") { + if (wci_line.injector_type_[0] == 'W') { if (!pu.phase_used[BlackoilPhases::Aqua]) { THROW("Water phase not used, yet found water-injecting well."); } cf[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; - } else if (wci_line.injector_type_ == "OIL") { + } else if (wci_line.injector_type_[0] == 'O') { if (!pu.phase_used[BlackoilPhases::Liquid]) { THROW("Oil phase not used, yet found oil-injecting well."); } cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; - } else if (wci_line.injector_type_ == "GAS") { + } else if (wci_line.injector_type_[0] == 'G') { if (!pu.phase_used[BlackoilPhases::Vapour]) { THROW("Water phase not used, yet found water-injecting well."); } From 58756f0ba5333af1c77b6fafafc760c57bd91d85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 09:24:34 +0200 Subject: [PATCH 16/31] Remove unneeded debug output. --- opm/core/WellsGroup.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index bfedccaf..903cc9b2 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -1076,11 +1076,6 @@ namespace Opm if (deck.hasField("WCONPROD")) { WCONPROD wconprod = deck.getWCONPROD(); - -#if THIS_STATEMENT_IS_REALLY_NEEDED - std::cout << wconprod.wconprod.size() << std::endl; -#endif - for (size_t i = 0; i < wconprod.wconprod.size(); i++) { if (wconprod.wconprod[i].well_ == name) { WconprodLine line = wconprod.wconprod[i]; From f619d1d909c7d5476b43b5a85c39293a338b7981 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 16 May 2012 11:13:38 +0200 Subject: [PATCH 17/31] Don't #include EclipseGridInspector for the side effect of --- opm/core/fluid/RockFromDeck.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 4e45b0d1..94777904 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -19,7 +19,8 @@ #include -#include + +#include namespace Opm { From 9379263646cb0f8af578c6c02945425528f94984 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 11:37:31 +0200 Subject: [PATCH 18/31] Renamed initStateTwophaseFromDeck() -> initStateFromDeck(). - Made initStateFromDeck() into a template taking arbitrary properties. Implementation detail: - initWaterOilContact() was also templatized on props. - initHydrostaticPressure() is now overloaded on prop interface types. --- examples/sim_wateroil.cpp | 8 +- examples/spu_2p.cpp | 2 +- examples/wells_example.cpp | 2 +- opm/core/TwophaseState.hpp | 25 ++-- opm/core/utility/initState.hpp | 18 +-- opm/core/utility/initState_impl.hpp | 191 ++++++++++++++++++++++++---- 6 files changed, 198 insertions(+), 48 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 9b63475e..94bdfa04 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -66,7 +66,7 @@ #include #include -#define COMPR_INIT_FIXED 0 +#define COMPR_INIT_FIXED 1 #define PRESSURE_SOLVER_FIXED 0 #define TRANSPORT_SOLVER_FIXED 0 @@ -199,9 +199,9 @@ main(int argc, char** argv) // Init state variables (saturation and pressure). #if COMPR_INIT_FIXED if (param.has("init_saturation")) { - initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); + initStateBlackoilBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { - initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } #endif // COMPR_INIT_FIXED } else { @@ -225,7 +225,7 @@ main(int argc, char** argv) gravity[2] = param.getDefault("gravity", 0.0); // Init state variables (saturation and pressure). #if COMPR_INIT_FIXED - initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); + initStateBlackoilBasic(*grid->c_grid(), *props, param, gravity[2], state); #endif // COMPR_INIT_FIXED } diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index e34d0441..366c1bb6 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -329,7 +329,7 @@ main(int argc, char** argv) if (param.has("init_saturation")) { initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { - initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } } else { // Grid init. diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp index 2c070991..18c051aa 100644 --- a/examples/wells_example.cpp +++ b/examples/wells_example.cpp @@ -58,7 +58,7 @@ int main(int argc, char** argv) Opm::TwophaseState state; - initStateTwophaseFromDeck(*grid.c_grid(), incomp_properties, parser, gravity[2], state); + initStateFromDeck(*grid.c_grid(), incomp_properties, parser, gravity[2], state); // Compute phase mobilities std::vector phase_mob; diff --git a/opm/core/TwophaseState.hpp b/opm/core/TwophaseState.hpp index 74e4b0ad..499a5b4f 100644 --- a/opm/core/TwophaseState.hpp +++ b/opm/core/TwophaseState.hpp @@ -32,38 +32,44 @@ namespace Opm { public: - void init(const UnstructuredGrid& g) + void init(const UnstructuredGrid& g, int num_phases) { + num_phases_ = num_phases; press_.resize(g.number_of_cells, 0.0); fpress_.resize(g.number_of_faces, 0.0); flux_.resize(g.number_of_faces, 0.0); - sat_.resize(2 * g.number_of_cells, 0.0); + sat_.resize(num_phases_ * g.number_of_cells, 0.0); for (int cell = 0; cell < g.number_of_cells; ++cell) { - sat_[2*cell + 1] = 1.0; // Defaulting oil saturations to 1.0. + // Defaulting the second saturation to 1.0. + // This will usually be oil in a water-oil case, + // gas in an oil-gas case. + // For proper initialization, one should not rely on this, + // but use available phase information instead. + sat_[num_phases_*cell + 1] = 1.0; } } enum ExtremalSat { MinSat, MaxSat }; - void setWaterSat(const std::vector& cells, + void setFirstSat(const std::vector& cells, const Opm::IncompPropertiesInterface& props, ExtremalSat es) { const int n = cells.size(); - std::vector smin(2*n); - std::vector smax(2*n); + std::vector smin(num_phases_*n); + std::vector smax(num_phases_*n); props.satRange(n, &cells[0], &smin[0], &smax[0]); const double* svals = (es == MinSat) ? &smin[0] : &smax[0]; for (int ci = 0; ci < n; ++ci) { const int cell = cells[ci]; - sat_[2*cell] = svals[2*ci]; - sat_[2*cell + 1] = 1.0 - sat_[2*cell]; + sat_[num_phases_*cell] = svals[num_phases_*ci]; + sat_[num_phases_*cell + 1] = 1.0 - sat_[num_phases_*cell]; } } int numPhases() const { - return 2; + return num_phases_; } std::vector& pressure () { return press_ ; } @@ -77,6 +83,7 @@ namespace Opm const std::vector& saturation () const { return sat_ ; } private: + int num_phases_; std::vector press_ ; std::vector fpress_; std::vector flux_ ; diff --git a/opm/core/utility/initState.hpp b/opm/core/utility/initState.hpp index a3f5a95a..981c9fad 100644 --- a/opm/core/utility/initState.hpp +++ b/opm/core/utility/initState.hpp @@ -28,8 +28,9 @@ namespace Opm namespace parameter { class ParameterGroup; } class EclipseGridParser; class IncompPropertiesInterface; + class BlackoilPropertiesInterface; - /// Initialize a state from parameters. + /// Initialize a two-phase state from parameters. /// The following parameters are accepted (defaults): /// convection_testcase (false) Water in the 'left' part of the grid. /// ref_pressure (100) Initial pressure in bar for all cells @@ -58,19 +59,20 @@ namespace Opm const double gravity, State& state); - /// Initialize a state from input deck. + /// Initialize a two-phase state from input deck. /// If EQUIL is present: /// - saturation is set according to the water-oil contact, /// - pressure is set to hydrostatic equilibrium. /// Otherwise: /// - saturation is set according to SWAT, /// - pressure is set according to PRESSURE. - template - void initStateTwophaseFromDeck(const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const EclipseGridParser& deck, - const double gravity, - State& state); + template + void initStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state); + } // namespace Opm diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 34f944f8..d61e27fe 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -25,8 +25,10 @@ #include #include #include +#include #include #include +#include namespace Opm @@ -62,9 +64,9 @@ namespace Opm // Initialize saturations so that there is water below woc, // and oil above. // If invert is true, water is instead above, oil below. - template + template void initWaterOilContact(const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, + const Props& props, const double woc, const WaterInit waterinit, State& state) @@ -73,10 +75,10 @@ namespace Opm std::vector water; std::vector oil; // Assuming that water should go below the woc, but warning if oil is heavier. - if (props.density()[0] < props.density()[1]) { - std::cout << "*** warning: water density is less than oil density, " - "but initialising water below woc." << std::endl; - } + // if (props.density()[0] < props.density()[1]) { + // std::cout << "*** warning: water density is less than oil density, " + // "but initialising water below woc." << std::endl; + // } switch (waterinit) { case WaterBelow: cellsBelowAbove(grid, woc, water, oil); @@ -85,10 +87,11 @@ namespace Opm cellsBelowAbove(grid, woc, oil, water); } // Set saturations. - state.setWaterSat(oil, props, State::MinSat); - state.setWaterSat(water, props, State::MaxSat); + state.setFirstSat(oil, props, State::MinSat); + state.setFirstSat(water, props, State::MaxSat); } + // Initialize hydrostatic pressures depending only on gravity, // (constant) phase densities and a water-oil contact depth. // The pressure difference between points is equal to @@ -123,11 +126,147 @@ namespace Opm } } + + // Facade to initHydrostaticPressure() taking a property object, + // for similarity to initHydrostaticPressure() for compressible fluids. + template + void initHydrostaticPressure(const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const double woc, + const double gravity, + const double datum_z, + const double datum_p, + State& state) + { + const double* densities = props.density(); + initHydrostaticPressure(grid, densities, woc, gravity, datum_z, datum_p, state); + } + + + + // Helper functor for initHydrostaticPressure() for compressible fluids. + struct Density + { + const BlackoilPropertiesInterface& props_; + Density(const BlackoilPropertiesInterface& props) : props_(props) {} + double operator()(const double pressure, const int phase) + { + ASSERT(props_.numPhases() == 2); + const double surfvol[2][2] = { { 1.0, 0.0 }, + { 0.0, 1.0 } }; + // We do not handle multi-region PVT/EQUIL at this point. + const int* cells = 0; + double A[4] = { 0.0 }; + props_.matrix(1, &pressure, surfvol[phase], cells, A, 0); + double rho[2] = { 0.0 }; + props_.density(1, A, rho); + return rho[phase]; + } + }; + + // Initialize hydrostatic pressures depending only on gravity, + // phase densities that may vary with pressure and a water-oil + // contact depth. The pressure ODE is given as + // \grad p = \rho g \grad z + // where rho is the oil density above the woc, water density + // below woc. Note that even if there is (immobile) water in + // the oil zone, it does not contribute to the pressure there. + template + void initHydrostaticPressure(const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const double woc, + const double gravity, + const double datum_z, + const double datum_p, + State& state) + { + ASSERT(props.numPhases() == 2); + + // Obtain max and min z for which we will need to compute p. + const int num_cells = grid.number_of_cells; + const int dim = grid.dimensions; + double zlim[2] = { 1e100, -1e100 }; + for (int c = 0; c < num_cells; ++c) { + const double z = grid.cell_centroids[dim*c + dim - 1]; + zlim[0] = std::min(zlim[0], z); + zlim[1] = std::max(zlim[1], z); + } + + // We'll use a minimum stepsize of 1/100 of the z range. + const double hmin = (zlim[1] - zlim[0])/100.0; + + // Store p(z) results in an associative array. + std::map press_by_z; + press_by_z[datum_z] = datum_p; + + // Set up density evaluator. + Density rho(props); + + // Solve the ODE from datum_z to woc. + int phase = (datum_z > woc) ? 0 : 1; + int num_steps = int(std::ceil(std::fabs(woc - datum_z)/hmin)); + double pval = datum_p; + double zval = datum_z; + double h = (woc - datum_z)/double(num_steps); + for (int i = 0; i < num_steps; ++i) { + zval += h; + const double dp = rho(pval, phase)*gravity; + pval += h*dp; + press_by_z[zval] = pval; + } + double woc_p = pval; + + // Solve the ODE from datum_z to the end of the interval. + double z_end = (datum_z > woc) ? zlim[1] : zlim[0]; + num_steps = int(std::ceil(std::fabs(z_end - datum_z)/hmin)); + pval = datum_p; + zval = datum_z; + h = (z_end - datum_z)/double(num_steps); + for (int i = 0; i < num_steps; ++i) { + zval += h; + const double dp = rho(pval, phase)*gravity; + pval += h*dp; + press_by_z[zval] = pval; + } + + // Solve the ODE from woc to the other end of the interval. + // Switching phase and z_end. + phase = (phase + 1) % 2; + z_end = (datum_z > woc) ? zlim[0] : zlim[1]; + pval = woc_p; + zval = woc; + h = (z_end - datum_z)/double(num_steps); + for (int i = 0; i < num_steps; ++i) { + zval += h; + const double dp = rho(pval, phase)*gravity; + pval += h*dp; + press_by_z[zval] = pval; + } + + // Create monotone spline for interpolating solution. + std::vector zv, pv; + zv.reserve(press_by_z.size()); + pv.reserve(press_by_z.size()); + for (std::map::const_iterator it = press_by_z.begin(); + it != press_by_z.end(); ++it) { + zv.push_back(it->first); + pv.push_back(it->second); + } + MonotCubicInterpolator press(zv, pv); + + // Evaluate pressure at each cell centroid. + std::vector& p = state.pressure(); + for (int c = 0; c < num_cells; ++c) { + const double z = grid.cell_centroids[dim*c + dim - 1]; + p[c] = press(z); + } + } + } // anonymous namespace - /// Initialize a state from parameters. + /// Initialize a twophase state from parameters. /// The following parameters are accepted (defaults): /// convection_testcase (false) Water in the 'left' part of the grid. /// ref_pressure (100) Initial pressure in bar for all cells @@ -156,17 +295,18 @@ namespace Opm const double gravity, State& state) { - if (state.numPhases() != 2) { - THROW("initStateTwophaseFromDeck(): state must have two phases."); + const int num_phases = props.numPhases(); + if (num_phases != 2) { + THROW("initStateTwophaseBasic(): currently handling only two-phase scenarios."); } - state.init(grid); + state.init(grid, num_phases); const int num_cells = props.numCells(); // By default: initialise water saturation to minimum everywhere. std::vector all_cells(num_cells); for (int i = 0; i < num_cells; ++i) { all_cells[i] = i; } - state.setWaterSat(all_cells, props, State::MinSat); + state.setFirstSat(all_cells, props, State::MinSat); const bool convection_testcase = param.getDefault("convection_testcase", false); const bool segregation_testcase = param.getDefault("segregation_testcase", false); if (convection_testcase) { @@ -182,7 +322,7 @@ namespace Opm left_cells.push_back(cell); } } - state.setWaterSat(left_cells, props, State::MaxSat); + state.setFirstSat(left_cells, props, State::MaxSat); const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa; std::fill(state.pressure().begin(), state.pressure().end(), init_p); } else if (segregation_testcase) { @@ -246,17 +386,18 @@ namespace Opm /// Otherwise: /// - saturation is set according to SWAT, /// - pressure is set according to PRESSURE. - template - void initStateTwophaseFromDeck(const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const EclipseGridParser& deck, - const double gravity, - State& state) + template + void initStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state) { - if (state.numPhases() != 2) { - THROW("initStateTwophaseFromDeck(): state must have two phases."); + const int num_phases = props.numPhases(); + if (num_phases != 2) { + THROW("initStateFromDeck(): currently handling only two-phase scenarios."); } - state.init(grid); + state.init(grid, num_phases); if (deck.hasField("EQUIL")) { // Set saturations depending on oil-water contact. const EQUIL& equil= deck.getEQUIL(); @@ -268,7 +409,7 @@ namespace Opm // Set pressure depending on densities and depths. const double datum_z = equil.equil[0].datum_depth_; const double datum_p = equil.equil[0].datum_depth_pressure_; - initHydrostaticPressure(grid, props.density(), woc, gravity, datum_z, datum_p, state); + initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); } else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) { // Set saturations from SWAT, pressure from PRESSURE. std::vector& s = state.saturation(); @@ -283,7 +424,7 @@ namespace Opm p[c] = p_deck[c_deck]; } } else { - THROW("initStateTwophaseFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); + THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); } } From 4c37676338fc79001fd212539e9e414d64f530bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 12:33:42 +0200 Subject: [PATCH 19/31] Renamed initTwophaseStateBasic() -> initStateBasic(). --- examples/spu_2p.cpp | 4 ++-- opm/core/utility/initState.hpp | 10 +++++----- opm/core/utility/initState_impl.hpp | 10 +++++----- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 366c1bb6..1f9fce7c 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -327,7 +327,7 @@ main(int argc, char** argv) gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; // Init state variables (saturation and pressure). if (param.has("init_saturation")) { - initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } @@ -351,7 +351,7 @@ main(int argc, char** argv) // Gravity. gravity[2] = param.getDefault("gravity", 0.0); // Init state variables (saturation and pressure). - initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); } // Extra fluid init for transport solver. diff --git a/opm/core/utility/initState.hpp b/opm/core/utility/initState.hpp index 981c9fad..b901ed69 100644 --- a/opm/core/utility/initState.hpp +++ b/opm/core/utility/initState.hpp @@ -53,11 +53,11 @@ namespace Opm /// In all three cases, pressure is initialised hydrostatically. /// In case 2) and 3), the depth of the first cell is used as reference depth. template - void initStateTwophaseBasic(const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const parameter::ParameterGroup& param, - const double gravity, - State& state); + void initStateBasic(const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state); /// Initialize a two-phase state from input deck. /// If EQUIL is present: diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index d61e27fe..e4dddb19 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -289,11 +289,11 @@ namespace Opm /// In all three cases, pressure is initialised hydrostatically. /// In case 2) and 3), the depth of the first cell is used as reference depth. template - void initStateTwophaseBasic(const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const parameter::ParameterGroup& param, - const double gravity, - State& state) + void initStateBasic(const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state) { const int num_phases = props.numPhases(); if (num_phases != 2) { From 99330b219e424dfc8fdaabb13b9156775872e34e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 12:52:58 +0200 Subject: [PATCH 20/31] Silence a warning. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 070cba05..44e14559 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -46,7 +46,7 @@ namespace Opm /// where P = s.size()/src.size(). /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const BlackoilPropertiesInterface& props, - const std::vector& p, + const std::vector& press, const std::vector& z, const std::vector& s, const std::vector& src, @@ -70,7 +70,7 @@ namespace Opm const double flux = -src[c]*dt; const double* sat = &s[np*c]; props.relperm(1, sat, &c, &mob[0], 0); - props.viscosity(1, &p[c], &z[np*c], &c, &visc[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]; @@ -94,14 +94,14 @@ namespace Opm /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, const std::vector& cells, - const std::vector& p, + const std::vector& press, const std::vector& z, const std::vector& s, std::vector& totmob) { std::vector pmobc; - computePhaseMobilities(props, cells, p, z, s, pmobc); + computePhaseMobilities(props, cells, press, z, s, pmobc); const std::size_t np = props.numPhases(); const std::vector::size_type nc = cells.size(); From 3d44427065cbb1497e57edd04bc021deabe2e887 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 12:54:12 +0200 Subject: [PATCH 21/31] Follow changes to TwophaseState class and initState*() functions. --- tutorials/tutorial3.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 0e9f4362..4ba293bc 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -372,11 +372,12 @@ int main () /// \page tutorial3 /// \details - /// We initialise water saturation to minimum everywhere. + /// We set up a two-phase state object, and + /// initialise water saturation to minimum everywhere. /// \code TwophaseState state; - state.init(*grid.c_grid()); - state.setWaterSat(allcells, props, TwophaseState::MinSat); + state.init(*grid.c_grid(), 2); + state.setFirstSat(allcells, props, TwophaseState::MinSat); /// \endcode /// \page tutorial3 From 79cdeaa5dff1f2f089a34c9a8c870854695fcc83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 12:54:48 +0200 Subject: [PATCH 22/31] Added initStateBasic() overload taking BlackoilPropertiesInterface props. --- examples/sim_wateroil.cpp | 4 +- opm/core/utility/initState.hpp | 23 +++++++++ opm/core/utility/initState_impl.hpp | 79 ++++++++++++++++++++++++++++- 3 files changed, 103 insertions(+), 3 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 94bdfa04..6833bead 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -199,7 +199,7 @@ main(int argc, char** argv) // Init state variables (saturation and pressure). #if COMPR_INIT_FIXED if (param.has("init_saturation")) { - initStateBlackoilBasic(*grid->c_grid(), *props, param, gravity[2], state); + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } @@ -225,7 +225,7 @@ main(int argc, char** argv) gravity[2] = param.getDefault("gravity", 0.0); // Init state variables (saturation and pressure). #if COMPR_INIT_FIXED - initStateBlackoilBasic(*grid->c_grid(), *props, param, gravity[2], state); + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); #endif // COMPR_INIT_FIXED } diff --git a/opm/core/utility/initState.hpp b/opm/core/utility/initState.hpp index b901ed69..e49a692c 100644 --- a/opm/core/utility/initState.hpp +++ b/opm/core/utility/initState.hpp @@ -59,6 +59,29 @@ namespace Opm const double gravity, State& state); + /// Initialize a blackoil state from parameters. + /// The following parameters are accepted (defaults): + /// convection_testcase (false) Water in the 'left' part of the grid. + /// ref_pressure (100) Initial pressure in bar for all cells + /// (if convection_testcase is true), + /// or pressure at woc depth. + /// water_oil_contact (none) Depth of water-oil contact (woc). + /// If convection_testcase is true, the saturation is initialised + /// as indicated, and pressure is initialised to a constant value + /// ('ref_pressure'). + /// Otherwise we have 2 cases: + /// 1) If 'water_oil_contact' is given, saturation is initialised + /// accordingly. + /// 2) Water saturation is set to minimum. + /// In both cases, pressure is initialised hydrostatically. + /// In case 2), the depth of the first cell is used as reference depth. + template + void initStateBasic(const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state); + /// Initialize a two-phase state from input deck. /// If EQUIL is present: /// - saturation is set according to the water-oil contact, diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index e4dddb19..7872ca98 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -265,7 +265,6 @@ namespace Opm } // anonymous namespace - /// Initialize a twophase state from parameters. /// The following parameters are accepted (defaults): /// convection_testcase (false) Water in the 'left' part of the grid. @@ -379,6 +378,84 @@ namespace Opm } + /// Initialize a blackoil state from parameters. + /// The following parameters are accepted (defaults): + /// convection_testcase (false) Water in the 'left' part of the grid. + /// ref_pressure (100) Initial pressure in bar for all cells + /// (if convection_testcase is true), + /// or pressure at woc depth. + /// water_oil_contact (none) Depth of water-oil contact (woc). + /// If convection_testcase is true, the saturation is initialised + /// as indicated, and pressure is initialised to a constant value + /// ('ref_pressure'). + /// Otherwise we have 2 cases: + /// 1) If 'water_oil_contact' is given, saturation is initialised + /// accordingly. + /// 2) Water saturation is set to minimum. + /// In both cases, pressure is initialised hydrostatically. + /// In case 2), the depth of the first cell is used as reference depth. + template + void initStateBasic(const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state) + { + // TODO: Refactor to exploit similarity with IncompProp* case. + const int num_phases = props.numPhases(); + if (num_phases != 2) { + THROW("initStateTwophaseBasic(): currently handling only two-phase scenarios."); + } + state.init(grid, num_phases); + const int num_cells = props.numCells(); + // By default: initialise water saturation to minimum everywhere. + std::vector all_cells(num_cells); + for (int i = 0; i < num_cells; ++i) { + all_cells[i] = i; + } + state.setFirstSat(all_cells, props, State::MinSat); + const bool convection_testcase = param.getDefault("convection_testcase", false); + if (convection_testcase) { + // Initialise water saturation to max in the 'left' part. + std::vector left_cells; + left_cells.reserve(num_cells/2); + const int *glob_cell = grid.global_cell; + const int* cd = grid.cartdims; + for (int cell = 0; cell < num_cells; ++cell) { + const int gc = glob_cell == 0 ? cell : glob_cell[cell]; + bool left = (gc % cd[0]) < cd[0]/2; + if (left) { + left_cells.push_back(cell); + } + } + state.setFirstSat(left_cells, props, State::MaxSat); + const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa; + std::fill(state.pressure().begin(), state.pressure().end(), init_p); + } else if (param.has("water_oil_contact")) { + // Warn against error-prone usage. + if (gravity == 0.0) { + std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl; + } + if (grid.cartdims[2] <= 1) { + std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl; + } + // Initialise water saturation to max below water-oil contact. + const double woc = param.get("water_oil_contact"); + initWaterOilContact(grid, props, woc, WaterBelow, state); + // Initialise pressure to hydrostatic state. + const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa; + initHydrostaticPressure(grid, props, woc, gravity, woc, ref_p, state); + } else { + // Use default: water saturation is minimum everywhere. + // Initialise pressure to hydrostatic state. + const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; + const double woc = -1e100; + initHydrostaticPressure(grid, props, woc, gravity, ref_z, ref_p, state); + } + } + + /// Initialize a state from input deck. /// If EQUIL is present: /// - saturation is set according to the water-oil contact, From cc0097d8b98ac222492ce1ef2a7dc1c5d0b059b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 14:15:50 +0200 Subject: [PATCH 23/31] Fix bug related to arithmetic if operator ( ? : ). --- opm/core/pressure/CompressibleTpfa.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index cdc4db54..564eeb8f 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -59,7 +59,7 @@ namespace Opm THROW("Inconsistent number of phases specified: " << wells_->number_of_phases << " != " << num_phases); } - const int num_dofs = g.number_of_cells + wells ? wells->number_of_wells : 0; + const int num_dofs = g.number_of_cells + (wells ? wells->number_of_wells : 0); pressure_increment_.resize(num_dofs); UnstructuredGrid* gg = const_cast(&grid_); tpfa_htrans_compute(gg, permeability, &htrans_[0]); From 3740242d7a511c31d5ba80ab3bc52ae28d166f75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 14:16:12 +0200 Subject: [PATCH 24/31] Removed macro disabling initialization code. --- examples/sim_wateroil.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 6833bead..d067d7b5 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -66,7 +66,6 @@ #include #include -#define COMPR_INIT_FIXED 1 #define PRESSURE_SOLVER_FIXED 0 #define TRANSPORT_SOLVER_FIXED 0 @@ -197,13 +196,11 @@ main(int argc, char** argv) // Gravity. gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; // Init state variables (saturation and pressure). -#if COMPR_INIT_FIXED if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } -#endif // COMPR_INIT_FIXED } else { // Grid init. const int nx = param.getDefault("nx", 100); @@ -224,9 +221,7 @@ main(int argc, char** argv) // Gravity. gravity[2] = param.getDefault("gravity", 0.0); // Init state variables (saturation and pressure). -#if COMPR_INIT_FIXED initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); -#endif // COMPR_INIT_FIXED } // Warn if gravity but no density difference. From bd5b8e9d29d474711e6eff1ca66e4fdc9e1c9df2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 14:37:39 +0200 Subject: [PATCH 25/31] Added necessary #include. --- opm/core/BlackoilState.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/BlackoilState.hpp b/opm/core/BlackoilState.hpp index 6a9e6118..6da857a7 100644 --- a/opm/core/BlackoilState.hpp +++ b/opm/core/BlackoilState.hpp @@ -21,7 +21,7 @@ #define OPM_BLACKOILSTATE_HEADER_INCLUDED #include -#include +#include #include namespace Opm From 7bb007f52643247c10e17986d4364e3d431b5257 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 14:37:55 +0200 Subject: [PATCH 26/31] Added WellState class. --- Makefile.am | 1 + opm/core/WellState.hpp | 53 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100644 opm/core/WellState.hpp diff --git a/Makefile.am b/Makefile.am index 7d365179..9823658a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -202,6 +202,7 @@ opm/core/GridAdapter.hpp \ opm/core/GridManager.hpp \ opm/core/TwophaseState.hpp \ opm/core/WellsManager.hpp \ +opm/core/WellState.hpp \ opm/core/pressure/fsh.h \ opm/core/pressure/HybridPressureSolver.hpp \ opm/core/pressure/IncompTpfa.hpp \ diff --git a/opm/core/WellState.hpp b/opm/core/WellState.hpp new file mode 100644 index 00000000..442b7533 --- /dev/null +++ b/opm/core/WellState.hpp @@ -0,0 +1,53 @@ +/* + 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_WELLSTATE_HEADER_INCLUDED +#define OPM_WELLSTATE_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + + class WellState + { + public: + void init(const Wells* wells) + { + if (wells) { + bhp_.resize(wells->number_of_wells); + perfrates_.resize(wells->well_connpos[wells->number_of_wells]); + } + } + + std::vector& bhp() { return bhp_; } + const std::vector& bhp() const { return bhp_; } + + std::vector& perfRates() { return perfrates_; } + const std::vector& perfRates() const { return perfrates_; } + + private: + std::vector bhp_; + std::vector perfrates_; + }; + +} // namespace Opm + +#endif // OPM_WELLSTATE_HEADER_INCLUDED From 9ad4272bc29a8486ac07efb181c0eaca914e71d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 14:38:17 +0200 Subject: [PATCH 27/31] Start using WellState class. --- examples/sim_wateroil.cpp | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index d067d7b5..0bcc22b0 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -45,6 +45,7 @@ #include #include +#include #include #include @@ -351,19 +352,17 @@ main(int argc, char** argv) Opm::Watercut watercut; watercut.push(0.0, 0.0, 0.0); Opm::WellReport wellreport; - std::vector well_bhp; - std::vector well_perfrates; + Opm::WellState well_state; std::vector fractional_flows; std::vector well_resflows_phase; int num_wells = 0; if (wells->c_wells()) { num_wells = wells->c_wells()->number_of_wells; - well_bhp.resize(num_wells, 0.0); - well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0); + well_state.init(wells->c_wells()); well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); wellreport.push(*props, *wells->c_wells(), state.pressure(), state.surfacevol(), state.saturation(), - 0.0, well_bhp, well_perfrates); + 0.0, well_state.bhp(), well_state.perfRates()); } for (; !simtimer.done(); ++simtimer) { // Report timestep and (optionally) write state to disk. @@ -407,7 +406,7 @@ main(int argc, char** argv) } computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); - std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells); + std::copy(well_state.bhp().begin(), well_state.bhp().end(), prev_pressure.begin() + num_cells); // prev_pressure = state.pressure(); // compute pressure increment @@ -421,7 +420,7 @@ main(int argc, char** argv) max_change = std::max(max_change, std::fabs(pressure_increment[cell])); } for (int well = 0; well < num_wells; ++well) { - well_bhp[well] += pressure_increment[num_cells + well]; + well_state.bhp()[well] += pressure_increment[num_cells + well]; max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well])); } @@ -431,10 +430,10 @@ main(int argc, char** argv) } } psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), - well_bhp, well_perfrates); + well_state.bhp(), well_state.perfRates()); } else { psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), - well_bhp, well_perfrates); + well_state.bhp(), well_state.perfRates()); } pressure_timer.stop(); double pt = pressure_timer.secsSinceStart(); @@ -444,11 +443,11 @@ main(int argc, char** argv) if (check_well_controls) { Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), fractional_flows, - well_perfrates, + well_state.perfRates(), well_resflows_phase); std::cout << "Checking well conditions." << std::endl; // For testing we set surface := reservoir - well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); + well_control_passed = wells->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."); @@ -464,7 +463,7 @@ main(int argc, char** argv) // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, - wells->c_wells(), well_perfrates, reorder_src); + wells->c_wells(), well_state.perfRates(), reorder_src); // Solve transport. transport_timer.start(); @@ -531,7 +530,7 @@ main(int argc, char** argv) wellreport.push(*props, *wells->c_wells(), state.pressure(), state.surfacevol(), state.saturation(), simtimer.currentTime() + simtimer.currentStepLength(), - well_bhp, well_perfrates); + well_state.bhp(), well_state.perfRates()); } } total_timer.stop(); From 30218d84ea018b72986d35d282ec577a553bf6f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 14:38:55 +0200 Subject: [PATCH 28/31] Started adding necessary parameters to pressure solver. Work in progress. --- opm/core/pressure/CompressibleTpfa.cpp | 16 +++++++++------- opm/core/pressure/CompressibleTpfa.hpp | 7 +++++-- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 564eeb8f..3fb78204 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include @@ -83,13 +84,14 @@ namespace Opm /// Solve pressure equation, by Newton iterations. - void CompressibleTpfa::solve() + void CompressibleTpfa::solve(const double dt, + BlackoilState& state) { // Set up dynamic data. computeDynamicData(); // Assemble J and F. - assemble(); + assemble(dt, state); bool residual_ok = false; // Replace with tolerance check. while (!residual_ok) { @@ -104,7 +106,7 @@ namespace Opm computeDynamicData(); // Assemble J and F. - assemble(); + assemble(dt, state); // Check for convergence. // Include both tolerance check for residual @@ -127,12 +129,12 @@ namespace Opm /// Compute the residual and Jacobian. - void CompressibleTpfa::assemble() + void CompressibleTpfa::assemble(const double dt, + const BlackoilState& state) { // Arguments or members? - const double dt = 0.0; - const double* z = NULL; - const double* cell_press = NULL; + const double* z = &state.surfacevol()[0]; + const double* cell_press = &state.pressure()[0]; const double* well_bhp = NULL; const double* porevol = NULL; diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index fdacf426..7141100d 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -31,6 +31,7 @@ struct FlowBoundaryConditions; namespace Opm { + class BlackoilState; class LinearSolverInterface; /// Encapsulating a tpfa pressure solver for the compressible-fluid case. @@ -60,11 +61,13 @@ namespace Opm /// Destructor. ~CompressibleTpfa(); - void solve(); + void solve(const double dt, + BlackoilState& state); private: void computeDynamicData(); - void assemble(); + void assemble(const double dt, + const BlackoilState& state); void solveIncrement(); void computeResults(std::vector& pressure, From e57ef29fae3de67dd27712284173fa9a752d4499 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 15:49:02 +0200 Subject: [PATCH 29/31] Work in progress on compressible pressure solver. --- examples/sim_wateroil.cpp | 3 +- opm/core/pressure/CompressibleTpfa.cpp | 57 +++++++++++++++----------- opm/core/pressure/CompressibleTpfa.hpp | 24 ++++++++--- 3 files changed, 55 insertions(+), 29 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 0bcc22b0..3cffd522 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -266,6 +266,7 @@ 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); @@ -296,7 +297,7 @@ main(int argc, char** argv) Opm::LinearSolverFactory linsolver(param); // Pressure solver. const double *grav = use_gravity ? &gravity[0] : 0; - Opm::CompressibleTpfa psolver(*grid->c_grid(), props->permeability(), grav, + Opm::CompressibleTpfa psolver(*grid->c_grid(), props->permeability(), &porevol[0], grav, linsolver, wells->c_wells(), props->numPhases()); // Reordering solver. const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 3fb78204..dc51f586 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -28,6 +28,7 @@ #include #include #include +#include #include @@ -46,11 +47,13 @@ namespace Opm /// is ignored if NULL CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& g, const double* permeability, + const double* porevol, const double* /* gravity */, // ??? const LinearSolverInterface& linsolver, const struct Wells* wells, const int num_phases) : grid_(g), + porevol_(porevol), linsolver_(linsolver), htrans_(g.cell_facepos[ g.number_of_cells ]), trans_ (g.number_of_faces), @@ -85,13 +88,15 @@ namespace Opm /// Solve pressure equation, by Newton iterations. void CompressibleTpfa::solve(const double dt, - BlackoilState& state) + BlackoilState& state, + WellState& well_state) { // Set up dynamic data. - computeDynamicData(); + computePerSolveDynamicData(); + computePerIterationDynamicData(); // Assemble J and F. - assemble(dt, state); + assemble(dt, state, well_state); bool residual_ok = false; // Replace with tolerance check. while (!residual_ok) { @@ -103,10 +108,10 @@ namespace Opm // Update pressure vars with increment. // Set up dynamic data. - computeDynamicData(); + computePerIterationDynamicData(); // Assemble J and F. - assemble(dt, state); + assemble(dt, state, well_state); // Check for convergence. // Include both tolerance check for residual @@ -120,8 +125,16 @@ namespace Opm - /// Solve pressure equation, by Newton iterations. - void CompressibleTpfa::computeDynamicData() + /// Compute per-solve dynamic properties. + void CompressibleTpfa::computePerSolveDynamicData() + { + } + + + + + /// Compute per-iteration dynamic properties. + void CompressibleTpfa::computePerIterationDynamicData() { } @@ -130,19 +143,18 @@ namespace Opm /// Compute the residual and Jacobian. void CompressibleTpfa::assemble(const double dt, - const BlackoilState& state) + const BlackoilState& state, + const WellState& well_state) { - // Arguments or members? const double* z = &state.surfacevol()[0]; const double* cell_press = &state.pressure()[0]; - const double* well_bhp = NULL; - const double* porevol = NULL; - + const double* well_bhp = well_state.bhp().empty() ? NULL : &well_state.bhp()[0]; UnstructuredGrid* gg = const_cast(&grid_); + CompletionData completion_data; - completion_data.gpot = 0; // TODO - completion_data.A = 0; // TODO - completion_data.phasemob = 0; // TODO + completion_data.gpot = &wellperf_gpot_[0]; + completion_data.A = &wellperf_A_[0]; + completion_data.phasemob = &wellperf_phasemob_[0]; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(wells_); wells_tmp.data = &completion_data; @@ -150,15 +162,14 @@ namespace Opm forces.wells = &wells_tmp; forces.src = NULL; // Check if it is legal to leave it as NULL. compr_quantities_gen cq; - cq.Ac = 0; // TODO - cq.dAc = 0; // TODO - cq.Af = 0; // TODO - cq.phasemobf = 0; // TODO - cq.voldiscr = 0; // TODO - // TODO: gravcapf_ must be set already. + cq.Ac = &cell_A_[0]; + cq.dAc = &cell_dA_[0]; + 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], - &gravcapf_[0], &cell_press[0], &well_bhp[0], - &porevol[0], h_); + &face_gravcap_[0], cell_press, well_bhp, + porevol_, h_); } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 7141100d..888b04ec 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -33,6 +33,7 @@ namespace Opm class BlackoilState; class LinearSolverInterface; + class WellState; /// Encapsulating a tpfa pressure solver for the compressible-fluid case. /// Supports gravity, wells and simple sources as driving forces. @@ -52,6 +53,7 @@ namespace Opm /// is ignored if NULL /// \param[in] num_phases Must be 2 or 3. CompressibleTpfa(const UnstructuredGrid& g, + const double* porevol, const double* permeability, const double* gravity, const LinearSolverInterface& linsolver, @@ -62,12 +64,15 @@ namespace Opm ~CompressibleTpfa(); void solve(const double dt, - BlackoilState& state); + BlackoilState& state, + WellState& well_state); private: - void computeDynamicData(); + void computePerSolveDynamicData(); + void computePerIterationDynamicData(); void assemble(const double dt, - const BlackoilState& state); + const BlackoilState& state, + const WellState& well_state); void solveIncrement(); void computeResults(std::vector& pressure, @@ -77,6 +82,7 @@ namespace Opm // ------ Data that will remain unmodified after construction. ------ const UnstructuredGrid& grid_; + const double* porevol_; const LinearSolverInterface& linsolver_; std::vector htrans_; std::vector trans_ ; @@ -86,15 +92,23 @@ namespace Opm struct cfs_tpfa_res_data* h_; // ------ Data that will be modified for every solve. ------ + std::vector wellperf_gpot_; // ------ Data that will be modified for every solver iteration. ------ + // Gravity and capillary contributions (per face). + std::vector face_gravcap_; + std::vector wellperf_A_; + std::vector wellperf_phasemob_; + std::vector cell_A_; + std::vector cell_dA_; + std::vector face_A_; + std::vector face_phasemob_; + std::vector cell_voldisc_; // The update to be applied to the pressures (cell and bhp). std::vector pressure_increment_; - // Gravity and capillary contributions (per face). - std::vector gravcapf_; }; From a1665d9902dfd0d63f9ae12bbf1414a65672b630 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 May 2012 15:56:04 +0200 Subject: [PATCH 30/31] Suppressed some warnings. --- examples/sim_wateroil.cpp | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 3cffd522..db71aaa5 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -143,7 +143,6 @@ main(int argc, char** argv) std::cout << "--------------- Reading parameters ---------------" << std::endl; // Reading various control parameters. - const bool guess_old_solution = param.getDefault("guess_old_solution", false); const bool use_reorder = param.getDefault("use_reorder", true); const bool output = param.getDefault("output", true); std::string output_dir; @@ -160,7 +159,7 @@ main(int argc, char** argv) } output_interval = param.getDefault("output_interval", output_interval); } - const int num_transport_substeps = param.getDefault("num_transport_substeps", 1); + // const int num_transport_substeps = param.getDefault("num_transport_substeps", 1); // If we have a "deck_filename", grid and props will be read from that. bool use_deck = param.has("deck_filename"); @@ -170,8 +169,8 @@ main(int argc, char** argv) boost::scoped_ptr rock_comp; Opm::SimulatorTimer simtimer; Opm::BlackoilState state; - bool check_well_controls = false; - int max_well_control_iterations = 0; + // 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"); @@ -184,8 +183,8 @@ main(int argc, char** argv) props.reset(new Opm::BlackoilPropertiesFromDeck(deck, global_cell)); // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); - check_well_controls = param.getDefault("check_well_controls", false); - max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Timer init. if (deck.hasField("TSTEP")) { simtimer.init(deck); @@ -246,15 +245,15 @@ main(int argc, char** argv) } // 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; + // 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 + // 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. @@ -300,9 +299,9 @@ main(int argc, char** argv) Opm::CompressibleTpfa psolver(*grid->c_grid(), props->permeability(), &porevol[0], grav, linsolver, wells->c_wells(), props->numPhases()); // Reordering solver. +#if TRANSPORT_SOLVER_FIXED const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const int nl_maxiter = param.getDefault("nl_maxiter", 30); -#if TRANSPORT_SOLVER_FIXED Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); if (use_gauss_seidel_gravity) { reorder_model.initGravity(grav); @@ -360,7 +359,7 @@ main(int argc, char** argv) if (wells->c_wells()) { num_wells = wells->c_wells()->number_of_wells; well_state.init(wells->c_wells()); - well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); + well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(num_wells), 0.0); wellreport.push(*props, *wells->c_wells(), state.pressure(), state.surfacevol(), state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); From b7139e07b5d69a4fa2f4dc6b1c29e66c09f57aee Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 16 May 2012 16:08:48 +0200 Subject: [PATCH 31/31] Added necessary #include. --- opm/core/utility/initState_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 7872ca98..6ccab85f 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -29,7 +29,7 @@ #include #include #include - +#include namespace Opm {