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