navier_stokes.cc 3.86 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   = DefaultTraitsMesh<Mesh, AMDIS_DIM, AMDIS_DOW, 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
50
    // define the differential operators
    using Op = StokesProblem::OperatorType;
51
52
    For<0,DOW>::loop([&](auto const _i)
    {
53
	// <1/tau * u_i, v_i>
54
55
56
57
58
59
        auto opTime = Op::zot( density );
        auto opTimeOld = Op::zot( valueOf(prob.getSolution(_i), density) );
        prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau());
        prob.addVectorOperator(opTimeOld,  _i, probInstat.getInvTau());

#if STRATEGY == 1
60
61
62
        // <(u * nabla)u_i^old, v_i>
        For<0, DOW>::loop([&](auto const _j)
        {
63
          auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) );
64
65
          prob.addMatrixOperator(opNonlin, _i, _j);
        });
66
#elif STRATEGY == 2
67
68
69
        // <(u^old * nabla)u_i, v_i>
        For<0, DOW>::loop([&](auto const _j)
        {
70
71
72
73
74
75
76
77
78
          auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
          prob.addMatrixOperator(opNonlin, _i, _i);
        });
#else
        // <(u^old * nabla)u_i^old, v_i>
        For<0, DOW>::loop([&](auto const _j)
        {
          auto opNonlin = Op::zot( density * valueOf(prob.getSolution(_j)) * derivativeOf(prob.getSolution(_i), _j) );
          prob.addVectorOperator(opNonlin, _i);
79
        });
80
81
82
83
84
85
86
87
88
89
90
#endif

	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>

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

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


97
98
99
    // 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; };
100
101
    auto corner = [](auto const& x) { return x[0] == 0.0 && x[1] == 0.0; };

102
103
104
    // 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]); };
105

106
107
108
    // set boundary conditions
    for (size_t i = 0; i < DOW; ++i)
	prob.addDirichletBC(box, i, i, zero);
109

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

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

116
    *prob.getSolution() = 0.0;
117

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

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

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