navier_stokes.cc 2.57 KB
Newer Older
1
2
#include <amdis/AMDiS.hpp>
#include <amdis/AdaptInstationary.hpp>
3
#include <amdis/LocalOperators.hpp>
4
5
#include <amdis/ProblemStat.hpp>
#include <amdis/ProblemInstat.hpp>
6
#include <amdis/localoperators/StokesOperator.hpp>
7

8
9
using namespace AMDiS;

Praetorius, Simon's avatar
Praetorius, Simon committed
10
using Grid = Dune::YaspGrid<2>;
11
using Basis = TaylorHoodBasis<Grid>;
12

13
14
int main(int argc, char** argv)
{
15
  Environment env(argc, argv);
16

Praetorius, Simon's avatar
Praetorius, Simon committed
17
  ProblemStat<Basis> prob("stokes");
18
19
  prob.initialize(INIT_ALL);

Praetorius, Simon's avatar
Praetorius, Simon committed
20
  ProblemInstat<Basis> probInstat("stokes", prob);
21
22
23
24
25
26
27
28
29
  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);

30
  AdaptInfo adaptInfo("adapt");
31
32

  // tree-paths for components
33
34
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
35

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

38
  // <1/tau * u, v>
Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
  auto opTime = makeOperator(tag::testvec_trialvec{},
    density * invTau);
41
42
43
  prob.addMatrixOperator(opTime, _v, _v);

  // <1/tau * u^old, v>
Praetorius, Simon's avatar
Praetorius, Simon committed
44
45
  auto opTimeOld = makeOperator(tag::testvec{},
    density * invTau * prob.solution(_v));
46
47
48
  prob.addVectorOperator(opTimeOld, _v);


49
50
51
  // 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());
52

53
  // <(u * nabla)u_i^old, v_i>
Praetorius, Simon's avatar
Praetorius, Simon committed
54
55
  auto opNonlin1 = makeOperator(tag::testvec_trialvec{},
    density * trans(gradientAtQP(prob.solution(_v))));
56
  prob.addMatrixOperator(opNonlin1, _v, _v);
57

58
  for (std::size_t i = 0; i < WORLDDIM; ++i) {
59
    // <(u^old * nabla)u_i, v_i>
Praetorius, Simon's avatar
Praetorius, Simon committed
60
61
    auto opNonlin2 = makeOperator(tag::test_gradtrial{},
      density * prob.solution(_v));
62
63
    prob.addMatrixOperator(opNonlin2, treepath(_v,i), treepath(_v,i));
  }
64

65
  // <(u^old * grad(u_i^old)), v_i>
Praetorius, Simon's avatar
Praetorius, Simon committed
66
67
  auto opNonlin3 = makeOperator(tag::testvec{},
    trans(gradientAtQP(prob.solution(_v))) * prob.solution(_v));
68
69
  prob.addVectorOperator(opNonlin3, _v);

70
  // define boundary values
Praetorius, Simon's avatar
Praetorius, Simon committed
71
  auto parabolic_y = [](auto const& x)
72
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
73
    return FieldVector<double,2>{0.0, x[1]*(1.0 - x[1])};
74
  };
75

Praetorius, Simon's avatar
Praetorius, Simon committed
76
  FieldVector<double,2> zero(0);
77

78
  // set boundary conditions for velocity
Praetorius, Simon's avatar
Praetorius, Simon committed
79
80
81
  prob.boundaryManager()->setBoxBoundary({1,2,2,2});
  prob.addDirichletBC(1, _v, _v, parabolic_y);
  prob.addDirichletBC(2, _v, _v, zero);
82

83
  // set initial conditions
Praetorius, Simon's avatar
Praetorius, Simon committed
84
  prob.solution(_v).interpolate(parabolic_y);
85
86
87
88
89
90
91
92

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

  // output solution
  prob.writeFiles(adaptInfo);

  return 0;
93
}