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