# Smoothing property of the Jacobi method

#### Content

Complete the implementation of the (weighted) Jacobi iteration.

## Code for exercise

#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

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

set terminal svg size 900, 500
set output "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 "smooth.data" index  10 using 1:2 with linespoints lt 2 lw 3, \
"smooth.data" index  20 using 1:2 with linespoints lt 4 lw 3

## How to compile and plot results

theon$g++ -Wall -I /home/numerik/pub/hpc/ws18/session25/ -std=c++17 -o direct_solver -o smooth smooth.cpp smooth.cpp: In instantiation of 'void jacobi_step(double, const VectorF&, const VectorU&, VectorS&) [with T = double; VectorF = hpc::matvec::DenseVector; VectorU = hpc::matvec::DenseVector; VectorS = hpc::matvec::DenseVector; typename std::common_type >::type::is_DenseVector:: value, bool>::type, typename std::enable_if >::type::is_DenseVector:: value, bool>::type, typename std::enable_if >::type::is_DenseVector:: value, bool>::type>::type = 1]': smooth.cpp:218:43: required from here smooth.cpp:178:10: warning: unused variable 'h2' [-Wunused-variable] auto h2 = h*h; ^~ theon$ ./smooth > smooth.data
theon\$

gives