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 : : 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 : : 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 : : 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 : : 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 : : 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( 1.0, ml );
112 : 0 : auto msr = splitmach_ausm( 1.0, mr );
113 : :
114 : : // Modified Mach polynomials
115 : : //// mtilde for 2nd order Mach polys (vanLeer)
116 : : //auto beta_l = -std::max(0.0, 1.0 - std::floor(std::abs(ml)));
117 : : //auto beta_r = -std::max(0.0, 1.0 - std::floor(std::abs(mr)));
118 : : //auto mtilde = 0.25*beta_l*beta_r*std::pow((std::sqrt(0.5*(ml*ml+mr*mr)) - 1.0), 2.0);
119 : : // mtilde for general-order Mach polys (AUSM+)
120 [ - - ][ - - ]: 0 : auto mtilde = 0.5*(msl[0] - std::max(0.0,ml) - msr[1] + std::min(0.0,mr));
121 [ - - ]: 0 : auto dp = pl - pr;
122 : : auto delta = 1.0;
123 : : auto mtilde_plus = mtilde*
124 [ - - ]: 0 : std::max( 0.0, 1.0 - (dp+delta*std::abs(dp))/(2.0*rhol*ac12*ac12) );
125 : : auto mtilde_minus = mtilde*
126 [ - - ]: 0 : std::max( 0.0, 1.0 + (dp-delta*std::abs(dp))/(2.0*rhor*ac12*ac12) );
127 : 0 : auto C_plus = msl[0] - mtilde_plus;
128 : 0 : auto C_minus = msr[1] + mtilde_minus;
129 : :
130 : : // Flux vector splitting
131 : 0 : auto l_plus = C_plus * ac12;
132 : 0 : auto l_minus = C_minus * ac12;
133 : 0 : auto vriem = l_plus + l_minus;
134 : :
135 : : // Riemann pressure
136 : 0 : auto p12 = 0.5*(pl+pr) + 0.5*(pl-pr)*(msl[2]-msr[3]) +
137 : 0 : rho12*ac12*ac12*(msl[2]+msr[3]-1.0);
138 : :
139 : : // Conservative fluxes
140 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
141 : : {
142 : 0 : flx[volfracIdx(nmat, k)] = l_plus*al_l[k] + l_minus*al_r[k];
143 : 0 : flx[densityIdx(nmat, k)] = l_plus*u[0][densityIdx(nmat, k)]
144 : 0 : + l_minus*u[1][densityIdx(nmat, k)];
145 : 0 : flx[energyIdx(nmat, k)] = l_plus*hml[k] + l_minus*hmr[k];
146 : : }
147 : :
148 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
149 : : {
150 : 0 : flx[momentumIdx(nmat, idir)] = l_plus*u[0][momentumIdx(nmat, idir)]
151 : 0 : + l_minus*u[1][momentumIdx(nmat, idir)]
152 : 0 : + p12*fn[idir];
153 : : }
154 : :
155 : 0 : l_plus = l_plus/( vnl + std::copysign(1.0e-12, vnl) );
156 : 0 : l_minus = l_minus/( vnr + std::copysign(1.0e-12, vnr) );
157 : :
158 : : // Store Riemann-advected partial pressures
159 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
160 [ - - ][ - - ]: 0 : flx.push_back( l_plus*pml[k] + l_minus*pmr[k] );
161 : :
162 : : // Store Riemann velocity
163 [ - - ]: 0 : flx.push_back( vriem );
164 : :
165 : : Assert( flx.size() == (3*nmat+3+nmat+1), "Size of multi-material flux "
166 : : "vector incorrect" );
167 : :
168 : 0 : return flx;
169 : : }
170 : :
171 : : //! Flux type accessor
172 : : //! \return Flux type
173 : : static ctr::FluxType type() noexcept { return ctr::FluxType::LDFSS; }
174 : : };
175 : :
176 : : } // inciter::
177 : :
178 : : #endif // LDFSS_h
|
