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 "Fields.hpp"
21 : : #include "FunctionPrototypes.hpp"
22 : : #include "Inciter/Options/Flux.hpp"
23 : : #include "SplitMachFns.hpp"
24 : : #include "EoS/EOS.hpp"
25 : : #include "MultiMat/MultiMatIndexing.hpp"
26 : :
27 : : namespace inciter {
28 : :
29 : : //! AUSM+up approximate Riemann solver
30 : : struct AUSM {
31 : :
32 : : //! AUSM+up approximate Riemann solver flux function
33 : : //! \param[in] fn Face/Surface normal
34 : : //! \param[in] u Left and right unknown/state vector
35 : : //! \return Riemann flux solution according to AUSM+up, appended by Riemann
36 : : //! velocities and volume-fractions.
37 : : //! \note The function signature must follow tk::RiemannFluxFn
38 : : static tk::RiemannFluxFn::result_type
39 : 6608351 : flux( const std::vector< EOS >& mat_blk,
40 : : const std::array< tk::real, 3 >& fn,
41 : : const std::array< std::vector< tk::real >, 2 >& u,
42 : : const std::vector< std::array< tk::real, 3 > >& = {} )
43 : : {
44 : 6608351 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
45 : :
46 : : // All-speed parameters
47 : : // These parameters control the amount of all-speed diffusion necessary for
48 : : // low-Mach flows. Setting k_u and k_p to zero does not add any all-speed
49 : : // diffusion, whereas setting k_u and k_p to 1 adds maximum recommended
50 : : // all-speed diffusion. See "Liou, M. S. (2006). A sequel to AUSM, Part II:
51 : : // AUSM+-up for all speeds. Journal of computational physics, 214(1),
52 : : // 137-170" for more mathematical explanation. k_u is the velocity diffusion
53 : : // term and k_p is the pressure diffusion term. These two terms reduce
54 : : // pressure-velocity decoupling (chequerboarding/odd-even oscillations).
55 : 6608351 : auto k_u = g_inputdeck.get< tag::lowspeed_ku >();
56 : 6608351 : auto k_p = g_inputdeck.get< tag::lowspeed_kp >();
57 : :
58 : 6608351 : auto ncomp = u[0].size()-(3+nmat);
59 : 6608351 : std::vector< tk::real > flx( ncomp, 0 );
60 : :
61 : : // Primitive variables
62 : : tk::real rhol(0.0), rhor(0.0);
63 [ + + ]: 22520119 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 : 15911768 : rhol += u[0][densityIdx(nmat, k)];
66 : 15911768 : rhor += u[1][densityIdx(nmat, k)];
67 : : }
68 : :
69 : : tk::real pl(0.0), pr(0.0), amatl(0.0), amatr(0.0);
70 [ + - ][ + - ]: 6608351 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
[ - - ][ - - ]
71 [ + - ][ + - ]: 6608351 : hml(nmat, 0.0), hmr(nmat, 0.0),
[ - - ][ - - ]
72 [ + - ][ + - ]: 6608351 : pml(nmat, 0.0), pmr(nmat, 0.0),
[ - - ][ - - ]
73 [ + - ][ - - ]: 6608351 : arhom12(nmat, 0.0),
74 [ + - ][ - - ]: 6608351 : amat12(nmat, 0.0);
75 [ + + ]: 22520119 : for (std::size_t k=0; k<nmat; ++k)
76 : : {
77 [ + - ]: 15911768 : al_l[k] = u[0][volfracIdx(nmat, k)];
78 : 15911768 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
79 [ + - ]: 15911768 : pl += pml[k];
80 : 15911768 : hml[k] = u[0][energyIdx(nmat, k)] + pml[k];
81 [ + - ]: 15911768 : amatl = mat_blk[k].compute< EOS::soundspeed >(
82 : : u[0][densityIdx(nmat, k)], pml[k], al_l[k], k );
83 : :
84 [ + - ]: 15911768 : al_r[k] = u[1][volfracIdx(nmat, k)];
85 : 15911768 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
86 [ + - ]: 15911768 : pr += pmr[k];
87 : 15911768 : hmr[k] = u[1][energyIdx(nmat, k)] + pmr[k];
88 [ + - ]: 15911768 : amatr = mat_blk[k].compute< EOS::soundspeed >(
89 : : u[1][densityIdx(nmat, k)], pmr[k], al_r[k], k );
90 : :
91 : : // Average states for mixture speed of sound
92 : 15911768 : arhom12[k] = 0.5*(u[0][densityIdx(nmat, k)] + u[1][densityIdx(nmat, k)]);
93 : 15911768 : amat12[k] = 0.5*(amatl+amatr);
94 : : }
95 : :
96 : 6608351 : auto rho12 = 0.5*(rhol+rhor);
97 : :
98 : : // mixture speed of sound
99 : : tk::real ac12(0.0);
100 [ + + ]: 22520119 : for (std::size_t k=0; k<nmat; ++k)
101 : : {
102 : 15911768 : ac12 += (arhom12[k]*amat12[k]*amat12[k]);
103 : : }
104 : 6608351 : ac12 = std::sqrt( ac12/rho12 );
105 : :
106 : : // Independently limited velocities for advection
107 : 6608351 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
108 : 6608351 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
109 : 6608351 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
110 : 6608351 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
111 : 6608351 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
112 : 6608351 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
113 : :
114 : : // Face-normal velocities from advective velocities
115 : 6608351 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
116 : 6608351 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
117 : :
118 : : // Mach numbers
119 : 6608351 : auto ml = vnl/ac12;
120 : 6608351 : auto mr = vnr/ac12;
121 : :
122 : : tk::real f_a(1.0);
123 : :
124 : : // Split Mach polynomials
125 : 6608351 : auto msl = splitmach_ausm( ml, f_a );
126 : 6608351 : auto msr = splitmach_ausm( mr, f_a );
127 : :
128 : : // Riemann Mach number
129 : 6608351 : auto m0 = 1.0 - (0.5*(vnl*vnl + vnr*vnr)/(ac12*ac12));
130 [ - + ]: 6608351 : auto mp = -k_p* std::max(m0, 0.0) * (pr-pl) / (f_a*rho12*ac12*ac12);
131 : 6608351 : auto m12 = msl[0] + msr[1] + mp;
132 : 6608351 : auto vriem = ac12 * m12;
133 : :
134 : : // Riemann pressure
135 : 6608351 : auto pu = -k_u* msl[2] * msr[3] * f_a * rho12 * ac12 * (vnr-vnl);
136 : 6608351 : auto p12 = msl[2]*pl + msr[3]*pr + pu;
137 : :
138 : : auto md = 0.0;
139 : : auto delta = 0.0;
140 : : // Velocity magnitudes
141 : : auto vmag_l = tk::dot( {{ul, vl, wl}}, {{ul, vl, wl}} );
142 : : auto vmag_r = tk::dot( {{ur, vr, wr}}, {{ur, vr, wr}} );
143 : 6608351 : auto m0_mod = 4.0
144 : 6608351 : * (0.25 - (0.5*(vmag_l*vmag_l + vmag_r*vmag_r)/(ac12*ac12)));
145 : : //// uncomment code below AND set k_u and k_p to zero for AUSM-2025u/p mods.
146 : : //// Additional diffusion
147 : : //delta = 4.0;
148 : : //md = std::max(m0_mod, 0.0) * delta * std::sqrt(std::abs(vnl - vnr) * ac12);
149 : : ////md = std::max(m0_mod, 0.0) * delta * std::sqrt(std::abs(pl - pr) / rho12);
150 : :
151 : : // Flux vector splitting
152 : 6608351 : auto l_plus = 0.5 * (vriem + std::fabs(vriem) + 2.0*md);
153 : 6608351 : auto l_minus = 0.5 * (vriem - std::fabs(vriem) - 2.0*md);
154 : :
155 : : // Conservative fluxes
156 [ + + ]: 22520119 : for (std::size_t k=0; k<nmat; ++k)
157 : : {
158 : 15911768 : flx[volfracIdx(nmat, k)] = l_plus*al_l[k] + l_minus*al_r[k];
159 : 15911768 : flx[densityIdx(nmat, k)] = l_plus*u[0][densityIdx(nmat, k)]
160 : 15911768 : + l_minus*u[1][densityIdx(nmat, k)];
161 : 15911768 : flx[energyIdx(nmat, k)] = l_plus*hml[k] + l_minus*hmr[k];
162 : : }
163 : :
164 [ + + ]: 26433404 : for (std::size_t idir=0; idir<3; ++idir)
165 : : {
166 : 19825053 : flx[momentumIdx(nmat, idir)] = l_plus*u[0][momentumIdx(nmat, idir)]
167 : 19825053 : + l_minus*u[1][momentumIdx(nmat, idir)]
168 : 19825053 : + p12*fn[idir];
169 : : }
170 : :
171 : : // Evaluate pressure work biasing in energy equation
172 : 13216702 : l_plus = 0.5 * (vriem + std::fabs(vriem))
173 : 6608351 : / ( vriem + std::copysign(1.0e-12,vriem) )
174 [ + + ]: 6608351 : + delta*std::max(m0_mod, 0.0);
175 : 13216702 : l_minus = 0.5 * (vriem - std::fabs(vriem))
176 : 6608351 : / ( vriem + std::copysign(1.0e-12,vriem) )
177 [ + + ]: 7561255 : - delta*std::max(m0_mod, 0.0);
178 : :
179 : : // Store Riemann-advected partial pressures
180 [ + + ]: 22520119 : for (std::size_t k=0; k<nmat; ++k)
181 [ + - ][ - - ]: 15911768 : flx.push_back( l_plus*pml[k] + l_minus*pmr[k] );
182 : :
183 : : // Store Riemann velocity
184 [ + - ]: 6608351 : flx.push_back( vriem );
185 : :
186 : : Assert( flx.size() == (3*nmat+3+nmat+1), "Size of multi-material flux "
187 : : "vector incorrect" );
188 : :
189 : 6608351 : return flx;
190 : : }
191 : :
192 : : //! Flux type accessor
193 : : //! \return Flux type
194 : : static ctr::FluxType type() noexcept { return ctr::FluxType::AUSM; }
195 : : };
196 : :
197 : : } // inciter::
198 : :
199 : : #endif // AUSM_h
|
