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

using namespace AMDiS;

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

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