Merge pull request #1735 from andlaus/fix_wells_Schur

well model: fix application of Schur complements
This commit is contained in:
Tor Harald Sandve 2019-02-25 14:21:17 +01:00 committed by GitHub
commit 5163131f59
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 20 additions and 23 deletions

View File

@ -246,7 +246,7 @@ namespace Opm {
report.total_linearizations = 1; report.total_linearizations = 1;
try { try {
report += assemble(timer, iteration); report += assembleReservoir(timer, iteration);
report.assemble_time += perfTimer.stop(); report.assemble_time += perfTimer.stop();
} }
catch (...) { catch (...) {
@ -293,6 +293,11 @@ namespace Opm {
const int nc = UgGridHelpers::numCells(grid_); const int nc = UgGridHelpers::numCells(grid_);
BVector x(nc); BVector x(nc);
// apply the Schur compliment of the well model to the reservoir linearized
// equations
wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
ebosSimulator().model().linearizer().residual());
try { try {
solveJacobianSystem(x); solveJacobianSystem(x);
report.linear_solve_time += perfTimer.stop(); report.linear_solve_time += perfTimer.stop();
@ -360,13 +365,13 @@ namespace Opm {
/// \param[in] reservoir_state reservoir state variables /// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables /// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep /// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
SimulatorReport assemble(const SimulatorTimerInterface& timer, SimulatorReport assembleReservoir(const SimulatorTimerInterface& timer,
const int iterationIdx) const int iterationIdx)
{ {
// -------- Mass balance equations -------- // -------- Mass balance equations --------
ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx); ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
ebosSimulator_.problem().beginIteration(); ebosSimulator_.problem().beginIteration();
ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.model().linearizer().linearizeDomain();
ebosSimulator_.problem().endIteration(); ebosSimulator_.problem().endIteration();
return wellModel().lastReport(); return wellModel().lastReport();
@ -463,19 +468,6 @@ namespace Opm {
auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
auto& ebosResid = ebosSimulator_.model().linearizer().residual(); auto& ebosResid = ebosSimulator_.model().linearizer().residual();
// J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well
// with the reservoir and D is the wells itself.
// The full system is reduced to a number of cells X number of cells system via Schur complement
// A -= B^T D^-1 C
// If matrix_add_well_contribution is false, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter
// instead of A.
// The residual is modified similarly.
// r = [r, r_well], where r is the residual and r_well the well residual.
// r -= B^T * D^-1 r_well
wellModel().apply(ebosResid);
if (param_.matrix_add_well_contributions_) {
wellModel().addWellContributions(ebosJac.istlMatrix());
}
// set initial guess // set initial guess
x = 0.0; x = 0.0;

View File

@ -101,18 +101,23 @@ namespace Opm {
if (!localWellsActive()) if (!localWellsActive())
return; return;
// we don't what to add the schur complement if (!param_.matrix_add_well_contributions_) {
// here since it affects the getConvergence method // 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_) {
// r = r - duneC_^T * invDuneD_ * resWell_
well->apply(res);
}
return;
}
for (const auto& well: well_container_) { for (const auto& well: well_container_) {
if (param_.matrix_add_well_contributions_) well->addWellContributions(mat.istlMatrix());
well->addWellContributions(mat);
// applying the well residual to reservoir residuals // applying the well residual to reservoir residuals
// r = r - duneC_^T * invDuneD_ * resWell_ // r = r - duneC_^T * invDuneD_ * resWell_
well->apply(res); well->apply(res);
} }
*/
} }