Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/HLLCMultiMat.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2023 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief HLLC Riemann flux function for solids
9 : : \details This file implements the HLLC Riemann solver for solids.
10 : : Ref. Ndanou, S., Favrie, N., & Gavrilyuk, S. (2015). Multi-solid
11 : : and multi-fluid diffuse interface model: Applications to dynamic
12 : : fracture and fragmentation. Journal of Computational Physics, 295,
13 : : 523-555.
14 : : */
15 : : // *****************************************************************************
16 : : #ifndef HLLCMultiMat_h
17 : : #define HLLCMultiMat_h
18 : :
19 : : #include <vector>
20 : :
21 : : #include "Fields.hpp"
22 : : #include "FunctionPrototypes.hpp"
23 : : #include "Inciter/Options/Flux.hpp"
24 : : #include "EoS/EOS.hpp"
25 : : #include "MultiMat/MultiMatIndexing.hpp"
26 : : #include "MultiMat/MiscMultiMatFns.hpp"
27 : :
28 : : namespace inciter {
29 : :
30 : : //! HLLC approximate Riemann solver for solids
31 : : struct HLLCMultiMat {
32 : :
33 : : //! HLLC approximate Riemann solver flux function
34 : : //! \param[in] fn Face/Surface normal
35 : : //! \param[in] u Left and right unknown/state vector
36 : : //! \return Riemann solution according to Harten-Lax-van Leer-Contact
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 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
46 : :
47 : 0 : auto nsld = numSolids(nmat, solidx);
48 : 0 : auto ncomp = u[0].size()-(3+nmat+nsld*6);
49 : 0 : std::vector< tk::real > flx(ncomp, 0);
50 : :
51 : : // Primitive variables
52 : : // -------------------------------------------------------------------------
53 : : tk::real rhol(0.0), rhor(0.0);
54 [ - - ]: 0 : for (size_t k=0; k<nmat; ++k) {
55 : 0 : rhol += u[0][densityIdx(nmat, k)];
56 : 0 : rhor += u[1][densityIdx(nmat, k)];
57 : : }
58 : :
59 [ - - ]: 0 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
60 : 0 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
61 : 0 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
62 : 0 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
63 : 0 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
64 : 0 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
65 : :
66 : : // Outer states
67 : : // -------------------------------------------------------------------------
68 : : tk::real pl(0.0), pr(0.0);
69 [ - - ][ - - ]: 0 : std::vector< tk::real > apl(nmat, 0.0), apr(nmat, 0.0);
[ - - ][ - - ]
70 : : tk::real acl(0.0), acr(0.0);
71 : :
72 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
73 : : // Left state
74 [ - - ]: 0 : apl[k] = u[0][ncomp+pressureIdx(nmat, k)];
75 : 0 : pl += apl[k];
76 [ - - ]: 0 : auto amatl = mat_blk[k].compute< EOS::soundspeed >(
77 : : u[0][densityIdx(nmat, k)], apl[k], u[0][volfracIdx(nmat, k)], k );
78 : :
79 : : // Right state
80 [ - - ]: 0 : apr[k] = u[1][ncomp+pressureIdx(nmat, k)];
81 : 0 : pr += apr[k];
82 [ - - ]: 0 : auto amatr = mat_blk[k].compute< EOS::soundspeed >(
83 : : u[1][densityIdx(nmat, k)], apr[k], u[1][volfracIdx(nmat, k)], k );
84 : :
85 : : // Mixture speed of sound
86 : 0 : acl += u[0][densityIdx(nmat, k)] * amatl * amatl;
87 : 0 : acr += u[1][densityIdx(nmat, k)] * amatr * amatr;
88 : : }
89 : 0 : acl = std::sqrt(acl/rhol);
90 : 0 : acr = std::sqrt(acr/rhor);
91 : :
92 : : // Face-normal velocities
93 : 0 : tk::real vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
94 : 0 : tk::real vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
95 : :
96 : : // Signal velocities
97 [ - - ]: 0 : auto Sl = std::min((vnl-acl), (vnr-acr));
98 [ - - ]: 0 : auto Sr = std::max((vnl+acl), (vnr+acr));
99 : 0 : auto Sm = ( rhor*vnr*(Sr-vnr) - rhol*vnl*(Sl-vnl) + pl-pr )
100 : 0 : /( rhor*(Sr-vnr) - rhol*(Sl-vnl) );
101 : :
102 : : // Middle-zone (star) variables
103 : : // -------------------------------------------------------------------------
104 : : tk::real pStar(0.0);
105 [ - - ][ - - ]: 0 : std::vector< tk::real > apStar(nmat, 0.0);
106 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
107 : 0 : apStar[k] = apl[k] + u[0][densityIdx(nmat, k)]*(vnl-Sl)*(vnl-Sm);
108 : 0 : pStar += apStar[k];
109 : : }
110 [ - - ]: 0 : auto uStar = u;
111 : :
112 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir) {
113 : 0 : uStar[0][momentumIdx(nmat, idir)] =
114 : 0 : ((Sl-vnl)*u[0][momentumIdx(nmat, idir)] + (pStar-pl)*fn[idir])/(Sl-Sm);
115 : 0 : uStar[1][momentumIdx(nmat, idir)] =
116 : 0 : ((Sr-vnr)*u[1][momentumIdx(nmat, idir)] + (pStar-pr)*fn[idir])/(Sr-Sm);
117 : : }
118 : :
119 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
120 : : // Left
121 : 0 : uStar[0][volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)];
122 : 0 : uStar[0][densityIdx(nmat, k)] =
123 : 0 : (Sl-vnl) * u[0][densityIdx(nmat, k)] / (Sl-Sm);
124 : 0 : uStar[0][energyIdx(nmat, k)] =
125 : 0 : ((Sl-vnl) * u[0][energyIdx(nmat, k)] - apl[k]*vnl + apStar[k]*Sm) / (Sl-Sm);
126 : :
127 : : // Right
128 : 0 : uStar[1][volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)];
129 : 0 : uStar[1][densityIdx(nmat, k)] =
130 : 0 : (Sr-vnr) * u[1][densityIdx(nmat, k)] / (Sr-Sm);
131 : 0 : uStar[1][energyIdx(nmat, k)] =
132 : 0 : ((Sr-vnr) * u[1][energyIdx(nmat, k)] - apr[k]*vnr + apStar[k]*Sm) / (Sr-Sm);
133 : : }
134 : :
135 : : // Numerical fluxes
136 : : // -------------------------------------------------------------------------
137 [ - - ]: 0 : if (Sl > 0.0) {
138 : :
139 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
140 : 0 : flx[momentumIdx(nmat, idir)] =
141 : 0 : u[0][momentumIdx(nmat, idir)] * vnl + pl*fn[idir];
142 : :
143 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
144 : 0 : flx[volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)] * vnl;
145 : 0 : flx[densityIdx(nmat, k)] = u[0][densityIdx(nmat, k)] * vnl;
146 : 0 : flx[energyIdx(nmat, k)] = (u[0][energyIdx(nmat, k)] + apl[k]) * vnl;
147 : : }
148 : :
149 : : // Quantities for non-conservative terms
150 : : // Store Riemann-advected partial pressures
151 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
152 [ - - ]: 0 : flx.push_back(apl[k]);
153 : : // Store Riemann velocity
154 [ - - ]: 0 : flx.push_back(vnl);
155 : : }
156 : :
157 [ - - ][ - - ]: 0 : else if (Sl <= 0.0 && Sm > 0.0) {
158 : :
159 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
160 : 0 : flx[momentumIdx(nmat, idir)] =
161 : 0 : uStar[0][momentumIdx(nmat, idir)] * Sm + pStar*fn[idir];
162 : :
163 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
164 : 0 : flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
165 : 0 : flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
166 : 0 : flx[energyIdx(nmat, k)] = (uStar[0][energyIdx(nmat, k)] + apStar[k]) * Sm;
167 : : }
168 : :
169 : : // Quantities for non-conservative terms
170 : : // Store Riemann-advected partial pressures
171 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
172 [ - - ]: 0 : flx.push_back(apStar[k]);
173 : : // Store Riemann velocity
174 [ - - ]: 0 : flx.push_back(Sm);
175 : : }
176 : :
177 [ - - ][ - - ]: 0 : else if (Sm <= 0.0 && Sr >= 0.0) {
178 : :
179 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
180 : 0 : flx[momentumIdx(nmat, idir)] =
181 : 0 : uStar[1][momentumIdx(nmat, idir)] * Sm + pStar*fn[idir];
182 : :
183 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
184 : 0 : flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
185 : 0 : flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
186 : 0 : flx[energyIdx(nmat, k)] = (uStar[1][energyIdx(nmat, k)] + apStar[k]) * Sm;
187 : : }
188 : :
189 : : // Quantities for non-conservative terms
190 : : // Store Riemann-advected partial pressures
191 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
192 [ - - ]: 0 : flx.push_back(apStar[k]);
193 : : // Store Riemann velocity
194 [ - - ]: 0 : flx.push_back(Sm);
195 : : }
196 : :
197 : : else {
198 : :
199 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
200 : 0 : flx[momentumIdx(nmat, idir)] =
201 : 0 : u[1][momentumIdx(nmat, idir)] * vnr + pr*fn[idir];
202 : :
203 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
204 : 0 : flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr;
205 : 0 : flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr;
206 : 0 : flx[energyIdx(nmat, k)] = (u[1][energyIdx(nmat, k)] + apr[k]) * vnr;
207 : : }
208 : :
209 : : // Quantities for non-conservative terms
210 : : // Store Riemann-advected partial pressures
211 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
212 [ - - ]: 0 : flx.push_back(apr[k]);
213 : : // Store Riemann velocity
214 [ - - ]: 0 : flx.push_back(vnr);
215 : : }
216 : :
217 : : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
218 : : "multi-material flux vector incorrect" );
219 : :
220 : 0 : return flx;
221 : : }
222 : :
223 : : ////! Flux type accessor
224 : : ////! \return Flux type
225 : : //static ctr::FluxType type() noexcept {
226 : : // return ctr::FluxType::HLLCMultiMat; }
227 : : };
228 : :
229 : : } // inciter::
230 : :
231 : : #endif // HLLCMultiMat_h
|
