- fixed extended communicator

- fixed scaling of well cpr
This commit is contained in:
hnil 2022-04-01 17:08:23 +02:00 committed by Atgeirr Flø Rasmussen
parent da572d1f60
commit 4975d5d9e7
2 changed files with 14 additions and 5 deletions

View File

@ -37,6 +37,7 @@ namespace Opm
IndexSet& indset_rw = commRW->indexSet();
const int max_nw = comm.communicator().max(nw);
const int num_proc = comm.communicator().size();
const int rank = comm.communicator().rank();
int glo_max = 0;
size_t loc_max = 0;
indset_rw.beginResize();
@ -53,7 +54,7 @@ namespace Opm
size_t local_ind = loc_max + 1;
for (int i = 0; i < nw; ++i) {
// need to set unique global number
const size_t v = global_max + max_nw * num_proc + i + 1;
const size_t v = global_max + max_nw * rank + i + 1;
// set to true to not have problems with higher levels if growing of domains is used
indset_rw.add(v, LocalIndex(local_ind, Dune::OwnerOverlapCopyAttributeSet::owner, true));
++local_ind;
@ -63,6 +64,7 @@ namespace Opm
// assume same communication pattern
commRW->remoteIndices().setNeighbours(comm.remoteIndices().getNeighbours());
commRW->remoteIndices().template rebuild<true>();
//commRW->clearInterfaces(); may need this for correct rebuild
}

View File

@ -2241,20 +2241,27 @@ namespace Opm
DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
// diag_block_transpose.solve(bweights, rhs);
//HACK due to template errors
double abs_max = 0;
for (size_t i = 0; i < blockSz; ++i) {
bweights[0][i] = 0;
for (size_t j = 0; j < blockSz; ++j) {
bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
}
abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
}
//inv_diag_block_transpose.mv(rhs[0], bweights[0]);
// NB how to scale to make it most symmetric
// double abs_max = *std::max_element(
// bweights.begin(), bweights.end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
// bweights /= std::fabs(abs_max);
//double abs_max = *std::max_element(
// bweights[0].begin(), bweights[0].end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
assert(abs_max>0.0);
for (size_t i = 0; i < blockSz; ++i) {
bweights[0][i] /= abs_max;
}
diagElem = 1.0/abs_max;
}
//
jacobian[welldof_ind][welldof_ind] = 1.0;
jacobian[welldof_ind][welldof_ind] = diagElem;
// set the matrix elements for well reservoir coupling
for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();