Home     Back     Next     1    2    3    4    5 Top of current page     Browse Files

(Simple) Multigrid solver

Content

This example a multigrid solver that uses the grids \(\Omega_{h}\), \(\Omega_{2h}\), \(\Omega_{4h}\), \(\Omega_{8h}\)

One iteration a multigrid V-cycle

Code

session31/ex01/smooth.cpp 
#include <cmath>

#include <hpc/matvec/axpy.hpp>
#include <hpc/matvec/copy.hpp>
#include <hpc/matvec/densevector.hpp>
#include <hpc/matvec/iterators.hpp>
#include <hpc/matvec/traits.hpp>
#include <hpc/matvec/print.hpp>

//------------------------------------------------------------------------------

using namespace hpc;
using namespace hpc::matvec;


template <typename T>
const T &
assertEqual(const T &a, const T &b)
{
    assert(a==b);
    return a;
}

//------------------------------------------------------------------------------

template <
    typename F, typename G,
    typename T, template<typename> class VectorF,
                template<typename> class VectorU,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>> > = true
>
void
discretizeProblem(const F &f, const G &g, VectorF<T> &f_h, VectorU<T> &u_h)
{
    auto N  = assertEqual(f_h.length(), u_h.length()) - 2;
    auto h  = T(1)/(N+1);

    u_h(0)   = g(0);
    u_h(N+1) = g(1);

    for (std::size_t i=1; i<N; ++i) {
        u_h(i) = 0;
        f_h(i) = f(i*h);
    }
    // actually we don't have to initialize f_h(0), f_h(N+1) so we store some
    // good and evil numbers there:
    f_h(0)   = 42;
    f_h(N+1) = 666;
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class Vector,
                Require< Dense<Vector<T>> > = true
>
T
lInfNorm(const Vector<T> &x_h)
{
    T res = 0;

    for (auto [i, xi] : x_h) {
        i=i;
        if (std::abs(xi) > res) {
            res = std::abs(xi);
        }
    }
    return res;
}

//------------------------------------------------------------------------------

template <
    typename U,
    typename T, template<typename> class VectorU,
                template<typename> class VectorE,
                Require< Dense<VectorU<T>>,
                         Dense<VectorE<T>> > = true
>
void
evalError(const U &u, const VectorU<T> &u_h, VectorE<T> &e_h)
{
    auto N  = assertEqual(u_h.length(), e_h.length()) - 2;
    auto h  = T(1)/(N+1);

    for (auto [i, ei] : e_h) {
        ei = u(i*h) - u_h(i);
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorU,
                template<typename> class VectorF,
                template<typename> class VectorR,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>>,
                         Dense<VectorR<T>> > = true
>
void
evalResidual(const VectorU<T> &u_h, const VectorF<T> &f_h, VectorR<T> &r_h)
{
    auto N  = assertEqual(u_h.length(), r_h.length()) - 2;
    auto h  = T(1)/(N+1);
    auto h2 = h*h;

    r_h(0) = r_h(N+1) = 0;
    for (std::size_t i=1; i<=N; ++i) {
        r_h(i) = f_h(i) - 1./h2 * (2*u_h(i) - u_h(i-1) - u_h(i+1));
    }
}


//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorF,
                template<typename> class VectorU,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>> > = true
>
void
direct_solver(const VectorF<T> &f_h, VectorU<T> &u_h)
{
    auto N  = assertEqual(f_h.length(), u_h.length()) - 2;
    auto h  = T(1)/(N+1);
    auto h2 = h*h;

    // forward substitution
    u_h(1) = h2*f_h(1) + u_h(0);
    for (std::size_t k=2; k<N; ++k) {
        u_h(k) = h2*f_h(k) + u_h(k-1)*(k-1)/k;
    }
    u_h(N) = h2*f_h(N) + u_h(N+1) + u_h(N-1)*(N-1)/N;

    // backward substitution
    u_h(N) = u_h(N)*N/(N+1);
    for (std::size_t k=N-1; k>=1; --k) {
        u_h(k) = (u_h(k)+u_h(k+1))*k/(k+1);
    }
}

//==============================================================================

template <
    typename T, template<typename> class Vector,
                Require< Dense<Vector<T>> > = true
>
void
addNoise(std::size_t k, Vector<T> &x_h)
{
    auto N  = x_h.length()-2;
    auto h  = T(1)/(N+1);

    for (std::size_t i=1; i<=N; ++i) {
        x_h(i) += sin(M_PI*k*i*h);
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorF,
                template<typename> class VectorU,
                template<typename> class VectorS,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>>,
                         Dense<VectorS<T>> > = true
>
void
jacobi_step(double omega, const VectorF<T> &f_h, const VectorU<T> &u_h,
            VectorS<T> &s_h)
{
    auto N  = assertEqual(f_h.length(), u_h.length()) - 2;
    auto h  = T(1)/(N+1);
    auto h2 = h*h;

    // TODO: Your code

    s_h(0)   = u_h(0);
    s_h(N+1) = u_h(N+1);
}

//------------------------------------------------------------------------------

int
main()
{
    /*
    auto f = [](double x) { return -2; };
    auto g = [](double x) { return 1 + 3*x; };
    auto u = [](double x) { return (x+1)*(x+1); };
    */
    auto f = [](double x) { return 0; };
    auto g = [](double x) { return 0; };
    auto u = [](double x) { return 0; };

    std::size_t N = 100;

    DenseVector<double>  f_h(N+2), u_h(N+2), tmp_h(N+2), r_h(N+2), e_h(N+2);

    discretizeProblem(f, g, f_h, u_h);

    // First initialize u_h with exact solution
    direct_solver(f_h, u_h);

    // then add some noise 
    addNoise(1, u_h);
    //addNoise(N/2, u_h);
    //addNoise(N-3, u_h);
    addNoise(N-1, u_h);

    double omega = 0.5;     // smoothing parameter

    for (std::size_t it=0; it<=20; ++it) {
        jacobi_step(omega, f_h, u_h, tmp_h);
        jacobi_step(omega, f_h, tmp_h, u_h);

        evalResidual(u_h, f_h, r_h);
        evalError(u, u_h, e_h);

        fmt::printf("#lInfNorm(r_h) = %.5e\n", lInfNorm(r_h));
        fmt::printf("#lInfNorm(e_h) = %.5e\n", lInfNorm(e_h));
        fmt::printf("it %3d\n", it);
        for (auto [i,ri] : r_h) {
            fmt::printf("%d %f\n", i, ri);
        }
        fmt::printf("\n\n");
    }
}

Code for gnuplot

session32/ex01/restrict_post_smooth.gnuplot 
set terminal svg size 900, 500
set output "restrict_post_smooth.svg"
set xlabel "Grid points"
set ylabel "Residual value"
set title "Smoothing"
set key outside
set key autotitle columnhead
set pointsize 0.5

plot "restrict_post_smooth.data" index  0 using 1:2 with linespoints lt 1 lw 1, \
     "restrict_post_smooth.data" index  1 using 1:2 with linespoints lt 2 lw 3

How to compile and plot results

theon$ g++ -Wall -I /home/numerik/pub/hpc/ws18/session25/ -std=c++17 -o restrict_post_smooth restrict_post_smooth.cpp
theon$ ./restrict_post_smooth > restrict_post_smooth.data
theon$ 

Plotting results

theon$ gnuplot restrict_post_smooth.gnuplot
theon$ gnuplot restrict_post_smooth2.gnuplot
theon$ gnuplot restrict_post_smooth3.gnuplot
theon$ gnuplot restrict_post_smooth4.gnuplot
theon$ gnuplot restrict_post_smooth5.gnuplot
theon$ gnuplot restrict_post_smooth6.gnuplot
theon$ gnuplot restrict_post_smooth7.gnuplot
theon$ gnuplot restrict_post_smooth8.gnuplot
theon$ 

Plot of smoothing properties

Restriction on \(\Omega_{2h}\) and post smoothing

Restriction on \(\Omega_{4h}\) and post smoothing

Restriction on \(\Omega_{8h}\) and post smoothing

Multigrid with 4 levels

Code

session32/ex01/mg_with_4_levels.cpp 
#include <cmath>
#include <sstream>
#include <memory>
#include <vector>

#include <hpc/matvec/axpy.hpp>
#include <hpc/matvec/copy.hpp>
#include <hpc/matvec/densevector.hpp>
#include <hpc/matvec/iterators.hpp>
#include <hpc/matvec/traits.hpp>
#include <hpc/matvec/print.hpp>

//------------------------------------------------------------------------------

using namespace hpc;
using namespace hpc::matvec;


template <typename T>
const T &
assertEqual(const T &a, const T &b)
{
    assert(a==b);
    return a;
}

//------------------------------------------------------------------------------

template <
    typename F, typename G,
    typename T, template<typename> class VectorF,
                template<typename> class VectorU,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>> > = true
>
void
discretizeProblem(const F &f, const G &g, VectorF<T> &f_h, VectorU<T> &u_h)
{
    auto N  = assertEqual(f_h.length(), u_h.length()) - 2;
    auto h  = T(1)/(N+1);

    u_h(0)   = g(0);
    u_h(N+1) = g(1);

    for (auto [i, fi] : f_h) {
        fi = f(i*h);
    }
}

//------------------------------------------------------------------------------

template <
    typename U,
    typename T, template<typename> class VectorU,
                template<typename> class VectorE,
                Require< Dense<VectorU<T>>,
                         Dense<VectorE<T>> > = true
>
void
evalError(const U &u, const VectorU<T> &u_h, VectorE<T> &e_h)
{
    auto N  = assertEqual(u_h.length(), e_h.length()) - 2;
    auto h  = T(1)/(N+1);

    for (auto [i, ei] : e_h) {
        ei = u(i*h) - u_h(i);
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorU,
                template<typename> class VectorF,
                template<typename> class VectorR,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>>,
                         Dense<VectorR<T>> > = true
>
void
evalResidual(const VectorU<T> &u_h, const VectorF<T> &f_h, VectorR<T> &r_h)
{
    auto N  = assertEqual(u_h.length(), r_h.length()) - 2;
    auto h  = T(1)/(N+1);
    auto h2 = h*h;

    r_h(0) = r_h(N+1) = 0;
    for (std::size_t i=1; i<=N; ++i) {
        r_h(i) = f_h(i) - 1./h2 * (2*u_h(i) - u_h(i-1) - u_h(i+1));
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorR1,
                template<typename> class VectorR2,
                Require< Dense<VectorR1<T>>,
                         Dense<VectorR2<T>> > = true
>
void
restriction(const VectorR1<T> &r_h, VectorR2<T> &r_2h)
{
    auto N  = r_h.length() -2;
    auto n  = r_2h.length() -2;

    assert(2*(n+1) == N+1);

    r_2h(0) = r_2h(n+1) = 0;
    for (std::size_t i=1; i<=n; ++i) {
        r_2h(i) = 0.25*r_h(2*i-1) + 0.5*r_h(2*i) + 0.25*r_h(2*i+1);
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorU1,
                template<typename> class VectorU2,
                Require< Dense<VectorU1<T>>,
                         Dense<VectorU2<T>> > = true
>
void
updateProlongation(const VectorU1<T> &u_2h, VectorU2<T> &u_h)
{
    auto N  = u_h.length() -2;
    auto n  = u_2h.length() -2;

    assert(2*(n+1) == N+1);

    for (std::size_t i=1; i<=n; ++i) {
        u_h(2*i-1) += 0.5 * u_2h(i);
        u_h(2*i  ) += 1.0 * u_2h(i);
        u_h(2*i+1) += 0.5 * u_2h(i);
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class Vector,
                Require< Dense<Vector<T>> > = true
>
void
addNoise(std::size_t k, Vector<T> &x_h)
{
    auto N  = x_h.length()-2;
    auto h  = T(1)/(N+1);

    for (std::size_t i=1; i<=N; ++i) {
        x_h(i) += sin(M_PI*k*i*h);
    }
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class Vector,
                Require< Dense<Vector<T>> > = true
>
T
lInfNorm(const Vector<T> &x_h)
{
    T res = 0;

    for (auto [i, xi] : x_h) {
        i=i;
        if (std::abs(xi) > res) {
            res = std::abs(xi);
        }
    }
    return res;
}

template <
    typename T, template<typename> class Vector1,
                template<typename> class Vector2,
                Require< Dense<Vector1<T>>,
                         Dense<Vector2<T>> > = true
>
T
lInfNorm(const Vector1<T> &x1_h, const Vector2<T> &x2_h)
{
    auto N  = assertEqual(x1_h.length(), x2_h.length());

    T res = 0;

    for (std::size_t i=0; i<N; ++i) {
        if (std::abs(x1_h(i) - x2_h(i)) > res) {
            res = std::abs(x1_h(i) - x2_h(i));
        }
    }
    return res;
}


//------------------------------------------------------------------------------

template <
    typename T, template<typename> class Vector,
                Require< Dense<Vector<T>> > = true
>
T
l1Norm(const Vector<T> &x_h)
{
    auto N  = x_h.length()-2;
    auto h  = T(1)/(N+1);

    T res = 0;

    for (std::size_t i=1; i<=N; ++i) {
        i=i;
        res += 0.5*(std::abs(x_h(i-1)) + std::abs(x_h(i)));
    }
    return h*res;
}

//------------------------------------------------------------------------------


template <
    typename T, template<typename> class VectorF,
                template<typename> class VectorU,
                template<typename> class VectorS,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>>,
                         Dense<VectorS<T>> > = true
>
void
jacobi_step(double omega, const VectorF<T> &f_h, const VectorU<T> &u_h,
            VectorS<T> &s_h)
{
    auto N  = assertEqual(f_h.length(), u_h.length()) - 2;
    auto h  = T(1)/(N+1);
    auto h2 = h*h;

    for (std::size_t i=1; i<=N; ++i) {
        s_h(i) = (1-omega)*u_h(i) + omega*0.5*(u_h(i-1) +u_h(i+1)+h2*f_h(i));
    }
    s_h(0)   = u_h(0);
    s_h(N+1) = u_h(N+1);
}

//------------------------------------------------------------------------------

template <
    typename T, template<typename> class VectorF,
                template<typename> class VectorU,
                Require< Dense<VectorF<T>>,
                         Dense<VectorU<T>> > = true
>
void
direct_solver(const VectorF<T> &f_h, VectorU<T> &u_h)
{
    auto N  = assertEqual(f_h.length(), u_h.length()) - 2;
    auto h  = T(1)/(N+1);
    auto h2 = h*h;

    // forward substitution
    u_h(1) = h2*f_h(1) + u_h(0);
    for (std::size_t k=2; k<N; ++k) {
        u_h(k) = h2*f_h(k) + u_h(k-1)*(k-1)/k;
    }
    u_h(N) = h2*f_h(N) + u_h(N+1) + u_h(N-1)*(N-1)/N;

    // backward substitution
    u_h(N) = u_h(N)*N/(N+1);
    for (std::size_t k=N-1; k>=1; --k) {
        u_h(k) = (u_h(k)+u_h(k+1))*k/(k+1);
    }
}

//------------------------------------------------------------------------------

template <typename T, template<typename> class Vector>
void
init_zero(Vector<T> &u_h)
{
    for (auto [i, ui] : u_h) {
        i = i;
        ui = T(0);
    }
}

//------------------------------------------------------------------------------

template <typename T, template<typename> class Vector>
void
print_record(const Vector<T> &r_h, const char *descr)
{
    auto N  = r_h.length() -2;

    fmt::printf("\"%s\"\n", descr);
    for (auto [i, ri] : r_h) {
        fmt::printf("%f %f\n", double(i)/(N+1), ri);
    }
    fmt::printf("\n\n");
}

int
main()
{
    /*
    auto f = [](double x) { return 0; };
    auto g = [](double x) { return 0; };
    auto u = [](double x) { return 0; };
    */

    /*
    auto f = [](double x) { return -2; };
    auto g = [](double x) { return 1 + 3*x; };
    auto u = [](double x) { return (x+1)*(x+1); };
    */

    auto f = [](double x) {
        return -8*M_PI*cos(2*M_PI*exp(2*x))
              +16*M_PI*M_PI*exp(2*x)*exp(2*x)*sin(2*M_PI*exp(2*x));
    };

    auto g = [](double x) {
        return x*sin(2*M_PI*exp(2.));
    };

    auto u = [](double x) {
        return sin(2*M_PI*exp(2*x));
    };

    std::size_t N0 = 512;
    std::size_t N1 = 256;
    std::size_t N2 = 128;
    std::size_t N3 = 64;

    DenseVector<double> f_h(N0+1), u_h(N0+1), tmp_h(N0+1), r_h(N0+1), e_h(N0+1);
    DenseVector<double> f_2h(N1+1), u_2h(N1+1), tmp_2h(N1+1), r_2h(N1+1);
    DenseVector<double> f_4h(N2+1), u_4h(N2+1), tmp_4h(N2+1), r_4h(N2+1);
    DenseVector<double> f_8h(N3+1), u_8h(N3+1), tmp_8h(N3+1), r_8h(N3+1);

    discretizeProblem(f, g, f_h, u_h);

    // First initialize u_h with exact solution
    direct_solver(f_h, u_h);

    // then add some noise 
    addNoise((N0-2)/16, u_h);
    addNoise((N0-2)/4, u_h);
    addNoise((N0-2)-3, u_h);
    addNoise((N0-2)-1, u_h);

    double omega = 0.5;     // smoothing parameter

    evalResidual(u_h, f_h, r_h);
    print_record(r_h, "Initial residual");

    for (int it=0; it<5; ++it) {
        // On Omega_{h}
        // ------------
        jacobi_step(omega, f_h, u_h, tmp_h);
        jacobi_step(omega, f_h, tmp_h, u_h);
        evalResidual(u_h, f_h, r_h);

        // On Omega_{2h}
        // -------------
        restriction(r_h, f_2h);
        init_zero(u_2h);
        evalResidual(u_2h, f_2h, r_2h);

        jacobi_step(omega, f_2h, u_2h, tmp_2h);
        jacobi_step(omega, f_2h, tmp_2h, u_2h);

        evalResidual(u_2h, f_2h, r_2h);

        // On Omega_{4h}
        // -------------
        restriction(r_2h, f_4h);
        init_zero(u_4h);
        evalResidual(u_4h, f_4h, r_4h);

        jacobi_step(omega, f_4h, u_4h, tmp_4h);
        jacobi_step(omega, f_4h, tmp_4h, u_4h);

        evalResidual(u_4h, f_4h, r_4h);

        // On Omega_{8h}
        // -------------
        restriction(r_4h, f_8h);
        direct_solver(f_8h, u_8h);

        // On Omega_{4h}
        // -------------
        updateProlongation(u_8h, u_4h);
        jacobi_step(omega, f_4h, u_4h, tmp_4h);
        jacobi_step(omega, f_4h, tmp_4h, u_4h);

        // On Omega_{2h}
        // -------------
        updateProlongation(u_4h, u_2h);
        jacobi_step(omega, f_2h, u_2h, tmp_2h);
        jacobi_step(omega, f_2h, tmp_2h, u_2h);


        // On Omega_{h}
        // ------------
        updateProlongation(u_2h, u_h);
        jacobi_step(omega, f_h, u_h, tmp_h);
        jacobi_step(omega, f_h, tmp_h, u_h);

        evalResidual(u_h, f_h, r_h);
        std::stringstream ss;
        ss << "r_{h} after " << (it+1) << " V-cycle(s)";
        print_record(r_h, ss.str().c_str());
    }

    evalError(u, u_h, e_h);
    print_record(e_h, "Error e_h = u - u_h after last V-cycle");
    print_record(u_h, "u_h after last V-cycle");
}

theon$ g++ -Wall -I /home/numerik/pub/hpc/ws18/session25/ -std=c++17 -o mg_with_4_levels mg_with_4_levels.cpp
theon$ ./mg_with_4_levels > mg_with_4_levels.data
theon$ 

Plotting results

theon$ gnuplot mg_with_4_levels.gnuplot
theon$ gnuplot mg_with_4_levels2.gnuplot
theon$ gnuplot mg_with_4_levels3.gnuplot
theon$ gnuplot mg_with_4_levels4.gnuplot
theon$ gnuplot mg_with_4_levels5.gnuplot
theon$ gnuplot mg_with_4_levels6.gnuplot
theon$ gnuplot mg_with_4_levels7.gnuplot
theon$ 

Plot

Home     Back     Next     1    2    3    4    5 Top of current page     Browse Files