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