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
12
13
14
15
16
17
18
19
20
21

#ifndef DIM
#define DIM 2
#endif
#ifndef DOW
#define DOW 2
#endif

using namespace AMDiS;

// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
using StokesParam   = ProblemStatTraits<DIM, DOW, 2, 2, 1>;
using StokesProblem = ProblemStat<StokesParam>;
22
using StokesProblemInstat = ProblemInstat<StokesParam>;
23
24
25
26
27
28
29

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