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