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 : 7050 : lhsLeastSq_P0P1( const inciter::FaceData& fd,
36 : : const Fields& geoElem,
37 : : const Fields& geoFace,
38 : : std::vector< std::array< std::array< real, 3 >, 3 > >& lhs_ls )
39 : : // *****************************************************************************
40 : : // Compute lhs matrix for the least-squares reconstruction
41 : : //! \param[in] fd Face connectivity and boundary conditions object
42 : : //! \param[in] geoElem Element geometry array
43 : : //! \param[in] geoFace Face geometry array
44 : : //! \param[in,out] lhs_ls LHS reconstruction matrix
45 : : //! \details This function computing the lhs matrix for reconstruction, is
46 : : //! common for primitive and conserved quantities.
47 : : // *****************************************************************************
48 : : {
49 : : const auto& esuf = fd.Esuf();
50 : 7050 : const auto nelem = fd.Esuel().size()/4;
51 : :
52 : : // Compute internal and boundary face contributions
53 [ + + ]: 4968450 : for (std::size_t f=0; f<esuf.size()/2; ++f)
54 : : {
55 : : Assert( esuf[2*f] > -1, "Left-side element detected as -1" );
56 : :
57 [ + + ]: 4961400 : auto el = static_cast< std::size_t >(esuf[2*f]);
58 : 4961400 : auto er = esuf[2*f+1];
59 : :
60 : : std::array< real, 3 > geoElemR;
61 : : std::size_t eR(0);
62 : :
63 : : // A second-order (piecewise linear) solution polynomial can be obtained
64 : : // from the first-order (piecewise constant) FV solutions by using a
65 : : // least-squares (LS) reconstruction process. LS uses the first-order
66 : : // solutions from the cell being processed, and the cells surrounding it.
67 : : // The LS system is obtaining by requiring the following to hold:
68 : : // 'Taylor expansions of solution from cell-i to the centroids of each of
69 : : // its neighboring cells should be equal to the cell average solution on
70 : : // that neighbor cell.'
71 : : // This gives a system of equations for the three second-order DOFs that are
72 : : // to be determined. In 3D tetrahedral meshes, this would give four
73 : : // equations (one for each neighbor )for the three unknown DOFs. This
74 : : // overdetermined system is solved in the least-squares sense using the
75 : : // normal equations approach. The normal equations approach involves
76 : : // pre-multiplying the overdetermined system by the transpose of the system
77 : : // matrix to obtain a square matrix (3x3 in this case).
78 : :
79 : : // get a 3x3 system by applying the normal equation approach to the
80 : : // least-squares overdetermined system
81 : :
82 [ + + ]: 4961400 : if (er > -1) {
83 : : // internal face contribution
84 : 3643500 : eR = static_cast< std::size_t >(er);
85 : : // Put in cell-centroid coordinates
86 : 3643500 : geoElemR = {{ geoElem(eR,1), geoElem(eR,2), geoElem(eR,3) }};
87 : : }
88 : : else {
89 : : // boundary face contribution
90 : : // Put in face-centroid coordinates
91 : 1317900 : geoElemR = {{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
92 : : }
93 : :
94 : 4961400 : std::array< real, 3 > wdeltax{{ geoElemR[0]-geoElem(el,1),
95 : 4961400 : geoElemR[1]-geoElem(el,2),
96 : 4961400 : geoElemR[2]-geoElem(el,3) }};
97 : :
98 : : // define a lambda for contributing to lhs matrix
99 : : auto lhs = [&]( std::size_t e ){
100 [ + + ][ + + ]: 33506400 : for (std::size_t idir=0; idir<3; ++idir)
101 [ + + ][ + + ]: 100519200 : for (std::size_t jdir=0; jdir<3; ++jdir)
102 : 75389400 : lhs_ls[e][idir][jdir] += wdeltax[idir] * wdeltax[jdir];
103 : : };
104 : :
105 : : // always add left element contribution (at a boundary face, the internal
106 : : // element is always the left element)
107 : : lhs(el);
108 : : // add right element contribution for internal faces only
109 [ + + ]: 4961400 : if (er > -1)
110 [ + + ]: 3643500 : if (eR < nelem) lhs(eR);
111 : :
112 : : }
113 : 7050 : }
114 : :
115 : : void
116 : 7050 : intLeastSq_P0P1( const std::size_t rdof,
117 : : const inciter::FaceData& fd,
118 : : const Fields& geoElem,
119 : : const Fields& W,
120 : : std::vector< std::vector< std::array< real, 3 > > >& rhs_ls,
121 : : const std::vector< std::size_t >& varList )
122 : : // *****************************************************************************
123 : : // \brief Compute internal surface contributions to rhs vector of the
124 : : // least-squares reconstruction
125 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
126 : : //! \param[in] fd Face connectivity and boundary conditions object
127 : : //! \param[in] geoElem Element geometry array
128 : : //! \param[in] W Solution vector to be reconstructed at recent time step
129 : : //! \param[in,out] rhs_ls RHS reconstruction vector
130 : : //! \param[in] varList List of indices in W, that need to be reconstructed
131 : : //! \details This function computing the internal face contributions to the rhs
132 : : //! vector for reconstruction, is common for primitive and conserved
133 : : //! quantities. If `W` == `U`, compute internal face contributions for the
134 : : //! conserved variables. If `W` == `P`, compute internal face contributions
135 : : //! for the primitive variables.
136 : : // *****************************************************************************
137 : : {
138 : : const auto& esuf = fd.Esuf();
139 : 7050 : const auto nelem = fd.Esuel().size()/4;
140 : :
141 : : // Compute internal face contributions
142 [ + + ]: 3650550 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
143 : : {
144 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
145 : : "as -1" );
146 : :
147 : 3643500 : auto el = static_cast< std::size_t >(esuf[2*f]);
148 : 3643500 : auto er = static_cast< std::size_t >(esuf[2*f+1]);
149 : :
150 : : // get a 3x3 system by applying the normal equation approach to the
151 : : // least-squares overdetermined system
152 : :
153 : : // 'wdeltax' is the distance vector between the centroids of this element
154 : : // and its neighbor
155 : 3643500 : std::array< real, 3 > wdeltax{{ geoElem(er,1)-geoElem(el,1),
156 : 3643500 : geoElem(er,2)-geoElem(el,2),
157 : 3643500 : geoElem(er,3)-geoElem(el,3) }};
158 : :
159 [ + + ]: 14574000 : for (std::size_t idir=0; idir<3; ++idir)
160 : : {
161 : : // rhs vector
162 [ + + ]: 31343400 : for (std::size_t i=0; i<varList.size(); ++i)
163 : : {
164 : 20412900 : auto c = varList[i];
165 : 20412900 : auto mark = c*rdof;
166 [ + + ]: 20412900 : rhs_ls[el][c][idir] +=
167 [ + + ]: 20412900 : wdeltax[idir] * (W(er,mark)-W(el,mark));
168 [ + + ]: 20412900 : if (er < nelem)
169 : 19728000 : rhs_ls[er][c][idir] +=
170 : 19728000 : wdeltax[idir] * (W(er,mark)-W(el,mark));
171 : : }
172 : : }
173 : : }
174 : 7050 : }
175 : :
176 : : void
177 : 42300 : bndLeastSqConservedVar_P0P1(
178 : : ncomp_t ncomp,
179 : : const std::vector< inciter::EOS >& mat_blk,
180 : : std::size_t rdof,
181 : : const std::vector< std::size_t >& bcconfig,
182 : : const inciter::FaceData& fd,
183 : : const Fields& geoFace,
184 : : const Fields& geoElem,
185 : : real t,
186 : : const StateFn& state,
187 : : const Fields& P,
188 : : const Fields& U,
189 : : std::vector< std::vector< std::array< real, 3 > > >& rhs_ls,
190 : : const std::vector< std::size_t >& varList,
191 : : std::size_t nprim )
192 : : // *****************************************************************************
193 : : // \brief Compute boundary surface contributions to rhs vector of the
194 : : // least-squares reconstruction of conserved quantities of the PDE system
195 : : //! \param[in] ncomp Number of scalar components in this PDE system
196 : : //! \param[in] mat_blk EOS material block
197 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
198 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
199 : : //! \param[in] fd Face connectivity and boundary conditions object
200 : : //! \param[in] geoFace Face geometry array
201 : : //! \param[in] geoElem Element geometry array
202 : : //! \param[in] t Physical time
203 : : //! \param[in] state Function to evaluate the left and right solution state at
204 : : //! boundaries
205 : : //! \param[in] P Primitive vector to be reconstructed at recent time step
206 : : //! \param[in] U Solution vector to be reconstructed at recent time step
207 : : //! \param[in,out] rhs_ls RHS reconstruction vector
208 : : //! \param[in] varList List of indices in W, that need to be reconstructed
209 : : //! \param[in] nprim This is the number of primitive quantities stored for this
210 : : //! PDE system. This is necessary to extend the state vector to the right
211 : : //! size, so that correct boundary conditions are obtained.
212 : : //! A default is set to 0, so that calling code for systems that do not store
213 : : //! primitive quantities does not need to specify this argument.
214 : : //! \details This function computing the boundary face contributions to the rhs
215 : : //! vector for reconstruction, is used for conserved quantities only.
216 : : // *****************************************************************************
217 : : {
218 : : const auto& bface = fd.Bface();
219 : : const auto& esuf = fd.Esuf();
220 : :
221 [ + + ]: 64350 : for (const auto& s : bcconfig) { // for all bc sidesets
222 : 22050 : auto bc = bface.find(static_cast<int>(s));// faces for side set
223 [ + + ]: 22050 : if (bc != end(bface))
224 : : {
225 : : // Compute boundary face contributions
226 [ + + ]: 1331550 : for (const auto& f : bc->second)
227 : : {
228 : : Assert( esuf[2*f+1] == -1, "physical boundary element not -1" );
229 : :
230 [ + - ]: 1317900 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
231 : :
232 : : // arrays for quadrature points
233 : : std::array< real, 3 >
234 : 1317900 : fc{{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
235 : : std::array< real, 3 >
236 : 1317900 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
237 : :
238 : : // Compute the state variables at the left element
239 [ + - ]: 1317900 : std::vector< real >B(1,1.0);
240 [ + - ]: 1317900 : auto ul = eval_state( ncomp, rdof, 1, el, U, B );
241 [ + - ]: 1317900 : auto uprim = eval_state( nprim, rdof, 1, el, P, B );
242 : :
243 : : // consolidate primitives into state vector
244 [ + - ][ - + ]: 1317900 : ul.insert(ul.end(), uprim.begin(), uprim.end());
[ - - ]
245 : :
246 : : Assert( ul.size() == ncomp+nprim, "Incorrect size for "
247 : : "appended state vector" );
248 : :
249 : : // Compute the state at the face-center using BC
250 : 1317900 : auto ustate = state( ncomp, mat_blk, ul, fc[0], fc[1], fc[2], t, fn );
251 : :
252 : 1317900 : std::array< real, 3 > wdeltax{{ fc[0]-geoElem(el,1),
253 : 1317900 : fc[1]-geoElem(el,2),
254 : 1317900 : fc[2]-geoElem(el,3) }};
255 : :
256 [ + + ]: 5271600 : for (std::size_t idir=0; idir<3; ++idir)
257 : : {
258 : : // rhs vector
259 [ + + ]: 10773000 : for (std::size_t i=0; i<varList.size(); ++i)
260 : : {
261 : 6819300 : auto c = varList[i];
262 : 6819300 : rhs_ls[el][c][idir] +=
263 : 6819300 : wdeltax[idir] * (ustate[1][c]-ustate[0][c]);
264 : : }
265 : : }
266 : : }
267 : : }
268 : : }
269 : 42300 : }
270 : :
271 : : void
272 : 7050 : solveLeastSq_P0P1(
273 : : const std::size_t rdof,
274 : : const std::vector< std::array< std::array< real, 3 >, 3 > >& lhs,
275 : : const std::vector< std::vector< std::array< real, 3 > > >& rhs,
276 : : Fields& W,
277 : : const std::vector< std::size_t >& varList )
278 : : // *****************************************************************************
279 : : // Solve the 3x3 linear system for least-squares reconstruction
280 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
281 : : //! \param[in] lhs LHS reconstruction matrix
282 : : //! \param[in] rhs RHS reconstruction vector
283 : : //! \param[in,out] W Solution vector to be reconstructed at recent time step
284 : : //! \param[in] varList List of indices in W, that need to be reconstructed
285 : : //! \details Solves the 3x3 linear system for each element, individually. For
286 : : //! systems that require reconstructions of primitive quantities, this should
287 : : //! be called twice, once with the argument 'W' as U (conserved), and again
288 : : //! with 'W' as P (primitive).
289 : : // *****************************************************************************
290 : : {
291 : 7050 : auto nelem = lhs.size();
292 : :
293 [ + + ]: 2101200 : for (std::size_t e=0; e<nelem; ++e)
294 : : {
295 [ + + ]: 6007500 : for (std::size_t i=0; i<varList.size(); ++i)
296 : : {
297 : 3913350 : auto mark = varList[i]*rdof;
298 : :
299 : : // solve system using Cramer's rule
300 : 3913350 : auto ux = tk::cramer( lhs[e], rhs[e][varList[i]] );
301 : :
302 : 3913350 : W(e,mark+1) = ux[0];
303 : 3913350 : W(e,mark+2) = ux[1];
304 : 3913350 : W(e,mark+3) = ux[2];
305 : : }
306 : : }
307 : 7050 : }
308 : :
309 : : void
310 : 1996120 : recoLeastSqExtStencil(
311 : : std::size_t rdof,
312 : : std::size_t e,
313 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
314 : : const std::vector< std::size_t >& inpoel,
315 : : const Fields& geoElem,
316 : : Fields& W,
317 : : const std::vector< std::size_t >& varList )
318 : : // *****************************************************************************
319 : : // \brief Reconstruct the second-order solution using least-squares approach
320 : : // from an extended stencil involving the node-neighbors
321 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
322 : : //! \param[in] e Element whoes solution is being reconstructed
323 : : //! \param[in] esup Elements surrounding points
324 : : //! \param[in] inpoel Element-node connectivity
325 : : //! \param[in] geoElem Element geometry array
326 : : //! \param[in,out] W Solution vector to be reconstructed at recent time step
327 : : //! \param[in] varList List of indices in W, that need to be reconstructed
328 : : //! \details A second-order (piecewise linear) solution polynomial is obtained
329 : : //! from the first-order (piecewise constant) FV solutions by using a
330 : : //! least-squares (LS) reconstruction process. This LS reconstruction function
331 : : //! using the nodal-neighbors of a cell, to get an overdetermined system of
332 : : //! equations for the derivatives of the solution. This overdetermined system
333 : : //! is solved in the least-squares sense using the normal equations approach.
334 : : // *****************************************************************************
335 : : {
336 : : // lhs matrix
337 : : std::array< std::array< tk::real, 3 >, 3 >
338 : 1996120 : lhs_ls( {{ {{0.0, 0.0, 0.0}},
339 : : {{0.0, 0.0, 0.0}},
340 : : {{0.0, 0.0, 0.0}} }} );
341 : : // rhs matrix
342 : : std::vector< std::array< tk::real, 3 > >
343 : 1996120 : rhs_ls( varList.size(), {{ 0.0, 0.0, 0.0 }} );
344 : :
345 : : // loop over all nodes of the element e
346 [ + + ]: 9980600 : for (std::size_t lp=0; lp<4; ++lp)
347 : : {
348 : 7984480 : auto p = inpoel[4*e+lp];
349 : : const auto& pesup = cref_find(esup, p);
350 : :
351 : : // loop over all the elements surrounding this node p
352 [ + + ]: 133443160 : for (auto er : pesup)
353 : : {
354 : : // centroid distance
355 : 125458680 : std::array< real, 3 > wdeltax{{ geoElem(er,1)-geoElem(e,1),
356 : 125458680 : geoElem(er,2)-geoElem(e,2),
357 : 125458680 : geoElem(er,3)-geoElem(e,3) }};
358 : :
359 : : // contribute to lhs matrix
360 [ + + ]: 501834720 : for (std::size_t idir=0; idir<3; ++idir)
361 [ + + ]: 1505504160 : for (std::size_t jdir=0; jdir<3; ++jdir)
362 : 1129128120 : lhs_ls[idir][jdir] += wdeltax[idir] * wdeltax[jdir];
363 : :
364 : : // compute rhs matrix
365 [ + + ]: 870737640 : for (std::size_t i=0; i<varList.size(); i++)
366 : : {
367 : 745278960 : auto mark = varList[i]*rdof;
368 [ + + ]: 2981115840 : for (std::size_t idir=0; idir<3; ++idir)
369 : 2235836880 : rhs_ls[i][idir] +=
370 : 2235836880 : wdeltax[idir] * (W(er,mark)-W(e,mark));
371 : :
372 : : }
373 : : }
374 : : }
375 : :
376 : : // solve least-square normal equation system using Cramer's rule
377 [ + + ]: 14010040 : for (std::size_t i=0; i<varList.size(); i++)
378 : : {
379 : 12013920 : auto mark = varList[i]*rdof;
380 : :
381 : 12013920 : auto ux = tk::cramer( lhs_ls, rhs_ls[i] );
382 : :
383 : : // Update the P1 dofs with the reconstructioned gradients.
384 : : // Since this reconstruction does not affect the cell-averaged solution,
385 : : // W(e,mark+0) is unchanged.
386 : 12013920 : W(e,mark+1) = ux[0];
387 : 12013920 : W(e,mark+2) = ux[1];
388 : 12013920 : W(e,mark+3) = ux[2];
389 : : }
390 : 1996120 : }
391 : :
392 : : void
393 : 18116 : transform_P0P1( std::size_t rdof,
394 : : std::size_t nelem,
395 : : const std::vector< std::size_t >& inpoel,
396 : : const UnsMesh::Coords& coord,
397 : : Fields& W,
398 : : const std::vector< std::size_t >& varList )
399 : : // *****************************************************************************
400 : : // Transform the reconstructed P1-derivatives to the Dubiner dofs
401 : : //! \param[in] rdof Total number of reconstructed dofs
402 : : //! \param[in] nelem Total number of elements
403 : : //! \param[in] inpoel Element-node connectivity
404 : : //! \param[in] coord Array of nodal coordinates
405 : : //! \param[in,out] W Second-order reconstructed vector which gets transformed to
406 : : //! the Dubiner reference space
407 : : //! \param[in] varList List of indices in W, that need to be reconstructed
408 : : //! \details Since the DG solution (and the primitive quantities) are assumed to
409 : : //! be stored in the Dubiner space, this transformation from Taylor
410 : : //! coefficients to Dubiner coefficients is necessary.
411 : : // *****************************************************************************
412 : : {
413 : : const auto& cx = coord[0];
414 : : const auto& cy = coord[1];
415 : : const auto& cz = coord[2];
416 : :
417 [ + + ]: 4108386 : for (std::size_t e=0; e<nelem; ++e)
418 : : {
419 : : // Extract the element coordinates
420 : : std::array< std::array< real, 3>, 4 > coordel {{
421 : 4090270 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
422 : 4090270 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
423 : 4090270 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
424 : 4090270 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
425 : 4090270 : }};
426 : :
427 : : auto jacInv =
428 : 4090270 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
429 : :
430 : : // Compute the derivatives of basis function for DG(P1)
431 : 4090270 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
432 : :
433 [ + + ]: 20017540 : for (std::size_t i=0; i<varList.size(); ++i)
434 : : {
435 : 15927270 : auto mark = varList[i]*rdof;
436 : :
437 : : // solve system using Cramer's rule
438 : 15927270 : auto ux = tk::cramer( {{ {{dBdx[0][1], dBdx[0][2], dBdx[0][3]}},
439 : 15927270 : {{dBdx[1][1], dBdx[1][2], dBdx[1][3]}},
440 : 15927270 : {{dBdx[2][1], dBdx[2][2], dBdx[2][3]}} }},
441 : 15927270 : {{ W(e,mark+1),
442 : 15927270 : W(e,mark+2),
443 : 15927270 : W(e,mark+3) }} );
444 : :
445 : : // replace physical derivatives with transformed dofs
446 : 15927270 : W(e,mark+1) = ux[0];
447 : 15927270 : W(e,mark+2) = ux[1];
448 : 15927270 : W(e,mark+3) = ux[2];
449 : : }
450 : : }
451 : 18116 : }
452 : :
453 : : void
454 : 6474300 : THINCReco( std::size_t rdof,
455 : : std::size_t nmat,
456 : : std::size_t e,
457 : : const std::vector< std::size_t >& inpoel,
458 : : const UnsMesh::Coords& coord,
459 : : const Fields& geoElem,
460 : : const std::array< real, 3 >& ref_xp,
461 : : const Fields& U,
462 : : const Fields& P,
463 : : bool intInd,
464 : : const std::vector< std::size_t >& matInt,
465 : : [[maybe_unused]] const std::vector< real >& vfmin,
466 : : [[maybe_unused]] const std::vector< real >& vfmax,
467 : : std::vector< real >& state )
468 : : // *****************************************************************************
469 : : // Compute THINC reconstructions at quadrature point for multi-material flows
470 : : //! \param[in] rdof Total number of reconstructed dofs
471 : : //! \param[in] nmat Total number of materials
472 : : //! \param[in] e Element for which interface reconstruction is being calculated
473 : : //! \param[in] inpoel Element-node connectivity
474 : : //! \param[in] coord Array of nodal coordinates
475 : : //! \param[in] geoElem Element geometry array
476 : : //! \param[in] ref_xp Quadrature point in reference space
477 : : //! \param[in] U Solution vector
478 : : //! \param[in] P Vector of primitives
479 : : //! \param[in] intInd Boolean which indicates if the element contains a
480 : : //! material interface
481 : : //! \param[in] matInt Array indicating which material has an interface
482 : : //! \param[in] vfmin Vector containing min volume fractions for each material
483 : : //! in this cell
484 : : //! \param[in] vfmax Vector containing max volume fractions for each material
485 : : //! in this cell
486 : : //! \param[in,out] state Unknown/state vector at quadrature point, modified
487 : : //! if near interfaces using THINC
488 : : //! \details This function is an interface for the multimat PDEs that use the
489 : : //! algebraic multi-material THINC reconstruction. This particular function
490 : : //! should only be called for multimat.
491 : : // *****************************************************************************
492 : : {
493 : : using inciter::volfracDofIdx;
494 : : using inciter::densityDofIdx;
495 : : using inciter::momentumDofIdx;
496 : : using inciter::energyDofIdx;
497 : : using inciter::pressureDofIdx;
498 : : using inciter::velocityDofIdx;
499 : : using inciter::deformDofIdx;
500 : : using inciter::stressDofIdx;
501 : : using inciter::volfracIdx;
502 : : using inciter::densityIdx;
503 : : using inciter::momentumIdx;
504 : : using inciter::energyIdx;
505 : : using inciter::pressureIdx;
506 : : using inciter::velocityIdx;
507 : : using inciter::deformIdx;
508 : : using inciter::stressIdx;
509 : :
510 : : auto bparam = inciter::g_inputdeck.get< tag::multimat,
511 : 6474300 : tag::intsharp_param >();
512 : 6474300 : const auto ncomp = U.nprop()/rdof;
513 : : const auto& solidx = inciter::g_inputdeck.get< tag::matidxmap,
514 : : tag::solidx >();
515 : :
516 : : // Step-1: Perform THINC reconstruction
517 : : // create a vector of volume-fractions and pass it to the THINC function
518 : 6474300 : std::vector< real > alSol(rdof*nmat, 0.0);
519 [ + - ][ - - ]: 6474300 : std::vector< real > alReco(nmat, 0.0);
520 [ + + ]: 19422900 : for (std::size_t k=0; k<nmat; ++k) {
521 : 12948600 : auto mark = k*rdof;
522 [ + + ]: 64743000 : for (std::size_t i=0; i<rdof; ++i) {
523 : 51794400 : alSol[mark+i] = U(e, volfracDofIdx(nmat,k,rdof,i));
524 : : }
525 : : // initialize with TVD reconstructions which will be modified if near
526 : : // material interface
527 : 12948600 : alReco[k] = state[volfracIdx(nmat,k)];
528 : : }
529 [ + - ]: 6474300 : THINCFunction(rdof, nmat, e, inpoel, coord, ref_xp, geoElem(e,0), bparam,
530 : : alSol, intInd, matInt, alReco);
531 : :
532 : : // check reconstructed volfracs for positivity
533 : : bool neg_vf = false;
534 [ + + ]: 19422900 : for (std::size_t k=0; k<nmat; ++k) {
535 [ - + ][ - - ]: 12948600 : if (alReco[k] < 1e-16 && matInt[k] > 0) neg_vf = true;
536 : : }
537 [ + + ]: 19422900 : for (std::size_t k=0; k<nmat; ++k) {
538 [ - + ]: 12948600 : if (neg_vf) {
539 : : std::cout << "Material-id: " << k << std::endl;
540 [ - - ]: 0 : std::cout << "Volume-fraction: " << std::setprecision(18) << alReco[k]
541 : : << std::endl;
542 : : std::cout << "Cell-avg vol-frac: " << std::setprecision(18) <<
543 [ - - ]: 0 : U(e,volfracDofIdx(nmat,k,rdof,0)) << std::endl;
544 : : std::cout << "Material-interface? " << intInd << std::endl;
545 [ - - ]: 0 : std::cout << "Mat-k-involved? " << matInt[k] << std::endl;
546 : : }
547 : : }
548 [ - + ][ - - ]: 6474300 : if (neg_vf) Throw("Material has negative volume fraction after THINC "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
549 : : "reconstruction.");
550 : :
551 : : // Step-2: Perform consistent reconstruction on other conserved quantities
552 [ + + ]: 6474300 : if (intInd)
553 : : {
554 : : auto rhobCC(0.0), rhobHO(0.0);
555 [ + + ]: 932616 : for (std::size_t k=0; k<nmat; ++k)
556 : : {
557 : 621744 : auto alCC = U(e, volfracDofIdx(nmat,k,rdof,0));
558 [ + - ]: 621744 : alCC = std::max(1e-14, alCC);
559 : :
560 [ + - ]: 621744 : if (matInt[k])
561 : : {
562 [ + + ]: 621744 : state[volfracIdx(nmat,k)] = alReco[k];
563 : 621744 : state[densityIdx(nmat,k)] = alReco[k]
564 [ + + ]: 621744 : * U(e, densityDofIdx(nmat,k,rdof,0))/alCC;
565 : 621744 : state[energyIdx(nmat,k)] = alReco[k]
566 [ + + ]: 621744 : * U(e, energyDofIdx(nmat,k,rdof,0))/alCC;
567 : 621744 : state[ncomp+pressureIdx(nmat,k)] = alReco[k]
568 : 621744 : * P(e, pressureDofIdx(nmat,k,rdof,0))/alCC;
569 [ + + ]: 621744 : if (solidx[k] > 0) {
570 [ + + ]: 146880 : for (std::size_t i=0; i<3; ++i)
571 [ + + ]: 440640 : for (std::size_t j=0; j<3; ++j)
572 : 330480 : state[deformIdx(nmat,solidx[k],i,j)] =
573 : 330480 : U(e, deformDofIdx(nmat,solidx[k],i,j,rdof,0));
574 : :
575 [ + + ]: 257040 : for (std::size_t i=0; i<6; ++i)
576 : 220320 : state[ncomp+stressIdx(nmat,solidx[k],i)] = alReco[k]
577 : 220320 : * P(e, stressDofIdx(nmat,solidx[k],i,rdof,0))/alCC;
578 : : }
579 : : }
580 : :
581 : 621744 : rhobCC += U(e, densityDofIdx(nmat,k,rdof,0));
582 : 621744 : rhobHO += state[densityIdx(nmat,k)];
583 : : }
584 : :
585 : : // consistent reconstruction for bulk momentum
586 [ + + ]: 1243488 : for (std::size_t i=0; i<3; ++i)
587 : : {
588 : 932616 : state[momentumIdx(nmat,i)] = rhobHO
589 : 932616 : * U(e, momentumDofIdx(nmat,i,rdof,0))/rhobCC;
590 : 932616 : state[ncomp+velocityIdx(nmat,i)] =
591 : 932616 : P(e, velocityDofIdx(nmat,i,rdof,0));
592 : : }
593 : : }
594 : 6474300 : }
595 : :
596 : : void
597 : 0 : THINCRecoTransport( std::size_t rdof,
598 : : std::size_t,
599 : : std::size_t e,
600 : : const std::vector< std::size_t >& inpoel,
601 : : const UnsMesh::Coords& coord,
602 : : const Fields& geoElem,
603 : : const std::array< real, 3 >& ref_xp,
604 : : const Fields& U,
605 : : const Fields&,
606 : : [[maybe_unused]] const std::vector< real >& vfmin,
607 : : [[maybe_unused]] const std::vector< real >& vfmax,
608 : : std::vector< real >& state )
609 : : // *****************************************************************************
610 : : // Compute THINC reconstructions at quadrature point for transport
611 : : //! \param[in] rdof Total number of reconstructed dofs
612 : : //! \param[in] e Element for which interface reconstruction is being calculated
613 : : //! \param[in] inpoel Element-node connectivity
614 : : //! \param[in] coord Array of nodal coordinates
615 : : //! \param[in] geoElem Element geometry array
616 : : //! \param[in] ref_xp Quadrature point in reference space
617 : : //! \param[in] U Solution vector
618 : : //! \param[in] vfmin Vector containing min volume fractions for each material
619 : : //! in this cell
620 : : //! \param[in] vfmax Vector containing max volume fractions for each material
621 : : //! in this cell
622 : : //! \param[in,out] state Unknown/state vector at quadrature point, modified
623 : : //! if near interfaces using THINC
624 : : //! \details This function is an interface for the transport PDEs that use the
625 : : //! algebraic multi-material THINC reconstruction. This particular function
626 : : //! should only be called for transport.
627 : : // *****************************************************************************
628 : : {
629 : : auto bparam = inciter::g_inputdeck.get< tag::transport,
630 : 0 : tag::intsharp_param >();
631 : 0 : auto ncomp = U.nprop()/rdof;
632 : :
633 : : // interface detection
634 : 0 : std::vector< std::size_t > matInt(ncomp, 0);
635 [ - - ][ - - ]: 0 : std::vector< tk::real > alAvg(ncomp, 0.0);
636 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k)
637 : 0 : alAvg[k] = U(e, k*rdof);
638 [ - - ]: 0 : auto intInd = inciter::interfaceIndicator(ncomp, alAvg, matInt);
639 : :
640 : : // create a vector of volume-fractions and pass it to the THINC function
641 [ - - ][ - - ]: 0 : std::vector< real > alSol(rdof*ncomp, 0.0);
642 : : // initialize with TVD reconstructions (modified if near interface)
643 [ - - ]: 0 : auto alReco = state;
644 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k) {
645 : 0 : auto mark = k*rdof;
646 [ - - ]: 0 : for (std::size_t i=0; i<rdof; ++i) {
647 : 0 : alSol[mark+i] = U(e,mark+i);
648 : : }
649 : : }
650 [ - - ]: 0 : THINCFunction(rdof, ncomp, e, inpoel, coord, ref_xp, geoElem(e,0), bparam,
651 : : alSol, intInd, matInt, alReco);
652 : :
653 [ - - ]: 0 : state = alReco;
654 : 0 : }
655 : :
656 : : void
657 : 6474300 : THINCFunction( std::size_t rdof,
658 : : std::size_t nmat,
659 : : std::size_t e,
660 : : const std::vector< std::size_t >& inpoel,
661 : : const UnsMesh::Coords& coord,
662 : : const std::array< real, 3 >& ref_xp,
663 : : real vol,
664 : : real bparam,
665 : : const std::vector< real >& alSol,
666 : : bool intInd,
667 : : const std::vector< std::size_t >& matInt,
668 : : std::vector< real >& alReco )
669 : : // *****************************************************************************
670 : : // Old version of the Multi-Medium THINC reconstruction function
671 : : //! \param[in] rdof Total number of reconstructed dofs
672 : : //! \param[in] nmat Total number of materials
673 : : //! \param[in] e Element for which interface reconstruction is being calculated
674 : : //! \param[in] inpoel Element-node connectivity
675 : : //! \param[in] coord Array of nodal coordinates
676 : : //! \param[in] ref_xp Quadrature point in reference space
677 : : //! \param[in] vol Element volume
678 : : //! \param[in] bparam User specified Beta for THINC, from the input file
679 : : //! \param[in] alSol Volume fraction solution vector for element e
680 : : //! \param[in] intInd Interface indicator, true if e is interface element
681 : : //! \param[in] matInt Vector indicating materials which constitute interface
682 : : //! \param[in,out] alReco Unknown/state vector at quadrature point, which gets
683 : : //! modified if near interface using MM-THINC
684 : : //! \details This function computes the interface reconstruction using the
685 : : //! algebraic multi-material THINC reconstruction for each material at the
686 : : //! given (ref_xp) quadrature point. This function is based on the following:
687 : : //! Pandare A. K., Waltz J., & Bakosi J. (2021) Multi-Material Hydrodynamics
688 : : //! with Algebraic Sharp Interface Capturing. Computers & Fluids,
689 : : //! doi: https://doi.org/10.1016/j.compfluid.2020.104804.
690 : : //! This function will be removed after the newer version (see
691 : : //! THINCFunction_new) is sufficiently tested.
692 : : // *****************************************************************************
693 : : {
694 : : // determine number of materials with interfaces in this cell
695 : 6474300 : auto epsl(1e-4), epsh(1e-1), bred(1.25), bmod(bparam);
696 : : std::size_t nIntMat(0);
697 [ + + ]: 19422900 : for (std::size_t k=0; k<nmat; ++k)
698 : : {
699 [ + + ]: 12948600 : auto alk = alSol[k*rdof];
700 [ + + ]: 12948600 : if (alk > epsl)
701 : : {
702 : 6727655 : ++nIntMat;
703 [ + + ]: 6727655 : if ((alk > epsl) && (alk < epsh))
704 : 191415 : bmod = std::min(bmod,
705 [ + + ]: 312410 : (alk-epsl)/(epsh-epsl) * (bred - bparam) + bparam);
706 [ + - ]: 6536240 : else if (alk > epsh)
707 : 6536240 : bmod = bred;
708 : : }
709 : : }
710 : :
711 [ - + ]: 6474300 : if (nIntMat > 2) bparam = bmod;
712 : :
713 : : // compression parameter
714 : 6474300 : auto beta = bparam/std::cbrt(6.0*vol);
715 : :
716 [ + + ]: 6474300 : if (intInd)
717 : : {
718 : : // 1. Get unit normals to material interface
719 : :
720 : : // Compute Jacobian matrix for converting Dubiner dofs to derivatives
721 : : const auto& cx = coord[0];
722 : : const auto& cy = coord[1];
723 : : const auto& cz = coord[2];
724 : :
725 : : std::array< std::array< real, 3>, 4 > coordel {{
726 : 310872 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
727 : 310872 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
728 : 310872 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
729 : 310872 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
730 : 310872 : }};
731 : :
732 : : auto jacInv =
733 : 310872 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
734 : :
735 [ + - ]: 310872 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
736 : :
737 : : std::array< real, 3 > nInt;
738 [ - + ][ - - ]: 310872 : std::vector< std::array< real, 3 > > ref_n(nmat, {{0.0, 0.0, 0.0}});
739 : :
740 : : // Get normals
741 [ + + ]: 932616 : for (std::size_t k=0; k<nmat; ++k)
742 : : {
743 : : // Get derivatives from moments in Dubiner space
744 [ + + ]: 2486976 : for (std::size_t i=0; i<3; ++i)
745 : 1865232 : nInt[i] = dBdx[i][1] * alSol[k*rdof+1]
746 : 1865232 : + dBdx[i][2] * alSol[k*rdof+2]
747 : 1865232 : + dBdx[i][3] * alSol[k*rdof+3];
748 : :
749 : 621744 : auto nMag = std::sqrt(tk::dot(nInt, nInt)) + 1e-14;
750 : :
751 [ + + ]: 2486976 : for (std::size_t i=0; i<3; ++i)
752 : 1865232 : nInt[i] /= nMag;
753 : :
754 : : // project interface normal onto local/reference coordinate system
755 [ + + ]: 2486976 : for (std::size_t i=0; i<3; ++i)
756 : : {
757 : : std::array< real, 3 > axis{
758 : 1865232 : coordel[i+1][0]-coordel[0][0],
759 : 1865232 : coordel[i+1][1]-coordel[0][1],
760 : 1865232 : coordel[i+1][2]-coordel[0][2] };
761 : 1865232 : ref_n[k][i] = tk::dot(nInt, axis);
762 : : }
763 : : }
764 : :
765 : : // 2. Reconstruct volume fractions using THINC
766 : 310872 : auto max_lim = 1.0 - (static_cast<tk::real>(nmat-1)*1.0e-12);
767 : 310872 : auto min_lim = 1e-12;
768 : : auto sum_inter(0.0), sum_non_inter(0.0);
769 [ + + ]: 932616 : for (std::size_t k=0; k<nmat; ++k)
770 : : {
771 [ + - ]: 621744 : if (matInt[k])
772 : : {
773 : : // get location of material interface (volume fraction 0.5) from the
774 : : // assumed tanh volume fraction distribution, and cell-averaged
775 : : // volume fraction
776 [ + + ]: 621744 : auto alCC(alSol[k*rdof]);
777 : : auto Ac(0.0), Bc(0.0), Qc(0.0);
778 [ + + ]: 621744 : if ((std::abs(ref_n[k][0]) > std::abs(ref_n[k][1]))
779 [ + + ][ + + ]: 621744 : && (std::abs(ref_n[k][0]) > std::abs(ref_n[k][2])))
780 : : {
781 : 251930 : Ac = std::exp(0.5*beta*ref_n[k][0]);
782 : 251930 : Bc = std::exp(0.5*beta*(ref_n[k][1]+ref_n[k][2]));
783 : 251930 : Qc = std::exp(0.5*beta*ref_n[k][0]*(2.0*alCC-1.0));
784 : : }
785 : : else if ((std::abs(ref_n[k][1]) > std::abs(ref_n[k][0]))
786 [ + + ][ + + ]: 369814 : && (std::abs(ref_n[k][1]) > std::abs(ref_n[k][2])))
787 : : {
788 : 270266 : Ac = std::exp(0.5*beta*ref_n[k][1]);
789 : 270266 : Bc = std::exp(0.5*beta*(ref_n[k][0]+ref_n[k][2]));
790 : 270266 : Qc = std::exp(0.5*beta*ref_n[k][1]*(2.0*alCC-1.0));
791 : : }
792 : : else
793 : : {
794 : 99548 : Ac = std::exp(0.5*beta*ref_n[k][2]);
795 : 99548 : Bc = std::exp(0.5*beta*(ref_n[k][0]+ref_n[k][1]));
796 : 99548 : Qc = std::exp(0.5*beta*ref_n[k][2]*(2.0*alCC-1.0));
797 : : }
798 : 621744 : auto d = std::log((1.0-Ac*Qc) / (Ac*Bc*(Qc-Ac))) / (2.0*beta);
799 : :
800 : : // THINC reconstruction
801 [ + - ]: 621744 : auto al_c = 0.5 * (1.0 + std::tanh(beta*(tk::dot(ref_n[k], ref_xp) + d)));
802 : :
803 : 621744 : alReco[k] = std::min(max_lim, std::max(min_lim, al_c));
804 : :
805 : 621744 : sum_inter += alReco[k];
806 : : } else
807 : : {
808 : 0 : sum_non_inter += alReco[k];
809 : : }
810 : : // else, if this material does not have an interface close-by, the TVD
811 : : // reconstructions must be used for state variables. This is ensured by
812 : : // initializing the alReco vector as the TVD state.
813 : : }
814 : :
815 : : // Rescale volume fractions of interface-materials to ensure unit sum
816 : 310872 : auto sum_rest = 1.0 - sum_non_inter;
817 [ + + ]: 932616 : for (std::size_t k=0; k<nmat; ++k)
818 [ + - ]: 621744 : if(matInt[k])
819 : 621744 : alReco[k] = alReco[k] * sum_rest / sum_inter;
820 : : }
821 : 6474300 : }
822 : :
823 : : void
824 : 0 : THINCFunction_new( std::size_t rdof,
825 : : std::size_t nmat,
826 : : std::size_t e,
827 : : const std::vector< std::size_t >& inpoel,
828 : : const UnsMesh::Coords& coord,
829 : : const std::array< real, 3 >& ref_xp,
830 : : real vol,
831 : : real bparam,
832 : : const std::vector< real >& alSol,
833 : : bool intInd,
834 : : const std::vector< std::size_t >& matInt,
835 : : std::vector< real >& alReco )
836 : : // *****************************************************************************
837 : : // New Multi-Medium THINC reconstruction function for volume fractions
838 : : //! \param[in] rdof Total number of reconstructed dofs
839 : : //! \param[in] nmat Total number of materials
840 : : //! \param[in] e Element for which interface reconstruction is being calculated
841 : : //! \param[in] inpoel Element-node connectivity
842 : : //! \param[in] coord Array of nodal coordinates
843 : : //! \param[in] ref_xp Quadrature point in reference space
844 : : //! \param[in] vol Element volume
845 : : //! \param[in] bparam User specified Beta for THINC, from the input file
846 : : //! \param[in] alSol Volume fraction solution vector for element e
847 : : //! \param[in] intInd Interface indicator, true if e is interface element
848 : : //! \param[in] matInt Vector indicating materials which constitute interface
849 : : //! \param[in,out] alReco Unknown/state vector at quadrature point, which gets
850 : : //! modified if near interface using MM-THINC
851 : : //! \details This function computes the interface reconstruction using the
852 : : //! algebraic multi-material THINC reconstruction for each material at the
853 : : //! given (ref_xp) quadrature point. This function succeeds the older version
854 : : //! of the mm-THINC (see THINCFunction), but is still under testing and is
855 : : //! currently experimental.
856 : : // *****************************************************************************
857 : : {
858 : : // compression parameter
859 : 0 : auto beta = bparam/std::cbrt(6.0*vol);
860 : :
861 : : // If the cell is not material interface, return this function
862 [ - - ]: 0 : if (not intInd) return;
863 : :
864 : : // If the cell is material interface, THINC reconstruction is applied
865 : : // Step 1. Get unit normals to material interface
866 : : // -------------------------------------------------------------------------
867 : :
868 : : // Compute Jacobian matrix for converting Dubiner dofs to derivatives
869 : : const auto& cx = coord[0];
870 : : const auto& cy = coord[1];
871 : : const auto& cz = coord[2];
872 : :
873 : : std::array< std::array< real, 3>, 4 > coordel {{
874 : 0 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
875 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
876 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
877 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
878 : 0 : }};
879 : :
880 : : auto jacInv =
881 : 0 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
882 : :
883 : 0 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
884 : :
885 : : std::array< real, 3 > nInt;
886 : 0 : std::array< real, 3 > ref_n{0.0, 0.0, 0.0};
887 : : auto almax(0.0);
888 : : std::size_t kmax(0);
889 : :
890 : : // Determine index of material present in majority
891 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
892 : : {
893 [ - - ]: 0 : auto alk = alSol[k*rdof];
894 [ - - ]: 0 : if (alk > almax)
895 : : {
896 : : almax = alk;
897 : : kmax = k;
898 : : }
899 : : }
900 : :
901 : : // Get normals of material present in majority
902 : : // Get derivatives from moments in Dubiner space
903 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
904 : 0 : nInt[i] = dBdx[i][1] * alSol[kmax*rdof+1]
905 : 0 : + dBdx[i][2] * alSol[kmax*rdof+2]
906 : 0 : + dBdx[i][3] * alSol[kmax*rdof+3];
907 : :
908 : 0 : auto nMag = std::sqrt(tk::dot(nInt, nInt)) + 1e-14;
909 : :
910 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
911 : 0 : nInt[i] /= nMag;
912 : :
913 : : // project interface normal onto local/reference coordinate system
914 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
915 : : {
916 : : std::array< real, 3 > axis{
917 : 0 : coordel[i+1][0]-coordel[0][0],
918 : 0 : coordel[i+1][1]-coordel[0][1],
919 : 0 : coordel[i+1][2]-coordel[0][2] };
920 : 0 : ref_n[i] = tk::dot(nInt, axis);
921 : : }
922 : :
923 : : // Step 2. Reconstruct volume fraction of majority material using THINC
924 : : // -------------------------------------------------------------------------
925 : :
926 : 0 : auto al_max = 1.0 - (static_cast<tk::real>(nmat-1)*1.0e-12);
927 : 0 : auto al_min = 1e-12;
928 : : auto alsum(0.0);
929 : : // get location of material interface (volume fraction 0.5) from the
930 : : // assumed tanh volume fraction distribution, and cell-averaged
931 : : // volume fraction
932 [ - - ]: 0 : auto alCC(alSol[kmax*rdof]);
933 : : auto Ac(0.0), Bc(0.0), Qc(0.0);
934 [ - - ]: 0 : if ((std::abs(ref_n[0]) > std::abs(ref_n[1]))
935 [ - - ][ - - ]: 0 : && (std::abs(ref_n[0]) > std::abs(ref_n[2])))
936 : : {
937 : 0 : Ac = std::exp(0.5*beta*ref_n[0]);
938 : 0 : Bc = std::exp(0.5*beta*(ref_n[1]+ref_n[2]));
939 : 0 : Qc = std::exp(0.5*beta*ref_n[0]*(2.0*alCC-1.0));
940 : : }
941 : : else if ((std::abs(ref_n[1]) > std::abs(ref_n[0]))
942 [ - - ][ - - ]: 0 : && (std::abs(ref_n[1]) > std::abs(ref_n[2])))
943 : : {
944 : 0 : Ac = std::exp(0.5*beta*ref_n[1]);
945 : 0 : Bc = std::exp(0.5*beta*(ref_n[0]+ref_n[2]));
946 : 0 : Qc = std::exp(0.5*beta*ref_n[1]*(2.0*alCC-1.0));
947 : : }
948 : : else
949 : : {
950 : 0 : Ac = std::exp(0.5*beta*ref_n[2]);
951 : 0 : Bc = std::exp(0.5*beta*(ref_n[0]+ref_n[1]));
952 : 0 : Qc = std::exp(0.5*beta*ref_n[2]*(2.0*alCC-1.0));
953 : : }
954 [ - - ]: 0 : auto d = std::log((1.0-Ac*Qc) / (Ac*Bc*(Qc-Ac))) / (2.0*beta);
955 : :
956 : : // THINC reconstruction
957 [ - - ]: 0 : auto al_c = 0.5 * (1.0 + std::tanh(beta*(tk::dot(ref_n, ref_xp) + d)));
958 : :
959 : 0 : alReco[kmax] = std::min(al_max, std::max(al_min, al_c));
960 : 0 : alsum += alReco[kmax];
961 : :
962 : : // if this material does not have an interface close-by, the TVD
963 : : // reconstructions must be used for state variables. This is ensured by
964 : : // initializing the alReco vector as the TVD state.
965 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
966 [ - - ]: 0 : if (!matInt[k]) {
967 : 0 : alsum += alReco[k];
968 : : }
969 : : }
970 : :
971 : : // Step 3. Do multimaterial cell corrections
972 : : // -------------------------------------------------------------------------
973 : :
974 : : // distribute remaining volume to rest of materials
975 : 0 : auto sum_left = 1.0 - alsum;
976 : : real den = 0.0;
977 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
978 [ - - ][ - - ]: 0 : if (matInt[k] && k != kmax) {
979 : 0 : auto mark = k * rdof;
980 : 0 : alReco[k] = sum_left * alSol[mark];
981 : 0 : den += alSol[mark];
982 : : }
983 : : }
984 : : // the distributed volfracs might be below al_min, correct that
985 : : real err = 0.0;
986 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
987 [ - - ][ - - ]: 0 : if (matInt[k] && k != kmax) {
988 : 0 : alReco[k] /= den;
989 [ - - ]: 0 : if (alReco[k] < al_min) {
990 : 0 : err += al_min - alReco[k];
991 : 0 : alReco[k] = al_min;
992 : : }
993 : : }
994 : : }
995 : :
996 : : // balance out errors
997 : 0 : alReco[kmax] -= err;
998 : : }
999 : :
1000 : : std::vector< tk::real >
1001 [ + - ]: 191576292 : evalPolynomialSol( const std::vector< inciter::EOS >& mat_blk,
1002 : : int intsharp,
1003 : : std::size_t ncomp,
1004 : : std::size_t nprim,
1005 : : std::size_t rdof,
1006 : : std::size_t nmat,
1007 : : std::size_t e,
1008 : : std::size_t dof_e,
1009 : : const std::vector< std::size_t >& inpoel,
1010 : : const UnsMesh::Coords& coord,
1011 : : const Fields& geoElem,
1012 : : const std::array< real, 3 >& ref_gp,
1013 : : const std::vector< real >& B,
1014 : : const Fields& U,
1015 : : const Fields& P )
1016 : : // *****************************************************************************
1017 : : // Evaluate polynomial solution at quadrature point
1018 : : //! \param[in] mat_blk EOS material block
1019 : : //! \param[in] intsharp Interface reconstruction indicator
1020 : : //! \param[in] ncomp Number of components in the PDE system
1021 : : //! \param[in] nprim Number of primitive quantities
1022 : : //! \param[in] rdof Total number of reconstructed dofs
1023 : : //! \param[in] nmat Total number of materials
1024 : : //! \param[in] e Element for which polynomial solution is being evaluated
1025 : : //! \param[in] dof_e Degrees of freedom for element
1026 : : //! \param[in] inpoel Element-node connectivity
1027 : : //! \param[in] coord Array of nodal coordinates
1028 : : //! \param[in] geoElem Element geometry array
1029 : : //! \param[in] ref_gp Quadrature point in reference space
1030 : : //! \param[in] B Basis function at given quadrature point
1031 : : //! \param[in] U Solution vector
1032 : : //! \param[in] P Vector of primitives
1033 : : //! \return High-order unknown/state vector at quadrature point, modified
1034 : : //! if near interfaces using THINC
1035 : : // *****************************************************************************
1036 : : {
1037 : : std::vector< real > state;
1038 : : std::vector< real > sprim;
1039 : :
1040 [ + - ]: 191576292 : state = eval_state( ncomp, rdof, dof_e, e, U, B );
1041 [ + - ]: 191576292 : sprim = eval_state( nprim, rdof, dof_e, e, P, B );
1042 : :
1043 : : // interface detection
1044 [ + - ][ - - ]: 191576292 : std::vector< std::size_t > matInt(nmat, 0);
1045 : : bool intInd(false);
1046 [ + + ]: 191576292 : if (nmat > 1) {
1047 [ + - ]: 27490830 : std::vector< tk::real > alAvg(nmat, 0.0);
1048 [ + + ]: 88022250 : for (std::size_t k=0; k<nmat; ++k)
1049 : 60531420 : alAvg[k] = U(e, inciter::volfracDofIdx(nmat,k,rdof,0));
1050 [ + - ]: 27490830 : intInd = inciter::interfaceIndicator(nmat, alAvg, matInt);
1051 : : }
1052 : :
1053 : : // consolidate primitives into state vector
1054 [ + - ]: 191576292 : state.insert(state.end(), sprim.begin(), sprim.end());
1055 : :
1056 [ + + ]: 191576292 : if (intsharp > 0)
1057 : : {
1058 [ + - ][ + - ]: 6474300 : std::vector< tk::real > vfmax(nmat, 0.0), vfmin(nmat, 0.0);
[ - - ][ - - ]
1059 : :
1060 : : // Until the appropriate setup for activating THINC with Transport
1061 : : // is ready, the following two chunks of code will need to be commented
1062 : : // for using THINC with Transport
1063 : : //for (std::size_t k=0; k<nmat; ++k) {
1064 : : // vfmin[k] = VolFracMax(el, 2*k, 0);
1065 : : // vfmax[k] = VolFracMax(el, 2*k+1, 0);
1066 : : //}
1067 [ + - ]: 6474300 : tk::THINCReco(rdof, nmat, e, inpoel, coord, geoElem,
1068 : : ref_gp, U, P, intInd, matInt, vfmin, vfmax, state);
1069 : :
1070 : : // Until the appropriate setup for activating THINC with Transport
1071 : : // is ready, the following lines will need to be uncommented for
1072 : : // using THINC with Transport
1073 : : //tk::THINCRecoTransport(rdof, nmat, el, inpoel, coord,
1074 : : // geoElem, ref_gp_l, U, P, vfmin, vfmax, state[0]);
1075 : : }
1076 : :
1077 : : // physical constraints
1078 [ + + ]: 191576292 : if (state.size() > ncomp) {
1079 : : using inciter::pressureIdx;
1080 : : using inciter::volfracIdx;
1081 : : using inciter::densityIdx;
1082 : :
1083 [ + + ]: 88022250 : for (std::size_t k=0; k<nmat; ++k) {
1084 [ + - ]: 60531420 : state[ncomp+pressureIdx(nmat,k)] = constrain_pressure( mat_blk,
1085 [ + - ]: 60531420 : state[ncomp+pressureIdx(nmat,k)], state[densityIdx(nmat,k)],
1086 [ + - ]: 60531420 : state[volfracIdx(nmat,k)], k );
1087 : : }
1088 : : }
1089 : :
1090 : 191576292 : return state;
1091 : : }
1092 : :
1093 : : std::vector< tk::real >
1094 [ + - ]: 975256 : evalFVSol( const std::vector< inciter::EOS >& mat_blk,
1095 : : int intsharp,
1096 : : std::size_t ncomp,
1097 : : std::size_t nprim,
1098 : : std::size_t rdof,
1099 : : std::size_t nmat,
1100 : : std::size_t e,
1101 : : const std::vector< std::size_t >& inpoel,
1102 : : const UnsMesh::Coords& coord,
1103 : : const Fields& geoElem,
1104 : : const std::array< real, 3 >& ref_gp,
1105 : : const std::vector< real >& B,
1106 : : const Fields& U,
1107 : : const Fields& P,
1108 : : int srcFlag )
1109 : : // *****************************************************************************
1110 : : // Evaluate second-order FV solution at quadrature point
1111 : : //! \param[in] mat_blk EOS material block
1112 : : //! \param[in] intsharp Interface reconstruction indicator
1113 : : //! \param[in] ncomp Number of components in the PDE system
1114 : : //! \param[in] nprim Number of primitive quantities
1115 : : //! \param[in] rdof Total number of reconstructed dofs
1116 : : //! \param[in] nmat Total number of materials
1117 : : //! \param[in] e Element for which polynomial solution is being evaluated
1118 : : //! \param[in] inpoel Element-node connectivity
1119 : : //! \param[in] coord Array of nodal coordinates
1120 : : //! \param[in] geoElem Element geometry array
1121 : : //! \param[in] ref_gp Quadrature point in reference space
1122 : : //! \param[in] B Basis function at given quadrature point
1123 : : //! \param[in] U Solution vector
1124 : : //! \param[in] P Vector of primitives
1125 : : //! \param[in] srcFlag Whether the energy source was added to element e
1126 : : //! \return High-order unknown/state vector at quadrature point, modified
1127 : : //! if near interfaces using THINC
1128 : : // *****************************************************************************
1129 : : {
1130 : : using inciter::pressureIdx;
1131 : : using inciter::velocityIdx;
1132 : : using inciter::volfracIdx;
1133 : : using inciter::densityIdx;
1134 : : using inciter::energyIdx;
1135 : : using inciter::momentumIdx;
1136 : :
1137 : : std::vector< real > state;
1138 : : std::vector< real > sprim;
1139 : :
1140 [ + - ]: 975256 : state = eval_state( ncomp, rdof, rdof, e, U, B );
1141 [ + - ]: 975256 : sprim = eval_state( nprim, rdof, rdof, e, P, B );
1142 : :
1143 : : // interface detection so that eos is called on the appropriate quantities
1144 [ + - ][ - - ]: 975256 : std::vector< std::size_t > matInt(nmat, 0);
1145 [ + - ][ - - ]: 975256 : std::vector< tk::real > alAvg(nmat, 0.0);
1146 [ + + ]: 3162824 : for (std::size_t k=0; k<nmat; ++k)
1147 : 2187568 : alAvg[k] = U(e, inciter::volfracDofIdx(nmat,k,rdof,0));
1148 [ + - ]: 975256 : auto intInd = inciter::interfaceIndicator(nmat, alAvg, matInt);
1149 : :
1150 : : // get mat-energy from reconstructed mat-pressure
1151 : : auto rhob(0.0);
1152 [ + + ]: 3162824 : for (std::size_t k=0; k<nmat; ++k) {
1153 [ + + ]: 2187568 : auto alk = state[volfracIdx(nmat,k)];
1154 [ + + ]: 2187568 : if (matInt[k]) {
1155 : 1061584 : alk = std::max(std::min(alk, 1.0-static_cast<tk::real>(nmat-1)*1e-12),
1156 [ - + ]: 1061584 : 1e-12);
1157 : : }
1158 : 2187568 : state[energyIdx(nmat,k)] = alk *
1159 [ + - ]: 2187568 : mat_blk[k].compute< inciter::EOS::totalenergy >(
1160 [ + - ]: 2187568 : state[densityIdx(nmat,k)]/alk, sprim[velocityIdx(nmat,0)],
1161 : : sprim[velocityIdx(nmat,1)], sprim[velocityIdx(nmat,2)],
1162 [ + - ]: 2187568 : sprim[pressureIdx(nmat,k)]/alk);
1163 : 2187568 : rhob += state[densityIdx(nmat,k)];
1164 : : }
1165 : : // get momentum from reconstructed velocity and bulk density
1166 [ + + ]: 3901024 : for (std::size_t i=0; i<3; ++i) {
1167 : 2925768 : state[momentumIdx(nmat,i)] = rhob * sprim[velocityIdx(nmat,i)];
1168 : : }
1169 : :
1170 : : // consolidate primitives into state vector
1171 [ + - ]: 975256 : state.insert(state.end(), sprim.begin(), sprim.end());
1172 : :
1173 [ - + ]: 975256 : if (intsharp > 0 && srcFlag == 0)
1174 : : {
1175 [ - - ][ - - ]: 0 : std::vector< tk::real > vfmax(nmat, 0.0), vfmin(nmat, 0.0);
[ - - ][ - - ]
1176 : :
1177 [ - - ]: 0 : tk::THINCReco(rdof, nmat, e, inpoel, coord, geoElem,
1178 : : ref_gp, U, P, intInd, matInt, vfmin, vfmax, state);
1179 : : }
1180 : :
1181 : : // physical constraints
1182 [ + - ]: 975256 : if (state.size() > ncomp) {
1183 [ + + ]: 3162824 : for (std::size_t k=0; k<nmat; ++k) {
1184 [ + - ]: 2187568 : state[ncomp+pressureIdx(nmat,k)] = constrain_pressure( mat_blk,
1185 [ + - ]: 2187568 : state[ncomp+pressureIdx(nmat,k)], state[densityIdx(nmat,k)],
1186 [ + - ]: 2187568 : state[volfracIdx(nmat,k)], k );
1187 : : }
1188 : : }
1189 : :
1190 : 975256 : return state;
1191 : : }
1192 : :
1193 : : void
1194 : 0 : safeReco( std::size_t rdof,
1195 : : std::size_t nmat,
1196 : : std::size_t el,
1197 : : int er,
1198 : : const Fields& U,
1199 : : std::array< std::vector< real >, 2 >& state )
1200 : : // *****************************************************************************
1201 : : // Compute safe reconstructions near material interfaces
1202 : : //! \param[in] rdof Total number of reconstructed dofs
1203 : : //! \param[in] nmat Total number of material is PDE system
1204 : : //! \param[in] el Element on the left-side of face
1205 : : //! \param[in] er Element on the right-side of face
1206 : : //! \param[in] U Solution vector at recent time-stage
1207 : : //! \param[in,out] state Second-order reconstructed state, at cell-face, that
1208 : : //! is being modified for safety
1209 : : //! \details When the consistent limiting is applied, there is a possibility
1210 : : //! that the material densities and energies violate TVD bounds. This function
1211 : : //! enforces the TVD bounds locally
1212 : : // *****************************************************************************
1213 : : {
1214 : : using inciter::densityIdx;
1215 : : using inciter::energyIdx;
1216 : : using inciter::densityDofIdx;
1217 : : using inciter::energyDofIdx;
1218 : :
1219 [ - - ][ - - ]: 0 : if (er < 0) Throw("safe limiting cannot be called for boundary cells");
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1220 : :
1221 : 0 : auto eR = static_cast< std::size_t >(er);
1222 : :
1223 : : // define a lambda for the safe limiting
1224 : 0 : auto safeLimit = [&]( std::size_t c, real ul, real ur )
1225 : : {
1226 : : // find min/max at the face
1227 [ - - ]: 0 : auto uMin = std::min(ul, ur);
1228 : 0 : auto uMax = std::max(ul, ur);
1229 : :
1230 : : // left-state limiting
1231 [ - - ]: 0 : state[0][c] = std::min(uMax, std::max(uMin, state[0][c]));
1232 : :
1233 : : // right-state limiting
1234 [ - - ]: 0 : state[1][c] = std::min(uMax, std::max(uMin, state[1][c]));
1235 : 0 : };
1236 : :
1237 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1238 : : {
1239 : : real ul(0.0), ur(0.0);
1240 : :
1241 : : // Density
1242 : : // establish left- and right-hand states
1243 : 0 : ul = U(el, densityDofIdx(nmat, k, rdof, 0));
1244 : 0 : ur = U(eR, densityDofIdx(nmat, k, rdof, 0));
1245 : :
1246 : : // limit reconstructed density
1247 : 0 : safeLimit(densityIdx(nmat,k), ul, ur);
1248 : :
1249 : : // Energy
1250 : : // establish left- and right-hand states
1251 : 0 : ul = U(el, energyDofIdx(nmat, k, rdof, 0));
1252 : 0 : ur = U(eR, energyDofIdx(nmat, k, rdof, 0));
1253 : :
1254 : : // limit reconstructed energy
1255 : 0 : safeLimit(energyIdx(nmat,k), ul, ur);
1256 : : }
1257 : 0 : }
1258 : :
1259 : : } // tk::
|
