navier_stokes.cc 3.49 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

10
11
#ifdef DOW
#undef DOW
12
#endif
13
#define DOW AMDIS_DOW
14
15
16
17

using namespace AMDiS;

// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
18
using StokesParam   = ProblemStatTraits<AMDIS_DIM, AMDIS_DOW, 2, 2, 1>;
19
using StokesProblem = ProblemStat<StokesParam>;
20
using StokesProblemInstat = ProblemInstat<StokesParam>;
21
22
23
24
25
26
27

int main(int argc, char** argv)
{
    AMDiS::init(argc, argv);
    
    StokesProblem prob("stokes");
    prob.initialize(INIT_ALL);
28
29
30
  
    StokesProblemInstat probInstat("stokes", prob);
    probInstat.initialize(INIT_UH_OLD);
31
32
33
34
35
36
37
38
    
    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);
    
39
40
41
42
    using DOFVectorU = std::decay_t< decltype(prob.getSolution<0>()) >;
    using DOFVectorP = std::decay_t< decltype(prob.getSolution<DOW>()) >;
    std::array<DOFVectorU*, DOW> velocity;
    For<0,DOW>::loop([&](auto const _i) { velocity[_i] = &prob.getSolution(_i); });
43
44
45
46
47
48
49
    
    AdaptInfo adaptInfo("adapt");
    
    // define the differential operators
    using Op = StokesProblem::OperatorType;
    For<0,DOW>::loop([&](auto const _i) 
    {    
50
	// <1/tau * u_i, v_i>
51
	Op *opTime = new Op, *opTimeOld = new Op;
52
53
	opTime->addZOT( density );
        opTimeOld->addZOT( valueOf(prob.getSolution(_i), density) );
54
55
        prob.addMatrixOperator(*opTime, _i, _i, probInstat.getInvTau());
        prob.addVectorOperator(*opTimeOld,  _i, probInstat.getInvTau());
56
57
58
59
60
        
#if 0
        // <(u * nabla)u_i^old, v_i>
        For<0, DOW>::loop([&](auto const _j)
        {
61
          auto opNonlin = new Op;
62
63
64
65
66
67
68
          opNonlin->addZOT( derivativeOf(prob.getSolution(_i), _j, density) );
          prob.addMatrixOperator(opNonlin, _i, _j);
        });
#else   
        // <(u^old * nabla)u_i, v_i>
        For<0, DOW>::loop([&](auto const _j)
        {
69
          Op *opNonlin = new Op;
70
          opNonlin->addFOT( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
71
          prob.addMatrixOperator(*opNonlin, _i, _i);
72
73
74
        });
#endif      
        // <viscosity*grad(u_i), grad(v_i)>
75
        Op *opL = new Op;
76
	opL->addSOT( viscosity ); 
77
	prob.addMatrixOperator(*opL, _i, _i);
78
79
	
	// <p, d_i(v_i)>
80
	Op *opP = new Op;
81
	opP->addFOT( 1.0, _i, GRD_PSI ); 
82
	prob.addMatrixOperator(*opP, _i, DOW); 
83
84
    
	// <d_i(u_i), q>
85
	Op *opDiv = new Op; 
86
	opDiv->addFOT( 1.0, _i, GRD_PHI );
87
	prob.addMatrixOperator(*opDiv, DOW, _i);
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
    });
    
    // 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; };
    
    // 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]); };
    
    // set boundary conditions
    for (size_t i = 0; i < DOW; ++i)
	prob.addDirichletBC(box, i, i, zero);
    
    prob.addDirichletBC(left, 0, 0, zero);
    prob.addDirichletBC(left, 1, 1, parabolic);
    
105
    *prob.getSolution() = 0.0;
106
    
107
108
    AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
    adapt.adapt();
109
110
    
    // output solution
111
    prob.writeFiles(adaptInfo);
112
113
114
115
    
    AMDiS::finalize();
    return 0;
}