diff --git a/CMakeLists.txt b/CMakeLists.txt index 78fd2da69..32bd28881 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -75,7 +75,7 @@ if (NOT EIGEN3_FOUND) include (ExternalProject) externalProject_Add (Eigen3 GIT_REPOSITORY git://github.com/OPM/eigen3 - UPDATE_COMMAND git checkout 9f6cc779c101b87184076322603f496e5fdd0432 + UPDATE_COMMAND git checkout 9e788db99d73f3199ade74bfda8d9f73fdfbbe4c CMAKE_ARGS -DEIGEN_TEST_NO_OPENGL=1 -DEIGEN_BUILD_PKGCONFIG=0 -DCMAKE_INSTALL_PREFIX=${CMAKE_BINARY_DIR}/eigen3-installed ) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 879cedefb..e9383485e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -28,6 +28,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/BlackoilPropsAd.cpp opm/autodiff/BlackoilPropsAdInterface.cpp + opm/autodiff/NewtonIterationBlackoilCPR.cpp opm/autodiff/NewtonIterationBlackoilSimple.cpp opm/autodiff/GridHelpers.cpp opm/autodiff/ImpesTPFAAD.cpp @@ -93,11 +94,13 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilPropsAd.hpp opm/autodiff/BlackoilPropsAdFromDeck.hpp opm/autodiff/BlackoilPropsAdInterface.hpp + opm/autodiff/CPRPreconditioner.hpp opm/autodiff/GeoProps.hpp opm/autodiff/GridHelpers.hpp opm/autodiff/ImpesTPFAAD.hpp opm/autodiff/FullyImplicitBlackoilSolver.hpp opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp + opm/autodiff/NewtonIterationBlackoilCPR.hpp opm/autodiff/NewtonIterationBlackoilInterface.hpp opm/autodiff/NewtonIterationBlackoilSimple.hpp opm/autodiff/LinearisedBlackoilResidual.hpp diff --git a/cmake/Modules/opm-autodiff-prereqs.cmake b/cmake/Modules/opm-autodiff-prereqs.cmake index 25972fa70..c08d036b8 100644 --- a/cmake/Modules/opm-autodiff-prereqs.cmake +++ b/cmake/Modules/opm-autodiff-prereqs.cmake @@ -21,5 +21,5 @@ set (opm-autodiff_DEPS dune-cornerpoint; opm-core REQUIRED" # Eigen - "Eigen3 3.1.2 " + "Eigen3 3.2.0 " ) diff --git a/examples/sim_2p_comp_ad.cpp b/examples/sim_2p_comp_ad.cpp index fd6183d0b..842517426 100644 --- a/examples/sim_2p_comp_ad.cpp +++ b/examples/sim_2p_comp_ad.cpp @@ -104,7 +104,7 @@ try if (use_deck) { std::string deck_filename = param.get("deck_filename"); Opm::ParserPtr parser(new Opm::Parser()); - Opm::DeckConstPtr deck = parser->parseFile( deck_filename ); + deck = parser->parseFile( deck_filename ); eclipseState.reset(new EclipseState(deck)); // Grid init diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index c275ef1a4..b8e9ca09a 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -38,6 +38,7 @@ #include #include +#include #include #include @@ -141,9 +142,13 @@ try bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); const double *grav = use_gravity ? &gravity[0] : 0; - // Linear solver. - LinearSolverFactory linsolver(param); - NewtonIterationBlackoilSimple fis_solver(linsolver); + // Solver for Newton iterations. + std::unique_ptr fis_solver; + if (param.getDefault("use_cpr", false)) { + fis_solver.reset(new NewtonIterationBlackoilCPR(param)); + } else { + fis_solver.reset(new NewtonIterationBlackoilSimple(param)); + } // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -212,7 +217,7 @@ try *new_props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - fis_solver, + *fis_solver, grav); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); diff --git a/examples/sim_fibo_ad_cp.cpp b/examples/sim_fibo_ad_cp.cpp index ecb4d0ace..300780df1 100644 --- a/examples/sim_fibo_ad_cp.cpp +++ b/examples/sim_fibo_ad_cp.cpp @@ -60,6 +60,7 @@ #include #include +#include #include #include @@ -194,9 +195,13 @@ try bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); const double *grav = use_gravity ? &gravity[0] : 0; - // Linear solver. - LinearSolverFactory linsolver(param); - NewtonIterationBlackoilSimple fis_solver(linsolver); + // Solver for Newton iterations. + std::unique_ptr fis_solver; + if (param.getDefault("use_cpr", false)) { + fis_solver.reset(new NewtonIterationBlackoilCPR(param)); + } else { + fis_solver.reset(new NewtonIterationBlackoilSimple(param)); + } // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -271,7 +276,7 @@ try *new_props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - fis_solver, + *fis_solver, grav); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 737793b92..842e981dd 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -102,8 +102,7 @@ try double grav[] = { 0.0, 0.0 }; Opm::DerivedGeology geo(*g, props, grav); - Opm::LinearSolverFactory linsolver(param); - Opm::NewtonIterationBlackoilSimple fis_solver(linsolver); + Opm::NewtonIterationBlackoilSimple fis_solver(param); Opm::FullyImplicitBlackoilSolver solver(param, *g, props, geo, 0, *wells, fis_solver); diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp new file mode 100644 index 000000000..327920c72 --- /dev/null +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -0,0 +1,182 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_CPRPRECONDITIONER_HEADER_INCLUDED +#define OPM_CPRPRECONDITIONER_HEADER_INCLUDED + + +#include "disable_warning_pragmas.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "reenable_warning_pragmas.h" + + +namespace Opm +{ + /*! + \brief Sequential CPR preconditioner. + + This is a two-stage preconditioner, combining an elliptic-type + partial solution with ILU0 for the whole system. + + \tparam M The matrix type to operate on + \tparam X Type of the update + \tparam Y Type of the defect + */ + template + class CPRPreconditioner : public Dune::Preconditioner { + public: + //! \brief The matrix type the preconditioner is for. + typedef typename Dune::remove_const::type matrix_type; + //! \brief The domain type of the preconditioner. + typedef X domain_type; + //! \brief The range type of the preconditioner. + typedef Y range_type; + //! \brief The field type of the preconditioner. + typedef typename X::field_type field_type; + + // define the category + enum { + //! \brief The category the preconditioner is part of. + category = Dune::SolverCategory::sequential + }; + + /*! \brief Constructor. + + Constructor gets all parameters to operate the prec. + \param A The matrix to operate on. + \param Ae The top-left elliptic part of A. + \param w The ILU0 relaxation factor. + */ + CPRPreconditioner (const M& A, const M& Ae, const field_type relax) + : A_(A), + ILU_(A), // copy A (will be overwritten by ILU decomp) + Ae_(Ae), + relax_(relax) + { + Dune::bilu0_decomposition(ILU_); + } + + /*! + \brief Prepare the preconditioner. + + \copydoc Preconditioner::pre(X&,Y&) + */ + virtual void pre (X& /*x*/, Y& /*b*/) + { + } + + /*! + \brief Apply the preconditoner. + + \copydoc Preconditioner::apply(X&,const Y&) + */ + virtual void apply (X& v, const Y& d) + { + // Extract part of d corresponding to elliptic part. + Y de(Ae_.N()); + // Note: Assumes that the elliptic part comes first. + std::copy_n(d.begin(), Ae_.N(), de.begin()); + + // Solve elliptic part, extend solution to full. + Y ve = solveElliptic(de); + Y vfull(ILU_.N()); + vfull = 0.0; + // Again assuming that the elliptic part comes first. + std::copy(ve.begin(), ve.end(), vfull.begin()); + + // Subtract elliptic residual from initial residual. + // dmodified = d - A * vfull + Y dmodified = d; + A_.mmv(vfull, dmodified); + + // Apply ILU0. + Y vilu(ILU_.N()); + Dune::bilu_backsolve(ILU_, vilu, dmodified); + v = vfull; + v += vilu; + v *= relax_; + } + + /*! + \brief Clean up. + + \copydoc Preconditioner::post(X&) + */ + virtual void post (X& /*x*/) + { + } + + private: + Y solveElliptic(Y& de) + { + // std::cout << "solveElliptic()" << std::endl; + // Construct operator, scalar product and vectors needed. + typedef Dune::MatrixAdapter Operator; + Operator opAe(Ae_); + Dune::SeqScalarProduct sp; + // Right hand side. + // System solution + X x(opAe.getmat().M()); + x = 0.0; + + // Construct preconditioner. + typedef typename Dune::SeqILU0 Preconditioner; + const double relax = 1.0; + Preconditioner precond(Ae_, relax); + + // Construct linear solver. + const double tolerance = 1e-4; + const int maxit = 5000; + const int verbosity = 0; + Dune::BiCGSTABSolver linsolve(opAe, sp, precond, tolerance, maxit, verbosity); + + // Solve system. + Dune::InverseOperatorResult result; + linsolve.apply(x, de, result); + if (result.converged) { + // std::cout << "solveElliptic() successful!" << std::endl; + } + return x; + } + + //! \brief The matrix for the full linear problem. + const matrix_type& A_; + //! \brief The ILU0 decomposition of the matrix. + matrix_type ILU_; + //! \brief The elliptic part of the matrix. + matrix_type Ae_; + //! \brief The relaxation factor to use. + field_type relax_; + }; + +} // namespace Opm + +#endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 0e93b90d3..89e98d7c4 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -283,9 +283,10 @@ namespace Opm { bool &oscillate, bool &stagnate ) const; void stablizeNewton(V &dx, V &dxOld, const double omega, const RelaxType relax_type) const; - const double dpMaxRel() const { return dp_max_rel_; } - const double dsMax() const { return ds_max_; } - const double drsMaxRel() const { return drs_max_rel_; } + + double dpMaxRel() const { return dp_max_rel_; } + double dsMax() const { return ds_max_; } + double drsMaxRel() const { return drs_max_rel_; } }; } // namespace Opm diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp new file mode 100644 index 000000000..3d0752d53 --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -0,0 +1,488 @@ +/* + Copyright 2014 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 . +*/ + +#include + +#include +#include +#include +#include +#include +#include + +#include "disable_warning_pragmas.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "reenable_warning_pragmas.h" + + +namespace Opm +{ + + + typedef AutoDiffBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector Vector; + + + 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 eliminateVariable(const std::vector& 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); + + /// Determine diagonality of a sparse matrix. + /// If there are off-diagonal elements in the sparse + /// structure, this function returns true if they are all + /// equal to zero. + /// \param[in] matrix the matrix under consideration + /// \return true if matrix is diagonal + bool isDiagonal(const M& matrix); + + /// 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& eqs, + Eigen::SparseMatrix& A, + V& b); + + /// Create a dune-istl matrix from an Eigen matrix. + /// \param[in] matrix input Eigen::SparseMatrix + /// \return output Dune::BCRSMatrix + Mat makeIstlMatrix(const Eigen::SparseMatrix& matrix); + + } // anonymous namespace + + + + + + /// Construct a system solver. + /// \param[in] linsolver linear solver to use + NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& /*param*/) + { + } + + + + + + /// 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 + NewtonIterationBlackoilCPR::SolutionVector + NewtonIterationBlackoilCPR::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const + { + // Build the vector of equations. + const int np = residual.material_balance_eq.size(); + std::vector eqs; + eqs.reserve(np + 2); + for (int phase = 0; phase < np; ++phase) { + eqs.push_back(residual.material_balance_eq[phase]); + } + eqs.push_back(residual.well_flux_eq); + eqs.push_back(residual.well_eq); + + // Eliminate the well-related unknowns, and corresponding equations. + std::vector elim_eqs; + 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]; + } + + // Add material balance equations (or other manipulations) to + // form pressure equation in top left of full system. + Eigen::SparseMatrix A; + V b; + formEllipticSystem(np, eqs, A, b); + + // Scale pressure equation. + const double pscale = 200*unit::barsa; + const int nc = residual.material_balance_eq[0].size(); + A.topRows(nc) *= pscale; + b.topRows(nc) *= pscale; + + // Solve reduced system. + SolutionVector dx(SolutionVector::Zero(b.size())); + + // Create ISTL matrix. + Mat istlA = makeIstlMatrix(A); + + // Create ISTL matrix for elliptic part. + Mat istlAe = makeIstlMatrix(A.topLeftCorner(nc, nc)); + + // Construct operator, scalar product and vectors needed. + typedef Dune::MatrixAdapter Operator; + Operator opA(istlA); + Dune::SeqScalarProduct sp; + // Right hand side. + Vector istlb(opA.getmat().N()); + std::copy_n(b.data(), istlb.size(), istlb.begin()); + // System solution + Vector x(opA.getmat().M()); + x = 0.0; + + // Construct preconditioner. + // typedef Dune::SeqILU0 Preconditioner; + typedef Opm::CPRPreconditioner Preconditioner; + const double relax = 1.0; + Preconditioner precond(istlA, istlAe, relax); + + // Construct linear solver. + const double tolerance = 1e-3; + const int maxit = 5000; + const int verbosity = 1; + const int restart = 40; + Dune::RestartedGMResSolver linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity); + + // Solve system. + Dune::InverseOperatorResult result; + linsolve.apply(x, istlb, result); + + // Check for failure of linear solver. + if (!result.converged) { + OPM_THROW(std::runtime_error, "Convergence failure for linear solver."); + } + + // Copy solver output to dx. + std::copy(x.begin(), x.end(), dx.data()); + + // 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; + } + + + + + namespace + { + + + std::vector eliminateVariable(const std::vector& 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. + // We require that D is diagonal. + const M& D = eqs[n].derivative()[n]; + if (!isDiagonal(D)) { + // std::cout << "++++++++++++++++++++++++++++++++++++++++++++\n" + // << D + // << "++++++++++++++++++++++++++++++++++++++++++++\n" << std::endl; + OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + } + V diag = D.diagonal(); + Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); + std::vector vals(num_eq); // Number n will remain empty. + std::vector> jacs(num_eq); // Number n will remain empty. + for (int eq = 0; eq < num_eq; ++eq) { + if (eq == n) { + continue; + } + jacs[eq].reserve(num_eq - 1); + const M& B = eqs[eq].derivative()[n]; + for (int var = 0; var < num_eq; ++var) { + if (var == n) { + continue; + } + // Create new jacobians. + M schur_jac = eqs[eq].derivative()[var] - B * (invD * eqs[n].derivative()[var]); + jacs[eq].push_back(schur_jac); + } + // Update right hand side. + vals[eq] = eqs[eq].value().matrix() - B * (invD * eqs[n].value().matrix()); + } + + // Create return value. + std::vector retval; + retval.reserve(num_eq - 1); + for (int eq = 0; eq < num_eq; ++eq) { + if (eq == n) { + continue; + } + retval.push_back(ADB::function(vals[eq], 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 + // y = inv(D) (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. + // We require that D is diagonal. + + // Find inv(D). + const M& D = equation.derivative()[n]; + if (!isDiagonal(D)) { + OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + } + V diag = D.diagonal(); + Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); + + // Build C. + std::vector C_jacs = equation.derivative(); + C_jacs.erase(C_jacs.begin() + n); + ADB eq_coll = collapseJacs(ADB::function(equation.value(), C_jacs)); + const M& C = eq_coll.derivative()[0]; + + // Compute value of eliminated variable. + V elim_var = invD * (equation.value().matrix() - C * partial_solution.matrix()); + + // 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; + } + + + + + + bool isDiagonal(const M& matr) + { + M matrix = matr; + matrix.makeCompressed(); + for (int k = 0; k < matrix.outerSize(); ++k) { + for (M::InnerIterator it(matrix, k); it; ++it) { + if (it.col() != it.row()) { + // Off-diagonal element. + if (it.value() != 0.0) { + // Nonzero off-diagonal element. + // std::cout << "off-diag: " << it.row() << ' ' << it.col() << std::endl; + return false; + } + } + } + } + return true; + } + + + + + /// 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& eqs_in, + Eigen::SparseMatrix& 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; + std::swap(eqs[0], eqs[1]); + + // Characterize the material balance equations. + const int n = eqs[0].size(); + const double ratio_limit = 0.01; + typedef Eigen::Array 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(); + } + + // 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 > 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 = eqs[0]; + for (int phase = 1; phase < num_phases; ++phase) { + total_residual = vertcat(total_residual, eqs[phase]); + } + total_residual = collapseJacs(total_residual); + + // Create output as product of L with equations. + A = L * total_residual.derivative()[0]; + b = L * total_residual.value().matrix(); + } + + + + + + Mat makeIstlMatrix(const Eigen::SparseMatrix& matrix) + { + // Create ISTL matrix. + const int size = matrix.rows(); + const int nonzeros = matrix.nonZeros(); + const int* ia = matrix.outerIndexPtr(); + const int* ja = matrix.innerIndexPtr(); + const double* sa = matrix.valuePtr(); + Mat A(size, size, nonzeros, Mat::row_wise); + for (Mat::CreateIterator row = A.createbegin(); row != A.createend(); ++row) { + const int ri = row.index(); + for (int i = ia[ri]; i < ia[ri + 1]; ++i) { + row.insert(ja[i]); + } + } + for (int ri = 0; ri < size; ++ri) { + for (int i = ia[ri]; i < ia[ri + 1]; ++i) { + A[ri][ja[i]] = sa[i]; + } + } + return A; + } + + + + } // anonymous namespace + + +} // namespace Opm + diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp new file mode 100644 index 000000000..c8101ccef --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilCPR.hpp @@ -0,0 +1,58 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED +#define OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED + + +#include +#include +#include +#include + +namespace Opm +{ + + /// This class solves the fully implicit black-oil system by + /// applying a Constrained Pressure Residual preconditioning + /// strategy. + /// The approach is similar to the one described in + /// "Preconditioning for Efficiently Applying Algebraic Multigrid + /// in Fully Implicit Reservoir Simulations" by Gries et al (SPE 163608). + class NewtonIterationBlackoilCPR : public NewtonIterationBlackoilInterface + { + public: + /// Construct a system solver. + /// \param[in] param parameters controlling the behaviour of + /// the preconditioning and choice of + /// linear solvers. + /// Note: parameters currently unused. + NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param); + + /// 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; + }; + +} // namespace Opm + +#endif // OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp index 92e885903..4ced9bcb5 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp +++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp @@ -22,15 +22,16 @@ #include #include #include +#include namespace Opm { /// Construct a system solver. /// \param[in] linsolver linear solver to use - NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver) - : linsolver_(linsolver) + NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param) { + linsolver_.reset(new LinearSolverFactory(param)); } /// Solve the linear system Ax = b, with A being the @@ -54,9 +55,9 @@ namespace Opm SolutionVector dx(SolutionVector::Zero(total_residual.size())); Opm::LinearSolverInterface::LinearSolverReport rep - = linsolver_.solve(matr.rows(), matr.nonZeros(), - matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), - total_residual.value().data(), dx.data()); + = linsolver_->solve(matr.rows(), matr.nonZeros(), + matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), + total_residual.value().data(), dx.data()); if (!rep.converged) { OPM_THROW(std::runtime_error, "FullyImplicitBlackoilSolver::solveJacobianSystem(): " diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.hpp b/opm/autodiff/NewtonIterationBlackoilSimple.hpp index ce17ef8b1..7fceaa28a 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.hpp +++ b/opm/autodiff/NewtonIterationBlackoilSimple.hpp @@ -22,7 +22,9 @@ #include +#include #include +#include namespace Opm { @@ -35,8 +37,9 @@ namespace Opm { public: /// Construct a system solver. - /// \param[in] linsolver linear solver to use - NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver); + /// \param[in] param parameters controlling the behaviour and + /// choice of linear solver. + NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param); /// Solve the system of linear equations Ax = b, with A being the /// combined derivative matrix of the residual and b @@ -46,7 +49,7 @@ namespace Opm virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const; private: - const LinearSolverInterface& linsolver_; + std::unique_ptr linsolver_; }; } // namespace Opm