From aa49d56d9cf87055ab1fdbe005238f40517ea1ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 10:36:57 +0200 Subject: [PATCH 01/31] Added function cfs_tpfa_res_comprock_assemble(). This is intended to handle cases with both fluid and rock compressibility. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 73 ++++++++++++++++++++++ opm/core/pressure/tpfa/cfs_tpfa_residual.h | 16 +++++ 2 files changed, 89 insertions(+) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 7a52db1e..ab79dbce 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1213,6 +1213,79 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , } +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_res_comprock_assemble( + struct UnstructuredGrid *G , + double dt , + struct cfs_tpfa_res_forces *forces , + const double *zc , + struct compr_quantities_gen *cq , + const double *trans , + const double *gravcap_f, + const double *cpress , + const double *wpress , + const double *porevol , + const double *porevol0 , + const double *rock_comp, + struct cfs_tpfa_res_data *h ) +/* ---------------------------------------------------------------------- */ +{ + /* We want to add this term to the usual residual: + * + * (porevol(pressure)-porevol(initial_pressure))/dt. + * + * Its derivative (for the diagonal term of the Jacobian) is: + * + * porevol(pressure)*rock_comp(pressure)/dt + */ + + int c, w, well_is_neumann, rock_is_incomp; + size_t j; + double dpvdt; + const struct Wells* W; + + /* Assemble usual system (without rock compressibility). */ + cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, + cpress, wpress, porevol, h); + + /* Check if we have only Neumann wells. */ + well_is_neumann = 1; + W = forces->wells->W; + for (w = 0; well_is_neumann && w < W->number_of_wells; w++) { + if ((W->ctrls[w]->current >= 0) && /* OPEN? */ + (W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) { + well_is_neumann = 0; + } + } + + /* If we made a singularity-removing adjustment in the + regular assembly, we undo it here. */ + if (well_is_neumann && h->pimpl->is_incomp) { + h->J->sa[0] /= 2; + } + + /* Add new terms to residual and Jacobian. */ + rock_is_incomp = 1; + for (c = 0; c < G->number_of_cells; c++) { + j = csrmatrix_elm_index(c, c, h->J); + + dpvdt = (porevol[c] - porevol0[c]) / dt; + if (dpvdt != 0.0 || rock_comp[c] != 0.0) { + rock_is_incomp = 0; + } + + h->J->sa[j] += porevol[c] * rock_comp[c] / dt; + h->F[c] -= dpvdt; + } + + /* Re-do the singularity-removing adjustment if necessary */ + if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) { + h->J->sa[0] *= 2; + } +} + + /* ---------------------------------------------------------------------- */ void cfs_tpfa_res_flux(struct UnstructuredGrid *G , diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 4fde3ba3..6a86c61b 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -72,6 +72,22 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G, const double *porevol, struct cfs_tpfa_res_data *h); +void +cfs_tpfa_res_comprock_assemble( + struct UnstructuredGrid *G, + double dt, + struct cfs_tpfa_res_forces *forces, + const double *zc, + struct compr_quantities_gen *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *wpress, + const double *porevol, + const double *porevol0, + const double *rock_comp, + struct cfs_tpfa_res_data *h); + void cfs_tpfa_res_flux(struct UnstructuredGrid *G , struct cfs_tpfa_res_forces *forces , From dec8ed5760abfab129d7706bc457ec6ca13ec116 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 10:38:41 +0200 Subject: [PATCH 02/31] Class CompressibleTpfa now handles rock compressibility. --- examples/sim_wateroil.cpp | 29 ++++++++----------- opm/core/pressure/CompressibleTpfa.cpp | 39 ++++++++++++++++++++++---- opm/core/pressure/CompressibleTpfa.hpp | 38 ++++++++++++++----------- 3 files changed, 67 insertions(+), 39 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 0ef94b05..fb43d98d 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -141,7 +141,6 @@ main(int argc, char** argv) std::cout << "--------------- Reading parameters ---------------" << std::endl; // Reading various control parameters. - const bool use_reorder = param.getDefault("use_reorder", true); const bool output = param.getDefault("output", true); std::string output_dir; int output_interval = 1; @@ -231,7 +230,7 @@ main(int argc, char** argv) bool use_segregation_split = false; bool use_column_solver = false; bool use_gauss_seidel_gravity = false; - if (use_gravity && use_reorder) { + if (use_gravity) { use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split); if (use_segregation_split) { use_column_solver = param.getDefault("use_column_solver", use_column_solver); @@ -241,18 +240,6 @@ 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; - 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; @@ -262,11 +249,11 @@ main(int argc, char** argv) // Extra rock init. std::vector porevol; if (rock_comp->isActive()) { - THROW("CompressibleTpfa solver does not handle this."); computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); } else { computePorevolume(*grid->c_grid(), props->porosity(), porevol); } + std::vector initial_porevol = porevol; double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); @@ -293,7 +280,7 @@ main(int argc, char** argv) const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0); const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20); const double *grav = use_gravity ? &gravity[0] : 0; - Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, linsolver, + Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver, nl_press_res_tol, nl_press_change_tol, nl_press_maxiter, grav, wells->c_wells()); // Reordering solver. @@ -409,6 +396,11 @@ main(int argc, char** argv) Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, wells->c_wells(), well_state.perfRates(), reorder_src); + // Compute new porevolumes after pressure solve, if necessary. + if (rock_comp->isActive()) { + initial_porevol = porevol; + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } // Solve transport. transport_timer.start(); double stepsize = simtimer.currentStepLength(); @@ -420,10 +412,11 @@ main(int argc, char** argv) // Note that for now we do not handle rock compressibility, // although the transport solver should be able to. reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0], - &porevol[0], &porevol[0], &reorder_src[0], stepsize, state.saturation()); + &porevol[0], &initial_porevol[0], &reorder_src[0], stepsize, state.saturation()); // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); if (use_segregation_split) { - reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation()); + reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0], + stepsize, grav, state.saturation()); } } transport_timer.stop(); diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 154f0633..dbd7ea47 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -58,6 +59,7 @@ namespace Opm /// to change. CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, const double residual_tol, const double change_tol, @@ -66,6 +68,7 @@ namespace Opm const struct Wells* wells) : grid_(grid), props_(props), + rock_comp_props_(rock_comp_props), linsolver_(linsolver), residual_tol_(residual_tol), change_tol_(change_tol), @@ -74,7 +77,6 @@ namespace Opm wells_(wells), htrans_(grid.cell_facepos[ grid.number_of_cells ]), trans_ (grid.number_of_faces), - porevol_(grid.number_of_cells), allcells_(grid.number_of_cells) { if (wells_ && (wells_->number_of_phases != props.numPhases())) { @@ -86,7 +88,12 @@ namespace Opm UnstructuredGrid* gg = const_cast(&grid_); tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]); tpfa_trans_compute(gg, &htrans_[0], &trans_[0]); - computePorevolume(grid_, props.porosity(), porevol_); + // If we have rock compressibility, pore volumes are updated + // in the compute*() methods, otherwise they are constant and + // hence may be computed here. + if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { + computePorevolume(grid_, props.porosity(), porevol_); + } for (int c = 0; c < grid.number_of_cells; ++c) { allcells_[c] = c; } @@ -230,6 +237,9 @@ namespace Opm const WellState& /*well_state*/) { computeWellPotentials(state); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_); + } } @@ -252,6 +262,8 @@ namespace Opm // std::vector face_gravcap_; // std::vector wellperf_A_; // std::vector wellperf_phasemob_; + // std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + // std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. computeCellDynamicData(dt, state, well_state); computeFaceDynamicData(dt, state, well_state); computeWellDynamicData(dt, state, well_state); @@ -273,6 +285,8 @@ namespace Opm // std::vector cell_viscosity_; // std::vector cell_phasemob_; // std::vector cell_voldisc_; + // std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + // std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. const int nc = grid_.number_of_cells; const int np = props_.numPhases(); const double* cell_p = &state.pressure()[0]; @@ -296,6 +310,14 @@ namespace Opm // TODO: Check this! cell_voldisc_.clear(); cell_voldisc_.resize(nc, 0.0); + + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_); + rock_comp_.resize(nc); + for (int cell = 0; cell < nc; ++cell) { + rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]); + } + } } @@ -465,9 +487,16 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; - cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], - &face_gravcap_[0], cell_press, well_bhp, - &porevol_[0], h_); + if (rock_comp_props_ == NULL) { + cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], h_); + } else { + cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], &initial_porevol_[0], + &rock_comp_[0], h_); + } } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 8233ee17..8a2857f2 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -33,6 +33,7 @@ namespace Opm class BlackoilState; class BlackoilPropertiesInterface; + class RockCompressibility; class LinearSolverInterface; class WellState; @@ -44,23 +45,25 @@ namespace Opm { public: /// Construct solver. - /// \param[in] grid A 2d or 3d grid. - /// \param[in] props Rock and fluid properties. - /// \param[in] linsolver Linear solver to use. - /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. - /// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller. - /// \param[in] maxiter Maximum acceptable number of iterations. - /// \param[in] gravity Gravity vector. If non-null, the array should - /// have D elements. - /// \param[in] wells The wells argument. Will be used in solution, - /// is ignored if NULL. - /// Note: this class observes the well object, and - /// makes the assumption that the well topology - /// and completions does not change during the - /// run. However, controls (only) are allowed - /// to change. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] props Rock and fluid properties. + /// \param[in] rock_comp_props Rock compressibility properties. May be null. + /// \param[in] linsolver Linear solver to use. + /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. + /// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller. + /// \param[in] maxiter Maximum acceptable number of iterations. + /// \param[in] gravity Gravity vector. If non-null, the array should + /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL. + /// Note: this class observes the well object, and + /// makes the assumption that the well topology + /// and completions does not change during the + /// run. However, controls (only) are allowed + /// to change. CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, const double residual_tol, const double change_tol, @@ -107,6 +110,7 @@ namespace Opm // ------ Data that will remain unmodified after construction. ------ const UnstructuredGrid& grid_; const BlackoilPropertiesInterface& props_; + const RockCompressibility* rock_comp_props_; const LinearSolverInterface& linsolver_; const double residual_tol_; const double change_tol_; @@ -115,7 +119,6 @@ namespace Opm const Wells* wells_; // May be NULL, outside may modify controls (only) between calls to solve(). std::vector htrans_; std::vector trans_ ; - std::vector porevol_; std::vector allcells_; // ------ Internal data for the cfs_tpfa_res solver. ------ @@ -123,6 +126,7 @@ namespace Opm // ------ Data that will be modified for every solve. ------ std::vector wellperf_gpot_; + std::vector initial_porevol_; // ------ Data that will be modified for every solver iteration. ------ std::vector cell_A_; @@ -135,6 +139,8 @@ namespace Opm std::vector face_gravcap_; std::vector wellperf_A_; std::vector wellperf_phasemob_; + std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. // The update to be applied to the pressures (cell and bhp). std::vector pressure_increment_; From 6e778b96414a142ad67c34a062b73956896fbbf8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 10:56:27 +0200 Subject: [PATCH 03/31] Bugfix: call *comprock* method only if active rock compressibility. --- 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 dbd7ea47..d28db5a3 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -487,7 +487,7 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; - if (rock_comp_props_ == NULL) { + if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], &face_gravcap_[0], cell_press, well_bhp, &porevol_[0], h_); From f15a6e03f09cbbda16604441dd5605d5a9829035 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 20 Aug 2012 15:14:01 +0200 Subject: [PATCH 04/31] Added common root name for vtk file names. --- opm/core/simulator/SimulatorTwophase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorTwophase.cpp index 84668b2a..d1c0ec3e 100644 --- a/opm/core/simulator/SimulatorTwophase.cpp +++ b/opm/core/simulator/SimulatorTwophase.cpp @@ -146,7 +146,7 @@ namespace Opm catch (...) { THROW("Creating directories failed: " << fpath); } - vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu"; + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); From 9e480cd997fc2ae1659d53ba10810f06d8540c3e Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 20 Aug 2012 15:18:24 +0200 Subject: [PATCH 05/31] Fix initialization bug. --- examples/sim_wateroil.cpp | 12 ++---------- .../reorder/TransportModelCompressibleTwophase.cpp | 1 + 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index fb43d98d..890ea981 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -228,16 +228,8 @@ main(int argc, char** argv) } } bool use_segregation_split = false; - bool use_column_solver = false; - bool use_gauss_seidel_gravity = false; if (use_gravity) { use_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); - } - } } // Source-related variables init. @@ -287,13 +279,13 @@ main(int argc, char** argv) const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const int nl_maxiter = param.getDefault("nl_maxiter", 30); Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); - if (use_gauss_seidel_gravity) { + if (use_segregation_split) { reorder_model.initGravity(); } // Column-based gravity segregation solver. std::vector > columns; - if (use_column_solver) { + if (use_segregation_split) { Opm::extractColumn(*grid->c_grid(), columns); } diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index f55f67ed..531aca4f 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -381,6 +381,7 @@ namespace Opm std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); const int nf = grid_.number_of_faces; trans_.resize(nf); + gravflux_.resize(nf); tpfa_htrans_compute(const_cast(&grid_), props_.permeability(), &htrans[0]); tpfa_trans_compute(const_cast(&grid_), &htrans[0], &trans_[0]); } From 6db144ec65d6d802d7230750bf1c2c22829e9e72 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 10:20:59 +0200 Subject: [PATCH 06/31] Fixed sign error and replaced porevol->porevol0. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index ab79dbce..b3ca1f9d 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1247,7 +1247,7 @@ cfs_tpfa_res_comprock_assemble( /* Assemble usual system (without rock compressibility). */ cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, - cpress, wpress, porevol, h); + cpress, wpress, porevol0, h); /* Check if we have only Neumann wells. */ well_is_neumann = 1; @@ -1276,7 +1276,7 @@ cfs_tpfa_res_comprock_assemble( } h->J->sa[j] += porevol[c] * rock_comp[c] / dt; - h->F[c] -= dpvdt; + h->F[c] += dpvdt; } /* Re-do the singularity-removing adjustment if necessary */ From c431df3d775c8afcd24703951c3f091160389de1 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 17:09:19 +0200 Subject: [PATCH 07/31] Corrected typo in comment. --- opm/core/pressure/tpfa/ifs_tpfa.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index 94e3d638..0e69bc49 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -771,7 +771,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G , assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok); /* We want to solve a Newton step for the residual - * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible + * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_incompressible * */ From 22fa63e6c4cc441e8df57fed4272c659ba8eced9 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 17:10:49 +0200 Subject: [PATCH 08/31] Changed way to get total number of cells. --- opm/core/fluid/RockFromDeck.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 6568e8cc..02f67556 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -87,7 +87,7 @@ namespace Opm double perm_threshold) { const int dim = 3; - const int num_global_cells = numGlobalCells(parser); + const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2]; const int nc = grid.number_of_cells; ASSERT (num_global_cells > 0); From c949319bf1fec7abb7498903868cf96ea9fee2fd Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 17:12:07 +0200 Subject: [PATCH 09/31] Corrected bug (residual should not be divided by dt). --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index b3ca1f9d..fb8ee879 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1242,7 +1242,7 @@ cfs_tpfa_res_comprock_assemble( int c, w, well_is_neumann, rock_is_incomp; size_t j; - double dpvdt; + double dpv; const struct Wells* W; /* Assemble usual system (without rock compressibility). */ @@ -1270,13 +1270,13 @@ cfs_tpfa_res_comprock_assemble( for (c = 0; c < G->number_of_cells; c++) { j = csrmatrix_elm_index(c, c, h->J); - dpvdt = (porevol[c] - porevol0[c]) / dt; - if (dpvdt != 0.0 || rock_comp[c] != 0.0) { + dpv = (porevol[c] - porevol0[c]); + if (dpv != 0.0 || rock_comp[c] != 0.0) { rock_is_incomp = 0; } - h->J->sa[j] += porevol[c] * rock_comp[c] / dt; - h->F[c] += dpvdt; + h->J->sa[j] += porevol[c] * rock_comp[c]; + h->F[c] += dpv; } /* Re-do the singularity-removing adjustment if necessary */ From d3641baf755803ebcc9f09dcee307c227058541f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 Aug 2012 09:05:54 +0200 Subject: [PATCH 10/31] Minor comment fix. --- opm/core/simulator/SimulatorIncompTwophase.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/simulator/SimulatorIncompTwophase.hpp b/opm/core/simulator/SimulatorIncompTwophase.hpp index ff4fbb68..10c46ff0 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.hpp +++ b/opm/core/simulator/SimulatorIncompTwophase.hpp @@ -96,4 +96,4 @@ namespace Opm } // namespace Opm -#endif // OPM_SIMULATORTWOPHASE_HEADER_INCLUDED +#endif // OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED From 6fd14bd76a40dc6410d4a3a5bf2e4dc4165ddce7 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 22 Aug 2012 10:11:32 +0200 Subject: [PATCH 11/31] Update surface volume in transport solver. --- examples/sim_wateroil.cpp | 5 +++-- .../TransportModelCompressibleTwophase.cpp | 15 ++++++++++++--- .../TransportModelCompressibleTwophase.hpp | 4 ++-- 3 files changed, 17 insertions(+), 7 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 890ea981..c1f26233 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -403,8 +403,9 @@ main(int argc, char** argv) for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { // Note that for now we do not handle rock compressibility, // although the transport solver should be able to. - reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0], - &porevol[0], &initial_porevol[0], &reorder_src[0], stepsize, state.saturation()); + reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], + &porevol[0], &initial_porevol[0], &reorder_src[0], stepsize, + state.saturation(), state.surfacevol()); // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); if (use_segregation_split) { reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0], diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 531aca4f..06708da7 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -75,15 +75,15 @@ namespace Opm void TransportModelCompressibleTwophase::solve(const double* darcyflux, const double* pressure, - const double* surfacevol0, const double* porevolume0, const double* porevolume, const double* source, const double dt, - std::vector& saturation) + std::vector& saturation, + std::vector& surfacevol) { darcyflux_ = darcyflux; - surfacevol0_ = surfacevol0; + surfacevol0_ = &surfacevol[0]; porevolume0_ = porevolume0; porevolume_ = porevolume; source_ = source; @@ -107,6 +107,15 @@ namespace Opm &ia_downw_[0], &ja_downw_[0]); reorderAndTransport(grid_, darcyflux); toBothSat(saturation_, saturation); + + // Compute surface volume as a postprocessing step from saturation and A_ + surfacevol = saturation; + const int np = props_.numPhases(); + for (int cell = 0; cell < grid_.number_of_cells; ++cell) { + for (int phase = 0; phase < np; ++phase) { + surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase]; + } + } } // Residual function r(s) for a single-cell implicit Euler transport diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 19c44b22..542363cd 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -54,12 +54,12 @@ namespace Opm /// \param[in, out] saturation Phase saturations. void solve(const double* darcyflux, const double* pressure, - const double* surfacevol0, const double* porevolume0, const double* porevolume, const double* source, const double dt, - std::vector& saturation); + std::vector& saturation, + std::vector& surfacevol); /// Initialise quantities needed by gravity solver. void initGravity(); From ae66fe74e31357fc3d453dca9cd0c52b67eef9e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 Aug 2012 10:37:52 +0200 Subject: [PATCH 12/31] New class SimulatorCompressibleTwophase. --- Makefile.am | 2 + .../SimulatorCompressibleTwophase.cpp | 564 ++++++++++++++++++ .../SimulatorCompressibleTwophase.hpp | 99 +++ 3 files changed, 665 insertions(+) create mode 100644 opm/core/simulator/SimulatorCompressibleTwophase.cpp create mode 100644 opm/core/simulator/SimulatorCompressibleTwophase.hpp diff --git a/Makefile.am b/Makefile.am index c1c85dec..712d7c65 100644 --- a/Makefile.am +++ b/Makefile.am @@ -92,6 +92,7 @@ opm/core/pressure/tpfa/compr_source.c \ opm/core/pressure/tpfa/ifs_tpfa.c \ opm/core/pressure/tpfa/trans_tpfa.c \ opm/core/pressure/well.c \ +opm/core/simulator/SimulatorCompressibleTwophase.cpp \ opm/core/simulator/SimulatorIncompTwophase.cpp \ opm/core/simulator/SimulatorReport.cpp \ opm/core/simulator/SimulatorTimer.cpp \ @@ -195,6 +196,7 @@ opm/core/pressure/tpfa/compr_source.h \ opm/core/pressure/tpfa/ifs_tpfa.h \ opm/core/pressure/tpfa/trans_tpfa.h \ opm/core/simulator/BlackoilState.hpp \ +opm/core/simulator/SimulatorCompressibleTwophase.hpp \ opm/core/simulator/SimulatorReport.hpp \ opm/core/simulator/SimulatorIncompTwophase.hpp \ opm/core/simulator/SimulatorTimer.hpp \ diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp new file mode 100644 index 00000000..52eb3e44 --- /dev/null +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -0,0 +1,564 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + + +namespace Opm +{ + + class SimulatorCompressibleTwophase::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity); + + SimulatorReport run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state); + + private: + // Data. + + // Parameters for output. + bool output_; + bool output_vtk_; + std::string output_dir_; + int output_interval_; + // Parameters for well control + bool check_well_controls_; + int max_well_control_iterations_; + // Parameters for transport solver. + int num_transport_substeps_; + bool use_segregation_split_; + // Observed objects. + const UnstructuredGrid& grid_; + const BlackoilPropertiesInterface& props_; + const RockCompressibility* rock_comp_; + WellsManager& wells_manager_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + const double* gravity_; + // Solvers + CompressibleTpfa psolver_; + TransportModelCompressibleTwophase tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorCompressibleTwophase::SimulatorCompressibleTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp, wells_manager, src, bcs, linsolver, gravity)); + } + + + + + + SimulatorReport SimulatorCompressibleTwophase::run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state) + { + return pimpl_->run(timer, state, well_state); + } + + + + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::BlackoilState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + } + + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::BlackoilState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + watercut.write(os); + } + + + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + wellreport.write(os); + } + + + static bool allNeumannBCs(const FlowBoundaryConditions* bcs) + { + if (bcs == NULL) { + return true; + } else { + return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE) + == bcs->type + bcs->nbc; + } + } + + + static bool allRateWells(const Wells* wells) + { + if (wells == NULL) { + return true; + } + const int nw = wells->number_of_wells; + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells->ctrls[w]; + if (wc->current >= 0) { + if (wc->type[wc->current] == BHP) { + return false; + } + } + } + return true; + } + + + + + // \TODO: make CompressibleTpfa take src and bcs. + SimulatorCompressibleTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_(rock_comp), + wells_manager_(wells_manager), + wells_(wells_manager.c_wells()), + src_(src), + bcs_(bcs), + gravity_(gravity), + psolver_(grid, props, /* rock_comp, */ linsolver, + param.getDefault("nl_pressure_residual_tolerance", 0.0), + param.getDefault("nl_pressure_change_tolerance", 1.0), + param.getDefault("nl_pressure_maxiter", 10), + gravity, wells_manager.c_wells() /*, src, bcs*/), + tsolver_(grid, props, + param.getDefault("nl_tolerance", 1e-9), + param.getDefault("nl_maxiter", 30)) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // Well control related init. + check_well_controls_ = param.getDefault("check_well_controls", false); + max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); + + // Transport related init. + num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); + use_segregation_split_ = param.getDefault("use_segregation_split", false); + if (gravity != 0 && use_segregation_split_){ + tsolver_.initGravity(); + extractColumn(grid_, columns_); + } + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + SimulatorReport SimulatorCompressibleTwophase::Impl::run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state) + { + std::vector transport_src; + + // Initialisation. + std::vector porevol; + if (rock_comp_ && rock_comp_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + std::vector initial_porevol = porevol; + + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + double init_satvol[2] = { 0.0 }; + double satvol[2] = { 0.0 }; + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); + std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init + << " " << init_satvol[1]/tot_porevol_init << std::endl; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + Opm::WellReport wellreport; + std::vector fractional_flows; + std::vector well_resflows_phase; + if (wells_) { + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + 0.0, well_state.bhp(), well_state.perfRates()); + } + std::fstream tstep_os; + if (output_) { + std::string filename = output_dir_ + "/step_timing.param"; + tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); + } + for (; !timer.done(); ++timer) { + // Report timestep and (optionally) write state to disk. + step_timer.start(); + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + } + + SimulatorReport sreport; + + // Solve pressure equation. + if (check_well_controls_) { + computeFractionalFlow(props_, allcells_, + state.pressure(), state.surfacevol(), state.saturation(), + fractional_flows); + wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver. + pressure_timer.start(); + std::vector initial_pressure = state.pressure(); + psolver_.solve(timer.currentStepLength(), state, well_state); + + // Renormalize pressure if rock is incompressible, and + // there are no pressure conditions (bcs or wells). + // It is deemed sufficient for now to renormalize + // using geometric volume instead of pore volume. + if ((rock_comp_ == NULL || !rock_comp_->isActive()) + && allNeumannBCs(bcs_) && allRateWells(wells_)) { + // Compute average pressures of previous and last + // step, and total volume. + double av_prev_press = 0.0; + double av_press = 0.0; + double tot_vol = 0.0; + const int num_cells = grid_.number_of_cells; + for (int cell = 0; cell < num_cells; ++cell) { + av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; + av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; + tot_vol += grid_.cell_volumes[cell]; + } + // Renormalization constant + const double ren_const = (av_prev_press - av_press)/tot_vol; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += ren_const; + } + const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += ren_const; + } + } + + // Stop timer and report. + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + sreport.pressure_time = pt; + + // Optionally, check if well controls are satisfied. + if (check_well_controls_) { + Opm::computePhaseFlowRatesPerWell(*wells_, + fractional_flows, + well_state.perfRates(), + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { + THROW("Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); + } + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; + } + } + } while (!well_control_passed); + + // Update pore volumes if rock is compressible. + if (rock_comp_ && rock_comp_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + } + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), transport_src); + + // Solve transport. + transport_timer.start(); + double stepsize = timer.currentStepLength(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0], + &porevol[0], &initial_porevol[0], &transport_src[0], stepsize, state.saturation()); + Opm::computeInjectedProduced(props_, + state.pressure(), state.surfacevol(), state.saturation(), + transport_src, stepsize, injected, produced); + if (gravity_ != 0 && use_segregation_split_) { + tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0], + stepsize, gravity_, state.saturation()); + } + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + sreport.transport_time = tt; + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.saturation(), satvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + std::cout.precision(5); + const int width = 18; + std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n"; + std::cout << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + std::cout << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + std::cout << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + std::cout << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + std::cout << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + std::cout << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; + std::cout << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::endl; + std::cout.precision(8); + + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + sreport.total_time = step_timer.secsSinceStart(); + if (output_) { + sreport.reportParam(tstep_os); + } + } + + if (output_) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } + tstep_os.close(); + } + + total_timer.stop(); + + SimulatorReport report; + report.pressure_time = ptime; + report.transport_time = ttime; + report.total_time = total_timer.secsSinceStart(); + return report; + } + + +} // namespace Opm diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.hpp b/opm/core/simulator/SimulatorCompressibleTwophase.hpp new file mode 100644 index 00000000..d38d74b2 --- /dev/null +++ b/opm/core/simulator/SimulatorCompressibleTwophase.hpp @@ -0,0 +1,99 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED +#define OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; +struct FlowBoundaryConditions; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class BlackoilPropertiesInterface; + class RockCompressibility; + class WellsManager; + class LinearSolverInterface; + class SimulatorTimer; + class BlackoilState; + class WellState; + struct SimulatorReport; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorCompressibleTwophase + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp if non-null, rock compressibility properties + /// \param[in] well_manager well manager, may manage no (null) wells + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorCompressibleTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED From f7dafe6c85cb1a319c18f2f402e263b8dbadd10d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 Aug 2012 11:16:51 +0200 Subject: [PATCH 13/31] Chase changes to some interfaces after merging. --- opm/core/simulator/SimulatorCompressibleTwophase.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index 52eb3e44..d0422edc 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -277,7 +277,7 @@ namespace Opm src_(src), bcs_(bcs), gravity_(gravity), - psolver_(grid, props, /* rock_comp, */ linsolver, + psolver_(grid, props, rock_comp, linsolver, param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_maxiter", 10), @@ -476,8 +476,9 @@ namespace Opm std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; } for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { - tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0], - &porevol[0], &initial_porevol[0], &transport_src[0], stepsize, state.saturation()); + tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], + &porevol[0], &initial_porevol[0], &transport_src[0], stepsize, + state.saturation(), state.surfacevol()); Opm::computeInjectedProduced(props_, state.pressure(), state.surfacevol(), state.saturation(), transport_src, stepsize, injected, produced); From a34fc775dd6a74709bbd615118ab44a7b295fa69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 Aug 2012 12:31:59 +0200 Subject: [PATCH 14/31] Document and check (primitive) for non-miscibility requirement. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 5 +++++ .../transport/reorder/TransportModelCompressibleTwophase.hpp | 2 ++ 2 files changed, 7 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 06708da7..98aeefdf 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -93,6 +93,11 @@ namespace Opm props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); + // Check non-miscibility requirement (only done for first cell). + if (A_[1] != 0.0 || A_[2] != 0.0) { + THROW("TransportModelCompressibleTwophase requires a property object without miscibility."); + } + std::vector seq(grid_.number_of_cells); std::vector comp(grid_.number_of_cells + 1); int ncomp; diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 542363cd..5b82ed62 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -30,6 +30,8 @@ namespace Opm class BlackoilPropertiesInterface; + /// Implements a reordering transport solver for compressible, + /// non-miscible two-phase flow. class TransportModelCompressibleTwophase : public TransportModelInterface { public: From 712f2ed16bb4321bb58cfec2bbf142f529f00871 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 22 Aug 2012 16:07:26 +0200 Subject: [PATCH 15/31] Added comments on how to run the script generate_doc_figures.py. --- generate_doc_figures.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/generate_doc_figures.py b/generate_doc_figures.py index ae42709d..8d4931ee 100755 --- a/generate_doc_figures.py +++ b/generate_doc_figures.py @@ -1,6 +1,24 @@ + +# This script generate the illustration pictures for the documentation. +# +# To run this script, you have to install paraview, see: +# +# http://www.paraview.org/paraview/resources/software.php +# +# Eventually, set up the paths (figure_path, tutorial_data_path, tutorial_path) according to your own installation. +# (The default values should be ok.) +# +# Make sure that pvpython is in your path of executables. +# +# Run the following command at the root of the directory where +# opm-core is installed (which also corresponds to the directory where +# generate_doc_figures is located): +# +# pvpython generate_doc_figures.py +# + from subprocess import call from paraview.simple import * -# from paraview import servermanager from os import remove, mkdir, curdir from os.path import join, isdir @@ -13,7 +31,6 @@ collected_garbage_file = [] if not isdir(figure_path): mkdir(figure_path) -# connection = servermanager.Connect() # tutorial 1 call(join(tutorial_path, "tutorial1")) From 234d125266a209d0360d5034f0603530492aa345 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 08:57:48 +0200 Subject: [PATCH 16/31] Added computeSurfacevol() method. The method is called by the reordering transport solver after computing new saturations in order to update the surface volumes. --- .../TransportModelCompressibleTwophase.cpp | 11 ++----- opm/core/utility/miscUtilitiesBlackoil.cpp | 30 +++++++++++++++++++ opm/core/utility/miscUtilitiesBlackoil.hpp | 19 ++++++++++++ 3 files changed, 52 insertions(+), 8 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 98aeefdf..b81a7e60 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -112,15 +113,9 @@ namespace Opm &ia_downw_[0], &ja_downw_[0]); reorderAndTransport(grid_, darcyflux); toBothSat(saturation_, saturation); - + // Compute surface volume as a postprocessing step from saturation and A_ - surfacevol = saturation; - const int np = props_.numPhases(); - for (int cell = 0; cell < grid_.number_of_cells; ++cell) { - for (int phase = 0; phase < np; ++phase) { - surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase]; - } - } + computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]); } // Residual function r(s) for a single-cell implicit Euler transport diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index af088e05..1c8bceee 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -220,5 +220,35 @@ namespace Opm } } + /// Computes the surface volume densities from saturations by the formula + /// z = A s + /// for a number of data points, where z is the surface volume density, + /// s is the saturation (both as column vectors) and A is the + /// phase-to-component relation matrix. + /// @param[in] n number of data points + /// @param[in] np number of phases, must be 2 or 3 + /// @param[in] A array containing n square matrices of size num_phases^2, + /// in Fortran ordering, typically the output of a call + /// to the matrix() method of a BlackoilProperties* class. + /// @param[in] saturation concatenated saturation values (for all P phases) + /// @param[out] surfacevol concatenated surface-volume values (for all P phases) + void computeSurfacevol(const int n, + const int np, + const double* A, + const double* saturation, + double* surfacevol) + { + // Note: since this is a simple matrix-vector product, it can + // be done by a BLAS call, but then we have to reorder the A + // matrix data. + std::fill(surfacevol, surfacevol + n*np, 0.0); + for (int i = 0; i < n; ++i) { + for (int row = 0; row < np; ++row) { + for (int col = 0; col < np; ++col) { + surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col]; + } + } + } + } } // namespace Opm diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 929d42ec..4a3d6448 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -112,6 +112,25 @@ namespace Opm const std::vector& s, std::vector& fractional_flows); + + /// Computes the surface volume densities from saturations by the formula + /// z = A s + /// for a number of data points, where z is the surface volume density, + /// s is the saturation (both as column vectors) and A is the + /// phase-to-component relation matrix. + /// @param[in] n number of data points + /// @param[in] np number of phases, must be 2 or 3 + /// @param[in] A array containing n square matrices of size num_phases, + /// in Fortran ordering, typically the output of a call + /// to the matrix() method of a BlackoilProperties* class. + /// @param[in] saturation concatenated saturation values (for all P phases) + /// @param[out] surfacevol concatenated surface-volume values (for all P phases) + void computeSurfacevol(const int n, + const int np, + const double* A, + const double* saturation, + double* surfacevol); + } // namespace Opm #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED From 82204f2b3b7c9645156984bee96db519d08e65d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 08:59:09 +0200 Subject: [PATCH 17/31] Whitespace cleanup. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 82 +++++++++++----------- opm/core/utility/miscUtilitiesBlackoil.hpp | 26 +++---- 2 files changed, 54 insertions(+), 54 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 1c8bceee..548aaa23 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -48,39 +48,39 @@ namespace Opm void computeInjectedProduced(const BlackoilPropertiesInterface& props, const std::vector& press, const std::vector& z, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced) + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); - } - std::fill(injected, injected + np, 0.0); - std::fill(produced, produced + np, 0.0); + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); std::vector visc(np); - std::vector mob(np); - for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; - const double* sat = &s[np*c]; - props.relperm(1, sat, &c, &mob[0], 0); + std::vector mob(np); + for (int c = 0; c < num_cells; ++c) { + if (src[c] > 0.0) { + injected[0] += src[c]*dt; + } else if (src[c] < 0.0) { + const double flux = -src[c]*dt; + const double* sat = &s[np*c]; + props.relperm(1, sat, &c, &mob[0], 0); props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); - double totmob = 0.0; - for (int p = 0; p < np; ++p) { - mob[p] /= visc[p]; - totmob += mob[p]; - } - for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; - } - } - } + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + mob[p] /= visc[p]; + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + } + } } @@ -93,11 +93,11 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& press, const std::vector& z, - const std::vector& s, - std::vector& totmob) + const std::vector& s, + std::vector& totmob) { std::vector pmobc; @@ -126,12 +126,12 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob, - std::vector& omega) + const std::vector& s, + std::vector& totmob, + std::vector& omega) { std::vector pmobc; @@ -185,10 +185,10 @@ namespace Opm props.relperm(nc, &s[0], &cells[0], &pmobc[0], dpmobc); - std::transform(pmobc.begin(), pmobc.end(), - mu.begin(), - pmobc.begin(), - std::divides()); + std::transform(pmobc.begin(), pmobc.end(), + mu.begin(), + pmobc.begin(), + std::divides()); } /// Computes the fractional flow for each cell in the cells argument diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 4a3d6448..565acc6e 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -46,11 +46,11 @@ namespace Opm void computeInjectedProduced(const BlackoilPropertiesInterface& props, const std::vector& p, const std::vector& z, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced); + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced); /// @brief Computes total mobility for a set of saturation values. /// @param[in] props rock and fluid properties @@ -60,11 +60,11 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob); + const std::vector& s, + std::vector& totmob); /// @brief Computes total mobility and omega for a set of saturation values. /// @param[in] props rock and fluid properties @@ -75,12 +75,12 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob, - std::vector& omega); + const std::vector& s, + std::vector& totmob, + std::vector& omega); /// @brief Computes phase mobilities for a set of saturation values. @@ -96,7 +96,7 @@ namespace Opm const std::vector& z, const std::vector& s, std::vector& pmobc); - + /// Computes the fractional flow for each cell in the cells argument /// @param[in] props rock and fluid properties From e5d0bb98c17e36ce2d06de0e14c1fecfa3cc376f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:00:04 +0200 Subject: [PATCH 18/31] Function cfs_tpfa_residual_assemble() and friends now return singularity flag. The singularity flag is true if there are no pressure conditions and no compressibility (so the absolute values of the pressure solution will be arbitrary). --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 42 ++++++++++------------ opm/core/pressure/tpfa/cfs_tpfa_residual.h | 13 +++++-- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index fb8ee879..eb8dc8f3 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , /* ---------------------------------------------------------------------- */ -void +int cfs_tpfa_res_assemble(struct UnstructuredGrid *G , double dt , struct cfs_tpfa_res_forces *forces , @@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , struct cfs_tpfa_res_data *h ) /* ---------------------------------------------------------------------- */ { - int res_is_neumann, well_is_neumann, c, np2; + int res_is_neumann, well_is_neumann, c, np2, singular; csrmatrix_zero( h->J); vector_zero (h->J->m, h->F); @@ -1207,14 +1207,17 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , assemble_sources(dt, forces->src, h); } - if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) { - h->J->sa[0] *= 2; + singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp; + if (singular) { + h->J->sa[0] *= 2.0; } + + return singular; } /* ---------------------------------------------------------------------- */ -void +int cfs_tpfa_res_comprock_assemble( struct UnstructuredGrid *G , double dt , @@ -1240,29 +1243,18 @@ cfs_tpfa_res_comprock_assemble( * porevol(pressure)*rock_comp(pressure)/dt */ - int c, w, well_is_neumann, rock_is_incomp; + int c, rock_is_incomp, singular; size_t j; - double dpv; - const struct Wells* W; + double dpv; /* Assemble usual system (without rock compressibility). */ - cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, - cpress, wpress, porevol0, h); - - /* Check if we have only Neumann wells. */ - well_is_neumann = 1; - W = forces->wells->W; - for (w = 0; well_is_neumann && w < W->number_of_wells; w++) { - if ((W->ctrls[w]->current >= 0) && /* OPEN? */ - (W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) { - well_is_neumann = 0; - } - } + singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, + cpress, wpress, porevol0, h); /* If we made a singularity-removing adjustment in the regular assembly, we undo it here. */ - if (well_is_neumann && h->pimpl->is_incomp) { - h->J->sa[0] /= 2; + if (singular) { + h->J->sa[0] /= 2.0; } /* Add new terms to residual and Jacobian. */ @@ -1280,9 +1272,11 @@ cfs_tpfa_res_comprock_assemble( } /* Re-do the singularity-removing adjustment if necessary */ - if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) { - h->J->sa[0] *= 2; + if (rock_is_incomp && singular) { + h->J->sa[0] *= 2.0; } + + return rock_is_incomp && singular; } diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 6a86c61b..5eddeed6 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , void cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); -void +/* Return value is 1 if the assembled matrix was adjusted to remove a + singularity. This happens if all fluids are incompressible and + there are no pressure conditions on wells or boundaries. + Otherwise return 0. */ +int cfs_tpfa_res_assemble(struct UnstructuredGrid *G, double dt, struct cfs_tpfa_res_forces *forces, @@ -72,7 +76,12 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G, const double *porevol, struct cfs_tpfa_res_data *h); -void +/* Return value is 1 if the assembled matrix was adjusted to remove a + singularity. This happens if all fluids are incompressible, the + rock is incompressible, and there are no pressure conditions on + wells or boundaries. + Otherwise return 0. */ +int cfs_tpfa_res_comprock_assemble( struct UnstructuredGrid *G, double dt, From b0522a499a03b0f10a3c7e376f866390a76237c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:03:20 +0200 Subject: [PATCH 19/31] New singularPressure() method propagates singularity information. --- opm/core/pressure/CompressibleTpfa.cpp | 36 ++++++++++++++++++++------ opm/core/pressure/CompressibleTpfa.hpp | 18 +++++++++---- 2 files changed, 41 insertions(+), 13 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index d28db5a3..67d3a651 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -77,7 +77,8 @@ namespace Opm wells_(wells), htrans_(grid.cell_facepos[ grid.number_of_cells ]), trans_ (grid.number_of_faces), - allcells_(grid.number_of_cells) + allcells_(grid.number_of_cells), + singular_(false) { if (wells_ && (wells_->number_of_phases != props.numPhases())) { THROW("Inconsistent number of phases specified (wells vs. props): " @@ -189,6 +190,21 @@ namespace Opm + /// @brief After solve(), was the resulting pressure singular. + /// Returns true if the pressure is singular in the following + /// sense: if everything is incompressible and there are no + /// pressure conditions, the absolute values of the pressure + /// solution are arbitrary. (But the differences in pressure + /// are significant.) + bool CompressibleTpfa::singularPressure() const + { + return singular_; + } + + + + + /// Compute well potentials. void CompressibleTpfa::computeWellPotentials(const BlackoilState& state) { @@ -487,16 +503,20 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; + int was_adjusted = 0; if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { - cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], - &face_gravcap_[0], cell_press, well_bhp, - &porevol_[0], h_); + was_adjusted = + cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], h_); } else { - cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], - &face_gravcap_[0], cell_press, well_bhp, - &porevol_[0], &initial_porevol_[0], - &rock_comp_[0], h_); + was_adjusted = + cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], &initial_porevol_[0], + &rock_comp_[0], h_); } + singular_ = (was_adjusted == 1); } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 8a2857f2..9fc5c43f 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -81,6 +81,14 @@ namespace Opm BlackoilState& state, WellState& well_state); + /// @brief After solve(), was the resulting pressure singular. + /// Returns true if the pressure is singular in the following + /// sense: if everything is incompressible and there are no + /// pressure conditions, the absolute values of the pressure + /// solution are arbitrary. (But the differences in pressure + /// are significant.) + bool singularPressure() const; + private: void computePerSolveDynamicData(const double dt, const BlackoilState& state, @@ -143,11 +151,11 @@ namespace Opm std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. // The update to be applied to the pressures (cell and bhp). std::vector pressure_increment_; - - - - - + // True if the matrix assembled would be singular but for the + // adjustment made in the cfs_*_assemble() calls. This happens + // if everything is incompressible and there are no pressure + // conditions. + bool singular_; }; } // namespace Opm From bab4f8cb57fc9d782cf21dcbfb2cfed7fe7c404a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:05:28 +0200 Subject: [PATCH 20/31] Fixed renormalization conditions. --- .../SimulatorCompressibleTwophase.cpp | 45 +++---------------- 1 file changed, 7 insertions(+), 38 deletions(-) diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index d0422edc..4c9d26c5 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -17,7 +17,6 @@ along with OPM. If not, see . */ - #if HAVE_CONFIG_H #include "config.h" #endif // HAVE_CONFIG_H @@ -153,7 +152,7 @@ namespace Opm catch (...) { THROW("Creating directories failed: " << fpath); } - vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu"; + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); @@ -228,36 +227,6 @@ namespace Opm } - static bool allNeumannBCs(const FlowBoundaryConditions* bcs) - { - if (bcs == NULL) { - return true; - } else { - return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE) - == bcs->type + bcs->nbc; - } - } - - - static bool allRateWells(const Wells* wells) - { - if (wells == NULL) { - return true; - } - const int nw = wells->number_of_wells; - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells->ctrls[w]; - if (wc->current >= 0) { - if (wc->type[wc->current] == BHP) { - return false; - } - } - } - return true; - } - - - // \TODO: make CompressibleTpfa take src and bcs. SimulatorCompressibleTwophase::Impl::Impl(const parameter::ParameterGroup& param, @@ -402,12 +371,12 @@ namespace Opm std::vector initial_pressure = state.pressure(); psolver_.solve(timer.currentStepLength(), state, well_state); - // Renormalize pressure if rock is incompressible, and - // there are no pressure conditions (bcs or wells). - // It is deemed sufficient for now to renormalize - // using geometric volume instead of pore volume. - if ((rock_comp_ == NULL || !rock_comp_->isActive()) - && allNeumannBCs(bcs_) && allRateWells(wells_)) { + // Renormalize pressure if both fluids and rock are + // incompressible, and there are no pressure + // conditions (bcs or wells). It is deemed sufficient + // for now to renormalize using geometric volume + // instead of pore volume. + if (psolver_.singularPressure()) { // Compute average pressures of previous and last // step, and total volume. double av_prev_press = 0.0; From 2030c6a735107f13111b20f89a80d77b5014e90b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:07:11 +0200 Subject: [PATCH 21/31] Added compressible 2-phase simulator, to replace sim_wateroil.cpp. --- examples/Makefile.am | 2 + examples/sim_2p_comp_reorder.cpp | 285 +++++++++++++++++++++++++++++++ 2 files changed, 287 insertions(+) create mode 100644 examples/sim_2p_comp_reorder.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 6d664bdf..b01f7674 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -17,6 +17,7 @@ $(BOOST_SYSTEM_LIB) noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ +sim_2p_comp_reorder \ sim_2p_incomp_reorder \ sim_wateroil \ wells_example @@ -29,6 +30,7 @@ wells_example # Please maintain sort order from "noinst_PROGRAMS". refine_wells_SOURCES = refine_wells.cpp +sim_2p_comp_reorder_SOURCES = sim_2p_comp_reorder.cpp sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp sim_wateroil_SOURCES = sim_wateroil.cpp wells_example_SOURCES = wells_example.cpp diff --git a/examples/sim_2p_comp_reorder.cpp b/examples/sim_2p_comp_reorder.cpp new file mode 100644 index 00000000..22e32c37 --- /dev/null +++ b/examples/sim_2p_comp_reorder.cpp @@ -0,0 +1,285 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for weakly compressible two-phase flow ===============\n\n"; + parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr rock_comp; + BlackoilState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + deck.reset(new EclipseGridParser(deck_filename)); + // Grid init + grid.reset(new GridManager(*deck)); + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid())); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(*deck)); + // Gravity. + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); + } + initBlackoilSurfvol(*grid->c_grid(), *props, state); + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + initBlackoilSurfvol(*grid->c_grid(), *props, state); + } + + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + const double *grav = use_gravity ? &gravity[0] : 0; + + // Initialising src + int num_cells = grid->c_grid()->number_of_cells; + std::vector src(num_cells, 0.0); + if (use_deck) { + // Do nothing, wells will be the driving force, not source terms. + } else { + // Compute pore volumes, in order to enable specifying injection rate + // terms of total pore volume. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + FlowBCManager bcs; + if (param.getDefault("use_pside", false)) { + int pside = param.get("pside"); + double pside_pressure = param.get("pside_pressure"); + bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + LinearSolverFactory linsolver(param); + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::ofstream epoch_os; + std::string output_dir; + if (output) { + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + std::string filename = output_dir + "/epoch_timing.param"; + epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out); + // open file to clean it. The file is appended to in SimulatorTwophase + filename = output_dir + "/step_timing.param"; + std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out); + step_os.close(); + param.writeParam(output_dir + "/simulation.param"); + } + + + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << " (number of epochs: " + << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + + SimulatorReport rep; + if (!use_deck) { + // Simple simulation without a deck. + WellsManager wells; // no wells. + SimulatorCompressibleTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); + SimulatorTimer simtimer; + simtimer.init(param); + warnIfUnusedParams(param); + WellState well_state; + well_state.init(0, state); + rep = simulator.run(simtimer, state, well_state); + } else { + // With a deck, we may have more epochs etc. + WellState well_state; + int step = 0; + SimulatorTimer simtimer; + // Use timer for last epoch to obtain total time. + deck->setCurrentEpoch(deck->numberOfEpochs() - 1); + simtimer.init(*deck); + const double total_time = simtimer.totalTime(); + for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { + // Set epoch index. + deck->setCurrentEpoch(epoch); + + // Update the timer. + if (deck->hasField("TSTEP")) { + simtimer.init(*deck); + } else { + if (epoch != 0) { + THROW("No TSTEP in deck for epoch " << epoch); + } + simtimer.init(param); + } + simtimer.setCurrentStepNum(step); + simtimer.setTotalTime(total_time); + + // Report on start of epoch. + std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" + << "\n (number of steps: " + << simtimer.numSteps() - step << ")\n\n" << std::flush; + + // Create new wells, well_state + WellsManager wells(*deck, *grid->c_grid(), props->permeability()); + // @@@ HACK: we should really make a new well state and + // properly transfer old well state to it every epoch, + // since number of wells may change etc. + if (epoch == 0) { + well_state.init(wells.c_wells(), state); + } + + // Create and run simulator. + SimulatorCompressibleTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); + if (epoch == 0) { + warnIfUnusedParams(param); + } + SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + if (output) { + epoch_rep.reportParam(epoch_os); + } + // Update total timing report and remember step number. + rep += epoch_rep; + step = simtimer.currentStepNum(); + } + } + + std::cout << "\n\n================ End of simulation ===============\n\n"; + rep.report(std::cout); + + if (output) { + std::string filename = output_dir + "/walltime.param"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + rep.reportParam(tot_os); + } + +} From 26149c30a2a88a95b56fe866ee62bdcd91b0fe05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:45:23 +0200 Subject: [PATCH 22/31] Fixed solveGravity(): now properly modifies surfacevolume. Also: - solveGravity() interface changed to take surface volume as a parameter, - gravity vector is now given in initGravity() instead of solveGravity(), for consistency with the incompressible solver. --- examples/sim_wateroil.cpp | 4 ++-- .../SimulatorCompressibleTwophase.cpp | 5 ++--- .../TransportModelCompressibleTwophase.cpp | 22 +++++++++++++------ .../TransportModelCompressibleTwophase.hpp | 13 ++++++----- 4 files changed, 27 insertions(+), 17 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index c1f26233..f7e0172f 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -280,7 +280,7 @@ main(int argc, char** argv) const int nl_maxiter = param.getDefault("nl_maxiter", 30); Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); if (use_segregation_split) { - reorder_model.initGravity(); + reorder_model.initGravity(grav); } // Column-based gravity segregation solver. @@ -409,7 +409,7 @@ main(int argc, char** argv) // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); if (use_segregation_split) { reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0], - stepsize, grav, state.saturation()); + stepsize, state.saturation(), state.surfacevol()); } } transport_timer.stop(); diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index 4c9d26c5..c5a77eb0 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -245,7 +245,6 @@ namespace Opm wells_(wells_manager.c_wells()), src_(src), bcs_(bcs), - gravity_(gravity), psolver_(grid, props, rock_comp, linsolver, param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), @@ -279,7 +278,7 @@ namespace Opm num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); use_segregation_split_ = param.getDefault("use_segregation_split", false); if (gravity != 0 && use_segregation_split_){ - tsolver_.initGravity(); + tsolver_.initGravity(gravity); extractColumn(grid_, columns_); } @@ -453,7 +452,7 @@ namespace Opm transport_src, stepsize, injected, produced); if (gravity_ != 0 && use_segregation_split_) { tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0], - stepsize, gravity_, state.saturation()); + stepsize, state.saturation(), state.surfacevol()); } } transport_timer.stop(); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index b81a7e60..c2d50327 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -39,7 +39,8 @@ namespace Opm typedef RegulaFalsi RootFinder; - TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid, + TransportModelCompressibleTwophase::TransportModelCompressibleTwophase( + const UnstructuredGrid& grid, const Opm::BlackoilPropertiesInterface& props, const double tol, const int maxit) @@ -52,6 +53,7 @@ namespace Opm dt_(0.0), saturation_(grid.number_of_cells, -1.0), fractionalflow_(grid.number_of_cells, -1.0), + gravity_(0), mob_(2*grid.number_of_cells, -1.0), ia_upw_(grid.number_of_cells + 1, -1), ja_upw_(grid.number_of_faces, -1), @@ -384,7 +386,7 @@ namespace Opm - void TransportModelCompressibleTwophase::initGravity() + void TransportModelCompressibleTwophase::initGravity(const double* grav) { // Set up transmissibilities. std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); @@ -393,12 +395,15 @@ namespace Opm gravflux_.resize(nf); tpfa_htrans_compute(const_cast(&grid_), props_.permeability(), &htrans[0]); tpfa_trans_compute(const_cast(&grid_), &htrans[0], &trans_[0]); + + // Remember gravity vector. + gravity_ = grav; } - void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav) + void TransportModelCompressibleTwophase::initGravityDynamic() { // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] @@ -420,7 +425,7 @@ namespace Opm for (int ci = 0; ci < 2; ++ci) { double gdz = 0.0; for (int d = 0; d < dim; ++d) { - gdz += grav[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]); + gdz += gravity_[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]); } gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]); } @@ -504,11 +509,11 @@ namespace Opm const double* pressure, const double* porevolume0, const double dt, - const double* grav, - std::vector& saturation) + std::vector& saturation, + std::vector& surfacevol) { // Assume that solve() has already been called, so that A_ is current. - initGravityDynamic(grav); + initGravityDynamic(); // Initialize mobilities. const int nc = grid_.number_of_cells; @@ -541,6 +546,9 @@ namespace Opm std::cout << "Gauss-Seidel column solver average iterations: " << double(num_iters)/double(columns.size()) << std::endl; toBothSat(saturation_, saturation); + + // Compute surface volume as a postprocessing step from saturation and A_ + computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]); } } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 5b82ed62..4ee93e0a 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -54,6 +54,7 @@ namespace Opm /// \param[in] source Transport source term. /// \param[in] dt Time step. /// \param[in, out] saturation Phase saturations. + /// \param[in, out] surfacevol Surface volume densities for each phase. void solve(const double* darcyflux, const double* pressure, const double* porevolume0, @@ -64,7 +65,8 @@ namespace Opm std::vector& surfacevol); /// Initialise quantities needed by gravity solver. - void initGravity(); + /// \param[in] grav Gravity vector + void initGravity(const double* grav); /// Solve for gravity segregation. /// This uses a column-wise nonlinear Gauss-Seidel approach. @@ -74,14 +76,14 @@ namespace Opm /// \param[in] columns Vector of cell-columns. /// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] dt Time step. - /// \param[in] grav Gravity vector. /// \param[in, out] saturation Phase saturations. + /// \param[in, out] surfacevol Surface volume densities for each phase. void solveGravity(const std::vector >& columns, const double* pressure, const double* porevolume0, const double dt, - const double* grav, - std::vector& saturation); + std::vector& saturation, + std::vector& surfacevol); private: virtual void solveSingleCell(const int cell); @@ -90,7 +92,7 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - void initGravityDynamic(const double* grav); + void initGravityDynamic(); private: const UnstructuredGrid& grid_; @@ -112,6 +114,7 @@ namespace Opm std::vector saturation_; // P (= num. phases) per cell std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell // For gravity segregation. + const double* gravity_; std::vector trans_; std::vector density_; std::vector gravflux_; From af5cc531ef0660240dea781214f4041b963616ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:58:32 +0200 Subject: [PATCH 23/31] Use the start-of-timestep pore volume in transport solver. This is to improve consistency with other solvers, and relates to the expression that is converted into a finite difference when discretising: (phi s) - (phi s)^0 = phi^0(s - s^0) + (phi - phi^0)s ^^^^^ The above marks the spot where we now use phi^0 instead of phi. --- opm/core/simulator/SimulatorIncompTwophase.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/opm/core/simulator/SimulatorIncompTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp index 76b2ac30..9de04f38 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -335,7 +335,7 @@ namespace Opm computePorevolume(grid_, props_.porosity(), porevol); } const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - + std::vector initial_porevol = porevol; // Main simulation loop. Opm::time::StopWatch pressure_timer; @@ -452,6 +452,7 @@ namespace Opm // Update pore volumes if rock is compressible. if (rock_comp_ && rock_comp_->isActive()) { + initial_porevol = porevol; computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); } @@ -467,11 +468,11 @@ namespace Opm std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; } for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { - tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], + tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], stepsize, state.saturation()); Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced); if (use_segregation_split_) { - tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation()); + tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation()); } } transport_timer.stop(); From 24a9b8b539540c8b183d8fa81e48fe6ddd7f944b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 08:18:59 +0200 Subject: [PATCH 24/31] Removed unneeded function numGlobalCells(). --- opm/core/fluid/RockFromDeck.cpp | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 02f67556..7ba8903f 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -36,8 +36,6 @@ namespace Opm PermeabilityKind fillTensor(const EclipseGridParser& parser, std::vector*>& tensor, std::tr1::array& kmap); - - int numGlobalCells(const EclipseGridParser& parser); } // anonymous namespace @@ -334,26 +332,6 @@ namespace Opm return kind; } - int numGlobalCells(const EclipseGridParser& parser) - { - int ngc = -1; - - if (parser.hasField("DIMENS")) { - const std::vector& - dims = parser.getIntegerValue("DIMENS"); - - ngc = dims[0] * dims[1] * dims[2]; - } - else if (parser.hasField("SPECGRID")) { - const SPECGRID& sgr = parser.getSPECGRID(); - - ngc = sgr.dimensions[ 0 ]; - ngc *= sgr.dimensions[ 1 ]; - ngc *= sgr.dimensions[ 2 ]; - } - - return ngc; - } } // anonymous namespace } // namespace Opm From 3f2d1773f5644c37521243d825992f0f6802edcd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 08:19:24 +0200 Subject: [PATCH 25/31] Require DIMENS or SPECGRID also for DXV, DYZ etc. grid specs. Also check that dimensions are consistent. --- opm/core/GridManager.cpp | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index f84d1c55..04d8cd04 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -166,6 +166,16 @@ namespace Opm // Construct tensor grid from deck. void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck) { + // Extract logical cartesian size. + std::vector dims; + if (deck.hasField("DIMENS")) { + dims = deck.getIntegerValue("DIMENS"); + } else if (deck.hasField("SPECGRID")) { + dims = deck.getSPECGRID().dimensions; + } else { + THROW("Deck must have either DIMENS or SPECGRID."); + } + // Extract coordinates (or offsets from top, in case of z). const std::vector& dxv = deck.getFloatingPointValue("DXV"); const std::vector& dyv = deck.getFloatingPointValue("DYV"); @@ -174,6 +184,17 @@ namespace Opm std::vector y = coordsFromDeltas(dyv); std::vector z = coordsFromDeltas(dzv); + // Check that number of cells given are consistent with DIMENS/SPECGRID. + if (dims[0] != int(dxv.size())) { + THROW("Number of DXV data points do not match DIMENS or SPECGRID."); + } + if (dims[1] != int(dyv.size())) { + THROW("Number of DYV data points do not match DIMENS or SPECGRID."); + } + if (dims[2] != int(dzv.size())) { + THROW("Number of DZV data points do not match DIMENS or SPECGRID."); + } + // Extract top corner depths, if available. const double* top_depths = 0; std::vector top_depths_vec; From 7d0b9fc3b039e905ab38046bac0e5485ed0b2cdb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 13:31:23 +0200 Subject: [PATCH 26/31] Minor code cleanup in TransportModelTwophase. --- .../transport/reorder/TransportModelTwophase.cpp | 14 +++++++++----- .../transport/reorder/TransportModelTwophase.hpp | 9 +++++---- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index b6451809..2740d065 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -52,8 +52,8 @@ namespace Opm source_(0), dt_(0.0), saturation_(grid.number_of_cells, -1.0), - reorder_iterations_(grid.number_of_cells, 0), fractionalflow_(grid.number_of_cells, -1.0), + reorder_iterations_(grid.number_of_cells, 0), mob_(2*grid.number_of_cells, -1.0) #ifdef EXPERIMENT_GAUSS_SEIDEL , ia_upw_(grid.number_of_cells + 1, -1), @@ -107,6 +107,13 @@ namespace Opm toBothSat(saturation_, saturation); } + + const std::vector& TransportModelTwophase::getReorderIterations() const + { + return reorder_iterations_; + } + + // Residual function r(s) for a single-cell implicit Euler transport // // r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) @@ -645,10 +652,7 @@ namespace Opm toBothSat(saturation_, saturation); } - void TransportModelTwophase::getReorderIterations() - { - return reorder_iterations_; - }; + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 358baaa3..3630ada8 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -74,10 +74,11 @@ namespace Opm const double* porevolume, const double dt, std::vector& saturation); - //// Return reorder iterations - //// + + //// Return the number of iterations used by the reordering solver. //// \param[out] vector of iteration per cell - const std::vector& getReorderIterations(){return reorder_iterations_;}; + const std::vector& getReorderIterations() const; + private: virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); @@ -86,7 +87,7 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - private: + private: const UnstructuredGrid& grid_; const IncompPropertiesInterface& props_; const double* visc_; From b7e9fbe9a50191bd3e5d8c03c90ead8870f2d036 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Fri, 24 Aug 2012 13:39:42 +0200 Subject: [PATCH 27/31] Corrected mistake from moving function calls. --- opm/core/transport/reorder/TransportModelTwophase.cpp | 2 +- opm/core/transport/reorder/TransportModelTwophase.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index b6451809..3346d28e 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -645,7 +645,7 @@ namespace Opm toBothSat(saturation_, saturation); } - void TransportModelTwophase::getReorderIterations() + const std::vector& TransportModelTwophase::getReorderIterations() { return reorder_iterations_; }; diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 358baaa3..64449f9f 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -77,7 +77,7 @@ namespace Opm //// Return reorder iterations //// //// \param[out] vector of iteration per cell - const std::vector& getReorderIterations(){return reorder_iterations_;}; + const std::vector& getReorderIterations(); private: virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); From c4608327cf9202b1755ecb4d2030b5087dc2b3f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 13:46:16 +0200 Subject: [PATCH 28/31] Improved README in some small ways. --- README | 72 ++++++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 52 insertions(+), 20 deletions(-) diff --git a/README b/README index 21eb98a2..54dbb531 100644 --- a/README +++ b/README @@ -1,41 +1,73 @@ -These are the release notes for opm-core - Open Porous Media Core Library +============================== + +These are release notes for opm-core. + CONTENT +------- + opm-core is the core library within OPM and contains the following -* Eclipse preprosessing +* Eclipse deck input and preprosessing * Fluid properties (basic PVT models and rock properties) -* Grid (generates different grids) -* Linear Algerbra (interface to different linear solvers) -* Pressure solvers (different discretization schemes, different flow models) -* Simulator (some basic examples of simulators based on sequential schemes) -* Transport solvers (different discretization schemes, different flow models) -* Utility (input and output processing, unit conversion) +* Grid handling (cornerpoint grids, unstructured grid interface) +* Linear Algebra (interface to different linear solvers) +* Pressure solvers (various discretization schemes, flow models) +* Simulators (some basic examples of simulators based on sequential splitting schemes) +* Transport solvers (various discretization schemes, flow models) +* Utilities (input and output processing, unit conversion) * Wells (basic well handling) -ON WHAT PLATFORMS DOES IT RUN? -The opm-core module is designed to run on linux platforms. No efforts have -been made to ensure that the code will compile and run on windows platforms. -DOCUMENTATION -Efforts have been made to document the code with Doxygen. +LICENSE +------- + +The library is distributed under the GNU General Public License, +version 3 or later (GPLv3+). + + +PLATFORMS +--------- + +The opm-core module is designed to run on Linux platforms. It is also +regularly run on Mac OS X. No efforts have been made to ensure that +the code will compile and run on windows platforms. + + +DOWNLOADING +----------- -DOWNLOAD opm-core git clone git://github.com/OPM/opm-core.git -BUILDING opm-core -cd ../opm-core - autoreconf -i + +BUILDING +-------- + + cd ../opm-core + autoreconf -i ./configure make sudo make install DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS +------------------------------------------- + +(to be updated) DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS -blas libblas3 lapack liblapack3 libboost libxml2 gcc automake autoconf git -doxygen umfpack +----------------------------------------- + +blas libblas3 lapack liblapack3 libboost libxml2 gcc automake autoconf +git doxygen umfpack + + +DOCUMENTATION +------------- + +Efforts have been made to document the code with Doxygen. +In order to build the documentation, enter the command +$ doxygen +in the topmost directory. From cd1edde45d3119a6825773831b37ad9243a79173 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 09:48:06 +0200 Subject: [PATCH 29/31] Minor revision, mostly whitespace cleanup and comments. --- examples/sim_2p_comp_reorder.cpp | 6 +++--- opm/core/pressure/CompressibleTpfa.cpp | 6 +++--- opm/core/pressure/CompressibleTpfa.hpp | 16 ++++++++-------- .../TransportModelCompressibleTwophase.cpp | 2 +- .../transport/reorder/TransportModelTwophase.hpp | 2 +- 5 files changed, 16 insertions(+), 16 deletions(-) diff --git a/examples/sim_2p_comp_reorder.cpp b/examples/sim_2p_comp_reorder.cpp index 22e32c37..9f2b516b 100644 --- a/examples/sim_2p_comp_reorder.cpp +++ b/examples/sim_2p_comp_reorder.cpp @@ -277,9 +277,9 @@ main(int argc, char** argv) rep.report(std::cout); if (output) { - std::string filename = output_dir + "/walltime.param"; - std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); - rep.reportParam(tot_os); + std::string filename = output_dir + "/walltime.param"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + rep.reportParam(tot_os); } } diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 67d3a651..954e07fa 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -504,13 +504,13 @@ namespace Opm cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; int was_adjusted = 0; - if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { - was_adjusted = + if (! (rock_comp_props_ && rock_comp_props_->isActive())) { + was_adjusted = cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], &face_gravcap_[0], cell_press, well_bhp, &porevol_[0], h_); } else { - was_adjusted = + was_adjusted = cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], &face_gravcap_[0], cell_press, well_bhp, &porevol_[0], &initial_porevol_[0], diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 9fc5c43f..054ac2e4 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -61,7 +61,7 @@ namespace Opm /// and completions does not change during the /// run. However, controls (only) are allowed /// to change. - CompressibleTpfa(const UnstructuredGrid& grid, + CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, @@ -71,8 +71,8 @@ namespace Opm const double* gravity, const Wells* wells); - /// Destructor. - ~CompressibleTpfa(); + /// Destructor. + ~CompressibleTpfa(); /// Solve the pressure equation by Newton-Raphson scheme. /// May throw an exception if the number of iterations @@ -112,11 +112,11 @@ namespace Opm void solveIncrement(); double residualNorm() const; double incrementNorm() const; - void computeResults(BlackoilState& state, + void computeResults(BlackoilState& state, WellState& well_state) const; // ------ Data that will remain unmodified after construction. ------ - const UnstructuredGrid& grid_; + const UnstructuredGrid& grid_; const BlackoilPropertiesInterface& props_; const RockCompressibility* rock_comp_props_; const LinearSolverInterface& linsolver_; @@ -125,12 +125,12 @@ namespace Opm const int maxiter_; const double* gravity_; // May be NULL const Wells* wells_; // May be NULL, outside may modify controls (only) between calls to solve(). - std::vector htrans_; - std::vector trans_ ; + std::vector htrans_; + std::vector trans_ ; std::vector allcells_; // ------ Internal data for the cfs_tpfa_res solver. ------ - struct cfs_tpfa_res_data* h_; + struct cfs_tpfa_res_data* h_; // ------ Data that will be modified for every solve. ------ std::vector wellperf_gpot_; diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index c2d50327..697f94a5 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -96,7 +96,7 @@ namespace Opm props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); - // Check non-miscibility requirement (only done for first cell). + // Check immiscibility requirement (only done for first cell). if (A_[1] != 0.0 || A_[2] != 0.0) { THROW("TransportModelCompressibleTwophase requires a property object without miscibility."); } diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 3630ada8..97386c8a 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -76,7 +76,7 @@ namespace Opm std::vector& saturation); //// Return the number of iterations used by the reordering solver. - //// \param[out] vector of iteration per cell + //// \return vector of iteration per cell const std::vector& getReorderIterations() const; private: From 8e5ef9ac0da5ebadccdc8c0396f3f43294efb928 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 11:19:22 +0200 Subject: [PATCH 30/31] Fixed bug in matrix multiplication (matrix has Fortran element order). --- opm/core/utility/miscUtilitiesBlackoil.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 548aaa23..74eb1c9c 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -245,7 +245,7 @@ namespace Opm for (int i = 0; i < n; ++i) { for (int row = 0; row < np; ++row) { for (int col = 0; col < np; ++col) { - surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col]; + surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col]; } } } From 80e54a93d07d11f62974276bc114649da3be4c52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 13:17:27 +0200 Subject: [PATCH 31/31] Switch loop ordering for better cache performance. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 74eb1c9c..6f4f7556 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -243,8 +243,8 @@ namespace Opm // matrix data. std::fill(surfacevol, surfacevol + n*np, 0.0); for (int i = 0; i < n; ++i) { - for (int row = 0; row < np; ++row) { - for (int col = 0; col < np; ++col) { + for (int col = 0; col < np; ++col) { + for (int row = 0; row < np; ++row) { surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col]; } }