navier_stokes.cc 3.03 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <iostream>
#include <ctime>
#include <cmath>

#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>

#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>;

int main(int argc, char** argv)
{
    AMDiS::init(argc, argv);
    
    StokesProblem prob("stokes");
    prob.initialize(INIT_ALL);
    
    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);
    
    
    AdaptInfo adaptInfo("adapt");
    double tau = adaptInfo.getTimestep();
    
    // define the differential operators
    using Op = StokesProblem::OperatorType;
    For<0,DOW>::loop([&](auto const _i) 
    {    
	// <1/tau * u_i, v_i> + <viscosity*grad(u_i), grad(v_i)>
	Op* opL = new Op;
	opL->addZOT( constant(density/tau) );
	opL->addSOT( constant(viscosity) ); 
	prob.addMatrixOperator(*opL, _i, _i);
	
	// <1/tau * u_i^old, v_i>
	Op* opRhs = new Op;
	opRhs->addZOT( valueOf(prob.getSolution(_i), density/tau) );
	prob.addVectorOperator(*opRhs, _i);
	
	// <p, d_i(v_i)>
	Op* opB = new Op;
	opB->addFOT( constant(1.0), _i, GRD_PSI ); 
	prob.addMatrixOperator(*opB, _i, DOW); 
    
	// <d_i(u_i), q>
	Op* opDiv = new Op; 
	opDiv->addFOT( constant(1.0), _i, GRD_PHI );
	prob.addMatrixOperator(*opDiv, DOW, _i); 
	
	// <(u^old * nabla)u_i, v_i>
	// NOTE: only simples linearization of nonlin-term, since gradientOf() not yet implemented
	For<0, DOW>::loop([&](auto const _j)
	{
	    Op* opNonlin = new Op;
	    opNonlin->addFOT( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
	    prob.addMatrixOperator(*opNonlin, _i, _i);
	});
    });
    
    // 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);
    
    // set initial values for the solver (maybe not necessary)
    prob.getSolution().getVector() = 0.0;
    
    double t = 0.0;
    for (t = adaptInfo.getStartTime();
	 t < adaptInfo.getEndTime();
	 t+= adaptInfo.getTimestep()) 
    {
	AMDIS_MSG("Timestep t = " << t);
	prob.writeFiles(adaptInfo, t);
    
	// assemble and solve system
	prob.buildAfterCoarsen(adaptInfo, Flag(0));
	prob.solve(adaptInfo);
    }
    
    // output solution
    prob.writeFiles(adaptInfo, t);
    
    AMDiS::finalize();
    return 0;
}