ellipt.cc 3.37 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>

6
7
#include <fmt/core.h>

8
9
10
11
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/common/Literals.hpp>
12
#include <amdis/gridfunctions/Integrate.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
13

14
using namespace AMDiS;
15
using namespace Dune::Indices;
16

Praetorius, Simon's avatar
Praetorius, Simon committed
17
// 1 component with polynomial degree 1
18
19
using Param   = YaspGridBasis<AMDIS_DIM, 1>;
using ElliptProblem = ProblemStat<Param>;
Praetorius, Simon's avatar
Praetorius, Simon committed
20
21
22
23

int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);
24

25
26
27
28
29
30
31
32
33
34
35
36
37
  int numLevels = AMDIS_DIM == 2 ? 8 : 5;
  if (argc > 2)
    numLevels = std::atoi(argv[2]);

  auto f = [](auto const& x) {
    double r2 = dot(x,x);
    double ux = std::exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * x.size()) * ux;
  };
  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
  auto grad_g = [](auto const& x) -> FieldMatrix<double,1,AMDIS_DIM> {
    return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
  };
38

Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
41
42
43
  using Grid = Dune::YaspGrid<AMDIS_DIM>;
  auto gridPtr = MeshCreator<Grid>::create("mesh");
  auto basis = makeBasis(gridPtr->leafGridView(), lagrange<1>());

  auto prob = makeProblemStat("ellipt", *gridPtr, basis);
Praetorius, Simon's avatar
Praetorius, Simon committed
44
  prob.initialize(INIT_ALL);
45

46
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
Praetorius, Simon's avatar
Praetorius, Simon committed
47
  prob.addMatrixOperator(opL);
48

49
  auto opForce = makeOperator(tag::test{}, f, 6);
50
  prob.addVectorOperator(opForce, _0);
51

Praetorius, Simon's avatar
Praetorius, Simon committed
52
  // set boundary condition
53
54
55
56
  auto boundary = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };
  prob.addDirichletBC(boundary, _0, _0, g);

  AdaptInfo adaptInfo("adapt");
57

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
  std::vector<double> errL2; errL2.reserve(numLevels);
  std::vector<double> errH1; errH1.reserve(numLevels);
  std::vector<double> widths; widths.reserve(numLevels);
  for (int i = 0; i < numLevels; ++i) {
    prob.getGrid()->globalRefine(1);
    auto gridView = prob.getGlobalBasis()->gridView();

    double h = 0;
    for (auto const& e : edges(gridView))
      h = std::max(h, e.geometry().volume());
    widths.push_back(h);

    prob.buildAfterCoarsen(adaptInfo, Flag(0));
    prob.solve(adaptInfo);

    double errorL2 = integrate(sqr(g - prob.getSolution(_0)), gridView, 6);
    errL2.push_back(std::sqrt(errorL2));
    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.getSolution(_0))), gridView, 6);
    errH1.push_back(std::sqrt(errorH1));

#if WRITE_FILES
    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(prob.getGlobalBasis()->gridView());
    vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
    vtkWriter.write("u_" + std::to_string(i));
#endif
  }
84

85
86
87
88
  std::cout << "\n";
  fmt::print("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}\n",
             "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
  fmt::print("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
89

90
91
92
93
  std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
  for (int i = 1; i < numLevels; ++i) {
    eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]);
    eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]);
Praetorius, Simon's avatar
Praetorius, Simon committed
94
  }
95

96
97
98
  for (int i = 0; i < numLevels; ++i)
    fmt::print("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}\n",
                i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
99

Praetorius, Simon's avatar
Praetorius, Simon committed
100
101
  AMDiS::finalize();
  return 0;
102
}