ebos: make linearization recycling potentially useful

the missing piece was determining if the wells have changed between
report steps. This patch adds a simple way to determine this, but it
relies on low-level properties of opm-parser it does not
guarantee. (concretely, these details are that the same well objects
are returned in the same order if nothing changes. Since IMO this is a
pretty reasonable assumption, we use this approach instead of a more
complicated one until opm-parser provides a "change determination API"
for wells...)

note that this patch may increase the number of iterations required
for the simulation because the linear system of equations which is
solved in the first iteration of a time step actually corresponds to
the second to last solution of the previous time step. This means that
that linearization recycling usually only works well if the tolerance
of the Newton-Raphson solver is "sufficiently" low. ("sufficiently"
means that the linearization errors made due to using the "wrong"
solution for the first iteration can be neglected compared to the
differences because of the change of the solution in this iteration.)

therefore, use this feature at your own risk...
This commit is contained in:
Andreas Lauser
2015-03-31 16:26:06 +02:00
parent a6740b83ee
commit 15ad669d5d
2 changed files with 51 additions and 3 deletions

View File

@@ -388,7 +388,14 @@ public:
{
wellManager_.endTimeStep();
// write the summary information after each time step
summaryWriter_.write(wellManager_);
#ifndef NDEBUG
// in debug mode, we don't care about performance, so we check if the model does
// the right thing (i.e., the mass change inside the whole reservoir must be
// equivalent to the fluxes over the grid's boundaries plus the source rates
// specified by the problem)
this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
#endif // NDEBUG
}
@@ -398,10 +405,15 @@ public:
*/
void endEpisode()
{
std::cout << "Episode " << this->simulator().episodeIndex() + 1 << " finished.\n";
const auto& simulator = this->simulator();
const auto& eclState = simulator.gridManager().eclState();
auto& linearizer = this->model().linearizer();
int episodeIdx = simulator.episodeIndex();
// first, write the summary information ...
summaryWriter_.write(wellManager_);
bool wellsWillChange = wellManager_.wellsChanged(eclState, episodeIdx + 1);
linearizer.setLinearizationReusable(!wellsWillChange);
std::cout << "Episode " << episodeIdx + 1 << " finished.\n";
}
/*!

View File

@@ -518,6 +518,38 @@ public:
beginEpisode(simulator_.gridManager().eclState(), /*wasRestarted=*/true);
}
/*!
* \brief Returns true if something in a well changed compared to the previous report
* step.
*
* "Something" can either be the well topology (i.e., which grid blocks are contained
* in which well) or it can be a well parameter like the bottom hole pressure...
*/
bool wellsChanged(Opm::EclipseStateConstPtr eclState, int reportStepIdx) const
{
if (wellTopologyChanged_(eclState, reportStepIdx))
return true;
// this is slightly hacky because it assumes that the object which stores the set
// of wells which are relevant for a report step does not change if there are no
// changed well parameters. opm-parser does not guarantee this, but so far it
// seems to adhere to it...
auto deckSchedule = eclState->getSchedule();
if (deckSchedule->getTimeMap()->numTimesteps() <= reportStepIdx)
// for the "until the universe dies" episode, the wells don't change
return false;
const auto& curDeckWells = deckSchedule->getWells(reportStepIdx);
const auto& prevDeckWells = deckSchedule->getWells(reportStepIdx - 1);
for (int i = 0; i < curDeckWells.size(); ++i) {
if (curDeckWells[i] != prevDeckWells[i])
return true;
}
return false;
}
protected:
bool wellTopologyChanged_(Opm::EclipseStateConstPtr eclState, int reportStepIdx) const
{
@@ -528,6 +560,10 @@ protected:
}
auto deckSchedule = eclState->getSchedule();
if (deckSchedule->getTimeMap()->numTimesteps() <= reportStepIdx)
// for the "until the universe dies" episode, the wells don't change
return false;
const auto& curDeckWells = deckSchedule->getWells(reportStepIdx);
const auto& prevDeckWells = deckSchedule->getWells(reportStepIdx - 1);