/
formulae_mpdata_hot_2d.hpp
173 lines (162 loc) · 5.75 KB
/
formulae_mpdata_hot_2d.hpp
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
161
162
163
164
165
166
167
168
169
170
171
172
173
/** @file
* @copyright University of Warsaw
* @section LICENSE
* GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
*/
#pragma once
#include <libmpdata++/formulae/mpdata/formulae_mpdata_common.hpp>
#include <libmpdata++/formulae/mpdata/formulae_mpdata_psi_2d.hpp>
#include <libmpdata++/formulae/mpdata/formulae_mpdata_gc_2d.hpp>
#include <libmpdata++/formulae/mpdata/formulae_mpdata_g_2d.hpp>
#include <boost/preprocessor/punctuation/comma.hpp>
namespace libmpdataxx
{
namespace formulae
{
namespace mpdata
{
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto ndxx_psi_coeff(
const arr_2d_t &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j
)
{
return return_helper<ix_t>(
(
3 * GC(pi<dim>(i+h, j)) * abs(GC(pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
- 2 * pow3(GC(pi<dim>(i+h, j))) / pow2(G_bar_x<opts, dim>(G, i, j))
- GC(pi<dim>(i+h, j))
) / 6
);
}
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto ndxy_psi_coeff(
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j
)
{
return return_helper<ix_t>(
(abs(GC[dim](pi<dim>(i+h, j))) - 2 * pow2(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j))
* GC1_bar_xy<dim>(GC[dim+1], i, j) / (2 * G_bar_x<opts, dim>(G, i, j))
);
}
// third order terms
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto TOT(
const arr_2d_t &psi,
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j,
typename std::enable_if<opts::isset(opts, opts::tot)>::type* = 0
)
{
return return_helper<ix_t>(
ndxx_psi<opts, dim>(psi, i, j) * ndxx_psi_coeff<opts, dim>(GC[dim], G, i, j)
+
ndxy_psi<opts, dim>(psi, i, j) * ndxy_psi_coeff<opts, dim>(GC, G, i, j)
);
}
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto TOT(
const arr_2d_t &psi,
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j,
typename std::enable_if<!opts::isset(opts, opts::tot)>::type* = 0
)
{
return 0;
}
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto ndxxx_psi_coeff(
const arr_2d_t &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j
)
{
return return_helper<ix_t>(
(
6 * abs(GC(pi<dim>(i+h, j))) * pow2(GC(pi<dim>(i+h, j))) / pow2(G_bar_x<opts, dim>(G, i, j))
- 3 * pow2(GC(pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
- 3 * pow4(GC(pi<dim>(i+h, j))) / pow3(G_bar_x<opts, dim>(G, i, j))
) / 2
);
}
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto ndxxy_psi_coeff(
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j
)
{
return return_helper<ix_t>(
3 *
(abs(GC[dim](pi<dim>(i+h, j))) - pow2(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j))
* GC[dim](pi<dim>(i+h, j)) * GC1_bar_xy<dim>(GC[dim-1], i, j)
/ pow2(G_bar_x<opts, dim>(G, i, j))
);
}
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto ndxyy_psi_coeff(
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j
)
{
return return_helper<ix_t>(
3 *
(
3 * abs(GC1_bar_xy<dim>(GC[dim-1], i, j)) * pow2(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
+ 3 * pow2(GC1_bar_xy<dim>(GC[dim-1], i, j)) * abs(GC[dim](pi<dim>(i+h, j))) / G_bar_x<opts, dim>(G, i, j)
- 2 * abs(GC1_bar_xy<dim>(GC[dim-1], i, j)) * abs(GC[dim](pi<dim>(i+h, j)))
- 9 * pow2(GC1_bar_xy<dim>(GC[dim-1], i, j)) * pow2(GC[dim](pi<dim>(i+h, j)))
/ pow2(G_bar_x<opts, dim>(G, i, j))
) / G_bar_x<opts, dim>(G, i, j) / 4
);
}
// fourth order terms
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto FOT(
const arr_2d_t &psi,
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j,
typename std::enable_if<opts::isset(opts, opts::fot)>::type* = 0
)
{
static_assert(opts::isset(opts, opts::tot) || opts::isset(opts, opts::div_3rd),
"adding fourth-order terms makes sense only when third-order terms are present (tot or div_3rd option)");
static_assert(opts::isset(opts, opts::iga), "fot option only available with iga");
return return_helper<ix_t>(
ndxxx_psi<opts, dim>(psi, i, j) * ndxxx_psi_coeff<opts, dim>(GC[dim], G, i, j)
+
ndxxy_psi<opts, dim>(psi, i, j) * ndxxy_psi_coeff<opts, dim>(GC, G, i, j)
+
ndxyy_psi<opts, dim>(psi, i, j) * ndxyy_psi_coeff<opts, dim>(GC, G, i, j)
);
}
template<opts_t opts, int dim, class arr_2d_t, class ix_t>
forceinline_macro auto FOT(
const arr_2d_t &psi,
const arrvec_t<arr_2d_t> &GC,
const arr_2d_t &G,
const ix_t &i,
const ix_t &j,
typename std::enable_if<!opts::isset(opts, opts::fot)>::type* = 0
)
{
return 0;
}
} // namespace mpdata
} // namespace formulae
} // namespcae libmpdataxx