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 01/13] 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 From 63285bf6f97881d968e9d645a245b971533292be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 10 Jun 2015 12:24:38 +0200 Subject: [PATCH 02/13] Remove usage of CPRPreconditioner. --- .../NewtonIterationBlackoilInterleaved.cpp | 73 +++++++------------ .../NewtonIterationBlackoilInterleaved.hpp | 13 ++-- 2 files changed, 30 insertions(+), 56 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 62b33f56a..dd37214fa 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -98,7 +98,7 @@ namespace Opm /// \return solution to complete system. V recoverVariable(const ADB& equation, const V& partial_solution, const int n); - /// Form an elliptic system of equations. + /// Form an interleaved 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 @@ -106,10 +106,10 @@ namespace Opm /// 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); + void formInterleavedSystem(const int num_phases, + const std::vector& eqs, + Eigen::SparseMatrix& A, + V& b); } // anonymous namespace @@ -120,8 +120,7 @@ namespace Opm /// Construct a system solver. NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param, const boost::any& parallelInformation) - : cpr_param_( param ), - iterations_( 0 ), + : iterations_( 0 ), parallelInformation_(parallelInformation), newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ), linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ), @@ -172,17 +171,16 @@ namespace Opm eqs[phase] = eqs[phase] * matbalscale[phase]; } - // Add material balance equations (or other manipulations) to - // form pressure equation in top left of full system. + // Form interleaved system. Eigen::SparseMatrix A; V b; - formEllipticSystem(np, eqs, A, b); + formInterleavedSystem(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; + // // 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())); @@ -191,7 +189,7 @@ namespace Opm DuneMatrix istlA( A ); // Create ISTL matrix for elliptic part. - DuneMatrix istlAe( A.topLeftCorner(nc, nc) ); + // DuneMatrix istlAe( A.topLeftCorner(nc, nc) ); // Right hand side. Vector istlb(istlA.N()); @@ -201,31 +199,11 @@ namespace Opm 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); - } + // 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; @@ -395,7 +373,7 @@ namespace Opm - /// Form an elliptic system of equations. + /// Form an interleaved 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 @@ -403,14 +381,13 @@ namespace Opm /// 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) + void formInterleavedSystem(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."); + OPM_THROW(std::logic_error, "formInterleavedSystem() requires 3 phases."); } // A concession to MRST, to obtain more similar behaviour: diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index fb609f9f4..0bab86c9a 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -101,21 +101,20 @@ namespace Opm /// \tparam P The type of the parallel information. /// \param parallelInformation the information about the parallelization. template - void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe, + void constructPreconditionerAndSolve(O& opA, 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; + typedef Dune::SeqILU0 Preconditioner; + // typedef Opm::CPRPreconditioner Preconditioner; parallelInformation.copyOwnerToAll(istlb, istlb); - Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation, - parallelInformationAe); + const double relax = 1.0; + Preconditioner precond(opA.getmat(), relax); // TODO: Revise when linear solvers interface opm-core is done // Construct linear solver. @@ -134,8 +133,6 @@ namespace Opm } } - CPRParameter cpr_param_; - mutable int iterations_; boost::any parallelInformation_; From 8af1ea1e4268c01e8907e19d5db094a61645414f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 15 Jun 2015 10:53:07 +0200 Subject: [PATCH 03/13] Add (disabled) experimental part, not changing the A matrix. This seems to cause a significant increase in iterations and failures, which is why it is disabled for now. --- opm/autodiff/NewtonIterationBlackoilInterleaved.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index dd37214fa..6466a9dcc 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -390,11 +390,11 @@ namespace Opm OPM_THROW(std::logic_error, "formInterleavedSystem() requires 3 phases."); } +#if 1 // 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; @@ -465,6 +465,11 @@ namespace Opm // Create output as product of L with equations. A = L * total_residual.derivative()[0]; b = L * total_residual.value().matrix(); +#else + ADB total_residual = vertcatCollapseJacs(eqs_in); + A = total_residual.derivative()[0]; + b = total_residual.value().matrix(); +#endif } } // anonymous namespace From 8cbce1bfdf6b7254ecd9b0e29c85968254bd06d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 15 Jun 2015 16:12:23 +0200 Subject: [PATCH 04/13] Working interleaved solver implemented. --- .../NewtonIterationBlackoilInterleaved.cpp | 88 +++++++++++++++++-- .../NewtonIterationBlackoilInterleaved.hpp | 4 +- 2 files changed, 84 insertions(+), 8 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 6466a9dcc..443df27c0 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -70,6 +70,17 @@ #include #endif +// namespace Dune +// { +// typedef Dune::FieldMatrix MatrixBlockType; +// typedef Dune::BCRSMatrix Mat; + +// template <> +// struct MatrixDimension : public MatrixDimension +// { +// }; +// } + namespace Opm { @@ -77,8 +88,8 @@ namespace Opm typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; - typedef Dune::FieldMatrix MatrixBlockType; - typedef Dune::BCRSMatrix Mat; + // typedef Dune::FieldMatrix MatrixBlockType; + // typedef Dune::BCRSMatrix Mat; namespace { @@ -171,11 +182,67 @@ namespace Opm eqs[phase] = eqs[phase] * matbalscale[phase]; } - // Form interleaved system. + // // Write matrices to disk for analysis. + // const std::string eqnames[] = { "wat", "oil", "gas" }; + // const std::string varnames[] = { "p", "sw", "x" }; + // for (int p1 = 0; p1 < np; ++p1) { + // for (int p2 = 0; p2 < np; ++p2) { + // DuneMatrix m(eqs[p1].derivative()[p2]); + // std::string filename = "mat-"; + // filename += eqnames[p1]; + // filename += "-"; + // filename += varnames[p2]; + // writeMatrixToMatlab(m, filename.c_str()); + // } + // } + + // Find sparsity structure as union of basic block sparsity structures. + // Use addition to get to the union structure. + Eigen::SparseMatrix structure = eqs[0].derivative()[0]; + for (int p1 = 0; p1 < np; ++p1) { + for (int p2 = 0; p2 < np; ++p2) { + structure += eqs[p1].derivative()[p2]; + } + } + // DuneMatrix ms(structure); + // writeMatrixToMatlab(ms, "structurematrix"); + // Get row major form. + Eigen::SparseMatrix s = structure; + + // Form modified system. Eigen::SparseMatrix A; V b; formInterleavedSystem(np, eqs, A, b); + // Create ISTL matrix. + 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; + } + } + + // // Scale pressure equation. // const double pscale = 200*unit::barsa; // const int nc = residual.material_balance_eq[0].size(); @@ -186,14 +253,19 @@ namespace Opm SolutionVector dx(SolutionVector::Zero(b.size())); // Create ISTL matrix. - DuneMatrix istlA( A ); + // 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()); + 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; @@ -214,7 +286,11 @@ namespace Opm } // Copy solver output to dx. - std::copy(x.begin(), x.end(), dx.data()); + 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 ) { diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index 0bab86c9a..1a8503bd9 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -61,8 +61,8 @@ namespace Opm /// 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::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector Vector; From d86de7bb79685193f2a2899e260a9a5dac2d8bbd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 10:08:19 +0200 Subject: [PATCH 05/13] Only use pressure jacobian to form sparsity pattern. Also clean up by eliminating commented-out debugging code. --- .../NewtonIterationBlackoilInterleaved.cpp | 26 +++---------------- 1 file changed, 4 insertions(+), 22 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 443df27c0..ff4366db7 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -182,31 +182,13 @@ namespace Opm eqs[phase] = eqs[phase] * matbalscale[phase]; } - // // Write matrices to disk for analysis. - // const std::string eqnames[] = { "wat", "oil", "gas" }; - // const std::string varnames[] = { "p", "sw", "x" }; - // for (int p1 = 0; p1 < np; ++p1) { - // for (int p2 = 0; p2 < np; ++p2) { - // DuneMatrix m(eqs[p1].derivative()[p2]); - // std::string filename = "mat-"; - // filename += eqnames[p1]; - // filename += "-"; - // filename += varnames[p2]; - // writeMatrixToMatlab(m, filename.c_str()); - // } - // } - - // Find sparsity structure as union of basic block sparsity structures. + // 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 p1 = 0; p1 < np; ++p1) { - for (int p2 = 0; p2 < np; ++p2) { - structure += eqs[p1].derivative()[p2]; - } + for (int phase = 0; phase < np; ++phase) { + structure += eqs[phase].derivative()[0]; } - // DuneMatrix ms(structure); - // writeMatrixToMatlab(ms, "structurematrix"); - // Get row major form. Eigen::SparseMatrix s = structure; // Form modified system. From c8cae85ea21c0e71c6eb6d77556c51c46c37700a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 10:10:05 +0200 Subject: [PATCH 06/13] Move functions needed by several NewtonIteration-classes to separate file. --- CMakeLists_files.cmake | 4 +- opm/autodiff/NewtonIterationUtilities.cpp | 275 ++++++++++++++++++++++ opm/autodiff/NewtonIterationUtilities.hpp | 64 +++++ 3 files changed, 342 insertions(+), 1 deletion(-) create mode 100644 opm/autodiff/NewtonIterationUtilities.cpp create mode 100644 opm/autodiff/NewtonIterationUtilities.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f3f6b0259..eb38e90cd 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -29,8 +29,9 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/BlackoilPropsAdInterface.cpp opm/autodiff/ExtractParallelGridInformationToISTL.cpp opm/autodiff/NewtonIterationBlackoilCPR.cpp - opm/autodiff/NewtonIterationBlackoilSimple.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 @@ -114,6 +115,7 @@ list (APPEND PUBLIC_HEADER_FILES 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/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 From 9e288579330f56a601d1ee9d6b7ce79ce402de86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 10:18:11 +0200 Subject: [PATCH 07/13] Remove functions that were moved to NewtonIterationUtilities. --- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 278 +--------------- .../NewtonIterationBlackoilInterleaved.cpp | 312 +----------------- 2 files changed, 19 insertions(+), 571 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index ba920acac..87294fc9a 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -25,12 +25,15 @@ #include #include +#include #include #include #include #include #include #include + + #include // #include @@ -62,39 +65,6 @@ namespace Opm 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 +84,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. @@ -229,6 +201,10 @@ namespace Opm return dx; } + + + + const boost::any& NewtonIterationBlackoilCPR::parallelInformation() const { return parallelInformation_; @@ -236,242 +212,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 index ff4366db7..834681277 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -44,6 +44,7 @@ #include #include +#include #include #include #include @@ -70,16 +71,6 @@ #include #endif -// namespace Dune -// { -// typedef Dune::FieldMatrix MatrixBlockType; -// typedef Dune::BCRSMatrix Mat; - -// template <> -// struct MatrixDimension : public MatrixDimension -// { -// }; -// } namespace Opm { @@ -88,41 +79,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 interleaved 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 formInterleavedSystem(const int num_phases, - const std::vector& eqs, - Eigen::SparseMatrix& A, - V& b); - - } // anonymous namespace @@ -143,6 +99,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. @@ -194,9 +152,9 @@ namespace Opm // Form modified system. Eigen::SparseMatrix A; V b; - formInterleavedSystem(np, eqs, A, b); + formEllipticSystem(np, eqs, A, b); - // Create ISTL matrix. + // 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(); @@ -224,22 +182,9 @@ namespace Opm } } - - // // 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()); for (int i = 0; i < size; ++i) { @@ -274,8 +219,7 @@ namespace Opm dx(2*size + i) = x[i][2]; } - if( hasWells ) - { + if ( hasWells ) { // Compute full solution using the eliminated equations. // Recovery in inverse order of elimination. dx = recoverVariable(elim_eqs[1], dx, np); @@ -284,6 +228,10 @@ namespace Opm return dx; } + + + + const boost::any& NewtonIterationBlackoilInterleaved::parallelInformation() const { return parallelInformation_; @@ -291,246 +239,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 interleaved 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 formInterleavedSystem(const int num_phases, - const std::vector& eqs_in, - Eigen::SparseMatrix& A, - V& b) - { - if (num_phases != 3) { - OPM_THROW(std::logic_error, "formInterleavedSystem() requires 3 phases."); - } - -#if 1 - // 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(); -#else - ADB total_residual = vertcatCollapseJacs(eqs_in); - A = total_residual.derivative()[0]; - b = total_residual.value().matrix(); -#endif - } - - } // anonymous namespace } // namespace Opm From 4eb77bebb49465a3b371033a8c328d7f04538d49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 10:42:50 +0200 Subject: [PATCH 08/13] Further cleanup: includes, copyright, whitespace. --- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 24 +----------- .../NewtonIterationBlackoilInterleaved.cpp | 39 ------------------- 2 files changed, 2 insertions(+), 61 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index 87294fc9a..721b260c1 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -27,33 +27,16 @@ #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 { @@ -61,8 +44,6 @@ namespace Opm typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; - typedef Dune::FieldMatrix MatrixBlockType; - typedef Dune::BCRSMatrix Mat; @@ -191,8 +172,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); diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 834681277..0758c3cf1 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -1,26 +1,5 @@ /* 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). @@ -46,24 +25,6 @@ #include #include #include -#include -#include -#include -#include -#include -#include - -// #include -#include -#include -#include -#include -#include -#include -#include -#include - -#include #if HAVE_UMFPACK #include From ca74b18784af3e9e247680e44d6cbce57f667f75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 11:34:31 +0200 Subject: [PATCH 09/13] Add missing include. --- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index 721b260c1..10664e107 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #if HAVE_UMFPACK From 5476b7a6ac15ac4770e21c6844f45efd5e90fadc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 11:35:01 +0200 Subject: [PATCH 10/13] Clean up headers, copyright. --- .../NewtonIterationBlackoilInterleaved.hpp | 39 +++++++------------ 1 file changed, 13 insertions(+), 26 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index 1a8503bd9..eb46ca461 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -1,27 +1,5 @@ /* 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 @@ -43,14 +21,23 @@ along with OPM. If not, see . */ -#include +#ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED +#define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED + + #include -#include #include -#include + +#include + #include #include -#include +#include +#include +#include + +#include + #include namespace Opm From a3d115ff22398be3cada56b3a3ebf85f33dfbcb0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 11:35:23 +0200 Subject: [PATCH 11/13] Add (disabled) parallel version. Disabled because the constructPreconditionerAndSolve() method does not have a way currently to construct a parallel preconditioner. --- .../NewtonIterationBlackoilInterleaved.cpp | 33 ++++++++++++++++--- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 0758c3cf1..6e3cc752f 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -1,5 +1,7 @@ /* 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). @@ -25,6 +27,8 @@ #include #include #include +#include +#include #if HAVE_UMFPACK #include @@ -159,11 +163,30 @@ namespace Opm x = 0.0; Dune::InverseOperatorResult result; - // Construct operator, scalar product and vectors needed. - typedef Dune::MatrixAdapter Operator; - Operator opA(istlA); - Dune::Amg::SequentialInformation info; - constructPreconditionerAndSolve(opA, x, istlb, info, result); +// Parallel version is deactivated until we figure out how to do it properly. +#if 0 // 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; From 5e513642d72580dd0019c7a8d80ab595f5694159 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 13:28:11 +0200 Subject: [PATCH 12/13] Add paralell preconditioner, enable parallel case. --- .../NewtonIterationBlackoilInterleaved.cpp | 2 +- .../NewtonIterationBlackoilInterleaved.hpp | 61 +++++++++++++++---- 2 files changed, 49 insertions(+), 14 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 6e3cc752f..d8851bc26 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -164,7 +164,7 @@ namespace Opm Dune::InverseOperatorResult result; // Parallel version is deactivated until we figure out how to do it properly. -#if 0 // HAVE_MPI +#if HAVE_MPI if (parallelInformation_.type() == typeid(ParallelISTLInformation)) { typedef Dune::OwnerOverlapCopyCommunication Comm; diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index eb46ca461..f6aacfb69 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -27,6 +27,7 @@ #include #include +#include #include @@ -34,6 +35,7 @@ #include #include #include +#include #include #include @@ -87,33 +89,66 @@ namespace Opm /// \brief construct the CPR preconditioner and the solver. /// \tparam P The type of the parallel information. /// \param parallelInformation the information about the parallelization. - template + template void constructPreconditionerAndSolve(O& opA, Vector& x, Vector& istlb, - const P& parallelInformation, + const POrComm& parallelInformation, 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); - const double relax = 1.0; - Preconditioner precond(opA.getmat(), relax); + // 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, + 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, + Dune::BiCGSTABSolver linsolve(opA, sp, precond, linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_); // Solve system. linsolve.apply(x, istlb, result); From 5002fe1e37c3b263ea0ee091de371bc0dd66bef9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 16 Jun 2015 13:33:18 +0200 Subject: [PATCH 13/13] Add use_interleaved option to flow_cp. --- examples/flow_cp.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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));