Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/AUSM.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Advection Upstream Splitting Method (AUSM+) Riemann flux function
9 : : \details This file implements the Advection Upstream Splitting Method (AUSM)
10 : : Riemann solver, with the all-speed corrections.
11 : : Ref. Liou, M. S. (2006). A sequel to AUSM, Part II: AUSM+-up for
12 : : all speeds. Journal of computational physics, 214(1), 137-170.
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef AUSM_h
16 : : #define AUSM_h
17 : :
18 : : #include <vector>
19 : :
20 : : #include "Types.hpp"
21 : : #include "Fields.hpp"
22 : : #include "Tags.hpp"
23 : : #include "FunctionPrototypes.hpp"
24 : : #include "Inciter/Options/Flux.hpp"
25 : : #include "EoS/EoS.hpp"
26 : : #include "MultiMat/MultiMatIndexing.hpp"
27 : :
28 : : namespace inciter {
29 : :
30 : : //! AUSM+up approximate Riemann solver
31 : : //! \details This class can be used polymorphically with inciter::RiemannSolver
32 : : struct AUSM {
33 : :
34 : : //! AUSM+up approximate Riemann solver flux function
35 : : //! \param[in] fn Face/Surface normal
36 : : //! \param[in] u Left and right unknown/state vector
37 : : //! \return Riemann flux solution according to AUSM+up, appended by Riemann
38 : : //! velocities and volume-fractions.
39 : : //! \note The function signature must follow tk::RiemannFluxFn
40 : : static tk::RiemannFluxFn::result_type
41 : 7366575 : flux( const std::array< tk::real, 3 >& fn,
42 : : const std::array< std::vector< tk::real >, 2 >& u,
43 : : const std::vector< std::array< tk::real, 3 > >& = {} )
44 : : {
45 : : const auto nmat =
46 : 7366575 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[0];
47 : :
48 : 7366575 : auto ncomp = u[0].size()-(3+nmat);
49 : 7366575 : std::vector< tk::real > flx( ncomp, 0 );
50 : :
51 : : // Primitive variables
52 : : tk::real rhol(0.0), rhor(0.0);
53 [ + + ]: 24661575 : for (std::size_t k=0; k<nmat; ++k)
54 : : {
55 : 17295000 : rhol += u[0][densityIdx(nmat, k)];
56 : 17295000 : rhor += u[1][densityIdx(nmat, k)];
57 : : }
58 : :
59 : : tk::real pl(0.0), pr(0.0), amatl(0.0), amatr(0.0);
60 [ + - ][ + - ]: 7366575 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
[ - - ][ - - ]
61 [ + - ][ + - ]: 7366575 : hml(nmat, 0.0), hmr(nmat, 0.0),
[ - - ][ - - ]
62 [ + - ][ + - ]: 7366575 : pml(nmat, 0.0), pmr(nmat, 0.0),
[ - - ][ - - ]
63 [ + - ][ - - ]: 7366575 : arhom12(nmat, 0.0),
64 [ + - ][ - - ]: 7366575 : amat12(nmat, 0.0);
65 [ + + ]: 24661575 : for (std::size_t k=0; k<nmat; ++k)
66 : : {
67 [ + - ]: 17295000 : al_l[k] = u[0][volfracIdx(nmat, k)];
68 : 17295000 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
69 [ + - ]: 17295000 : pl += pml[k];
70 : 17295000 : hml[k] = u[0][energyIdx(nmat, k)] + pml[k];
71 [ + - ]: 17295000 : amatl = eos_soundspeed< tag::multimat >( 0,
72 : : u[0][densityIdx(nmat, k)],
73 : : pml[k], al_l[k], k );
74 : :
75 [ + - ]: 17295000 : al_r[k] = u[1][volfracIdx(nmat, k)];
76 : 17295000 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
77 : 17295000 : pr += pmr[k];
78 : 17295000 : hmr[k] = u[1][energyIdx(nmat, k)] + pmr[k];
79 [ + - ]: 17295000 : amatr = eos_soundspeed< tag::multimat >( 0,
80 : : u[1][densityIdx(nmat, k)],
81 : : pmr[k], al_r[k], k );
82 : :
83 : : // Average states for mixture speed of sound
84 : 17295000 : arhom12[k] = 0.5*(u[0][densityIdx(nmat, k)] + u[1][densityIdx(nmat, k)]);
85 : 17295000 : amat12[k] = 0.5*(amatl+amatr);
86 : : }
87 : :
88 : 7366575 : auto rho12 = 0.5*(rhol+rhor);
89 : :
90 : : // mixture speed of sound
91 : : tk::real ac12(0.0);
92 [ + + ]: 24661575 : for (std::size_t k=0; k<nmat; ++k)
93 : : {
94 : 17295000 : ac12 += (arhom12[k]*amat12[k]*amat12[k]);
95 : : }
96 : 7366575 : ac12 = std::sqrt( ac12/rho12 );
97 : :
98 : : // Independently limited velocities for advection
99 : 7366575 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
100 : 7366575 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
101 : 7366575 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
102 : 7366575 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
103 : 7366575 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
104 : 7366575 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
105 : :
106 : : // Face-normal velocities from advective velocities
107 : 7366575 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
108 : 7366575 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
109 : :
110 : : // Mach numbers
111 : 7366575 : auto ml = vnl/ac12;
112 : 7366575 : auto mr = vnr/ac12;
113 : :
114 : : // All-speed parameters
115 : : // These parameters control the amount of all-speed diffusion necessary for
116 : : // low-Mach flows. Setting k_u and k_p to zero does not add any all-speed
117 : : // diffusion, whereas setting k_u and k_p to 1 adds maximum recommended
118 : : // all-speed diffusion. See "Liou, M. S. (2006). A sequel to AUSM, Part II:
119 : : // AUSM+-up for all speeds. Journal of computational physics, 214(1),
120 : : // 137-170" for more mathematical explanation. k_u is the velocity diffusion
121 : : // term and k_p is the pressure diffusion term. These two terms reduce
122 : : // pressure-velocity decoupling (chequerboarding/odd-even oscillations).
123 : : tk::real k_u(1.0), k_p(1.0), f_a(1.0);
124 : :
125 : : // Split Mach polynomials
126 : 7366575 : auto msl = splitmach_ausm( f_a, ml );
127 : 7366575 : auto msr = splitmach_ausm( f_a, mr );
128 : :
129 : : // Riemann Mach number
130 : 7366575 : auto m0 = 1.0 - (0.5*(vnl*vnl + vnr*vnr)/(ac12*ac12));
131 [ + + ]: 7366575 : auto mp = -k_p* std::max(m0, 0.0) * (pr-pl) / (f_a*rho12*ac12*ac12);
132 : 7366575 : auto m12 = msl[0] + msr[1] + mp;
133 : 7366575 : auto vriem = ac12 * m12;
134 : :
135 : : // Riemann pressure
136 : 7366575 : auto pu = -k_u* msl[2] * msr[3] * f_a * rho12 * ac12 * (vnr-vnl);
137 : 7366575 : auto p12 = msl[2]*pl + msr[3]*pr + pu;
138 : :
139 : : // Flux vector splitting
140 : 7366575 : auto l_plus = 0.5 * (vriem + std::fabs(vriem));
141 : 7366575 : auto l_minus = 0.5 * (vriem - std::fabs(vriem));
142 : :
143 : : // Conservative fluxes
144 [ + + ]: 24661575 : for (std::size_t k=0; k<nmat; ++k)
145 : : {
146 : 17295000 : flx[volfracIdx(nmat, k)] = l_plus*al_l[k] + l_minus*al_r[k];
147 : 17295000 : flx[densityIdx(nmat, k)] = l_plus*u[0][densityIdx(nmat, k)]
148 : 17295000 : + l_minus*u[1][densityIdx(nmat, k)];
149 : 17295000 : flx[energyIdx(nmat, k)] = l_plus*hml[k] + l_minus*hmr[k];
150 : : }
151 : :
152 [ + + ]: 29466300 : for (std::size_t idir=0; idir<3; ++idir)
153 : : {
154 : 22099725 : flx[momentumIdx(nmat, idir)] = l_plus*u[0][momentumIdx(nmat, idir)]
155 : 22099725 : + l_minus*u[1][momentumIdx(nmat, idir)]
156 : 22099725 : + p12*fn[idir];
157 : : }
158 : :
159 : 7366575 : l_plus = l_plus/( std::fabs(vriem) + 1.0e-16 );
160 : 7366575 : l_minus = l_minus/( std::fabs(vriem) + 1.0e-16 );
161 : :
162 : : // Store Riemann-advected partial pressures
163 [ + + ]: 7366575 : if (std::fabs(l_plus) > 1.0e-10)
164 : : {
165 [ + + ]: 8250369 : for (std::size_t k=0; k<nmat; ++k)
166 [ + - ]: 5818771 : flx.push_back( pml[k] );
167 : : }
168 [ + + ]: 4934977 : else if (std::fabs(l_minus) > 1.0e-10)
169 : : {
170 [ + + ]: 8307320 : for (std::size_t k=0; k<nmat; ++k)
171 [ + - ]: 5865218 : flx.push_back( pmr[k] );
172 : : }
173 : : else
174 : : {
175 [ + + ]: 8103886 : for (std::size_t k=0; k<nmat; ++k)
176 [ + - ][ - - ]: 5611011 : flx.push_back( 0.5*(pml[k] + pmr[k]) );
177 : : }
178 : :
179 : : // Store Riemann velocity
180 [ + - ]: 7366575 : flx.push_back( vriem );
181 : :
182 : : Assert( flx.size() == (3*nmat+3+nmat+1), "Size of multi-material flux "
183 : : "vector incorrect" );
184 : :
185 : 7366575 : return flx;
186 : : }
187 : :
188 : : //! Flux type accessor
189 : : //! \return Flux type
190 : : static ctr::FluxType type() noexcept { return ctr::FluxType::AUSM; }
191 : :
192 : : private:
193 : :
194 : : //! Split Mach polynomials for AUSM+ flux
195 : : //! \param[in] fa All-speed parameter
196 : : //! \param[in] mach Local Mach numner
197 : : //! \return Values of the positive and negative split Mach and pressure
198 : : //! polynomials.
199 : : //! \details This function returns a vector with positive and negative Mach
200 : : //! and pressure polynomials, as:
201 : : //! ms[0] = M_4(+),
202 : : //! ms[1] = M_4(-),
203 : : //! ms[2] = P_5(+), and
204 : : //! ms[3] = P_5(-).
205 : : //! For more details, ref. Liou, M. S. (2006). A sequel to AUSM, Part II:
206 : : //! AUSM+-up for all speeds. J. Comp. Phys., 214(1), 137-170.
207 : 14733150 : static std::array< tk::real, 4 > splitmach_ausm( tk::real fa,
208 : : tk::real mach )
209 : : {
210 : : std::array< tk::real, 4 > ms;
211 : :
212 : : std::array< tk::real, 3 > msplus, msminus;
213 : : tk::real psplus, psminus;
214 : :
215 : 14733150 : msplus[0] = 0.5*(mach + std::fabs(mach));
216 : 14733150 : msminus[0]= 0.5*(mach - std::fabs(mach));
217 : :
218 : 14733150 : msplus[1] = +0.25*(mach + 1.0)*(mach + 1.0);
219 : 14733150 : msminus[1]= -0.25*(mach - 1.0)*(mach - 1.0);
220 : :
221 : 14733150 : auto alph_fa = (3.0/16.0) * (-4.0 + 5.0*fa*fa);
222 : :
223 [ + + ]: 14733150 : if (std::fabs(mach) >= 1.0)
224 : : {
225 : : msplus[2] = msplus[0];
226 : : msminus[2]= msminus[0];
227 : 168 : psplus = msplus[0]/mach;
228 : 168 : psminus = msminus[0]/mach;
229 : : }
230 : : else
231 : : {
232 : 14732982 : msplus[2] = msplus[1]* (1.0 - 2.0*msminus[1]);
233 : 14732982 : msminus[2]= msminus[1]* (1.0 + 2.0*msplus[1]);
234 : 14732982 : psplus = msplus[1]*
235 : 14732982 : ((+2.0 - mach) - (16.0 * alph_fa)*mach*msminus[1]);
236 : 14732982 : psminus = msminus[1]*
237 : 14732982 : ((-2.0 - mach) + (16.0 * alph_fa)*mach*msplus[1]);
238 : : }
239 : :
240 : 14733150 : ms[0] = msplus[2];
241 : 14733150 : ms[1] = msminus[2];
242 : 14733150 : ms[2] = psplus;
243 : 14733150 : ms[3] = psminus;
244 : :
245 : 14733150 : return ms;
246 : : }
247 : : };
248 : :
249 : : } // inciter::
250 : :
251 : : #endif // AUSM_h
|
