Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/LaxFriedrichsSolids.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 Lax-Friedrichs Riemann flux function for solids
9 : : \details This file implements the Lax-Friedrichs Riemann solver for solids.
10 : : */
11 : : // *****************************************************************************
12 : : #ifndef LaxFriedrichsSolids_h
13 : : #define LaxFriedrichsSolids_h
14 : :
15 : : #include <vector>
16 : :
17 : : #include "Fields.hpp"
18 : : #include "FunctionPrototypes.hpp"
19 : : #include "Inciter/Options/Flux.hpp"
20 : : #include "EoS/EOS.hpp"
21 : : #include "MultiMat/MultiMatIndexing.hpp"
22 : : #include "MultiMat/MiscMultiMatFns.hpp"
23 : :
24 : : namespace inciter {
25 : :
26 : : //! Lax-Friedrichs approximate Riemann solver for solids
27 : : struct LaxFriedrichsSolids {
28 : :
29 : : //! Lax-Friedrichs approximate Riemann solver flux function
30 : : //! \param[in] fn Face/Surface normal
31 : : //! \param[in] u Left and right unknown/state vector
32 : : //! \return Riemann flux solution according to Lax-Friedrichs, appended by
33 : : //! Riemann velocities and volume-fractions.
34 : : //! \note The function signature must follow tk::RiemannFluxFn
35 : : static tk::RiemannFluxFn::result_type
36 : 0 : flux( const std::vector< EOS >& mat_blk,
37 : : const std::array< tk::real, 3 >& fn,
38 : : const std::array< std::vector< tk::real >, 2 >& u,
39 : : const std::vector< std::array< tk::real, 3 > >& = {} )
40 : : {
41 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
42 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
43 : :
44 : 0 : auto nsld = numSolids(nmat, solidx);
45 : 0 : auto ncomp = u[0].size()-(3+nmat+nsld*6);
46 [ - - ][ - - ]: 0 : std::vector< tk::real > flx( ncomp, 0 ), fluxl(ncomp, 0), fluxr(ncomp,0);
[ - - ][ - - ]
47 : :
48 : : // Primitive variables
49 : : tk::real rhol(0.0), rhor(0.0);
50 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
51 : : {
52 : 0 : rhol += u[0][densityIdx(nmat, k)];
53 : 0 : rhor += u[1][densityIdx(nmat, k)];
54 : : }
55 : :
56 : : // Independently limited velocities for advection
57 [ - - ]: 0 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
58 : 0 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
59 : 0 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
60 : 0 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
61 : 0 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
62 : 0 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
63 : :
64 [ - - ][ - - ]: 0 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
[ - - ][ - - ]
65 [ - - ][ - - ]: 0 : pml(nmat, 0.0), pmr(nmat, 0.0),
[ - - ][ - - ]
66 [ - - ][ - - ]: 0 : am_l(nmat, 0.0),
67 [ - - ][ - - ]: 0 : am_r(nmat, 0.0);
68 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > > g_l, g_r,
69 : : asig_l, asig_r;
70 : : std::vector< std::array< tk::real, 3 > > asign_l, asign_r;
71 : : std::array< std::array< tk::real, 3 >, 3 > gn_l, gn_r;
72 : 0 : std::array< tk::real, 3 > sign_l {{0, 0, 0}}, sign_r {{0, 0, 0}};
73 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
74 : : {
75 : : // Left state
76 [ - - ]: 0 : al_l[k] = u[0][volfracIdx(nmat, k)];
77 [ - - ]: 0 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
78 : :
79 : : // inv deformation gradient and Cauchy stress tensors
80 [ - - ]: 0 : g_l.push_back(getDeformGrad(nmat, k, u[0]));
81 [ - - ]: 0 : asig_l.push_back(getCauchyStress(nmat, k, ncomp, u[0]));
82 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) asig_l[k][i][i] -= pml[k];
83 : :
84 : : // normal stress (traction) vector
85 : 0 : asign_l.push_back(tk::matvec(asig_l[k], fn));
86 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
87 : 0 : sign_l[i] += asign_l[k][i];
88 : :
89 : : // rotate deformation gradient tensor for speed of sound in normal dir
90 [ - - ]: 0 : gn_l = tk::rotateTensor(g_l[k], fn);
91 [ - - ][ - - ]: 0 : am_l[k] = mat_blk[k].compute< EOS::soundspeed >(
92 [ - - ]: 0 : u[0][densityIdx(nmat, k)], pml[k], al_l[k], k, gn_l );
93 : :
94 : : // Right state
95 [ - - ]: 0 : al_r[k] = u[1][volfracIdx(nmat, k)];
96 [ - - ]: 0 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
97 : :
98 : : // inv deformation gradient and Cauchy stress tensors
99 [ - - ]: 0 : g_r.push_back(getDeformGrad(nmat, k, u[1]));
100 [ - - ]: 0 : asig_r.push_back(getCauchyStress(nmat, k, ncomp, u[1]));
101 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) asig_r[k][i][i] -= pmr[k];
102 : :
103 : : // normal stress (traction) vector
104 : 0 : asign_r.push_back(tk::matvec(asig_r[k], fn));
105 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
106 : 0 : sign_r[i] += asign_r[k][i];
107 : :
108 : : // rotate deformation gradient tensor for speed of sound in normal dir
109 [ - - ]: 0 : gn_r = tk::rotateTensor(g_r[k], fn);
110 [ - - ]: 0 : am_r[k] = mat_blk[k].compute< EOS::soundspeed >(
111 [ - - ]: 0 : u[1][densityIdx(nmat, k)], pmr[k], al_r[k], k, gn_r );
112 : : }
113 : :
114 : : // Mixture speed of sound
115 : 0 : tk::real ac_l(0.0), ac_r(0.0);
116 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
117 : : {
118 : 0 : ac_l += (u[0][densityIdx(nmat,k)]*am_l[k]*am_l[k]);
119 : 0 : ac_r += (u[1][densityIdx(nmat,k)]*am_r[k]*am_r[k]);
120 : : }
121 : 0 : ac_l = std::sqrt( ac_l/rhol );
122 : 0 : ac_r = std::sqrt( ac_r/rhor );
123 : :
124 : : // Face-normal velocities from advective velocities
125 : 0 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
126 [ - - ]: 0 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
127 : :
128 : : // Maximum eignevalue and Riemann velocity
129 [ - - ][ - - ]: 0 : auto lambda = std::max(std::abs(vnl), std::abs(vnr)) + std::max(ac_l, ac_r);
130 : 0 : auto vriem = 0.5 * (vnl + vnr);
131 : :
132 : : // Conservative fluxes
133 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
134 : : {
135 : : // Left fluxes
136 : 0 : fluxl[volfracIdx(nmat, k)] = vnl*al_l[k];
137 : 0 : fluxl[densityIdx(nmat, k)] = vnl*u[0][densityIdx(nmat, k)];
138 : 0 : fluxl[energyIdx(nmat, k)] = vnl*u[0][energyIdx(nmat, k)];
139 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
140 : 0 : fluxl[energyIdx(nmat, k)] -= u[0][ncomp+velocityIdx(nmat,i)] *
141 : 0 : asign_l[k][i];
142 : : }
143 : :
144 : : // fluxes for inv deformation gradient tensor
145 [ - - ]: 0 : if (solidx[k] > 0) {
146 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
147 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
148 : 0 : fluxl[deformIdx(nmat,solidx[k],i,j)] = (
149 : 0 : g_l[k][i][0] * ul +
150 : 0 : g_l[k][i][1] * vl +
151 : 0 : g_l[k][i][2] * wl ) * fn[j];
152 : : }
153 : :
154 : : // Right fluxes
155 : 0 : fluxr[volfracIdx(nmat, k)] = vnr*al_r[k];
156 : 0 : fluxr[densityIdx(nmat, k)] = vnr*u[1][densityIdx(nmat, k)];
157 : 0 : fluxr[energyIdx(nmat, k)] = vnr*u[1][energyIdx(nmat, k)];
158 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
159 : 0 : fluxr[energyIdx(nmat, k)] -= u[1][ncomp+velocityIdx(nmat,i)] *
160 : 0 : asign_r[k][i];
161 : : }
162 : :
163 : : // fluxes for inv deformation gradient tensor
164 [ - - ]: 0 : if (solidx[k] > 0) {
165 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
166 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
167 : 0 : fluxr[deformIdx(nmat,solidx[k],i,j)] = (
168 : 0 : g_r[k][i][0] * ur +
169 : 0 : g_r[k][i][1] * vr +
170 : 0 : g_r[k][i][2] * wr ) * fn[j];
171 : : }
172 : : }
173 : :
174 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
175 : : {
176 : 0 : fluxl[momentumIdx(nmat, idir)] = vnl*u[0][momentumIdx(nmat, idir)]
177 : 0 : - sign_l[idir];
178 : 0 : fluxr[momentumIdx(nmat, idir)] = vnr*u[1][momentumIdx(nmat, idir)]
179 : 0 : - sign_r[idir];
180 : : }
181 : :
182 : : // Numerical flux function
183 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
184 : 0 : flx[c] = 0.5 * (fluxl[c] + fluxr[c] - lambda*(u[1][c] - u[0][c]));
185 : :
186 : : // Store Riemann-advected partial pressures (nmat)
187 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
188 [ - - ]: 0 : flx.push_back( 0.5*(pml[k]+pmr[k]) );
189 : :
190 : : // Store Riemann velocity (1)
191 [ - - ]: 0 : flx.push_back( vriem );
192 : :
193 : : // Flux vector splitting
194 : : auto l_plus = 0.5 * (vriem + std::fabs(vriem));
195 : : auto l_minus = 0.5 * (vriem - std::fabs(vriem));
196 : :
197 : : l_plus = l_plus/( std::fabs(vriem) + 1.0e-12 );
198 : : l_minus = l_minus/( std::fabs(vriem) + 1.0e-12 );
199 : :
200 : : // Store Riemann asign_ij (3*nsld)
201 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
202 [ - - ]: 0 : if (solidx[k] > 0) {
203 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
204 [ - - ][ - - ]: 0 : flx.push_back( 0.5 * (asign_l[k][i] + asign_r[k][i]) );
205 : : }
206 : : }
207 : :
208 : : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
209 : : "multi-material flux vector incorrect" );
210 : :
211 : 0 : return flx;
212 : : }
213 : :
214 : : ////! Flux type accessor
215 : : ////! \return Flux type
216 : : //static ctr::FluxType type() noexcept {
217 : : // return ctr::FluxType::LaxFriedrichs; }
218 : : };
219 : :
220 : : } // inciter::
221 : :
222 : : #endif // LaxFriedrichsSolids_h
|
