Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

stokes1.cc 2.38 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
10
11
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.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
}