navier_stokes.cc 3.29 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/common/FieldMatVec.hpp>
11
12
13
14
15
#include <amdis/AdaptInstationary.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/ProblemInstat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/assembler/StokesOperator.hpp>
16

17
18
using namespace AMDiS;

19
20
struct NavierStokesBasis
{
21
  using Grid = Dune::YaspGrid<AMDIS_DIM>;
22
23
  using GlobalBasis = typename TaylorHoodBasis<Grid::LeafGridView>::GlobalBasis;
};
24

25
26
using StokesProblem = ProblemStat<NavierStokesBasis>;
using StokesProblemInstat = ProblemInstat<NavierStokesBasis>;
27

28
29
int main(int argc, char** argv)
{
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
  AMDiS::init(argc, argv);

  StokesProblem prob("stokes");
  prob.initialize(INIT_ALL);

  StokesProblemInstat probInstat("stokes", prob);
  probInstat.initialize(INIT_UH_OLD);

  double viscosity = 1.0;
  double density = 1.0;
  double vel = 1.0;
  Parameters::get("stokes->viscosity", viscosity);
  Parameters::get("stokes->density", density);
  Parameters::get("stokes->boundary velocity", vel);

45
  AdaptInfo adaptInfo("adapt");
46
47

  // tree-paths for components
48
49
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
50

Praetorius, Simon's avatar
Praetorius, Simon committed
51
52
  auto invTau = std::ref(probInstat.invTau());

53
  // <1/tau * u, v>
Praetorius, Simon's avatar
Praetorius, Simon committed
54
55
  auto opTime = makeOperator(tag::testvec_trialvec{},
    density * invTau);
56
57
58
  prob.addMatrixOperator(opTime, _v, _v);

  // <1/tau * u^old, v>
Praetorius, Simon's avatar
Praetorius, Simon committed
59
60
  auto opTimeOld = makeOperator(tag::testvec{},
    density * invTau * prob.solution(_v));
61
62
63
  prob.addVectorOperator(opTimeOld, _v);


64
65
66
  // sum_i <grad(u_i),grad(v_i)> + <p, div(v)> + <div(u), q>
  auto opStokes = makeOperator(tag::stokes{}, viscosity);
  prob.addMatrixOperator(opStokes, treepath(), treepath());
67

68
  // <(u * nabla)u_i^old, v_i>
Praetorius, Simon's avatar
Praetorius, Simon committed
69
70
  auto opNonlin1 = makeOperator(tag::testvec_trialvec{},
    density * trans(gradientAtQP(prob.solution(_v))));
71
  prob.addMatrixOperator(opNonlin1, _v, _v);
72

73
  for (std::size_t i = 0; i < AMDIS_DOW; ++i) {
74
    // <(u^old * nabla)u_i, v_i>
Praetorius, Simon's avatar
Praetorius, Simon committed
75
76
    auto opNonlin2 = makeOperator(tag::test_gradtrial{},
      density * prob.solution(_v));
77
78
    prob.addMatrixOperator(opNonlin2, treepath(_v,i), treepath(_v,i));
  }
79

80
  // <(u^old * grad(u_i^old)), v_i>
Praetorius, Simon's avatar
Praetorius, Simon committed
81
82
  auto opNonlin3 = makeOperator(tag::testvec{},
    trans(gradientAtQP(prob.solution(_v))) * prob.solution(_v));
83
84
  prob.addVectorOperator(opNonlin3, _v);

85
86
87
  // 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; };
88

89
90
91
92
93
  // define boundary values
  auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,2>
  {
    return {0.0, x[1]*(1.0 - x[1])};
  };
94

95
96
97
98
  auto zero = [](auto const& x) -> Dune::FieldVector<double,2>
  {
    return {0.0, 0.0};
  };
99

100
101
102
  // set boundary conditions for velocity
  prob.addDirichletBC(left, _v, _v, parabolic_y);
  prob.addDirichletBC(not_left, _v, _v, zero);
103

104
105
106
  // set point constraint for pressure
  prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);

107
  // set initial conditions
Praetorius, Simon's avatar
Praetorius, Simon committed
108
109
  prob.solution(_v).interpolate(parabolic_y);
  prob.solution(_p).interpolate(0.0);
110
111
112
113
114
115
116
117
118

  AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
  adapt.adapt();

  // output solution
  prob.writeFiles(adaptInfo);

  AMDiS::finalize();
  return 0;
119
}