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