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