navier_stokes.cc 3.4 KB
Newer Older
1
2
3
4
5
#include <iostream>
#include <ctime>
#include <cmath>

#include <dune/amdis/AMDiS.hpp>
6
#include <dune/amdis/AdaptInstationary.hpp>
7
#include <dune/amdis/ProblemStat.hpp>
8
#include <dune/amdis/ProblemInstat.hpp>
9
#include <dune/amdis/Terms.hpp>
10

11
12
#ifdef DOW
#undef DOW
13
#endif
14
#define DOW AMDIS_DOW
15

16
17
18
19
#ifndef STRATEGY
#define STRATEGY 1
#endif

20
21
using namespace AMDiS;

22
23
using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
24
using StokesProblem = ProblemStat<StokesParam>;
25
using StokesProblemInstat = ProblemInstat<StokesParam>;
26

27
28
using Op = StokesProblem::ElementOperator;

29
30
int main(int argc, char** argv)
{
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
  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);

  AdaptInfo adaptInfo("adapt", prob.getNumComponents());

  // tree-paths for components
  auto _v = 0_c;
  auto _p = 1_c;

  // <1/tau * u, v>
  Op opTime;
  opTime.addZOT(density);
  prob.addMatrixOperator(opTime, _v, _v);

  // <1/tau * u^old, v>
  Op opTimeOld;
  opTimeOld.addZOT(valueOf(prob.getSolution(_v), density));
  prob.addVectorOperator(opTimeOld, _v);

  // sum_i <grad(u_i),grad(v_i)>
  Op opL;
  opL.addSOT(viscosity);
  prob.addMatrixOperator(opL, _v, _v);

  // <p, div(v)>
  Op opP;
  opP.addFOT(1.0, GRD_PSI);
  prob.addMatrixOperator(opP, _v, _p);

  // <div(u), q>
  Op opDiv;
  opDiv.addFOT(1.0, GRD_PHI);
  prob.addMatrixOperator(opDiv, _p, _v);


  for (std::size_t i = 0; i < DOW; ++i) {
    // <(u * nabla)u_i^old, v_i>
    for (std::size_t j = 0; j < DOW; ++j) {
      Op opNonlin1;
      opNonlin1.addZOT(derivativeOf(prob.getSolution(treepath(_v,i)), j, density));
      prob.addMatrixOperator(opNonlin1, treepath(_v,i), treepath(_v,j));
    }

    // <(u^old * nabla)u_i, v_i>
    Op opNonlin2;
    opNonlin2.addFOT(valueOf(prob.getSolution(_v)), GRD_PHI);
    prob.addMatrixOperator(opNonlin2, treepath(_v,i), treepath(_v,i));

    // <(u^old * grad(u_i^old)), v_i>
    Op opNonlin3;
    opNonlin3.addZOT( func([density](Dune::FieldVector<double,DOW> const& u, Dune::FieldVector<double,DOW> const& grad_ui)
    {
      return density * (u * grad_ui);
    }, valueOf(prob.getSolution(_v)), gradientOf(prob.getSolution(treepath(_v,i)))) );
    prob.addVectorOperator(opNonlin, treepath(_v,i));
  }
99

100
101
102
  // 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; };
103

104
105
106
107
108
  // define boundary values
  auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,2>
  {
    return {0.0, x[1]*(1.0 - x[1])};
  };
109

110
111
112
113
  auto zero = [](auto const& x) -> Dune::FieldVector<double,2>
  {
    return {0.0, 0.0};
  };
114

115
116
117
  // set boundary conditions for velocity
  prob.addDirichletBC(left, _v, _v, parabolic_y);
  prob.addDirichletBC(not_left, _v, _v, zero);
118

119
120
121
122
123
124
125
126
127
128
129
130
  // 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);


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

  // output solution
  prob.writeFiles(adaptInfo);

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