From 5d0f654443268c62e31162f734c8692d28246c83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 1 Jun 2015 13:33:37 +0200 Subject: [PATCH] Add class NewtonIterationBlackoilInterleaved. Initially it is just a copy of the NewtonIterationBlackoilCPR class. Also, add use_interleaved parameter to use the class. --- CMakeLists_files.cmake | 2 + examples/flow.cpp | 5 +- .../NewtonIterationBlackoilInterleaved.cpp | 497 ++++++++++++++++++ .../NewtonIterationBlackoilInterleaved.hpp | 152 ++++++ 4 files changed, 655 insertions(+), 1 deletion(-) create mode 100644 opm/autodiff/NewtonIterationBlackoilInterleaved.cpp create mode 100644 opm/autodiff/NewtonIterationBlackoilInterleaved.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 6f131a692..f3f6b0259 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -30,6 +30,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/ExtractParallelGridInformationToISTL.cpp opm/autodiff/NewtonIterationBlackoilCPR.cpp opm/autodiff/NewtonIterationBlackoilSimple.cpp + opm/autodiff/NewtonIterationBlackoilInterleaved.cpp opm/autodiff/GridHelpers.cpp opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -111,6 +112,7 @@ 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/NewtonSolver.hpp opm/autodiff/NewtonSolver_impl.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/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp new file mode 100644 index 000000000..62b33f56a --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -0,0 +1,497 @@ +/* + 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 . +*/ + +/* + Copyright 2014 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 +#include +#include +#include + +// #include +#include +#include +#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; + 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 + + + + + + /// Construct a system solver. + NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param, + const boost::any& parallelInformation) + : cpr_param_( param ), + 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]; + } + + // 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. + DuneMatrix istlA( A ); + + // Create ISTL matrix for elliptic part. + DuneMatrix istlAe( A.topLeftCorner(nc, nc) ); + + // Right hand side. + Vector istlb(istlA.N()); + std::copy_n(b.data(), istlb.size(), istlb.begin()); + // System solution + Vector x(istlA.M()); + x = 0.0; + + Dune::InverseOperatorResult result; +#if HAVE_MPI + if(parallelInformation_.type()==typeid(ParallelISTLInformation)) + { + typedef Dune::OwnerOverlapCopyCommunication Comm; + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation_); + Comm istlComm(info.communicator()); + Comm istlAeComm(info.communicator()); + info.copyValuesTo(istlAeComm.indexSet(), istlAeComm.remoteIndices()); + info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(), + istlAe.N(), istlA.N()/istlAe.N()); + // Construct operator, scalar product and vectors needed. + typedef Dune::OverlappingSchwarzOperator Operator; + Operator opA(istlA, istlComm); + constructPreconditionerAndSolve(opA, istlAe, x, istlb, istlComm, istlAeComm, result); + } + else +#endif + { + // Construct operator, scalar product and vectors needed. + typedef Dune::MatrixAdapter Operator; + Operator opA(istlA); + Dune::Amg::SequentialInformation info; + constructPreconditionerAndSolve(opA, istlAe, x, istlb, info, 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. + std::copy(x.begin(), x.end(), dx.data()); + + 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 + { + + + 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.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp new file mode 100644 index 000000000..fb609f9f4 --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -0,0 +1,152 @@ +/* + 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_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED +#define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED + +/* + Copyright 2014 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 . +*/ + +#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, DuneMatrix& istlAe, + Vector& x, Vector& istlb, + const P& parallelInformation, + const P& parallelInformationAe, + Dune::InverseOperatorResult& result) const + { + typedef Dune::ScalarProductChooser ScalarProductChooser; + std::unique_ptr + sp(ScalarProductChooser::construct(parallelInformation)); + // Construct preconditioner. + // typedef Dune::SeqILU0 Preconditioner; + typedef Opm::CPRPreconditioner Preconditioner; + parallelInformation.copyOwnerToAll(istlb, istlb); + Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation, + parallelInformationAe); + + // 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); + } + } + + CPRParameter cpr_param_; + + 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