stokes0.cc 2.54 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
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

#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

using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;

int main(int argc, char** argv)
{
    AMDiS::init(argc, argv);

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

    double viscosity = 1.0;
    Parameters::get("stokes->viscosity", viscosity);

    // tree-paths for components
37
38
    auto _v = Dune::Indices::_0;
    auto _p = Dune::Indices::_1;
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

    // <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));

      // <d_i(v_i), p>
      auto opP = makeOperator(tag::partialtest_trial{i}, 1.0);
      prob.addMatrixOperator(opP, treepath(_v,i), _p);

      // <q, d_i(u_i)>
      auto opDiv = makeOperator(tag::test_partialtrial{i}, 1.0);
      prob.addMatrixOperator(opDiv, _p, treepath(_v,i));
    }

54
55
    auto opZero = makeOperator(tag::test_trial{}, 0.0);
    prob.addMatrixOperator(opZero, _p, _p);
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

    // 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; };

    // define boundary values
    auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,DOW>
    {
      return {0.0, x[1]*(1.0 - x[1])};
    };

    auto zero = [](auto const& x) -> Dune::FieldVector<double,DOW>
    {
      return {0.0, 0.0};
    };

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

    prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);

78
    AdaptInfo adaptInfo("adapt");
79
80

    // assemble and solve system
81
    prob.buildAfterAdapt(adaptInfo, Flag(0));
82

83
#ifdef DEBUG_MTL
84
85
    // write matrix to file
    mtl::io::matrix_market_ostream out("matrix_stokes0.mtx");
86
87
    out << prob.getSystemMatrix().matrix();
    std::cout << prob.getSystemMatrix().matrix() << '\n';
88
#endif
89

90
91
92
93
94
95
96
97
    prob.solve(adaptInfo);

    // output solution
    prob.writeFiles(adaptInfo);

    AMDiS::finalize();
    return 0;
}