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