From f5d33f4602f4a43e7f53efbfe8e3777f4dd204f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 26 Mar 2012 13:35:19 +0200 Subject: [PATCH 1/4] Do not require ACTNUM presence to construct grid. --- opm/core/GridManager.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index 721e0731..d941f393 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -36,7 +36,10 @@ namespace Opm // Extract data from deck. const std::vector& zcorn = deck.getFloatingPointValue("ZCORN"); const std::vector& coord = deck.getFloatingPointValue("COORD"); - const std::vector& actnum = deck.getIntegerValue("ACTNUM"); + const int* actnum = 0; + if (deck.hasField("ACTNUM")) { + actnum = &(deck.getIntegerValue("ACTNUM")[0]); + } std::vector dims; if (deck.hasField("DIMENS")) { dims = deck.getIntegerValue("DIMENS"); @@ -50,7 +53,7 @@ namespace Opm struct grdecl grdecl; grdecl.zcorn = &zcorn[0]; grdecl.coord = &coord[0]; - grdecl.actnum = actnum.empty() ? 0 : &actnum[0]; + grdecl.actnum = actnum; grdecl.dims[0] = dims[0]; grdecl.dims[1] = dims[1]; grdecl.dims[2] = dims[2]; From 97c3a9d9e3904cd6952764d313006a1de332a2c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 26 Mar 2012 14:19:30 +0200 Subject: [PATCH 2/4] Gravity column solver output: make it clear where the output comes from. --- opm/core/transport/GravityColumnSolver_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/GravityColumnSolver_impl.hpp b/opm/core/transport/GravityColumnSolver_impl.hpp index 0bb32863..941da0b4 100644 --- a/opm/core/transport/GravityColumnSolver_impl.hpp +++ b/opm/core/transport/GravityColumnSolver_impl.hpp @@ -96,7 +96,7 @@ namespace Opm const double maxelem = *std::max_element(increment.begin(), increment.end()); const double minelem = *std::min_element(increment.begin(), increment.end()); max_delta = std::max(maxelem, -minelem); - std::cout << "Iteration " << iter << " max_delta = " << max_delta << std::endl; + std::cout << "Gravity column solver iteration " << iter << " max_delta = " << max_delta << std::endl; if (max_delta < tol_) { break; } From f406487f72f3fbbc7a877013ed803f806d65c995 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 26 Mar 2012 14:20:00 +0200 Subject: [PATCH 3/4] Bugfix in solveMultiCell(). --- opm/core/transport/reorder/TransportModelTwophase.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index 7d9aadac..70d9876e 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -362,7 +362,9 @@ namespace Opm for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { const int downwind_cell = ja_downw_[j]; int ci = pos[downwind_cell]; - needs_update[ci] = 1; + if (ci != -1) { + needs_update[ci] = 1; + } // ++needs_update[ci]; // if (needs_update[ci] == num_upstream[ci]) { // fully_marked_stack.push_back(ci); From bb953a5fba7170250cbaac68c23ff99b2189af16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 26 Mar 2012 16:08:57 +0200 Subject: [PATCH 4/4] Added controls to nonlinear pressure iteration. Fixed segregation scenario init. --- examples/spu_2p.cpp | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index f5d7508c..7920dc20 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -453,6 +453,8 @@ 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()) { if (!use_reorder) { THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet."); @@ -462,6 +464,8 @@ main(int argc, char** argv) THROW("For gravity segregation splitting, only use_gauss_seidel_gravity=true supports rock compressibility."); } } + nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10); + nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal } // State-related and source-related variables init. @@ -554,16 +558,18 @@ main(int argc, char** argv) if (grid->c_grid()->cartdims[2] <= 1) { std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl; } - std::vector& sat = state.saturation(); + std::vector top_cells; const int *glob_cell = grid->c_grid()->global_cell; // Water on top for (int cell = 0; cell < num_cells; ++cell) { const int* cd = grid->c_grid()->cartdims; const int gc = glob_cell == 0 ? cell : glob_cell[cell]; bool top = (gc / cd[0] / cd[1]) < cd[2]/2; - sat[2*cell] = top ? 1.0 : 0.0; - sat[2*cell + 1 ] = 1.0 - sat[2*cell]; + if (top) { + top_cells.push_back(cell); + } } + state.setWaterSat(top_cells, *props, ReservoirState::MaxSat); break; } case 4: @@ -618,9 +624,9 @@ main(int argc, char** argv) const double *grav = use_gravity ? &gravity[0] : 0; Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver); // Reordering solver. - const double nltol = param.getDefault("nl_tolerance", 1e-9); - const int maxit = param.getDefault("nl_maxiter", 30); - Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nltol, maxit); + const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); + const int nl_maxiter = param.getDefault("nl_maxiter", 30); + Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); if (use_gauss_seidel_gravity) { reorder_model.initGravity(grav); } @@ -636,7 +642,7 @@ main(int argc, char** argv) if (use_column_solver) { Opm::extractColumn(*grid->c_grid(), columns); } - Opm::GravityColumnSolver colsolver(model, *grid->c_grid(), nltol, maxit); + Opm::GravityColumnSolver colsolver(model, *grid->c_grid(), nl_tolerance, nl_maxiter); // Control init. Opm::ImplicitTransportDetails::NRReport rpt; @@ -707,8 +713,7 @@ main(int argc, char** argv) rc.resize(num_cells); std::vector initial_pressure = state.pressure(); std::vector prev_pressure; - const int num_pressure_iter = 10; - for (int iter = 0; iter < num_pressure_iter; ++iter) { + for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { prev_pressure = state.pressure(); for (int cell = 0; cell < num_cells; ++cell) { rc[cell] = rock_comp->rockComp(state.pressure()[cell]); @@ -720,7 +725,10 @@ main(int argc, char** argv) for (int cell = 0; cell < num_cells; ++cell) { max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell])); } - std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance) { + break; + } } computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); } else {