OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master' into PR/merge-flow-and-flow_cp

Browse Source
This commit is contained in:
Robert Kloefkorn 2015-06-30 13:39:03 +02:00
commit 5b51a5b7b1
9 changed files with 852 additions and 350 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@ -29,7 +29,9 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/BlackoilPropsAdInterface.cpp
opm/autodiff/ExtractParallelGridInformationToISTL.cpp
opm/autodiff/NewtonIterationBlackoilCPR.cpp
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
opm/autodiff/NewtonIterationBlackoilSimple.cpp
opm/autodiff/NewtonIterationUtilities.cpp
opm/autodiff/GridHelpers.cpp
opm/autodiff/ImpesTPFAAD.cpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
@ -111,7 +113,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/ImpesTPFAAD.hpp
opm/autodiff/NewtonIterationBlackoilCPR.hpp
opm/autodiff/NewtonIterationBlackoilInterface.hpp
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
opm/autodiff/NewtonIterationBlackoilSimple.hpp
opm/autodiff/NewtonIterationUtilities.hpp
opm/autodiff/NewtonSolver.hpp
opm/autodiff/NewtonSolver_impl.hpp
opm/autodiff/LinearisedBlackoilResidual.hpp

7
examples/flow.cpp
View File

@ -70,6 +70,7 @@
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@ -345,8 +346,10 @@ try
// Solver for Newton iterations.
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
if (param.getDefault("use_cpr", true)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information));
if (param.getDefault("use_interleaved", false)) {
fis_solver.reset(new NewtonIterationBlackoilInterleaved(param));
} else if (param.getDefault("use_cpr", true)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
}

42
opm/autodiff/BlackoilModelBase.hpp
View File

@ -307,9 +307,6 @@ namespace Opm {
std::vector<int>
variableWellStateIndices() const;
void
addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s);
SolutionState
variableStateExtractVars(const ReservoirState& x,
const std::vector<int>& indices,
@ -330,11 +327,6 @@ namespace Opm {
void
assembleMassBalanceEq(const SolutionState& state);
void
addWellControlEq(const SolutionState& state,
const WellState& xw,
const V& aliveWells);
void
solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
@ -342,20 +334,30 @@ namespace Opm {
WellState& well_state);
void
addWellEq(const SolutionState& state,
WellState& xw,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s);
computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s);
void
extraAddWellEq(const SolutionState& state,
const WellState& xw,
const std::vector<ADB>& cq_ps,
const std::vector<ADB>& cmix_s,
const ADB& cqt_is,
const std::vector<int>& well_cells);
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw);
void
addWellFluxEq(const std::vector<ADB>& cq_s,
const SolutionState& state);
void
addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
const SolutionState& state,
const WellState& xw);
void
addWellControlEq(const SolutionState& state,
const WellState& xw,
const V& aliveWells);
void updateWellControls(WellState& xw) const;

90
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -810,8 +810,10 @@ namespace detail {
solveWellEq(mob_perfcells, b_perfcells, state, well_state);
}
asImpl().addWellEq(state, well_state, mob_perfcells, b_perfcells, aliveWells, cq_s);
addWellContributionToMassBalanceEq(cq_s);
asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
asImpl().addWellFluxEq(cq_s, state);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
addWellControlEq(state, well_state, aliveWells);
}
@ -875,7 +877,9 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s)
BlackoilModelBase<Grid, Implementation>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
const SolutionState&,
const WellState&)
{
// Add well contributions to mass balance equations
const int nc = Opm::AutoDiffGrid::numCells(grid_);
@ -894,12 +898,11 @@ namespace detail {
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::addWellEq(const SolutionState& state,
WellState& xw,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s)
BlackoilModelBase<Grid, Implementation>::computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s)
{
if( ! wellsActive() ) return ;
@ -919,8 +922,6 @@ namespace detail {
// Perforation pressure
const ADB perfpressure = (wops_.w2p * state.bhp) + cdp;
std::vector<double> perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf);
xw.perfPress() = perfpressure_d;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;
@ -1013,13 +1014,6 @@ namespace detail {
cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
}
// WELL EQUATIONS
ADB qs = state.qs;
for (int phase = 0; phase < np; ++phase) {
qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
}
// check for dead wells (used in the well controll equations)
aliveWells = V::Constant(nw, 1.0);
for (int w = 0; w < nw; ++w) {
@ -1027,19 +1021,6 @@ namespace detail {
aliveWells[w] = 0.0;
}
}
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore)
V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
}
std::vector<double> cq_d(cq.data(), cq.data() + nperf*np);
xw.perfPhaseRates() = cq_d;
residual_.well_flux_eq = qs;
asImpl().extraAddWellEq(state, xw, cq_ps, cmix_s, cqt_is, well_cells);
}
@ -1047,14 +1028,43 @@ namespace detail {
template <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::extraAddWellEq(const SolutionState& /* state */,
const WellState& /* xw */,
const std::vector<ADB>& /* cq_ps */,
const std::vector<ADB>& /* cmix_s */,
const ADB& /* cqt_is */,
const std::vector<int>& /* well_cells */)
void BlackoilModelBase<Grid, Implementation>::updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw)
{
// Does nothing in this model.
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
}
xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
// Update the perforation pressures.
const V& cdp = well_perforation_pressure_diffs_;
const V perfpressure = (wops_.w2p * state.bhp.value().matrix()).array() + cdp;
xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
}
template <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::addWellFluxEq(const std::vector<ADB>& cq_s,
const SolutionState& state)
{
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
ADB qs = state.qs;
for (int phase = 0; phase < np; ++phase) {
qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
}
residual_.well_flux_eq = qs;
}
@ -1247,7 +1257,9 @@ namespace detail {
SolutionState wellSolutionState = state0;
variableStateExtractWellsVars(indices, vars, wellSolutionState);
asImpl().addWellEq(wellSolutionState, well_state, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
asImpl().addWellFluxEq(cq_s, wellSolutionState);
addWellControlEq(wellSolutionState, well_state, aliveWells);
converged = getWellConvergence(it);

299
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@ -25,25 +25,11 @@
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/utility/platform_dependent/disable_warnings.h>
// #include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/io.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
#if HAVE_UMFPACK
#include <Eigen/UmfPackSupport>
@ -51,6 +37,7 @@
#include <Eigen/SparseLU>
#endif
namespace Opm
{
@ -58,41 +45,6 @@ namespace Opm
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
namespace {
/// Eliminate a variable via Schur complement.
/// \param[in] eqs set of equations with Jacobians
/// \param[in] n index of equation/variable to eliminate.
/// \return new set of equations, one smaller than eqs.
/// Note: this method requires the eliminated variable to have the same size
/// as the equation in the corresponding position (that also will be eliminated).
/// It also required the jacobian block n of equation n to be diagonal.
std::vector<ADB> eliminateVariable(const std::vector<ADB>& eqs, const int n);
/// Recover that value of a variable previously eliminated.
/// \param[in] equation previously eliminated equation.
/// \param[in] partial_solution solution to the remainder system after elimination.
/// \param[in] n index of equation/variable that was eliminated.
/// \return solution to complete system.
V recoverVariable(const ADB& equation, const V& partial_solution, const int n);
/// Form an elliptic system of equations.
/// \param[in] num_phases the number of fluid phases
/// \param[in] eqs the equations
/// \param[out] A the resulting full system matrix
/// \param[out] b the right hand side
/// This function will deal with the first num_phases
/// equations in eqs, and return a matrix A for the full
/// system that has a elliptic upper left corner, if possible.
void formEllipticSystem(const int num_phases,
const std::vector<ADB>& eqs,
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
V& b);
} // anonymous namespace
@ -114,6 +66,8 @@ namespace Opm
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
@ -219,8 +173,7 @@ namespace Opm
// Copy solver output to dx.
std::copy(x.begin(), x.end(), dx.data());
if( hasWells )
{
if ( hasWells ) {
// Compute full solution using the eliminated equations.
// Recovery in inverse order of elimination.
dx = recoverVariable(elim_eqs[1], dx, np);
@ -229,6 +182,10 @@ namespace Opm
return dx;
}
const boost::any& NewtonIterationBlackoilCPR::parallelInformation() const
{
return parallelInformation_;
@ -236,242 +193,6 @@ namespace Opm
namespace
{
std::vector<ADB> eliminateVariable(const std::vector<ADB>& eqs, const int n)
{
// Check that the variable index to eliminate is within bounds.
const int num_eq = eqs.size();
const int num_vars = eqs[0].derivative().size();
if (num_eq != num_vars) {
OPM_THROW(std::logic_error, "eliminateVariable() requires the same number of variables and equations.");
}
if (n >= num_eq) {
OPM_THROW(std::logic_error, "Trying to eliminate variable from too small set of equations.");
}
// Schur complement of (A B ; C D) wrt. D is A - B*inv(D)*C.
// This is applied to all 2x2 block submatrices
// The right hand side is modified accordingly. bi = bi - B * inv(D)* bn;
// We do not explicitly compute inv(D) instead Du = C is solved
// Extract the submatrix
const std::vector<M>& Jn = eqs[n].derivative();
// Use sparse LU to solve the block submatrices i.e compute inv(D)
#if HAVE_UMFPACK
const Eigen::UmfPackLU< M > solver(Jn[n]);
#else
const Eigen::SparseLU< M > solver(Jn[n]);
#endif
M id(Jn[n].rows(), Jn[n].cols());
id.setIdentity();
const Eigen::SparseMatrix<M::Scalar, Eigen::ColMajor> Di = solver.solve(id);
// compute inv(D)*bn for the update of the right hand side
const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());
std::vector<V> vals(num_eq); // Number n will remain empty.
std::vector<std::vector<M>> jacs(num_eq); // Number n will remain empty.
for (int eq = 0; eq < num_eq; ++eq) {
jacs[eq].reserve(num_eq - 1);
const std::vector<M>& Je = eqs[eq].derivative();
const M& B = Je[n];
// Update right hand side.
vals[eq] = eqs[eq].value().matrix() - B * Dibn;
}
for (int var = 0; var < num_eq; ++var) {
if (var == n) {
continue;
}
// solve Du = C
// const M u = Di * Jn[var]; // solver.solve(Jn[var]);
M u;
fastSparseProduct(Di, Jn[var], u); // solver.solve(Jn[var]);
for (int eq = 0; eq < num_eq; ++eq) {
if (eq == n) {
continue;
}
const std::vector<M>& Je = eqs[eq].derivative();
const M& B = Je[n];
// Create new jacobians.
// Add A
jacs[eq].push_back(Je[var]);
M& J = jacs[eq].back();
// Subtract Bu (B*inv(D)*C)
M Bu;
fastSparseProduct(B, u, Bu);
J -= Bu;
}
}
// Create return value.
std::vector<ADB> retval;
retval.reserve(num_eq - 1);
for (int eq = 0; eq < num_eq; ++eq) {
if (eq == n) {
continue;
}
retval.push_back(ADB::function(std::move(vals[eq]), std::move(jacs[eq])));
}
return retval;
}
V recoverVariable(const ADB& equation, const V& partial_solution, const int n)
{
// The equation to solve for the unknown y (to be recovered) is
// Cx + Dy = b
// Dy = (b - Cx)
// where D is the eliminated block, C is the jacobian of
// the eliminated equation with respect to the
// non-eliminated unknowms, b is the right-hand side of
// the eliminated equation, and x is the partial solution
// of the non-eliminated unknowns.
const M& D = equation.derivative()[n];
// Build C.
std::vector<M> C_jacs = equation.derivative();
C_jacs.erase(C_jacs.begin() + n);
V equation_value = equation.value();
ADB eq_coll = collapseJacs(ADB::function(std::move(equation_value), std::move(C_jacs)));
const M& C = eq_coll.derivative()[0];
// Use sparse LU to solve the block submatrices
#if HAVE_UMFPACK
const Eigen::UmfPackLU< M > solver(D);
#else
const Eigen::SparseLU< M > solver(D);
#endif
// Compute value of eliminated variable.
const Eigen::VectorXd b = (equation.value().matrix() - C * partial_solution.matrix());
const Eigen::VectorXd elim_var = solver.solve(b);
// Find the relevant sizes to use when reconstructing the full solution.
const int nelim = equation.size();
const int npart = partial_solution.size();
assert(C.cols() == npart);
const int full_size = nelim + npart;
int start = 0;
for (int i = 0; i < n; ++i) {
start += equation.derivative()[i].cols();
}
assert(start < full_size);
// Reconstruct complete solution vector.
V sol(full_size);
std::copy_n(partial_solution.data(), start, sol.data());
std::copy_n(elim_var.data(), nelim, sol.data() + start);
std::copy_n(partial_solution.data() + start, npart - start, sol.data() + start + nelim);
return sol;
}
/// Form an elliptic system of equations.
/// \param[in] num_phases the number of fluid phases
/// \param[in] eqs the equations
/// \param[out] A the resulting full system matrix
/// \param[out] b the right hand side
/// This function will deal with the first num_phases
/// equations in eqs, and return a matrix A for the full
/// system that has a elliptic upper left corner, if possible.
void formEllipticSystem(const int num_phases,
const std::vector<ADB>& eqs_in,
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
// M& A,
V& b)
{
if (num_phases != 3) {
OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases.");
}
// A concession to MRST, to obtain more similar behaviour:
// swap the first two equations, so that oil is first, then water.
auto eqs = eqs_in;
eqs[0].swap(eqs[1]);
// Characterize the material balance equations.
const int n = eqs[0].size();
const double ratio_limit = 0.01;
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> Block;
// The l1 block indicates if the equation for a given cell and phase is
// sufficiently strong on the diagonal.
Block l1 = Block::Zero(n, num_phases);
for (int phase = 0; phase < num_phases; ++phase) {
const M& J = eqs[phase].derivative()[0];
V dj = J.diagonal().cwiseAbs();
V sod = V::Zero(n);
for (int elem = 0; elem < n; ++elem) {
sod(elem) = J.col(elem).cwiseAbs().sum() - dj(elem);
}
l1.col(phase) = (dj/sod > ratio_limit).cast<double>();
}
// By default, replace first equation with sum of all phase equations.
// Build helper vectors.
V l21 = V::Zero(n);
V l22 = V::Ones(n);
V l31 = V::Zero(n);
V l33 = V::Ones(n);
// If the first phase diagonal is not strong enough, we need further treatment.
// Then the first equation will be the sum of the remaining equations,
// and we swap the first equation into one of their slots.
for (int elem = 0; elem < n; ++elem) {
if (!l1(elem, 0)) {
const double l12x = l1(elem, 1);
const double l13x = l1(elem, 2);
const bool allzero = (l12x + l13x == 0);
if (allzero) {
l1(elem, 0) = 1;
} else {
if (l12x >= l13x) {
l21(elem) = 1;
l22(elem) = 0;
} else {
l31(elem) = 1;
l33(elem) = 0;
}
}
}
}
// Construct the sparse matrix L that does the swaps and sums.
Span i1(n, 1, 0);
Span i2(n, 1, n);
Span i3(n, 1, 2*n);
std::vector< Eigen::Triplet<double> > t;
t.reserve(7*n);
for (int ii = 0; ii < n; ++ii) {
t.emplace_back(i1[ii], i1[ii], l1(ii));
t.emplace_back(i1[ii], i2[ii], l1(ii+n));
t.emplace_back(i1[ii], i3[ii], l1(ii+2*n));
t.emplace_back(i2[ii], i1[ii], l21(ii));
t.emplace_back(i2[ii], i2[ii], l22(ii));
t.emplace_back(i3[ii], i1[ii], l31(ii));
t.emplace_back(i3[ii], i3[ii], l33(ii));
}
M L(3*n, 3*n);
L.setFromTriplets(t.begin(), t.end());
// Combine in single block.
ADB total_residual = vertcatCollapseJacs(eqs);
// Create output as product of L with equations.
A = L * total_residual.derivative()[0];
b = L * total_residual.value().matrix();
}
} // anonymous namespace
} // namespace Opm

246
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp Normal file
View File

@ -0,0 +1,246 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2015 Statoil AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#if HAVE_UMFPACK
#include <Eigen/UmfPackSupport>
#else
#include <Eigen/SparseLU>
#endif
namespace Opm
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
/// Construct a system solver.
NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param,
const boost::any& parallelInformation)
: iterations_( 0 ),
parallelInformation_(parallelInformation),
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
linear_solver_restart_( param.getDefault("linear_solver_restart", 40 ) ),
linear_solver_verbosity_( param.getDefault("linear_solver_verbosity", 0 ))
{
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] residual residual object containing A and b.
/// \return the solution x
NewtonIterationBlackoilInterleaved::SolutionVector
NewtonIterationBlackoilInterleaved::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
{
// Build the vector of equations.
const int np = residual.material_balance_eq.size();
std::vector<ADB> eqs;
eqs.reserve(np + 2);
for (int phase = 0; phase < np; ++phase) {
eqs.push_back(residual.material_balance_eq[phase]);
}
// check if wells are present
const bool hasWells = residual.well_flux_eq.size() > 0 ;
std::vector<ADB> elim_eqs;
if( hasWells )
{
eqs.push_back(residual.well_flux_eq);
eqs.push_back(residual.well_eq);
// Eliminate the well-related unknowns, and corresponding equations.
elim_eqs.reserve(2);
elim_eqs.push_back(eqs[np]);
eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns.
elim_eqs.push_back(eqs[np]);
eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
assert(int(eqs.size()) == np);
}
// Scale material balance equations.
const double matbalscale[3] = { 1.1169, 1.0031, 0.0031 }; // HACK hardcoded instead of computed.
for (int phase = 0; phase < np; ++phase) {
eqs[phase] = eqs[phase] * matbalscale[phase];
}
// Form modified system.
Eigen::SparseMatrix<double, Eigen::RowMajor> A;
V b;
formEllipticSystem(np, eqs, A, b);
// Create ISTL matrix with interleaved rows and columns (block structured).
Mat istlA;
formInterleavedSystem(eqs, A, istlA);
// Solve reduced system.
SolutionVector dx(SolutionVector::Zero(b.size()));
// Right hand side.
const int size = istlA.N();
Vector istlb(size);
for (int i = 0; i < size; ++i) {
istlb[i][0] = b(i);
istlb[i][1] = b(size + i);
istlb[i][2] = b(2*size + i);
}
// System solution
Vector x(istlA.M());
x = 0.0;
Dune::InverseOperatorResult result;
// Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
Comm istlComm(info.communicator());
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
size, np);
// Construct operator, scalar product and vectors needed.
typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
Operator opA(istlA, istlComm);
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, istlb, istlComm, result);
}
else
#endif
{
// Construct operator, scalar product and vectors needed.
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
Operator opA(istlA);
Dune::Amg::SequentialInformation info;
constructPreconditionerAndSolve(opA, x, istlb, info, result);
}
// store number of iterations
iterations_ = result.iterations;
// Check for failure of linear solver.
if (!result.converged) {
OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver.");
}
// Copy solver output to dx.
for (int i = 0; i < size; ++i) {
dx(i) = x[i][0];
dx(size + i) = x[i][1];
dx(2*size + i) = x[i][2];
}
if ( hasWells ) {
// Compute full solution using the eliminated equations.
// Recovery in inverse order of elimination.
dx = recoverVariable(elim_eqs[1], dx, np);
dx = recoverVariable(elim_eqs[0], dx, np);
}
return dx;
}
void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs,
const Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
Mat& istlA) const
{
const int np = eqs.size();
// Find sparsity structure as union of basic block sparsity structures,
// corresponding to the jacobians with respect to pressure.
// Use addition to get to the union structure.
Eigen::SparseMatrix<double> structure = eqs[0].derivative()[0];
for (int phase = 0; phase < np; ++phase) {
structure += eqs[phase].derivative()[0];
}
Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;
// Create ISTL matrix with interleaved rows and columns (block structured).
assert(np == 3);
istlA.setSize(s.rows(), s.cols(), s.nonZeros());
istlA.setBuildMode(Mat::row_wise);
const int* ia = s.outerIndexPtr();
const int* ja = s.innerIndexPtr();
for (Mat::CreateIterator row = istlA.createbegin(); row != istlA.createend(); ++row) {
int ri = row.index();
for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
row.insert(ja[i]);
}
}
const int size = s.rows();
Span span[3] = { Span(size, 1, 0),
Span(size, 1, size),
Span(size, 1, 2*size) };
for (int row = 0; row < size; ++row) {
for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
const int col = ja[col_ix];
MatrixBlockType block;
for (int p1 = 0; p1 < np; ++p1) {
for (int p2 = 0; p2 < np; ++p2) {
block[p1][p2] = A.coeff(span[p1][row], span[p2][col]);
}
}
istlA[row][col] = block;
}
}
}
const boost::any& NewtonIterationBlackoilInterleaved::parallelInformation() const
{
return parallelInformation_;
}
} // namespace Opm

175
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp Normal file
View File

@ -0,0 +1,175 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
Copyright 2015 IRIS AS
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2015 Statoil AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
#define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/utility/platform_dependent/disable_warnings.h>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
#include <memory>
namespace Opm
{
/// This class solves the fully implicit black-oil system by
/// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables).
class NewtonIterationBlackoilInterleaved : public NewtonIterationBlackoilInterface
{
typedef Dune::FieldVector<double, 3 > VectorBlockType;
typedef Dune::FieldMatrix<double, 3, 3> MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector;
public:
/// Construct a system solver.
/// \param[in] param parameters controlling the behaviour of
/// the preconditioning and choice of
/// linear solvers.
/// Parameters:
/// cpr_relax (default 1.0) relaxation for the preconditioner
/// cpr_ilu_n (default 0) use ILU(n) for preconditioning of the linear system
/// cpr_use_amg (default false) if true, use AMG preconditioner for elliptic part
/// cpr_use_bicgstab (default true) if true, use BiCGStab (else use CG) for elliptic part
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param,
const boost::any& parallelInformation=boost::any());
/// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] residual residual object containing A and b.
/// \return the solution x
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
/// \copydoc NewtonIterationBlackoilInterface::iterations
virtual int iterations () const { return iterations_; }
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
virtual const boost::any& parallelInformation() const;
private:
/// \brief construct the CPR preconditioner and the solver.
/// \tparam P The type of the parallel information.
/// \param parallelInformation the information about the parallelization.
template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
void constructPreconditionerAndSolve(O& opA,
Vector& x, Vector& istlb,
const POrComm& parallelInformation,
Dune::InverseOperatorResult& result) const
{
// Construct scalar product.
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
auto sp = ScalarProductChooser::construct(parallelInformation);
// Construct preconditioner.
auto precond = constructPrecond(opA, parallelInformation);
// Communicate if parallel.
parallelInformation.copyOwnerToAll(istlb, istlb);
// Solve.
solve(opA, x, istlb, *sp, precond, result);
}
typedef Dune::SeqILU0<Mat, Vector, Vector> SeqPreconditioner;
template <class Operator>
SeqPreconditioner constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
{
const double relax = 1.0;
SeqPreconditioner precond(opA.getmat(), relax);
return precond;
}
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
typedef Dune::BlockPreconditioner<Vector, Vector, Comm, SeqPreconditioner> ParPreconditioner;
template <class Operator>
ParPreconditioner constructPrecond(Operator& opA, const Comm& comm) const
{
const double relax = 1.0;
SeqPreconditioner seq_precond(opA.getmat(), relax);
ParPreconditioner precond(seq_precond, comm);
return precond;
}
#endif
/// \brief Solve the system using the given preconditioner and scalar product.
template <class Operator, class ScalarProd, class Precond>
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
{
// TODO: Revise when linear solvers interface opm-core is done
// Construct linear solver.
// GMRes solver
if ( newton_use_gmres_ ) {
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
else { // BiCGstab solver
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
}
void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs,
const Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
Mat& istlA) const;
mutable int iterations_;
boost::any parallelInformation_;
const bool newton_use_gmres_;
const double linear_solver_reduction_;
const int linear_solver_maxiter_;
const int linear_solver_restart_;
const int linear_solver_verbosity_;
};
} // namespace Opm
#endif // OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED

275
opm/autodiff/NewtonIterationUtilities.cpp Normal file
View File

@ -0,0 +1,275 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#if HAVE_UMFPACK
#include <Eigen/UmfPackSupport>
#else
#include <Eigen/SparseLU>
#endif
namespace Opm
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
std::vector<ADB> eliminateVariable(const std::vector<ADB>& eqs, const int n)
{
// Check that the variable index to eliminate is within bounds.
const int num_eq = eqs.size();
const int num_vars = eqs[0].derivative().size();
if (num_eq != num_vars) {
OPM_THROW(std::logic_error, "eliminateVariable() requires the same number of variables and equations.");
}
if (n >= num_eq) {
OPM_THROW(std::logic_error, "Trying to eliminate variable from too small set of equations.");
}
// Schur complement of (A B ; C D) wrt. D is A - B*inv(D)*C.
// This is applied to all 2x2 block submatrices
// The right hand side is modified accordingly. bi = bi - B * inv(D)* bn;
// We do not explicitly compute inv(D) instead Du = C is solved
// Extract the submatrix
const std::vector<M>& Jn = eqs[n].derivative();
// Use sparse LU to solve the block submatrices i.e compute inv(D)
#if HAVE_UMFPACK
const Eigen::UmfPackLU< M > solver(Jn[n]);
#else
const Eigen::SparseLU< M > solver(Jn[n]);
#endif
M id(Jn[n].rows(), Jn[n].cols());
id.setIdentity();
const Eigen::SparseMatrix<M::Scalar, Eigen::ColMajor> Di = solver.solve(id);
// compute inv(D)*bn for the update of the right hand side
const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());
std::vector<V> vals(num_eq); // Number n will remain empty.
std::vector<std::vector<M>> jacs(num_eq); // Number n will remain empty.
for (int eq = 0; eq < num_eq; ++eq) {
jacs[eq].reserve(num_eq - 1);
const std::vector<M>& Je = eqs[eq].derivative();
const M& B = Je[n];
// Update right hand side.
vals[eq] = eqs[eq].value().matrix() - B * Dibn;
}
for (int var = 0; var < num_eq; ++var) {
if (var == n) {
continue;
}
// solve Du = C
// const M u = Di * Jn[var]; // solver.solve(Jn[var]);
M u;
fastSparseProduct(Di, Jn[var], u); // solver.solve(Jn[var]);
for (int eq = 0; eq < num_eq; ++eq) {
if (eq == n) {
continue;
}
const std::vector<M>& Je = eqs[eq].derivative();
const M& B = Je[n];
// Create new jacobians.
// Add A
jacs[eq].push_back(Je[var]);
M& J = jacs[eq].back();
// Subtract Bu (B*inv(D)*C)
M Bu;
fastSparseProduct(B, u, Bu);
J -= Bu;
}
}
// Create return value.
std::vector<ADB> retval;
retval.reserve(num_eq - 1);
for (int eq = 0; eq < num_eq; ++eq) {
if (eq == n) {
continue;
}
retval.push_back(ADB::function(std::move(vals[eq]), std::move(jacs[eq])));
}
return retval;
}
V recoverVariable(const ADB& equation, const V& partial_solution, const int n)
{
// The equation to solve for the unknown y (to be recovered) is
// Cx + Dy = b
// Dy = (b - Cx)
// where D is the eliminated block, C is the jacobian of
// the eliminated equation with respect to the
// non-eliminated unknowms, b is the right-hand side of
// the eliminated equation, and x is the partial solution
// of the non-eliminated unknowns.
const M& D = equation.derivative()[n];
// Build C.
std::vector<M> C_jacs = equation.derivative();
C_jacs.erase(C_jacs.begin() + n);
V equation_value = equation.value();
ADB eq_coll = collapseJacs(ADB::function(std::move(equation_value), std::move(C_jacs)));
const M& C = eq_coll.derivative()[0];
// Use sparse LU to solve the block submatrices
#if HAVE_UMFPACK
const Eigen::UmfPackLU< M > solver(D);
#else
const Eigen::SparseLU< M > solver(D);
#endif
// Compute value of eliminated variable.
const Eigen::VectorXd b = (equation.value().matrix() - C * partial_solution.matrix());
const Eigen::VectorXd elim_var = solver.solve(b);
// Find the relevant sizes to use when reconstructing the full solution.
const int nelim = equation.size();
const int npart = partial_solution.size();
assert(C.cols() == npart);
const int full_size = nelim + npart;
int start = 0;
for (int i = 0; i < n; ++i) {
start += equation.derivative()[i].cols();
}
assert(start < full_size);
// Reconstruct complete solution vector.
V sol(full_size);
std::copy_n(partial_solution.data(), start, sol.data());
std::copy_n(elim_var.data(), nelim, sol.data() + start);
std::copy_n(partial_solution.data() + start, npart - start, sol.data() + start + nelim);
return sol;
}
/// Form an elliptic system of equations.
/// \param[in] num_phases the number of fluid phases
/// \param[in] eqs the equations
/// \param[out] A the resulting full system matrix
/// \param[out] b the right hand side
/// This function will deal with the first num_phases
/// equations in eqs, and return a matrix A for the full
/// system that has a elliptic upper left corner, if possible.
void formEllipticSystem(const int num_phases,
const std::vector<ADB>& eqs_in,
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
V& b)
{
if (num_phases != 3) {
OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases.");
}
// A concession to MRST, to obtain more similar behaviour:
// swap the first two equations, so that oil is first, then water.
auto eqs = eqs_in;
eqs[0].swap(eqs[1]);
// Characterize the material balance equations.
const int n = eqs[0].size();
const double ratio_limit = 0.01;
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> Block;
// The l1 block indicates if the equation for a given cell and phase is
// sufficiently strong on the diagonal.
Block l1 = Block::Zero(n, num_phases);
for (int phase = 0; phase < num_phases; ++phase) {
const M& J = eqs[phase].derivative()[0];
V dj = J.diagonal().cwiseAbs();
V sod = V::Zero(n);
for (int elem = 0; elem < n; ++elem) {
sod(elem) = J.col(elem).cwiseAbs().sum() - dj(elem);
}
l1.col(phase) = (dj/sod > ratio_limit).cast<double>();
}
// By default, replace first equation with sum of all phase equations.
// Build helper vectors.
V l21 = V::Zero(n);
V l22 = V::Ones(n);
V l31 = V::Zero(n);
V l33 = V::Ones(n);
// If the first phase diagonal is not strong enough, we need further treatment.
// Then the first equation will be the sum of the remaining equations,
// and we swap the first equation into one of their slots.
for (int elem = 0; elem < n; ++elem) {
if (!l1(elem, 0)) {
const double l12x = l1(elem, 1);
const double l13x = l1(elem, 2);
const bool allzero = (l12x + l13x == 0);
if (allzero) {
l1(elem, 0) = 1;
} else {
if (l12x >= l13x) {
l21(elem) = 1;
l22(elem) = 0;
} else {
l31(elem) = 1;
l33(elem) = 0;
}
}
}
}
// Construct the sparse matrix L that does the swaps and sums.
Span i1(n, 1, 0);
Span i2(n, 1, n);
Span i3(n, 1, 2*n);
std::vector< Eigen::Triplet<double> > t;
t.reserve(7*n);
for (int ii = 0; ii < n; ++ii) {
t.emplace_back(i1[ii], i1[ii], l1(ii));
t.emplace_back(i1[ii], i2[ii], l1(ii+n));
t.emplace_back(i1[ii], i3[ii], l1(ii+2*n));
t.emplace_back(i2[ii], i1[ii], l21(ii));
t.emplace_back(i2[ii], i2[ii], l22(ii));
t.emplace_back(i3[ii], i1[ii], l31(ii));
t.emplace_back(i3[ii], i3[ii], l33(ii));
}
M L(3*n, 3*n);
L.setFromTriplets(t.begin(), t.end());
// Combine in single block.
ADB total_residual = vertcatCollapseJacs(eqs);
// Create output as product of L with equations.
A = L * total_residual.derivative()[0];
b = L * total_residual.value().matrix();
}
} // namespace Opm

64
opm/autodiff/NewtonIterationUtilities.hpp Normal file
View File

@ -0,0 +1,64 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED
#define OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <vector>
namespace Opm
{
/// Eliminate a variable via Schur complement.
/// \param[in] eqs set of equations with Jacobians
/// \param[in] n index of equation/variable to eliminate.
/// \return new set of equations, one smaller than eqs.
/// Note: this method requires the eliminated variable to have the same size
/// as the equation in the corresponding position (that also will be eliminated).
std::vector< AutoDiffBlock<double> >
eliminateVariable(const std::vector< AutoDiffBlock<double> >& eqs,
const int n);
/// Recover that value of a variable previously eliminated.
/// \param[in] equation previously eliminated equation.
/// \param[in] partial_solution solution to the remainder system after elimination.
/// \param[in] n index of equation/variable that was eliminated.
/// \return solution to complete system.
AutoDiffBlock<double>::V recoverVariable(const AutoDiffBlock<double>& equation,
const AutoDiffBlock<double>::V& partial_solution,
const int n);
/// Form an elliptic system of equations.
/// \param[in] num_phases the number of fluid phases
/// \param[in] eqs the equations
/// \param[out] A the resulting full system matrix
/// \param[out] b the right hand side
/// This function will deal with the first num_phases
/// equations in eqs, and return a matrix A for the full
/// system that has a elliptic upper left corner, if possible.
void formEllipticSystem(const int num_phases,
const std::vector< AutoDiffBlock<double> >& eqs,
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
AutoDiffBlock<double>::V& b);
} // namespace Opm
#endif // OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED