From f4b25ea27ffaab84df146ebfdad37d27bf976a5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 5 Jan 2023 15:13:19 +0100 Subject: [PATCH 01/20] Add SubDomain class, aspinPartition() function and matrix utilities. --- CMakeLists_files.cmake | 6 + opm/simulators/flow/SubDomain.hpp | 45 +++++++ opm/simulators/flow/aspinPartition.cpp | 93 +++++++++++++++ opm/simulators/flow/aspinPartition.hpp | 34 ++++++ opm/simulators/linalg/extractMatrix.hpp | 151 ++++++++++++++++++++++++ tests/partition.txt | 10 ++ tests/test_aspinPartition.cpp | 35 ++++++ 7 files changed, 374 insertions(+) create mode 100644 opm/simulators/flow/SubDomain.hpp create mode 100644 opm/simulators/flow/aspinPartition.cpp create mode 100644 opm/simulators/flow/aspinPartition.hpp create mode 100644 opm/simulators/linalg/extractMatrix.hpp create mode 100644 tests/partition.txt create mode 100644 tests/test_aspinPartition.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2ee0be1c437..4ead8403e2f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -50,6 +50,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/NonlinearSolverEbos.cpp opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp opm/simulators/flow/ValidationFunctions.cpp + opm/simulators/flow/aspinPartition.cpp opm/simulators/linalg/bda/WellContributions.cpp opm/simulators/linalg/bda/MultisegmentWellContribution.cpp opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp @@ -226,6 +227,7 @@ endif() # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort list (APPEND TEST_SOURCE_FILES tests/test_ALQState.cpp + tests/test_aspinPartition.cpp tests/test_blackoil_amg.cpp tests/test_convergenceoutputconfiguration.cpp tests/test_convergencereport.cpp @@ -343,6 +345,7 @@ list (APPEND TEST_DATA_FILES tests/include/summary.inc tests/include/test1_20x30x10.grdecl tests/include/well_vfp.ecl + tests/partition.txt ) @@ -362,6 +365,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp opm/simulators/flow/KeywordValidation.hpp opm/simulators/flow/ValidationFunctions.hpp + opm/simulators/flow/aspinPartition.hpp + opm/simulators/flow/SubDomain.hpp opm/core/props/BlackoilPhases.hpp opm/core/props/phaseUsageFromDeck.hpp opm/core/props/satfunc/RelpermDiagnostics.hpp @@ -421,6 +426,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/SmallDenseMatrixUtils.hpp opm/simulators/linalg/WellOperators.hpp opm/simulators/linalg/WriteSystemMatrixHelper.hpp + opm/simulators/linalg/extractMatrix.hpp opm/simulators/linalg/findOverlapRowsAndColumns.hpp opm/simulators/linalg/getQuasiImpesWeights.hpp opm/simulators/linalg/setupPropertyTree.hpp diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp new file mode 100644 index 00000000000..f7aeb21eb6a --- /dev/null +++ b/opm/simulators/flow/SubDomain.hpp @@ -0,0 +1,45 @@ +/* + Copyright 2021 Total SE + + 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_SUBDOMAIN_HEADER_INCLUDED +#define OPM_SUBDOMAIN_HEADER_INCLUDED + +#include + +#include + +namespace Opm +{ + + template + struct SubDomain + { + int index; + std::vector cells; + std::vector interior; + Dune::SubGridView view; + SubDomain(const int i, std::vector&& c, std::vector&& in, Dune::SubGridView&& v) + : index(i), cells(std::move(c)), interior(std::move(in)), view(std::move(v)) + {} + }; + +} // namespace Opm + + +#endif // OPM_SUBDOMAIN_HEADER_INCLUDED diff --git a/opm/simulators/flow/aspinPartition.cpp b/opm/simulators/flow/aspinPartition.cpp new file mode 100644 index 00000000000..05dfb8096ab --- /dev/null +++ b/opm/simulators/flow/aspinPartition.cpp @@ -0,0 +1,93 @@ +/* + Copyright 2021 Total SE + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + +namespace +{ + + // Read from file, containing one number per cell, from [0, ... , num_domains - 1]. + std::pair, int> partitionCellsFromFile([[maybe_unused]] const int num_cells) + { + // TODO: refactor to make more flexible. + // Read file into single vector. + const std::string filename = "partition.txt"; + std::ifstream is(filename); + const std::vector cellpart{std::istream_iterator(is), std::istream_iterator()}; + if (cellpart.size() != size_t(num_cells)) { + auto msg = fmt::format("Partition file contains {} entries, but there are {} cells.", + cellpart.size(), num_cells); + throw std::runtime_error(msg); + } + + // Create and return the output domain vector. + const int num_domains = (*std::max_element(cellpart.begin(), cellpart.end())) + 1; + return { cellpart, num_domains }; + } + + + // Trivially simple partitioner + std::pair, int> partitionCellsSimple(const int num_cells, const int num_domains) + { + // Build the partitions. + std::vector bounds(num_domains + 1); + bounds[0] = 0; + const int dom_sz = num_cells / num_domains; + for (int i = 1; i < num_domains; ++i) { + bounds[i] = dom_sz * i; + } + bounds[num_domains] = num_cells; + std::vector part(num_cells); + for (int i = 0; i < num_domains; ++i) { + std::fill(part.begin() + bounds[i], part.begin() + bounds[i + 1], i); + } + return { part, num_domains }; + } + +} // anonymous namespace + + +std::pair, int> +partitionCells(const int num_cells) +{ + // const std::string method = "simple"; + const std::string method = "file"; + if (method == "simple") { + return partitionCellsSimple(num_cells, 1); + } else if (method == "file") { + return partitionCellsFromFile(num_cells); + } else { + return {}; + } +} + +} // namespace Opm diff --git a/opm/simulators/flow/aspinPartition.hpp b/opm/simulators/flow/aspinPartition.hpp new file mode 100644 index 00000000000..f1781a08696 --- /dev/null +++ b/opm/simulators/flow/aspinPartition.hpp @@ -0,0 +1,34 @@ +/* + Copyright 2021 Total SE + + 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_ASPINPARTITION_HEADER_INCLUDED +#define OPM_ASPINPARTITION_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + +std::pair, int> partitionCells(const int num_cells); + +} // namespace Opm + + +#endif // OPM_ASPINPARTITION_HEADER_INCLUDED diff --git a/opm/simulators/linalg/extractMatrix.hpp b/opm/simulators/linalg/extractMatrix.hpp new file mode 100644 index 00000000000..b659f2e28b3 --- /dev/null +++ b/opm/simulators/linalg/extractMatrix.hpp @@ -0,0 +1,151 @@ +/* + Copyright 2021 SINTEF Digital, Mathematics and Cybernetics. + + 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_EXTRACT_MATRIX_HEADER_INCLUDED +#define OPM_EXTRACT_MATRIX_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + +namespace Details +{ + + template + void copySubMatrix(const Matrix& A, const std::vector& indices, Matrix& B) + { + // Copy elements so that B_{i, j} = A_{indices[i], indices[j]}. + for (auto row = B.begin(); row != B.end(); ++row) { + for (auto col = row->begin(); col != row->end(); ++col) { + *col = A[indices[row.index()]][indices[col.index()]]; + } + } + } + + + template + Matrix extractMatrix(const Matrix& m, const std::vector& indices) + { + assert(std::is_sorted(indices.begin(), indices.end())); + + // Set up reverse index map. + const size_t n = indices.size(); + size_t newrow = 0; + enum { NotIncluded = -1 }; + std::vector index_map(m.N(), NotIncluded); + for (auto row = m.begin(); row != m.end(); ++row) { + if (row.index() == indices[newrow]) { + index_map[row.index()] = newrow; + ++newrow; + if (newrow == n) { + break; + } + } + } + assert(newrow == n); + + // Count nonzeroes. + size_t nnz = 0; + for (auto row = m.begin(); row != m.end(); ++row) { + if (index_map[row.index()] != NotIncluded) { + for (auto col = row->begin(); col != row->end(); ++col) { + if (index_map[col.index()] != NotIncluded) { + ++nnz; + } + } + } + } + + // Create the matrix structure. + Matrix res(n, n, nnz, Matrix::row_wise); + auto from_row = m.begin(); + for (auto row = res.createbegin(); row != res.createend(); ++row) { + // Move from_row to point to the row to be extracted. + while (from_row.index() < indices[row.index()]) { + ++from_row; + } + assert(from_row.index() == indices[row.index()]); + // Insert nonzeros for row. + for (auto from_col = from_row->begin(); from_col != from_row->end(); ++from_col) { + const int new_col = index_map[from_col.index()]; + if (new_col != NotIncluded) { + row.insert(new_col); + } + } + } + + copySubMatrix(m, indices, res); + return res; + } + + + template + Vector extractVector(const Vector& x, const std::vector& indices) + { + Vector res(indices.size()); + for (size_t ii = 0; ii < indices.size(); ++ii) { + res[ii] = x[indices[ii]]; + } + return res; + } + + + template + void setGlobal(const Vector& x, const std::vector& indices, Vector& global_x) + { + for (size_t ii = 0; ii < indices.size(); ++ii) { + global_x[indices[ii]] = x[ii]; + } + } + + + + template + bool matrixEqual(const Matrix& m1, const Matrix& m2) + { + // Compare size and nonzeroes. + if (m1.N() != m2.N()) return false; + if (m1.M() != m2.M()) return false; + if (m1.nonzeroes() != m2.nonzeroes()) return false; + + auto row1 = m1.begin(); + auto row2 = m2.begin(); + for (; row1 != m1.end(); ++row1, ++row2) { + if (row2 == m2.end()) return false; + if (row1.index() != row2.index()) return false; + auto col1 = row1->begin(); + auto col2 = row2->begin(); + for (; col1 != row1->end(); ++col1, ++col2) { + if (col2 == row2->end()) return false; + if (col1.index() != col2.index()) return false; + if (*col1 != *col2) return false; + } + } + return true; + } + + +} // namespace Details + + +} // namespace Opm + +#endif // OPM_EXTRACT_MATRIX_HEADER_INCLUDED diff --git a/tests/partition.txt b/tests/partition.txt new file mode 100644 index 00000000000..86d2c9c5f94 --- /dev/null +++ b/tests/partition.txt @@ -0,0 +1,10 @@ +0 +0 +1 +1 +2 +2 +1 +1 +0 +0 diff --git a/tests/test_aspinPartition.cpp b/tests/test_aspinPartition.cpp new file mode 100644 index 00000000000..de542d81c52 --- /dev/null +++ b/tests/test_aspinPartition.cpp @@ -0,0 +1,35 @@ +/* + Copyright 2021 SINTEF Digital, Mathematics and Cybernetics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#define BOOST_TEST_MODULE OPM_test_aspinPartition +#include + +#include + + +BOOST_AUTO_TEST_CASE(fileBased) +{ + auto [part, num_part] = Opm::partitionCells(10); + BOOST_CHECK_EQUAL(num_part, 3); + std::vector expected = { 0, 0, 1, 1, 2, 2, 1, 1, 0, 0 }; + BOOST_CHECK_EQUAL_COLLECTIONS(expected.begin(), expected.end(), part.begin(), part.end()); +} + From fc933240e797b763e1235b8d93a344566aa427ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 5 Jan 2023 15:15:14 +0100 Subject: [PATCH 02/20] Add new parameters for aspin usage. --- .../flow/BlackoilModelParametersEbos.hpp | 78 +++++++++++++++++-- 1 file changed, 73 insertions(+), 5 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 1da083a476c..5e8b3e7234d 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -168,7 +168,30 @@ template struct UseAverageDensityMsWells { using type = UndefinedProperty; }; - +template +struct NonlinearSolver { + using type = UndefinedProperty; +}; +template +struct LocalSolveApproach { + using type = UndefinedProperty; +}; +template +struct OuterAspinTolerance { + using type = UndefinedProperty; +}; +template +struct MaxLocalSolveIterations { + using type = UndefinedProperty; +}; +template +struct LocalToleranceScalingMb { + using type = UndefinedProperty; +}; +template +struct LocalToleranceScalingCnv { + using type = UndefinedProperty; +}; template struct DbhpMaxRel { using type = GetPropType; @@ -316,10 +339,34 @@ template struct UseAverageDensityMsWells { static constexpr bool value = false; }; - - - - +template +struct NonlinearSolver { + static constexpr auto value = "newton"; +}; +template +struct LocalSolveApproach { + static constexpr auto value = "jacobi"; +}; +template +struct OuterAspinTolerance { + using type = GetPropType; + static constexpr type value = 1e-3; +}; +template +struct MaxLocalSolveIterations { + static constexpr int value = 20; +}; +template +struct LocalToleranceScalingMb { + using type = GetPropType; + static constexpr type value = 1.0; +}; +template +struct LocalToleranceScalingCnv { + using type = GetPropType; + static constexpr type value = 0.01; +}; +>>>>>>> d1d200f5b (Add new parameters for aspin usage.) // if openMP is available, determine the number threads per process automatically. #if _OPENMP @@ -436,7 +483,16 @@ namespace Opm /// Whether to approximate segment densities by averaging over segment and its outlet bool use_average_density_ms_wells_; + // Choose nonlinear solver type: newton, aspin, nldd. + std::string nonlinear_solver_; + std::string local_solve_approach_; + + double outer_aspin_tolerance_; + + int max_local_solve_iterations_; + double local_tolerance_scaling_mb_; + double local_tolerance_scaling_cnv_; /// Construct from user parameters or defaults. BlackoilModelParametersEbos() @@ -473,6 +529,12 @@ namespace Opm check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter); max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches); use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells); + nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver); + local_solve_approach_ = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach); + outer_aspin_tolerance_ = EWOMS_GET_PARAM(TypeTag, double, OuterAspinTolerance); + max_local_solve_iterations_ = EWOMS_GET_PARAM(TypeTag, int, MaxLocalSolveIterations); + local_tolerance_scaling_mb_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingMb); + local_tolerance_scaling_cnv_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingCnv); deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); } @@ -512,6 +574,12 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations"); EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, NonlinearSolver, "Choose nonlinear solver. Valid choices are newton, aspin or nldd."); + EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalSolveApproach, "Choose local solve approach. Valid choices are jacobi and gauss-seidel"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, OuterAspinTolerance, "Tolerance for ASPIN residual."); + EWOMS_REGISTER_PARAM(TypeTag, int, MaxLocalSolveIterations, "Max iterations for local solves with ASPIN or NLDD."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingMb, "Set lower than 1.0 to use stricter convergence tolerance for local solves."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingCnv, "Set lower than 1.0 to use stricter convergence tolerance for local solves."); } }; } // namespace Opm From 497b9cf86bb4b59b690a28c1c0b2ee9aa244e56b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 5 Jan 2023 15:17:19 +0100 Subject: [PATCH 03/20] Store solver class instead of recreating. --- .../SimulatorFullyImplicitBlackoilEbos.hpp | 20 ++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp index 162f57e2131..efbc724b3c8 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp @@ -362,7 +362,9 @@ class SimulatorFullyImplicitBlackoilEbos // Run a multiple steps of the solver depending on the time step control. solverTimer_->start(); - auto solver = createSolver(wellModel_()); + if (!solver_) { + solver_ = createSolver(wellModel_()); + } ebosSimulator_.startNextEpisode( ebosSimulator_.startTime() @@ -375,8 +377,7 @@ class SimulatorFullyImplicitBlackoilEbos loadStep_ = -1; ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); } - solver->model().beginReportStep(); - + solver_->model().beginReportStep(); bool enableTUNING = EWOMS_GET_PARAM(TypeTag, bool, EnableTuning); // If sub stepping is enabled allow the solver to sub cycle @@ -398,13 +399,13 @@ class SimulatorFullyImplicitBlackoilEbos events.hasEvent(ScheduleEvents::INJECTION_TYPE_CHANGED) || events.hasEvent(ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER) || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE); - auto stepReport = adaptiveTimeStepping_->step(timer, *solver, event, nullptr); + auto stepReport = adaptiveTimeStepping_->step(timer, *solver_, event, nullptr); report_ += stepReport; //Pass simulation report to eclwriter for summary output ebosSimulator_.problem().setSimulationReport(report_); } else { // solve for complete report step - auto stepReport = solver->step(timer); + auto stepReport = solver_->step(timer); report_ += stepReport; if (terminalOutput_) { std::ostringstream ss; @@ -421,7 +422,7 @@ class SimulatorFullyImplicitBlackoilEbos ebosSimulator_.problem().writeOutput(); report_.success.output_write_time += perfTimer.stop(); - solver->model().endReportStep(); + solver_->model().endReportStep(); // take time that was used to solve system for this reportStep solverTimer_->stop(); @@ -433,7 +434,7 @@ class SimulatorFullyImplicitBlackoilEbos // Destructively grab the step convergence reports. The solver // object and the model object contained therein are about to go // out of scope. - this->writeConvergenceOutput(solver->model().getStepReportsDestructively()); + this->writeConvergenceOutput(solver_->model().getStepReportsDestructively()); } // Increment timer, remember well state. @@ -486,6 +487,9 @@ class SimulatorFullyImplicitBlackoilEbos serializer(adaptiveTimeStepping_); } + const Model& model() const + { return solver_->model(); } + protected: std::unique_ptr createSolver(WellModel& wellModel) @@ -702,6 +706,8 @@ class SimulatorFullyImplicitBlackoilEbos ModelParameters modelParam_; SolverParameters solverParam_; + std::unique_ptr solver_; + // Observed objects. PhaseUsage phaseUsage_; // Misc. data From 2e733ca0d8ccff71b6f3601dc730b33d31376ec7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 6 Jan 2023 11:36:30 +0100 Subject: [PATCH 04/20] HACK: ensure we rebuild UMFPACK solver always. Should probably be handled with fixing update() for preconditioner wrapper instead, making a new wrapper for those preconditioners that need to rebuild. --- opm/simulators/linalg/FlexibleSolver.hpp | 1 + opm/simulators/linalg/FlexibleSolver_impl.hpp | 16 ++++++++++++++-- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp index 22849fd0f71..3096bafe380 100644 --- a/opm/simulators/linalg/FlexibleSolver.hpp +++ b/opm/simulators/linalg/FlexibleSolver.hpp @@ -101,6 +101,7 @@ class FlexibleSolver : public Dune::InverseOperator preconditioner_; std::shared_ptr scalarproduct_; std::shared_ptr linsolver_; + bool HACK_is_umfpack_ = false; }; } // namespace Dune diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp index 3002c2e25b9..938825c207c 100644 --- a/opm/simulators/linalg/FlexibleSolver_impl.hpp +++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp @@ -73,6 +73,12 @@ namespace Dune FlexibleSolver:: apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) { +#if HAVE_SUITESPARSE_UMFPACK + if (HACK_is_umfpack_) { + using MatrixType = std::remove_const_tgetmat())>>; + linsolver_.reset(new Dune::UMFPack(linearoperator_for_solver_->getmat(), 0, false)); + } +#endif linsolver_->apply(x, rhs, res); } @@ -81,6 +87,12 @@ namespace Dune FlexibleSolver:: apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res) { +#if HAVE_SUITESPARSE_UMFPACK + if (HACK_is_umfpack_) { + using MatrixType = std::remove_const_tgetmat())>>; + linsolver_.reset(new Dune::UMFPack(linearoperator_for_solver_->getmat(), 0, false)); + } +#endif linsolver_->apply(x, rhs, reduction, res); } @@ -185,9 +197,9 @@ namespace Dune verbosity); #if HAVE_SUITESPARSE_UMFPACK } else if (solver_type == "umfpack") { - bool dummy = false; using MatrixType = std::remove_const_tgetmat())>>; - linsolver_ = std::make_shared>(linearoperator_for_solver_->getmat(), verbosity, dummy); + linsolver_ = std::make_shared>(linearoperator_for_solver_->getmat(), verbosity, false); + HACK_is_umfpack_ = true; #endif #if HAVE_CUDA } else if (solver_type == "cubicgstab") { From 34d36aa52fd3e15ffa6a9c11351662ded567de92 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 6 Jan 2023 11:40:25 +0100 Subject: [PATCH 05/20] Changes to well classes. --- opm/simulators/wells/BlackoilWellModel.hpp | 38 ++ .../wells/BlackoilWellModel_impl.hpp | 341 +++++++++++++++++- opm/simulators/wells/StandardWell.hpp | 19 + .../wells/StandardWellPrimaryVariables.hpp | 6 +- opm/simulators/wells/WellInterface.hpp | 10 + opm/simulators/wells/WellInterfaceGeneric.cpp | 5 + opm/simulators/wells/WellInterfaceGeneric.hpp | 3 + 7 files changed, 417 insertions(+), 5 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 0abbaad9d2a..6b731040e19 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -45,6 +45,7 @@ #include #include +#include #include #include #include @@ -137,6 +138,8 @@ namespace Opm { using AverageRegionalPressureType = RegionAverageCalculator:: AverageRegionalPressure >; + using Domain = SubDomain; + BlackoilWellModel(Simulator& ebosSimulator); void init(); @@ -155,12 +158,18 @@ namespace Opm { {} void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) override; + void linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res); void postSolve(GlobalEqVector& deltaX) override { recoverWellSolutionAndUpdateWellState(deltaX); } + void postSolveLocal(GlobalEqVector& deltaX, const Domain& domain) + { + recoverWellSolutionAndUpdateWellStateLocal(deltaX, domain); + } + ///////////// // ///////////// @@ -270,6 +279,11 @@ namespace Opm { // Check if well equations is converged. ConvergenceReport getWellConvergence(const std::vector& B_avg, const bool checkWellGroupControls = false) const; + // Check if well equations are converged locally. + ConvergenceReport getLocalWellConvergence(const Domain& domain, + const std::vector& B_avg, + const bool checkGroupConvergence = false) const; + const SimulatorReportSingle& lastReport() const; void addWellContributions(SparseMatrixAdapter& jacobian) const; @@ -292,6 +306,7 @@ namespace Opm { // at the beginning of each time step (Not report step) void prepareTimeStep(DeferredLogger& deferred_logger); void initPrimaryVariablesEvaluation() const; + void initPrimaryVariablesEvaluationLocal(const Domain& domain) const; std::pair updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false); @@ -328,6 +343,20 @@ namespace Opm { int numLocalNonshutWells() const; + // prototype for assemble function for ASPIN solveLocal() + // will try to merge back to assemble() when done prototyping + void assembleLocal(const int iterationIdx, + const double dt, + const Domain& domain); + void updateWellControlsLocal(DeferredLogger& deferred_logger, const bool checkGroupControls, + const Domain& domain); + + void logPrimaryVars() const; + std::vector getPrimaryVarsDomain(const Domain& domain) const; + void setPrimaryVarsDomain(const Domain& domain, const std::vector& vars); + + void setupDomains(const std::vector& domains); + protected: Simulator& ebosSimulator_; @@ -377,6 +406,9 @@ namespace Opm { std::vector B_avg_{}; + // Keep track of the domain of each well, if using subdomains. + std::map well_domain_; + const Grid& grid() const { return ebosSimulator_.vanguard().grid(); } @@ -404,6 +436,7 @@ namespace Opm { const double dt, DeferredLogger& local_deferredLogger); + // called at the end of a time step void timeStepSucceeded(const double& simulationTime, const double dt); @@ -414,6 +447,10 @@ namespace Opm { // xw to update Well State void recoverWellSolutionAndUpdateWellState(const BVector& x); + // using the solution x to recover the solution xw for wells and applying + // xw to update Well State + void recoverWellSolutionAndUpdateWellStateLocal(const BVector& x, const Domain& domain); + // setting the well_solutions_ based on well_state. void updatePrimaryVariables(DeferredLogger& deferred_logger); @@ -438,6 +475,7 @@ namespace Opm { int reportStepIndex() const; void assembleWellEq(const double dt, DeferredLogger& deferred_logger); + void assembleWellEqLocal(const double dt, const Domain& domain, DeferredLogger& deferred_logger); void prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index aa47b599e59..aa5012ce7c2 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -38,6 +38,7 @@ #endif #include +#include #include #include @@ -169,6 +170,36 @@ namespace Opm { } + template + void + BlackoilWellModel:: + linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res) + { + if (!param_.matrix_add_well_contributions_) { + // if the well contributions are not supposed to be included explicitly in + // the matrix, we only apply the vector part of the Schur complement here. + for (const auto& well: well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + // r = r - duneC_^T * invDuneD_ * resWell_ + well->apply(res); + } + } + return; + } + + for (const auto& well: well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->addWellContributions(jacobian); + // applying the well residual to reservoir residuals + // r = r - duneC_^T * invDuneD_ * resWell_ + well->apply(res); + } + } + } + + + + template void BlackoilWellModel:: @@ -1026,6 +1057,66 @@ namespace Opm { + template + void + BlackoilWellModel:: + assembleLocal(const int iterationIdx, + const double dt, + const Domain& domain) + { + + DeferredLogger local_deferredLogger; + // TODO: not handling gas lift for now + /* if (this->glift_debug) { + const std::string msg = fmt::format( + "assemble() : iteration {}" , iterationIdx); + gliftDebug(msg, local_deferredLogger); + } */ + last_report_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + perfTimer.start(); + + /* if ( ! wellsActive() ) { + return; + } */ + + + // TODO: this function remains to be optimized, because we do not need to handle all the wells + updatePerforationIntensiveQuantities(); + + auto exc_type = ExceptionType::NONE; + std::string exc_msg; + try { + // if (iterationIdx == 0) { + // calculateExplicitQuantities(local_deferredLogger); + // prepareTimeStep(local_deferredLogger); + // } + updateWellControlsLocal(local_deferredLogger, /* check group controls */ false, domain); + + // Set the well primary variables based on the value of well solutions + initPrimaryVariablesEvaluationLocal(domain); + + // maybeDoGasLiftOptimize(local_deferredLogger); + assembleWellEqLocal(dt, domain, local_deferredLogger); + } catch (const std::runtime_error& e) { + exc_type = ExceptionType::RUNTIME_ERROR; + exc_msg = e.what(); + } catch (const std::invalid_argument& e) { + exc_type = ExceptionType::INVALID_ARGUMENT; + exc_msg = e.what(); + } catch (const std::logic_error& e) { + exc_type = ExceptionType::LOGIC_ERROR; + exc_msg = e.what(); + } catch (const std::exception& e) { + exc_type = ExceptionType::DEFAULT; + exc_msg = e.what(); + } + // @TODO errors here must be caught higher up, as this method is not called in parallel. + logAndCheckForExceptionsAndThrow(local_deferredLogger, exc_type, "assemble() failed: " + exc_msg, terminal_output_, Parallel::Communication()); + last_report_.converged = true; + last_report_.assemble_time_well += perfTimer.stop(); + } + template bool @@ -1237,6 +1328,17 @@ namespace Opm { well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger); } } + template + void + BlackoilWellModel:: + assembleWellEqLocal(const double dt, const Domain& domain, DeferredLogger& deferred_logger) + { + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger); + } + } + } template void @@ -1466,7 +1568,6 @@ namespace Opm { BlackoilWellModel:: recoverWellSolutionAndUpdateWellState(const BVector& x) { - DeferredLogger local_deferredLogger; OPM_BEGIN_PARALLEL_TRY_CATCH(); { @@ -1474,12 +1575,34 @@ namespace Opm { const auto& summary_state = ebosSimulator_.vanguard().summaryState(); well->recoverWellSolutionAndUpdateWellState(summary_state, x, this->wellState(), local_deferredLogger); } - } OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "recoverWellSolutionAndUpdateWellState() failed: ", terminal_output_, ebosSimulator_.vanguard().grid().comm()); + } + + + template + void + BlackoilWellModel:: + recoverWellSolutionAndUpdateWellStateLocal(const BVector& x, const Domain& domain) + { + DeferredLogger local_deferredLogger; + OPM_BEGIN_PARALLEL_TRY_CATCH(); + { + const auto& summary_state = this->ebosSimulator_.vanguard().summaryState(); + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->recoverWellSolutionAndUpdateWellState(summary_state, x, + this->wellState(), + local_deferredLogger); + } + } + } + OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, + "recoverWellSolutionAndUpdateWellState() failed: ", + terminal_output_, ebosSimulator_.vanguard().grid().comm()); } @@ -1495,6 +1618,81 @@ namespace Opm { } } + template + void + BlackoilWellModel:: + initPrimaryVariablesEvaluationLocal(const Domain& domain) const + { + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->initPrimaryVariablesEvaluation(); + } + } + } + + + + + + + template + ConvergenceReport + BlackoilWellModel:: + getLocalWellConvergence(const Domain& domain, + const std::vector& B_avg, + const bool checkWellGroupControls) const + { + const auto& summary_state = ebosSimulator_.vanguard().summaryState(); + + const auto relax_tolerance = ebosSimulator_.model().newtonMethod().numIterations() + > param_.strict_outer_iter_wells_; + + Opm::DeferredLogger local_deferredLogger; + + // Get global (from all processes) convergence report. + ConvergenceReport local_report; + + for (const auto& well : well_container_) { + if ((well_domain_.at(well->name()) == domain.index) && + well->isOperableAndSolvable()) + { + local_report += well->getWellConvergence(summary_state, this->wellState(), B_avg, + local_deferredLogger, relax_tolerance); + } + } + + // This function will be called once for each domain on a process, + // and the number of domains on a process will vary, so there is + // no way to communicate here. There is also no need, as a domain + // is local to a single process in our current approach. + // Therefore there is no call to gatherDeferredLogger() or to + // gatherConvergenceReport() below. + Opm::DeferredLogger global_deferredLogger = local_deferredLogger; + ConvergenceReport report = local_report; + + // the well_group_control_changed info is already communicated + // TODO: Check what the implications are here for local solves. + if (checkWellGroupControls) { + report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed); + } + + if (terminal_output_) { + global_deferredLogger.logMessages(); + } + + + // Log debug messages for NaN or too large residuals. + if (terminal_output_) { + for (const auto& f : report.wellFailures()) { + if (f.severity() == ConvergenceReport::Severity::NotANumber) { + OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName()); + } else if (f.severity() == ConvergenceReport::Severity::TooLarge) { + OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName()); + } + } + } + return report; + } @@ -1641,6 +1839,55 @@ namespace Opm { return { changed_well_group, more_network_update }; } + template + void + BlackoilWellModel:: + updateWellControlsLocal(DeferredLogger& deferred_logger, const bool checkGroupControls, const Domain& domain) + { + // Even if there are no wells active locally, we cannot + // return as the DeferredLogger uses global communication. + // For no well active globally we simply return. + if( !wellsActive() ) return ; + +/* updateAndCommunicateGroupData(); + + updateNetworkPressures(); */ + + std::set switched_wells; + std::set switched_groups; + +/* if (checkGroupControls) { + // Check group individual constraints. + updateGroupIndividualControls(deferred_logger, switched_groups); + + // Check group's constraints from higher levels. + updateGroupHigherControls(deferred_logger, switched_groups); + + updateAndCommunicateGroupData(); + + // Check wells' group constraints and communicate. + for (const auto& well : well_container_) { + const auto mode = WellInterface::IndividualOrGroup::Group; + const bool changed = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger); + if (changed) { + switched_wells.insert(well->name()); + } + } + updateAndCommunicateGroupData(); + } */ + + // Check individual well constraints and communicate. + for (const auto& well : well_container_) { + if (switched_wells.count(well->name()) || well_domain_.at(well->name()) != domain.index) { + continue; + } + const auto mode = WellInterface::IndividualOrGroup::Individual; + well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger); + } + // updateAndCommunicateGroupData(); + } + + template void @@ -1949,6 +2196,7 @@ namespace Opm { BlackoilWellModel:: updatePerforationIntensiveQuantities() { +#if 0 ElementContext elemCtx(ebosSimulator_); const auto& gridView = ebosSimulator_.gridView(); @@ -1963,6 +2211,7 @@ namespace Opm { elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); } OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updatePerforationIntensiveQuantities() failed: ", ebosSimulator_.vanguard().grid().comm()); +#endif } @@ -2094,6 +2343,59 @@ namespace Opm { } + template + void + BlackoilWellModel:: + logPrimaryVars() const + { + std::ostringstream os; + for (const auto& w : well_container_) { + os << w->name() << ":"; + auto pv = w->getPrimaryVars(); + for (const double v : pv) { + os << ' ' << v; + } + os << '\n'; + } + OpmLog::debug(os.str()); + } + + + + template + std::vector + BlackoilWellModel:: + getPrimaryVarsDomain(const Domain& domain) const + { + std::vector ret; + for (const auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + const auto& pv = well->getPrimaryVars(); + for (const double v : pv) { + ret.push_back(v); + } + } + } + return ret; + } + + + + template + void + BlackoilWellModel:: + setPrimaryVarsDomain(const Domain& domain, const std::vector& vars) + { + int offset = 0; + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + int num_pri_vars = well->setPrimaryVars(vars.begin() + offset); + offset += num_pri_vars; + } + } + } + + template void @@ -2118,7 +2420,38 @@ namespace Opm { } - - + template + void + BlackoilWellModel:: + setupDomains(const std::vector& domains) + { + for (const auto& wellPtr : this->well_container_) { + const int first_well_cell = wellPtr->cells()[0]; + for (const auto& domain : domains) { + const bool found = std::binary_search(domain.cells.begin(), domain.cells.end(), first_well_cell); + if (found) { + // Assuming that if the first well cell is found in a domain, + // then all of that well's cells are in that same domain. + // TODO verify assumption. + well_domain_[wellPtr->name()] = domain.index; + + // Verify that all of that well's cells are in that same domain. + for (int well_cell : wellPtr->cells()) { + if (!std::binary_search(domain.cells.begin(), domain.cells.end(), well_cell)) { + // TODO: ensure proper exit in parallel. + OPM_THROW(std::runtime_error, "Well found on multiple domains."); + } + } + } + } + } + // TODO: ensure output on rank 0 only. + std::ostringstream os; + os << "Well name Domain\n"; + for (const auto& [wname, domain] : well_domain_) { + os << wname << std::setw(21 - wname.size()) << domain << '\n'; + } + OpmLog::debug(os.str()); + } } // namespace Opm diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index d2febff4af1..4c48699a497 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -254,6 +254,25 @@ namespace Opm DeferredLogger& deferred_logger) const; + std::vector getPrimaryVars() const override + { + const int num_pri_vars = this->primary_variables_.numWellEq(); + std::vector retval(num_pri_vars); + for (int ii = 0; ii < num_pri_vars; ++ii) { + retval[ii] = this->primary_variables_.value(ii); + } + return retval; + } + + int setPrimaryVars(std::vector::const_iterator it) override + { + const int num_pri_vars = this->primary_variables_.numWellEq(); + for (int ii = 0; ii < num_pri_vars; ++ii) { + this->primary_variables_.setValue(ii, it[ii]); + } + return num_pri_vars; + } + protected: bool regularize_; diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index 31951536ede..10037f6771a 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -135,7 +135,7 @@ class StandardWellPrimaryVariables { //! \brief Returns scaled rate for a component. EvalWell getQs(const int compIdx) const; - //! \brief Returns a const ref to an evaluation. + //! \brief Returns a value. Scalar value(const int idx) const { return value_[idx]; } @@ -143,6 +143,10 @@ class StandardWellPrimaryVariables { const EvalWell& eval(const int idx) const { return evaluation_[idx]; } + //! \brief Set a value. Note that this does not also set the corresponding evaluation. + void setValue(const int idx, const Scalar val) + { value_[idx] = val; } + private: //! \brief Calculate a relaxation factor for producers. //! \details To avoid overshoot of the fractions which might result in negative rates. diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index d3a2fa336ad..cbb49126fed 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -328,6 +328,16 @@ class WellInterface : public WellInterfaceIndices getPrimaryVars() const + { + return {}; + } + + virtual int setPrimaryVars(std::vector::const_iterator) + { + return 0; + } + protected: // simulation parameters const ModelParameters& param_; diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 754372b8c04..88c6e5fc135 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -397,6 +397,11 @@ bool WellInterfaceGeneric::isVFPActive(DeferredLogger& deferred_logger) const } } +const Well& WellInterfaceGeneric::scheduleWell() const +{ + return well_ecl_; +} + bool WellInterfaceGeneric::isOperableAndSolvable() const { return operability_status_.isOperableAndSolvable(); diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 674e50fff27..852b0956de8 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -197,6 +197,9 @@ class WellInterfaceGeneric { bool stopppedOrZeroRateTarget(const SummaryState& summary_state, const WellState& well_state) const; + + const Well& scheduleWell() const; + protected: bool getAllowCrossFlow() const; From 3083eecba206bfb613206a4adeef112e9e9db36a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 9 Jan 2023 10:46:19 +0100 Subject: [PATCH 06/20] Cleanup file and add LinearSolverPrintJsonDefinition. 1. The order of options has been revised slightly to be more logical, and the order is now the same in all 6 places they appear: Declaration of parameters, definition of defaults, struct definition, init(), registerParameters() and reset(). 2. The reset() function now sets the members to the corresponding parameter defaults. 3. Removed some unused members and parameters. --- .../linalg/FlowLinearSolverParameters.hpp | 90 ++++++++++--------- 1 file changed, 50 insertions(+), 40 deletions(-) diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 425c41a2c3e..23803d16352 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -51,10 +51,6 @@ struct RelaxedLinearSolverReduction { using type = UndefinedProperty; }; -template -struct IluRelaxation { - using type = UndefinedProperty; -}; template struct LinearSolverMaxIter { using type = UndefinedProperty; @@ -67,6 +63,10 @@ struct LinearSolverRestart { // LinearSolverVerbosity defined in opm-models // template +struct IluRelaxation { + using type = UndefinedProperty; +}; +template struct IluFillinLevel { using type = UndefinedProperty; }; @@ -91,23 +91,23 @@ struct LinearSolverIgnoreConvergenceFailure{ using type = UndefinedProperty; }; template -struct PreconditionerAddWellContributions { +struct ScaleLinearSystem { using type = UndefinedProperty; }; template -struct ScaleLinearSystem { +struct LinearSolver { using type = UndefinedProperty; }; template -struct CprReuseSetup { +struct LinearSolverPrintJsonDefinition { using type = UndefinedProperty; }; template -struct CprReuseInterval { +struct CprReuseSetup { using type = UndefinedProperty; }; template -struct LinearSolver { +struct CprReuseInterval { using type = UndefinedProperty; }; template @@ -137,11 +137,6 @@ struct RelaxedLinearSolverReduction { static constexpr type value = 1e-2; }; template -struct IluRelaxation { - using type = GetPropType; - static constexpr type value = 0.9; -}; -template struct LinearSolverMaxIter { static constexpr int value = 200; }; @@ -154,6 +149,11 @@ struct LinearSolverVerbosity { static constexpr int value = 0; }; template +struct IluRelaxation { + using type = GetPropType; + static constexpr type value = 0.9; +}; +template struct IluFillinLevel { static constexpr int value = 0; }; @@ -178,16 +178,16 @@ struct LinearSolverIgnoreConvergenceFailure static constexpr bool value = false; }; template -struct LinearSolverBackend { - using type = ISTLSolverEbos; +struct ScaleLinearSystem { + static constexpr bool value = false; }; template -struct PreconditionerAddWellContributions { - static constexpr bool value = false; +struct LinearSolver { + static constexpr auto value = "ilu0"; }; template -struct ScaleLinearSystem { - static constexpr bool value = false; +struct LinearSolverPrintJsonDefinition { + static constexpr auto value = true; }; template struct CprReuseSetup { @@ -198,10 +198,6 @@ struct CprReuseInterval { static constexpr int value = 30; }; template -struct LinearSolver { - static constexpr auto value = "ilu0"; -}; -template struct AcceleratorMode { static constexpr auto value = "none"; }; @@ -218,6 +214,12 @@ struct OpenclIluParallel { static constexpr bool value = true; // note: false should only be used in debug }; +// Set the backend to be used. +template +struct LinearSolverBackend { + using type = ISTLSolverEbos; +}; + } // namespace Opm::Properties namespace Opm @@ -229,24 +231,24 @@ namespace Opm { double linear_solver_reduction_; double relaxed_linear_solver_reduction_; - double ilu_relaxation_; int linear_solver_maxiter_; int linear_solver_restart_; int linear_solver_verbosity_; + double ilu_relaxation_; int ilu_fillin_level_; MILU_VARIANT ilu_milu_; bool ilu_redblack_; bool ilu_reorder_sphere_; bool newton_use_gmres_; - bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; bool scale_linear_system_; std::string linsolver_; + bool linear_solver_print_json_definition_; + int cpr_reuse_setup_; + int cpr_reuse_interval_; std::string accelerator_mode_; int bda_device_id_; int opencl_platform_id_; - int cpr_reuse_setup_; - int cpr_reuse_interval_; bool opencl_ilu_parallel_; template @@ -255,10 +257,10 @@ namespace Opm // TODO: these parameters have undocumented non-trivial dependencies linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction); relaxed_linear_solver_reduction_ = EWOMS_GET_PARAM(TypeTag, double, RelaxedLinearSolverReduction); - ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation); linear_solver_maxiter_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter); linear_solver_restart_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverRestart); linear_solver_verbosity_ = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity); + ilu_relaxation_ = EWOMS_GET_PARAM(TypeTag, double, IluRelaxation); ilu_fillin_level_ = EWOMS_GET_PARAM(TypeTag, int, IluFillinLevel); ilu_milu_ = convertString2Milu(EWOMS_GET_PARAM(TypeTag, std::string, MiluVariant)); ilu_redblack_ = EWOMS_GET_PARAM(TypeTag, bool, IluRedblack); @@ -266,6 +268,8 @@ namespace Opm newton_use_gmres_ = EWOMS_GET_PARAM(TypeTag, bool, UseGmres); ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure); scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem); + linsolver_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver); + linear_solver_print_json_definition_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition); cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); cpr_reuse_interval_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseInterval); @@ -284,12 +288,12 @@ namespace Opm template static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve should try to achive"); + EWOMS_REGISTER_PARAM(TypeTag, double, LinearSolverReduction, "The minimum reduction of the residual which the linear solver must achieve"); EWOMS_REGISTER_PARAM(TypeTag, double, RelaxedLinearSolverReduction, "The minimum reduction of the residual which the linear solver need to achieve for the solution to be accepted"); - EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIter, "The maximum number of iterations of the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverRestart, "The number of iterations after which GMRES is restarted"); EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, "The verbosity level of the linear solver (0: off, 2: all)"); + EWOMS_REGISTER_PARAM(TypeTag, double, IluRelaxation, "The relaxation factor of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, int, IluFillinLevel, "The fill-in level of the linear solver's ILU preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, std::string, MiluVariant, "Specify which variant of the modified-ILU preconditioner ought to be used. Possible variants are: ILU (default, plain ILU), MILU_1 (lump diagonal with dropped row entries), MILU_2 (lump diagonal with the sum of the absolute values of the dropped row entries), MILU_3 (if diagonal is positive add sum of dropped row entrires. Otherwise subtract them), MILU_4 (if diagonal is positive add sum of dropped row entrires. Otherwise do nothing"); EWOMS_REGISTER_PARAM(TypeTag, bool, IluRedblack, "Use red-black partitioning for the ILU preconditioner"); @@ -297,6 +301,8 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge"); EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cpr (an alias for cpr_trueimpes), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); + EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition, "Write the JSON definition of the linear solver setup to the DBG wfle."); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate, 4: recreated every CprReuseInterval"); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseInterval, "Reuse preconditioner interval. Used when CprReuseSetup is set to 4, then the preconditioner will be fully recreated instead of reused every N linear solve, where N is this parameter."); EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); @@ -311,19 +317,23 @@ namespace Opm // set default values void reset() { - newton_use_gmres_ = false; - linear_solver_reduction_ = 1e-2; relaxed_linear_solver_reduction_ = 1e-2; - linear_solver_maxiter_ = 150; - linear_solver_restart_ = 40; - linear_solver_verbosity_ = 0; - require_full_sparsity_pattern_ = false; - ignoreConvergenceFailure_ = false; - ilu_fillin_level_ = 0; + linear_solver_reduction_ = 1e-2; + linear_solver_maxiter_ = 200; + linear_solver_restart_ = 40; + linear_solver_verbosity_ = 0; ilu_relaxation_ = 0.9; + ilu_fillin_level_ = 0; ilu_milu_ = MILU_VARIANT::ILU; ilu_redblack_ = false; - ilu_reorder_sphere_ = true; + ilu_reorder_sphere_ = false; + newton_use_gmres_ = false; + ignoreConvergenceFailure_ = false; + scale_linear_system_ = false; + linsolver_ = "ilu0"; + linear_solver_print_json_definition_ = true; + cpr_reuse_setup_ = 4; + cpr_reuse_interval_ = 30; accelerator_mode_ = "none"; bda_device_id_ = 0; opencl_platform_id_ = 0; From 9b6d6ad2ca6c07f954e3860133b78d3bf53690ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 9 Jan 2023 11:03:09 +0100 Subject: [PATCH 07/20] Add two new functions to ISTLSolverEbos. 1. A new constructor taking the FlowLinearSolverParameters object directly. 2. A new overload of prepare() taking the matrix directly rather than the SparseMatrixAdapter class. --- opm/simulators/linalg/ISTLSolverEbos.hpp | 47 ++++++++++++++++++------ 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 068a3ce93e4..c7c2656db63 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -231,8 +231,22 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } /// Construct a system solver. - /// \param[in] parallelInformation In the case of a parallel run - /// with dune-istl the information about the parallelization. + /// \param[in] simulator The opm-models simulator object + /// \param[in] parameters Explicit parameters for solver setup, do not + /// read them from command line parameters. + ISTLSolverEbos(const Simulator& simulator, const FlowLinearSolverParameters& parameters) + : simulator_(simulator), + iterations_( 0 ), + calls_( 0 ), + converged_(false), + matrix_(), + parameters_(parameters) + { + initialize(); + } + + /// Construct a system solver. + /// \param[in] simulator The opm-models simulator object explicit ISTLSolverEbos(const Simulator& simulator) : simulator_(simulator), iterations_( 0 ), @@ -240,16 +254,22 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, converged_(false), matrix_() { - OPM_TIMEBLOCK(IstlSolverEbos); - const bool on_io_rank = (simulator.gridView().comm().rank() == 0); -#if HAVE_MPI - comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); -#endif parameters_.template init(simulator_.vanguard().eclState().getSimulationConfig().useCPR()); + initialize(); + } + + void initialize() + { + OPM_TIMEBLOCK(IstlSolverEbos); prm_ = setupPropertyTree(parameters_, EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter), EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction)); + const bool on_io_rank = (simulator_.gridView().comm().rank() == 0); +#if HAVE_MPI + comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); +#endif + #if COMPILE_BDA_BRIDGE { std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode); @@ -300,7 +320,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true); // Print parameters to PRT/DBG logs. - if (on_io_rank) { + if (on_io_rank && parameters_.linear_solver_print_json_definition_) { std::ostringstream os; os << "Property tree for linear solver:\n"; prm_.write_json(os, true); @@ -314,12 +334,17 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } void prepare(const SparseMatrixAdapter& M, Vector& b) + { + prepare(M.istlMatrix(), b); + } + + void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverEbosPrepare); static bool firstcall = true; #if HAVE_MPI if (firstcall) { - const size_t size = M.istlMatrix().N(); + const size_t size = M.N(); detail::copyParValues(parallelInformation_, size, *comm_); } #endif @@ -329,7 +354,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, // ebos will not change the matrix object. Hence simply store a pointer // to the original one with a deleter that does nothing. // Outch! We need to be able to scale the linear system! Hence const_cast - matrix_ = const_cast(&M.istlMatrix()); + matrix_ = const_cast(&M); useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver) @@ -343,7 +368,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #endif } else { // Pointers should not change - if ( &(M.istlMatrix()) != matrix_ ) { + if ( &M != matrix_ ) { OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!"); } From 90cb66a8427c35351fcaea2a6ff5cf1f934ebd70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 6 Jan 2023 11:41:22 +0100 Subject: [PATCH 08/20] Main implementation parts for aspin. --- opm/simulators/flow/BlackoilModelEbos.hpp | 1337 ++++++++++++++++++++- opm/simulators/flow/FlowMainEbos.hpp | 7 + 2 files changed, 1314 insertions(+), 30 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index fcd6bd1edf2..02e07e8972b 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -27,6 +27,24 @@ #include #include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + #include #include #include @@ -43,6 +61,9 @@ #include #include #include + + + #include #include #include @@ -52,11 +73,17 @@ #include #include + +#include +#include + #include #include #include #include +#include + #include #include #include @@ -200,6 +227,8 @@ namespace Opm { typedef typename SparseMatrixAdapter::IstlMatrix Mat; typedef Dune::BlockVector BVector; + using Domain = SubDomain; + typedef ISTLSolverEbos ISTLSolverType; class ComponentName @@ -295,14 +324,96 @@ namespace Opm { // compute global sum of number of cells global_nc_ = detail::countGlobalCells(grid_); convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves. + // TODO: remember to fix! + if (param_.nonlinear_solver_ == "aspin" || param_.nonlinear_solver_ == "nldd" || param_.nonlinear_solver_ == "purelocal") { + setupAspinDomains(); + } } + bool isParallel() const { return grid_.comm().size() > 1; } + const EclipseState& eclState() const { return ebosSimulator_.vanguard().eclState(); } + + + void setupAspinDomains() + { + // Create partitions. + const int num_cells = detail::countLocalInteriorCells(grid_); + auto [partition_vector, num_partitions] = partitionCells(num_cells); + + // Scan through partitioning to get correct size for each. + std::vector sizes(num_partitions, 0); + for (int p : partition_vector) { + ++sizes[p]; + } + + // Set up correctly sized vectors of entity seeds and of indices for each partition. + using EntitySeed = typename Grid::template Codim<0>::EntitySeed; + std::vector> seeds(num_partitions); + std::vector> partitions(num_partitions); + for (int part = 0; part < num_partitions; ++part) { + seeds[part].resize(sizes[part]); + partitions[part].resize(sizes[part]); + } + + // Iterate through grid once, setting the seeds of all partitions. + std::vector count(num_partitions, 0); + const auto beg = grid_.template leafbegin<0>(); + const auto end = grid_.template leafend<0>(); + int cell = 0; + for (auto it = beg; it != end; ++it, ++cell) { + const int p = partition_vector[cell]; + seeds[p][count[p]] = it->seed(); + partitions[p][count[p]] = cell; + ++count[p]; + } + assert(cell == num_cells); + assert(count == sizes); + + // Create the domains. + const int num_domains = partitions.size(); + for (int index = 0; index < num_domains; ++index) { + std::vector interior(num_cells, false); + for (int ix : partitions[index]) { + interior[ix] = true; + } + Dune::SubGridView view(ebosSimulator_.vanguard().grid(), std::move(seeds[index])); + domains_.emplace_back(index, std::move(partitions[index]), std::move(interior), std::move(view)); + } + + // Set up container for the local system matrices. + domain_matrices_.resize(num_domains); + + // Set up container for the local linear solvers. + for (int index = 0; index < num_domains; ++index) { + // TODO: The ISTLSolverEbos constructor will make + // parallel structures appropriate for the full grid + // only. This must be addressed before going parallel. + FlowLinearSolverParameters param; + param.template init(); + // Override solver type with umfpack if small domain. + // Otherwise hardcode to ILU0 + if (domains_[index].cells.size() < 200) { + param.linsolver_ = "umfpack.json"; + } else { + param.linsolver_ = "ilu0"; + param.linear_solver_reduction_ = 1e-2; + } + param.linear_solver_print_json_definition_ = false; + domain_linsolvers_.emplace_back(ebosSimulator_, param); + } + + assert(int(domains_.size()) == num_domains); + } + + + + /// Called once before each time step. /// \param[in] timer simulation timer SimulatorReportSingle prepareStep(const SimulatorTimerInterface& timer) @@ -335,6 +446,10 @@ namespace Opm { std::cout << "equation scaling not supported yet" << std::endl; //updateEquationsScaling(); } + + // Setup domain->well mapping. + wellModel().setupDomains(domains_); + report.pre_post_time += perfTimer.stop(); return report; @@ -348,18 +463,11 @@ namespace Opm { /// \param[in] iteration should be 0 for the first call of a new timestep /// \param[in] timer simulation timer /// \param[in] nonlinear_solver nonlinear solver used (for oscillation/relaxation control) - /// \param[in, out] reservoir_state reservoir state variables - /// \param[in, out] well_state well state variables template SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface& timer, NonlinearSolverType& nonlinear_solver) { - SimulatorReportSingle report; - failureReport_ = SimulatorReportSingle(); - Dune::Timer perfTimer; - - perfTimer.start(); if (iteration == 0) { // For each iteration we store in a vector the norms of the residual of // the mass balance for each active phase, the well flux and the well equations. @@ -370,8 +478,29 @@ namespace Opm { convergence_reports_.back().report.reserve(11); } + if (param_.nonlinear_solver_ == "aspin" || param_.nonlinear_solver_ == "nldd" || param_.nonlinear_solver_ == "purelocal") { + return nonlinearIterationAspin(iteration, timer, nonlinear_solver); + } else { + return nonlinearIterationNewton(iteration, timer, nonlinear_solver); + } + } + + + template + SimulatorReportSingle nonlinearIterationNewton(const int iteration, + const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + + // ----------- Set up reports and timer ----------- + SimulatorReportSingle report; + failureReport_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + + perfTimer.start(); report.total_linearizations = 1; + // ----------- Assemble ----------- try { report += assembleReservoir(timer, iteration); report.assemble_time += perfTimer.stop(); @@ -383,12 +512,41 @@ namespace Opm { throw; // continue throwing the stick } +#if 0 + // DEBUG HACK + // Check if domain-diagonal parts are identical to stored subdomain matrices. + if (param_.nonlinear_solver_ == "nldd") { + const Mat& main_matrix = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); + for (const auto& domain : domains_) { + if (domain_matrices_[domain.index]) { + auto domain_diag_part = Details::extractMatrix(main_matrix, domain.cells); + if (Details::matrixEqual(domain_diag_part, *domain_matrices_[domain.index])) { + OpmLog::debug(" *** matrices IDENTICAL for domain " + std::to_string(domain.index)); + } else { + OpmLog::debug(" *** matrices DIFFER for domain " + std::to_string(domain.index)); + } + } else { + OpmLog::debug(" *** no local matrix for domain " + std::to_string(domain.index)); + } + } + } +#endif + + // ----------- Check if converged ----------- std::vector residual_norms; perfTimer.reset(); perfTimer.start(); // the step is not considered converged until at least minIter iterations is done { - auto convrep = getConvergence(timer, iteration,residual_norms); + auto convrep = getConvergence(timer, iteration, residual_norms); + // if (!convrep.converged()) { + // std::ostringstream os; + // os << "Convergence data for Newton iteration " << iteration << "\n" + // << convrep; + // OpmLog::debug(os.str()); + // } + ////wellModel().logPrimaryVars(); + report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();; ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); convergence_reports_.back().report.push_back(std::move(convrep)); @@ -400,8 +558,11 @@ namespace Opm { OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); } } - report.update_time += perfTimer.stop(); + //report.convergence_check_time += perfTimer.stop(); residual_norms_history_.push_back(residual_norms); + // Dune::writeMatrixToMatlab(ebosSimulator_.model().linearizer().jacobian().istlMatrix(), "beforeWellModelLinearize.matlab"); + + // ----------- If not converged, solve linear system ----------- if (!report.converged) { perfTimer.reset(); perfTimer.start(); @@ -422,6 +583,7 @@ namespace Opm { // Note that linearize may throw for MSwells. wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), ebosSimulator().model().linearizer().residual()); + //wellModel().logPrimaryVars(); solveJacobianSystem(x); report.linear_solve_setup_time += linear_solve_setup_time_; @@ -444,12 +606,13 @@ namespace Opm { // on observation to avoid some big performance degeneration under some circumstances. // there is no theorectical explanation which way is better for sure. wellModel().postSolve(x); + //wellModel().logPrimaryVars(); if (param_.use_update_stabilization_) { // Stabilize the nonlinear update. bool isOscillate = false; bool isStagnate = false; - nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate); + nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate); if (isOscillate) { current_relaxation_ -= nonlinear_solver.relaxIncrement(); current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); @@ -472,6 +635,731 @@ namespace Opm { return report; } + + + + template + SimulatorReportSingle nonlinearIterationAspin(const int iteration, + const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + // ----------- Set up reports and timer ----------- + SimulatorReportSingle report; + failureReport_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + + //wellModel().logPrimaryVars(); + + perfTimer.start(); + report.total_linearizations = 1; + + // ----------- Assemble ----------- + try { + report += assembleReservoir(timer, iteration); + report.assemble_time += perfTimer.stop(); + } + catch (...) { + report.assemble_time += perfTimer.stop(); + failureReport_ += report; + // todo (?): make the report an attribute of the class + throw; // continue throwing the stick + } + //wellModel().logPrimaryVars(); + + // ----------- Check if converged ----------- + std::vector residual_norms; + perfTimer.reset(); + perfTimer.start(); + // the step is not considered converged until at least minIter iterations is done + { + auto convrep = getConvergence(timer, iteration, residual_norms); + if (!convrep.converged()) { + // std::ostringstream os; + // os << "Convergence data for first global iteration:\n" + // << convrep; + // OpmLog::debug(os.str()); + } + report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();; + ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); + convergence_reports_.back().report.push_back(std::move(convrep)); + + // Throw if any NaN or too large residual found. + if (severity == ConvergenceReport::Severity::NotANumber) { + OPM_THROW(NumericalProblem, "NaN residual found!"); + } else if (severity == ConvergenceReport::Severity::TooLarge) { + OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); + } + } + //report.convergence_check_time += perfTimer.stop(); + residual_norms_history_.push_back(residual_norms); + + if (report.converged) { + return report; + } + + // ----------- If not converged, do an ASPIN or NLDD iteration ----------- + + // apply the Schur compliment of the well model to the reservoir linearized + // equations + // wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), + // ebosSimulator().model().linearizer().residual()); + + // Take a copy of the FI residual. + // auto fully_implicit_residual = ebosSimulator().model().linearizer().residual(); + // ... and the Jacobian. + // auto fully_implicit_jacobian = ebosSimulator().model().linearizer().jacobian().istlMatrix(); + // ... and the solution state that generated it. + auto& solution = ebosSimulator().model().solution(0); + auto initial_solution = solution; + // auto initial_wellstate = wellModel().wellState(); + auto locally_solved = initial_solution; + + // ----------- Decide on an ordering for the domains ----------- + std::vector domain_order(domains_.size()); + if (param_.local_solve_approach_ == "gauss-seidel") { + if (true) { + // Use average pressures to order domains. + std::vector> avgpress_per_domain(domains_.size()); + for (const auto& domain : domains_) { + double press_sum = 0.0; + for (const int c : domain.cells) { + press_sum += solution[c][Indices::pressureSwitchIdx]; + } + const double avgpress = press_sum / domain.cells.size(); + avgpress_per_domain[domain.index] = std::make_pair(avgpress, domain.index); + } + // Lexicographical sort by pressure, then index. + std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end()); + // Reverse since we want high-pressure regions solved first. + std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end()); + for (size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = avgpress_per_domain[ii].second; + } + } else { + // Use maximum residual to order domains. + const auto& residual = ebosSimulator_.model().linearizer().residual(); + const int num_vars = residual[0].size(); + std::vector> maxres_per_domain(domains_.size()); + for (const auto& domain : domains_) { + double maxres = 0.0; + for (const int c : domain.cells) { + for (int ii = 0; ii < num_vars; ++ii) { + maxres = std::max(maxres, std::fabs(residual[c][ii])); + } + } + maxres_per_domain[domain.index] = std::make_pair(maxres, domain.index); + } + // Lexicographical sort by pressure, then index. + std::sort(maxres_per_domain.begin(), maxres_per_domain.end()); + // Reverse since we want high-pressure regions solved first. + std::reverse(maxres_per_domain.begin(), maxres_per_domain.end()); + for (size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = maxres_per_domain[ii].second; + } + } + } else { + std::iota(domain_order.begin(), domain_order.end(), 0); + } + + // ----------- Solve each domain separately ----------- + std::vector domain_reports(domains_.size()); + for (const int domain_index : domain_order) { + const auto& domain = domains_[domain_index]; + SimulatorReportSingle local_report; + if (param_.local_solve_approach_ == "jacobi") { + auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; +#if 0 + auto local_solution = Details::extractVector(solution, domain.cells); + Details::setGlobal(local_solution, domain.cells, locally_solved); + Details::setGlobal(initial_local_solution, domain.cells, solution); + // ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); +#else + if (/*!local_report.converged*/ false) { + // Try again with a less strict tolerance. + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + // auto param_orig = param_; + // param_.local_tolerance_scaling_cnv_ *= 10.0; + // param_.local_tolerance_scaling_mb_ *= 10.0; + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; + // param_ = param_orig; + } + if (local_report.converged) { + auto local_solution = Details::extractVector(solution, domain.cells); + Details::setGlobal(local_solution, domain.cells, locally_solved); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } else { + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } +#endif + } else { + assert(param_.local_solve_approach_ == "gauss-seidel"); + auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; + if (/*!local_report.converged*/ false) { + // Try again with a less strict tolerance. + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + // auto param_orig = param_; + // param_.local_tolerance_scaling_cnv_ = 50.0; + // param_.local_tolerance_scaling_mb_ = 50.0; + auto res = solveLocal(domain, timer, iteration); + local_report = res.first; + // param_ = param_orig; + } + if (!local_report.converged) { + // We look at the detailed convergence report to evaluate + // if we should accept the unconverged solution. + const auto& convrep = res.second; + // We do not accept a solution if the wells are unconverged. + if (!convrep.wellFailed()) { + // Calculare the sums of the mb and cnv failures. + double mb_sum = 0.0; + double cnv_sum = 0.0; + for (const auto& rf : convrep.reservoirFailures()) { + if (rf.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) { + // mb_sum += rf.magnitude(); + } else if (rf.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) { + // cnv_sum += rf.magnitude(); + } + } + // If not too high, we overrule the convergence failure. + if (mb_sum < 1e-3 && cnv_sum < 1.0) { + local_report.converged = true; + OpmLog::debug("Accepting solution in unconverged domain " + std::to_string(domain.index)); + } + } + } + if (local_report.converged) { + auto local_solution = Details::extractVector(solution, domain.cells); + Details::setGlobal(local_solution, domain.cells, locally_solved); + } else { + wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } + } + // This should have updated the global matrix to be + // dR_i/du_j evaluated at new local solutions for + // i == j, at old solution for i != j. + if (!local_report.converged) { + // TODO: more proper treatment, including in parallel. + OpmLog::debug("Convergence failure in domain " + std::to_string(domain.index)); + } + domain_reports[domain.index] = local_report; + } + +#if 0 + // DEBUG write solution vectors. + { + const int nc = UgGridHelpers::numCells(grid_); + BVector sol(nc); + for (int c = 0; c < nc; ++c) { + sol[c] = solution[c]; + } + BVector loc_sol(nc); + for (int c = 0; c < nc; ++c) { + loc_sol[c] = locally_solved[c]; + } + Dune::storeMatrixMarket(sol, param_.local_solve_approach_ + "-" + std::to_string(iteration) + ".mm"); + Dune::storeMatrixMarket(loc_sol, param_.local_solve_approach_ + "-local-" + std::to_string(iteration) + ".mm"); + } +#endif + + // Print summary of local solve convergence. + { + int num_converged = 0; + SimulatorReportSingle rep; + for (const auto& dr : domain_reports) { + if (dr.converged) { + ++num_converged; + } + rep += dr; + } + std::ostringstream os; + os << fmt::format("Local solves finished. Converged for {}/{} domains.\n", + num_converged, domain_reports.size()); + rep.reportFullyImplicit(os, nullptr); + OpmLog::debug(os.str()); + local_reports_accumulated_ += rep; + } + +#define ASPIN_EXTRA_WELL_OUTPUT 1 +#if ASPIN_EXTRA_WELL_OUTPUT + // Extra debug output. + //wellModel().logPrimaryVars(); + { + std::ostringstream os; + os << " ** Well rates:"; + for (int w = 0; w < wellModel().wellState().numWells(); ++w) { + os << " | "; + const auto& wr = wellModel().wellState().wellRates(w); + for (const double r : wr) { + os << ' ' << r; + } + } + OpmLog::debug(os.str()); + } +#endif + + if (param_.nonlinear_solver_ == "purelocal") { + if (param_.local_solve_approach_ == "jacobi") { + solution = locally_solved; + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + return report; + } + else if (param_.nonlinear_solver_ == "nldd") { + if (param_.local_solve_approach_ == "jacobi") { + solution = locally_solved; + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + // Finish with a Newton step. + // ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + auto rep = nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver); + report += rep; + if (rep.converged) { + report.converged = true; + } + return report; + } + report.total_newton_iterations = 1; + + /* + // HACK to check FI convergence + // Take a copy of the FI residual. + auto res1 = ebosSimulator().model().linearizer().residual(); + // ... and the Jacobian. + auto jac1 = ebosSimulator().model().linearizer().jacobian().istlMatrix(); + // ... and the solution state that generated it. + auto sol1 = solution; + auto ws1 = wellModel().wellState(); + + // ----------- Assemble ----------- + try { + report += assembleReservoir(timer, iteration); + report.assemble_time += perfTimer.stop(); + } + catch (...) { + report.assemble_time += perfTimer.stop(); + failureReport_ += report; + // todo (?): make the report an attribute of the class + throw; // continue throwing the stick + } + //wellModel().logPrimaryVars(); + + // ----------- Check if converged ----------- + perfTimer.reset(); + perfTimer.start(); + // the step is not considered converged until at least minIter iterations is done + { + auto convrep = getConvergence(timer, iteration, residual_norms); + report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();; + ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); + convergence_reports_.back().report.push_back(std::move(convrep)); + + // Throw if any NaN or too large residual found. + if (severity == ConvergenceReport::Severity::NotANumber) { + OPM_THROW(NumericalProblem, "NaN residual found!"); + } else if (severity == ConvergenceReport::Severity::TooLarge) { + OPM_THROW(NumericalProblem, "Too large residual found!"); + } + } + report.update_time += perfTimer.stop(); + residual_norms_history_.push_back(residual_norms); + + if (report.converged) { + return report; + } + // HACK reinstate everything + ebosSimulator().model().linearizer().residual() = res1; + ebosSimulator().model().linearizer().jacobian().istlMatrix() = jac1; + solution = sol1; + wellModel().wellState() = ws1; + */ + + // ----------- Compute ASPIN residual, check convergence ----------- + const int nc = grid_.size(0); // Includes both owned cells and overlap. + BVector aspin_residual(nc); + const int num_vars = aspin_residual[0].size(); + std::vector max_sol(num_vars, 0.0); + std::vector max_diff(num_vars, 0.0); + for (int c = 0; c < nc; ++c) { + aspin_residual[c] = initial_solution[c] - locally_solved[c]; + for (int var = 0; var < num_vars; ++var) { + max_sol[var] = std::max(max_sol[var], std::fabs(locally_solved[c][var])); + max_diff[var] = std::max(max_diff[var], std::fabs(aspin_residual[c][var])); + } + } + double max_scaled_diff = 0.0; + for (int var = 0; var < num_vars; ++var) { + max_scaled_diff = std::max(max_scaled_diff, max_diff[var]/max_sol[var]); + } + OpmLog::debug(" Scaled ASPIN residual: " + std::to_string(max_scaled_diff)); + if (max_scaled_diff < param_.outer_aspin_tolerance_) { + report.converged = true; + solution = locally_solved; + return report; + } + + // ----------- Solve global linear system ----------- + // The matrix is kept as is for now. + // TODO: test with proper ASPIN (in what we have now, columns are updated outside the diagonal + // domain-domain interaction blocks in a Gauss-Seidel-like manner), and other variations. + // NOTE: the above should no longer be true, the matrix is ASPIN now (should be). + // +#if 0 + // Compute rhs. + // The modified residual is D * aspin_residual in each domain. + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + for (const auto& domain : domains_) { + auto aspin_res_local = Details::extractVector(aspin_residual, domain.cells); + BVector aspin_res_mod(aspin_res_local.size()); + aspin_res_mod = 0; + if (domain_reports[domain.index].total_newton_iterations > 0) { + (*domain_matrices_[domain.index]).mv(aspin_res_local, aspin_res_mod); + } + Details::setGlobal(aspin_res_mod, domain.cells, ebosResid); + } + // Add J_aspin * aspin_residual to the modified residual, so we solve for + // the difference from the local solution rather than from the very start + // of the iteration. + { + const auto& ebosJac = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); + BVector yy(nc); + ebosJac.mv(aspin_residual, yy); + ebosResid -= yy; + } +#else + // Compute rhs. + { + auto Jcopy = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); + Dune::storeMatrixMarket(Jcopy, "aspinjac.mm"); + // Eliminate the domain diagonal blocks from Jcopy. Jcopy is then (J - D) + for (const auto& domain : domains_) { + for (const int cell : domain.cells) { + auto& row = Jcopy[cell]; + for (auto iter = row.begin(); iter != row.end(); ++iter) { + int col = iter.index(); + auto lb = std::find(domain.cells.begin(), domain.cells.end(), col); + if (lb != domain.cells.end()) { + //std::cerr << "row = " << cell << " col = " << col <<'\n'; + *iter = 0.0; + } + } + } + } + Dune::storeMatrixMarket(Jcopy, "aspinjacoffdiag.mm"); + // Set the residual to -Jcopy * aspin_residual, which will then be (D - J)*(x^k - x^{k+1}) + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + ebosResid = 0.0; + Jcopy.mmv(aspin_residual, ebosResid); + } +#endif + + // Check convergence before solving. + { + auto convrep = getConvergence(timer, -iteration, residual_norms); + // report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();; + // ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); + // convergence_reports_.back().report.push_back(std::move(convrep)); + + // // Throw if any NaN or too large residual found. + // if (severity == ConvergenceReport::Severity::NotANumber) { + // OPM_THROW(NumericalProblem, "NaN residual found!"); + // } else if (severity == ConvergenceReport::Severity::TooLarge) { + // OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); + // } + } + // report.convergence_check_time += perfTimer.stop(); + // residual_norms_history_.push_back(residual_norms); + + // if (report.converged) { + // return report; + // } + + + // Solve the linear system. + linear_solve_setup_time_ = 0.0; + BVector x(nc); + try { + solveJacobianSystem(x); + report.linear_solve_setup_time += linear_solve_setup_time_; + report.linear_solve_time += perfTimer.stop(); + report.total_linear_iterations += linearIterationsLastSolve(); + } + catch (...) { + report.linear_solve_setup_time += linear_solve_setup_time_; + report.linear_solve_time += perfTimer.stop(); + report.total_linear_iterations += linearIterationsLastSolve(); + + failureReport_ += report; + throw; // re-throw up + } + + perfTimer.reset(); + perfTimer.start(); + + // ----------- Update solution ----------- + + // Compute the nonlinear update. + + // apply the Schur compliment of the well model to the reservoir linearized + // equations + // wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), + // ebosSimulator().model().linearizer().residual()); + + + // handling well state update before oscillation treatment is a decision based + // on observation to avoid some big performance degeneration under some circumstances. + // there is no theorectical explanation which way is better for sure. + // wellModel().wellState() = initial_wellstate; + // ebosSimulator_.model().newtonMethod().setIterationIndex(iteration); + // wellModel().beginIteration(); + + // Before we can do postSolve() we must take into account that the current + // state in the wells are from after the local solves, not from the start + // of the iteration. + // Instead of passing x (where state_{i+1} = state_i + x, i.e. x is the + // difference from one global iteration to the next) we should pass + // y, with y = state_{i+1} - state_after_local_solve. + // BVector y(nc); + // for (int c = 0; c < nc; ++c) { + // y[c] = solution[c] - locally_solved[c]; + // } + // wellModel().postSolve(y); + + if (param_.use_update_stabilization_) { + // Stabilize the nonlinear update. + bool isOscillate = false; + bool isStagnate = false; + nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate); + if (isOscillate) { + current_relaxation_ -= nonlinear_solver.relaxIncrement(); + current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); + if (terminalOutputEnabled()) { + std::string msg = " Oscillating behavior detected: Relaxation set to " + + std::to_string(current_relaxation_); + OpmLog::info(msg); + } + } + nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_); + } + + // Apply the update, with considering model-dependent limitations and + // chopping of the update. This update must be based on the initial solution, + // not the solution after the local solve. + // ebosSimulator().model().solution(0) = initial_solution; + updateSolution(x); + + + + + // Before we can do postSolve() we must take into account that the current + // state in the wells are from after the local solves, not from the start + // of the iteration. + // Instead of passing x (where state_{i+1} = state_i - x, i.e. x is the + // difference from one global iteration to the next) we should pass + // y, with y = state_after_local_solve - state_{i+1}. + // BVector y(nc); + // for (int c = 0; c < nc; ++c) { + // y[c] = locally_solved[c] - solution[c]; + // } + // wellModel().postSolve(y); + wellModel().postSolve(x); + + + + + report.update_time += perfTimer.stop(); + + //wellModel().logPrimaryVars(); + + // Extra debug output. + { + std::ostringstream os; + os << " ** Well rates:"; + for (int w = 0; w < wellModel().wellState().numWells(); ++w) { + os << " | "; + const auto& wr = wellModel().wellState().wellRates(w); + for (const double r : wr) { + os << ' ' << r; + } + } + OpmLog::debug(os.str()); + } + + + return report; + } + + + + + std::pair + solveLocal(const Domain& domain, + const SimulatorTimerInterface& timer, + const int global_iteration, + const bool initial_assembly_required = false) + { + SimulatorReportSingle report; + Dune::Timer solveTimer; + solveTimer.start(); + Dune::Timer detailTimer; + + // TODO: the following assertion will fail, does it indicate a programming error? + // assert(ebosSimulator_.model().newtonMethod().numIterations() == 0); + ebosSimulator_.model().newtonMethod().setIterationIndex(0); + //ebosSimulator_.problem().beginIteration(); + + // When called, if assembly has already been performed + // with the initial values, we only need to check + // for local convergence. Otherwise, we must do a local + // assembly. + int iter = 0; + if (initial_assembly_required) { + OpmLog::debug("solveLocal() initial assembly START"); + //wellModel().logPrimaryVars(); + detailTimer.start(); + ebosSimulator_.model().newtonMethod().setIterationIndex(iter); + // TODO: we should have a beginIterationLocal function() + // only handling the well model for now + // ebosSimulator_.problem().beginIteration(); + ebosSimulator_.problem().wellModel().assembleLocal(ebosSimulator_.model().newtonMethod().numIterations(), + ebosSimulator_.timeStepSize(), + domain); + // Assemble reservoir locally. + report += assembleReservoirLocal(domain, iter); + report.assemble_time += detailTimer.stop(); + OpmLog::debug("solveLocal() initial assembly END"); + //wellModel().logPrimaryVars(); + } + detailTimer.reset(); + detailTimer.start(); + std::vector resnorms; + auto convreport = getLocalConvergence(domain, timer, 0, resnorms); + if (convreport.converged()) { + // TODO: set more info, timing etc. + report.converged = true; + // Exit, but not if we are on the first global iteration + // (if so, we require at least one iteration). + //if (global_iteration > 0) { + return { report, convreport }; + //} + } else { + // std::ostringstream os; + // os << "Convergence data for first local iteration:\n" + // << convreport; + // OpmLog::debug(os.str()); + } + //report.convergence_check_time += detailTimer.stop(); + + // We have already assembled for the first iteration, + // but not done the Schur complement for the wells yet. + detailTimer.reset(); + detailTimer.start(); + wellModel().linearizeDomain(domain, + ebosSimulator().model().linearizer().jacobian(), + ebosSimulator().model().linearizer().residual()); + const double tt1 = detailTimer.stop(); + report.assemble_time += tt1; + report.assemble_time_well += tt1; + //wellModel().logPrimaryVars(); + + // Local Newton loop. + const int max_iter = param_.max_local_solve_iterations_; + bool converged = false; + do { + // Solve local linear system. + // Note that x has full size, we expect it to be nonzero only for in-domain cells. + const int nc = grid_.size(0); + BVector x(nc); + detailTimer.reset(); + detailTimer.start(); + solveLocalJacobianSystem(domain, x); + wellModel().postSolveLocal(x, domain); + report.linear_solve_time += detailTimer.stop(); + report.linear_solve_setup_time += linear_solve_setup_time_; + report.total_linear_iterations = linearIterationsLastSolve(); + //wellModel().logPrimaryVars(); + + // Update local solution. // TODO: x is still full size, should we optimize it? + detailTimer.reset(); + detailTimer.start(); + updateDomainSolution(domain, x); + report.update_time += detailTimer.stop(); + // Note: endIteration does not do anything + // ebosSimulator_.problem().endIteration(); + + // Assemble well and reservoir. + detailTimer.reset(); + detailTimer.start(); + ++iter; + ebosSimulator_.model().newtonMethod().setIterationIndex(iter); + // TODO: we should have a beginIterationLocal function() + // only handling the well model for now + // ebosSimulator_.problem().beginIteration(); + // Assemble reservoir locally. + ebosSimulator_.problem().wellModel().assembleLocal(ebosSimulator_.model().newtonMethod().numIterations(), + ebosSimulator_.timeStepSize(), + domain); + report += assembleReservoirLocal(domain, iter); + report.assemble_time += detailTimer.stop(); + + + // Check for local convergence. + detailTimer.reset(); + detailTimer.start(); + convreport = getLocalConvergence(domain, timer, iter, resnorms); + // report.convergence_check_time += detailTimer.stop(); + //wellModel().logPrimaryVars(); + // Dune::writeMatrixToMatlab(ebosSimulator_.model().linearizer().jacobian().istlMatrix(), "beforeWellModelLinearize.matlab"); + + // apply the Schur compliment of the well model to the reservoir linearized + // equations + detailTimer.reset(); + detailTimer.start(); + wellModel().linearizeDomain(domain, + ebosSimulator().model().linearizer().jacobian(), + ebosSimulator().model().linearizer().residual()); + const double tt2 = detailTimer.stop(); + report.assemble_time += tt2; + report.assemble_time_well += tt2; + //wellModel().logPrimaryVars(); + + } while (!convreport.converged() && iter <= max_iter); + + ebosSimulator_.problem().endIteration(); + + if (!convreport.converged()) { + // std::ostringstream os; + // os << "Convergence data for unconverged local solution:\n" + // << convreport; + // OpmLog::debug(os.str()); + } + + report.converged = convreport.converged(); + report.total_newton_iterations = iter; + report.total_linearizations = iter; + report.total_time = solveTimer.stop(); + // TODO: set more info, timing etc. + return { report, convreport }; + } + + + + void printIf(int c, double x, double y, double eps, std::string type) { if (std::abs(x-y) > eps) { std::cout << type << " " <(Details::extractMatrix(main_matrix, domain.cells)); + } + auto& jac = *domain_matrices_[domain.index]; + auto res = Details::extractVector(ebosSimulator_.model().linearizer().residual(), domain.cells); + auto x = res; + + // set initial guess + global_x = 0.0; + x = 0.0; + + auto& linsolver = domain_linsolvers_[domain.index]; + + linsolver.prepare(jac, res); + linear_solve_setup_time_ = perfTimer.stop(); + linsolver.setResidual(res); + //linsolver.setMatrix(jac); + linsolver.solve(x); + + Details::setGlobal(x, domain.cells, global_x); + } + /// Solve the Jacobian system Jx = r where J is the Jacobian and /// r is the residual. void solveJacobianSystem(BVector& x) { - auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); + auto& ebosJac = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); // set initial guess @@ -619,12 +1569,33 @@ namespace Opm { // account for parallelization properly. since the residual of ECFV // discretizations does not need to be synchronized across processes to be // consistent, this is not relevant for OPM-flow... - ebosSolver.setMatrix(ebosJac); + //ebosSolver.setMatrix(ebosJac); ebosSolver.solve(x); } + /// Apply an update to the primary variables. + void updateDomainSolution(const Domain& domain, const BVector& dx) + { + auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod(); + SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0); + + ebosNewtonMethod.update_(/*nextSolution=*/solution, + /*curSolution=*/solution, + /*update=*/dx, + /*resid=*/dx, + domain.cells); // the update routines of the black + // oil model do not care about the + // residual + + // if the solution is updated, the intensive quantities need to be recalculated + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } + + + + /// Apply an update to the primary variables. void updateSolution(const BVector& dx) { @@ -719,7 +1690,8 @@ namespace Opm { /// of the cells associated with a numerical aquifer. std::tuple localConvergenceData(std::vector& R_sum, std::vector& maxCoeff, - std::vector& B_avg) + std::vector& B_avg, + std::vector& maxCoeffCell) { OPM_TIMEBLOCK(localConvergenceData); double pvSumLocal = 0.0; @@ -734,9 +1706,11 @@ namespace Opm { OPM_BEGIN_PARALLEL_TRY_CATCH(); for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + // elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + // const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); @@ -759,7 +1733,140 @@ namespace Opm { const auto R2 = ebosResid[cell_idx][compIdx]; R_sum[ compIdx ] += R2; - maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue ); + const double Rval = std::abs( R2 ) / pvValue; + if (Rval > maxCoeff[ compIdx ]) { + maxCoeff[ compIdx ] = Rval; + maxCoeffCell[ compIdx ] = cell_idx; + } + } + + if constexpr (has_solvent_) { + B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + const auto R2 = ebosResid[cell_idx][contiSolventEqIdx]; + R_sum[ contiSolventEqIdx ] += R2; + maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_extbo_) { + B_avg[ contiZfracEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiZfracEqIdx]; + R_sum[ contiZfracEqIdx ] += R2; + maxCoeff[ contiZfracEqIdx ] = std::max( maxCoeff[ contiZfracEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_polymer_) { + B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx]; + R_sum[ contiPolymerEqIdx ] += R2; + maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_foam_) { + B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiFoamEqIdx]; + R_sum[ contiFoamEqIdx ] += R2; + maxCoeff[ contiFoamEqIdx ] = std::max( maxCoeff[ contiFoamEqIdx ], std::abs( R2 ) / pvValue ); + } + if constexpr (has_brine_) { + B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiBrineEqIdx]; + R_sum[ contiBrineEqIdx ] += R2; + maxCoeff[ contiBrineEqIdx ] = std::max( maxCoeff[ contiBrineEqIdx ], std::abs( R2 ) / pvValue ); + } + + if constexpr (has_polymermw_) { + static_assert(has_polymer_); + + B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight + // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations + // TODO: there should be a more general way to determine the scaling-down coefficient + const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.; + R_sum[contiPolymerMWEqIdx] += R2; + maxCoeff[contiPolymerMWEqIdx] = std::max( maxCoeff[contiPolymerMWEqIdx], std::abs( R2 ) / pvValue ); + } + + if constexpr (has_energy_) { + B_avg[ contiEnergyEqIdx ] += 1.0; + const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx]; + R_sum[ contiEnergyEqIdx ] += R2; + maxCoeff[ contiEnergyEqIdx ] = std::max( maxCoeff[ contiEnergyEqIdx ], std::abs( R2 ) / pvValue ); + } + + } + + OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm()); + + // compute local average in terms of global number of elements + const int bSize = B_avg.size(); + for ( int i = 0; i& R_sum, + std::vector& maxCoeff, + std::vector& B_avg, + std::vector& maxCoeffCell) + { + double pvSumLocal = 0.0; + const auto& ebosModel = ebosSimulator_.model(); + const auto& ebosProblem = ebosSimulator_.problem(); + + const auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + + ElementContext elemCtx(ebosSimulator_); + const auto& gridView = domain.view; + const auto& elemEndIt = gridView.template end(); + OPM_BEGIN_PARALLEL_TRY_CATCH(); + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; + ++elemIt) + { + if (elemIt->partitionType() != Dune::InteriorEntity) { + continue; + } + const auto& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + // elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + + // const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + auto* intQuantsPtr = ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0); + // const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + if (intQuantsPtr == nullptr) { + OpmLog::debug("** no cached intensive quantities for cell " + std::to_string(cell_idx)); + elemCtx.updatePrimaryIntensiveQuantities(0); + intQuantsPtr = ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0); + assert(intQuantsPtr != nullptr); + } + const auto& intQuants = *intQuantsPtr; + + const auto& fs = intQuants.fluidState(); + + const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); + pvSumLocal += pvValue; + + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) + { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + + B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value(); + const auto R2 = ebosResid[cell_idx][compIdx]; + + R_sum[ compIdx ] += R2; + const double Rval = std::abs( R2 ) / pvValue; + if (Rval > maxCoeff[ compIdx ]) { + maxCoeff[ compIdx ] = Rval; + maxCoeffCell[ compIdx ] = cell_idx; + } } if constexpr (has_solvent_) { @@ -842,12 +1949,43 @@ namespace Opm { const int bSize = B_avg.size(); for ( int i = 0; i& B_avg, double dt) + { + double errorPV{}; + const auto& ebosModel = ebosSimulator_.model(); + const auto& ebosProblem = ebosSimulator_.problem(); + const auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + + for (const int cell_idx : domain.cells) + { + const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); + const auto& cellResidual = ebosResid[cell_idx]; + bool cnvViolated = false; + + for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) + { + using std::fabs; + Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue; + cnvViolated = cnvViolated || (fabs(CNV) > param_.tolerance_cnv_); + } + + if (cnvViolated) + { + errorPV += pvValue; + } + } + return errorPV; } + /// \brief Compute the total pore volume of cells violating CNV that are not part /// of a numerical aquifer. double computeCnvErrorPv(const std::vector& B_avg, double dt) @@ -870,7 +2008,7 @@ namespace Opm { continue; } elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + // elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); const auto& cellResidual = ebosResid[cell_idx]; @@ -894,6 +2032,119 @@ namespace Opm { return grid_.comm().sum(errorPV); } + ConvergenceReport getLocalReservoirConvergence(const double reportTime, + const double dt, + const int iteration, + const Domain& domain, + std::vector& B_avg, + std::vector& residual_norms) + { + typedef std::vector< Scalar > Vector; + + const int numComp = numEq; + Vector R_sum(numComp, 0.0 ); + Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() ); + std::vector maxCoeffCell(numComp, -1); + const double pvSum = localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell); + + // compute global sum and max of quantities + // const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal, + // R_sum, maxCoeff, B_avg); + + auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt); + cnvErrorPvFraction /= pvSum; + + const double tol_mb = param_.local_tolerance_scaling_mb_ * param_.tolerance_mb_; + // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0. + // For each iteration, we need to determine whether to use the relaxed CNV tolerance. + // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0. + const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.min_strict_cnv_iter_; + // Tighter bound for local convergence should increase the + // likelyhood of: local convergence => global convergence + const double tol_cnv = param_.local_tolerance_scaling_cnv_ + * (use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_); + + // Finish computation + std::vector CNV(numComp); + std::vector mass_balance_residual(numComp); + for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + { + CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx]; + mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum; + residual_norms.push_back(CNV[compIdx]); + } + + // Create convergence report. + ConvergenceReport report; + using CR = ConvergenceReport; + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] }; + CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance, + CR::ReservoirFailure::Type::Cnv }; + double tol[2] = { tol_mb, tol_cnv }; + int cell[2] = { -1, maxCoeffCell[compIdx] }; // No cell associated with MB failures. + for (int ii : {0, 1}) { + if (std::isnan(res[ii])) { + report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});//, cell[ii], res[ii]}); + if ( terminal_output_ ) { + OpmLog::debug("NaN residual for " + compNames_.name(compIdx) + " equation."); + } + } else if (res[ii] > maxResidualAllowed()) { + report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});//, cell[ii], res[ii]}); + if ( terminal_output_ ) { + OpmLog::debug("Too large residual for " + compNames_.name(compIdx) + " equation."); + } + } else if (res[ii] < 0.0) { + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});//, cell[ii], res[ii]}); + if ( terminal_output_ ) { + OpmLog::debug("Negative residual for " + compNames_.name(compIdx) + " equation."); + } + } else if (res[ii] > tol[ii]) { + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});//, cell[ii], res[ii]}); + } + } + } + + // Output of residuals. + if ( terminal_output_ ) + { + // Only rank 0 does print to std::cout + if (iteration == 0) { + std::string msg = fmt::format("Domain {}, size {}, containing cell {}\n| Iter", + domain.index, domain.cells.size(), domain.cells[0]); + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + msg += " MB("; + msg += compNames_.name(compIdx)[0]; + msg += ") "; + } + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + msg += " CNV("; + msg += compNames_.name(compIdx)[0]; + msg += ") "; + } + OpmLog::debug(msg); + } + std::ostringstream ss; + ss << "| "; + const std::streamsize oprec = ss.precision(3); + const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); + ss << std::setw(4) << iteration; + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + ss << std::setw(11) << mass_balance_residual[compIdx]; + } + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + ss << std::setw(11) << CNV[compIdx]; + } + ss.precision(oprec); + ss.flags(oflags); + OpmLog::debug(ss.str()); + } + + return report; + } + + + ConvergenceReport getReservoirConvergence(const double reportTime, const double dt, const int iteration, @@ -906,7 +2157,8 @@ namespace Opm { const int numComp = numEq; Vector R_sum(numComp, 0.0 ); Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() ); - const auto [ pvSumLocal, numAquiferPvSumLocal] = localConvergenceData(R_sum, maxCoeff, B_avg); + std::vector maxCoeffCell(numComp, -1); + const auto [ pvSumLocal, numAquiferPvSumLocal] = localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell); // compute global sum and max of quantities const auto [ pvSum, numAquiferPvSum ] = @@ -942,24 +2194,25 @@ namespace Opm { CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance, CR::ReservoirFailure::Type::Cnv }; double tol[2] = { tol_mb, tol_cnv }; + int cell[2] = { -1, maxCoeffCell[compIdx] }; // No cell associated with MB failures. for (int ii : {0, 1}) { if (std::isnan(res[ii])) { - report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});//, cell[ii], res[ii]}); if ( terminal_output_ ) { OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation."); } } else if (res[ii] > maxResidualAllowed()) { - report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});//, cell[ii], res[ii]}); if ( terminal_output_ ) { OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation."); } } else if (res[ii] < 0.0) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});//, cell[ii], res[ii]}); if ( terminal_output_ ) { OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation."); } } else if (res[ii] > tol[ii]) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});//, cell[ii], res[ii]}); } report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]); } @@ -1001,6 +2254,22 @@ namespace Opm { return report; } + ConvergenceReport getLocalConvergence(const Domain& domain, + const SimulatorTimerInterface& timer, + const int iteration, + std::vector& residual_norms) + { + std::vector B_avg(numEq, 0.0); + auto report = getLocalReservoirConvergence(timer.simulationTimeElapsed(), + timer.currentStepLength(), + iteration, + domain, + B_avg, + residual_norms); + report += wellModel().getLocalWellConvergence(domain, B_avg, false); + return report; + } + /// Compute convergence based on total mass balance (tol_mb) and maximum /// residual mass balance (tol_cnv). /// \param[in] timer simulation timer @@ -1059,6 +2328,10 @@ namespace Opm { const SimulatorReportSingle& failureReport() const { return failureReport_; } + /// return the statistics if the nonlinearIteration() method failed + const SimulatorReportSingle& localAccumulatedReports() const + { return local_reports_accumulated_; } + const std::vector& stepReports() const { return convergence_reports_; @@ -1086,6 +2359,7 @@ namespace Opm { ModelParameters param_; SimulatorReportSingle failureReport_; + SimulatorReportSingle local_reports_accumulated_; // Well Model BlackoilWellModel& well_model_; @@ -1101,6 +2375,9 @@ namespace Opm { std::vector convergence_reports_; ComponentName compNames_{}; + std::vector domains_; + std::vector> domain_matrices_; + std::vector domain_linsolvers_; public: /// return the StandardWells object diff --git a/opm/simulators/flow/FlowMainEbos.hpp b/opm/simulators/flow/FlowMainEbos.hpp index 66c6d2a2057..b24f11892e8 100644 --- a/opm/simulators/flow/FlowMainEbos.hpp +++ b/opm/simulators/flow/FlowMainEbos.hpp @@ -460,6 +460,13 @@ void handleExtraConvergenceOutput(SimulatorReport& report, { SimulatorReport report = simulator_->run(*simtimer_); runSimulatorAfterSim_(report); + { + SimulatorReportSingle rep_local = simulator_->model().localAccumulatedReports(); + std::ostringstream os; + os << "=== Accumulated local solve data ===\n"; + rep_local.reportFullyImplicit(os); + OpmLog::info(os.str()); + } return report.success.exit_status; } From ccc75268688d56d9a3c595720b47706054f2a067 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 11 Jan 2023 17:26:54 +0100 Subject: [PATCH 09/20] Use a single partition, to verify against base runs. --- opm/simulators/flow/aspinPartition.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/simulators/flow/aspinPartition.cpp b/opm/simulators/flow/aspinPartition.cpp index 05dfb8096ab..d30886b73a9 100644 --- a/opm/simulators/flow/aspinPartition.cpp +++ b/opm/simulators/flow/aspinPartition.cpp @@ -79,8 +79,8 @@ namespace std::pair, int> partitionCells(const int num_cells) { - // const std::string method = "simple"; - const std::string method = "file"; + const std::string method = "simple"; + // const std::string method = "file"; if (method == "simple") { return partitionCellsSimple(num_cells, 1); } else if (method == "file") { From 7a06bc7c17e58b2a2d73603c280c233398ecfb13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 13 Jan 2023 09:27:53 +0100 Subject: [PATCH 10/20] Bugfix: avoid static variable, fixing case with multiple linear solvers. --- opm/simulators/linalg/ISTLSolverEbos.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index c7c2656db63..f4d9194ea97 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -239,7 +239,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, iterations_( 0 ), calls_( 0 ), converged_(false), - matrix_(), + matrix_(nullptr), parameters_(parameters) { initialize(); @@ -252,7 +252,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, iterations_( 0 ), calls_( 0 ), converged_(false), - matrix_() + matrix_(nullptr) { parameters_.template init(simulator_.vanguard().eclState().getSimulationConfig().useCPR()); initialize(); @@ -341,7 +341,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverEbosPrepare); - static bool firstcall = true; + const bool firstcall = (matrix_ == nullptr); #if HAVE_MPI if (firstcall) { const size_t size = M.N(); @@ -385,8 +385,6 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #else prepareFlexibleSolver(); #endif - - firstcall = false; } From 27e5509c5c872da60dc1b2322065541b4ce4f755 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 16 Jan 2023 09:01:28 +0100 Subject: [PATCH 11/20] [WIP] Use 25 partitions, i.e. for the layer 10 test. --- opm/simulators/flow/aspinPartition.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/aspinPartition.cpp b/opm/simulators/flow/aspinPartition.cpp index d30886b73a9..168770eda0e 100644 --- a/opm/simulators/flow/aspinPartition.cpp +++ b/opm/simulators/flow/aspinPartition.cpp @@ -82,7 +82,7 @@ partitionCells(const int num_cells) const std::string method = "simple"; // const std::string method = "file"; if (method == "simple") { - return partitionCellsSimple(num_cells, 1); + return partitionCellsSimple(num_cells, 25); } else if (method == "file") { return partitionCellsFromFile(num_cells); } else { From 9e3ac6748a93788fd65b2a8ad04814857b807bc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 16 Jan 2023 11:27:23 +0100 Subject: [PATCH 12/20] Avoid clearing derivatives when pressure diff is 0. --- ebos/eclfluxmodule.hh | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 51d50cb6113..566cb851910 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -389,14 +389,16 @@ public: // of threshold pressure is a quite big hack that only makes sense for ECL // datasets. (and even there, its physical justification is quite // questionable IMO.) - if (std::abs(Toolbox::value(pressureDifference)) > thpres) { - if (pressureDifference < 0.0) - pressureDifference += thpres; - else - pressureDifference -= thpres; - } - else { - pressureDifference = 0.0; + if (thpres > 0.0) { + if (std::abs(Toolbox::value(pressureDifference)) > thpres) { + if (pressureDifference < 0.0) + pressureDifference += thpres; + else + pressureDifference -= thpres; + } + else { + pressureDifference = 0.0; + } } } From 1806ee5ef2193dc68e887bfb68f7adbe02a3dec3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 18 Jan 2023 13:26:38 +0100 Subject: [PATCH 13/20] Use simplified interface for domain-specific linearizer methods. --- opm/simulators/flow/BlackoilModelEbos.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 02e07e8972b..046e0ee4974 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -1419,9 +1419,9 @@ namespace Opm { // ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx); // ebosSimulator_.problem().beginIteration(); // Need to set residual and jacobian to zero in the domain. - ebosSimulator_.model().linearizer().resetSystem(domain.view, domain.interior); + ebosSimulator_.model().linearizer().resetSystem(domain); // Call the domain-dependent linearization. - ebosSimulator_.model().linearizer().linearizeDomain(domain.view, domain.interior); + ebosSimulator_.model().linearizer().linearizeDomain(domain); // ebosSimulator_.problem().endIteration(); return wellModel().lastReport(); From 3ddfc2cacd8ed722181d676ca3c6ac848d6e6faa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 18 Jan 2023 14:21:15 +0100 Subject: [PATCH 14/20] HACK to do init and restart, but incorrect with drsdt etc. --- ebos/eclproblem.hh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index e41959d0068..9aea91b8dac 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -1698,6 +1698,9 @@ public: */ bool recycleFirstIterationStorage() const { + if (this->simulator().timeStepIndex() == 0) + return true; + int episodeIdx = this->episodeIndex(); return !this->drsdtActive_(episodeIdx) && !this->drvdtActive_(episodeIdx) && From 789b459815146eb110d1557eaa64a10899e21c20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 24 Jan 2023 08:48:26 +0100 Subject: [PATCH 15/20] [WIP] Use Newton-first, and partition from file. --- opm/simulators/flow/BlackoilModelEbos.hpp | 3 +++ opm/simulators/flow/aspinPartition.cpp | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 046e0ee4974..291258fbd61 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -478,6 +478,9 @@ namespace Opm { convergence_reports_.back().report.reserve(11); } + if (iteration == 0) { + return nonlinearIterationNewton(iteration, timer, nonlinear_solver); + } if (param_.nonlinear_solver_ == "aspin" || param_.nonlinear_solver_ == "nldd" || param_.nonlinear_solver_ == "purelocal") { return nonlinearIterationAspin(iteration, timer, nonlinear_solver); } else { diff --git a/opm/simulators/flow/aspinPartition.cpp b/opm/simulators/flow/aspinPartition.cpp index 168770eda0e..6927b721342 100644 --- a/opm/simulators/flow/aspinPartition.cpp +++ b/opm/simulators/flow/aspinPartition.cpp @@ -79,8 +79,8 @@ namespace std::pair, int> partitionCells(const int num_cells) { - const std::string method = "simple"; - // const std::string method = "file"; + //const std::string method = "simple"; + const std::string method = "file"; if (method == "simple") { return partitionCellsSimple(num_cells, 25); } else if (method == "file") { From 1066f8ab89c70958696e837f3f82d89995fe76f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 6 Feb 2023 12:58:58 +0100 Subject: [PATCH 16/20] Follow API update of init() call for linear solver params. --- opm/simulators/flow/BlackoilModelEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 291258fbd61..aaf33596a01 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -395,7 +395,7 @@ namespace Opm { // parallel structures appropriate for the full grid // only. This must be addressed before going parallel. FlowLinearSolverParameters param; - param.template init(); + param.template init(ebosSimulator_.vanguard().eclState().getSimulationConfig().useCPR()); // Override solver type with umfpack if small domain. // Otherwise hardcode to ILU0 if (domains_[index].cells.size() < 200) { From 9a5338d3154be0e420b3961d29434fdb8f653b6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 6 Feb 2023 13:15:42 +0100 Subject: [PATCH 17/20] Move LinearSolver parameter and fix duplication. --- opm/simulators/linalg/FlowLinearSolverParameters.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 23803d16352..7c7acfa8cfc 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -301,11 +301,10 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge"); EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types"); - EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cpr (an alias for cpr_trueimpes), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition, "Write the JSON definition of the linear solver setup to the DBG wfle."); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate, 4: recreated every CprReuseInterval"); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseInterval, "Reuse preconditioner interval. Used when CprReuseSetup is set to 4, then the preconditioner will be fully recreated instead of reused every N linear solve, where N is this parameter."); - EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'"); EWOMS_REGISTER_PARAM(TypeTag, std::string, AcceleratorMode, "Choose a linear solver, usage: '--accelerator-mode=[none|cusparse|opencl|amgcl|rocalution]'"); EWOMS_REGISTER_PARAM(TypeTag, int, BdaDeviceId, "Choose device ID for cusparseSolver or openclSolver, use 'nvidia-smi' or 'clinfo' to determine valid IDs"); EWOMS_REGISTER_PARAM(TypeTag, int, OpenclPlatformId, "Choose platform ID for openclSolver, use 'clinfo' to determine valid platform IDs"); From 9ed89a916add92340457e9a9bdfaa2f3fbd40898 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 May 2023 11:31:16 +0200 Subject: [PATCH 18/20] Account for renaming SubGridView -> SubGridPart. --- opm/simulators/flow/BlackoilModelEbos.hpp | 4 ++-- opm/simulators/flow/SubDomain.hpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index aaf33596a01..eccad057f39 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -40,7 +40,7 @@ #include #include -#include +#include #include #include #include @@ -382,7 +382,7 @@ namespace Opm { for (int ix : partitions[index]) { interior[ix] = true; } - Dune::SubGridView view(ebosSimulator_.vanguard().grid(), std::move(seeds[index])); + Dune::SubGridPart view(ebosSimulator_.vanguard().grid(), std::move(seeds[index])); domains_.emplace_back(index, std::move(partitions[index]), std::move(interior), std::move(view)); } diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp index f7aeb21eb6a..db09d823bf4 100644 --- a/opm/simulators/flow/SubDomain.hpp +++ b/opm/simulators/flow/SubDomain.hpp @@ -20,7 +20,7 @@ #ifndef OPM_SUBDOMAIN_HEADER_INCLUDED #define OPM_SUBDOMAIN_HEADER_INCLUDED -#include +#include #include @@ -33,8 +33,8 @@ namespace Opm int index; std::vector cells; std::vector interior; - Dune::SubGridView view; - SubDomain(const int i, std::vector&& c, std::vector&& in, Dune::SubGridView&& v) + Dune::SubGridPart view; + SubDomain(const int i, std::vector&& c, std::vector&& in, Dune::SubGridPart&& v) : index(i), cells(std::move(c)), interior(std::move(in)), view(std::move(v)) {} }; From 96ee59973628f2e254ae518e73a83e3c2fee1274 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 11 May 2023 14:30:44 +0200 Subject: [PATCH 19/20] No longer necessary to call resetSystem() explicitly. --- opm/simulators/flow/BlackoilModelEbos.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index eccad057f39..168b0221c9f 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -1397,7 +1397,6 @@ namespace Opm { // if (param_.enable_aspin_) { // ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx); // ebosSimulator_.problem().beginIteration(); - // ebosSimulator_.model().linearizer().resetSystem(); // for (const auto& domain : domains_) { // assembleReservoirLocal(domain, iterationIdx); // } @@ -1421,8 +1420,6 @@ namespace Opm { // -------- Mass balance equations -------- // ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx); // ebosSimulator_.problem().beginIteration(); - // Need to set residual and jacobian to zero in the domain. - ebosSimulator_.model().linearizer().resetSystem(domain); // Call the domain-dependent linearization. ebosSimulator_.model().linearizer().linearizeDomain(domain); // ebosSimulator_.problem().endIteration(); From 935f4cbb32acc6c902e953d61e4854243b0a8e56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 13 Mar 2023 13:40:36 +0100 Subject: [PATCH 20/20] Enable Zoltan-Based In-Process Subdomain Partitioning Tested only for sequential running (mpirun -np 1) at the moment. Intended to replace the original partitioning scheme based on a user-provided input file named partition.txt. There are two new user-controlled runtime parameters NumLocalDomains (--num-local-domains): Positive integer. The partitioner will divide the domain into this many non-overlapping subdomains and run the non-linear domain decomposition method on those LocalDomainsPartitioningImbalance (--local-domains-partitioning-imbalance`): Floating-point value >= 1 which controls Zoltan's imbalance tolerance. A value of 1.03 corresponds to a tolerance of three percent domain imbalance These parameters are not queried unless the run uses the NLDD non-linear solver. --- ebos/eclgenericcpgridvanguard.cc | 16 +++- ebos/eclgenericcpgridvanguard.hh | 10 ++ opm/simulators/flow/BlackoilModelEbos.hpp | 92 ++++++++++++++++--- .../flow/BlackoilModelParametersEbos.hpp | 26 +++++- 4 files changed, 127 insertions(+), 17 deletions(-) diff --git a/ebos/eclgenericcpgridvanguard.cc b/ebos/eclgenericcpgridvanguard.cc index 2fc206bf5ac..3e23eb41724 100644 --- a/ebos/eclgenericcpgridvanguard.cc +++ b/ebos/eclgenericcpgridvanguard.cc @@ -49,7 +49,7 @@ #include #include #include -#endif //HAVE_DUNE_FEM +#endif // HAVE_DUNE_FEM #include #include @@ -84,6 +84,20 @@ void EclGenericCpGridVanguard::releaseEquilGrid() } #if HAVE_MPI + +template +std::vector +EclGenericCpGridVanguard:: +partitionCells(const int numDomains, + const double imbalanceTolerance, + const std::vector& wells) const +{ + const auto* faceTrans = static_cast(nullptr); + + return this->grid_-> + zoltanPartitionWithoutScatter(&wells, faceTrans, numDomains, imbalanceTolerance); +} + template void EclGenericCpGridVanguard:: doLoadBalance_(const Dune::EdgeWeightMethod edgeWeightsMethod, diff --git a/ebos/eclgenericcpgridvanguard.hh b/ebos/eclgenericcpgridvanguard.hh index c873b410b83..fa9f98ce598 100644 --- a/ebos/eclgenericcpgridvanguard.hh +++ b/ebos/eclgenericcpgridvanguard.hh @@ -28,9 +28,12 @@ #define EWOMS_ECL_CP_GRID_GENERIC_VANGUARD_HH #include + #include +#include #include +#include #include #include @@ -115,6 +118,13 @@ public: return this->cell_part_; } +#if HAVE_MPI + std::vector + partitionCells(const int numDomains, + const double imbalanceTolerance, + const std::vector& wells) const; +#endif + protected: /*! * \brief Distribute the simulation grid over multiple processes diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 168b0221c9f..c32f6ea393f 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -90,6 +90,7 @@ #include #include #include +#include #include #include @@ -343,26 +344,26 @@ namespace Opm { void setupAspinDomains() { // Create partitions. - const int num_cells = detail::countLocalInteriorCells(grid_); - auto [partition_vector, num_partitions] = partitionCells(num_cells); + const auto& [partition_vector, num_domains] = + this->partitionCells(); // Scan through partitioning to get correct size for each. - std::vector sizes(num_partitions, 0); - for (int p : partition_vector) { + std::vector sizes(num_domains, 0); + for (const auto& p : partition_vector) { ++sizes[p]; } // Set up correctly sized vectors of entity seeds and of indices for each partition. using EntitySeed = typename Grid::template Codim<0>::EntitySeed; - std::vector> seeds(num_partitions); - std::vector> partitions(num_partitions); - for (int part = 0; part < num_partitions; ++part) { - seeds[part].resize(sizes[part]); - partitions[part].resize(sizes[part]); + std::vector> seeds(num_domains); + std::vector> partitions(num_domains); + for (int domain = 0; domain < num_domains; ++domain) { + seeds[domain].resize(sizes[domain]); + partitions[domain].resize(sizes[domain]); } // Iterate through grid once, setting the seeds of all partitions. - std::vector count(num_partitions, 0); + std::vector count(num_domains, 0); const auto beg = grid_.template leafbegin<0>(); const auto end = grid_.template leafend<0>(); int cell = 0; @@ -372,18 +373,26 @@ namespace Opm { partitions[p][count[p]] = cell; ++count[p]; } - assert(cell == num_cells); + //assert(cell == num_cells); assert(count == sizes); // Create the domains. - const int num_domains = partitions.size(); + //const int num_domains = partitions.size(); for (int index = 0; index < num_domains; ++index) { - std::vector interior(num_cells, false); + std::vector interior(partition_vector.size(), false); for (int ix : partitions[index]) { interior[ix] = true; } - Dune::SubGridPart view(ebosSimulator_.vanguard().grid(), std::move(seeds[index])); - domains_.emplace_back(index, std::move(partitions[index]), std::move(interior), std::move(view)); + + Dune::SubGridPart view { + ebosSimulator_.vanguard().grid(), + std::move(seeds[index]) + }; + + this->domains_.emplace_back(index, + std::move(partitions[index]), + std::move(interior), + std::move(view)); } // Set up container for the local system matrices. @@ -2398,6 +2407,19 @@ namespace Opm { } private: + template + struct HasZoltanPartitioning : public std::false_type {}; + + template + struct HasZoltanPartitioning< + GridType, + std::void_t().zoltanPartitionWithoutScatter + (std::declval*>(), + std::declval(), + std::declval(), + std::declval()))> + > : public std::true_type {}; + template bool isNumericalAquiferCell(const Dune::CpGrid& grid, const T& elem) { @@ -2424,6 +2446,46 @@ namespace Opm { double maxResidualAllowed() const { return param_.max_residual_allowed_; } double linear_solve_setup_time_; + std::pair, int> partitionCells() const + { + if constexpr (HasZoltanPartitioning::value) { + return this->partitionCellsZoltan(); + } + else { + return this->partitionCellsFallback(); + } + } + + std::pair, int> partitionCellsZoltan() const + { + const auto wells = this->ebosSimulator_.vanguard().schedule().getWellsatEnd(); + + auto partition_vector = this->grid_.zoltanPartitionWithoutScatter + (&wells, nullptr, this->param_.num_local_domains_, + this->param_.local_domain_partition_imbalance_); + + return this->countDomains(std::move(partition_vector)); + } + + std::pair, int> partitionCellsFallback() const + { + const int num_cells = detail::countLocalInteriorCells(this->grid_); + + return this->countDomains(partitionCells(num_cells)); + } + + std::pair, int> + countDomains(std::vector partition_vector) const + { + auto maxPos = std::max_element(partition_vector.begin(), + partition_vector.end()); + + const auto num_domains = (maxPos == partition_vector.end()) + ? 0 : *maxPos + 1; + + return { std::move(partition_vector), num_domains }; + } + public: std::vector wasSwitched_; }; diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 5e8b3e7234d..62969f9d788 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -192,6 +192,14 @@ template struct LocalToleranceScalingCnv { using type = UndefinedProperty; }; +template +struct NumLocalDomains { + using type = UndefinedProperty; +}; +template +struct LocalDomainsPartitioningImbalance { + using type = UndefinedProperty; +}; template struct DbhpMaxRel { using type = GetPropType; @@ -366,7 +374,16 @@ struct LocalToleranceScalingCnv { using type = GetPropType; static constexpr type value = 0.01; }; ->>>>>>> d1d200f5b (Add new parameters for aspin usage.) +template +struct NumLocalDomains { + using type = int; + static constexpr auto value = 0; +}; +template +struct LocalDomainsPartitioningImbalance { + using type = GetPropType; + static constexpr auto value = type{1.03}; +}; // if openMP is available, determine the number threads per process automatically. #if _OPENMP @@ -494,6 +511,9 @@ namespace Opm double local_tolerance_scaling_mb_; double local_tolerance_scaling_cnv_; + int num_local_domains_{0}; + double local_domain_partition_imbalance_{1.03}; + /// Construct from user parameters or defaults. BlackoilModelParametersEbos() { @@ -535,6 +555,8 @@ namespace Opm max_local_solve_iterations_ = EWOMS_GET_PARAM(TypeTag, int, MaxLocalSolveIterations); local_tolerance_scaling_mb_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingMb); local_tolerance_scaling_cnv_ = EWOMS_GET_PARAM(TypeTag, double, LocalToleranceScalingCnv); + num_local_domains_ = EWOMS_GET_PARAM(TypeTag, int, NumLocalDomains); + local_domain_partition_imbalance_ = std::max(1.0, EWOMS_GET_PARAM(TypeTag, double, LocalDomainsPartitioningImbalance)); deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); } @@ -580,6 +602,8 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, int, MaxLocalSolveIterations, "Max iterations for local solves with ASPIN or NLDD."); EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingMb, "Set lower than 1.0 to use stricter convergence tolerance for local solves."); EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalToleranceScalingCnv, "Set lower than 1.0 to use stricter convergence tolerance for local solves."); + EWOMS_REGISTER_PARAM(TypeTag, int, NumLocalDomains, "Number of local domains for ASPIN/NLDD solves."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LocalDomainsPartitioningImbalance, "Subdomain partitioning imbalance tolerance. 1.03 is 3 percent imbalance."); } }; } // namespace Opm