stokes1.cc 1.52 KB
Newer Older
1
#include <amdis/AMDiS.hpp>
2
#include <amdis/LocalOperators.hpp>
3
#include <amdis/ProblemStat.hpp>
4
5
6

using namespace AMDiS;

Praetorius, Simon's avatar
Praetorius, Simon committed
7
using Grid = Dune::YaspGrid<2>;
8
using StokesParam = TaylorHoodBasis<Grid>;
9
10
11
12
using StokesProblem = ProblemStat<StokesParam>;

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

Praetorius, Simon's avatar
Praetorius, Simon committed
15
16
  StokesProblem prob("stokes");
  prob.initialize(INIT_ALL);
17

Praetorius, Simon's avatar
Praetorius, Simon committed
18
19
  double viscosity = 1.0;
  Parameters::get("stokes->viscosity", viscosity);
20

Praetorius, Simon's avatar
Praetorius, Simon committed
21
22
23
  // tree-paths for components
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
24

Praetorius, Simon's avatar
Praetorius, Simon committed
25
26
  for (std::size_t i = 0; i < WORLDDIM; ++i) {
    // <viscosity*grad(u_i), grad(v_i)>
Praetorius, Simon's avatar
Praetorius, Simon committed
27
28
29
    auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
    prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
  }
30

Praetorius, Simon's avatar
Praetorius, Simon committed
31
32
33
  // <d_i(v_i), p>
  auto opP = makeOperator(tag::divtestvec_trial{}, 1.0);
  prob.addMatrixOperator(opP, _v, _p);
34

Praetorius, Simon's avatar
Praetorius, Simon committed
35
36
37
  // <q, d_i(u_i)>
  auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0);
  prob.addMatrixOperator(opDiv, _p, _v);
38

Praetorius, Simon's avatar
Praetorius, Simon committed
39
  // define boundary values
Praetorius, Simon's avatar
Praetorius, Simon committed
40
  auto parabolic_y = [](auto const& x)
Praetorius, Simon's avatar
Praetorius, Simon committed
41
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
42
    return FieldVector<double,2>{0.0, x[1]*(1.0 - x[1])};
Praetorius, Simon's avatar
Praetorius, Simon committed
43
  };
44

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

Praetorius, Simon's avatar
Praetorius, Simon committed
47
  // set boundary conditions for velocity
Praetorius, Simon's avatar
Praetorius, Simon committed
48
49
50
  prob.boundaryManager()->setBoxBoundary({1,2,2,2});
  prob.addDirichletBC(1, _v, _v, parabolic_y);
  prob.addDirichletBC(2, _v, _v, zero);
51

Praetorius, Simon's avatar
Praetorius, Simon committed
52
  AdaptInfo adaptInfo("adapt");
53

Praetorius, Simon's avatar
Praetorius, Simon committed
54
55
56
  // assemble and solve system
  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);
57

Praetorius, Simon's avatar
Praetorius, Simon committed
58
59
  // output solution
  prob.writeFiles(adaptInfo);
60

Praetorius, Simon's avatar
Praetorius, Simon committed
61
62
  AMDiS::finalize();
  return 0;
63
}