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