stokes0.cc 2.48 KB
Newer Older
1
2
3
4
#include <iostream>
#include <ctime>
#include <cmath>

5
6
7
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
8
9
10
11
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
37
38
39
40
41
42
43
44
45
46
47
48
49

#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
    auto _v = 0_c;
    auto _p = 1_c;

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

50
51
    auto opZero = makeOperator(tag::test_trial{}, 0.0);
    prob.addMatrixOperator(opZero, _p, _p);
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

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

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

    // assemble and solve system
    prob.buildAfterCoarsen(adaptInfo, Flag(0));

79
80
81
82
83
    // write matrix to file
    mtl::io::matrix_market_ostream out("matrix_stokes0.mtx");
    out << prob.getSystemMatrix()->getMatrix();
    std::cout << prob.getSystemMatrix()->getMatrix() << '\n';

84
85
86
87
88
89
90
91
    prob.solve(adaptInfo);

    // output solution
    prob.writeFiles(adaptInfo);

    AMDiS::finalize();
    return 0;
}