diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 6f131a692..eb38e90cd 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/examples/flow.cpp b/examples/flow.cpp index 606149c37..58f1cacbf 100644 --- a/examples/flow.cpp +++ b/examples/flow.cpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -238,7 +239,9 @@ try // Solver for Newton iterations. std::unique_ptr fis_solver; - if (param.getDefault("use_cpr", true)) { + 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)); diff --git a/examples/flow_cp.cpp b/examples/flow_cp.cpp index 370ea1828..225bc1d83 100644 --- a/examples/flow_cp.cpp +++ b/examples/flow_cp.cpp @@ -61,6 +61,7 @@ #include #include #include +#include #include #include @@ -338,7 +339,9 @@ try boost::any parallel_information; Opm::extractParallelGridInformationToISTL(*grid, parallel_information); - if (param.getDefault("use_cpr", true)) { + if (param.getDefault("use_interleaved", false)) { + fis_solver.reset(new NewtonIterationBlackoilInterleaved(param, parallel_information)); + } else if (param.getDefault("use_cpr", true)) { fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information)); } else { fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information)); diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index ba920acac..10664e107 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -25,25 +25,11 @@ #include #include +#include #include -#include -#include #include -#include +#include #include -#include - -// #include -#include -#include -#include -#include -#include -#include -#include -#include - -#include #if HAVE_UMFPACK #include @@ -51,6 +37,7 @@ #include #endif + namespace Opm { @@ -58,41 +45,6 @@ namespace Opm typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; - typedef Dune::FieldMatrix MatrixBlockType; - typedef Dune::BCRSMatrix 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 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); - - /// 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); - - } // 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 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 - // 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& 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 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 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) { - jacs[eq].reserve(num_eq - 1); - const std::vector& 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& 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 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 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& eqs_in, - Eigen::SparseMatrix& 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 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 = 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 diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp new file mode 100644 index 000000000..d8851bc26 --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -0,0 +1,229 @@ +/* + 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 . +*/ + +#include + +#include + +#include +#include +#include +#include +#include + +#if HAVE_UMFPACK +#include +#else +#include +#endif + + +namespace Opm +{ + + + typedef AutoDiffBlock 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 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 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]; + } + + // 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 structure = eqs[0].derivative()[0]; + for (int phase = 0; phase < np; ++phase) { + structure += eqs[phase].derivative()[0]; + } + Eigen::SparseMatrix s = structure; + + // Form modified system. + Eigen::SparseMatrix A; + V b; + formEllipticSystem(np, eqs, A, b); + + // Create ISTL matrix with interleaved rows and columns (block structured). + assert(np == 3); + Mat istlA(s.rows(), s.cols(), s.nonZeros(), 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; + } + } + + // Solve reduced system. + SolutionVector dx(SolutionVector::Zero(b.size())); + + // Right hand side. + Vector istlb(istlA.N()); + 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 Comm; + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation_); + Comm istlComm(info.communicator()); + info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(), + size, np); + // Construct operator, scalar product and vectors needed. + typedef Dune::OverlappingSchwarzOperator Operator; + Operator opA(istlA, istlComm); + constructPreconditionerAndSolve(opA, x, istlb, istlComm, result); + } + else +#endif + { + // Construct operator, scalar product and vectors needed. + typedef Dune::MatrixAdapter 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; + } + + + + + + const boost::any& NewtonIterationBlackoilInterleaved::parallelInformation() const + { + return parallelInformation_; + } + + + + + +} // namespace Opm + diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp new file mode 100644 index 000000000..f6aacfb69 --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -0,0 +1,171 @@ +/* + 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 . +*/ + +#ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED +#define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED + + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include + +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 VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector 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 + void constructPreconditionerAndSolve(O& opA, + Vector& x, Vector& istlb, + const POrComm& parallelInformation, + Dune::InverseOperatorResult& result) const + { + // Construct scalar product. + typedef Dune::ScalarProductChooser 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 SeqPreconditioner; + typedef Dune::OwnerOverlapCopyCommunication Comm; + typedef Dune::BlockPreconditioner ParPreconditioner; + + + template + SeqPreconditioner constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const + { + const double relax = 1.0; + SeqPreconditioner precond(opA.getmat(), relax); + return precond; + } + + + template + 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; + } + + + /// \brief Solve the system using the given preconditioner and scalar product. + template + 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 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 linsolve(opA, sp, precond, + linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_); + // Solve system. + linsolve.apply(x, istlb, result); + } + } + + 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 diff --git a/opm/autodiff/NewtonIterationUtilities.cpp b/opm/autodiff/NewtonIterationUtilities.cpp new file mode 100644 index 000000000..aa78eaee6 --- /dev/null +++ b/opm/autodiff/NewtonIterationUtilities.cpp @@ -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 . +*/ + +#include + +#include +#include +#include + +#if HAVE_UMFPACK +#include +#else +#include +#endif + +namespace Opm +{ + + + typedef AutoDiffBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + + + 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 + // 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& 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 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 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) { + jacs[eq].reserve(num_eq - 1); + const std::vector& 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& 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 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 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& 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; + 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 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 = vertcatCollapseJacs(eqs); + + // Create output as product of L with equations. + A = L * total_residual.derivative()[0]; + b = L * total_residual.value().matrix(); + } + + + +} // namespace Opm + diff --git a/opm/autodiff/NewtonIterationUtilities.hpp b/opm/autodiff/NewtonIterationUtilities.hpp new file mode 100644 index 000000000..1e383d477 --- /dev/null +++ b/opm/autodiff/NewtonIterationUtilities.hpp @@ -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 . +*/ + +#ifndef OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED +#define OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED + +#include +#include + +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 > + eliminateVariable(const std::vector< AutoDiffBlock >& 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::V recoverVariable(const AutoDiffBlock& equation, + const AutoDiffBlock::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 >& eqs, + Eigen::SparseMatrix& A, + AutoDiffBlock::V& b); + + +} // namespace Opm + +#endif // OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED