Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Reconstruction.cpp
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 Reconstruction for reconstructed discontinuous Galerkin methods
9 : : \details This file contains functions that reconstruct an "n"th order
10 : : polynomial to an "n+1"th order polynomial using a least-squares
11 : : reconstruction procedure.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <array>
16 : : #include <vector>
17 : : #include <iostream>
18 : : #include <iomanip>
19 : :
20 : : #include "Vector.hpp"
21 : : #include "Around.hpp"
22 : : #include "Base/HashMapReducer.hpp"
23 : : #include "Reconstruction.hpp"
24 : : #include "Inciter/Options/PDE.hpp"
25 : : #include "MultiMat/MultiMatIndexing.hpp"
26 : : #include "Inciter/InputDeck/InputDeck.hpp"
27 : : #include "Limiter.hpp"
28 : : #include "Integrate/Quadrature.hpp"
29 : :
30 : : namespace inciter {
31 : : extern ctr::InputDeck g_inputdeck;
32 : : }
33 : :
34 : : namespace tk {
35 : :
36 : : void
37 : 2894560 : recoLeastSqExtStencil(
38 : : std::size_t rdof,
39 : : std::size_t e,
40 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
41 : : const std::vector< std::size_t >& inpoel,
42 : : const Fields& geoElem,
43 : : Fields& W,
44 : : const std::vector< std::size_t >& varList )
45 : : // *****************************************************************************
46 : : // \brief Reconstruct the second-order solution using least-squares approach
47 : : // from an extended stencil involving the node-neighbors
48 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
49 : : //! \param[in] e Element whoes solution is being reconstructed
50 : : //! \param[in] esup Elements surrounding points
51 : : //! \param[in] inpoel Element-node connectivity
52 : : //! \param[in] geoElem Element geometry array
53 : : //! \param[in,out] W Solution vector to be reconstructed at recent time step
54 : : //! \param[in] varList List of indices in W, that need to be reconstructed
55 : : //! \details A second-order (piecewise linear) solution polynomial is obtained
56 : : //! from the first-order (piecewise constant) FV solutions by using a
57 : : //! least-squares (LS) reconstruction process. This LS reconstruction function
58 : : //! using the nodal-neighbors of a cell, to get an overdetermined system of
59 : : //! equations for the derivatives of the solution. This overdetermined system
60 : : //! is solved in the least-squares sense using the normal equations approach.
61 : : // *****************************************************************************
62 : : {
63 : : // Element centroid
64 [ + - ]: 2894560 : auto ex = geoElem(e,1);
65 [ + - ]: 2894560 : auto ey = geoElem(e,2);
66 [ + - ]: 2894560 : auto ez = geoElem(e,3);
67 : :
68 : 2894560 : tk::real m_xx = 0.0, m_xy = 0.0, m_xz = 0.0,
69 : 2894560 : m_yy = 0.0, m_yz = 0.0, m_zz = 0.0;
70 : :
71 : : // RHS for each variable in varList
72 : 2894560 : auto nvar = varList.size();
73 [ + - ]: 5789120 : std::vector< std::array< tk::real,3 > > rhs(nvar);
74 [ + + ]: 15281860 : for (auto& r : rhs) r = {{0.0, 0.0, 0.0}};
75 : :
76 : : // Cache P0 of the target cell for all variables (used for all neighbors)
77 [ + - ]: 5789120 : std::vector< tk::real > w0(nvar);
78 [ + + ]: 15281860 : for (std::size_t iv=0; iv<nvar; ++iv) {
79 : 12387300 : auto c = varList[iv];
80 : 12387300 : auto mark = c*rdof;
81 [ + - ]: 12387300 : w0[iv] = W(e, mark + 0);
82 : : }
83 : :
84 : : // Build LHS and accumulate RHS (neighbor outer loop)
85 [ + + ]: 14472800 : for (std::size_t lp=0; lp<4; ++lp) {
86 : 11578240 : auto p = inpoel[4*e + lp];
87 [ + - ]: 11578240 : auto& pesup = cref_find(esup, p);
88 : :
89 [ + + ]: 192571840 : for (auto er : pesup) {
90 [ + - ]: 180993600 : auto dx = geoElem(er,1) - ex;
91 [ + - ]: 180993600 : auto dy = geoElem(er,2) - ey;
92 [ + - ]: 180993600 : auto dz = geoElem(er,3) - ez;
93 : :
94 : : // LHS (use symmetry)
95 : 180993600 : m_xx += dx*dx;
96 : 180993600 : m_xy += dx*dy;
97 : 180993600 : m_xz += dx*dz;
98 : 180993600 : m_yy += dy*dy;
99 : 180993600 : m_yz += dy*dz;
100 : 180993600 : m_zz += dz*dz;
101 : :
102 : : // RHS for all variables
103 [ + + ]: 957705300 : for (std::size_t iv=0; iv<nvar; ++iv) {
104 : 776711700 : auto c = varList[iv];
105 : 776711700 : auto mark = c*rdof;
106 [ + - ]: 776711700 : auto dW = W(er, mark + 0) - w0[iv];
107 : 776711700 : rhs[iv][0] += dx * dW;
108 : 776711700 : rhs[iv][1] += dy * dW;
109 : 776711700 : rhs[iv][2] += dz * dW;
110 : : }
111 : : }
112 : : }
113 : :
114 : : // Form full symmetric 3×3 and factor once (reused for all RHS)
115 : : std::array<std::array<tk::real,3>,3> A {{
116 : : {{ m_xx, m_xy, m_xz }},
117 : : {{ m_xy, m_yy, m_yz }},
118 : : {{ m_xz, m_yz, m_zz }}
119 : 2894560 : }};
120 : : std::array<std::array<tk::real,3>,3> L;
121 [ + - ]: 2894560 : chol3x3(A, L);
122 : :
123 : : // Solve for and update the P1 dofs with the reconstructioned gradients.
124 : : // Since this reconstruction does not affect the cell-averaged solution,
125 : : // W(e,mark+0) is unchanged.
126 [ + + ]: 15281860 : for (std::size_t iv=0; iv<nvar; ++iv) {
127 : 12387300 : auto c = varList[iv];
128 : 12387300 : auto mark = c*rdof;
129 : 12387300 : auto ux = solve_chol3x3(L, rhs[iv]);
130 : :
131 [ + - ]: 12387300 : W(e, mark + 1) = ux[0];
132 [ + - ]: 12387300 : W(e, mark + 2) = ux[1];
133 [ + - ]: 12387300 : W(e, mark + 3) = ux[2];
134 : : }
135 : 2894560 : }
136 : :
137 : : void
138 : 2894560 : transform_P0P1( std::size_t rdof,
139 : : std::size_t e,
140 : : const std::vector< std::size_t >& inpoel,
141 : : const UnsMesh::Coords& coord,
142 : : Fields& W,
143 : : const std::vector< std::size_t >& varList )
144 : : // *****************************************************************************
145 : : // Transform the reconstructed P1-derivatives to the Dubiner dofs
146 : : //! \param[in] rdof Total number of reconstructed dofs
147 : : //! \param[in] e Element for which reconstruction is being calculated
148 : : //! \param[in] inpoel Element-node connectivity
149 : : //! \param[in] coord Array of nodal coordinates
150 : : //! \param[in,out] W Second-order reconstructed vector which gets transformed to
151 : : //! the Dubiner reference space
152 : : //! \param[in] varList List of indices in W, that need to be reconstructed
153 : : //! \details Since the DG solution (and the primitive quantities) are assumed to
154 : : //! be stored in the Dubiner space, this transformation from Taylor
155 : : //! coefficients to Dubiner coefficients is necessary.
156 : : // *****************************************************************************
157 : : {
158 : 2894560 : const auto& cx = coord[0];
159 : 2894560 : const auto& cy = coord[1];
160 : 2894560 : const auto& cz = coord[2];
161 : :
162 : : // Extract the element coordinates
163 : : std::array< std::array< real, 3>, 4 > coordel {{
164 : 8683680 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
165 : 8683680 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
166 : 8683680 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
167 : 8683680 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
168 : 31840160 : }};
169 : :
170 : : auto jacInv =
171 : 2894560 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
172 : :
173 : : // Compute the derivatives of basis function for DG(P1)
174 [ + - ]: 5789120 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
175 : :
176 [ + + ]: 15281860 : for (std::size_t i=0; i<varList.size(); ++i)
177 : : {
178 : 12387300 : auto mark = varList[i]*rdof;
179 : :
180 : : // solve system using Cramer's rule
181 : 37161900 : auto ux = tk::cramer( {{ {{dBdx[0][1], dBdx[0][2], dBdx[0][3]}},
182 : 37161900 : {{dBdx[1][1], dBdx[1][2], dBdx[1][3]}},
183 : 37161900 : {{dBdx[2][1], dBdx[2][2], dBdx[2][3]}} }},
184 [ + - ]: 12387300 : {{ W(e,mark+1),
185 : 24774600 : W(e,mark+2),
186 [ + - ][ + - ]: 123873000 : W(e,mark+3) }} );
187 : :
188 : : // replace physical derivatives with transformed dofs
189 [ + - ]: 12387300 : W(e,mark+1) = ux[0];
190 [ + - ]: 12387300 : W(e,mark+2) = ux[1];
191 [ + - ]: 12387300 : W(e,mark+3) = ux[2];
192 : : }
193 : 2894560 : }
194 : :
195 : : void
196 : 6757100 : THINCReco( std::size_t rdof,
197 : : std::size_t nmat,
198 : : std::size_t e,
199 : : const std::vector< std::size_t >& inpoel,
200 : : const UnsMesh::Coords& coord,
201 : : const Fields& geoElem,
202 : : const std::array< real, 3 >& ref_xp,
203 : : const Fields& U,
204 : : const Fields& P,
205 : : bool intInd,
206 : : const std::vector< std::size_t >& matInt,
207 : : [[maybe_unused]] const std::vector< real >& vfmin,
208 : : [[maybe_unused]] const std::vector< real >& vfmax,
209 : : std::vector< real >& state )
210 : : // *****************************************************************************
211 : : // Compute THINC reconstructions at quadrature point for multi-material flows
212 : : //! \param[in] rdof Total number of reconstructed dofs
213 : : //! \param[in] nmat Total number of materials
214 : : //! \param[in] e Element for which interface reconstruction is being calculated
215 : : //! \param[in] inpoel Element-node connectivity
216 : : //! \param[in] coord Array of nodal coordinates
217 : : //! \param[in] geoElem Element geometry array
218 : : //! \param[in] ref_xp Quadrature point in reference space
219 : : //! \param[in] U Solution vector
220 : : //! \param[in] P Vector of primitives
221 : : //! \param[in] intInd Boolean which indicates if the element contains a
222 : : //! material interface
223 : : //! \param[in] matInt Array indicating which material has an interface
224 : : //! \param[in] vfmin Vector containing min volume fractions for each material
225 : : //! in this cell
226 : : //! \param[in] vfmax Vector containing max volume fractions for each material
227 : : //! in this cell
228 : : //! \param[in,out] state Unknown/state vector at quadrature point, modified
229 : : //! if near interfaces using THINC
230 : : //! \details This function is an interface for the multimat PDEs that use the
231 : : //! algebraic multi-material THINC reconstruction. This particular function
232 : : //! should only be called for multimat.
233 : : // *****************************************************************************
234 : : {
235 : : using inciter::volfracDofIdx;
236 : : using inciter::densityDofIdx;
237 : : using inciter::momentumDofIdx;
238 : : using inciter::energyDofIdx;
239 : : using inciter::pressureDofIdx;
240 : : using inciter::velocityDofIdx;
241 : : using inciter::deformDofIdx;
242 : : using inciter::stressDofIdx;
243 : : using inciter::volfracIdx;
244 : : using inciter::densityIdx;
245 : : using inciter::momentumIdx;
246 : : using inciter::energyIdx;
247 : : using inciter::pressureIdx;
248 : : using inciter::velocityIdx;
249 : : using inciter::deformIdx;
250 : : using inciter::stressIdx;
251 : :
252 : : auto bparam = inciter::g_inputdeck.get< tag::multimat,
253 : 6757100 : tag::intsharp_param >();
254 : 6757100 : const auto ncomp = U.nprop()/rdof;
255 : : const auto& solidx = inciter::g_inputdeck.get< tag::matidxmap,
256 : 6757100 : tag::solidx >();
257 : :
258 : : // Step-1: Perform THINC reconstruction
259 : : // create a vector of volume-fractions and pass it to the THINC function
260 [ + - ]: 13514200 : std::vector< real > alSol(rdof*nmat, 0.0);
261 [ + - ]: 13514200 : std::vector< real > alReco(nmat, 0.0);
262 [ + + ]: 20271300 : for (std::size_t k=0; k<nmat; ++k) {
263 : 13514200 : auto mark = k*rdof;
264 [ + + ]: 67571000 : for (std::size_t i=0; i<rdof; ++i) {
265 [ + - ]: 54056800 : alSol[mark+i] = U(e, volfracDofIdx(nmat,k,rdof,i));
266 : : }
267 : : // initialize with TVD reconstructions which will be modified if near
268 : : // material interface
269 : 13514200 : alReco[k] = state[volfracIdx(nmat,k)];
270 : : }
271 [ + - ][ + - ]: 6757100 : THINCFunction(rdof, nmat, e, inpoel, coord, ref_xp, geoElem(e,0), bparam,
272 : : alSol, intInd, matInt, alReco);
273 : :
274 : : // check reconstructed volfracs for positivity
275 : 6757100 : bool neg_vf = false;
276 [ + + ]: 20271300 : for (std::size_t k=0; k<nmat; ++k) {
277 [ - + ][ - - ]: 13514200 : if (alReco[k] < 1e-16 && matInt[k] > 0) neg_vf = true;
[ - + ]
278 : : }
279 [ + + ]: 20271300 : for (std::size_t k=0; k<nmat; ++k) {
280 [ - + ]: 13514200 : if (neg_vf) {
281 [ - - ][ - - ]: 0 : std::cout << "Material-id: " << k << std::endl;
[ - - ]
282 [ - - ][ - - ]: 0 : std::cout << "Volume-fraction: " << std::setprecision(18) << alReco[k]
[ - - ]
283 [ - - ]: 0 : << std::endl;
284 [ - - ][ - - ]: 0 : std::cout << "Cell-avg vol-frac: " << std::setprecision(18) <<
285 [ - - ][ - - ]: 0 : U(e,volfracDofIdx(nmat,k,rdof,0)) << std::endl;
[ - - ]
286 [ - - ][ - - ]: 0 : std::cout << "Material-interface? " << intInd << std::endl;
[ - - ]
287 [ - - ][ - - ]: 0 : std::cout << "Mat-k-involved? " << matInt[k] << std::endl;
[ - - ]
288 : : }
289 : : }
290 [ - + ][ - - ]: 6757100 : if (neg_vf) Throw("Material has negative volume fraction after THINC "
[ - - ][ - - ]
291 : : "reconstruction.");
292 : :
293 : : // Step-2: Perform consistent reconstruction on other conserved quantities
294 [ + + ]: 6757100 : if (intInd)
295 : : {
296 : 239692 : auto rhobCC(0.0), rhobHO(0.0);
297 [ + + ]: 719076 : for (std::size_t k=0; k<nmat; ++k)
298 : : {
299 [ + - ]: 479384 : auto alCC = U(e, volfracDofIdx(nmat,k,rdof,0));
300 : 479384 : alCC = std::max(1e-14, alCC);
301 : :
302 [ + - ]: 479384 : if (matInt[k])
303 : : {
304 : 479384 : state[volfracIdx(nmat,k)] = alReco[k];
305 : 958768 : state[densityIdx(nmat,k)] = alReco[k]
306 [ + - ]: 479384 : * U(e, densityDofIdx(nmat,k,rdof,0))/alCC;
307 : 958768 : state[energyIdx(nmat,k)] = alReco[k]
308 [ + - ]: 479384 : * U(e, energyDofIdx(nmat,k,rdof,0))/alCC;
309 : 958768 : state[ncomp+pressureIdx(nmat,k)] = alReco[k]
310 [ + - ]: 479384 : * P(e, pressureDofIdx(nmat,k,rdof,0))/alCC;
311 [ - + ]: 479384 : if (solidx[k] > 0) {
312 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
313 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
314 : 0 : state[deformIdx(nmat,solidx[k],i,j)] =
315 [ - - ]: 0 : U(e, deformDofIdx(nmat,solidx[k],i,j,rdof,0));
316 : :
317 [ - - ]: 0 : for (std::size_t i=0; i<6; ++i)
318 : 0 : state[ncomp+stressIdx(nmat,solidx[k],i)] = alReco[k]
319 [ - - ]: 0 : * P(e, stressDofIdx(nmat,solidx[k],i,rdof,0))/alCC;
320 : : }
321 : : }
322 : :
323 [ + - ]: 479384 : rhobCC += U(e, densityDofIdx(nmat,k,rdof,0));
324 : 479384 : rhobHO += state[densityIdx(nmat,k)];
325 : : }
326 : :
327 : : // consistent reconstruction for bulk momentum
328 [ + + ]: 958768 : for (std::size_t i=0; i<3; ++i)
329 : : {
330 : 719076 : state[momentumIdx(nmat,i)] = rhobHO
331 [ + - ]: 719076 : * U(e, momentumDofIdx(nmat,i,rdof,0))/rhobCC;
332 : 719076 : state[ncomp+velocityIdx(nmat,i)] =
333 [ + - ]: 719076 : P(e, velocityDofIdx(nmat,i,rdof,0));
334 : : }
335 : : }
336 : 6757100 : }
337 : :
338 : : void
339 : 0 : THINCRecoTransport( std::size_t rdof,
340 : : std::size_t,
341 : : std::size_t e,
342 : : const std::vector< std::size_t >& inpoel,
343 : : const UnsMesh::Coords& coord,
344 : : const Fields& geoElem,
345 : : const std::array< real, 3 >& ref_xp,
346 : : const Fields& U,
347 : : const Fields&,
348 : : [[maybe_unused]] const std::vector< real >& vfmin,
349 : : [[maybe_unused]] const std::vector< real >& vfmax,
350 : : std::vector< real >& state )
351 : : // *****************************************************************************
352 : : // Compute THINC reconstructions at quadrature point for transport
353 : : //! \param[in] rdof Total number of reconstructed dofs
354 : : //! \param[in] e Element for which interface reconstruction is being calculated
355 : : //! \param[in] inpoel Element-node connectivity
356 : : //! \param[in] coord Array of nodal coordinates
357 : : //! \param[in] geoElem Element geometry array
358 : : //! \param[in] ref_xp Quadrature point in reference space
359 : : //! \param[in] U Solution vector
360 : : //! \param[in] vfmin Vector containing min volume fractions for each material
361 : : //! in this cell
362 : : //! \param[in] vfmax Vector containing max volume fractions for each material
363 : : //! in this cell
364 : : //! \param[in,out] state Unknown/state vector at quadrature point, modified
365 : : //! if near interfaces using THINC
366 : : //! \details This function is an interface for the transport PDEs that use the
367 : : //! algebraic multi-material THINC reconstruction. This particular function
368 : : //! should only be called for transport.
369 : : // *****************************************************************************
370 : : {
371 : : auto bparam = inciter::g_inputdeck.get< tag::transport,
372 : 0 : tag::intsharp_param >();
373 : 0 : auto ncomp = U.nprop()/rdof;
374 : :
375 : : // interface detection
376 [ - - ]: 0 : std::vector< std::size_t > matInt(ncomp, 0);
377 [ - - ]: 0 : std::vector< tk::real > alAvg(ncomp, 0.0);
378 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k)
379 [ - - ]: 0 : alAvg[k] = U(e, k*rdof);
380 [ - - ]: 0 : auto intInd = inciter::interfaceIndicator(ncomp, alAvg, matInt);
381 : :
382 : : // create a vector of volume-fractions and pass it to the THINC function
383 [ - - ]: 0 : std::vector< real > alSol(rdof*ncomp, 0.0);
384 : : // initialize with TVD reconstructions (modified if near interface)
385 [ - - ]: 0 : auto alReco = state;
386 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k) {
387 : 0 : auto mark = k*rdof;
388 [ - - ]: 0 : for (std::size_t i=0; i<rdof; ++i) {
389 [ - - ]: 0 : alSol[mark+i] = U(e,mark+i);
390 : : }
391 : : }
392 [ - - ][ - - ]: 0 : THINCFunction(rdof, ncomp, e, inpoel, coord, ref_xp, geoElem(e,0), bparam,
393 : : alSol, intInd, matInt, alReco);
394 : :
395 [ - - ]: 0 : state = alReco;
396 : 0 : }
397 : :
398 : : void
399 : 0 : THINCFunction_old( std::size_t rdof,
400 : : std::size_t nmat,
401 : : std::size_t e,
402 : : const std::vector< std::size_t >& inpoel,
403 : : const UnsMesh::Coords& coord,
404 : : const std::array< real, 3 >& ref_xp,
405 : : real vol,
406 : : real bparam,
407 : : const std::vector< real >& alSol,
408 : : bool intInd,
409 : : const std::vector< std::size_t >& matInt,
410 : : std::vector< real >& alReco )
411 : : // *****************************************************************************
412 : : // Old version of the Multi-Medium THINC reconstruction function
413 : : //! \param[in] rdof Total number of reconstructed dofs
414 : : //! \param[in] nmat Total number of materials
415 : : //! \param[in] e Element for which interface reconstruction is being calculated
416 : : //! \param[in] inpoel Element-node connectivity
417 : : //! \param[in] coord Array of nodal coordinates
418 : : //! \param[in] ref_xp Quadrature point in reference space
419 : : //! \param[in] vol Element volume
420 : : //! \param[in] bparam User specified Beta for THINC, from the input file
421 : : //! \param[in] alSol Volume fraction solution vector for element e
422 : : //! \param[in] intInd Interface indicator, true if e is interface element
423 : : //! \param[in] matInt Vector indicating materials which constitute interface
424 : : //! \param[in,out] alReco Unknown/state vector at quadrature point, which gets
425 : : //! modified if near interface using MM-THINC
426 : : //! \details This function computes the interface reconstruction using the
427 : : //! algebraic multi-material THINC reconstruction for each material at the
428 : : //! given (ref_xp) quadrature point. This function is based on the following:
429 : : //! Pandare A. K., Waltz J., & Bakosi J. (2021) Multi-Material Hydrodynamics
430 : : //! with Algebraic Sharp Interface Capturing. Computers & Fluids,
431 : : //! doi: https://doi.org/10.1016/j.compfluid.2020.104804.
432 : : //! This function will be removed after the newer version (see
433 : : //! THINCFunction_new) is sufficiently tested.
434 : : // *****************************************************************************
435 : : {
436 : : auto min_al = inciter::g_inputdeck.get< tag::multimat,
437 : 0 : tag::min_volumefrac >();
438 : :
439 : : // determine number of materials with interfaces in this cell
440 : 0 : auto epsl(1e-4), epsh(1e-1), bred(1.25), bmod(bparam);
441 : 0 : std::size_t nIntMat(0);
442 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
443 : : {
444 : 0 : auto alk = alSol[k*rdof];
445 [ - - ]: 0 : if (alk > epsl)
446 : : {
447 : 0 : ++nIntMat;
448 [ - - ][ - - ]: 0 : if ((alk > epsl) && (alk < epsh))
449 : 0 : bmod = std::min(bmod,
450 : 0 : (alk-epsl)/(epsh-epsl) * (bred - bparam) + bparam);
451 [ - - ]: 0 : else if (alk > epsh)
452 : 0 : bmod = bred;
453 : : }
454 : : }
455 : :
456 [ - - ]: 0 : if (nIntMat > 2) bparam = bmod;
457 : :
458 : : // compression parameter
459 : 0 : auto beta = bparam/std::cbrt(6.0*vol);
460 : :
461 [ - - ]: 0 : if (intInd)
462 : : {
463 : : // 1. Get unit normals to material interface
464 : :
465 : : // Compute Jacobian matrix for converting Dubiner dofs to derivatives
466 : 0 : const auto& cx = coord[0];
467 : 0 : const auto& cy = coord[1];
468 : 0 : const auto& cz = coord[2];
469 : :
470 : : std::array< std::array< real, 3>, 4 > coordel {{
471 : 0 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
472 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
473 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
474 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
475 : 0 : }};
476 : :
477 : : auto jacInv =
478 : 0 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
479 : :
480 [ - - ]: 0 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
481 : :
482 : : std::array< real, 3 > nInt;
483 [ - - ]: 0 : std::vector< std::array< real, 3 > > ref_n(nmat, {{0.0, 0.0, 0.0}});
484 : :
485 : : // Get normals
486 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
487 : : {
488 : : // Get derivatives from moments in Dubiner space
489 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
490 : 0 : nInt[i] = dBdx[i][1] * alSol[k*rdof+1]
491 : 0 : + dBdx[i][2] * alSol[k*rdof+2]
492 : 0 : + dBdx[i][3] * alSol[k*rdof+3];
493 : :
494 : 0 : auto nMag = std::sqrt(tk::dot(nInt, nInt)) + 1e-14;
495 : :
496 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
497 : 0 : nInt[i] /= nMag;
498 : :
499 : : // project interface normal onto local/reference coordinate system
500 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
501 : : {
502 : : std::array< real, 3 > axis{
503 : 0 : coordel[i+1][0]-coordel[0][0],
504 : 0 : coordel[i+1][1]-coordel[0][1],
505 : 0 : coordel[i+1][2]-coordel[0][2] };
506 : 0 : ref_n[k][i] = tk::dot(nInt, axis);
507 : : }
508 : : }
509 : :
510 : : // 2. Reconstruct volume fractions using THINC
511 : 0 : auto max_lim = 1.0 - (static_cast<tk::real>(nmat-1)*min_al);
512 : 0 : auto min_lim = min_al;
513 : 0 : auto sum_inter(0.0), sum_non_inter(0.0);
514 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
515 : : {
516 [ - - ]: 0 : if (matInt[k])
517 : : {
518 : : // get location of material interface (volume fraction 0.5) from the
519 : : // assumed tanh volume fraction distribution, and cell-averaged
520 : : // volume fraction
521 : 0 : auto alCC(alSol[k*rdof]);
522 : 0 : auto Ac(0.0), Bc(0.0), Qc(0.0);
523 : 0 : if ((std::abs(ref_n[k][0]) > std::abs(ref_n[k][1]))
524 [ - - ][ - - ]: 0 : && (std::abs(ref_n[k][0]) > std::abs(ref_n[k][2])))
[ - - ]
525 : : {
526 : 0 : Ac = std::exp(0.5*beta*ref_n[k][0]);
527 : 0 : Bc = std::exp(0.5*beta*(ref_n[k][1]+ref_n[k][2]));
528 : 0 : Qc = std::exp(0.5*beta*ref_n[k][0]*(2.0*alCC-1.0));
529 : : }
530 : 0 : else if ((std::abs(ref_n[k][1]) > std::abs(ref_n[k][0]))
531 [ - - ][ - - ]: 0 : && (std::abs(ref_n[k][1]) > std::abs(ref_n[k][2])))
[ - - ]
532 : : {
533 : 0 : Ac = std::exp(0.5*beta*ref_n[k][1]);
534 : 0 : Bc = std::exp(0.5*beta*(ref_n[k][0]+ref_n[k][2]));
535 : 0 : Qc = std::exp(0.5*beta*ref_n[k][1]*(2.0*alCC-1.0));
536 : : }
537 : : else
538 : : {
539 : 0 : Ac = std::exp(0.5*beta*ref_n[k][2]);
540 : 0 : Bc = std::exp(0.5*beta*(ref_n[k][0]+ref_n[k][1]));
541 : 0 : Qc = std::exp(0.5*beta*ref_n[k][2]*(2.0*alCC-1.0));
542 : : }
543 : 0 : auto d = std::log((1.0-Ac*Qc) / (Ac*Bc*(Qc-Ac))) / (2.0*beta);
544 : :
545 : : // THINC reconstruction
546 : 0 : auto al_c = 0.5 * (1.0 + std::tanh(beta*(tk::dot(ref_n[k], ref_xp) + d)));
547 : :
548 : 0 : alReco[k] = std::min(max_lim, std::max(min_lim, al_c));
549 : :
550 : 0 : sum_inter += alReco[k];
551 : : } else
552 : : {
553 : 0 : sum_non_inter += alReco[k];
554 : : }
555 : : // else, if this material does not have an interface close-by, the TVD
556 : : // reconstructions must be used for state variables. This is ensured by
557 : : // initializing the alReco vector as the TVD state.
558 : : }
559 : :
560 : : // Rescale volume fractions of interface-materials to ensure unit sum
561 : 0 : auto sum_rest = 1.0 - sum_non_inter;
562 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
563 [ - - ]: 0 : if(matInt[k])
564 : 0 : alReco[k] = alReco[k] * sum_rest / sum_inter;
565 : : }
566 : 0 : }
567 : :
568 : : void
569 : 6757100 : THINCFunction( std::size_t rdof,
570 : : std::size_t nmat,
571 : : std::size_t e,
572 : : const std::vector< std::size_t >& inpoel,
573 : : const UnsMesh::Coords& coord,
574 : : const std::array< real, 3 >& ref_xp,
575 : : real vol,
576 : : real bparam,
577 : : const std::vector< real >& alSol,
578 : : bool intInd,
579 : : const std::vector< std::size_t >& matInt,
580 : : std::vector< real >& alReco )
581 : : // *****************************************************************************
582 : : // Multi-Medium THINC reconstruction function for volume fractions
583 : : //! \param[in] rdof Total number of reconstructed dofs
584 : : //! \param[in] nmat Total number of materials
585 : : //! \param[in] e Element for which interface reconstruction is being calculated
586 : : //! \param[in] inpoel Element-node connectivity
587 : : //! \param[in] coord Array of nodal coordinates
588 : : //! \param[in] ref_xp Quadrature point in reference space
589 : : //! \param[in] vol Element volume
590 : : //! \param[in] bparam User specified Beta for THINC, from the input file
591 : : //! \param[in] alSol Volume fraction solution vector for element e
592 : : //! \param[in] intInd Interface indicator, true if e is interface element
593 : : //! \param[in] matInt Vector indicating materials which constitute interface
594 : : //! \param[in,out] alReco Unknown/state vector at quadrature point, which gets
595 : : //! modified if near interface using MM-THINC
596 : : //! \details This function computes the interface reconstruction using the
597 : : //! algebraic multi-material THINC reconstruction for each material at the
598 : : //! given (ref_xp) quadrature point. This function succeeds the older version
599 : : //! of the mm-THINC (see THINCFunction_old).
600 : : // *****************************************************************************
601 : : {
602 : : auto min_al = inciter::g_inputdeck.get< tag::multimat,
603 : 6757100 : tag::min_volumefrac >();
604 : :
605 : : // compression parameter
606 : 6757100 : auto beta = bparam/std::cbrt(6.0*vol);
607 : :
608 : : // If the cell is not material interface, return this function
609 [ + + ]: 6757100 : if (not intInd) return;
610 : :
611 : : // If the cell is material interface, THINC reconstruction is applied
612 : : // Step 1. Get unit normals to material interface
613 : : // -------------------------------------------------------------------------
614 : :
615 : : // Compute Jacobian matrix for converting Dubiner dofs to derivatives
616 : 239692 : const auto& cx = coord[0];
617 : 239692 : const auto& cy = coord[1];
618 : 239692 : const auto& cz = coord[2];
619 : :
620 : : std::array< std::array< real, 3>, 4 > coordel {{
621 : 719076 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
622 : 719076 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
623 : 719076 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
624 : 719076 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
625 : 2636612 : }};
626 : :
627 : : auto jacInv =
628 : 239692 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
629 : :
630 [ + - ]: 239692 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
631 : :
632 : : std::array< real, 3 > nInt;
633 : 239692 : std::array< real, 3 > ref_n{0.0, 0.0, 0.0};
634 : 239692 : auto almax(0.0);
635 : 239692 : std::size_t kmax(0);
636 : :
637 : : // Determine index of material present in majority
638 [ + + ]: 719076 : for (std::size_t k=0; k<nmat; ++k)
639 : : {
640 : 479384 : auto alk = alSol[k*rdof];
641 [ + - ]: 479384 : if (alk > almax)
642 : : {
643 : 479384 : almax = alk;
644 : 479384 : kmax = k;
645 : : }
646 : : }
647 : :
648 : : // Get normals of material present in majority
649 : : // Get derivatives from moments in Dubiner space
650 [ + + ]: 958768 : for (std::size_t i=0; i<3; ++i)
651 : 1438152 : nInt[i] = dBdx[i][1] * alSol[kmax*rdof+1]
652 : 719076 : + dBdx[i][2] * alSol[kmax*rdof+2]
653 : 719076 : + dBdx[i][3] * alSol[kmax*rdof+3];
654 : :
655 : 239692 : auto nMag = std::sqrt(tk::dot(nInt, nInt)) + 1e-14;
656 : :
657 [ + + ]: 958768 : for (std::size_t i=0; i<3; ++i)
658 : 719076 : nInt[i] /= nMag;
659 : :
660 : : // project interface normal onto local/reference coordinate system
661 [ + + ]: 958768 : for (std::size_t i=0; i<3; ++i)
662 : : {
663 : : std::array< real, 3 > axis{
664 : 719076 : coordel[i+1][0]-coordel[0][0],
665 : 719076 : coordel[i+1][1]-coordel[0][1],
666 : 1438152 : coordel[i+1][2]-coordel[0][2] };
667 : 719076 : ref_n[i] = tk::dot(nInt, axis);
668 : : }
669 : :
670 : : // Step 2. Reconstruct volume fraction of majority material using THINC
671 : : // -------------------------------------------------------------------------
672 : :
673 : 239692 : auto al_max = 1.0 - (static_cast<tk::real>(nmat-1)*min_al);
674 : 239692 : auto al_min = min_al;
675 : 239692 : auto alsum(0.0);
676 : : // get location of material interface (volume fraction 0.5) from the
677 : : // assumed tanh volume fraction distribution, and cell-averaged
678 : : // volume fraction
679 : 239692 : auto alCC(alSol[kmax*rdof]);
680 : 239692 : auto Ac(0.0), Bc(0.0), Qc(0.0);
681 : 239692 : if ((std::abs(ref_n[0]) > std::abs(ref_n[1]))
682 [ + + ][ + + ]: 239692 : && (std::abs(ref_n[0]) > std::abs(ref_n[2])))
[ + + ]
683 : : {
684 : 96141 : Ac = std::exp(0.5*beta*ref_n[0]);
685 : 96141 : Bc = std::exp(0.5*beta*(ref_n[1]+ref_n[2]));
686 : 96141 : Qc = std::exp(0.5*beta*ref_n[0]*(2.0*alCC-1.0));
687 : : }
688 : 143551 : else if ((std::abs(ref_n[1]) > std::abs(ref_n[0]))
689 [ + + ][ + + ]: 143551 : && (std::abs(ref_n[1]) > std::abs(ref_n[2])))
[ + + ]
690 : : {
691 : 106223 : Ac = std::exp(0.5*beta*ref_n[1]);
692 : 106223 : Bc = std::exp(0.5*beta*(ref_n[0]+ref_n[2]));
693 : 106223 : Qc = std::exp(0.5*beta*ref_n[1]*(2.0*alCC-1.0));
694 : : }
695 : : else
696 : : {
697 : 37328 : Ac = std::exp(0.5*beta*ref_n[2]);
698 : 37328 : Bc = std::exp(0.5*beta*(ref_n[0]+ref_n[1]));
699 : 37328 : Qc = std::exp(0.5*beta*ref_n[2]*(2.0*alCC-1.0));
700 : : }
701 : 239692 : auto d = std::log((1.0-Ac*Qc) / (Ac*Bc*(Qc-Ac))) / (2.0*beta);
702 : :
703 : : // THINC reconstruction
704 : 239692 : auto al_c = 0.5 * (1.0 + std::tanh(beta*(tk::dot(ref_n, ref_xp) + d)));
705 : :
706 : 239692 : alReco[kmax] = std::min(al_max, std::max(al_min, al_c));
707 : 239692 : alsum += alReco[kmax];
708 : :
709 : : // if this material does not have an interface close-by, the TVD
710 : : // reconstructions must be used for state variables. This is ensured by
711 : : // initializing the alReco vector as the TVD state.
712 [ + + ]: 719076 : for (std::size_t k=0; k<nmat; ++k) {
713 [ - + ]: 479384 : if (!matInt[k]) {
714 : 0 : alsum += alReco[k];
715 : : }
716 : : }
717 : :
718 : : // Step 3. Do multimaterial cell corrections
719 : : // -------------------------------------------------------------------------
720 : :
721 : : // distribute remaining volume to rest of materials
722 : 239692 : auto sum_left = 1.0 - alsum;
723 : 239692 : real den = 0.0;
724 [ + + ]: 719076 : for (std::size_t k=0; k<nmat; ++k) {
725 [ + - ][ + + ]: 479384 : if (matInt[k] && k != kmax) {
[ + + ]
726 : 239692 : auto mark = k * rdof;
727 : 239692 : alReco[k] = sum_left * alSol[mark];
728 : 239692 : den += alSol[mark];
729 : : }
730 : : }
731 : : // the distributed volfracs might be below al_min, correct that
732 : 239692 : real err = 0.0;
733 [ + + ]: 719076 : for (std::size_t k=0; k<nmat; ++k) {
734 [ + - ][ + + ]: 479384 : if (matInt[k] && k != kmax) {
[ + + ]
735 : 239692 : alReco[k] /= den;
736 [ - + ]: 239692 : if (alReco[k] < al_min) {
737 : 0 : err += al_min - alReco[k];
738 : 0 : alReco[k] = al_min;
739 : : }
740 : : }
741 : : }
742 : :
743 : : // balance out errors
744 : 239692 : alReco[kmax] -= err;
745 : : }
746 : :
747 : : void
748 : 0 : computeTemperaturesFV(
749 : : const std::vector< inciter::EOS >& mat_blk,
750 : : std::size_t nmat,
751 : : const std::vector< std::size_t >& inpoel,
752 : : const tk::UnsMesh::Coords& coord,
753 : : const tk::Fields& geoElem,
754 : : const tk::Fields& unk,
755 : : const tk::Fields& prim,
756 : : const std::vector< int >& srcFlag,
757 : : tk::Fields& T )
758 : : // *****************************************************************************
759 : : // Compute the temperatures based on FV conserved quantities
760 : : //! \param[in] mat_blk EOS material block
761 : : //! \param[in] nmat Number of materials in this PDE system
762 : : //! \param[in] inpoel Element-node connectivity
763 : : //! \param[in] coord Array of nodal coordinates
764 : : //! \param[in] geoElem Element geometry array
765 : : //! \param[in] unk Array of conservative variables
766 : : //! \param[in] prim Array of primitive variables
767 : : //! \param[in] srcFlag Whether the energy source was added
768 : : //! \param[in,out] T Array of material temperature dofs
769 : : //! \details This function computes the dofs of material temperatures based on
770 : : //! conservative quantities from an FV scheme, using EOS calls. It uses the
771 : : //! weak form m_{ij} T_i = \int T_{EOS}(rho, rhoE, u) B_j.
772 : : // *****************************************************************************
773 : : {
774 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
775 : 0 : std::size_t ncomp = unk.nprop()/rdof;
776 : 0 : std::size_t nprim = prim.nprop()/rdof;
777 : : const auto intsharp = inciter::g_inputdeck.get< tag::multimat,
778 : 0 : tag::intsharp >();
779 : 0 : auto nelem = unk.nunk();
780 : :
781 : 0 : auto L = tk::massMatrixDubiner();
782 : :
783 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e) {
784 [ - - ]: 0 : auto vole = geoElem(e,0);
785 : :
786 [ - - ]: 0 : std::vector< tk::real > R(nmat*rdof, 0.0);
787 : :
788 : 0 : std::size_t ng = 11;
789 : :
790 : : // Arrays for quadrature points
791 : 0 : std::array< std::vector< tk::real >, 3 > coordgp;
792 : 0 : std::vector< tk::real > wgp;
793 : :
794 [ - - ]: 0 : coordgp[0].resize( ng );
795 [ - - ]: 0 : coordgp[1].resize( ng );
796 [ - - ]: 0 : coordgp[2].resize( ng );
797 [ - - ]: 0 : wgp.resize( ng );
798 : :
799 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
800 : :
801 : : // Loop over quadrature points in element e
802 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp) {
803 : : // Compute the basis function
804 : 0 : auto B = tk::eval_basis( rdof, coordgp[0][igp], coordgp[1][igp],
805 [ - - ]: 0 : coordgp[2][igp] );
806 : :
807 : 0 : auto w = wgp[igp] * vole;
808 : :
809 : : // Evaluate the solution at quadrature point
810 : : auto state = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
811 : : nmat, e, inpoel, coord, geoElem,
812 : 0 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, unk, prim,
813 [ - - ]: 0 : srcFlag[e]);
814 : :
815 : : // Velocity vector at quadrature point
816 : : std::array< tk::real, 3 >
817 : 0 : vel{ state[ncomp+inciter::velocityIdx(nmat, 0)],
818 : 0 : state[ncomp+inciter::velocityIdx(nmat, 1)],
819 : 0 : state[ncomp+inciter::velocityIdx(nmat, 2)] };
820 : :
821 : : // Evaluate the right-hand-side vector (for temperature)
822 [ - - ]: 0 : for(std::size_t k=0; k<nmat; k++) {
823 [ - - ]: 0 : auto tk = mat_blk[k].compute< inciter::EOS::temperature >(
824 : 0 : state[inciter::densityIdx(nmat, k)], vel[0], vel[1], vel[2],
825 : 0 : state[inciter::energyIdx(nmat, k)],
826 : 0 : state[inciter::volfracIdx(nmat, k)] );
827 : 0 : auto mark = k * rdof;
828 [ - - ]: 0 : for(std::size_t idof=0; idof<rdof; idof++)
829 : 0 : R[mark+idof] += w * tk * B[idof];
830 : : }
831 : : }
832 : :
833 : : // Update the high order dofs of the temperature
834 [ - - ]: 0 : for(std::size_t k=0; k<nmat; k++) {
835 : 0 : auto mark = k * rdof;
836 [ - - ]: 0 : for(std::size_t idof=1; idof<rdof; idof++)
837 [ - - ]: 0 : T(e, mark+idof) = R[mark+idof] / (vole*L[idof]);
838 : : }
839 : : }
840 : 0 : }
841 : :
842 : : std::vector< tk::real >
843 : 92378342 : evalPolynomialSol( const std::vector< inciter::EOS >& mat_blk,
844 : : int intsharp,
845 : : std::size_t ncomp,
846 : : std::size_t nprim,
847 : : std::size_t rdof,
848 : : std::size_t nmat,
849 : : std::size_t e,
850 : : std::size_t dof_e,
851 : : const std::vector< std::size_t >& inpoel,
852 : : const UnsMesh::Coords& coord,
853 : : const Fields& geoElem,
854 : : const std::array< real, 3 >& ref_gp,
855 : : const std::vector< real >& B,
856 : : const Fields& U,
857 : : const Fields& P )
858 : : // *****************************************************************************
859 : : // Evaluate polynomial solution at quadrature point
860 : : //! \param[in] mat_blk EOS material block
861 : : //! \param[in] intsharp Interface reconstruction indicator
862 : : //! \param[in] ncomp Number of components in the PDE system
863 : : //! \param[in] nprim Number of primitive quantities
864 : : //! \param[in] rdof Total number of reconstructed dofs
865 : : //! \param[in] nmat Total number of materials
866 : : //! \param[in] e Element for which polynomial solution is being evaluated
867 : : //! \param[in] dof_e Degrees of freedom for element
868 : : //! \param[in] inpoel Element-node connectivity
869 : : //! \param[in] coord Array of nodal coordinates
870 : : //! \param[in] geoElem Element geometry array
871 : : //! \param[in] ref_gp Quadrature point in reference space
872 : : //! \param[in] B Basis function at given quadrature point
873 : : //! \param[in] U Solution vector
874 : : //! \param[in] P Vector of primitives
875 : : //! \return High-order unknown/state vector at quadrature point, modified
876 : : //! if near interfaces using THINC
877 : : // *****************************************************************************
878 : : {
879 : 92378342 : std::vector< real > state;
880 : 184756684 : std::vector< real > sprim;
881 : :
882 [ + - ]: 92378342 : state = eval_state( ncomp, rdof, dof_e, e, U, B );
883 [ + - ]: 92378342 : sprim = eval_state( nprim, rdof, dof_e, e, P, B );
884 : :
885 : : // interface detection
886 [ + - ]: 184756684 : std::vector< std::size_t > matInt(nmat, 0);
887 : 92378342 : bool intInd(false);
888 [ + + ]: 92378342 : if (nmat > 1) {
889 [ + - ]: 21687110 : std::vector< tk::real > alAvg(nmat, 0.0);
890 [ + + ]: 70611090 : for (std::size_t k=0; k<nmat; ++k)
891 [ + - ]: 48923980 : alAvg[k] = U(e, inciter::volfracDofIdx(nmat,k,rdof,0));
892 [ + - ]: 21687110 : intInd = inciter::interfaceIndicator(nmat, alAvg, matInt);
893 : : }
894 : :
895 : : // consolidate primitives into state vector
896 [ + - ]: 92378342 : state.insert(state.end(), sprim.begin(), sprim.end());
897 : :
898 [ + + ]: 92378342 : if (intsharp > 0)
899 : : {
900 [ + - ][ + - ]: 12411000 : std::vector< tk::real > vfmax(nmat, 0.0), vfmin(nmat, 0.0);
901 : :
902 : : // Until the appropriate setup for activating THINC with Transport
903 : : // is ready, the following two chunks of code will need to be commented
904 : : // for using THINC with Transport
905 : : //for (std::size_t k=0; k<nmat; ++k) {
906 : : // vfmin[k] = VolFracMax(el, 2*k, 0);
907 : : // vfmax[k] = VolFracMax(el, 2*k+1, 0);
908 : : //}
909 [ + - ]: 6205500 : tk::THINCReco(rdof, nmat, e, inpoel, coord, geoElem,
910 : : ref_gp, U, P, intInd, matInt, vfmin, vfmax, state);
911 : :
912 : : // Until the appropriate setup for activating THINC with Transport
913 : : // is ready, the following lines will need to be uncommented for
914 : : // using THINC with Transport
915 : : //tk::THINCRecoTransport(rdof, nmat, el, inpoel, coord,
916 : : // geoElem, ref_gp_l, U, P, vfmin, vfmax, state[0]);
917 : : }
918 : :
919 : : // physical constraints
920 [ + - ]: 92378342 : enforcePhysicalConstraints(mat_blk, ncomp, nmat, state);
921 : :
922 : 184756684 : return state;
923 : : }
924 : :
925 : : std::vector< tk::real >
926 : 1091856 : evalFVSol( const std::vector< inciter::EOS >& mat_blk,
927 : : int intsharp,
928 : : std::size_t ncomp,
929 : : std::size_t nprim,
930 : : std::size_t rdof,
931 : : std::size_t nmat,
932 : : std::size_t e,
933 : : const std::vector< std::size_t >& inpoel,
934 : : const UnsMesh::Coords& coord,
935 : : const Fields& geoElem,
936 : : const std::array< real, 3 >& ref_gp,
937 : : const std::vector< real >& B,
938 : : const Fields& U,
939 : : const Fields& P,
940 : : int srcFlag )
941 : : // *****************************************************************************
942 : : // Evaluate second-order FV solution at quadrature point
943 : : //! \param[in] mat_blk EOS material block
944 : : //! \param[in] intsharp Interface reconstruction indicator
945 : : //! \param[in] ncomp Number of components in the PDE system
946 : : //! \param[in] nprim Number of primitive quantities
947 : : //! \param[in] rdof Total number of reconstructed dofs
948 : : //! \param[in] nmat Total number of materials
949 : : //! \param[in] e Element for which polynomial solution is being evaluated
950 : : //! \param[in] inpoel Element-node connectivity
951 : : //! \param[in] coord Array of nodal coordinates
952 : : //! \param[in] geoElem Element geometry array
953 : : //! \param[in] ref_gp Quadrature point in reference space
954 : : //! \param[in] B Basis function at given quadrature point
955 : : //! \param[in] U Solution vector
956 : : //! \param[in] P Vector of primitives
957 : : //! \param[in] srcFlag Whether the energy source was added to element e
958 : : //! \return High-order unknown/state vector at quadrature point, modified
959 : : //! if near interfaces using THINC
960 : : // *****************************************************************************
961 : : {
962 : : auto min_al = inciter::g_inputdeck.get< tag::multimat,
963 : 1091856 : tag::min_volumefrac >();
964 : :
965 : : using inciter::pressureIdx;
966 : : using inciter::velocityIdx;
967 : : using inciter::volfracIdx;
968 : : using inciter::densityIdx;
969 : : using inciter::energyIdx;
970 : : using inciter::momentumIdx;
971 : :
972 : 1091856 : std::vector< real > state;
973 : 2183712 : std::vector< real > sprim;
974 : :
975 [ + - ]: 1091856 : state = eval_state( ncomp, rdof, rdof, e, U, B );
976 [ + - ]: 1091856 : sprim = eval_state( nprim, rdof, rdof, e, P, B );
977 : :
978 : : // interface detection so that eos is called on the appropriate quantities
979 [ + - ]: 2183712 : std::vector< std::size_t > matInt(nmat, 0);
980 [ + - ]: 2183712 : std::vector< tk::real > alAvg(nmat, 0.0);
981 [ + + ]: 3512624 : for (std::size_t k=0; k<nmat; ++k)
982 [ + - ]: 2420768 : alAvg[k] = U(e, inciter::volfracDofIdx(nmat,k,rdof,0));
983 [ + - ]: 1091856 : auto intInd = inciter::interfaceIndicator(nmat, alAvg, matInt);
984 : :
985 : : std::array< tk::real, 3 > vel{{
986 : 1091856 : sprim[velocityIdx(nmat,0)],
987 : 1091856 : sprim[velocityIdx(nmat,1)],
988 : 2183712 : sprim[velocityIdx(nmat,2)] }};
989 : :
990 : : // get mat-energy from reconstructed mat-pressure
991 : 1091856 : auto rhob(0.0);
992 [ + + ]: 3512624 : for (std::size_t k=0; k<nmat; ++k) {
993 : 2420768 : auto alk = state[volfracIdx(nmat,k)];
994 [ + + ]: 2420768 : if (matInt[k]) {
995 : 95028 : alk = std::max(std::min(alk, 1.0-static_cast<tk::real>(nmat-1)*min_al),
996 : 95028 : min_al);
997 : : }
998 : 2420768 : state[energyIdx(nmat,k)] =
999 [ + - ]: 4841536 : mat_blk[k].compute< inciter::EOS::totalenergy >(
1000 : 2420768 : state[densityIdx(nmat,k)], vel[0],
1001 : 2420768 : vel[1], vel[2],
1002 : 2420768 : sprim[pressureIdx(nmat,k)], alk);
1003 : 2420768 : rhob += state[densityIdx(nmat,k)];
1004 : : }
1005 : : // get momentum from reconstructed velocity and bulk density
1006 [ + + ]: 4367424 : for (std::size_t i=0; i<3; ++i) {
1007 : 3275568 : state[momentumIdx(nmat,i)] = rhob * vel[i];
1008 : : }
1009 : :
1010 : : // consolidate primitives into state vector
1011 [ + - ]: 1091856 : state.insert(state.end(), sprim.begin(), sprim.end());
1012 : :
1013 [ + + ][ + - ]: 1091856 : if (intsharp > 0 && srcFlag == 0)
1014 : : {
1015 [ + - ][ + - ]: 1103200 : std::vector< tk::real > vfmax(nmat, 0.0), vfmin(nmat, 0.0);
1016 : :
1017 [ + - ]: 551600 : tk::THINCReco(rdof, nmat, e, inpoel, coord, geoElem,
1018 : : ref_gp, U, P, intInd, matInt, vfmin, vfmax, state);
1019 : : }
1020 : :
1021 : : // physical constraints
1022 [ + - ]: 1091856 : enforcePhysicalConstraints(mat_blk, ncomp, nmat, state);
1023 : :
1024 : 2183712 : return state;
1025 : : }
1026 : :
1027 : : void
1028 : 93470198 : enforcePhysicalConstraints(
1029 : : const std::vector< inciter::EOS >& mat_blk,
1030 : : std::size_t ncomp,
1031 : : std::size_t nmat,
1032 : : std::vector< tk::real >& state )
1033 : : // *****************************************************************************
1034 : : // Enforce physical constraints on state at quadrature point
1035 : : //! \param[in] mat_blk EOS material block
1036 : : //! \param[in] ncomp Number of components in the PDE system
1037 : : //! \param[in] nmat Total number of materials
1038 : : //! \param[in,out] state state at quadrature point
1039 : : // *****************************************************************************
1040 : : {
1041 : 93470198 : auto myPDE = inciter::g_inputdeck.get< tag::pde >();
1042 : :
1043 : : // unfortunately have to query PDEType here. alternative will potentially
1044 : : // require refactor that passes PDEType from DGPDE to this level.
1045 [ + + ]: 93470198 : if (myPDE == inciter::ctr::PDEType::MULTIMAT) {
1046 : : using inciter::pressureIdx;
1047 : : using inciter::volfracIdx;
1048 : : using inciter::densityIdx;
1049 : :
1050 [ + + ]: 74123714 : for (std::size_t k=0; k<nmat; ++k) {
1051 : 51344748 : state[ncomp+pressureIdx(nmat,k)] = constrain_pressure( mat_blk,
1052 : 51344748 : state[ncomp+pressureIdx(nmat,k)], state[densityIdx(nmat,k)],
1053 : 51344748 : state[volfracIdx(nmat,k)], k );
1054 : : }
1055 : : }
1056 : : else if (myPDE == inciter::ctr::PDEType::MULTISPECIES) {
1057 : : // TODO: consider clipping temperature here
1058 : : }
1059 : 93470198 : }
1060 : :
1061 : : void
1062 : 0 : safeReco( std::size_t rdof,
1063 : : std::size_t nmat,
1064 : : std::size_t el,
1065 : : int er,
1066 : : const Fields& U,
1067 : : std::array< std::vector< real >, 2 >& state )
1068 : : // *****************************************************************************
1069 : : // Compute safe reconstructions near material interfaces
1070 : : //! \param[in] rdof Total number of reconstructed dofs
1071 : : //! \param[in] nmat Total number of material is PDE system
1072 : : //! \param[in] el Element on the left-side of face
1073 : : //! \param[in] er Element on the right-side of face
1074 : : //! \param[in] U Solution vector at recent time-stage
1075 : : //! \param[in,out] state Second-order reconstructed state, at cell-face, that
1076 : : //! is being modified for safety
1077 : : //! \details When the consistent limiting is applied, there is a possibility
1078 : : //! that the material densities and energies violate TVD bounds. This function
1079 : : //! enforces the TVD bounds locally
1080 : : // *****************************************************************************
1081 : : {
1082 : : using inciter::densityIdx;
1083 : : using inciter::energyIdx;
1084 : : using inciter::densityDofIdx;
1085 : : using inciter::energyDofIdx;
1086 : :
1087 [ - - ][ - - ]: 0 : if (er < 0) Throw("safe limiting cannot be called for boundary cells");
[ - - ][ - - ]
1088 : :
1089 : 0 : auto eR = static_cast< std::size_t >(er);
1090 : :
1091 : : // define a lambda for the safe limiting
1092 : 0 : auto safeLimit = [&]( std::size_t c, real ul, real ur )
1093 : : {
1094 : : // find min/max at the face
1095 : 0 : auto uMin = std::min(ul, ur);
1096 : 0 : auto uMax = std::max(ul, ur);
1097 : :
1098 : : // left-state limiting
1099 : 0 : state[0][c] = std::min(uMax, std::max(uMin, state[0][c]));
1100 : :
1101 : : // right-state limiting
1102 : 0 : state[1][c] = std::min(uMax, std::max(uMin, state[1][c]));
1103 : 0 : };
1104 : :
1105 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1106 : : {
1107 : 0 : real ul(0.0), ur(0.0);
1108 : :
1109 : : // Density
1110 : : // establish left- and right-hand states
1111 [ - - ]: 0 : ul = U(el, densityDofIdx(nmat, k, rdof, 0));
1112 [ - - ]: 0 : ur = U(eR, densityDofIdx(nmat, k, rdof, 0));
1113 : :
1114 : : // limit reconstructed density
1115 [ - - ]: 0 : safeLimit(densityIdx(nmat,k), ul, ur);
1116 : :
1117 : : // Energy
1118 : : // establish left- and right-hand states
1119 [ - - ]: 0 : ul = U(el, energyDofIdx(nmat, k, rdof, 0));
1120 [ - - ]: 0 : ur = U(eR, energyDofIdx(nmat, k, rdof, 0));
1121 : :
1122 : : // limit reconstructed energy
1123 [ - - ]: 0 : safeLimit(energyIdx(nmat,k), ul, ur);
1124 : : }
1125 : 0 : }
1126 : :
1127 : : } // tk::
|
