diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index d1442441c..66c9b8aa1 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -242,9 +242,6 @@ namespace Opm { convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves. // TODO: remember to fix! if (param_.nonlinear_solver_ == "nldd") { - if (!param_.matrix_add_well_contributions_) { - OPM_THROW(std::runtime_error, "The --nonlinear-solver=nldd option can only be used with --matrix-add-well-contributions=true"); - } if (terminal_output) { OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd)."); } diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index b43f0a6b6..51615cda5 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -441,6 +441,9 @@ template class WellContributions; // Keep track of the domain of each well, if using subdomains. std::map well_domain_; + // Store the global index of the cells in the domain, if using sumdomains + std::vector> domains_cells_; + const Grid& grid() const { return simulator_.vanguard().grid(); } diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index f48641712..dc9048c23 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1745,6 +1745,24 @@ namespace Opm { auto cells = well->cells(); r_local.resize(cells.size()); + if (this->param_.nonlinear_solver_ == "nldd") { + // transfer global cells index to local subdomain cells index + + auto domain_cells = domains_cells_[well_domain_.at(well->name())]; + + // Assuming domain_cells is sorted + for (size_t i = 0; i < cells.size(); ++i) { + auto it = std::lower_bound(domain_cells.begin(), domain_cells.end(), cells[i]); + if (it != domain_cells.end() && *it == cells[i]) { + // Found the cell, get the index + auto local_index = std::distance(domain_cells.begin(), it); + cells[i] = local_index; + } else { + std::cerr << "Cell value " << cells[i] << " not found in domain_cells." << std::endl; + } + } + } + for (size_t i = 0; i < cells.size(); ++i) { r_local[i] = r[cells[i]]; } @@ -1773,6 +1791,24 @@ namespace Opm { x_local.resize(cells.size()); Ax_local.resize(cells.size()); + if (this->param_.nonlinear_solver_ == "nldd") { + // transfer global cells index to local subdomain cells index + + auto domain_cells = domains_cells_[well_domain_.at(well->name())]; + + // Assuming domain_cells is sorted + for (size_t i = 0; i < cells.size(); ++i) { + auto it = std::lower_bound(domain_cells.begin(), domain_cells.end(), cells[i]); + if (it != domain_cells.end() && *it == cells[i]) { + // Found the cell, get the index + auto local_index = std::distance(domain_cells.begin(), it); + cells[i] = local_index; + } else { + std::cerr << "Cell value " << cells[i] << " not found in domain_cells." << std::endl; + } + } + } + for (size_t i = 0; i < cells.size(); ++i) { x_local[i] = x[cells[i]]; Ax_local[i] = Ax[cells[i]]; @@ -2945,6 +2981,10 @@ namespace Opm { if (this->terminal_output_) { global_log.logMessages(); } - } + // Store the global index of the cells in the domain + for (const auto& domain : domains) { + domains_cells_.push_back(domain.cells); + } + } } // namespace Opm