Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Limiter.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 Limiters for discontiunous Galerkin methods
9 : : \details This file contains functions that provide limiter function
10 : : calculations for maintaining monotonicity near solution discontinuities
11 : : for the DG discretization.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <array>
16 : : #include <vector>
17 : :
18 : : #include "FaceData.hpp"
19 : : #include "Vector.hpp"
20 : : #include "Limiter.hpp"
21 : : #include "DerivedData.hpp"
22 : : #include "Integrate/Quadrature.hpp"
23 : : #include "Integrate/Basis.hpp"
24 : : #include "Inciter/InputDeck/InputDeck.hpp"
25 : : #include "PrefIndicator.hpp"
26 : : #include "Reconstruction.hpp"
27 : : #include "MultiMat/MiscMultiMatFns.hpp"
28 : : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
29 : : #include "MultiSpecies/Mixture/Mixture.hpp"
30 : :
31 : : namespace inciter {
32 : :
33 : : extern ctr::InputDeck g_inputdeck;
34 : :
35 : : void
36 : 0 : WENO_P1( const std::vector< int >& esuel,
37 : : tk::Fields& U )
38 : : // *****************************************************************************
39 : : // Weighted Essentially Non-Oscillatory (WENO) limiter for DGP1
40 : : //! \param[in] esuel Elements surrounding elements
41 : : //! \param[in,out] U High-order solution vector which gets limited
42 : : //! \details This WENO function should be called for transport and compflow
43 : : //! \note This limiter function is experimental and untested. Use with caution.
44 : : // *****************************************************************************
45 : : {
46 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
47 : 0 : const auto cweight = inciter::g_inputdeck.get< tag::cweight >();
48 : 0 : auto nelem = esuel.size()/4;
49 : : std::array< std::vector< tk::real >, 3 >
50 : : limU {{ std::vector< tk::real >(nelem),
51 : : std::vector< tk::real >(nelem),
52 [ - - ][ - - ]: 0 : std::vector< tk::real >(nelem) }};
[ - - ]
53 : :
54 : 0 : std::size_t ncomp = U.nprop()/rdof;
55 : :
56 [ - - ]: 0 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
57 : : {
58 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
59 : : {
60 [ - - ]: 0 : WENOLimiting(U, esuel, e, c, rdof, cweight, limU);
61 : : }
62 : :
63 : 0 : auto mark = c*rdof;
64 : :
65 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
66 : : {
67 [ - - ]: 0 : U(e, mark+1) = limU[0][e];
68 [ - - ]: 0 : U(e, mark+2) = limU[1][e];
69 [ - - ]: 0 : U(e, mark+3) = limU[2][e];
70 : : }
71 : : }
72 : 0 : }
73 : :
74 : : void
75 : 1650 : Superbee_P1( const std::vector< int >& esuel,
76 : : const std::vector< std::size_t >& inpoel,
77 : : const std::vector< std::size_t >& ndofel,
78 : : const tk::UnsMesh::Coords& coord,
79 : : tk::Fields& U )
80 : : // *****************************************************************************
81 : : // Superbee limiter for DGP1
82 : : //! \param[in] esuel Elements surrounding elements
83 : : //! \param[in] inpoel Element connectivity
84 : : //! \param[in] ndofel Vector of local number of degrees of freedom
85 : : //! \param[in] coord Array of nodal coordinates
86 : : //! \param[in,out] U High-order solution vector which gets limited
87 : : //! \details This Superbee function should be called for transport and compflow
88 : : // *****************************************************************************
89 : : {
90 : 1650 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
91 : 1650 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
92 : 1650 : std::size_t ncomp = U.nprop()/rdof;
93 : :
94 : 1650 : auto beta_lim = 2.0;
95 : :
96 [ + + ]: 766680 : for (std::size_t e=0; e<esuel.size()/4; ++e)
97 : : {
98 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
99 : : // basis functions and solutions by default. This implies that P0P1 is
100 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
101 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
102 : : // element for pDG.
103 : : std::size_t dof_el;
104 [ - + ]: 765030 : if (rdof > ndof)
105 : : {
106 : 0 : dof_el = rdof;
107 : : }
108 : : else
109 : : {
110 : 765030 : dof_el = ndofel[e];
111 : : }
112 : :
113 [ + + ]: 765030 : if (dof_el > 1)
114 : : {
115 : : auto phi = SuperbeeLimiting(U, esuel, inpoel, coord, e, ndof, rdof,
116 [ + - ]: 952638 : dof_el, ncomp, beta_lim);
117 : :
118 : : // apply limiter function
119 [ + + ]: 2857914 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
120 : : {
121 : 2381595 : auto mark = c*rdof;
122 [ + - ][ + - ]: 2381595 : U(e, mark+1) = phi[c] * U(e, mark+1);
123 [ + - ][ + - ]: 2381595 : U(e, mark+2) = phi[c] * U(e, mark+2);
124 [ + - ][ + - ]: 2381595 : U(e, mark+3) = phi[c] * U(e, mark+3);
125 : : }
126 : : }
127 : : }
128 : 1650 : }
129 : :
130 : : void
131 : 0 : SuperbeeMultiMat_P1(
132 : : const std::vector< int >& esuel,
133 : : const std::vector< std::size_t >& inpoel,
134 : : const std::vector< std::size_t >& ndofel,
135 : : const tk::UnsMesh::Coords& coord,
136 : : const std::vector< std::size_t >& solidx,
137 : : tk::Fields& U,
138 : : tk::Fields& P,
139 : : std::size_t nmat )
140 : : // *****************************************************************************
141 : : // Superbee limiter for multi-material DGP1
142 : : //! \param[in] esuel Elements surrounding elements
143 : : //! \param[in] inpoel Element connectivity
144 : : //! \param[in] ndofel Vector of local number of degrees of freedom
145 : : //! \param[in] coord Array of nodal coordinates
146 : : //! \param[in] solidx Solid material index indicator
147 : : //! \param[in,out] U High-order solution vector which gets limited
148 : : //! \param[in,out] P High-order vector of primitives which gets limited
149 : : //! \param[in] nmat Number of materials in this PDE system
150 : : //! \details This Superbee function should be called for multimat
151 : : // *****************************************************************************
152 : : {
153 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
154 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
155 : : const auto intsharp = inciter::g_inputdeck.get< tag::multimat,
156 : 0 : tag::intsharp >();
157 : 0 : std::size_t ncomp = U.nprop()/rdof;
158 : 0 : std::size_t nprim = P.nprop()/rdof;
159 : :
160 : 0 : auto beta_lim = 2.0;
161 : :
162 [ - - ]: 0 : for (std::size_t e=0; e<esuel.size()/4; ++e)
163 : : {
164 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
165 : : // basis functions and solutions by default. This implies that P0P1 is
166 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
167 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
168 : : // element for pDG.
169 : : std::size_t dof_el;
170 [ - - ]: 0 : if (rdof > ndof)
171 : : {
172 : 0 : dof_el = rdof;
173 : : }
174 : : else
175 : : {
176 : 0 : dof_el = ndofel[e];
177 : : }
178 : :
179 [ - - ]: 0 : if (dof_el > 1)
180 : : {
181 : : // limit conserved quantities
182 : : auto phic = SuperbeeLimiting(U, esuel, inpoel, coord, e, ndof, rdof,
183 [ - - ]: 0 : dof_el, ncomp, beta_lim);
184 : : // limit primitive quantities
185 : : auto phip = SuperbeeLimiting(P, esuel, inpoel, coord, e, ndof, rdof,
186 [ - - ]: 0 : dof_el, nprim, beta_lim);
187 : :
188 : 0 : std::vector< tk::real > phic_p2;
189 : 0 : std::vector< std::vector< tk::real > > unk, prim;
190 [ - - ]: 0 : if(ndof > 1)
191 [ - - ]: 0 : BoundPreservingLimiting(nmat, ndof, e, inpoel, coord, U, phic,
192 : : phic_p2);
193 : :
194 : : // limits under which compression is to be performed
195 [ - - ]: 0 : std::vector< std::size_t > matInt(nmat, 0);
196 [ - - ]: 0 : std::vector< tk::real > alAvg(nmat, 0.0);
197 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
198 [ - - ]: 0 : alAvg[k] = U(e, volfracDofIdx(nmat,k,rdof,0));
199 [ - - ]: 0 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
200 [ - - ][ - - ]: 0 : if ((intsharp > 0) && intInd)
201 : : {
202 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
203 : : {
204 [ - - ]: 0 : if (matInt[k])
205 : 0 : phic[volfracIdx(nmat,k)] = 1.0;
206 : 0 : }
207 : : }
208 : : else
209 : : {
210 [ - - ]: 0 : if (!g_inputdeck.get< tag::accuracy_test >())
211 [ - - ]: 0 : consistentMultiMatLimiting_P1(nmat, rdof, e, solidx, U, P, phic,
212 : : phic_p2);
213 : : }
214 : :
215 : : // apply limiter function
216 [ - - ]: 0 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
217 : : {
218 : 0 : auto mark = c*rdof;
219 [ - - ][ - - ]: 0 : U(e, mark+1) = phic[c] * U(e, mark+1);
220 [ - - ][ - - ]: 0 : U(e, mark+2) = phic[c] * U(e, mark+2);
221 [ - - ][ - - ]: 0 : U(e, mark+3) = phic[c] * U(e, mark+3);
222 : : }
223 [ - - ]: 0 : for (inciter::ncomp_t c=0; c<nprim; ++c)
224 : : {
225 : 0 : auto mark = c*rdof;
226 [ - - ][ - - ]: 0 : P(e, mark+1) = phip[c] * P(e, mark+1);
227 [ - - ][ - - ]: 0 : P(e, mark+2) = phip[c] * P(e, mark+2);
228 [ - - ][ - - ]: 0 : P(e, mark+3) = phip[c] * P(e, mark+3);
229 : : }
230 : : }
231 : : }
232 : 0 : }
233 : :
234 : : void
235 : 0 : VertexBasedTransport_P1(
236 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
237 : : const std::vector< std::size_t >& inpoel,
238 : : const std::vector< std::size_t >& ndofel,
239 : : std::size_t nelem,
240 : : const tk::UnsMesh::Coords& coord,
241 : : tk::Fields& U )
242 : : // *****************************************************************************
243 : : // Kuzmin's vertex-based limiter for transport DGP1
244 : : //! \param[in] esup Elements surrounding points
245 : : //! \param[in] inpoel Element connectivity
246 : : //! \param[in] ndofel Vector of local number of degrees of freedom
247 : : //! \param[in] nelem Number of elements
248 : : //! \param[in] coord Array of nodal coordinates
249 : : //! \param[in,out] U High-order solution vector which gets limited
250 : : //! \details This vertex-based limiter function should be called for transport.
251 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
252 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
253 : : //! computational and applied mathematics, 233(12), 3077-3085.
254 : : // *****************************************************************************
255 : : {
256 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
257 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
258 : : const auto intsharp = inciter::g_inputdeck.get< tag::transport,
259 : 0 : tag::intsharp >();
260 : 0 : std::size_t ncomp = U.nprop()/rdof;
261 : :
262 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
263 : : {
264 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
265 : : // basis functions and solutions by default. This implies that P0P1 is
266 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
267 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
268 : : // element for pDG.
269 : : std::size_t dof_el;
270 [ - - ]: 0 : if (rdof > ndof)
271 : : {
272 : 0 : dof_el = rdof;
273 : : }
274 : : else
275 : : {
276 : 0 : dof_el = ndofel[e];
277 : : }
278 : :
279 [ - - ]: 0 : if (dof_el > 1)
280 : : {
281 [ - - ]: 0 : std::vector< tk::real > phi(ncomp, 1.0);
282 : 0 : std::vector< std::size_t > var;
283 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) var.push_back(c);
284 : : // limit conserved quantities
285 [ - - ]: 0 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
286 : : ncomp, phi, var);
287 : :
288 : : // limits under which compression is to be performed
289 [ - - ]: 0 : std::vector< std::size_t > matInt(ncomp, 0);
290 [ - - ]: 0 : std::vector< tk::real > alAvg(ncomp, 0.0);
291 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k)
292 [ - - ]: 0 : alAvg[k] = U(e,k*rdof);
293 [ - - ]: 0 : auto intInd = interfaceIndicator(ncomp, alAvg, matInt);
294 [ - - ][ - - ]: 0 : if ((intsharp > 0) && intInd)
295 : : {
296 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k)
297 : : {
298 [ - - ]: 0 : if (matInt[k]) phi[k] = 1.0;
299 : : }
300 : : }
301 : :
302 : : // apply limiter function
303 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
304 : : {
305 : 0 : auto mark = c*rdof;
306 [ - - ][ - - ]: 0 : U(e, mark+1) = phi[c] * U(e, mark+1);
307 [ - - ][ - - ]: 0 : U(e, mark+2) = phi[c] * U(e, mark+2);
308 [ - - ][ - - ]: 0 : U(e, mark+3) = phi[c] * U(e, mark+3);
309 : : }
310 : : }
311 : : }
312 : 0 : }
313 : :
314 : : void
315 : 0 : VertexBasedCompflow_P1(
316 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
317 : : const std::vector< std::size_t >& inpoel,
318 : : const std::vector< std::size_t >& ndofel,
319 : : std::size_t nelem,
320 : : const std::vector< inciter::EOS >& mat_blk,
321 : : const inciter::FaceData& fd,
322 : : const tk::Fields& geoFace,
323 : : const tk::Fields& geoElem,
324 : : const tk::UnsMesh::Coords& coord,
325 : : const tk::FluxFn& flux,
326 : : const std::vector< std::size_t >& solidx,
327 : : tk::Fields& U,
328 : : std::vector< std::size_t >& shockmarker )
329 : : // *****************************************************************************
330 : : // Kuzmin's vertex-based limiter for single-material DGP1
331 : : //! \param[in] esup Elements surrounding points
332 : : //! \param[in] inpoel Element connectivity
333 : : //! \param[in] ndofel Vector of local number of degrees of freedom
334 : : //! \param[in] nelem Number of elements
335 : : //! \param[in] mat_blk EOS material block
336 : : //! \param[in] fd Face connectivity and boundary conditions object
337 : : //! \param[in] geoFace Face geometry array
338 : : // //! \param[in] geoElem Element geometry array
339 : : //! \param[in] coord Array of nodal coordinates
340 : : //! \param[in] flux Riemann flux function to use
341 : : //! \param[in] solidx Solid material index indicator
342 : : //! \param[in,out] U High-order solution vector which gets limited
343 : : //! \param[in,out] shockmarker Shock detection marker array
344 : : //! \details This vertex-based limiter function should be called for compflow.
345 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
346 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
347 : : //! computational and applied mathematics, 233(12), 3077-3085.
348 : : // *****************************************************************************
349 : : {
350 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
351 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
352 : 0 : std::size_t ncomp = U.nprop()/rdof;
353 : :
354 : : // Null field for MarkShockCells argument
355 : 0 : tk::Fields P;
356 : :
357 [ - - ]: 0 : if (inciter::g_inputdeck.get< tag::shock_detector_coeff >() > 1e-6) {
358 : : // Indicator based on momentum flux jump
359 : 0 : std::set< std::size_t > vars;
360 [ - - ][ - - ]: 0 : for (std::size_t i=1; i<=3; ++i) vars.insert(i);
361 [ - - ]: 0 : MarkShockCells(false, nelem, 1, ndof, rdof, mat_blk, ndofel,
362 : : inpoel, coord, fd, geoFace, geoElem, flux, solidx, U, P,
363 : : vars, shockmarker);
364 : : }
365 : :
366 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
367 : : {
368 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
369 : : // basis functions and solutions by default. This implies that P0P1 is
370 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
371 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
372 : : // element for pDG.
373 : : std::size_t dof_el;
374 [ - - ]: 0 : if (rdof > ndof)
375 : : {
376 : 0 : dof_el = rdof;
377 : : }
378 : : else
379 : : {
380 : 0 : dof_el = ndofel[e];
381 : : }
382 : :
383 [ - - ][ - - ]: 0 : if (dof_el > 1 && shockmarker[e])
[ - - ]
384 : : {
385 [ - - ]: 0 : std::vector< tk::real > phi(ncomp, 1.0);
386 : 0 : std::vector< std::size_t > var;
387 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) var.push_back(c);
388 : : // limit conserved quantities
389 [ - - ]: 0 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
390 : : ncomp, phi, var);
391 : :
392 : : // apply limiter function
393 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
394 : : {
395 : 0 : auto mark = c*rdof;
396 [ - - ][ - - ]: 0 : U(e, mark+1) = phi[c] * U(e, mark+1);
397 [ - - ][ - - ]: 0 : U(e, mark+2) = phi[c] * U(e, mark+2);
398 [ - - ][ - - ]: 0 : U(e, mark+3) = phi[c] * U(e, mark+3);
399 : : }
400 : : }
401 : : }
402 : 0 : }
403 : :
404 : : void
405 : 2580 : VertexBasedCompflow_P2(
406 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
407 : : const std::vector< std::size_t >& inpoel,
408 : : const std::vector< std::size_t >& ndofel,
409 : : std::size_t nelem,
410 : : const std::vector< inciter::EOS >& mat_blk,
411 : : const inciter::FaceData& fd,
412 : : const tk::Fields& geoFace,
413 : : const tk::Fields& geoElem,
414 : : const tk::UnsMesh::Coords& coord,
415 : : [[maybe_unused]] const std::vector< std::size_t >& gid,
416 : : [[maybe_unused]] const std::unordered_map< std::size_t, std::size_t >& bid,
417 : : [[maybe_unused]] const std::vector< std::vector<tk::real> >& uNodalExtrm,
418 : : [[maybe_unused]] const std::vector< std::vector<tk::real> >& mtInv,
419 : : const tk::FluxFn& flux,
420 : : const std::vector< std::size_t >& solidx,
421 : : tk::Fields& U,
422 : : std::vector< std::size_t >& shockmarker )
423 : : // *****************************************************************************
424 : : // Kuzmin's vertex-based limiter on reference element for single-material DGP2
425 : : //! \param[in] esup Elements surrounding points
426 : : //! \param[in] inpoel Element connectivity
427 : : //! \param[in] ndofel Vector of local number of degrees of freedom
428 : : //! \param[in] nelem Number of elements
429 : : //! \param[in] mat_blk EOS material block
430 : : //! \param[in] fd Face connectivity and boundary conditions object
431 : : //! \param[in] geoFace Face geometry array
432 : : // //! \param[in] geoElem Element geometry array
433 : : //! \param[in] coord Array of nodal coordinates
434 : : //! \param[in] gid Local->global node id map
435 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
436 : : //! global node ids (key)
437 : : //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
438 : : //! variables
439 : : //! \param[in] mtInv Inverse of Taylor mass matrix
440 : : //! \param[in] flux Riemann flux function to use
441 : : //! \param[in] solidx Solid material index indicator
442 : : //! \param[in,out] U High-order solution vector which gets limited
443 : : //! \param[in,out] shockmarker Shock detection marker array
444 : : //! \details This vertex-based limiter function should be called for compflow.
445 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
446 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
447 : : //! computational and applied mathematics, 233(12), 3077-3085.
448 : : // *****************************************************************************
449 : : {
450 : 2580 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
451 : 2580 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
452 : 2580 : std::size_t ncomp = U.nprop()/rdof;
453 : :
454 : : // Null field for MarkShockCells argument
455 : 5160 : tk::Fields P;
456 : :
457 [ + - ]: 2580 : if (inciter::g_inputdeck.get< tag::shock_detector_coeff >() > 1e-6) {
458 : : // Indicator based on momentum flux jump
459 : 5160 : std::set< std::size_t > vars;
460 [ + + ][ + - ]: 10320 : for (std::size_t i=1; i<=3; ++i) vars.insert(i);
461 [ + - ]: 2580 : MarkShockCells(false, nelem, 1, ndof, rdof, mat_blk, ndofel,
462 : : inpoel, coord, fd, geoFace, geoElem, flux, solidx, U, P,
463 : : vars, shockmarker);
464 : : }
465 : :
466 [ + + ]: 366420 : for (std::size_t e=0; e<nelem; ++e)
467 : : {
468 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
469 : : // basis functions and solutions by default. This implies that P0P1 is
470 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
471 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
472 : : // element for pDG.
473 : : std::size_t dof_el;
474 [ - + ]: 363840 : if (rdof > ndof)
475 : : {
476 : 0 : dof_el = rdof;
477 : : }
478 : : else
479 : : {
480 : 363840 : dof_el = ndofel[e];
481 : : }
482 : :
483 [ + + ][ + + ]: 363840 : if (dof_el > 1 && shockmarker[e])
[ + + ]
484 : : {
485 : : // The vector of limiting coefficients for P1 and P2 coefficients
486 [ + - ]: 33128 : std::vector< tk::real > phic_p1(ncomp, 1.0);
487 : :
488 : : // Removing 3rd order DOFs if discontinuity is detected, and applying
489 : : // limiting to the 2nd order/P1 solution
490 [ + + ]: 99384 : for (std::size_t c=0; c<ncomp; ++c) {
491 [ + + ]: 579740 : for(std::size_t idof = 4; idof < rdof; idof++) {
492 : 496920 : auto mark = c * rdof + idof;
493 [ + - ]: 496920 : U(e, mark) = 0.0;
494 : : }
495 : : }
496 : :
497 : : // Obtain limiting coefficient for P1 coefficients
498 : 33128 : std::vector< std::size_t > var;
499 [ + + ][ + - ]: 99384 : for (std::size_t c=0; c<ncomp; ++c) var.push_back(c);
500 [ + - ]: 16564 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
501 : : ncomp, phic_p1, var);
502 : :
503 : : // apply limiter function to the solution with Taylor basis
504 [ + + ]: 99384 : for (std::size_t c=0; c<ncomp; ++c) {
505 : 82820 : auto mark = c * rdof;
506 [ + + ]: 331280 : for(std::size_t idof=1; idof<4; idof++)
507 [ + - ][ + - ]: 248460 : U(e, mark+idof) = phic_p1[c] * U(e, mark+idof);
508 : : }
509 : : }
510 : : }
511 : 2580 : }
512 : :
513 : : void
514 : 4080 : VertexBasedMultiMat_P1(
515 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
516 : : const std::vector< std::size_t >& inpoel,
517 : : const std::vector< std::size_t >& ndofel,
518 : : std::size_t nelem,
519 : : const std::vector< inciter::EOS >& mat_blk,
520 : : const inciter::FaceData& fd,
521 : : const tk::Fields& geoFace,
522 : : const tk::Fields& geoElem,
523 : : const tk::UnsMesh::Coords& coord,
524 : : const tk::FluxFn& flux,
525 : : const std::vector< std::size_t >& solidx,
526 : : tk::Fields& U,
527 : : tk::Fields& P,
528 : : std::size_t nmat,
529 : : std::vector< std::size_t >& shockmarker )
530 : : // *****************************************************************************
531 : : // Kuzmin's vertex-based limiter for multi-material DGP1
532 : : //! \param[in] esup Elements surrounding points
533 : : //! \param[in] inpoel Element connectivity
534 : : //! \param[in] ndofel Vector of local number of degrees of freedom
535 : : //! \param[in] nelem Number of elements
536 : : //! \param[in] mat_blk EOS material block
537 : : //! \param[in] fd Face connectivity and boundary conditions object
538 : : //! \param[in] geoFace Face geometry array
539 : : // //! \param[in] geoElem Element geometry array
540 : : //! \param[in] coord Array of nodal coordinates
541 : : //! \param[in] flux Riemann flux function to use
542 : : //! \param[in] solidx Solid material index indicator
543 : : //! \param[in,out] U High-order solution vector which gets limited
544 : : //! \param[in,out] P High-order vector of primitives which gets limited
545 : : //! \param[in] nmat Number of materials in this PDE system
546 : : //! \param[in,out] shockmarker Shock detection marker array
547 : : //! \details This vertex-based limiter function should be called for multimat.
548 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
549 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
550 : : //! computational and applied mathematics, 233(12), 3077-3085.
551 : : // *****************************************************************************
552 : : {
553 : 4080 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
554 : 4080 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
555 : : const auto intsharp = inciter::g_inputdeck.get< tag::multimat,
556 : 4080 : tag::intsharp >();
557 : 4080 : std::size_t ncomp = U.nprop()/rdof;
558 : 4080 : std::size_t nprim = P.nprop()/rdof;
559 : :
560 : : // Evaluate the interface condition and mark the shock cells
561 : 4080 : if (inciter::g_inputdeck.get< tag::shock_detector_coeff >()
562 [ + + ][ + - ]: 4080 : > 1e-6 && ndof > 1) {
[ + + ]
563 : : // Indicator based on momentum flux jump
564 : 1500 : std::set< std::size_t > vars;
565 [ + + ][ + - ]: 3000 : for (std::size_t i=0; i<3; ++i) vars.insert(momentumIdx(nmat, i));
566 [ + - ]: 750 : MarkShockCells(false, nelem, nmat, ndof, rdof, mat_blk, ndofel,
567 : : inpoel, coord, fd, geoFace, geoElem, flux, solidx, U, P,
568 : : vars, shockmarker);
569 : : }
570 : :
571 [ + + ]: 845460 : for (std::size_t e=0; e<nelem; ++e)
572 : : {
573 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
574 : : // basis functions and solutions by default. This implies that P0P1 is
575 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
576 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
577 : : // element for pDG.
578 : : std::size_t dof_el;
579 [ + + ]: 841380 : if (rdof > ndof)
580 : : {
581 : 341100 : dof_el = rdof;
582 : : }
583 : : else
584 : : {
585 : 500280 : dof_el = ndofel[e];
586 : : }
587 : :
588 [ + - ]: 841380 : if (dof_el > 1)
589 : : {
590 [ + - ]: 1682760 : std::vector< tk::real > phic(ncomp, 1.0);
591 [ + - ]: 1682760 : std::vector< tk::real > phip(nprim, 1.0);
592 [ + + ]: 841380 : if(shockmarker[e]) {
593 : : // When shockmarker is 1, there is discontinuity within the element.
594 : : // Hence, the vertex-based limiter will be applied.
595 : :
596 : : // limit conserved quantities
597 : 840640 : std::vector< std::size_t > varc;
598 [ + + ][ + - ]: 4203200 : for (std::size_t c=0; c<ncomp; ++c) varc.push_back(c);
599 [ + - ]: 420320 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
600 : : ncomp, phic, varc);
601 : : // limit primitive quantities
602 : 840640 : std::vector< std::size_t > varp;
603 [ + + ][ + - ]: 2521920 : for (std::size_t c=0; c<nprim; ++c) varp.push_back(c);
604 [ + - ]: 420320 : VertexBasedLimiting(P, esup, inpoel, coord, e, rdof, dof_el,
605 : : nprim, phip, varp);
606 : : } else {
607 : : // When shockmarker is 0, the volume fraction, density and energy
608 : : // of minor material will still be limited to ensure a stable solution.
609 : 842120 : std::vector< std::size_t > vars;
610 [ + + ][ + - ]: 1263180 : for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat,k));
611 [ + - ]: 421060 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
612 : : ncomp, phic, vars);
613 : :
614 [ + + ]: 1263180 : for(std::size_t k=0; k<nmat; ++k) {
615 [ + - ][ + + ]: 842120 : if(U(e, volfracDofIdx(nmat,k,rdof,0)) < 1e-4) {
616 : : // limit the density and energy of minor materials
617 : 421060 : vars.clear();
618 [ + - ]: 421060 : vars.push_back(densityIdx(nmat, k));
619 [ + - ]: 421060 : vars.push_back(energyIdx(nmat, k));
620 [ - + ]: 421060 : if (solidx[k] > 0) {
621 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
622 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
623 [ - - ]: 0 : vars.push_back(deformIdx(nmat, solidx[k], i, j));
624 : : }
625 [ + - ]: 421060 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
626 : : ncomp, phic, vars);
627 : :
628 : : // limit the pressure of minor materials
629 [ + - ]: 421060 : VertexBasedLimiting(P, esup, inpoel, coord, e, rdof, dof_el,
630 [ + - ]: 842120 : nprim, phip, std::vector< std::size_t >{pressureIdx(nmat, k)});
631 : : }
632 : : }
633 : : }
634 : :
635 : 1682760 : std::vector< tk::real > phic_p2, phip_p2;
636 : :
637 [ + - ]: 841380 : PositivityLimiting(nmat, 1, mat_blk, rdof, dof_el, ndofel, e,
638 : : inpoel, coord, fd.Esuel(), U, P, phic, phic_p2, phip, phip_p2);
639 : :
640 : : // limits under which compression is to be performed
641 [ + - ]: 1682760 : std::vector< std::size_t > matInt(nmat, 0);
642 [ + - ]: 1682760 : std::vector< tk::real > alAvg(nmat, 0.0);
643 [ + + ]: 2524140 : for (std::size_t k=0; k<nmat; ++k)
644 [ + - ]: 1682760 : alAvg[k] = U(e, volfracDofIdx(nmat,k,rdof,0));
645 [ + - ]: 841380 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
646 [ + + ][ + + ]: 841380 : if ((intsharp > 0) && intInd) {
647 [ + + ]: 32214 : for (std::size_t k=0; k<nmat; ++k) {
648 [ + - ]: 21476 : if (matInt[k]) {
649 : 21476 : phic[volfracIdx(nmat,k)] = 1.0;
650 : : }
651 : 10738 : }
652 : : }
653 : : else {
654 [ + - ]: 830642 : if(nmat > 1)
655 [ + - ]: 830642 : BoundPreservingLimiting(nmat, ndof, e, inpoel, coord, U, phic,
656 : : phic_p2);
657 : :
658 [ + - ]: 830642 : if (!g_inputdeck.get< tag::accuracy_test >())
659 [ + - ]: 830642 : consistentMultiMatLimiting_P1(nmat, rdof, e, solidx, U, P, phic,
660 : : phic_p2);
661 : : }
662 : :
663 : : // apply limiter function
664 [ + + ]: 8413800 : for (std::size_t c=0; c<ncomp; ++c)
665 : : {
666 : 7572420 : auto mark = c*rdof;
667 [ + - ][ + - ]: 7572420 : U(e, mark+1) = phic[c] * U(e, mark+1);
668 [ + - ][ + - ]: 7572420 : U(e, mark+2) = phic[c] * U(e, mark+2);
669 [ + - ][ + - ]: 7572420 : U(e, mark+3) = phic[c] * U(e, mark+3);
670 : : }
671 [ + + ]: 5048280 : for (std::size_t c=0; c<nprim; ++c)
672 : : {
673 : 4206900 : auto mark = c*rdof;
674 [ + - ][ + - ]: 4206900 : P(e, mark+1) = phip[c] * P(e, mark+1);
675 [ + - ][ + - ]: 4206900 : P(e, mark+2) = phip[c] * P(e, mark+2);
676 [ + - ][ + - ]: 4206900 : P(e, mark+3) = phip[c] * P(e, mark+3);
677 : : }
678 : : }
679 : : }
680 : 4080 : }
681 : :
682 : : void
683 : 0 : VertexBasedMultiMat_P2(
684 : : const bool pref,
685 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
686 : : const std::vector< std::size_t >& inpoel,
687 : : const std::vector< std::size_t >& ndofel,
688 : : std::size_t nelem,
689 : : const std::vector< inciter::EOS >& mat_blk,
690 : : const inciter::FaceData& fd,
691 : : const tk::Fields& geoFace,
692 : : const tk::Fields& geoElem,
693 : : const tk::UnsMesh::Coords& coord,
694 : : [[maybe_unused]] const std::vector< std::size_t >& gid,
695 : : [[maybe_unused]] const std::unordered_map< std::size_t, std::size_t >& bid,
696 : : [[maybe_unused]] const std::vector< std::vector<tk::real> >& uNodalExtrm,
697 : : [[maybe_unused]] const std::vector< std::vector<tk::real> >& pNodalExtrm,
698 : : [[maybe_unused]] const std::vector< std::vector<tk::real> >& mtInv,
699 : : const tk::FluxFn& flux,
700 : : const std::vector< std::size_t >& solidx,
701 : : tk::Fields& U,
702 : : tk::Fields& P,
703 : : std::size_t nmat,
704 : : std::vector< std::size_t >& shockmarker )
705 : : // *****************************************************************************
706 : : // Kuzmin's vertex-based limiter for multi-material DGP2
707 : : //! \param[in] pref Indicator for p-adaptive algorithm
708 : : //! \param[in] esup Elements surrounding points
709 : : //! \param[in] inpoel Element connectivity
710 : : //! \param[in] ndofel Vector of local number of degrees of freedom
711 : : //! \param[in] nelem Number of elements
712 : : //! \param[in] mat_blk EOS material block
713 : : //! \param[in] fd Face connectivity and boundary conditions object
714 : : //! \param[in] geoFace Face geometry array
715 : : // //! \param[in] geoElem Element geometry array
716 : : //! \param[in] coord Array of nodal coordinates
717 : : //! \param[in] gid Local->global node id map
718 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
719 : : //! global node ids (key)
720 : : //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
721 : : //! variables
722 : : //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
723 : : //! variables
724 : : //! \param[in] mtInv Inverse of Taylor mass matrix
725 : : //! \param[in] flux Riemann flux function to use
726 : : //! \param[in] solidx Solid material index indicator
727 : : //! \param[in,out] U High-order solution vector which gets limited
728 : : //! \param[in,out] P High-order vector of primitives which gets limited
729 : : //! \param[in] nmat Number of materials in this PDE system
730 : : //! \param[in,out] shockmarker Shock detection marker array
731 : : //! \details This vertex-based limiter function should be called for multimat.
732 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
733 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
734 : : //! computational and applied mathematics, 233(12), 3077-3085.
735 : : // *****************************************************************************
736 : : {
737 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
738 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
739 : : const auto intsharp = inciter::g_inputdeck.get< tag::multimat,
740 : 0 : tag::intsharp >();
741 : 0 : std::size_t ncomp = U.nprop()/rdof;
742 : 0 : std::size_t nprim = P.nprop()/rdof;
743 : :
744 : : // Evaluate the interface condition and mark the shock cells
745 : 0 : if (inciter::g_inputdeck.get< tag::shock_detector_coeff >()
746 [ - - ]: 0 : > 1e-6) {
747 : : // Indicator based on momentum flux jump
748 : 0 : std::set< std::size_t > vars;
749 [ - - ][ - - ]: 0 : for (std::size_t i=0; i<3; ++i) vars.insert(momentumIdx(nmat, i));
750 [ - - ]: 0 : MarkShockCells(pref, nelem, nmat, ndof, rdof, mat_blk, ndofel,
751 : : inpoel, coord, fd, geoFace, geoElem, flux, solidx, U, P,
752 : : vars, shockmarker);
753 : : }
754 : :
755 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
756 : : {
757 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
758 : : // basis functions and solutions by default. This implies that P0P1 is
759 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
760 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
761 : : // element for pDG.
762 : : std::size_t dof_el;
763 [ - - ]: 0 : if (rdof > ndof)
764 : : {
765 : 0 : dof_el = rdof;
766 : : }
767 : : else
768 : : {
769 : 0 : dof_el = ndofel[e];
770 : : }
771 : :
772 : : // For multi-material simulation, when dofel = 1, p0p1 is applied and ndof
773 : : // for solution evaluation should be 4
774 [ - - ][ - - ]: 0 : if(ncomp > 5 && dof_el == 1)
775 : 0 : dof_el = 4;
776 : :
777 [ - - ]: 0 : if (dof_el > 1)
778 : : {
779 : : // The vector of limiting coefficients for P1
780 [ - - ][ - - ]: 0 : std::vector< tk::real > phic_p1(ncomp, 1.0), phic_p2(ncomp, 1.0);
781 [ - - ][ - - ]: 0 : std::vector< tk::real > phip_p1(nprim, 1.0), phip_p2(nprim, 1.0);
782 : :
783 : : // Only when the cell is marked with discontinuous solution or P0P1 scheme
784 : : // is used, the vertex-based slope limiter will be applied.
785 [ - - ][ - - ]: 0 : if(shockmarker[e] || dof_el == 4) {
[ - - ]
786 : : // Removing 3rd order DOFs if discontinuity is detected, and applying
787 : : // limiting to the 2nd order/P1 solution
788 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
789 : 0 : auto mark = c * rdof;
790 [ - - ]: 0 : for(std::size_t idof = 4; idof < rdof; idof++)
791 [ - - ]: 0 : U(e, mark+idof) = 0.0;
792 : : }
793 [ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c) {
794 : 0 : auto mark = c * rdof;
795 [ - - ]: 0 : for(std::size_t idof = 4; idof < rdof; idof++)
796 [ - - ]: 0 : P(e, mark+idof) = 0.0;
797 : : }
798 : :
799 : : // Obtain limiter coefficient for P1 conserved quantities
800 : 0 : std::vector< std::size_t > varc;
801 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) varc.push_back(c);
802 [ - - ]: 0 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
803 : : ncomp, phic_p1, varc);
804 : : // Obtain limiter coefficient for P1 primitive quantities
805 : 0 : std::vector< std::size_t > varp;
806 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c) varp.push_back(c);
807 [ - - ]: 0 : VertexBasedLimiting(P, esup, inpoel, coord, e, rdof, dof_el,
808 : : nprim, phip_p1, varp);
809 : : } else {
810 : : // When shockmarker is 0, the volume fraction will still be limited to
811 : : // ensure a stable solution. Since the limiting strategy for third order
812 : : // solution will downgrade the accuracy to second order, the density,
813 : : // energy and pressure of minor material will not be limited.
814 : 0 : std::vector< std::size_t > vars;
815 [ - - ][ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat,k));
816 [ - - ]: 0 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
817 : : ncomp, phic_p1, vars);
818 : :
819 : : //for(std::size_t k=0; k<nmat; ++k) {
820 : : // if(U(e, volfracDofIdx(nmat,k,rdof,0)) < 1e-4) {
821 : : // // limit the density of minor materials
822 : : // VertexBasedLimiting(unk, U, esup, inpoel, coord, e, rdof, dof_el,
823 : : // ncomp, phic_p1, std::vector< std::size_t >{densityIdx(nmat,k)});
824 : :
825 : : // // limit the pressure of minor materials
826 : : // VertexBasedLimiting(prim, P, esup, inpoel, coord, e, rdof, dof_el,
827 : : // nprim, phip_p1, std::vector< std::size_t >{pressureIdx(nmat,k)});
828 : : // }
829 : : //}
830 : : }
831 : :
832 [ - - ]: 0 : PositivityLimiting(nmat, 1, mat_blk, ndof, dof_el, ndofel, e,
833 : : inpoel, coord, fd.Esuel(), U, P, phic_p1, phic_p2, phip_p1, phic_p2);
834 : :
835 : : // limits under which compression is to be performed
836 [ - - ]: 0 : std::vector< std::size_t > matInt(nmat, 0);
837 [ - - ]: 0 : std::vector< tk::real > alAvg(nmat, 0.0);
838 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
839 [ - - ]: 0 : alAvg[k] = U(e, volfracDofIdx(nmat,k,rdof,0));
840 [ - - ]: 0 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
841 [ - - ][ - - ]: 0 : if ((intsharp > 0) && intInd) {
842 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
843 [ - - ]: 0 : if (matInt[k]) {
844 : 0 : phic_p1[volfracIdx(nmat,k)] = 1.0;
845 : : }
846 : 0 : }
847 : : }
848 : : else {
849 [ - - ]: 0 : if(nmat > 1)
850 [ - - ]: 0 : BoundPreservingLimiting(nmat, ndof, e, inpoel, coord, U,
851 : : phic_p1, phic_p2);
852 : :
853 [ - - ]: 0 : if (!g_inputdeck.get< tag::accuracy_test >())
854 [ - - ]: 0 : consistentMultiMatLimiting_P1(nmat, rdof, e, solidx, U, P, phic_p1,
855 : : phic_p2);
856 : : }
857 : :
858 : : // apply limiing coefficient
859 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
860 : : {
861 : 0 : auto mark = c * rdof;
862 [ - - ]: 0 : for(std::size_t idof=1; idof<4; idof++)
863 [ - - ][ - - ]: 0 : U(e, mark+idof) = phic_p1[c] * U(e, mark+idof);
864 [ - - ]: 0 : for(std::size_t idof=4; idof<rdof; idof++)
865 [ - - ][ - - ]: 0 : U(e, mark+idof) = phic_p2[c] * U(e, mark+idof);
866 : : }
867 [ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c)
868 : : {
869 : 0 : auto mark = c * rdof;
870 [ - - ]: 0 : for(std::size_t idof=1; idof<4; idof++)
871 [ - - ][ - - ]: 0 : P(e, mark+idof) = phip_p1[c] * P(e, mark+idof);
872 [ - - ]: 0 : for(std::size_t idof=4; idof<rdof; idof++)
873 [ - - ][ - - ]: 0 : P(e, mark+idof) = phip_p2[c] * P(e, mark+idof);
874 : : }
875 : : }
876 : : }
877 : 0 : }
878 : :
879 : : void
880 : 1368 : VertexBasedMultiMat_FV(
881 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
882 : : const std::vector< std::size_t >& inpoel,
883 : : std::size_t nelem,
884 : : const tk::UnsMesh::Coords& coord,
885 : : const std::vector< int >& srcFlag,
886 : : const std::vector< std::size_t >& solidx,
887 : : tk::Fields& U,
888 : : tk::Fields& P,
889 : : std::size_t nmat )
890 : : // *****************************************************************************
891 : : // Kuzmin's vertex-based limiter for multi-material FV
892 : : //! \param[in] esup Elements surrounding points
893 : : //! \param[in] inpoel Element connectivity
894 : : //! \param[in] nelem Number of elements
895 : : //! \param[in] coord Array of nodal coordinates
896 : : //! \param[in] srcFlag Whether the energy source was added
897 : : //! \param[in] solidx Solid material index indicator
898 : : //! \param[in,out] U High-order solution vector which gets limited
899 : : //! \param[in,out] P High-order vector of primitives which gets limited
900 : : //! \param[in] nmat Number of materials in this PDE system
901 : : //! \details This vertex-based limiter function should be called for multimat.
902 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
903 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
904 : : //! computational and applied mathematics, 233(12), 3077-3085.
905 : : // *****************************************************************************
906 : : {
907 : 1368 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
908 : : const auto intsharp = inciter::g_inputdeck.get< tag::multimat,
909 : 1368 : tag::intsharp >();
910 : 1368 : std::size_t ncomp = U.nprop()/rdof;
911 : 1368 : std::size_t nprim = P.nprop()/rdof;
912 : :
913 [ + + ]: 266168 : for (std::size_t e=0; e<nelem; ++e)
914 : : {
915 [ + - ]: 529600 : std::vector< tk::real > phic(ncomp, 1.0);
916 [ + - ]: 529600 : std::vector< tk::real > phip(nprim, 1.0);
917 : : // limit conserved quantities
918 : 529600 : std::vector< std::size_t > var;
919 [ + + ]: 846960 : for (std::size_t k=0; k<nmat; ++k) {
920 [ + - ]: 582160 : var.push_back(volfracIdx(nmat,k));
921 [ + - ]: 582160 : var.push_back(densityIdx(nmat,k));
922 : : }
923 [ + - ]: 264800 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, rdof, ncomp,
924 : : phic, var);
925 : : // limit primitive quantities
926 : 264800 : var.clear();
927 [ + + ][ + - ]: 1641360 : for (std::size_t c=0; c<nprim; ++c) var.push_back(c);
928 [ + - ]: 264800 : VertexBasedLimiting(P, esup, inpoel, coord, e, rdof, rdof, nprim,
929 : : phip, var);
930 : :
931 : : // limits under which compression is to be performed
932 [ + - ]: 529600 : std::vector< std::size_t > matInt(nmat, 0);
933 [ + - ]: 529600 : std::vector< tk::real > alAvg(nmat, 0.0);
934 [ + + ]: 846960 : for (std::size_t k=0; k<nmat; ++k)
935 [ + - ]: 582160 : alAvg[k] = U(e, volfracDofIdx(nmat,k,rdof,0));
936 [ + - ]: 264800 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
937 [ + + ][ + + ]: 264800 : if ((intsharp > 0) && intInd && srcFlag[e] == 0)
[ + - ][ + + ]
938 : : {
939 [ + + ]: 7791 : for (std::size_t k=0; k<nmat; ++k)
940 : : {
941 [ + - ]: 5194 : if (matInt[k])
942 : 5194 : phic[volfracIdx(nmat,k)] = 1.0;
943 : : }
944 : : }
945 : : else
946 : : {
947 [ + - ]: 262203 : if (!g_inputdeck.get< tag::accuracy_test >()) {
948 [ + - ]: 524406 : std::vector< tk::real > phic_p2(ncomp, 1.0);
949 [ + - ]: 262203 : consistentMultiMatLimiting_P1(nmat, rdof, e, solidx, U, P, phic,
950 : : phic_p2);
951 : : }
952 : : }
953 : :
954 : : // apply limiter function
955 [ + + ]: 2805680 : for (std::size_t c=0; c<ncomp; ++c)
956 : : {
957 : 2540880 : auto mark = c*rdof;
958 [ + - ][ + - ]: 2540880 : U(e, mark+1) = phic[c] * U(e, mark+1);
959 [ + - ][ + - ]: 2540880 : U(e, mark+2) = phic[c] * U(e, mark+2);
960 [ + - ][ + - ]: 2540880 : U(e, mark+3) = phic[c] * U(e, mark+3);
961 : : }
962 [ + + ]: 1641360 : for (std::size_t c=0; c<nprim; ++c)
963 : : {
964 : 1376560 : auto mark = c*rdof;
965 [ + - ][ + - ]: 1376560 : P(e, mark+1) = phip[c] * P(e, mark+1);
966 [ + - ][ + - ]: 1376560 : P(e, mark+2) = phip[c] * P(e, mark+2);
967 [ + - ][ + - ]: 1376560 : P(e, mark+3) = phip[c] * P(e, mark+3);
968 : : }
969 : : }
970 : 1368 : }
971 : :
972 : : void
973 : 6600 : VertexBasedMultiSpecies_P1(
974 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
975 : : const std::vector< std::size_t >& inpoel,
976 : : const std::vector< std::size_t >& ndofel,
977 : : std::size_t nelem,
978 : : const std::vector< inciter::EOS >& mat_blk,
979 : : const inciter::FaceData& fd,
980 : : const tk::Fields& geoFace,
981 : : const tk::Fields& geoElem,
982 : : const tk::UnsMesh::Coords& coord,
983 : : const tk::FluxFn& flux,
984 : : const std::vector< std::size_t >& solidx,
985 : : tk::Fields& U,
986 : : tk::Fields& P,
987 : : std::size_t nspec,
988 : : std::vector< std::size_t >& shockmarker )
989 : : // *****************************************************************************
990 : : // Kuzmin's vertex-based limiter for multi-species DGP1
991 : : //! \param[in] esup Elements surrounding points
992 : : //! \param[in] inpoel Element connectivity
993 : : //! \param[in] ndofel Vector of local number of degrees of freedom
994 : : //! \param[in] nelem Number of elements
995 : : //! \param[in] mat_blk EOS material block
996 : : //! \param[in] fd Face connectivity and boundary conditions object
997 : : //! \param[in] geoFace Face geometry array
998 : : //! \param[in] geoElem Element geometry array
999 : : //! \param[in] coord Array of nodal coordinates
1000 : : //! \param[in] flux Riemann flux function to use
1001 : : //! \param[in] solidx Solid material index indicator
1002 : : //! \param[in,out] U High-order solution vector which gets limited
1003 : : //! \param[in,out] P High-order primitive vector which gets limited
1004 : : //! \param[in] nspec Number of species in this PDE system
1005 : : //! \param[in,out] shockmarker Shock detection marker array
1006 : : //! \details This vertex-based limiter function should be called for
1007 : : //! multispecies. For details see: Kuzmin, D. (2010). A vertex-based
1008 : : //! hierarchical slope limiter for p-adaptive discontinuous Galerkin methods.
1009 : : //! Journal of computational and applied mathematics, 233(12), 3077-3085.
1010 : : // *****************************************************************************
1011 : : {
1012 : 6600 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
1013 : 6600 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
1014 : 6600 : std::size_t ncomp = U.nprop()/rdof;
1015 : 6600 : std::size_t nprim = P.nprop()/rdof;
1016 : :
1017 : : // Evaluate the interface condition and mark the shock cells
1018 : 6600 : if (inciter::g_inputdeck.get< tag::shock_detector_coeff >()
1019 [ + - ][ - + ]: 6600 : > 1e-6 && ndof > 1) {
[ - + ]
1020 : : // Indicator based on momentum flux jump
1021 : 0 : std::set< std::size_t > vars;
1022 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1023 [ - - ]: 0 : vars.insert(multispecies::momentumIdx(nspec, i));
1024 [ - - ]: 0 : MarkShockCells(false, nelem, 1, ndof, rdof, mat_blk,
1025 : : ndofel, inpoel, coord, fd, geoFace, geoElem, flux, solidx, U, P,
1026 : : vars, shockmarker);
1027 : : }
1028 : :
1029 [ + + ]: 688800 : for (std::size_t e=0; e<nelem; ++e)
1030 : : {
1031 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
1032 : : // basis functions and solutions by default. This implies that P0P1 is
1033 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
1034 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
1035 : : // element for pDG.
1036 : : std::size_t dof_el;
1037 [ + - ]: 682200 : if (rdof > ndof)
1038 : : {
1039 : 682200 : dof_el = rdof;
1040 : : }
1041 : : else
1042 : : {
1043 : 0 : dof_el = ndofel[e];
1044 : : }
1045 : :
1046 [ + - ]: 682200 : if (dof_el > 1)
1047 : : {
1048 : 1364400 : std::vector< std::size_t > varc, varp;
1049 [ + - ]: 682200 : if(shockmarker[e]) {
1050 : : // When shockmarker is 1, there is discontinuity within the element.
1051 : : // Hence, the vertex-based limiter will be applied.
1052 : :
1053 [ + + ][ + - ]: 4434300 : for (std::size_t c=0; c<ncomp; ++c) varc.push_back(c);
1054 [ + + ][ + - ]: 1364400 : for (std::size_t c=0; c<nprim; ++c) varp.push_back(c);
1055 : : } else {
1056 : : // When shockmarker is 0, the density of minor species will still be
1057 : : // limited to ensure a stable solution.
1058 : :
1059 : 0 : tk::real rhob(0.0);
1060 [ - - ]: 0 : for(std::size_t k=0; k<nspec; ++k)
1061 [ - - ]: 0 : rhob += U(e, multispecies::densityDofIdx(nspec,k,rdof,0));
1062 [ - - ]: 0 : for(std::size_t k=0; k<nspec; ++k) {
1063 [ - - ][ - - ]: 0 : if (U(e, multispecies::densityDofIdx(nspec,k,rdof,0))/rhob < 1e-4) {
1064 : : // limit the density of minor species
1065 [ - - ]: 0 : varc.push_back(multispecies::densityIdx(nspec, k));
1066 : : }
1067 : : }
1068 : : }
1069 : :
1070 [ + - ][ + - ]: 1364400 : std::vector< tk::real > phic(ncomp, 1.0), phip(nprim, 1.0);
1071 [ + - ]: 682200 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
1072 : : ncomp, phic, varc);
1073 [ + - ]: 682200 : if (!varp.empty())
1074 [ + - ]: 682200 : VertexBasedLimiting(P, esup, inpoel, coord, e, rdof, dof_el,
1075 : : nprim, phip, varp);
1076 : :
1077 : 1364400 : std::vector< tk::real > phic_p2, phip_p2;
1078 [ + - ]: 682200 : PositivityLimiting(1, nspec, mat_blk, rdof, dof_el, ndofel, e,
1079 : : inpoel, coord, fd.Esuel(), U, P, phic, phic_p2, phip, phip_p2);
1080 : :
1081 : : // TODO: Unit sum of mass fractions is maintained by using common limiter
1082 : : // for all species densities. Investigate better approaches.
1083 [ + - ]: 682200 : if (!g_inputdeck.get< tag::accuracy_test >()) {
1084 : 682200 : tk::real phi_rhos_p1(1.0);
1085 [ + + ]: 1705500 : for (std::size_t k=0; k<nspec; ++k)
1086 : 1023300 : phi_rhos_p1 = std::min( phi_rhos_p1,
1087 : 1023300 : phic[multispecies::densityIdx(nspec, k)] );
1088 : : // same limiter for all densities
1089 [ + + ]: 1705500 : for (std::size_t k=0; k<nspec; ++k)
1090 : 1023300 : phic[multispecies::densityIdx(nspec, k)] = phi_rhos_p1;
1091 : : }
1092 : :
1093 : : // apply limiter function
1094 [ + + ]: 4434300 : for (std::size_t c=0; c<ncomp; ++c)
1095 : : {
1096 : 3752100 : auto mark = c*rdof;
1097 [ + - ][ + - ]: 3752100 : U(e, mark+1) = phic[c] * U(e, mark+1);
1098 [ + - ][ + - ]: 3752100 : U(e, mark+2) = phic[c] * U(e, mark+2);
1099 [ + - ][ + - ]: 3752100 : U(e, mark+3) = phic[c] * U(e, mark+3);
1100 : : }
1101 [ + + ]: 1364400 : for (std::size_t c=0; c<nprim; ++c)
1102 : : {
1103 : 682200 : auto mark = c*rdof;
1104 [ + - ][ + - ]: 682200 : P(e, mark+1) = phip[c] * P(e, mark+1);
1105 [ + - ][ + - ]: 682200 : P(e, mark+2) = phip[c] * P(e, mark+2);
1106 [ + - ][ + - ]: 682200 : P(e, mark+3) = phip[c] * P(e, mark+3);
1107 : : }
1108 : : }
1109 : : }
1110 : 6600 : }
1111 : :
1112 : : void
1113 : 0 : VertexBasedMultiSpecies_P2(
1114 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
1115 : : const std::vector< std::size_t >& inpoel,
1116 : : const std::vector< std::size_t >& ndofel,
1117 : : std::size_t nelem,
1118 : : const std::vector< inciter::EOS >& mat_blk,
1119 : : const inciter::FaceData& fd,
1120 : : const tk::Fields& geoFace,
1121 : : const tk::Fields& geoElem,
1122 : : const tk::UnsMesh::Coords& coord,
1123 : : const tk::FluxFn& flux,
1124 : : const std::vector< std::size_t >& solidx,
1125 : : tk::Fields& U,
1126 : : tk::Fields& P,
1127 : : std::size_t nspec,
1128 : : std::vector< std::size_t >& shockmarker )
1129 : : // *****************************************************************************
1130 : : // Kuzmin's vertex-based limiter for multi-species DGP2
1131 : : //! \param[in] esup Elements surrounding points
1132 : : //! \param[in] inpoel Element connectivity
1133 : : //! \param[in] ndofel Vector of local number of degrees of freedom
1134 : : //! \param[in] nelem Number of elements
1135 : : //! \param[in] mat_blk EOS material block
1136 : : //! \param[in] fd Face connectivity and boundary conditions object
1137 : : //! \param[in] geoFace Face geometry array
1138 : : //! \param[in] geoElem Element geometry array
1139 : : //! \param[in] coord Array of nodal coordinates
1140 : : //! \param[in] flux Riemann flux function to use
1141 : : //! \param[in] solidx Solid material index indicator
1142 : : //! \param[in,out] U High-order solution vector which gets limited
1143 : : //! \param[in,out] P High-order primitive vector which gets limited
1144 : : //! \param[in] nspec Number of species in this PDE system
1145 : : //! \param[in,out] shockmarker Shock detection marker array
1146 : : //! \details This vertex-based limiter function should be called for
1147 : : //! multispecies. For details see: Kuzmin, D. (2010). A vertex-based
1148 : : //! hierarchical slope limiter for p-adaptive discontinuous Galerkin methods.
1149 : : //! Journal of computational and applied mathematics, 233(12), 3077-3085.
1150 : : // *****************************************************************************
1151 : : {
1152 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
1153 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
1154 : 0 : std::size_t ncomp = U.nprop()/rdof;
1155 : 0 : std::size_t nprim = P.nprop()/rdof;
1156 : :
1157 : : // Evaluate the interface condition and mark the shock cells
1158 [ - - ]: 0 : if (inciter::g_inputdeck.get< tag::shock_detector_coeff >() > 1e-6) {
1159 : 0 : std::set< std::size_t > vars;
1160 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1161 [ - - ]: 0 : vars.insert(multispecies::momentumIdx(nspec, i));
1162 [ - - ]: 0 : MarkShockCells(false, nelem, 1, ndof, rdof, mat_blk,
1163 : : ndofel, inpoel, coord, fd, geoFace, geoElem, flux, solidx, U, P,
1164 : : vars, shockmarker);
1165 : : }
1166 : :
1167 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
1168 : : {
1169 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
1170 : : // basis functions and solutions by default. This implies that P0P1 is
1171 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
1172 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
1173 : : // element for pDG.
1174 : : std::size_t dof_el;
1175 [ - - ]: 0 : if (rdof > ndof)
1176 : : {
1177 : 0 : dof_el = rdof;
1178 : : }
1179 : : else
1180 : : {
1181 : 0 : dof_el = ndofel[e];
1182 : : }
1183 : :
1184 [ - - ]: 0 : if (dof_el > 1)
1185 : : {
1186 : 0 : std::vector< std::size_t > varc, varp;
1187 [ - - ]: 0 : if(shockmarker[e]) {
1188 : : // When shockmarker is 1, there is discontinuity within the element.
1189 : : // Hence, the vertex-based limiter will be applied.
1190 : :
1191 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) varc.push_back(c);
1192 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c) varp.push_back(c);
1193 : : }
1194 : : // When shockmarker is 0, the density of minor species is not limited
1195 : : // since it will degrade the accuracy to second order.
1196 : :
1197 : : // Removing 3rd order DOFs if discontinuity is detected, and applying
1198 : : // limiting to the 2nd order/P1 solution
1199 [ - - ]: 0 : for (std::size_t c=0; c<varc.size(); c++) {
1200 [ - - ]: 0 : for(std::size_t idof = 4; idof < rdof; idof++) {
1201 : 0 : auto mark = varc[c] * rdof + idof;
1202 [ - - ]: 0 : U(e, mark) = 0.0;
1203 : : }
1204 : : }
1205 [ - - ]: 0 : for (std::size_t c=0; c<varp.size(); c++) {
1206 [ - - ]: 0 : for(std::size_t idof = 4; idof < rdof; idof++) {
1207 : 0 : auto mark = varp[c] * rdof + idof;
1208 [ - - ]: 0 : P(e, mark) = 0.0;
1209 : : }
1210 : : }
1211 : :
1212 [ - - ][ - - ]: 0 : std::vector< tk::real > phic_p1(ncomp, 1.0), phip_p1(nprim, 1.0);
1213 [ - - ]: 0 : if (!varc.empty())
1214 [ - - ]: 0 : VertexBasedLimiting(U, esup, inpoel, coord, e, rdof, dof_el,
1215 : : ncomp, phic_p1, varc);
1216 [ - - ]: 0 : if (!varp.empty())
1217 [ - - ]: 0 : VertexBasedLimiting(P, esup, inpoel, coord, e, rdof, dof_el,
1218 : : nprim, phip_p1, varp);
1219 : :
1220 [ - - ][ - - ]: 0 : std::vector< tk::real > phic_p2(ncomp, 1.0), phip_p2(nprim, 1.0);
1221 [ - - ]: 0 : PositivityLimiting(1, nspec, mat_blk, ndof, dof_el, ndofel, e,
1222 : : inpoel, coord, fd.Esuel(), U, P, phic_p1, phic_p2, phip_p1, phip_p2);
1223 : :
1224 : : // TODO: Unit sum of mass fractions is maintained by using common limiter
1225 : : // for all species densities. Investigate better approaches.
1226 [ - - ]: 0 : if (!g_inputdeck.get< tag::accuracy_test >()) {
1227 : 0 : tk::real phi_rhos_p1(1.0), phi_rhos_p2(1.0);
1228 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
1229 : 0 : phi_rhos_p1 = std::min( phi_rhos_p1,
1230 : 0 : phic_p1[multispecies::densityIdx(nspec, k)] );
1231 : 0 : phi_rhos_p2 = std::min( phi_rhos_p2,
1232 : 0 : phic_p2[multispecies::densityIdx(nspec, k)] );
1233 : : }
1234 : : // same limiter for all densities
1235 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
1236 : 0 : phic_p1[multispecies::densityIdx(nspec, k)] = phi_rhos_p1;
1237 : 0 : phic_p2[multispecies::densityIdx(nspec, k)] = phi_rhos_p2;
1238 : : }
1239 : : }
1240 : :
1241 : : // apply limiter function to the solution
1242 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
1243 : : {
1244 : 0 : auto mark = c * rdof;
1245 [ - - ]: 0 : for(std::size_t idof=1; idof<4; idof++)
1246 [ - - ][ - - ]: 0 : U(e, mark+idof) = phic_p1[c] * U(e, mark+idof);
1247 [ - - ]: 0 : for(std::size_t idof=4; idof<rdof; idof++)
1248 [ - - ][ - - ]: 0 : U(e, mark+idof) = phic_p2[c] * U(e, mark+idof);
1249 : : }
1250 [ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c)
1251 : : {
1252 : 0 : auto mark = c * rdof;
1253 [ - - ]: 0 : for(std::size_t idof=1; idof<4; idof++)
1254 [ - - ][ - - ]: 0 : P(e, mark+idof) = phip_p1[c] * P(e, mark+idof);
1255 [ - - ]: 0 : for(std::size_t idof=4; idof<rdof; idof++)
1256 [ - - ][ - - ]: 0 : P(e, mark+idof) = phip_p2[c] * P(e, mark+idof);
1257 : : }
1258 : : }
1259 : : }
1260 : 0 : }
1261 : :
1262 : : void
1263 : 0 : WENOLimiting( const tk::Fields& U,
1264 : : const std::vector< int >& esuel,
1265 : : std::size_t e,
1266 : : inciter::ncomp_t c,
1267 : : std::size_t rdof,
1268 : : tk::real cweight,
1269 : : std::array< std::vector< tk::real >, 3 >& limU )
1270 : : // *****************************************************************************
1271 : : // WENO limiter function calculation for P1 dofs
1272 : : //! \param[in] U High-order solution vector which is to be limited
1273 : : //! \param[in] esuel Elements surrounding elements
1274 : : //! \param[in] e Id of element whose solution is to be limited
1275 : : //! \param[in] c Index of component which is to be limited
1276 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
1277 : : //! \param[in] cweight Weight of the central stencil
1278 : : //! \param[in,out] limU Limited gradients of component c
1279 : : // *****************************************************************************
1280 : : {
1281 : : std::array< std::array< tk::real, 3 >, 5 > gradu;
1282 : : std::array< tk::real, 5 > wtStencil, osc, wtDof;
1283 : :
1284 : 0 : auto mark = c*rdof;
1285 : :
1286 : : // reset all stencil values to zero
1287 [ - - ][ - - ]: 0 : for (auto& g : gradu) g.fill(0.0);
1288 [ - - ]: 0 : osc.fill(0);
1289 [ - - ]: 0 : wtDof.fill(0);
1290 [ - - ]: 0 : wtStencil.fill(0);
1291 : :
1292 : : // The WENO limiter uses solution data from the neighborhood in the form
1293 : : // of stencils to enforce non-oscillatory conditions. The immediate
1294 : : // (Von Neumann) neighborhood of a tetrahedral cell consists of the 4
1295 : : // cells that share faces with it. These are the 4 neighborhood-stencils
1296 : : // for the tetrahedron. The primary stencil is the tet itself. Weights are
1297 : : // assigned to these stencils, with the primary stencil usually assigned
1298 : : // the highest weight. The lower the primary/central weight, the more
1299 : : // dissipative the limiting effect. This central weight is usually problem
1300 : : // dependent. It is set higher for relatively weaker discontinuities, and
1301 : : // lower for stronger discontinuities.
1302 : :
1303 : : // primary stencil
1304 [ - - ]: 0 : gradu[0][0] = U(e, mark+1);
1305 [ - - ]: 0 : gradu[0][1] = U(e, mark+2);
1306 [ - - ]: 0 : gradu[0][2] = U(e, mark+3);
1307 : 0 : wtStencil[0] = cweight;
1308 : :
1309 : : // stencils from the neighborhood
1310 [ - - ]: 0 : for (std::size_t is=1; is<5; ++is)
1311 : : {
1312 : 0 : auto nel = esuel[ 4*e+(is-1) ];
1313 : :
1314 : : // ignore physical domain ghosts
1315 [ - - ]: 0 : if (nel == -1)
1316 : : {
1317 [ - - ]: 0 : gradu[is].fill(0.0);
1318 : 0 : wtStencil[is] = 0.0;
1319 : 0 : continue;
1320 : : }
1321 : :
1322 : 0 : std::size_t n = static_cast< std::size_t >( nel );
1323 [ - - ]: 0 : gradu[is][0] = U(n, mark+1);
1324 [ - - ]: 0 : gradu[is][1] = U(n, mark+2);
1325 [ - - ]: 0 : gradu[is][2] = U(n, mark+3);
1326 : 0 : wtStencil[is] = 1.0;
1327 : : }
1328 : :
1329 : : // From these stencils, an oscillation indicator is calculated, which
1330 : : // determines the effective weights for the high-order solution DOFs.
1331 : : // These effective weights determine the contribution of each of the
1332 : : // stencils to the high-order solution DOFs of the current cell which are
1333 : : // being limited. If this indicator detects a large oscillation in the
1334 : : // solution of the current cell, it reduces the effective weight for the
1335 : : // central stencil contribution to its high-order DOFs. This results in
1336 : : // a more dissipative and well-behaved solution in the troubled cell.
1337 : :
1338 : : // oscillation indicators
1339 [ - - ]: 0 : for (std::size_t is=0; is<5; ++is)
1340 : 0 : osc[is] = std::sqrt( tk::dot(gradu[is], gradu[is]) );
1341 : :
1342 : 0 : tk::real wtotal = 0;
1343 : :
1344 : : // effective weights for dofs
1345 [ - - ]: 0 : for (std::size_t is=0; is<5; ++is)
1346 : : {
1347 : : // A small number (1.0e-8) is needed here to avoid dividing by a zero in
1348 : : // the case of a constant solution, where osc would be zero. The number
1349 : : // is not set to machine zero because it is squared, and a number
1350 : : // between 1.0e-8 to 1.0e-6 is needed.
1351 : 0 : wtDof[is] = wtStencil[is] * pow( (1.0e-8 + osc[is]), -2 );
1352 : 0 : wtotal += wtDof[is];
1353 : : }
1354 : :
1355 [ - - ]: 0 : for (std::size_t is=0; is<5; ++is)
1356 : : {
1357 : 0 : wtDof[is] = wtDof[is]/wtotal;
1358 : : }
1359 : :
1360 : 0 : limU[0][e] = 0.0;
1361 : 0 : limU[1][e] = 0.0;
1362 : 0 : limU[2][e] = 0.0;
1363 : :
1364 : : // limiter function
1365 [ - - ]: 0 : for (std::size_t is=0; is<5; ++is)
1366 : : {
1367 : 0 : limU[0][e] += wtDof[is]*gradu[is][0];
1368 : 0 : limU[1][e] += wtDof[is]*gradu[is][1];
1369 : 0 : limU[2][e] += wtDof[is]*gradu[is][2];
1370 : : }
1371 : 0 : }
1372 : :
1373 : : std::vector< tk::real >
1374 : 476319 : SuperbeeLimiting( const tk::Fields& U,
1375 : : const std::vector< int >& esuel,
1376 : : const std::vector< std::size_t >& inpoel,
1377 : : const tk::UnsMesh::Coords& coord,
1378 : : std::size_t e,
1379 : : std::size_t ndof,
1380 : : std::size_t rdof,
1381 : : std::size_t dof_el,
1382 : : inciter:: ncomp_t ncomp,
1383 : : tk::real beta_lim )
1384 : : // *****************************************************************************
1385 : : // Superbee limiter function calculation for P1 dofs
1386 : : //! \param[in] U High-order solution vector which is to be limited
1387 : : //! \param[in] esuel Elements surrounding elements
1388 : : //! \param[in] inpoel Element connectivity
1389 : : //! \param[in] coord Array of nodal coordinates
1390 : : //! \param[in] e Id of element whose solution is to be limited
1391 : : //! \param[in] ndof Maximum number of degrees of freedom
1392 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
1393 : : //! \param[in] dof_el Local number of degrees of freedom
1394 : : //! \param[in] ncomp Number of scalar components in this PDE system
1395 : : //! \param[in] beta_lim Parameter which is equal to 2 for Superbee and 1 for
1396 : : //! minmod limiter
1397 : : //! \return phi Limiter function for solution in element e
1398 : : // *****************************************************************************
1399 : : {
1400 : : // Superbee is a TVD limiter, which uses min-max bounds that the
1401 : : // high-order solution should satisfy, to ensure TVD properties. For a
1402 : : // high-order method like DG, this involves the following steps:
1403 : : // 1. Find min-max bounds in the immediate neighborhood of cell.
1404 : : // 2. Calculate the Superbee function for all the points where solution
1405 : : // needs to be reconstructed to (all quadrature points). From these,
1406 : : // use the minimum value of the limiter function.
1407 : :
1408 [ + - ][ + - ]: 952638 : std::vector< tk::real > uMin(ncomp, 0.0), uMax(ncomp, 0.0);
1409 : :
1410 [ + + ]: 2857914 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
1411 : : {
1412 : 2381595 : auto mark = c*rdof;
1413 [ + - ]: 2381595 : uMin[c] = U(e, mark);
1414 [ + - ]: 2381595 : uMax[c] = U(e, mark);
1415 : : }
1416 : :
1417 : : // ----- Step-1: find min/max in the neighborhood
1418 [ + + ]: 2381595 : for (std::size_t is=0; is<4; ++is)
1419 : : {
1420 : 1905276 : auto nel = esuel[ 4*e+is ];
1421 : :
1422 : : // ignore physical domain ghosts
1423 [ + + ]: 1905276 : if (nel == -1) continue;
1424 : :
1425 : 1591101 : auto n = static_cast< std::size_t >( nel );
1426 [ + + ]: 9546606 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
1427 : : {
1428 : 7955505 : auto mark = c*rdof;
1429 [ + - ]: 7955505 : uMin[c] = std::min(uMin[c], U(n, mark));
1430 [ + - ]: 7955505 : uMax[c] = std::max(uMax[c], U(n, mark));
1431 : : }
1432 : : }
1433 : :
1434 : : // ----- Step-2: loop over all quadrature points to get limiter function
1435 : :
1436 : : // to loop over all the quadrature points of all faces of element e,
1437 : : // coordinates of the quadrature points are needed.
1438 : : // Number of quadrature points for face integration
1439 [ + - ]: 476319 : auto ng = tk::NGfa(ndof);
1440 : :
1441 : : // arrays for quadrature points
1442 : 952638 : std::array< std::vector< tk::real >, 2 > coordgp;
1443 : 952638 : std::vector< tk::real > wgp;
1444 : :
1445 [ + - ]: 476319 : coordgp[0].resize( ng );
1446 [ + - ]: 476319 : coordgp[1].resize( ng );
1447 [ + - ]: 476319 : wgp.resize( ng );
1448 : :
1449 : : // get quadrature point weights and coordinates for triangle
1450 [ + - ]: 476319 : tk::GaussQuadratureTri( ng, coordgp, wgp );
1451 : :
1452 : 476319 : const auto& cx = coord[0];
1453 : 476319 : const auto& cy = coord[1];
1454 : 476319 : const auto& cz = coord[2];
1455 : :
1456 : : // Extract the element coordinates
1457 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
1458 : 1428957 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
1459 : 1428957 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
1460 : 1428957 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
1461 : 4286871 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
1462 : :
1463 : : // Compute the determinant of Jacobian matrix
1464 : : auto detT =
1465 : 476319 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
1466 : :
1467 : : // initialize limiter function
1468 [ + - ]: 476319 : std::vector< tk::real > phi(ncomp, 1.0);
1469 [ + + ]: 2381595 : for (std::size_t lf=0; lf<4; ++lf)
1470 : : {
1471 : : // Extract the face coordinates
1472 : 1905276 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
1473 : 1905276 : inpoel[4*e+tk::lpofa[lf][1]],
1474 : 3810552 : inpoel[4*e+tk::lpofa[lf][2]] }};
1475 : :
1476 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
1477 : 5715828 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
1478 : 5715828 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
1479 : 11431656 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
1480 : :
1481 : : // Gaussian quadrature
1482 [ + + ]: 8091012 : for (std::size_t igp=0; igp<ng; ++igp)
1483 : : {
1484 : : // Compute the coordinates of quadrature point at physical domain
1485 [ + - ]: 6185736 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
1486 : :
1487 : : //Compute the basis functions
1488 : : auto B_l = tk::eval_basis( rdof,
1489 : 6185736 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
1490 : 6185736 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
1491 [ + - ]: 24742944 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
1492 : :
1493 : : auto state =
1494 [ + - ]: 12371472 : tk::eval_state(ncomp, rdof, dof_el, e, U, B_l);
1495 : :
1496 [ - + ][ - - ]: 6185736 : Assert( state.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
1497 : :
1498 : : // compute the limiter function
1499 [ + + ]: 37114416 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
1500 : : {
1501 : 30928680 : auto phi_gp = 1.0;
1502 : 30928680 : auto mark = c*rdof;
1503 [ + - ]: 30928680 : auto uNeg = state[c] - U(e, mark);
1504 [ + + ]: 30928680 : if (uNeg > 1.0e-14)
1505 : : {
1506 : 468201 : uNeg = std::max(uNeg, 1.0e-08);
1507 [ + - ]: 468201 : phi_gp = std::min( 1.0, (uMax[c]-U(e, mark))/(2.0*uNeg) );
1508 : : }
1509 [ + + ]: 30460479 : else if (uNeg < -1.0e-14)
1510 : : {
1511 : 438960 : uNeg = std::min(uNeg, -1.0e-08);
1512 [ + - ]: 438960 : phi_gp = std::min( 1.0, (uMin[c]-U(e, mark))/(2.0*uNeg) );
1513 : : }
1514 : : else
1515 : : {
1516 : 30021519 : phi_gp = 1.0;
1517 : : }
1518 : 61857360 : phi_gp = std::max( 0.0,
1519 : 61857360 : std::max( std::min(beta_lim*phi_gp, 1.0),
1520 : 30928680 : std::min(phi_gp, beta_lim) ) );
1521 : 30928680 : phi[c] = std::min( phi[c], phi_gp );
1522 : : }
1523 : : }
1524 : : }
1525 : :
1526 : 952638 : return phi;
1527 : : }
1528 : :
1529 : : void
1530 : 4014384 : VertexBasedLimiting(
1531 : : const tk::Fields& U,
1532 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
1533 : : const std::vector< std::size_t >& inpoel,
1534 : : const tk::UnsMesh::Coords& coord,
1535 : : std::size_t e,
1536 : : std::size_t rdof,
1537 : : std::size_t dof_el,
1538 : : std::size_t ncomp,
1539 : : std::vector< tk::real >& phi,
1540 : : const std::vector< std::size_t >& VarList )
1541 : : // *****************************************************************************
1542 : : // Kuzmin's vertex-based limiter function calculation for P1 dofs
1543 : : //! \param[in] U High-order solution vector which is to be limited
1544 : : //! \param[in] esup Elements surrounding points
1545 : : //! \param[in] inpoel Element connectivity
1546 : : //! \param[in] coord Array of nodal coordinates
1547 : : //! \param[in] e Id of element whose solution is to be limited
1548 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
1549 : : //! \param[in] dof_el Local number of degrees of freedom
1550 : : //! \param[in] ncomp Number of scalar components in this PDE system
1551 : : //! \param[in,out] phi Limiter function for solution in element e
1552 : : //! \param[in] VarList List of variable indices to be limited
1553 : : // *****************************************************************************
1554 : : {
1555 : : // Kuzmin's vertex-based TVD limiter uses min-max bounds that the
1556 : : // high-order solution should satisfy, to ensure TVD properties. For a
1557 : : // high-order method like DG, this involves the following steps:
1558 : : // 1. Find min-max bounds in the nodal-neighborhood of cell.
1559 : : // 2. Calculate the limiter function (Superbee) for all the vertices of cell.
1560 : : // From these, use the minimum value of the limiter function.
1561 : :
1562 : : // Prepare for calculating Basis functions
1563 : 4014384 : const auto& cx = coord[0];
1564 : 4014384 : const auto& cy = coord[1];
1565 : 4014384 : const auto& cz = coord[2];
1566 : :
1567 : : // Extract the element coordinates
1568 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
1569 : 12043152 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
1570 : 12043152 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
1571 : 12043152 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
1572 : 36129456 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
1573 : :
1574 : : // Compute the determinant of Jacobian matrix
1575 : : auto detT =
1576 : 4014384 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
1577 : :
1578 [ + - ]: 8028768 : std::vector< tk::real > uMin(VarList.size(), 0.0),
1579 [ + - ]: 8028768 : uMax(VarList.size(), 0.0);
1580 : :
1581 : : // loop over all nodes of the element e
1582 [ + + ]: 20071920 : for (std::size_t lp=0; lp<4; ++lp)
1583 : : {
1584 : : // reset min/max
1585 [ + + ]: 76248656 : for (std::size_t i=0; i<VarList.size(); ++i)
1586 : : {
1587 : 60191120 : auto mark = VarList[i]*rdof;
1588 [ + - ]: 60191120 : uMin[i] = U(e, mark);
1589 [ + - ]: 60191120 : uMax[i] = U(e, mark);
1590 : : }
1591 : 16057536 : auto p = inpoel[4*e+lp];
1592 [ + - ]: 16057536 : const auto& pesup = tk::cref_find(esup, p);
1593 : :
1594 : : // ----- Step-1: find min/max in the neighborhood of node p
1595 : : // loop over all the internal elements surrounding this node p
1596 [ + + ]: 266573880 : for (auto er : pesup)
1597 : : {
1598 [ + + ]: 1192790316 : for (std::size_t i=0; i<VarList.size(); ++i)
1599 : : {
1600 : 942273972 : auto mark = VarList[i]*rdof;
1601 [ + - ]: 942273972 : uMin[i] = std::min(uMin[i], U(er, mark));
1602 [ + - ]: 942273972 : uMax[i] = std::max(uMax[i], U(er, mark));
1603 : : }
1604 : : }
1605 : :
1606 : : // ----- Step-2: compute the limiter function at this node
1607 : : // find high-order solution
1608 : 32115072 : std::vector< tk::real > state;
1609 : 16057536 : std::array< tk::real, 3 > gp{cx[p], cy[p], cz[p]};
1610 : : auto B_p = tk::eval_basis( rdof,
1611 : 16057536 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
1612 : 16057536 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
1613 [ + - ]: 64230144 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
1614 [ + - ]: 16057536 : state = tk::eval_state(ncomp, rdof, dof_el, e, U, B_p);
1615 : :
1616 [ - + ][ - - ]: 16057536 : Assert( state.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
1617 : :
1618 : : // compute the limiter function
1619 [ + + ]: 76248656 : for (std::size_t i=0; i<VarList.size(); ++i)
1620 : : {
1621 : 60191120 : auto c = VarList[i];
1622 : 60191120 : auto phi_gp = 1.0;
1623 : 60191120 : auto mark = c*rdof;
1624 [ + - ]: 60191120 : auto uNeg = state[c] - U(e, mark);
1625 [ + - ]: 60191120 : auto uref = std::max(std::fabs(U(e,mark)), 1e-14);
1626 [ + + ]: 60191120 : if (uNeg > 1.0e-06*uref)
1627 : : {
1628 [ + - ]: 9064131 : phi_gp = std::min( 1.0, (uMax[i]-U(e, mark))/uNeg );
1629 : : }
1630 [ + + ]: 51126989 : else if (uNeg < -1.0e-06*uref)
1631 : : {
1632 [ + - ]: 9097787 : phi_gp = std::min( 1.0, (uMin[i]-U(e, mark))/uNeg );
1633 : : }
1634 : : else
1635 : : {
1636 : 42029202 : phi_gp = 1.0;
1637 : : }
1638 : :
1639 : : // ----- Step-3: take the minimum of the nodal-limiter functions
1640 : 60191120 : phi[c] = std::min( phi[c], phi_gp );
1641 : : }
1642 : : }
1643 : 4014384 : }
1644 : :
1645 : : void
1646 : 0 : VertexBasedLimiting_P2( const std::vector< std::vector< tk::real > >& unk,
1647 : : const tk::Fields& U,
1648 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
1649 : : const std::vector< std::size_t >& inpoel,
1650 : : std::size_t e,
1651 : : std::size_t rdof,
1652 : : [[maybe_unused]] std::size_t dof_el,
1653 : : std::size_t ncomp,
1654 : : const std::vector< std::size_t >& gid,
1655 : : const std::unordered_map< std::size_t, std::size_t >& bid,
1656 : : const std::vector< std::vector<tk::real> >& NodalExtrm,
1657 : : const std::vector< std::size_t >& VarList,
1658 : : std::vector< tk::real >& phi )
1659 : : // *****************************************************************************
1660 : : // Kuzmin's vertex-based limiter function calculation for P2 dofs
1661 : : //! \param[in] U High-order solution vector which is to be limited
1662 : : //! \param[in] esup Elements surrounding points
1663 : : //! \param[in] inpoel Element connectivity
1664 : : //! \param[in] e Id of element whose solution is to be limited
1665 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
1666 : : //! \param[in] dof_el Local number of degrees of freedom
1667 : : //! \param[in] ncomp Number of scalar components in this PDE system
1668 : : //! \param[in] gid Local->global node id map
1669 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
1670 : : //! global node ids (key)
1671 : : //! \param[in] NodalExtrm Chare-boundary nodal extrema
1672 : : //! \param[in] VarList List of variable indices that need to be limited
1673 : : //! \param[out] phi Limiter function for solution in element e
1674 : : //! \details This function limits the P2 dofs of P2 solution in a hierachical
1675 : : //! way to P1 dof limiting. Here we treat the first order derivatives the same
1676 : : //! way as cell average while second order derivatives represent the gradients
1677 : : //! to be limited in the P1 limiting procedure.
1678 : : // *****************************************************************************
1679 : : {
1680 : 0 : const auto nelem = inpoel.size() / 4;
1681 : :
1682 : 0 : std::vector< std::vector< tk::real > > uMin, uMax;
1683 [ - - ][ - - ]: 0 : uMin.resize( VarList.size(), std::vector<tk::real>(3, 0.0) );
1684 [ - - ][ - - ]: 0 : uMax.resize( VarList.size(), std::vector<tk::real>(3, 0.0) );
1685 : :
1686 : : // The coordinates of centroid in the reference domain
1687 : 0 : std::array< std::vector< tk::real >, 3 > center;
1688 [ - - ]: 0 : center[0].resize(1, 0.25);
1689 [ - - ]: 0 : center[1].resize(1, 0.25);
1690 [ - - ]: 0 : center[2].resize(1, 0.25);
1691 : :
1692 : 0 : std::array< std::array< tk::real, 4 >, 3 > cnodes{{
1693 : : {{0, 1, 0, 0}},
1694 : : {{0, 0, 1, 0}},
1695 : : {{0, 0, 0, 1}} }};
1696 : :
1697 : : // loop over all nodes of the element e
1698 [ - - ]: 0 : for (std::size_t lp=0; lp<4; ++lp)
1699 : : {
1700 : : // Find the max/min first-order derivatives for internal element
1701 [ - - ]: 0 : for (std::size_t i=0; i<VarList.size(); ++i)
1702 : : {
1703 [ - - ]: 0 : for (std::size_t idir=1; idir < 4; ++idir)
1704 : : {
1705 : 0 : uMin[i][idir-1] = unk[VarList[i]][idir];
1706 : 0 : uMax[i][idir-1] = unk[VarList[i]][idir];
1707 : : }
1708 : : }
1709 : :
1710 : 0 : auto p = inpoel[4*e+lp];
1711 [ - - ]: 0 : const auto& pesup = tk::cref_find(esup, p);
1712 : :
1713 : : // Step-1: find min/max first order derivative at the centroid in the
1714 : : // neighborhood of node p
1715 [ - - ]: 0 : for (auto er : pesup)
1716 : : {
1717 [ - - ]: 0 : if(er < nelem) // If this is internal element
1718 : : {
1719 : : // Compute the derivatives of basis function in the reference domain
1720 : : auto dBdxi_er = tk::eval_dBdxi(rdof,
1721 [ - - ]: 0 : {{center[0][0], center[1][0], center[2][0]}});
1722 : :
1723 [ - - ]: 0 : for (std::size_t i=0; i<VarList.size(); ++i)
1724 : : {
1725 : 0 : auto mark = VarList[i]*rdof;
1726 [ - - ]: 0 : for (std::size_t idir = 0; idir < 3; ++idir)
1727 : : {
1728 : : // The first order derivative at the centroid of element er
1729 : 0 : tk::real slope_er(0.0);
1730 [ - - ]: 0 : for(std::size_t idof = 1; idof < rdof; idof++)
1731 [ - - ]: 0 : slope_er += U(er, mark+idof) * dBdxi_er[idir][idof];
1732 : :
1733 : 0 : uMin[i][idir] = std::min(uMin[i][idir], slope_er);
1734 : 0 : uMax[i][idir] = std::max(uMax[i][idir], slope_er);
1735 : :
1736 : : }
1737 : : }
1738 : : }
1739 : : }
1740 : : // If node p is the chare-boundary node, find min/max by comparing with
1741 : : // the chare-boundary nodal extrema from vector NodalExtrm
1742 [ - - ]: 0 : auto gip = bid.find( gid[p] );
1743 [ - - ]: 0 : if(gip != end(bid))
1744 : : {
1745 : 0 : auto ndof_NodalExtrm = NodalExtrm[0].size() / (ncomp * 2);
1746 [ - - ]: 0 : for (std::size_t i=0; i<VarList.size(); ++i)
1747 : : {
1748 [ - - ]: 0 : for (std::size_t idir = 0; idir < 3; idir++)
1749 : : {
1750 : 0 : auto max_mark = 2*VarList[i]*ndof_NodalExtrm + 2*idir;
1751 : 0 : auto min_mark = max_mark + 1;
1752 : 0 : const auto& ex = NodalExtrm[gip->second];
1753 : 0 : uMax[i][idir] = std::max(ex[max_mark], uMax[i][idir]);
1754 : 0 : uMin[i][idir] = std::min(ex[min_mark], uMin[i][idir]);
1755 : : }
1756 : : }
1757 : : }
1758 : :
1759 : : //Step-2: compute the limiter function at this node
1760 : 0 : std::array< tk::real, 3 > node{cnodes[0][lp], cnodes[1][lp], cnodes[2][lp]};
1761 : :
1762 : : // find high-order solution
1763 : 0 : std::vector< std::array< tk::real, 3 > > state;
1764 [ - - ]: 0 : state.resize(VarList.size());
1765 : :
1766 [ - - ]: 0 : for (std::size_t i=0; i<VarList.size(); ++i)
1767 : : {
1768 : 0 : auto dx = node[0] - center[0][0];
1769 : 0 : auto dy = node[1] - center[1][0];
1770 : 0 : auto dz = node[2] - center[2][0];
1771 : :
1772 : 0 : auto c = VarList[i];
1773 : :
1774 : 0 : state[i][0] = unk[c][1] + unk[c][4]*dx + unk[c][7]*dy + unk[c][8]*dz;
1775 : 0 : state[i][1] = unk[c][2] + unk[c][5]*dy + unk[c][7]*dx + unk[c][9]*dz;
1776 : 0 : state[i][2] = unk[c][3] + unk[c][6]*dz + unk[c][8]*dx + unk[c][9]*dy;
1777 : : }
1778 : :
1779 : : // compute the limiter function
1780 [ - - ]: 0 : for (std::size_t i=0; i<VarList.size(); ++i)
1781 : : {
1782 : 0 : auto c = VarList[i];
1783 : 0 : tk::real phi_dir(1.0);
1784 [ - - ]: 0 : for (std::size_t idir = 1; idir <= 3; ++idir)
1785 : : {
1786 : 0 : phi_dir = 1.0;
1787 : 0 : auto uNeg = state[i][idir-1] - unk[c][idir];
1788 : 0 : auto uref = std::max(std::fabs(unk[c][idir]), 1e-14);
1789 [ - - ]: 0 : if (uNeg > 1.0e-6*uref)
1790 : : {
1791 : 0 : phi_dir =
1792 : 0 : std::min( 1.0, ( uMax[i][idir-1] - unk[c][idir])/uNeg );
1793 : : }
1794 [ - - ]: 0 : else if (uNeg < -1.0e-6*uref)
1795 : : {
1796 : 0 : phi_dir =
1797 : 0 : std::min( 1.0, ( uMin[i][idir-1] - unk[c][idir])/uNeg );
1798 : : }
1799 : : else
1800 : : {
1801 : 0 : phi_dir = 1.0;
1802 : : }
1803 : :
1804 : 0 : phi[c] = std::min( phi[c], phi_dir );
1805 : : }
1806 : : }
1807 : : }
1808 : 0 : }
1809 : :
1810 : 1092845 : void consistentMultiMatLimiting_P1(
1811 : : std::size_t nmat,
1812 : : std::size_t rdof,
1813 : : std::size_t e,
1814 : : const std::vector< std::size_t >& solidx,
1815 : : tk::Fields& U,
1816 : : [[maybe_unused]] tk::Fields& P,
1817 : : std::vector< tk::real >& phic_p1,
1818 : : std::vector< tk::real >& phic_p2 )
1819 : : // *****************************************************************************
1820 : : // Consistent limiter modifications for conservative variables
1821 : : //! \param[in] nmat Number of materials in this PDE system
1822 : : //! \param[in] rdof Total number of reconstructed dofs
1823 : : //! \param[in] e Element being checked for consistency
1824 : : //! \param[in] solidx Solid material index indicator
1825 : : //! \param[in] U Vector of conservative variables
1826 : : //! \param[in] P Vector of primitive variables
1827 : : //! \param[in,out] phic_p1 Vector of limiter functions for P1 dofs of the
1828 : : //! conserved quantities
1829 : : //! \param[in,out] phip_p2 Vector of limiter functions for P2 dofs of the
1830 : : //! conserved quantities
1831 : : // *****************************************************************************
1832 : : {
1833 : : // find the limiter-function for volume-fractions
1834 : 1092845 : auto phi_al_p1(1.0), phi_al_p2(1.0), almax(0.0), dalmax(0.0);
1835 : : //std::size_t nmax(0);
1836 [ + + ]: 3331095 : for (std::size_t k=0; k<nmat; ++k)
1837 : : {
1838 : 2238250 : phi_al_p1 = std::min( phi_al_p1, phic_p1[volfracIdx(nmat, k)] );
1839 [ - + ]: 2238250 : if(rdof > 4)
1840 : 0 : phi_al_p2 = std::min( phi_al_p2, phic_p2[volfracIdx(nmat, k)] );
1841 [ + - ][ + + ]: 2238250 : if (almax < U(e,volfracDofIdx(nmat, k, rdof, 0)))
1842 : : {
1843 : : //nmax = k;
1844 [ + - ]: 1573036 : almax = U(e,volfracDofIdx(nmat, k, rdof, 0));
1845 : : }
1846 : 2238250 : tk::real dmax(0.0);
1847 : 2238250 : dmax = std::max(
1848 : : std::max(
1849 [ + - ]: 2238250 : std::abs(U(e,volfracDofIdx(nmat, k, rdof, 1))),
1850 [ + - ]: 2238250 : std::abs(U(e,volfracDofIdx(nmat, k, rdof, 2))) ),
1851 [ + - ]: 4476500 : std::abs(U(e,volfracDofIdx(nmat, k, rdof, 3))) );
1852 : 2238250 : dalmax = std::max( dalmax, dmax );
1853 : : }
1854 : :
1855 : 1092845 : auto al_band = 1e-4;
1856 : :
1857 : : //phi_al = phic[nmax];
1858 : :
1859 : : // determine if cell is a material-interface cell based on ad-hoc tolerances.
1860 : : // if interface-cell, then modify high-order dofs of conserved unknowns
1861 : : // consistently and use same limiter for all equations.
1862 : : // Slopes of solution variables \alpha_k \rho_k and \alpha_k \rho_k E_k need
1863 : : // to be modified in interface cells, such that slopes in the \rho_k and
1864 : : // \rho_k E_k part are ignored and only slopes in \alpha_k are considered.
1865 : : // Ideally, we would like to not do this, but this is a necessity to avoid
1866 : : // limiter-limiter interactions in multiphase CFD (see "K.-M. Shyue, F. Xiao,
1867 : : // An Eulerian interface sharpening algorithm for compressible two-phase flow:
1868 : : // the algebraic THINC approach, Journal of Computational Physics 268, 2014,
1869 : : // 326–354. doi:10.1016/j.jcp.2014.03.010." and "A. Chiapolino, R. Saurel,
1870 : : // B. Nkonga, Sharpening diffuse interfaces with compressible fluids on
1871 : : // unstructured meshes, Journal of Computational Physics 340 (2017) 389–417.
1872 : : // doi:10.1016/j.jcp.2017.03.042."). This approximation should be applied in
1873 : : // as narrow a band of interface-cells as possible. The following if-test
1874 : : // defines this band of interface-cells. This tests checks the value of the
1875 : : // maximum volume-fraction in the cell (almax) and the maximum change in
1876 : : // volume-fraction in the cell (dalmax, calculated from second-order DOFs),
1877 : : // to determine the band of interface-cells where the aforementioned fix needs
1878 : : // to be applied. This if-test says that, the fix is applied when the change
1879 : : // in volume-fraction across a cell is greater than 0.1, *and* the
1880 : : // volume-fraction is between 0.1 and 0.9.
1881 [ + - ]: 1092845 : if ( //dalmax > al_band &&
1882 [ + + ]: 1092845 : (almax > al_band && almax < (1.0-al_band)) )
1883 : : {
1884 : : // 1. consistent high-order dofs
1885 [ + + ]: 78426 : for (std::size_t k=0; k<nmat; ++k)
1886 : : {
1887 : : auto alk =
1888 [ + - ]: 52350 : std::max( 1.0e-14, U(e,volfracDofIdx(nmat, k, rdof, 0)) );
1889 [ + - ]: 52350 : auto rhok = U(e,densityDofIdx(nmat, k, rdof, 0)) / alk;
1890 [ + - ]: 52350 : auto rhoE = U(e,energyDofIdx(nmat, k, rdof, 0)) / alk;
1891 [ + + ]: 209400 : for (std::size_t idof=1; idof<rdof; ++idof)
1892 : : {
1893 [ + - ]: 157050 : U(e,densityDofIdx(nmat, k, rdof, idof)) = rhok *
1894 [ + - ]: 157050 : U(e,volfracDofIdx(nmat, k, rdof, idof));
1895 [ + - ]: 157050 : U(e,energyDofIdx(nmat, k, rdof, idof)) = rhoE *
1896 [ + - ]: 157050 : U(e,volfracDofIdx(nmat, k, rdof, idof));
1897 : : }
1898 [ - + ]: 52350 : if (solidx[k] > 0)
1899 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1900 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1901 : : {
1902 [ - - ]: 0 : for (std::size_t idof=1; idof<rdof; ++idof)
1903 [ - - ]: 0 : U(e,deformDofIdx(nmat,solidx[k],i,j,rdof,idof)) = 0.0;
1904 : : }
1905 : : }
1906 : :
1907 : : // 2. same limiter for all volume-fractions and densities
1908 [ + + ]: 78426 : for (std::size_t k=0; k<nmat; ++k)
1909 : : {
1910 : 52350 : phic_p1[volfracIdx(nmat, k)] = phi_al_p1;
1911 : 52350 : phic_p1[densityIdx(nmat, k)] = phi_al_p1;
1912 : 52350 : phic_p1[energyIdx(nmat, k)] = phi_al_p1;
1913 [ - + ]: 52350 : if (solidx[k] > 0)
1914 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1915 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1916 : 0 : phic_p1[deformIdx(nmat,solidx[k],i,j)] = phi_al_p1;
1917 : : }
1918 [ - + ]: 26076 : if(rdof > 4)
1919 : : {
1920 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1921 : : {
1922 : 0 : phic_p2[volfracIdx(nmat, k)] = phi_al_p2;
1923 : 0 : phic_p2[densityIdx(nmat, k)] = phi_al_p2;
1924 : 0 : phic_p2[energyIdx(nmat, k)] = phi_al_p2;
1925 [ - - ]: 0 : if (solidx[k] > 0)
1926 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1927 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1928 : 0 : phic_p2[deformIdx(nmat,solidx[k],i,j)] = phi_al_p2;
1929 : : }
1930 : 26076 : }
1931 : : }
1932 : : else
1933 : : {
1934 : : // same limiter for all volume-fractions
1935 [ + + ]: 3252669 : for (std::size_t k=0; k<nmat; ++k)
1936 : 2185900 : phic_p1[volfracIdx(nmat, k)] = phi_al_p1;
1937 [ - + ]: 1066769 : if(rdof > 4)
1938 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1939 : 0 : phic_p2[volfracIdx(nmat, k)] = phi_al_p2;
1940 : : }
1941 : 1092845 : }
1942 : :
1943 : 830642 : void BoundPreservingLimiting( std::size_t nmat,
1944 : : std::size_t ndof,
1945 : : std::size_t e,
1946 : : const std::vector< std::size_t >& inpoel,
1947 : : const tk::UnsMesh::Coords& coord,
1948 : : const tk::Fields& U,
1949 : : std::vector< tk::real >& phic_p1,
1950 : : std::vector< tk::real >& phic_p2 )
1951 : : // *****************************************************************************
1952 : : // Bound preserving limiter for volume fractions when MulMat scheme is selected
1953 : : //! \param[in] nmat Number of materials in this PDE system
1954 : : //! \param[in] ndof Total number of reconstructed dofs
1955 : : //! \param[in] e Element being checked for consistency
1956 : : //! \param[in] inpoel Element connectivity
1957 : : //! \param[in] coord Array of nodal coordinates
1958 : : //! \param[in,out] U Second-order solution vector which gets modified near
1959 : : //! material interfaces for consistency
1960 : : //! \param[in] unk Vector of conservative variables based on Taylor basis
1961 : : //! \param[in,out] phic_p1 Vector of limiter functions for P1 dofs of the
1962 : : //! conserved quantities
1963 : : //! \param[in,out] phic_p2 Vector of limiter functions for P2 dofs of the
1964 : : //! conserved quantities
1965 : : //! \details This bound-preserving limiter is specifically meant to enforce
1966 : : //! bounds [0,1], but it does not suppress oscillations like the other 'TVD'
1967 : : //! limiters. TVD limiters on the other hand, do not preserve such bounds. A
1968 : : //! combination of oscillation-suppressing and bound-preserving limiters can
1969 : : //! obtain a non-oscillatory and bounded solution.
1970 : : // *****************************************************************************
1971 : : {
1972 : 830642 : const auto& cx = coord[0];
1973 : 830642 : const auto& cy = coord[1];
1974 : 830642 : const auto& cz = coord[2];
1975 : :
1976 : : // Extract the element coordinates
1977 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
1978 : 2491926 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
1979 : 2491926 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
1980 : 2491926 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
1981 : 7475778 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
1982 : :
1983 : : // Compute the determinant of Jacobian matrix
1984 : : auto detT =
1985 : 830642 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
1986 : :
1987 [ + - ]: 1661284 : std::vector< tk::real > phi_bound(nmat, 1.0);
1988 : :
1989 : : // Compute the upper and lower bound for volume fraction
1990 : 830642 : const tk::real min = 1e-14;
1991 : 830642 : const tk::real max = 1.0 - min * static_cast<tk::real>(nmat - 1);
1992 : :
1993 : : // loop over all faces of the element e
1994 [ + + ]: 4153210 : for (std::size_t lf=0; lf<4; ++lf)
1995 : : {
1996 : : // Extract the face coordinates
1997 : 3322568 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
1998 : 3322568 : inpoel[4*e+tk::lpofa[lf][1]],
1999 : 6645136 : inpoel[4*e+tk::lpofa[lf][2]] }};
2000 : :
2001 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
2002 : 9967704 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
2003 : 9967704 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
2004 : 19935408 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
2005 : :
2006 [ + - ]: 3322568 : auto ng = tk::NGfa(ndof);
2007 : :
2008 : : // arrays for quadrature points
2009 : 6645136 : std::array< std::vector< tk::real >, 2 > coordgp;
2010 : 6645136 : std::vector< tk::real > wgp;
2011 : :
2012 [ + - ]: 3322568 : coordgp[0].resize( ng );
2013 [ + - ]: 3322568 : coordgp[1].resize( ng );
2014 [ + - ]: 3322568 : wgp.resize( ng );
2015 : :
2016 : : // get quadrature point weights and coordinates for triangle
2017 [ + - ]: 3322568 : tk::GaussQuadratureTri( ng, coordgp, wgp );
2018 : :
2019 : : // Gaussian quadrature
2020 [ + + ]: 10561472 : for (std::size_t igp=0; igp<ng; ++igp)
2021 : : {
2022 : : // Compute the coordinates of quadrature point at physical domain
2023 [ + - ]: 7238904 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
2024 : :
2025 : : //Compute the basis functions
2026 : : auto B = tk::eval_basis( ndof,
2027 : 7238904 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
2028 : 7238904 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
2029 [ + - ]: 28955616 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
2030 : :
2031 [ + - ]: 14477808 : auto state = eval_state( U.nprop()/ndof, ndof, ndof, e, U, B );
2032 : :
2033 [ + + ]: 21716712 : for(std::size_t imat = 0; imat < nmat; imat++)
2034 : : {
2035 [ + - ]: 14477808 : auto phi = BoundPreservingLimitingFunction( min, max,
2036 : 14477808 : state[volfracIdx(nmat, imat)],
2037 [ + - ]: 14477808 : U(e,volfracDofIdx(nmat, imat, ndof, 0)) );
2038 : 14477808 : phi_bound[imat] = std::min( phi_bound[imat], phi );
2039 : : }
2040 : : }
2041 : : }
2042 : :
2043 : : // If DG(P2), the bound-preserving limiter should also be applied to the gauss
2044 : : // point within the element
2045 [ - + ]: 830642 : if(ndof > 4)
2046 : : {
2047 [ - - ]: 0 : auto ng = tk::NGvol(ndof);
2048 : :
2049 : : // arrays for quadrature points
2050 : 0 : std::array< std::vector< tk::real >, 3 > coordgp;
2051 : 0 : std::vector< tk::real > wgp;
2052 : :
2053 [ - - ]: 0 : coordgp[0].resize( ng );
2054 [ - - ]: 0 : coordgp[1].resize( ng );
2055 [ - - ]: 0 : coordgp[2].resize( ng );
2056 [ - - ]: 0 : wgp.resize( ng );
2057 : :
2058 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
2059 : :
2060 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
2061 : : {
2062 : : // Compute the basis function
2063 : 0 : auto B = tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp],
2064 [ - - ]: 0 : coordgp[2][igp] );
2065 : :
2066 [ - - ]: 0 : auto state = tk::eval_state(U.nprop()/ndof, ndof, ndof, e, U, B);
2067 : :
2068 [ - - ]: 0 : for(std::size_t imat = 0; imat < nmat; imat++)
2069 : : {
2070 [ - - ]: 0 : auto phi = BoundPreservingLimitingFunction(min, max,
2071 : 0 : state[volfracIdx(nmat, imat)],
2072 [ - - ]: 0 : U(e,volfracDofIdx(nmat, imat, ndof, 0)) );
2073 : 0 : phi_bound[imat] = std::min( phi_bound[imat], phi );
2074 : : }
2075 : : }
2076 : : }
2077 : :
2078 [ + + ]: 2491926 : for(std::size_t k=0; k<nmat; k++)
2079 : 3322568 : phic_p1[volfracIdx(nmat, k)] = std::min(phi_bound[k],
2080 : 3322568 : phic_p1[volfracIdx(nmat, k)]);
2081 : :
2082 [ - + ]: 830642 : if(ndof > 4)
2083 [ - - ]: 0 : for(std::size_t k=0; k<nmat; k++)
2084 : 0 : phic_p2[volfracIdx(nmat, k)] = std::min(phi_bound[k],
2085 : 0 : phic_p2[volfracIdx(nmat, k)]);
2086 : 830642 : }
2087 : :
2088 : : tk::real
2089 : 14477808 : BoundPreservingLimitingFunction( const tk::real min,
2090 : : const tk::real max,
2091 : : const tk::real al_gp,
2092 : : const tk::real al_avg )
2093 : : // *****************************************************************************
2094 : : // Bound-preserving limiter function for the volume fractions
2095 : : //! \param[in] min Minimum bound for volume fraction
2096 : : //! \param[in] max Maximum bound for volume fraction
2097 : : //! \param[in] al_gp Volume fraction at the quadrature point
2098 : : //! \param[in] al_avg Cell-average volume fraction
2099 : : //! \return The limiting coefficient from the bound-preserving limiter function
2100 : : // *****************************************************************************
2101 : : {
2102 : 14477808 : tk::real phi(1.0), al_diff(0.0);
2103 : 14477808 : al_diff = al_gp - al_avg;
2104 [ + + ][ + - ]: 14477808 : if(al_gp > max && fabs(al_diff) > 1e-15)
2105 : 269653 : phi = std::fabs( (max - al_avg) / al_diff );
2106 [ + + ][ + + ]: 14208155 : else if(al_gp < min && fabs(al_diff) > 1e-15)
2107 : 269653 : phi = std::fabs( (min - al_avg) / al_diff );
2108 : 14477808 : return phi;
2109 : : }
2110 : :
2111 : 1523580 : void PositivityLimiting( std::size_t nmat,
2112 : : std::size_t nspec,
2113 : : const std::vector< inciter::EOS >& mat_blk,
2114 : : std::size_t rdof,
2115 : : std::size_t ndof_el,
2116 : : const std::vector< std::size_t >& ndofel,
2117 : : std::size_t e,
2118 : : const std::vector< std::size_t >& inpoel,
2119 : : const tk::UnsMesh::Coords& coord,
2120 : : const std::vector< int >& esuel,
2121 : : const tk::Fields& U,
2122 : : const tk::Fields& P,
2123 : : std::vector< tk::real >& phic_p1,
2124 : : std::vector< tk::real >& phic_p2,
2125 : : std::vector< tk::real >& phip_p1,
2126 : : std::vector< tk::real >& phip_p2 )
2127 : : // *****************************************************************************
2128 : : // Positivity preserving limiter for multi-material and multspecies solver
2129 : : //! \param[in] nmat Number of materials in this PDE system
2130 : : //! \param[in] nspec Number of species in this PDE system
2131 : : //! \param[in] mat_blk EOS material block
2132 : : //! \param[in] rdof Total number of reconstructed dofs
2133 : : //! \param[in] ndof_el Number of dofs for element e
2134 : : //! \param[in] ndofel Vector of local number of degrees of freedome
2135 : : //! \param[in] e Element being checked for consistency
2136 : : //! \param[in] inpoel Element connectivity
2137 : : //! \param[in] coord Array of nodal coordinates
2138 : : //! \param[in] esuel Elements surrounding elements
2139 : : //! \param[in] U Vector of conservative variables
2140 : : //! \param[in] P Vector of primitive variables
2141 : : //! \param[in,out] phic_p1 Vector of limiter functions for P1 dofs of the
2142 : : //! conserved quantities
2143 : : //! \param[in,out] phic_p2 Vector of limiter functions for P2 dofs of the
2144 : : //! conserved quantities
2145 : : //! \param[in,out] phip_p1 Vector of limiter functions for P1 dofs of the
2146 : : //! primitive quantities
2147 : : //! \param[in,out] phip_p2 Vector of limiter functions for P2 dofs of the
2148 : : //! primitive quantities
2149 : : // *****************************************************************************
2150 : : {
2151 : 1523580 : const auto ncomp = U.nprop() / rdof;
2152 : 1523580 : const auto nprim = P.nprop() / rdof;
2153 : :
2154 : 1523580 : const auto& cx = coord[0];
2155 : 1523580 : const auto& cy = coord[1];
2156 : 1523580 : const auto& cz = coord[2];
2157 : :
2158 : : // Extract the element coordinates
2159 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
2160 : 4570740 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
2161 : 4570740 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
2162 : 4570740 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
2163 : 13712220 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
2164 : :
2165 : : // Compute the determinant of Jacobian matrix
2166 : : auto detT =
2167 : 1523580 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
2168 : :
2169 [ + - ]: 3047160 : std::vector< tk::real > phic_bound(ncomp, 1.0);
2170 [ + - ]: 3047160 : std::vector< tk::real > phip_bound(nprim, 1.0);
2171 : :
2172 [ + + ]: 7617900 : for (std::size_t lf=0; lf<4; ++lf)
2173 : : {
2174 : 6094320 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
2175 : 6094320 : inpoel[4*e+tk::lpofa[lf][1]],
2176 : 12188640 : inpoel[4*e+tk::lpofa[lf][2]] }};
2177 : :
2178 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
2179 : 18282960 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
2180 : 18282960 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
2181 : 36565920 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
2182 : :
2183 : : std::size_t nel;
2184 [ + + ]: 6094320 : if (esuel[4*e+lf] == -1) nel = e;
2185 : 5294340 : else nel = static_cast< std::size_t >(esuel[4*e+lf]);
2186 : :
2187 [ + - ]: 6094320 : auto ng = tk::NGfa(std::max(ndofel[e], ndofel[nel]));
2188 : :
2189 : 12188640 : std::array< std::vector< tk::real >, 2 > coordgp;
2190 : 12188640 : std::vector< tk::real > wgp;
2191 : :
2192 [ + - ]: 6094320 : coordgp[0].resize( ng );
2193 [ + - ]: 6094320 : coordgp[1].resize( ng );
2194 [ + - ]: 6094320 : wgp.resize( ng );
2195 : :
2196 [ + - ]: 6094320 : tk::GaussQuadratureTri( ng, coordgp, wgp );
2197 : :
2198 [ + + ]: 16190880 : for (std::size_t igp=0; igp<ng; ++igp)
2199 : : {
2200 [ + - ]: 10096560 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
2201 : : auto B = tk::eval_basis( ndof_el,
2202 : 10096560 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
2203 : 10096560 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
2204 [ + - ]: 40386240 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
2205 : :
2206 [ + - ]: 20193120 : auto state = eval_state(ncomp, rdof, ndof_el, e, U, B);
2207 [ + - ]: 20193120 : auto sprim = eval_state(nprim, rdof, ndof_el, e, P, B);
2208 : :
2209 [ + + ]: 10096560 : if (nmat > 1) {
2210 : : // multi-material PDE bounds
2211 [ + - ]: 7367760 : PositivityBoundsMultiMat(nmat, mat_blk, rdof, e, U, P, state, sprim,
2212 : : phic_bound, phip_bound);
2213 : : }
2214 : : else {
2215 : : // multispecies PDE bounds
2216 [ + - ]: 2728800 : PositivityBoundsMultiSpecies(nspec, mat_blk, rdof, e, U, P, state,
2217 : : sprim, phic_bound, phip_bound);
2218 : : }
2219 : : }
2220 : : }
2221 : :
2222 [ - + ]: 1523580 : if(ndofel[e] > 4)
2223 : : {
2224 [ - - ]: 0 : auto ng = tk::NGvol(ndof_el);
2225 : 0 : std::array< std::vector< tk::real >, 3 > coordgp;
2226 : 0 : std::vector< tk::real > wgp;
2227 : :
2228 [ - - ]: 0 : coordgp[0].resize( ng );
2229 [ - - ]: 0 : coordgp[1].resize( ng );
2230 [ - - ]: 0 : coordgp[2].resize( ng );
2231 [ - - ]: 0 : wgp.resize( ng );
2232 : :
2233 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
2234 : :
2235 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
2236 : : {
2237 : 0 : auto B = tk::eval_basis( ndof_el, coordgp[0][igp], coordgp[1][igp],
2238 [ - - ]: 0 : coordgp[2][igp] );
2239 : :
2240 [ - - ]: 0 : auto state = eval_state(ncomp, rdof, ndof_el, e, U, B);
2241 [ - - ]: 0 : auto sprim = eval_state(nprim, rdof, ndof_el, e, P, B);
2242 : :
2243 [ - - ]: 0 : if (nmat > 1) {
2244 : : // multi-material PDE bounds
2245 [ - - ]: 0 : PositivityBoundsMultiMat(nmat, mat_blk, rdof, e, U, P, state, sprim,
2246 : : phic_bound, phip_bound);
2247 : : }
2248 : : else {
2249 : : // multispecies PDE bounds
2250 [ - - ]: 0 : PositivityBoundsMultiSpecies(nspec, mat_blk, rdof, e, U, P, state,
2251 : : sprim, phic_bound, phip_bound);
2252 : : }
2253 : : }
2254 : : }
2255 : :
2256 : : // apply new bounds
2257 [ + + ]: 12848100 : for (std::size_t c=0; c<ncomp; ++c) {
2258 : 11324520 : phic_p1[c] = std::min( phic_bound[c], phic_p1[c] );
2259 [ - + ]: 11324520 : if (ndof_el > 4) phic_p2[c] = std::min( phic_bound[c], phic_p2[c] );
2260 : : }
2261 [ + + ]: 6412680 : for (std::size_t c=0; c<nprim; ++c) {
2262 : 4889100 : phip_p1[c] = std::min( phip_bound[c], phip_p1[c] );
2263 [ - + ]: 4889100 : if (ndof_el > 4) phip_p2[c] = std::min( phip_bound[c], phip_p2[c] );
2264 : : }
2265 : 1523580 : }
2266 : :
2267 : 7367760 : void PositivityBoundsMultiMat(
2268 : : std::size_t nmat,
2269 : : const std::vector< inciter::EOS >& mat_blk,
2270 : : std::size_t rdof,
2271 : : std::size_t e,
2272 : : const tk::Fields& U,
2273 : : const tk::Fields& P,
2274 : : const std::vector< tk::real >& state,
2275 : : const std::vector< tk::real >& sprim,
2276 : : std::vector< tk::real >& phic_bound,
2277 : : std::vector< tk::real >& phip_bound )
2278 : : // *****************************************************************************
2279 : : // Positivity bounds for multi-material PDE system
2280 : : //! \param[in] nmat Number of materials in this system
2281 : : //! \param[in] mat_blk EOS material block
2282 : : //! \param[in] rdof Total number of reconstructed dofs
2283 : : //! \param[in] e Element for which bounds are being calculated
2284 : : //! \param[in] U Vector of conservative variables
2285 : : //! \param[in] P Vector of primitive variables
2286 : : //! \param[in] state Vector of state of conserved quantities at quadrature point
2287 : : //! \param[in] sprim Vector of state of primitive quantities at quadrature point
2288 : : //! \param[in,out] phic_bound Vector of bounding-limiter functions for dofs of
2289 : : //! the conserved quantities
2290 : : //! \param[in,out] phip_bound Vector of bounding-limiter functions for dofs of
2291 : : //! the primitive quantities
2292 : : // *****************************************************************************
2293 : : {
2294 : 7367760 : const tk::real min = 1e-15;
2295 : :
2296 [ + + ]: 22103280 : for(std::size_t k = 0; k < nmat; k++)
2297 : : {
2298 : 14735520 : tk::real phi_rho(1.0), phi_rhoe(1.0), phi_pre(1.0);
2299 : : // Evaluate the limiting coefficient for material density
2300 : 14735520 : auto rho = state[densityIdx(nmat, k)];
2301 [ + - ]: 14735520 : auto rho_avg = U(e, densityDofIdx(nmat, k, rdof, 0));
2302 [ + - ]: 14735520 : phi_rho = PositivityFunction(min, rho, rho_avg);
2303 : 14735520 : phic_bound[densityIdx(nmat, k)] =
2304 : 14735520 : std::min(phic_bound[densityIdx(nmat, k)], phi_rho);
2305 : : // Evaluate the limiting coefficient for material energy
2306 : 14735520 : auto rhoe = state[energyIdx(nmat, k)];
2307 [ + - ]: 14735520 : auto rhoe_avg = U(e, energyDofIdx(nmat, k, rdof, 0));
2308 [ + - ]: 14735520 : phi_rhoe = PositivityFunction(min, rhoe, rhoe_avg);
2309 : 14735520 : phic_bound[energyIdx(nmat, k)] =
2310 : 14735520 : std::min(phic_bound[energyIdx(nmat, k)], phi_rhoe);
2311 : : // Evaluate the limiting coefficient for material pressure
2312 : 14735520 : auto min_pre = std::max(min, state[volfracIdx(nmat, k)] *
2313 [ + - ]: 29471040 : mat_blk[k].compute< EOS::min_eff_pressure >(min, rho,
2314 : 29471040 : state[volfracIdx(nmat, k)]));
2315 : 14735520 : auto pre = sprim[pressureIdx(nmat, k)];
2316 [ + - ]: 14735520 : auto pre_avg = P(e, pressureDofIdx(nmat, k, rdof, 0));
2317 [ + - ]: 14735520 : phi_pre = PositivityFunction(min_pre, pre, pre_avg);
2318 : 14735520 : phip_bound[pressureIdx(nmat, k)] =
2319 : 14735520 : std::min(phip_bound[pressureIdx(nmat, k)], phi_pre);
2320 : : }
2321 : 7367760 : }
2322 : :
2323 : 2728800 : void PositivityBoundsMultiSpecies(
2324 : : std::size_t nspec,
2325 : : const std::vector< inciter::EOS >& /*mat_blk*/,
2326 : : std::size_t rdof,
2327 : : std::size_t e,
2328 : : const tk::Fields& /*U*/,
2329 : : const tk::Fields& P,
2330 : : const std::vector< tk::real >& /*state*/,
2331 : : const std::vector< tk::real >& sprim,
2332 : : std::vector< tk::real >& /*phic_bound*/,
2333 : : std::vector< tk::real >& phip_bound )
2334 : : // *****************************************************************************
2335 : : // Positivity bounds for multispecies PDE system
2336 : : //! \param[in] nmat Number of materials in this system
2337 : : // //! \param[in] mat_blk EOS material block
2338 : : //! \param[in] rdof Total number of reconstructed dofs
2339 : : //! \param[in] e Element for which bounds are being calculated
2340 : : // //! \param[in] U Vector of conservative variables
2341 : : //! \param[in] P Vector of primitive variables
2342 : : // //! \param[in] state Vector of state of conserved quantities at quadrature point
2343 : : //! \param[in] sprim Vector of state of primitive quantities at quadrature point
2344 : : // //! \param[in,out] phic_bound Vector of bounding-limiter functions for dofs of
2345 : : // //! the conserved quantities
2346 : : //! \param[in,out] phip_bound Vector of bounding-limiter functions for dofs of
2347 : : //! the primitive quantities
2348 : : // *****************************************************************************
2349 : : {
2350 : 2728800 : const tk::real min = 1e-8;
2351 : :
2352 : 2728800 : tk::real phi_T(1.0);
2353 : : // Evaluate the limiting coefficient for mixture temperature
2354 : 2728800 : auto T = sprim[multispecies::temperatureIdx(nspec, 0)];
2355 [ + - ]: 2728800 : auto T_avg = P(e, multispecies::temperatureDofIdx(nspec, 0, rdof, 0));
2356 [ + - ]: 2728800 : phi_T = PositivityFunction(min, T, T_avg);
2357 : 2728800 : phip_bound[multispecies::temperatureIdx(nspec, 0)] =
2358 : 2728800 : std::min(phip_bound[multispecies::temperatureIdx(nspec, 0)], phi_T);
2359 : 2728800 : }
2360 : :
2361 : 1368 : void PositivityPreservingMultiMat_FV(
2362 : : const std::vector< std::size_t >& inpoel,
2363 : : std::size_t nelem,
2364 : : std::size_t nmat,
2365 : : const std::vector< inciter::EOS >& mat_blk,
2366 : : const tk::UnsMesh::Coords& coord,
2367 : : const tk::Fields& /*geoFace*/,
2368 : : tk::Fields& U,
2369 : : tk::Fields& P )
2370 : : // *****************************************************************************
2371 : : // Positivity preserving limiter for the FV multi-material solver
2372 : : //! \param[in] inpoel Element connectivity
2373 : : //! \param[in] nelem Number of elements
2374 : : //! \param[in] nmat Number of materials in this PDE system
2375 : : //! \param[in] mat_blk Material EOS block
2376 : : //! \param[in] coord Array of nodal coordinates
2377 : : ////! \param[in] geoFace Face geometry array
2378 : : //! \param[in,out] U High-order solution vector which gets limited
2379 : : //! \param[in,out] P High-order vector of primitives which gets limited
2380 : : //! \details This positivity preserving limiter function should be called for
2381 : : //! FV multimat.
2382 : : // *****************************************************************************
2383 : : {
2384 : 1368 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
2385 : 1368 : const auto ncomp = U.nprop() / rdof;
2386 : 1368 : const auto nprim = P.nprop() / rdof;
2387 : :
2388 : 1368 : const auto& cx = coord[0];
2389 : 1368 : const auto& cy = coord[1];
2390 : 1368 : const auto& cz = coord[2];
2391 : :
2392 [ + + ]: 266168 : for (std::size_t e=0; e<nelem; ++e)
2393 : : {
2394 : : // Extract the element coordinates
2395 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
2396 : 794400 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
2397 : 794400 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
2398 : 794400 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
2399 : 2383200 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
2400 : :
2401 : : // Compute the determinant of Jacobian matrix
2402 : : auto detT =
2403 : 264800 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
2404 : :
2405 [ + - ]: 529600 : std::vector< tk::real > phic(ncomp, 1.0);
2406 [ + - ]: 529600 : std::vector< tk::real > phip(nprim, 1.0);
2407 : :
2408 : 264800 : const tk::real min = 1e-15;
2409 : :
2410 : : // 1. Enforce positive density (total energy will be positive if pressure
2411 : : // and density are positive)
2412 [ + + ]: 1324000 : for (std::size_t lf=0; lf<4; ++lf)
2413 : : {
2414 : 1059200 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
2415 : 1059200 : inpoel[4*e+tk::lpofa[lf][1]],
2416 : 2118400 : inpoel[4*e+tk::lpofa[lf][2]] }};
2417 : :
2418 : : // face coordinates
2419 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
2420 : 3177600 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
2421 : 3177600 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
2422 : 6355200 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
2423 : :
2424 : : // face centroid
2425 : : std::array< tk::real, 3 > fc{{
2426 : 1059200 : (coordfa[0][0]+coordfa[1][0]+coordfa[2][0])/3.0 ,
2427 : 1059200 : (coordfa[0][1]+coordfa[1][1]+coordfa[2][1])/3.0 ,
2428 : 2118400 : (coordfa[0][2]+coordfa[1][2]+coordfa[2][2])/3.0 }};
2429 : :
2430 : : auto B = tk::eval_basis( rdof,
2431 : 1059200 : tk::Jacobian( coordel[0], fc, coordel[2], coordel[3] ) / detT,
2432 : 1059200 : tk::Jacobian( coordel[0], coordel[1], fc, coordel[3] ) / detT,
2433 [ + - ]: 4236800 : tk::Jacobian( coordel[0], coordel[1], coordel[2], fc ) / detT );
2434 [ + - ]: 2118400 : auto state = eval_state(ncomp, rdof, rdof, e, U, B);
2435 : :
2436 [ + + ]: 3387840 : for(std::size_t i=0; i<nmat; i++)
2437 : : {
2438 : : // Evaluate the limiting coefficient for material density
2439 : 2328640 : auto rho = state[densityIdx(nmat, i)];
2440 [ + - ]: 2328640 : auto rho_avg = U(e, densityDofIdx(nmat, i, rdof, 0));
2441 [ + - ]: 2328640 : auto phi_rho = PositivityFunction(min, rho, rho_avg);
2442 : 2328640 : phic[densityIdx(nmat, i)] =
2443 : 2328640 : std::min(phic[densityIdx(nmat, i)], phi_rho);
2444 : : }
2445 : : }
2446 : : // apply limiter coefficient
2447 [ + + ]: 846960 : for(std::size_t i=0; i<nmat; i++)
2448 : : {
2449 [ + - ]: 582160 : U(e, densityDofIdx(nmat,i,rdof,1)) *= phic[densityIdx(nmat,i)];
2450 [ + - ]: 582160 : U(e, densityDofIdx(nmat,i,rdof,2)) *= phic[densityIdx(nmat,i)];
2451 [ + - ]: 582160 : U(e, densityDofIdx(nmat,i,rdof,3)) *= phic[densityIdx(nmat,i)];
2452 : : }
2453 : :
2454 : : // 2. Enforce positive pressure (assuming density is positive)
2455 [ + + ]: 1324000 : for (std::size_t lf=0; lf<4; ++lf)
2456 : : {
2457 : 1059200 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
2458 : 1059200 : inpoel[4*e+tk::lpofa[lf][1]],
2459 : 2118400 : inpoel[4*e+tk::lpofa[lf][2]] }};
2460 : :
2461 : : // face coordinates
2462 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
2463 : 3177600 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
2464 : 3177600 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
2465 : 6355200 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
2466 : :
2467 : : // face centroid
2468 : : std::array< tk::real, 3 > fc{{
2469 : 1059200 : (coordfa[0][0]+coordfa[1][0]+coordfa[2][0])/3.0 ,
2470 : 1059200 : (coordfa[0][1]+coordfa[1][1]+coordfa[2][1])/3.0 ,
2471 : 2118400 : (coordfa[0][2]+coordfa[1][2]+coordfa[2][2])/3.0 }};
2472 : :
2473 : : auto B = tk::eval_basis( rdof,
2474 : 1059200 : tk::Jacobian( coordel[0], fc, coordel[2], coordel[3] ) / detT,
2475 : 1059200 : tk::Jacobian( coordel[0], coordel[1], fc, coordel[3] ) / detT,
2476 [ + - ]: 4236800 : tk::Jacobian( coordel[0], coordel[1], coordel[2], fc ) / detT );
2477 [ + - ]: 2118400 : auto state = eval_state(ncomp, rdof, rdof, e, U, B);
2478 [ + - ]: 2118400 : auto sprim = eval_state(nprim, rdof, rdof, e, P, B);
2479 : :
2480 [ + + ]: 3387840 : for(std::size_t i=0; i<nmat; i++)
2481 : : {
2482 : 2328640 : tk::real phi_pre(1.0);
2483 : : // Evaluate the limiting coefficient for material pressure
2484 : 2328640 : auto rho = state[densityIdx(nmat, i)];
2485 [ + - ]: 2328640 : auto min_pre = std::max(min, U(e,volfracDofIdx(nmat,i,rdof,0)) *
2486 : 2328640 : mat_blk[i].compute< EOS::min_eff_pressure >(min, rho,
2487 [ + - ][ + - ]: 2328640 : U(e,volfracDofIdx(nmat,i,rdof,0))));
2488 : 2328640 : auto pre = sprim[pressureIdx(nmat, i)];
2489 [ + - ]: 2328640 : auto pre_avg = P(e, pressureDofIdx(nmat, i, rdof, 0));
2490 [ + - ]: 2328640 : phi_pre = PositivityFunction(min_pre, pre, pre_avg);
2491 : 2328640 : phip[pressureIdx(nmat, i)] =
2492 : 2328640 : std::min(phip[pressureIdx(nmat, i)], phi_pre);
2493 : : }
2494 : : }
2495 : : // apply limiter coefficient
2496 [ + + ]: 846960 : for(std::size_t i=0; i<nmat; i++)
2497 : : {
2498 [ + - ]: 582160 : P(e, pressureDofIdx(nmat,i,rdof,1)) *= phip[pressureIdx(nmat,i)];
2499 [ + - ]: 582160 : P(e, pressureDofIdx(nmat,i,rdof,2)) *= phip[pressureIdx(nmat,i)];
2500 [ + - ]: 582160 : P(e, pressureDofIdx(nmat,i,rdof,3)) *= phip[pressureIdx(nmat,i)];
2501 : : }
2502 : : }
2503 : 1368 : }
2504 : :
2505 : : tk::real
2506 : 51592640 : PositivityFunction( const tk::real min,
2507 : : const tk::real u_gp,
2508 : : const tk::real u_avg )
2509 : : // *****************************************************************************
2510 : : // Positivity-preserving limiter function
2511 : : //! \param[in] min Minimum bound for volume fraction
2512 : : //! \param[in] u_gp Variable quantity at the quadrature point
2513 : : //! \param[in] u_avg Cell-average variable quantitiy
2514 : : //! \return The limiting coefficient from the positivity-preserving limiter
2515 : : //! function
2516 : : // *****************************************************************************
2517 : : {
2518 : 51592640 : tk::real phi(1.0);
2519 : 51592640 : tk::real diff = u_gp - u_avg;
2520 : : // Only when u_gp is less than minimum threshold and the high order
2521 : : // contribution is not zero, the limiting function will be applied
2522 [ + + ]: 51592640 : if(u_gp < min)
2523 : 322004 : phi = std::fabs( (min - u_avg) / (diff+std::copysign(1e-15,diff)) );
2524 : 51592640 : return phi;
2525 : : }
2526 : :
2527 : : bool
2528 : 30363106 : interfaceIndicator( std::size_t nmat,
2529 : : const std::vector< tk::real >& al,
2530 : : std::vector< std::size_t >& matInt )
2531 : : // *****************************************************************************
2532 : : // Interface indicator function, which checks element for material interface
2533 : : //! \param[in] nmat Number of materials in this PDE system
2534 : : //! \param[in] al Cell-averaged volume fractions
2535 : : //! \param[in] matInt Array indicating which material has an interface
2536 : : //! \return Boolean which indicates if the element contains a material interface
2537 : : // *****************************************************************************
2538 : : {
2539 : : auto alphamin =
2540 : 30363106 : inciter::g_inputdeck.get< tag::multimat, tag::min_volumefrac >();
2541 : :
2542 : 30363106 : bool intInd = false;
2543 : :
2544 : : // limits under which compression is to be performed
2545 : 30363106 : auto al_eps = std::min(1e-08, 1e4*alphamin); // limit this value at 1e-8
2546 : 30363106 : auto loLim = 2.0 * al_eps;
2547 : 30363106 : auto hiLim = 1.0 - loLim;
2548 : :
2549 : 30363106 : auto almax = 0.0;
2550 [ + + ]: 96928694 : for (std::size_t k=0; k<nmat; ++k)
2551 : : {
2552 : 66565588 : almax = std::max(almax, al[k]);
2553 : 66565588 : matInt[k] = 0;
2554 [ + + ][ + + ]: 66565588 : if ((al[k] > loLim) && (al[k] < hiLim)) matInt[k] = 1;
[ + + ]
2555 : : }
2556 : :
2557 [ + - ][ + + ]: 30363106 : if ((almax > loLim) && (almax < hiLim)) intInd = true;
2558 : :
2559 : 30363106 : return intInd;
2560 : : }
2561 : :
2562 : 3330 : void MarkShockCells ( const bool pref,
2563 : : const std::size_t nelem,
2564 : : const std::size_t nmat,
2565 : : const std::size_t ndof,
2566 : : const std::size_t rdof,
2567 : : const std::vector< inciter::EOS >& mat_blk,
2568 : : const std::vector< std::size_t >& ndofel,
2569 : : const std::vector< std::size_t >& inpoel,
2570 : : const tk::UnsMesh::Coords& coord,
2571 : : const inciter::FaceData& fd,
2572 : : [[maybe_unused]] const tk::Fields& geoFace,
2573 : : const tk::Fields& geoElem,
2574 : : const tk::FluxFn& flux,
2575 : : const std::vector< std::size_t >& solidx,
2576 : : const tk::Fields& U,
2577 : : const tk::Fields& P,
2578 : : const std::set< std::size_t >& vars,
2579 : : std::vector< std::size_t >& shockmarker )
2580 : : // *****************************************************************************
2581 : : // Mark the cells that contain discontinuity according to the interface
2582 : : // condition
2583 : : //! \param[in] pref Indicator for p-adaptive algorithm
2584 : : //! \param[in] nelem Number of elements
2585 : : //! \param[in] nmat Number of materials in this PDE system
2586 : : //! \param[in] ndof Maximum number of degrees of freedom
2587 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
2588 : : //! \param[in] mat_blk EOS material block
2589 : : //! \param[in] ndofel Vector of local number of degrees of freedome
2590 : : //! \param[in] inpoel Element-node connectivity
2591 : : //! \param[in] coord Array of nodal coordinates
2592 : : //! \param[in] fd Face connectivity and boundary conditions object
2593 : : //! \param[in] geoFace Face geometry array
2594 : : //! \param[in] geoElem Element geometry array
2595 : : //! \param[in] flux Flux function to use
2596 : : //! \param[in] solidx Solid material index indicator
2597 : : //! \param[in] U Solution vector at recent time step
2598 : : //! \param[in] P Vector of primitives at recent time step
2599 : : //! \param[in] vars Vector of variable indices to evaluate flux jump
2600 : : //! \param[in, out] shockmarker Vector of the shock indicator
2601 : : //! \details This function computes the discontinuity indicator based on
2602 : : //! interface conditon. It is based on the following paper:
2603 : : //! Hong L., Gianni A., Robert N. (2021) A moving discontinuous Galerkin
2604 : : //! finite element method with interface condition enforcement for
2605 : : //! compressible flows. Journal of Computational Physics,
2606 : : //! doi: https://doi.org/10.1016/j.jcp.2021.110618
2607 : : // *****************************************************************************
2608 : : {
2609 : 3330 : const auto coeff = g_inputdeck.get< tag::shock_detector_coeff >();
2610 : :
2611 [ + - ]: 6660 : std::vector< tk::real > IC(U.nunk(), 0.0);
2612 : 3330 : const auto& esuf = fd.Esuf();
2613 : 3330 : const auto& inpofa = fd.Inpofa();
2614 : :
2615 : 3330 : const auto& cx = coord[0];
2616 : 3330 : const auto& cy = coord[1];
2617 : 3330 : const auto& cz = coord[2];
2618 : :
2619 : 3330 : auto ncomp = U.nprop()/rdof;
2620 : 3330 : auto nprim = P.nprop()/rdof;
2621 : :
2622 : : // Loop over faces
2623 [ + - ][ + + ]: 1558710 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f) {
2624 [ + - ][ + - ]: 1555380 : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
[ - - ][ - - ]
[ - - ]
2625 : : "as -1" );
2626 : :
2627 : 1555380 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
2628 : 1555380 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
2629 : :
2630 : : // When the number of gauss points for the left and right element are
2631 : : // different, choose the larger ng
2632 [ + - ]: 1555380 : auto ng_l = tk::NGfa(ndofel[el]);
2633 [ + - ]: 1555380 : auto ng_r = tk::NGfa(ndofel[er]);
2634 : :
2635 : 1555380 : auto ng = std::max( ng_l, ng_r );
2636 : :
2637 : : std::array< std::vector< tk::real >, 2 > coordgp
2638 [ + - ][ + - ]: 3110760 : { std::vector<tk::real>(ng), std::vector<tk::real>(ng) };
2639 [ + - ]: 3110760 : std::vector< tk::real > wgp( ng );
2640 : :
2641 [ + - ]: 1555380 : tk::GaussQuadratureTri( ng, coordgp, wgp );
2642 : :
2643 : : // Extract the element coordinates
2644 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
2645 : 4666140 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
2646 : 4666140 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
2647 : 4666140 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
2648 : 13998420 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
2649 : :
2650 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
2651 : 4666140 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
2652 : 4666140 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
2653 : 4666140 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
2654 : 13998420 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
2655 : :
2656 : : // Compute the determinant of Jacobian matrix
2657 : : auto detT_l =
2658 : 1555380 : tk::Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
2659 : : auto detT_r =
2660 : 1555380 : tk::Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
2661 : :
2662 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
2663 : 4666140 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
2664 : 4666140 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
2665 : 9332280 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
2666 : :
2667 : : std::array< tk::real, 3 >
2668 [ + - ][ + - ]: 1555380 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ + - ]
2669 : :
2670 : : // Numerator and denominator of the shock indicator
2671 : 1555380 : tk::real numer(0.0), denom(0.0);
2672 : 3110760 : std::vector< tk::real > fl_jump, fl_avg;
2673 [ + - ]: 1555380 : fl_jump.resize(3, 0.0);
2674 [ + - ]: 1555380 : fl_avg.resize(3, 0.0);
2675 : :
2676 [ + + ]: 5394186 : for (std::size_t igp=0; igp<ng; ++igp) {
2677 [ + - ]: 3838806 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
2678 : : std::size_t dof_el, dof_er;
2679 [ - + ]: 3838806 : if (rdof > ndof)
2680 : : {
2681 : 0 : dof_el = rdof;
2682 : 0 : dof_er = rdof;
2683 : : }
2684 : : else
2685 : : {
2686 : 3838806 : dof_el = ndofel[el];
2687 : 3838806 : dof_er = ndofel[er];
2688 : : }
2689 : : std::array< tk::real, 3> ref_gp_l{
2690 : 3838806 : tk::Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
2691 : 3838806 : tk::Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
2692 : 7677612 : tk::Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
2693 : : std::array< tk::real, 3> ref_gp_r{
2694 : 3838806 : tk::Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
2695 : 3838806 : tk::Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
2696 : 7677612 : tk::Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
2697 [ + - ]: 7677612 : auto B_l = tk::eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
2698 [ + - ]: 7677612 : auto B_r = tk::eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
2699 : :
2700 : : // Evaluate the high order solution at the qudrature point
2701 : 7677612 : std::array< std::vector< tk::real >, 2 > state;
2702 [ + - ]: 7677612 : state[0] = tk::evalPolynomialSol(mat_blk, 0, ncomp, nprim, rdof,
2703 : 3838806 : nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
2704 [ + - ]: 7677612 : state[1] = tk::evalPolynomialSol(mat_blk, 0, ncomp, nprim, rdof,
2705 : 3838806 : nmat, er, dof_er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P);
2706 : :
2707 [ - + ][ - - ]: 3838806 : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
2708 : : "appended boundary state vector" );
2709 [ - + ][ - - ]: 3838806 : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
2710 : : "appended boundary state vector" );
2711 : :
2712 : : // Force deformation unknown to first order
2713 [ + + ]: 10113912 : for (std::size_t k=0; k<nmat; ++k)
2714 [ - + ]: 6275106 : if (solidx[k] > 0)
2715 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
2716 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
2717 : : {
2718 : 0 : state[0][deformIdx(nmat, solidx[k], i, j)] = U(el,deformDofIdx(
2719 [ - - ]: 0 : nmat, solidx[k], i, j, rdof, 0));
2720 : 0 : state[1][deformIdx(nmat, solidx[k], i, j)] = U(er,deformDofIdx(
2721 [ - - ]: 0 : nmat, solidx[k], i, j, rdof, 0));
2722 : : }
2723 : :
2724 : : // Evaluate the flux
2725 [ + - ]: 7677612 : auto fl = flux( ncomp, mat_blk, state[0], {} );
2726 [ + - ]: 7677612 : auto fr = flux( ncomp, mat_blk, state[1], {} );
2727 : :
2728 : 3838806 : std::size_t i(0);
2729 [ + + ]: 15355224 : for (const auto& c : vars) {
2730 : 11516418 : tk::real fn_l(0.0), fn_r(0.0);
2731 [ + + ]: 46065672 : for(std::size_t idir = 0; idir < 3; idir++) {
2732 : 34549254 : fn_l += fl[c][idir] * fn[idir];
2733 : 34549254 : fn_r += fr[c][idir] * fn[idir];
2734 : : }
2735 : 11516418 : fl_jump[i] += wgp[igp] * (fn_l - fn_r) * (fn_l - fn_r);
2736 : 11516418 : fl_avg[i] += wgp[igp] * (fn_l + fn_r) * (fn_l + fn_r) * 0.25;
2737 : 11516418 : ++i;
2738 : : }
2739 : : }
2740 : :
2741 : : // Evaluate the numerator and denominator
2742 [ + + ]: 6221520 : for(std::size_t idir = 0; idir < 3; idir++) {
2743 : 4666140 : numer += std::sqrt(fl_jump[idir]);
2744 : 4666140 : denom += std::sqrt(fl_avg[idir]);
2745 : : }
2746 : :
2747 : 1555380 : tk::real Ind(0.0);
2748 [ + - ]: 1555380 : if(denom > 1e-8)
2749 : 1555380 : Ind = numer / denom;
2750 : 1555380 : IC[el] = std::max(IC[el], Ind);
2751 : 1555380 : IC[er] = std::max(IC[er], Ind);
2752 : : }
2753 : :
2754 : : // Loop over element to mark shock cell
2755 [ + + ]: 821970 : for (std::size_t e=0; e<nelem; ++e) {
2756 [ - + ]: 818640 : std::size_t dof_el = pref ? ndofel[e] : rdof;
2757 : :
2758 : 818640 : tk::real power = 0.0;
2759 [ + + ]: 818640 : if(dof_el == 10) power = 1.5;
2760 : 454800 : else power = 1.0;
2761 : :
2762 : : // Evaluate the threshold
2763 [ + - ]: 818640 : auto thres = coeff * std::pow(geoElem(e, 4), power);
2764 [ + + ]: 818640 : if(IC[e] > thres)
2765 : 50304 : shockmarker[e] = 1;
2766 : : else
2767 : 768336 : shockmarker[e] = 0;
2768 : : }
2769 : 3330 : }
2770 : :
2771 : : void
2772 : 30 : correctLimConservMultiMat(
2773 : : std::size_t nelem,
2774 : : const std::vector< EOS >& mat_blk,
2775 : : std::size_t nmat,
2776 : : const std::vector< std::size_t >& inpoel,
2777 : : const tk::UnsMesh::Coords& coord,
2778 : : const tk::Fields& geoElem,
2779 : : const tk::Fields& prim,
2780 : : tk::Fields& unk )
2781 : : // *****************************************************************************
2782 : : // Update the conservative quantities after limiting for multi-material systems
2783 : : //! \param[in] nelem Number of internal elements
2784 : : //! \param[in] mat_blk EOS material block
2785 : : //! \param[in] nmat Number of materials in this PDE system
2786 : : //! \param[in] inpoel Element-node connectivity
2787 : : //! \param[in] coord Array of nodal coordinates
2788 : : //! \param[in] geoElem Element geometry array
2789 : : //! \param[in] prim Array of primitive variables
2790 : : //! \param[in,out] unk Array of conservative variables
2791 : : //! \details This function computes the updated dofs for conservative
2792 : : //! quantities based on the limited primitive quantities, to re-instate
2793 : : //! consistency between the limited primitive and evolved quantities. For
2794 : : //! further details, see Pandare et al. (2023). On the Design of Stable,
2795 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
2796 : : //! Hydrodynamics. J Comp Phys, 112313.
2797 : : // *****************************************************************************
2798 : : {
2799 : 30 : const auto rdof = g_inputdeck.get< tag::rdof >();
2800 : 30 : std::size_t ncomp = unk.nprop()/rdof;
2801 : 30 : std::size_t nprim = prim.nprop()/rdof;
2802 : : const auto intsharp = inciter::g_inputdeck.get< tag::multimat,
2803 : 30 : tag::intsharp >();
2804 : :
2805 : 30 : auto L = tk::massMatrixDubiner();
2806 : :
2807 [ + + ]: 45510 : for (std::size_t e=0; e<nelem; ++e) {
2808 [ + - ]: 45480 : auto vole = geoElem(e,0);
2809 : :
2810 : : // The right-hand side vector is sized as nprim, i.e. the primitive quantity
2811 : : // vector. However, it stores the consistently obtained values of evolved
2812 : : // quantities, since nprim is the number of evolved quantities that need to
2813 : : // be evaluated consistently. For this reason, accessing R will require
2814 : : // the primitive quantity accessors. But this access is intended to give
2815 : : // the corresponding evolved quantites, as follows:
2816 : : // pressureIdx() - mat. total energy
2817 : : // velocityIdx() - bulk momentum components
2818 : : // stressIdx() - mat. inverse deformation gradient tensor components
2819 [ + - ]: 90960 : std::vector< tk::real > R(nprim*rdof, 0.0);
2820 : :
2821 [ + - ]: 45480 : auto ng = tk::NGvol(rdof);
2822 : :
2823 : : // Arrays for quadrature points
2824 : 90960 : std::array< std::vector< tk::real >, 3 > coordgp;
2825 : 90960 : std::vector< tk::real > wgp;
2826 : :
2827 [ + - ]: 45480 : coordgp[0].resize( ng );
2828 [ + - ]: 45480 : coordgp[1].resize( ng );
2829 [ + - ]: 45480 : coordgp[2].resize( ng );
2830 [ + - ]: 45480 : wgp.resize( ng );
2831 : :
2832 [ + - ]: 45480 : tk::GaussQuadratureTet( ng, coordgp, wgp );
2833 : :
2834 : : // Loop over quadrature points in element e
2835 [ + + ]: 272880 : for (std::size_t igp=0; igp<ng; ++igp) {
2836 : : // Compute the basis function
2837 : 454800 : auto B = tk::eval_basis( rdof, coordgp[0][igp], coordgp[1][igp],
2838 [ + - ]: 909600 : coordgp[2][igp] );
2839 : :
2840 : 227400 : auto w = wgp[igp] * vole;
2841 : :
2842 : : // Evaluate the solution at quadrature point
2843 : : auto state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
2844 : : rdof, nmat, e, rdof, inpoel, coord, geoElem,
2845 [ + - ]: 454800 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, unk, prim);
2846 : :
2847 : : // Solution vector that stores the material energy and bulk momentum
2848 [ + - ]: 454800 : std::vector< tk::real > s(nprim, 0.0);
2849 : :
2850 : : // Bulk density at quadrature point
2851 : 227400 : tk::real rhob(0.0);
2852 [ + + ]: 682200 : for (std::size_t k=0; k<nmat; ++k)
2853 : 454800 : rhob += state[densityIdx(nmat, k)];
2854 : :
2855 : : // Velocity vector at quadrature point
2856 : : std::array< tk::real, 3 >
2857 : 227400 : vel{ state[ncomp+velocityIdx(nmat, 0)],
2858 : 227400 : state[ncomp+velocityIdx(nmat, 1)],
2859 : 454800 : state[ncomp+velocityIdx(nmat, 2)] };
2860 : :
2861 : : // Compute and store the bulk momentum
2862 [ + + ]: 909600 : for(std::size_t idir = 0; idir < 3; idir++)
2863 : 682200 : s[velocityIdx(nmat, idir)] = rhob * vel[idir];
2864 : :
2865 : : // Compute and store material energy at quadrature point
2866 [ + + ]: 682200 : for(std::size_t imat = 0; imat < nmat; imat++) {
2867 : 454800 : auto alphamat = state[volfracIdx(nmat, imat)];
2868 : 454800 : auto arhomat = state[densityIdx(nmat, imat)];
2869 : 454800 : auto apremat = state[ncomp+pressureIdx(nmat, imat)];
2870 [ + - ]: 454800 : auto gmat = getDeformGrad(nmat, imat, state);
2871 : 454800 : s[pressureIdx(nmat,imat)] =
2872 [ + - ]: 909600 : mat_blk[imat].compute< EOS::totalenergy >( arhomat, vel[0], vel[1],
2873 : 454800 : vel[2], apremat, alphamat, gmat );
2874 : : }
2875 : :
2876 : : // Evaluate the righ-hand-side vector
2877 [ + + ]: 1364400 : for(std::size_t k = 0; k < nprim; k++) {
2878 : 1137000 : auto mark = k * rdof;
2879 [ + + ]: 5685000 : for(std::size_t idof = 0; idof < rdof; idof++)
2880 : 4548000 : R[mark+idof] += w * s[k] * B[idof];
2881 : : }
2882 : : }
2883 : :
2884 : : // Update the high order dofs of the material energy
2885 [ + + ]: 136440 : for(std::size_t imat = 0; imat < nmat; imat++) {
2886 [ + + ]: 363840 : for(std::size_t idof = 1; idof < rdof; idof++)
2887 [ + - ]: 272880 : unk(e, energyDofIdx(nmat, imat, rdof, idof)) =
2888 : 272880 : R[pressureDofIdx(nmat,imat,rdof,idof)] / (vole*L[idof]);
2889 : : }
2890 : :
2891 : : // Update the high order dofs of the bulk momentum
2892 [ + + ]: 181920 : for(std::size_t idir = 0; idir < 3; idir++) {
2893 [ + + ]: 545760 : for(std::size_t idof = 1; idof < rdof; idof++)
2894 [ + - ]: 409320 : unk(e, momentumDofIdx(nmat, idir, rdof, idof)) =
2895 : 409320 : R[velocityDofIdx(nmat,idir,rdof,idof)] / (vole*L[idof]);
2896 : : }
2897 : : }
2898 : 30 : }
2899 : :
2900 : : void
2901 : 6600 : correctLimConservMultiSpecies(
2902 : : std::size_t nelem,
2903 : : const std::vector< EOS >& mat_blk,
2904 : : std::size_t nspec,
2905 : : const std::vector< std::size_t >& inpoel,
2906 : : const tk::UnsMesh::Coords& coord,
2907 : : const tk::Fields& geoElem,
2908 : : const tk::Fields& prim,
2909 : : tk::Fields& unk )
2910 : : // *****************************************************************************
2911 : : // Update the conservative quantities after limiting for multispecies systems
2912 : : //! \param[in] nelem Number of internal elements
2913 : : //! \param[in] mat_blk EOS material block
2914 : : //! \param[in] nspec Number of species in this PDE system
2915 : : //! \param[in] inpoel Element-node connectivity
2916 : : //! \param[in] coord Array of nodal coordinates
2917 : : //! \param[in] geoElem Element geometry array
2918 : : //! \param[in] prim Array of primitive variables
2919 : : //! \param[in,out] unk Array of conservative variables
2920 : : //! \details This function computes the updated dofs for conservative
2921 : : //! quantities based on the limited primitive quantities, to re-instate
2922 : : //! consistency between the limited primitive and evolved quantities. For
2923 : : //! further details, see Pandare et al. (2023). On the Design of Stable,
2924 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
2925 : : //! Hydrodynamics. J Comp Phys, 112313.
2926 : : // *****************************************************************************
2927 : : {
2928 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
2929 : 6600 : std::size_t ncomp = unk.nprop()/rdof;
2930 : 6600 : std::size_t nprim = prim.nprop()/rdof;
2931 : :
2932 : 6600 : auto L = tk::massMatrixDubiner();
2933 : :
2934 [ + + ]: 688800 : for (std::size_t e=0; e<nelem; ++e) {
2935 [ + - ]: 682200 : auto vole = geoElem(e,0);
2936 : :
2937 : : // The right-hand side vector is sized as nprim, i.e. the primitive quantity
2938 : : // vector. However, it stores the consistently obtained values of evolved
2939 : : // quantities, since nprim is the number of evolved quantities that need to
2940 : : // be evaluated consistently. For this reason, accessing R will require
2941 : : // the primitive quantity accessors. But this access is intended to give
2942 : : // the corresponding evolved quantites, as follows:
2943 : : // multispecies::tempratureIdx() - total energy
2944 [ + - ]: 1364400 : std::vector< tk::real > R(nprim*rdof, 0.0);
2945 : :
2946 [ + - ]: 682200 : auto ng = tk::NGvol(rdof);
2947 : :
2948 : : // Arrays for quadrature points
2949 : 1364400 : std::array< std::vector< tk::real >, 3 > coordgp;
2950 : 1364400 : std::vector< tk::real > wgp;
2951 : :
2952 [ + - ]: 682200 : coordgp[0].resize( ng );
2953 [ + - ]: 682200 : coordgp[1].resize( ng );
2954 [ + - ]: 682200 : coordgp[2].resize( ng );
2955 [ + - ]: 682200 : wgp.resize( ng );
2956 : :
2957 [ + - ]: 682200 : tk::GaussQuadratureTet( ng, coordgp, wgp );
2958 : :
2959 : : // Loop over quadrature points in element e
2960 [ + + ]: 4093200 : for (std::size_t igp=0; igp<ng; ++igp) {
2961 : : // Compute the basis function
2962 : 6822000 : auto B = tk::eval_basis( rdof, coordgp[0][igp], coordgp[1][igp],
2963 [ + - ]: 13644000 : coordgp[2][igp] );
2964 : :
2965 : 3411000 : auto w = wgp[igp] * vole;
2966 : :
2967 : : // Evaluate the solution at quadrature point
2968 : : auto state = evalPolynomialSol(mat_blk, 0, ncomp, nprim,
2969 : : rdof, 1, e, rdof, inpoel, coord, geoElem,
2970 [ + - ]: 6822000 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, unk, prim);
2971 : :
2972 : : // Solution vector that stores the material energy and bulk momentum
2973 [ + - ]: 6822000 : std::vector< tk::real > s(nprim, 0.0);
2974 : :
2975 : : // Mixture state at quadrature point
2976 [ + - ]: 6822000 : Mixture mixgp(nspec, state, mat_blk);
2977 : :
2978 : : // Mixture density at quadrature point
2979 : 3411000 : tk::real rhob = mixgp.get_mix_density();
2980 : :
2981 : : // velocity vector at quadrature point
2982 : : std::array< tk::real, 3 >
2983 : 3411000 : vel{ state[multispecies::momentumIdx(nspec,0)]/rhob,
2984 : 3411000 : state[multispecies::momentumIdx(nspec,1)]/rhob,
2985 : 6822000 : state[multispecies::momentumIdx(nspec,2)]/rhob };
2986 : :
2987 : : // Compute and store total energy at quadrature point
2988 [ + - ]: 3411000 : s[multispecies::temperatureIdx(nspec,0)] = mixgp.totalenergy(rhob,
2989 : 3411000 : vel[0], vel[1], vel[2],
2990 : 3411000 : state[ncomp+multispecies::temperatureIdx(nspec,0)], mat_blk);
2991 : :
2992 : : // Evaluate the right-hand-side vector
2993 [ + + ]: 6822000 : for(std::size_t k = 0; k < nprim; k++) {
2994 : 3411000 : auto mark = k * rdof;
2995 [ + + ]: 17055000 : for(std::size_t idof = 0; idof < rdof; idof++)
2996 : 13644000 : R[mark+idof] += w * s[k] * B[idof];
2997 : : }
2998 : : }
2999 : :
3000 : : // Update the high order dofs of the total energy
3001 [ + + ]: 2728800 : for(std::size_t idof = 1; idof < rdof; idof++)
3002 [ + - ]: 2046600 : unk(e, multispecies::energyDofIdx(nspec,0,rdof,idof)) =
3003 : 2046600 : R[multispecies::temperatureDofIdx(nspec,0,rdof,idof)] / (vole*L[idof]);
3004 : : }
3005 : 6600 : }
3006 : :
3007 : : tk::real
3008 : 74684686 : constrain_pressure( const std::vector< EOS >& mat_blk,
3009 : : tk::real apr,
3010 : : tk::real arho,
3011 : : tk::real alpha=1.0,
3012 : : std::size_t imat=0 )
3013 : : // *****************************************************************************
3014 : : // Constrain material partial pressure (alpha_k * p_k)
3015 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
3016 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
3017 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for the
3018 : : //! single-material system, this argument can be left unspecified by the
3019 : : //! calling code
3020 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
3021 : : //! for the single-material system, this argument can be left unspecified by
3022 : : //! the calling code
3023 : : //! \return Constrained material partial pressure (alpha_k * p_k)
3024 : : // *****************************************************************************
3025 : : {
3026 : 74684686 : return std::max(apr, alpha*mat_blk[imat].compute<
3027 [ + - ]: 74684686 : EOS::min_eff_pressure >(1e-12, arho, alpha));
3028 : : }
3029 : :
3030 : :
3031 : : } // inciter::
|
