0%

AcousticSolver求解器代码阅读

发表于 分类于 Valine：
本文字数： 2.8k 阅读时长 ≈ 5 分钟

EquationSystems->AcousticSystem.cpp

  • AcousticSystem::v_InitObject()
    m_bfNames背景流场变量名
    1
    2
    3
    4
    5
    m_bfNames.push_back("c0sq");
    m_bfNames.push_back("rho0");
    m_bfNames.push_back("u0");
    m_bfNames.push_back("v0");
    m_bfNames.push_back("w0");

EquationSystems->APE.cpp

下面一段代码求解了APE方程的通量项，发现这里求连续方程通量是把c0sq直接乘进去了，而不是求完通量后再乘以c0sq，这意味着求解器只能求c0sq空间导数很小的流动。好心累。。。为什么发布的版本这里都不该，ppt里说是改了。。。

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);
            }
        }
    }
}