#include #include #include #include #include #include #include #include #include //------------------------------------------------------------------------------ using namespace hpc; using namespace hpc::matvec; template const T & assertEqual(const T &a, const T &b) { assert(a==b); return a; } //------------------------------------------------------------------------------ template < typename F, typename G, typename T, template class VectorF, template class VectorU, Require< Dense>, Dense> > = true > void discretizeProblem(const F &f, const G &g, VectorF &f_h, VectorU &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 class VectorU, template class VectorE, Require< Dense>, Dense> > = true > void evalError(const U &u, const VectorU &u_h, VectorE &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 class VectorU, template class VectorF, template class VectorR, Require< Dense>, Dense>, Dense> > = true > void evalResidual(const VectorU &u_h, const VectorF &f_h, VectorR &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 class VectorR1, template class VectorR2, Require< Dense>, Dense> > = true > void restriction(const VectorR1 &r_h, VectorR2 &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 class VectorU1, template class VectorU2, Require< Dense>, Dense> > = true > void updateProlongation(const VectorU1 &u_2h, VectorU2 &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 class Vector, Require< Dense> > = true > void addNoise(std::size_t k, Vector &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 class Vector, Require< Dense> > = true > T lInfNorm(const Vector &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 class Vector1, template class Vector2, Require< Dense>, Dense> > = true > T lInfNorm(const Vector1 &x1_h, const Vector2 &x2_h) { auto N = assertEqual(x1_h.length(), x2_h.length()); T res = 0; for (std::size_t i=0; i res) { res = std::abs(x1_h(i) - x2_h(i)); } } return res; } //------------------------------------------------------------------------------ template < typename T, template class Vector, Require< Dense> > = true > T l1Norm(const Vector &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 class VectorF, template class VectorU, template class VectorS, Require< Dense>, Dense>, Dense> > = true > void jacobi_step(double omega, const VectorF &f_h, const VectorU &u_h, VectorS &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 class VectorF, template class VectorU, Require< Dense>, Dense> > = true > void direct_solver(const VectorF &f_h, VectorU &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=1; --k) { u_h(k) = (u_h(k)+u_h(k+1))*k/(k+1); } } //------------------------------------------------------------------------------ template class Vector> void init_zero(Vector &u_h) { for (auto [i, ui] : u_h) { i = i; ui = T(0); } } //------------------------------------------------------------------------------ template class Vector> void print_record(const Vector &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 -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; DenseVector f_h(N0+1), u_h(N0+1), tmp_h(N0+1), r_h(N0+1); DenseVector f_2h(N1+1), u_2h(N1+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, "Resifual r_h of initial guess"); for (std::size_t it=1; it<=3; ++it) { // On Omega_{h} // ------------ jacobi_step(omega, f_h, u_h, tmp_h); jacobi_step(omega, f_h, tmp_h, u_h); // On Omega_{2h} // ------------- restriction(r_h, f_2h); direct_solver(f_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); // usual we would compute the norm of r_h and stop if is below a // given threshold (and continue otherwise) std::stringstream ss; ss << "r_{h} after " << it << " iteration(s)"; print_record(r_h, ss.str().c_str()); } print_record(u_h, "u_h"); }