Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/LDFSS.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 Low Diffusion Flux Splitting Scheme (LDFSS+) Riemann flux function
9 : : \details This file implements the modified Low Diffusion Flux Splitting
10 : : Scheme (LDFSS) Riemann solver. See Edwards, J. (2001, June).
11 : : Towards unified CFD simulations of real fluid flows. In 15th AIAA
12 : : computational fluid dynamics conference (p. 2524).
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef LDFSS_h
16 : : #define LDFSS_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 : : //! LDFSS approximate Riemann solver
30 : : struct LDFSS {
31 : :
32 : : //! LDFSS 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 LDFSS, appended by Riemann
36 : : //! velocities and volume-fractions.
37 : : //! \note The function signature must follow tk::RiemannFluxFn
38 : : static tk::RiemannFluxFn::result_type
39 : 0 : 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 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
45 : :
46 : 0 : auto ncomp = u[0].size()-(3+nmat);
47 [ - - ]: 0 : std::vector< tk::real > flx( ncomp, 0 );
48 : :
49 : : // Primitive variables
50 : 0 : tk::real rhol(0.0), rhor(0.0);
51 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
52 : : {
53 : 0 : rhol += u[0][densityIdx(nmat, k)];
54 : 0 : rhor += u[1][densityIdx(nmat, k)];
55 : : }
56 : :
57 : 0 : tk::real pl(0.0), pr(0.0), amatl(0.0), amatr(0.0);
58 [ - - ][ - - ]: 0 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
59 [ - - ][ - - ]: 0 : hml(nmat, 0.0), hmr(nmat, 0.0),
60 [ - - ][ - - ]: 0 : pml(nmat, 0.0), pmr(nmat, 0.0),
61 [ - - ]: 0 : arhom12(nmat, 0.0),
62 [ - - ]: 0 : amat12(nmat, 0.0);
63 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 : 0 : al_l[k] = u[0][volfracIdx(nmat, k)];
66 : 0 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
67 : 0 : pl += pml[k];
68 : 0 : hml[k] = u[0][energyIdx(nmat, k)] + pml[k];
69 [ - - ]: 0 : amatl = mat_blk[k].compute< EOS::soundspeed >(
70 : 0 : u[0][densityIdx(nmat, k)], pml[k], al_l[k], k );
71 : :
72 : 0 : al_r[k] = u[1][volfracIdx(nmat, k)];
73 : 0 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
74 : 0 : pr += pmr[k];
75 : 0 : hmr[k] = u[1][energyIdx(nmat, k)] + pmr[k];
76 [ - - ]: 0 : amatr = mat_blk[k].compute< EOS::soundspeed >(
77 : 0 : u[1][densityIdx(nmat, k)], pmr[k], al_r[k], k );
78 : :
79 : : // Average states for mixture speed of sound
80 : 0 : arhom12[k] = 0.5*(u[0][densityIdx(nmat, k)] + u[1][densityIdx(nmat, k)]);
81 : 0 : amat12[k] = 0.5*(amatl+amatr);
82 : : }
83 : :
84 : 0 : auto rho12 = 0.5*(rhol+rhor);
85 : :
86 : : // mixture speed of sound
87 : 0 : tk::real ac12(0.0);
88 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
89 : : {
90 : 0 : ac12 += (arhom12[k]*amat12[k]*amat12[k]);
91 : : }
92 : 0 : ac12 = std::sqrt( ac12/rho12 );
93 : :
94 : : // Independently limited velocities for advection
95 : 0 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
96 : 0 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
97 : 0 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
98 : 0 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
99 : 0 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
100 : 0 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
101 : :
102 : : // Face-normal velocities from advective velocities
103 : 0 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
104 : 0 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
105 : :
106 : : // Mach numbers
107 : 0 : auto ml = vnl/ac12;
108 : 0 : auto mr = vnr/ac12;
109 : :
110 : : // Split Mach polynomials
111 : 0 : auto msl = splitmach_ausm( ml, 1.0, 0.0, 0.0 );
112 : 0 : auto msr = splitmach_ausm( mr, 1.0, 0.0, 0.0 );
113 : :
114 : : // Modified Mach polynomials
115 : : // mtilde for general-order Mach polys (AUSM+)
116 : 0 : auto mtilde = 0.5*(msl[0] - std::max(0.0,ml) - msr[1] + std::min(0.0,mr));
117 : : //// Uncomment for LDFSS-2025M
118 : : //// Modified mtilde for 2nd order Mach polys (vanLeer)
119 : : //auto beta_l = -std::max(0.0, 1.0 - std::floor(std::abs(ml)));
120 : : //auto beta_r = -std::max(0.0, 1.0 - std::floor(std::abs(mr)));
121 : : //auto mtilde = 0.25*beta_l*beta_r*std::pow((std::pow(0.5*(ml*ml+mr*mr),0.25) - 1.0), 2.0);
122 : :
123 : 0 : auto dp = pl - pr;
124 : 0 : auto delta = 4.0;
125 : :
126 : : // LDFSS-2001
127 : : // Additional diffusion
128 : 0 : auto dpplus = std::abs(dp);
129 : : //dpplus += 8.0*rhol*ac12*std::sqrt(std::abs(vnl - vnr)*ac12); // 2025u mod
130 : : //dpplus += 8.0*rhol*ac12*std::sqrt(std::abs(dp)/rho12); // 2025p mod
131 : 0 : auto dpminus = std::abs(dp);
132 : : //dpminus += 8.0*rhor*ac12*std::sqrt(std::abs(vnl - vnr)*ac12); // 2025u mod
133 : : //dpminus += 8.0*rhor*ac12*std::sqrt(std::abs(dp)/rho12); // 2025p mod
134 : : auto mtilde_plus = mtilde*
135 : 0 : std::max( 0.0, 1.0 - (dp+delta*dpplus)/(2.0*rhol*ac12*ac12) );
136 : : auto mtilde_minus = mtilde*
137 : 0 : std::max( 0.0, 1.0 + (dp-delta*dpminus)/(2.0*rhor*ac12*ac12) );
138 : :
139 : : //// LDFSS(2) with the 2025p mod
140 : : //auto dpplus = std::abs(dp)
141 : : // + 4.0*pl*std::sqrt(2.0*std::abs(dp)/(pl+pr));
142 : : //auto dpminus = std::abs(dp)
143 : : // + 4.0*pr*std::sqrt(2.0*std::abs(dp)/(pl+pr));
144 : : //auto mtilde_plus = mtilde*
145 : : // std::max( 0.0, 1.0 - dp/(pl+pr) - delta*dpplus/pl );
146 : : //auto mtilde_minus = mtilde*
147 : : // std::max( 0.0, 1.0 + dp/(pl+pr) - delta*dpminus/pr );
148 : :
149 : 0 : auto C_plus = msl[0] - mtilde_plus;
150 : 0 : auto C_minus = msr[1] + mtilde_minus;
151 : :
152 : : // Flux vector splitting
153 : 0 : auto l_plus = C_plus * ac12;
154 : 0 : auto l_minus = C_minus * ac12;
155 : 0 : auto vriem = l_plus + l_minus;
156 : :
157 : : // Riemann pressure
158 : 0 : auto p12 = 0.5*(pl+pr) + 0.5*(pl-pr)*(msl[2]-msr[3]) +
159 : 0 : rho12*ac12*ac12*(msl[2]+msr[3]-1.0);
160 : :
161 : : // Conservative fluxes
162 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
163 : : {
164 : 0 : flx[volfracIdx(nmat, k)] = l_plus*al_l[k] + l_minus*al_r[k];
165 : 0 : flx[densityIdx(nmat, k)] = l_plus*u[0][densityIdx(nmat, k)]
166 : 0 : + l_minus*u[1][densityIdx(nmat, k)];
167 : 0 : flx[energyIdx(nmat, k)] = l_plus*hml[k] + l_minus*hmr[k];
168 : : }
169 : :
170 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
171 : : {
172 : 0 : flx[momentumIdx(nmat, idir)] = l_plus*u[0][momentumIdx(nmat, idir)]
173 : 0 : + l_minus*u[1][momentumIdx(nmat, idir)]
174 : 0 : + p12*fn[idir];
175 : : }
176 : :
177 : 0 : l_plus = l_plus/( vnl + std::copysign(1.0e-12, vnl) );
178 : 0 : l_minus = l_minus/( vnr + std::copysign(1.0e-12, vnr) );
179 : :
180 : : // Store Riemann-advected partial pressures
181 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
182 [ - - ]: 0 : flx.push_back( l_plus*pml[k] + l_minus*pmr[k] );
183 : :
184 : : // Store Riemann velocity
185 [ - - ]: 0 : flx.push_back( vriem );
186 : :
187 [ - - ][ - - ]: 0 : Assert( flx.size() == (3*nmat+3+nmat+1), "Size of multi-material flux "
[ - - ][ - - ]
188 : : "vector incorrect" );
189 : :
190 : 0 : return flx;
191 : : }
192 : :
193 : : //! Flux type accessor
194 : : //! \return Flux type
195 : : static ctr::FluxType type() noexcept { return ctr::FluxType::LDFSS; }
196 : : };
197 : :
198 : : } // inciter::
199 : :
200 : : #endif // LDFSS_h
|
