stokes1.cc 2.39 KB
Newer Older
1
2
3
4
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

5
6
7
8
#include <iostream>
#include <ctime>
#include <cmath>

9
#include <amdis/AMDiS.hpp>
10
#include <amdis/LocalOperators.hpp>
11
#include <amdis/ProblemStat.hpp>
12
13
14
15
16
17
18
19
20
21

#ifdef DOW
#undef DOW
#endif
#define DOW AMDIS_DOW

using namespace AMDiS;

// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1

22
using Grid = Dune::YaspGrid<AMDIS_DIM>;
23
24
25
26
27
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;

int main(int argc, char** argv)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
28
  AMDiS::init(argc, argv);
29

Praetorius, Simon's avatar
Praetorius, Simon committed
30
31
  StokesProblem prob("stokes");
  prob.initialize(INIT_ALL);
32

Praetorius, Simon's avatar
Praetorius, Simon committed
33
34
  double viscosity = 1.0;
  Parameters::get("stokes->viscosity", viscosity);
35

Praetorius, Simon's avatar
Praetorius, Simon committed
36
37
38
  // tree-paths for components
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
39

Praetorius, Simon's avatar
Praetorius, Simon committed
40
41
42
43
44
  // <viscosity*grad(u_i), grad(v_i)>
  for (std::size_t i = 0; i < DOW; ++i) {
    auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
    prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
  }
45

Praetorius, Simon's avatar
Praetorius, Simon committed
46
47
48
  // <d_i(v_i), p>
  auto opP = makeOperator(tag::divtestvec_trial{}, 1.0);
  prob.addMatrixOperator(opP, _v, _p);
49

Praetorius, Simon's avatar
Praetorius, Simon committed
50
51
52
  // <q, d_i(u_i)>
  auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0);
  prob.addMatrixOperator(opDiv, _p, _v);
53

Praetorius, Simon's avatar
Praetorius, Simon committed
54
55
  auto opZero = makeOperator(tag::test_trial{}, 0.0);
  prob.addMatrixOperator(opZero, _p, _p);
56
57


Praetorius, Simon's avatar
Praetorius, Simon committed
58
59
60
  // define boundary regions
  auto left = [](auto const& x) { return x[0] < 1.e-8; };
  auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };
61

Praetorius, Simon's avatar
Praetorius, Simon committed
62
63
64
65
66
  // define boundary values
  auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,DOW>
  {
    return {0.0, x[1]*(1.0 - x[1])};
  };
67

Praetorius, Simon's avatar
Praetorius, Simon committed
68
69
70
71
  auto zero = [](auto const& x) -> Dune::FieldVector<double,DOW>
  {
    return {0.0, 0.0};
  };
72

Praetorius, Simon's avatar
Praetorius, Simon committed
73
74
75
  // set boundary conditions for velocity
  prob.addDirichletBC(left, _v, _v, parabolic_y);
  prob.addDirichletBC(not_left, _v, _v, zero);
76

Praetorius, Simon's avatar
Praetorius, Simon committed
77
  prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);
78

Praetorius, Simon's avatar
Praetorius, Simon committed
79
  AdaptInfo adaptInfo("adapt");
80

Praetorius, Simon's avatar
Praetorius, Simon committed
81
82
  // assemble and solve system
  prob.assemble(adaptInfo);
83

84
#ifdef DEBUG_MTL
Praetorius, Simon's avatar
Praetorius, Simon committed
85
86
87
88
  // write matrix to file
  mtl::io::matrix_market_ostream out("matrix_stokes1.mtx");
  out << prob.systemMatrix().matrix();
  std::cout << prob.systemMatrix().matrix() << '\n';
89
#endif
90

Praetorius, Simon's avatar
Praetorius, Simon committed
91
  prob.solve(adaptInfo);
92

Praetorius, Simon's avatar
Praetorius, Simon committed
93
94
  // output solution
  prob.writeFiles(adaptInfo);
95

Praetorius, Simon's avatar
Praetorius, Simon committed
96
97
  AMDiS::finalize();
  return 0;
98
}