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] 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