Map global index to local subdomain index in well operator

This commit is contained in:
jakobtorben 2024-05-08 13:01:33 +02:00
parent 5eea7055c3
commit 8ba83d4c7b
3 changed files with 44 additions and 4 deletions

View File

@ -244,9 +244,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).");
}

View File

@ -423,6 +423,9 @@ class WellContributions;
// Keep track of the domain of each well, if using subdomains.
std::map<std::string, int> well_domain_;
// Store the global index of the cells in the domain, if using sumdomains
std::vector<std::vector<int>> domains_cells_;
const Grid& grid() const
{ return simulator_.vanguard().grid(); }

View File

@ -1575,6 +1575,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]];
}
@ -1603,6 +1621,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]];
@ -2772,6 +2808,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