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