Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

stokes1.cc 2.5 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

#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

22
using Grid = Dune::YaspGrid<AMDIS_DIM>;
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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
    // <viscosity*grad(u_i), grad(v_i)>
41
42
    for (std::size_t i = 0; i < DOW; ++i) {
      auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
43
44
45
      prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
    }

46
47
48
    // <d_i(v_i), p>
    auto opP = makeOperator(tag::divtestvec_trial{}, 1.0);
    prob.addMatrixOperator(opP, _v, _p);
49

50
51
    // <q, d_i(u_i)>
    auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0);
52
53
    prob.addMatrixOperator(opDiv, _p, _v);

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


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

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

79
    AdaptInfo adaptInfo("adapt");
80
81

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

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

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

    // output solution
    prob.writeFiles(adaptInfo);

    AMDiS::finalize();
    return 0;
}