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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
struct Dissipation {
scalar dis;
vector u;
face vector mu;
#if COMPRESSIBLE
face vector lambdav;
#endif
};
#if AXI
# define gamma (1./y)
#else
# define gamma (1.)
#endif
void dissipation (struct Dissipation p)
{
vector u = p.u;
scalar dis = p.dis;
(const) face vector mu = p.mu;
#if COMPRESSIBLE
(const) face vector lambdav = p.lambdav;
#endif
#if AXI
scalar ur = u.y;
#endif
#if TREE
/* conservative coarse/fine discretisation (2nd order) */
foreach () {
#if AXI
dis[] = ((mu.x[] + mu.x[1] + mu.y[] + mu.y[0,1])/2.*gamma*sq(ur[]/y)
#if COMPRESSIBLE
+ (lambdav.x[] + lambdav.x[1] +lambdav.y[]+lambdav.y[0,1])/4.*
(ur[]/y + (u.x[1] -u.x[-1] + u.y[0,1] -u.y[0,-1])/2./Delta)*ur[]/y
#endif
);
#else
dis[] = 0.;
#endif
}
foreach_dimension() {
face vector taux[];
#if COMPRESSIBLE
face vector tauc[];
#if AXI
face vector axic[];
#endif
#endif
foreach_face(x) {
#if COMPRESSIBLE
#if AXI
axic.x[] = lambdav.x[]*(ur[]+ur[-1])
*(u.x[] - u.x[-1])/Delta/2.;
#endif
tauc.x[] = lambdav.x[]*(u.x[] - u.x[-1]
#if dimension > 1
+ (u.y[0,1] + u.y[-1,1])/4
- (u.y[0,-1] + u.y[-1,-1])/4.
#endif
#if dimension > 2
+ (u.z[0,0,1] + u.z[-1,0,1])/4
- (u.z[0,0,-1] + u.z[-1,0,-1])/4.
#endif
)*(u.x[] - u.x[-1])/sq(Delta);
#endif
taux.x[] = 2.*mu.x[]*sq((u.x[] - u.x[-1])/Delta);
}
#if dimension > 1
foreach_face(y)
taux.y[] = mu.y[]*(u.x[] - u.x[0,-1] +
(u.y[1,-1] + u.y[1,0])/4. -
(u.y[-1,-1] + u.y[-1,0])/4.)
*(u.x[] - u.x[0,-1])/sq(Delta);
#endif
#if dimension > 2
foreach_face(z)
taux.z[] = mu.z[]*(u.x[] - u.x[0,0,-1] +
(u.z[1,0,-1] + u.z[1,0,0])/4. -
(u.z[-1,0,-1] + u.z[-1,0,0])/4.)
*(u.x[] - u.x[0,0,-1])/sq(Delta);
#endif
foreach () {
double d = 0;
foreach_dimension()
d += (taux.x[1] + taux.x[])*gamma;
dis[] += (d
#if COMPRESSIBLE
+ tauc.x[1] + tauc.x[]
#if AXI
+ (axic.x[1] + axic.x[])/y
#endif
#endif
)/2.;
}
}
#else
/* /\* "naive" discretisation (only 1st order on trees) *\/ */
foreach () {
#if AXI
dis[] = ((mu.x[] + mu.x[1] + mu.y[] + mu.y[0,1])/2.*gamma*sq(ur[]/y)
#if COMPRESSIBLE
+ (lambdav.x[] + lambdav.x[1] +lambdav.y[]+lambdav.y[0,1])/4.*
(ur[]/y + (u.x[1] -u.x[-1] + u.y[0,1] -u.y[0,-1])/2./Delta)*ur[]/y
#endif
);
#else
dis[] = 0.;
#endif
foreach_dimension()
dis[] += ((mu.x[1,0]*sq(u.x[1] - u.x[])
+ mu.x[]*sq(u.x[] - u.x[-1])
#if dimension > 1
+ mu.y[0,1]*(u.x[0,1] - u.x[] +
(u.y[1,0] + u.y[1,1])/4. -
(u.y[-1,0] + u.y[-1,1])/4.)*(u.x[0,1] - u.x[])/2.
+ mu.y[]*(u.x[] - u.x[0,-1] +
(u.y[1,-1] + u.y[1,0])/4. -
(u.y[-1,-1] + u.y[-1,0])/4.)*(u.x[] - u.x[0,-1])/2.
#endif
#if dimension > 2
+ mu.z[0,0,1]*(u.x[0,0,1] - u.x[] +
(u.z[1,0,0] + u.z[1,0,1])/4. -
(u.z[-1,0,0] + u.z[-1,0,1])/4.)*(u.x[0,0,1] - u.x[])/2.
- mu.z[]*(u.x[] - u.x[0,0,-1] +
(u.z[1,0,-1] + u.z[1,0,0])/4. -
(u.z[-1,0,-1] + u.z[-1,0,0])/4.)*(u.x[] - u.x[0,0,-1])/2.
#endif
)*gamma
#if COMPRESSIBLE
+ lambdav.x[1]*sq(u.x[1] - u.x[])/2.
+ lambdav.x[]*sq(u.x[] - u.x[-1])/2.
#if dimension > 1
+ lambdav.x[1]*((u.y[1,1] + u.y[0,1])/4 -
(u.y[1,-1] + u.y[0,-1])/4.)*(u.x[1] - u.x[])/2.
+ lambdav.x[]*((u.y[0,1] + u.y[-1,1])/4 -
(u.y[0,-1] + u.y[-1,-1])/4.)*(u.x[] - u.x[-1])/2.
#if AXI
+ lambdav.x[1]*(ur[1] + ur[])*(u.x[1] - u.x[])/4./y*Delta
+ lambdav.x[]*(ur[-1] + ur[])*(u.x[] - u.x[-1])/4./y*Delta
#endif
#endif
#if dimension > 2
+ lambdav.x[1]*((u.z[1,0,1] + u.z[0,0,1])/4 -
(u.z[1,0,-1] + u.z[0,0,-1])/4.)*(u.x[1] - u.x[])/2.
+ lambdav.x[]*((u.z[0,0,1] + u.z[-1,0,1])/4 -
(u.z[0,0,-1] + u.z[-1,0,-1])/4.)*(u.x[] - u.x[-1])/2.
#endif
#endif
)/sq(Delta);
}
#endif
}
#undef gamma
|
