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
|
/**
* @brief Return the flux vector for the APE equations.
*
* @param physfield Fields.
* @param flux Resulting flux. flux[eq][dir][pt]
*/
void APE::v_GetFluxVector(
const Array<OneD, Array<OneD, NekDouble>> &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &flux)
{
int nq = physfield[0].num_elements();
Array<OneD, NekDouble> tmp1(nq);
Array<OneD, NekDouble> tmp2(nq);
ASSERTL1(flux[0].num_elements() == m_spacedim,
"Dimension of flux array and velocity array do not match");
// F_{adv,p',j} = bar{rho} bar{c^2} u'_j + p' bar{u}_j
for (int j = 0; j < m_spacedim; ++j)
{
Vmath::Zero(nq, flux[0][j], 1);
// construct bar{rho} bar{c^2} u'_j
Vmath::Vmul(nq, m_bf[0], 1, m_bf[1], 1, tmp1, 1);
Vmath::Vmul(nq, tmp1, 1, physfield[j + 1], 1, tmp1, 1);
// construct p' bar{u}_j term
Vmath::Vmul(nq, physfield[0], 1, m_bf[j + 2], 1, tmp2, 1);
// add both terms
Vmath::Vadd(nq, tmp1, 1, tmp2, 1, flux[0][j], 1);
}
for (int i = 1; i < flux.num_elements(); ++i)
{
ASSERTL1(flux[i].num_elements() == m_spacedim,
"Dimension of flux array and velocity array do not match");
// F_{adv,u'_i,j} = (p'/ bar{rho} + bar{u}_k u'_k) delta_{ij}
for (int j = 0; j < m_spacedim; ++j)
{
Vmath::Zero(nq, flux[i][j], 1);
if (i - 1 == j)
{
// contruct p'/ bar{rho} term
Vmath::Vdiv(nq, physfield[0], 1, m_bf[1], 1, flux[i][j], 1);
// construct bar{u}_k u'_k term
Vmath::Zero(nq, tmp1, 1);
for (int k = 0; k < m_spacedim; ++k)
{
Vmath::Vvtvp(nq, physfield[k + 1], 1, m_bf[k + 2], 1, tmp1,
1, tmp1, 1);
}
// add terms
Vmath::Vadd(nq, flux[i][j], 1, tmp1, 1, flux[i][j], 1);
}
}
}
}
|
近期评论