From e22ade91642c25fcdd4567f7d1ddd9a29dfef747 Mon Sep 17 00:00:00 2001 From: weiwongg Date: Wed, 18 Sep 2024 11:00:39 +0200 Subject: [PATCH] add multigrid module --- Bembel/MultiGrid | 23 +++ Bembel/src/MultiGrid/MultiGridSolver.hpp | 73 +++++++++ Bembel/src/MultiGrid/Prolongation.hpp | 60 ++++++++ Bembel/src/MultiGrid/ProlongationBsplines.hpp | 140 +++++++++++++++++ examples/LaplaceBeltramiMultiGrid.cpp | 144 ++++++++++++++++++ 5 files changed, 440 insertions(+) create mode 100644 Bembel/MultiGrid create mode 100644 Bembel/src/MultiGrid/MultiGridSolver.hpp create mode 100644 Bembel/src/MultiGrid/Prolongation.hpp create mode 100644 Bembel/src/MultiGrid/ProlongationBsplines.hpp create mode 100644 examples/LaplaceBeltramiMultiGrid.cpp diff --git a/Bembel/MultiGrid b/Bembel/MultiGrid new file mode 100644 index 000000000..c6dd02465 --- /dev/null +++ b/Bembel/MultiGrid @@ -0,0 +1,23 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#ifndef BEMBEL_MULTIGRID_MODULE_ +#define BEMBEL_MULTIGRID_MODULE_ + +#include +#include +#include + +#include "AnsatzSpace" +#include "Spline" + +#include "src/MultiGrid/ProlongationBsplines.hpp" +#include "src/MultiGrid/Prolongation.hpp" +#include "src/MultiGrid/MultiGridSolver.hpp" + +#endif diff --git a/Bembel/src/MultiGrid/MultiGridSolver.hpp b/Bembel/src/MultiGrid/MultiGridSolver.hpp new file mode 100644 index 000000000..e0f6481dc --- /dev/null +++ b/Bembel/src/MultiGrid/MultiGridSolver.hpp @@ -0,0 +1,73 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#ifndef BEMBEL_SRC_MULTIGRID_MULTIGRIDSOLVER_HPP_ +#define BEMBEL_SRC_MULTIGRID_MULTIGRIDSOLVER_HPP_ + +namespace Bembel { +namespace MG { + +#define MULTIGRID_minLvl 1 +#define MULTIGRID_cycleType 2 +#define MULTIGRID_preSmooth 3 +#define MULTIGRID_postSmooth 3 + +/** + * \ingroup MultiGrid + * \brief Add Gauss Seidel smoother for use with multigrid. + **/ +template +void GaussSeidelSmoother(const Derived &S, const Eigen::MatrixBase &b, + Eigen::MatrixBase *x, + const unsigned int steps = 3) { + for (auto i = 0; i < steps; ++i) + x->derived() = S.template triangularView().solve( + b - S.template triangularView() * x->derived()); +} + +/** + * \ingroup MultiGrid + * \brief Implements a single step of multiplicative multigrid. + **/ +template +void multiplicativeMultiGrid(const Derived &Ss, + const Eigen::MatrixBase &f, + Eigen::MatrixBase *x, const Derived4 &Ps, + const unsigned int lvl = MULTIGRID_minLvl) { + typedef typename Derived3::Scalar Scalar; + typedef Eigen::Matrix Matrix; + if (lvl <= MULTIGRID_minLvl) { + Matrix A = Matrix(Ss[lvl]); + x->derived() = A.fullPivHouseholderQr().solve(f); + } else { + // Pre-smooth + GaussSeidelSmoother(Ss[lvl], f, x, MULTIGRID_preSmooth); + Matrix r = f - Ss[lvl] * x->derived(); + // Restrict the residual + r = Ps[lvl - 1].transpose() * r; + Matrix y = r; + y.setZero(); + // Recursively calling multiplicativeMultiGrid for computing the error + for (auto i = 0; i < MULTIGRID_cycleType; ++i) + multiplicativeMultiGrid(Ss, r, &y, Ps, lvl - 1); + // Prolongate the error + y = Ps[lvl - 1] * y; + // Correction + x->derived() += y; + // Post-smooth + GaussSeidelSmoother(Ss[lvl], f, x, MULTIGRID_postSmooth); + } +} + +} // namespace MG +} // namespace Bembel +#endif // BEMBEL_SRC_MULTIGRID_MULTIGRIDSOLVER_HPP_ diff --git a/Bembel/src/MultiGrid/Prolongation.hpp b/Bembel/src/MultiGrid/Prolongation.hpp new file mode 100644 index 000000000..17b265f11 --- /dev/null +++ b/Bembel/src/MultiGrid/Prolongation.hpp @@ -0,0 +1,60 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#ifndef BEMBEL_SRC_MULTIGRID_PROLONGATION_HPP_ +#define BEMBEL_SRC_MULTIGRID_PROLONGATION_HPP_ + +namespace Bembel { +namespace MG { + +/** + * \ingroup MultiGrid + * \brief Return the prolongation matrix from level l-1 to level l. + **/ +template +Eigen::SparseMatrix prolongationMatrix( + const AnsatzSpace &ansatz_space, unsigned int level) { + assert(level > 0 && + "Level must be 1 or larger for build prolongation matrix from level-1 " + "to level"); + unsigned int knot_repetition = ansatz_space.get_knot_repetition(); + unsigned int polynomial_degree = ansatz_space.get_polynomial_degree(); + unsigned int num_patches = ansatz_space.get_number_of_patches(); + + SuperSpace coarse_super_space; + coarse_super_space.init_SuperSpace( + Geometry(ansatz_space.get_superspace().get_geometry()), level - 1, + polynomial_degree); + Projector coarse_proj(coarse_super_space, knot_repetition); + Glue coarse_glue(coarse_super_space, coarse_proj); + + SuperSpace fine_super_space; + fine_super_space.init_SuperSpace( + Geometry(ansatz_space.get_superspace().get_geometry()), level, + polynomial_degree); + Projector fine_proj(fine_super_space, knot_repetition); + Glue fine_glue(fine_super_space, fine_proj); + + Eigen::VectorXd degree_vector = + fine_glue.get_glue_matrix().transpose() * + Eigen::VectorXd::Ones(fine_glue.get_glue_matrix().rows()); + // the reason we multiply 0.5 in the begining is because we use the scaled + // B-spline basis function in the Bembel + return 0.5 * degree_vector.cwiseInverse().asDiagonal() * + fine_glue.get_glue_matrix().transpose() * + MG::prolongateBsplinesMultiPatches2D(polynomial_degree, level, + num_patches, knot_repetition) * + coarse_glue.get_glue_matrix(); +} + +} // namespace MG +} // namespace Bembel +#endif // BEMBEL_SRC_MULTIGRID_PROLONGATION_HPP_ diff --git a/Bembel/src/MultiGrid/ProlongationBsplines.hpp b/Bembel/src/MultiGrid/ProlongationBsplines.hpp new file mode 100644 index 000000000..4da80af10 --- /dev/null +++ b/Bembel/src/MultiGrid/ProlongationBsplines.hpp @@ -0,0 +1,140 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#ifndef BEMBEL_SRC_MULTIGRID_PROLONGATIONBSPLINES_HPP_ +#define BEMBEL_SRC_MULTIGRID_PROLONGATIONBSPLINES_HPP_ + +namespace Bembel { +namespace MG { + +template +std::vector> toEigenTriplets( + Eigen::SparseMatrix &M) { + std::vector> v; + for (int i = 0; i < M.outerSize(); i++) + for (typename Eigen::SparseMatrix::InnerIterator it(M, i); it; ++it) + v.emplace_back(it.row(), it.col(), it.value()); + return v; +} + +/** + * \ingroup MultiGrid + * \brief Return the nonzero entries info of the row id of the prolongation + * matrix using Oslo-Algorithm + * https://www.uio.no/studier/emner/matnat/ifi/nedlagte-emner/INF-MAT5340/v07/undervisningsmateriale/kap4.pdf + */ +std::tuple nonZerosInfoRow( + const std::vector &knots, + const std::vector &uniform_refined_knots, + unsigned int polynomial_degree, unsigned int row_id, + unsigned int guess = 0) { + unsigned int d = guess; + while (uniform_refined_knots[row_id] >= knots[d]) ++d; + --d; + Eigen::MatrixXd non_zero_alpha = + Eigen::MatrixXd::Ones(1, polynomial_degree + 1); + for (auto i = 0; i < polynomial_degree; ++i) { + Eigen::MatrixXd R = Eigen::MatrixXd::Zero(i + 1, i + 2); + for (auto j = 0; j < i + 1; ++j) { + R(j, j) = (knots[d + j + 1] - uniform_refined_knots[row_id + i + 1]) / + (knots[d + j + 1] - knots[d + j - i]); + R(j, j + 1) = 1.0 - R(j, j); + } + non_zero_alpha.block(0, 0, 1, i + 2) = + non_zero_alpha.block(0, 0, 1, i + 1) * R; + } + + return std::make_tuple(non_zero_alpha, d); +} + +/** + * \ingroup MultiGrid + * \brief Return prolongation matrix from level-1 to level in 1D. + */ +Eigen::SparseMatrix prolongateBsplines1D( + unsigned int polynomial_degree, unsigned int level, + unsigned int knotrepetition = 1) { + assert(level > 0 && + "Level must be 1 or larger for build prolongation matrix from level-1 " + "to level"); + std::vector knots = Bembel::Spl::MakeUniformKnotVector( + polynomial_degree + 1, ((1 << (level - 1)) - 1), knotrepetition); + std::vector uniform_refined_knots = Spl::MakeUniformKnotVector( + polynomial_degree + 1, ((1 << level) - 1), knotrepetition); + unsigned int num_basis = knots.size() - polynomial_degree - 1; + unsigned int num_fine_basis = + uniform_refined_knots.size() - polynomial_degree - 1; + typedef typename Eigen::Triplet Triplet; + std::vector tripletList; + tripletList.reserve(num_fine_basis * (polynomial_degree + 1)); + int d = 0; + for (int i = 0; i < num_fine_basis; ++i) { + for (int j = 0; j < polynomial_degree + 1; ++j) { + Eigen::MatrixXd non_zero_alpha; + std::tie(non_zero_alpha, d) = nonZerosInfoRow( + knots, uniform_refined_knots, polynomial_degree, i, d); + tripletList.push_back( + Triplet(i, d - polynomial_degree + j, non_zero_alpha(j))); + } + } + Eigen::SparseMatrix P; + P.resize(num_fine_basis, num_basis); + P.setFromTriplets(tripletList.begin(), tripletList.end()); + return P; +} + +/** + * \ingroup MultiGrid + * \brief Return prolongation matrix from level-1 to level in 2D. + */ +Eigen::SparseMatrix prolongateBsplines2D( + unsigned int polynomial_degree, unsigned int level, + unsigned int knotrepetition = 1) { + assert(level > 0 && + "Level must be 1 or larger for build prolongation matrix from level-1 " + "to level"); + Eigen::SparseMatrix P_ = + prolongateBsplines1D(polynomial_degree, level, knotrepetition); + Eigen::SparseMatrix P = Eigen::kroneckerProduct(P_, P_); + return P; +} + +/** + * \ingroup MultiGrid + * \brief Return prolongation matrix from level-1 to level in 2D in case of + * multiple patches. + */ +Eigen::SparseMatrix prolongateBsplinesMultiPatches2D( + unsigned int polynomial_degree, unsigned int level, + unsigned int num_patches, unsigned int knotrepetition = 1) { + assert(level > 0 && + "Level must be 1 or larger for build prolongation matrix from level-1 " + "to level"); + Eigen::SparseMatrix P_ = + prolongateBsplines2D(polynomial_degree, level, knotrepetition); + std::vector> triplets_ = toEigenTriplets(P_); + std::vector> triplets; + Eigen::SparseMatrix P(num_patches * P_.rows(), + num_patches * P_.cols()); + for (auto i = 0; i < num_patches; ++i) { + for (auto entry : triplets_) { + triplets.push_back(Eigen::Triplet(entry.row() + i * P_.rows(), + entry.col() + i * P_.cols(), + entry.value())); + } + } + P.setFromTriplets(triplets.begin(), triplets.end()); + return P; +} + +} // namespace MG +} // namespace Bembel +#endif // BEMBEL_SRC_MULTIGRID_PROLONGATIONBSPLINES_HPP_ diff --git a/examples/LaplaceBeltramiMultiGrid.cpp b/examples/LaplaceBeltramiMultiGrid.cpp new file mode 100644 index 000000000..89249a8b5 --- /dev/null +++ b/examples/LaplaceBeltramiMultiGrid.cpp @@ -0,0 +1,144 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +// MG solver for the Laplace Beltrami with Neumann condition. +int mmg(const Eigen::VectorXd &constant_one, + const Eigen::SparseMatrix &S, + const Eigen::Matrix &L, + const Eigen::Matrix &f, + Eigen::Matrix &x, + const std::vector> &Ps, int lvl, + double tol = 1e-6) { + // stiffness matrices + std::vector> Ss(lvl + 1); + // the stiffness matrices on each level + Ss[lvl] = S; + for (int i = lvl - 1; i >= 0; --i) + Ss[i] = Ps[i].transpose() * Ss[i + 1] * Ps[i]; + + // Multi Grid Part + int iter = 0; + // initialize the x with zeros + x.resize(f.size()); + x.setZero(); + Eigen::VectorXd res; + do { + // repetitively call multiplicativeMultiGrid() + Bembel::MG::multiplicativeMultiGrid(Ss, f, &x, Ps, lvl); + // minus regularization terms such that the integral of function represented + // by x is 0. L.dot(x) is integral of x and L.dot(constant_one) is integral + // of constant one function. + x = x - L.dot(x) / L.dot(constant_one) * constant_one; + res = f - S * x; + ++iter; + } while (res.norm() > tol); + return iter; +} + +int main() { + using namespace Bembel; + using namespace Eigen; + IO::Stopwatch sw; + std::function fun = [](const Vector3d &in) { + // Spherical harmonics + return sqrt(3 / BEMBEL_PI) * in(2); + }; + std::function refsol = [](const Vector3d &in) { + return 0.5 * sqrt(3 / BEMBEL_PI) * in(2); + }; + + std::function const_fun = [](const Vector3d &in) { + return 1.0; + }; + Geometry geometry("sphere.dat"); + std::cout << "\n=============================================================" + "==========\n"; + // Iterate over polynomial degree. + for (auto polynomial_degree : {1, 2, 3, 4}) { + std::cout << "Degree " << polynomial_degree << std::endl; + int MAX_LVL = 7; + + std::vector> Ps; + + // Iterate over refinement levels + std::cout << std::left << std::setw(8) << "Level" << std::left + << std::setw(8) << "Iters" << std::left << std::setw(15) + << "L2 norm" << std::left << std::setw(15) << "Time/sec" + << std::endl; + + for (int refinement_level = 0; refinement_level <= MAX_LVL; + ++refinement_level) { + // Build ansatz space + AnsatzSpace ansatz_space( + geometry, refinement_level, polynomial_degree); + if (refinement_level != 0) + Ps.push_back(MG::prolongationMatrix(ansatz_space, refinement_level)); + // Set up and compute discrete operator + DiscreteLocalOperator disc_op(ansatz_space); + disc_op.compute(); + + DiscreteLinearForm, LaplaceBeltramiOperator> + disc_lf_const(ansatz_space); + disc_lf_const.get_linear_form().set_function(const_fun); + disc_lf_const.compute(); + + DiscreteLinearForm, LaplaceBeltramiOperator> + disc_lf(ansatz_space); + disc_lf.get_linear_form().set_function(fun); + disc_lf.compute(); + + AnsatzSpace temp_ansatz_space( + geometry, refinement_level, polynomial_degree); + DiscreteLocalOperator disc_identity_op( + temp_ansatz_space); + disc_identity_op.compute(); + disc_identity_op.get_discrete_operator(); + + // prepare solver for the mass matrix + Eigen::SimplicialLLT> solver; + // Compute the numerical factorization + solver.compute(disc_identity_op.get_discrete_operator()); + Eigen::VectorXd constant_one = + solver.solve(disc_lf_const.get_discrete_linear_form()); + + // solve linear system of equation with MG + Eigen::Matrix x; + sw.tic(); + int iters = mmg(constant_one, disc_op.get_discrete_operator(), + disc_lf_const.get_discrete_linear_form(), + disc_lf.get_discrete_linear_form(), x, Ps, + refinement_level, 1e-10); + auto difft = sw.toc(); + // print INFO + auto err = surfaceL2error(ansatz_space, x, refsol); + std::cout << std::left << std::setw(8) << refinement_level << std::left + << std::setw(8) << iters << std::left << std::setw(15) << err + << std::left << std::setw(15) << difft << std::endl; + } + } + std::cout << "=============================================================" + "==========" + << std::endl; + + return 0; +}