navier_stokes.cc 3.94 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
24
using Mesh = Dune::YaspGrid<AMDIS_DIM>;
// using Mesh = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;

25
// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
26
using StokesParam   = LagrangeTraits<Mesh, 2, 2, 1>;
27
using StokesProblem = ProblemStat<StokesParam>;
28
using StokesProblemInstat = ProblemInstat<StokesParam>;
29
30
31
32

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

34
35
    StokesProblem prob("stokes");
    prob.initialize(INIT_ALL);
36

37
38
    StokesProblemInstat probInstat("stokes", prob);
    probInstat.initialize(INIT_UH_OLD);
39

40
41
42
43
44
45
    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);
46

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

49
    // define the differential operators
50
    using Op = StokesProblem::ElementOperator;
51
    forEach(range_<0,DOW>, [&](auto const _i)
52
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
53
	      // <1/tau * u_i, v_i>
54
55
        auto opTime = Op::zot( density );
        auto opTimeOld = Op::zot( valueOf(prob.getSolution(_i), density) );
Praetorius, Simon's avatar
Praetorius, Simon committed
56
57
        prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau());
        prob.addVectorOperator(opTimeOld,  _i, probInstat.getInvTau());
58

59
        // <(u * nabla)u_i^old, v_i>
60
        forEach(range_<0, DOW>, [&](auto const _j)
61
        {
62
          auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) );
Praetorius, Simon's avatar
Praetorius, Simon committed
63
          prob.addMatrixOperator(opNonlin, _i, _j);
64
        });
Praetorius, Simon's avatar
Praetorius, Simon committed
65

66
        // <(u^old * nabla)u_i, v_i>
67
        forEach(range_<0, DOW>, [&](auto const _j)
68
        {
69
          auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
Praetorius, Simon's avatar
Praetorius, Simon committed
70
          prob.addMatrixOperator(opNonlin, _i, _i);
71
        });
Praetorius, Simon's avatar
Praetorius, Simon committed
72

73
        // <(u^old * nabla)u_i^old, v_i>
74
        forEach(range_<0, DOW>, [&](auto const _j)
75
        {
Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
78
79
80
          auto opNonlin = Op::zot( func([density](double uj, double dj_ui)
            {
              return density * uj * dj_ui;
            }, valueOf(prob.getSolution(_j)), derivativeOf(prob.getSolution(_i), _j)) );
          prob.addVectorOperator(opNonlin, _i);
81
        });
82

Praetorius, Simon's avatar
Praetorius, Simon committed
83
84
85
        auto opL = std::make_pair( Op::sot( viscosity ), false ); // <viscosity*grad(u_i), grad(v_i)>
        auto opP = std::make_pair( Op::fot( 1.0, _i, GRD_PSI ), false); // <p, d_i(v_i)>
        auto opDiv = std::make_pair( Op::fot( 1.0, _i, GRD_PHI ), false); // <d_i(u_i), q>
86
87
88
89
90
91

        const int i = _i;
        prob.addMatrixOperator({{
          {{i,  i}, opL}, {{i,DOW}, opP},
          {{DOW,i}, opDiv}
        }});
92
    });
93
94
95
96
97

    auto opZero = std::make_pair( Op::zot( 0.0 ), false );
    prob.addMatrixOperator({{ {{DOW,  DOW}, opZero} }});


98
99
100
    // define boundary regions
    auto left = [](auto const& x) { return x[0] < 1.e-8; };
    auto box  = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };
101
102
    auto corner = [](auto const& x) { return x[0] == 0.0 && x[1] == 0.0; };

103
104
105
    // define boundary values
    auto zero      = [](auto const& x) { return 0.0; };
    auto parabolic = [vel](auto const& x) { return vel*4.0*x[1]*(1.0 - x[1]); };
106

107
108
    // set boundary conditions
    for (size_t i = 0; i < DOW; ++i)
Praetorius, Simon's avatar
Praetorius, Simon committed
109
	    prob.addDirichletBC(box, i, i, zero);
110

111
112
    prob.addDirichletBC(left, 0, 0, zero);
    prob.addDirichletBC(left, 1, 1, parabolic);
113
114
115
116

    // pressure fixation
    prob.addDirichletBC(corner, DOW, DOW, zero);

117
    *prob.getSolution() = 0.0;
118

119
120
    AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
    adapt.adapt();
121

122
    // output solution
123
    prob.writeFiles(adaptInfo);
124

125
126
127
    AMDiS::finalize();
    return 0;
}