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 : :
28 : : namespace inciter {
29 : :
30 : : extern ctr::InputDeck g_inputdeck;
31 : :
32 : : void
33 : 6000 : WENO_P1( const std::vector< int >& esuel,
34 : : inciter::ncomp_t offset,
35 : : tk::Fields& U )
36 : : // *****************************************************************************
37 : : // Weighted Essentially Non-Oscillatory (WENO) limiter for DGP1
38 : : //! \param[in] esuel Elements surrounding elements
39 : : //! \param[in] offset Index for equation systems
40 : : //! \param[in,out] U High-order solution vector which gets limited
41 : : //! \details This WENO function should be called for transport and compflow
42 : : //! \note This limiter function is experimental and untested. Use with caution.
43 : : // *****************************************************************************
44 : : {
45 : 6000 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
46 : 6000 : const auto cweight = inciter::g_inputdeck.get< tag::discr, tag::cweight >();
47 : 6000 : auto nelem = esuel.size()/4;
48 : : std::array< std::vector< tk::real >, 3 >
49 : : limU {{ std::vector< tk::real >(nelem),
50 : : std::vector< tk::real >(nelem),
51 [ + - ][ + - ]: 12000 : std::vector< tk::real >(nelem) }};
[ + - ]
52 : :
53 : 6000 : std::size_t ncomp = U.nprop()/rdof;
54 : :
55 [ + + ]: 12000 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
56 : : {
57 [ + + ]: 552450 : for (std::size_t e=0; e<nelem; ++e)
58 : : {
59 [ + - ]: 546450 : WENOLimiting(U, esuel, e, c, rdof, offset, cweight, limU);
60 : : }
61 : :
62 : 6000 : auto mark = c*rdof;
63 : :
64 [ + + ]: 552450 : for (std::size_t e=0; e<nelem; ++e)
65 : : {
66 [ + - ]: 546450 : U(e, mark+1, offset) = limU[0][e];
67 [ + - ]: 546450 : U(e, mark+2, offset) = limU[1][e];
68 [ + - ]: 546450 : U(e, mark+3, offset) = limU[2][e];
69 : : }
70 : : }
71 : 6000 : }
72 : :
73 : : void
74 : 24030 : Superbee_P1( const std::vector< int >& esuel,
75 : : const std::vector< std::size_t >& inpoel,
76 : : const std::vector< std::size_t >& ndofel,
77 : : inciter::ncomp_t offset,
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] offset Index for equation systems
86 : : //! \param[in] coord Array of nodal coordinates
87 : : //! \param[in,out] U High-order solution vector which gets limited
88 : : //! \details This Superbee function should be called for transport and compflow
89 : : // *****************************************************************************
90 : : {
91 : 24030 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
92 : 24030 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
93 : 24030 : std::size_t ncomp = U.nprop()/rdof;
94 : :
95 : 24030 : auto beta_lim = 2.0;
96 : :
97 [ + + ]: 5432850 : for (std::size_t e=0; e<esuel.size()/4; ++e)
98 : : {
99 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
100 : : // basis functions and solutions by default. This implies that P0P1 is
101 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
102 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
103 : : // element for pDG.
104 : : std::size_t dof_el;
105 [ + + ]: 5408820 : if (rdof > ndof)
106 : : {
107 : 2094150 : dof_el = rdof;
108 : : }
109 : : else
110 : : {
111 : 3314670 : dof_el = ndofel[e];
112 : : }
113 : :
114 [ + + ]: 5408820 : if (dof_el > 1)
115 : : {
116 : : auto phi = SuperbeeLimiting(U, esuel, inpoel, coord, e, ndof, rdof,
117 [ + - ]: 9583206 : dof_el, offset, ncomp, beta_lim);
118 : :
119 : : // apply limiter function
120 [ + + ]: 13449018 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
121 : : {
122 : 8657415 : auto mark = c*rdof;
123 [ + - ][ + - ]: 8657415 : U(e, mark+1, offset) = phi[c] * U(e, mark+1, offset);
124 [ + - ][ + - ]: 8657415 : U(e, mark+2, offset) = phi[c] * U(e, mark+2, offset);
125 [ + - ][ + - ]: 8657415 : U(e, mark+3, offset) = phi[c] * U(e, mark+3, offset);
126 : : }
127 : : }
128 : : }
129 : 24030 : }
130 : :
131 : : void
132 : 0 : SuperbeeMultiMat_P1(
133 : : const std::vector< int >& esuel,
134 : : const std::vector< std::size_t >& inpoel,
135 : : const std::vector< std::size_t >& ndofel,
136 : : std::size_t system,
137 : : inciter::ncomp_t offset,
138 : : const tk::UnsMesh::Coords& coord,
139 : : tk::Fields& U,
140 : : tk::Fields& P,
141 : : std::size_t nmat )
142 : : // *****************************************************************************
143 : : // Superbee limiter for multi-material DGP1
144 : : //! \param[in] esuel Elements surrounding elements
145 : : //! \param[in] inpoel Element connectivity
146 : : //! \param[in] ndofel Vector of local number of degrees of freedom
147 : : //! \param[in] system Index for equation systems
148 : : //! \param[in] offset Offset this PDE system operates from
149 : : //! \param[in] coord Array of nodal coordinates
150 : : //! \param[in,out] U High-order solution vector which gets limited
151 : : //! \param[in,out] P High-order vector of primitives which gets limited
152 : : //! \param[in] nmat Number of materials in this PDE system
153 : : //! \details This Superbee function should be called for multimat
154 : : // *****************************************************************************
155 : : {
156 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
157 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
158 : : const auto intsharp = inciter::g_inputdeck.get< tag::param, tag::multimat,
159 : 0 : tag::intsharp >()[system];
160 : 0 : std::size_t ncomp = U.nprop()/rdof;
161 : 0 : std::size_t nprim = P.nprop()/rdof;
162 : :
163 : 0 : auto beta_lim = 2.0;
164 : :
165 [ - - ]: 0 : for (std::size_t e=0; e<esuel.size()/4; ++e)
166 : : {
167 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
168 : : // basis functions and solutions by default. This implies that P0P1 is
169 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
170 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
171 : : // element for pDG.
172 : : std::size_t dof_el;
173 [ - - ]: 0 : if (rdof > ndof)
174 : : {
175 : 0 : dof_el = rdof;
176 : : }
177 : : else
178 : : {
179 : 0 : dof_el = ndofel[e];
180 : : }
181 : :
182 [ - - ]: 0 : if (dof_el > 1)
183 : : {
184 : : // limit conserved quantities
185 : : auto phic = SuperbeeLimiting(U, esuel, inpoel, coord, e, ndof, rdof,
186 [ - - ]: 0 : dof_el, offset, ncomp, beta_lim);
187 : : // limit primitive quantities
188 : : auto phip = SuperbeeLimiting(P, esuel, inpoel, coord, e, ndof, rdof,
189 [ - - ]: 0 : dof_el, offset, nprim, beta_lim);
190 : :
191 [ - - ]: 0 : if(ndof > 1)
192 [ - - ]: 0 : BoundPreservingLimiting(nmat, offset, ndof, e, inpoel, coord, U, phic);
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), offset);
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 : consistentMultiMatLimiting_P1(nmat, offset, rdof, e, U, P, phic, phip);
211 : : }
212 : :
213 : : // apply limiter function
214 [ - - ]: 0 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
215 : : {
216 : 0 : auto mark = c*rdof;
217 [ - - ][ - - ]: 0 : U(e, mark+1, offset) = phic[c] * U(e, mark+1, offset);
218 [ - - ][ - - ]: 0 : U(e, mark+2, offset) = phic[c] * U(e, mark+2, offset);
219 [ - - ][ - - ]: 0 : U(e, mark+3, offset) = phic[c] * U(e, mark+3, offset);
220 : : }
221 [ - - ]: 0 : for (inciter::ncomp_t c=0; c<nprim; ++c)
222 : : {
223 : 0 : auto mark = c*rdof;
224 [ - - ][ - - ]: 0 : P(e, mark+1, offset) = phip[c] * P(e, mark+1, offset);
225 [ - - ][ - - ]: 0 : P(e, mark+2, offset) = phip[c] * P(e, mark+2, offset);
226 [ - - ][ - - ]: 0 : P(e, mark+3, offset) = phip[c] * P(e, mark+3, offset);
227 : : }
228 : : }
229 : : }
230 : 0 : }
231 : :
232 : : void
233 : 0 : VertexBasedTransport_P1(
234 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
235 : : const std::vector< std::size_t >& inpoel,
236 : : const std::vector< std::size_t >& ndofel,
237 : : std::size_t nelem,
238 : : std::size_t system,
239 : : std::size_t offset,
240 : : const tk::Fields& geoElem,
241 : : const tk::UnsMesh::Coords& coord,
242 : : tk::Fields& U )
243 : : // *****************************************************************************
244 : : // Kuzmin's vertex-based limiter for transport DGP1
245 : : //! \param[in] esup Elements surrounding points
246 : : //! \param[in] inpoel Element connectivity
247 : : //! \param[in] ndofel Vector of local number of degrees of freedom
248 : : //! \param[in] nelem Number of elements
249 : : //! \param[in] system Index for equation systems
250 : : //! \param[in] offset Index for equation systems
251 : : //! \param[in] geoElem Element geometry array
252 : : //! \param[in] coord Array of nodal coordinates
253 : : //! \param[in,out] U High-order solution vector which gets limited
254 : : //! \details This vertex-based limiter function should be called for transport.
255 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
256 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
257 : : //! computational and applied mathematics, 233(12), 3077-3085.
258 : : // *****************************************************************************
259 : : {
260 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
261 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
262 : : const auto intsharp = inciter::g_inputdeck.get< tag::param, tag::transport,
263 : 0 : tag::intsharp >()[system];
264 : 0 : std::size_t ncomp = U.nprop()/rdof;
265 : :
266 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
267 : : {
268 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
269 : : // basis functions and solutions by default. This implies that P0P1 is
270 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
271 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
272 : : // element for pDG.
273 : : std::size_t dof_el;
274 [ - - ]: 0 : if (rdof > ndof)
275 : : {
276 : 0 : dof_el = rdof;
277 : : }
278 : : else
279 : : {
280 : 0 : dof_el = ndofel[e];
281 : : }
282 : :
283 [ - - ]: 0 : if (dof_el > 1)
284 : : {
285 : 0 : std::vector< std::vector< tk::real > > unk;
286 [ - - ]: 0 : std::vector< tk::real > phi(ncomp, 1.0);
287 : : // limit conserved quantities
288 : 0 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
289 [ - - ]: 0 : dof_el, offset, ncomp, phi, {0, ncomp-1});
290 : :
291 : : // limits under which compression is to be performed
292 [ - - ]: 0 : std::vector< std::size_t > matInt(ncomp, 0);
293 [ - - ]: 0 : std::vector< tk::real > alAvg(ncomp, 0.0);
294 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k)
295 [ - - ]: 0 : alAvg[k] = U(e,k*rdof,offset);
296 [ - - ]: 0 : auto intInd = interfaceIndicator(ncomp, alAvg, matInt);
297 [ - - ][ - - ]: 0 : if ((intsharp > 0) && intInd)
298 : : {
299 [ - - ]: 0 : for (std::size_t k=0; k<ncomp; ++k)
300 : : {
301 [ - - ]: 0 : if (matInt[k]) phi[k] = 1.0;
302 : : }
303 : : }
304 : :
305 : : // apply limiter function
306 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
307 : : {
308 : 0 : auto mark = c*rdof;
309 [ - - ][ - - ]: 0 : U(e, mark+1, offset) = phi[c] * U(e, mark+1, offset);
310 [ - - ][ - - ]: 0 : U(e, mark+2, offset) = phi[c] * U(e, mark+2, offset);
311 [ - - ][ - - ]: 0 : U(e, mark+3, offset) = phi[c] * U(e, mark+3, offset);
312 : : }
313 : : }
314 : : }
315 : 0 : }
316 : :
317 : : void
318 : 0 : VertexBasedCompflow_P1(
319 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
320 : : const std::vector< std::size_t >& inpoel,
321 : : const std::vector< std::size_t >& ndofel,
322 : : std::size_t nelem,
323 : : std::size_t offset,
324 : : const tk::Fields& geoElem,
325 : : const tk::UnsMesh::Coords& coord,
326 : : tk::Fields& U )
327 : : // *****************************************************************************
328 : : // Kuzmin's vertex-based limiter for single-material DGP1
329 : : //! \param[in] esup Elements surrounding points
330 : : //! \param[in] inpoel Element connectivity
331 : : //! \param[in] ndofel Vector of local number of degrees of freedom
332 : : //! \param[in] nelem Number of elements
333 : : //! \param[in] offset Index for equation systems
334 : : //! \param[in] geoElem Element geometry array
335 : : //! \param[in] coord Array of nodal coordinates
336 : : //! \param[in,out] U High-order solution vector which gets limited
337 : : //! \details This vertex-based limiter function should be called for compflow.
338 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
339 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
340 : : //! computational and applied mathematics, 233(12), 3077-3085.
341 : : // *****************************************************************************
342 : : {
343 : 0 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
344 : 0 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
345 : 0 : std::size_t ncomp = U.nprop()/rdof;
346 : :
347 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
348 : : {
349 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
350 : : // basis functions and solutions by default. This implies that P0P1 is
351 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
352 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
353 : : // element for pDG.
354 : : std::size_t dof_el;
355 [ - - ]: 0 : if (rdof > ndof)
356 : : {
357 : 0 : dof_el = rdof;
358 : : }
359 : : else
360 : : {
361 : 0 : dof_el = ndofel[e];
362 : : }
363 : :
364 [ - - ]: 0 : if (dof_el > 1)
365 : : {
366 : 0 : std::vector< std::vector< tk::real > > unk;
367 [ - - ]: 0 : std::vector< tk::real > phi(ncomp, 1.0);
368 : : // limit conserved quantities
369 : 0 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
370 [ - - ]: 0 : dof_el, offset, ncomp, phi, {0, ncomp-1});
371 : :
372 : : // apply limiter function
373 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
374 : : {
375 : 0 : auto mark = c*rdof;
376 [ - - ][ - - ]: 0 : U(e, mark+1, offset) = phi[c] * U(e, mark+1, offset);
377 [ - - ][ - - ]: 0 : U(e, mark+2, offset) = phi[c] * U(e, mark+2, offset);
378 [ - - ][ - - ]: 0 : U(e, mark+3, offset) = phi[c] * U(e, mark+3, offset);
379 : : }
380 : : }
381 : : }
382 : 0 : }
383 : :
384 : : void
385 : 300 : VertexBasedCompflow_P2(
386 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
387 : : const std::vector< std::size_t >& inpoel,
388 : : const std::vector< std::size_t >& ndofel,
389 : : std::size_t nelem,
390 : : std::size_t offset,
391 : : const tk::Fields& geoElem,
392 : : const tk::UnsMesh::Coords& coord,
393 : : const std::vector< std::size_t >& gid,
394 : : const std::unordered_map< std::size_t, std::size_t >& bid,
395 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
396 : : tk::Fields& U )
397 : : // *****************************************************************************
398 : : // Kuzmin's vertex-based limiter for single-material DGP2
399 : : //! \param[in] esup Elements surrounding points
400 : : //! \param[in] inpoel Element connectivity
401 : : //! \param[in] ndofel Vector of local number of degrees of freedom
402 : : //! \param[in] nelem Number of elements
403 : : //! \param[in] offset Index for equation systems
404 : : //! \param[in] geoElem Element geometry array
405 : : //! \param[in] coord Array of nodal coordinates
406 : : //! \param[in] gid Local->global node id map
407 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
408 : : //! global node ids (key)
409 : : //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
410 : : //! variables
411 : : //! \param[in,out] U High-order solution vector which gets limited
412 : : //! \details This vertex-based limiter function should be called for compflow.
413 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
414 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
415 : : //! computational and applied mathematics, 233(12), 3077-3085.
416 : : // *****************************************************************************
417 : : {
418 : 300 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
419 : 300 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
420 : 300 : std::size_t ncomp = U.nprop()/rdof;
421 : :
422 : : // Copy field data U to U_lim. U_lim will store the limited solution
423 : : // temporarily, so that the limited solution is NOT used to find the
424 : : // min/max bounds for the limiting function
425 [ + - ]: 600 : auto U_lim = U;
426 : :
427 [ + + ]: 437460 : for (std::size_t e=0; e<nelem; ++e)
428 : : {
429 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
430 : : // basis functions and solutions by default. This implies that P0P1 is
431 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
432 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
433 : : // element for pDG.
434 : : std::size_t dof_el;
435 [ - + ]: 437160 : if (rdof > ndof)
436 : : {
437 : 0 : dof_el = rdof;
438 : : }
439 : : else
440 : : {
441 : 437160 : dof_el = ndofel[e];
442 : : }
443 : :
444 : 437160 : bool shock_detec(false);
445 : :
446 : : // Evaluate the shock detection indicator
447 [ + - ]: 437160 : auto Ind = evalDiscIndicator_CompFlow(e, ncomp, dof_el, ndofel[e], U);
448 [ + + ]: 437160 : if(Ind > 1e-6)
449 : 3604 : shock_detec = true;
450 : :
451 [ + - ][ + + ]: 437160 : if (dof_el > 1 && shock_detec)
452 : : {
453 : : // Transform the solution with Dubiner basis to Taylor basis so that the
454 : : // limiting function could be applied to physical derivatives in a
455 : : // hierarchical manner
456 : : auto unk =
457 [ + - ]: 7208 : tk::DubinerToTaylor(ncomp, offset, e, dof_el, U, inpoel, coord);
458 : :
459 : : // The vector of limiting coefficients for P1 and P2 coefficients
460 [ + - ]: 7208 : std::vector< tk::real > phic_p1(ncomp, 1.0);
461 [ + - ]: 7208 : std::vector< tk::real > phic_p2(ncomp, 1.0);
462 : :
463 : : // If DGP2 is applied, apply the limiter function to the first derivative
464 : : // to obtain the limiting coefficient for P2 coefficients
465 [ + - ]: 3604 : if(dof_el > 4)
466 [ + - ]: 7208 : phic_p2 = VertexBasedLimiting_P2(unk, U, esup, inpoel, coord, geoElem,
467 : 3604 : e, rdof, dof_el, offset, ncomp, gid, bid, uNodalExtrm);
468 : :
469 : : // limit conserved quantities
470 : 0 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
471 [ + - ]: 3604 : dof_el, offset, ncomp, phic_p1, {0, ncomp-1});
472 : :
473 [ + - ]: 3604 : if(dof_el > 4)
474 [ + + ]: 21624 : for (std::size_t c=0; c<ncomp; ++c)
475 : 18020 : phic_p1[c] = std::max(phic_p1[c], phic_p2[c]);
476 : :
477 : : // apply limiter function to the solution with Taylor basis
478 [ + + ]: 21624 : for (std::size_t c=0; c<ncomp; ++c)
479 : : {
480 : 18020 : unk[c][1] = phic_p1[c] * unk[c][1];
481 : 18020 : unk[c][2] = phic_p1[c] * unk[c][2];
482 : 18020 : unk[c][3] = phic_p1[c] * unk[c][3];
483 : : }
484 [ + - ]: 3604 : if(dof_el > 4)
485 : : {
486 [ + + ]: 21624 : for (std::size_t c=0; c<ncomp; ++c)
487 : : {
488 : 18020 : unk[c][4] = phic_p2[c] * unk[c][4];
489 : 18020 : unk[c][5] = phic_p2[c] * unk[c][5];
490 : 18020 : unk[c][6] = phic_p2[c] * unk[c][6];
491 : 18020 : unk[c][7] = phic_p2[c] * unk[c][7];
492 : 18020 : unk[c][8] = phic_p2[c] * unk[c][8];
493 : 18020 : unk[c][9] = phic_p2[c] * unk[c][9];
494 : : }
495 : : }
496 : :
497 : : // Convert the solution with Taylor basis to the solution with Dubiner
498 : : // basis
499 [ + - ]: 3604 : tk::TaylorToDubiner( ncomp, e, dof_el, inpoel, coord, geoElem, unk );
500 : :
501 : : // Store the limited solution in U_lim
502 [ + + ]: 21624 : for(std::size_t c=0; c<ncomp; ++c)
503 : : {
504 : 18020 : auto mark = c*rdof;
505 [ + + ]: 180200 : for(std::size_t idof = 1; idof < rdof; idof++)
506 [ + - ]: 162180 : U_lim(e, mark+idof, offset) = unk[c][idof];
507 : : }
508 : : }
509 : : }
510 : :
511 : : // Store the limited solution with Dubiner basis
512 [ + + ]: 437460 : for (std::size_t e=0; e<nelem; ++e)
513 : : {
514 [ + + ]: 2622960 : for (std::size_t c=0; c<ncomp; ++c)
515 : : {
516 : 2185800 : auto mark = c*rdof;
517 [ + - ][ + - ]: 2185800 : U(e, mark+1, offset) = U_lim(e, mark+1, offset);
518 [ + - ][ + - ]: 2185800 : U(e, mark+2, offset) = U_lim(e, mark+2, offset);
519 [ + - ][ + - ]: 2185800 : U(e, mark+3, offset) = U_lim(e, mark+3, offset);
520 : :
521 [ + - ]: 2185800 : if(ndof > 4)
522 : : {
523 [ + - ][ + - ]: 2185800 : U(e, mark+4, offset) = U_lim(e, mark+4, offset);
524 [ + - ][ + - ]: 2185800 : U(e, mark+5, offset) = U_lim(e, mark+5, offset);
525 [ + - ][ + - ]: 2185800 : U(e, mark+6, offset) = U_lim(e, mark+6, offset);
526 [ + - ][ + - ]: 2185800 : U(e, mark+7, offset) = U_lim(e, mark+7, offset);
527 [ + - ][ + - ]: 2185800 : U(e, mark+8, offset) = U_lim(e, mark+8, offset);
528 [ + - ][ + - ]: 2185800 : U(e, mark+9, offset) = U_lim(e, mark+9, offset);
529 : : }
530 : : }
531 : : }
532 : 300 : }
533 : :
534 : : void
535 : 4125 : VertexBasedMultiMat_P1(
536 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
537 : : const std::vector< std::size_t >& inpoel,
538 : : const std::vector< std::size_t >& ndofel,
539 : : std::size_t nelem,
540 : : std::size_t system,
541 : : std::size_t offset,
542 : : [[maybe_unused]] const inciter::FaceData& fd,
543 : : [[maybe_unused]] const tk::Fields& geoFace,
544 : : const tk::Fields& geoElem,
545 : : const tk::UnsMesh::Coords& coord,
546 : : tk::Fields& U,
547 : : tk::Fields& P,
548 : : std::size_t nmat,
549 : : std::vector< std::size_t >& shockmarker )
550 : : // *****************************************************************************
551 : : // Kuzmin's vertex-based limiter for multi-material DGP1
552 : : //! \param[in] esup Elements surrounding points
553 : : //! \param[in] inpoel Element connectivity
554 : : //! \param[in] ndofel Vector of local number of degrees of freedom
555 : : //! \param[in] nelem Number of elements
556 : : //! \param[in] system Index for equation systems
557 : : //! \param[in] offset Offset this PDE system operates from
558 : : //! \param[in] geoElem Element geometry array
559 : : //! \param[in] coord Array of nodal coordinates
560 : : //! \param[in,out] U High-order solution vector which gets limited
561 : : //! \param[in,out] P High-order vector of primitives which gets limited
562 : : //! \param[in] nmat Number of materials in this PDE system
563 : : //! \param[in,out] shockmarker Shock detection marker array
564 : : //! \details This vertex-based limiter function should be called for multimat.
565 : : //! For details see: Kuzmin, D. (2010). A vertex-based hierarchical slope
566 : : //! limiter for p-adaptive discontinuous Galerkin methods. Journal of
567 : : //! computational and applied mathematics, 233(12), 3077-3085.
568 : : // *****************************************************************************
569 : : {
570 : 4125 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
571 : 4125 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
572 : : const auto intsharp = inciter::g_inputdeck.get< tag::param, tag::multimat,
573 : 4125 : tag::intsharp >()[system];
574 : 4125 : std::size_t ncomp = U.nprop()/rdof;
575 : 4125 : std::size_t nprim = P.nprop()/rdof;
576 : :
577 : : // Evaluate the interface condition and mark the shock cells
578 : : //MarkShockCells(nelem, nmat, system, offset, ndof, rdof, ndofel, inpoel, coord,
579 : : // fd, geoFace, geoElem, U, P, shockmarker);
580 : :
581 : : // Threshold for shock detection indicator
582 : 4125 : auto threshold = pow(10, -5.7);
583 : :
584 [ + + ]: 913725 : for (std::size_t e=0; e<nelem; ++e)
585 : : {
586 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
587 : : // basis functions and solutions by default. This implies that P0P1 is
588 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
589 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
590 : : // element for pDG.
591 : : std::size_t dof_el;
592 [ + + ]: 909600 : if (rdof > ndof)
593 : : {
594 : 454800 : dof_el = rdof;
595 : : }
596 : : else
597 : : {
598 : 454800 : dof_el = ndofel[e];
599 : : }
600 : :
601 [ + + ]: 909600 : if(ndofel[e] > 1) {
602 : : // Evaluate the shock detection indicator to determine whether the limiter
603 : : // is applied or not
604 : 454800 : auto Ind = evalDiscIndicator_MultiMat(e, nmat, ncomp, dof_el, ndofel[e], U);
605 [ + + ]: 454800 : if(Ind > threshold)
606 : 19010 : shockmarker[e] = 1;
607 : : else
608 : 435790 : shockmarker[e] = 0;
609 : : } else { // If P0P1, the limiter is always applied
610 : 454800 : shockmarker[e] = 1;
611 : : }
612 : :
613 [ + - ]: 909600 : if (dof_el > 1)
614 : : {
615 : 1819200 : std::vector< std::vector< tk::real > > unk;
616 [ + - ]: 1819200 : std::vector< tk::real > phic(ncomp, 1.0);
617 [ + - ]: 1819200 : std::vector< tk::real > phip(nprim, 1.0);
618 [ + + ]: 909600 : if(shockmarker[e]) {
619 : : // When shockmarker is 1, there is discontinuity within the element.
620 : : // Hence, the vertex-based limiter will be applied.
621 : :
622 : : // limit conserved quantities
623 : 0 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
624 [ + - ]: 473810 : dof_el, offset, ncomp, phic, {0, ncomp-1});
625 : : // limit primitive quantities
626 : 473810 : VertexBasedLimiting(unk, P, esup, inpoel, coord, geoElem, e, rdof,
627 [ + - ]: 473810 : dof_el, offset, nprim, phip, {0, nprim-1});
628 : : } else {
629 : : // When shockmarker is 0, the volume fraction, density and energy
630 : : // of minor material will still be limited to ensure a stable solution.
631 [ + - ]: 435790 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
632 : : dof_el, offset, ncomp, phic,
633 : 435790 : {volfracIdx(nmat,0), volfracIdx(nmat,nmat-1)});
634 : :
635 [ + + ]: 1307370 : for(std::size_t k=0; k<nmat; ++k) {
636 [ + - ][ + + ]: 871580 : if(U(e, volfracDofIdx(nmat,k,rdof,0), offset) < 1e-4) {
637 : : // Vector to store the range of limited variables
638 : : std::array< std::size_t, 2 > VarRange;
639 : :
640 : : // limit the density of minor materials
641 : 427000 : VarRange[0] = densityIdx(nmat, k);
642 : 427000 : VarRange[1] = VarRange[0];
643 [ + - ]: 427000 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
644 : : dof_el, offset, ncomp, phic, VarRange);
645 : :
646 : : // limit the energy of minor materials
647 : 427000 : VarRange[0] = energyIdx(nmat, k);
648 : 427000 : VarRange[1] = VarRange[0];
649 [ + - ]: 427000 : VertexBasedLimiting(unk, U, esup, inpoel, coord, geoElem, e, rdof,
650 : : dof_el, offset, ncomp, phic, VarRange);
651 : :
652 : : // limit the pressure of minor materials
653 : 427000 : VarRange[0] = pressureIdx(nmat, k);
654 : 427000 : VarRange[1] = VarRange[0];
655 [ + - ]: 427000 : VertexBasedLimiting(unk, P, esup, inpoel, coord, geoElem, e, rdof,
656 : : dof_el, offset, nprim, phip, VarRange);
657 : : }
658 : : }
659 : : }
660 : :
661 [ + + ][ + + ]: 909600 : if(ndof > 1 && intsharp == 0)
662 [ + - ]: 227400 : BoundPreservingLimiting(nmat, offset, ndof, e, inpoel, coord, U, phic);
663 : :
664 : : // limits under which compression is to be performed
665 [ + - ]: 1819200 : std::vector< std::size_t > matInt(nmat, 0);
666 [ + - ]: 1819200 : std::vector< tk::real > alAvg(nmat, 0.0);
667 [ + + ]: 2728800 : for (std::size_t k=0; k<nmat; ++k)
668 [ + - ]: 1819200 : alAvg[k] = U(e, volfracDofIdx(nmat,k,rdof,0), offset);
669 [ + - ]: 909600 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
670 [ + + ][ + + ]: 909600 : if ((intsharp > 0) && intInd)
671 : : {
672 [ + + ]: 32190 : for (std::size_t k=0; k<nmat; ++k)
673 : : {
674 [ + - ]: 21460 : if (matInt[k])
675 : 21460 : phic[volfracIdx(nmat,k)] = 1.0;
676 : 10730 : }
677 : : }
678 : : else
679 : : {
680 [ + - ]: 898870 : consistentMultiMatLimiting_P1(nmat, offset, rdof, e, U, P, phic, phip);
681 : : }
682 : :
683 : : // apply limiter function
684 [ + + ]: 9096000 : for (std::size_t c=0; c<ncomp; ++c)
685 : : {
686 : 8186400 : auto mark = c*rdof;
687 [ + - ][ + - ]: 8186400 : U(e, mark+1, offset) = phic[c] * U(e, mark+1, offset);
688 [ + - ][ + - ]: 8186400 : U(e, mark+2, offset) = phic[c] * U(e, mark+2, offset);
689 [ + - ][ + - ]: 8186400 : U(e, mark+3, offset) = phic[c] * U(e, mark+3, offset);
690 : : }
691 [ + + ]: 5457600 : for (std::size_t c=0; c<nprim; ++c)
692 : : {
693 : 4548000 : auto mark = c*rdof;
694 [ + - ][ + - ]: 4548000 : P(e, mark+1, offset) = phip[c] * P(e, mark+1, offset);
695 [ + - ][ + - ]: 4548000 : P(e, mark+2, offset) = phip[c] * P(e, mark+2, offset);
696 [ + - ][ + - ]: 4548000 : P(e, mark+3, offset) = phip[c] * P(e, mark+3, offset);
697 : : }
698 : : }
699 : : }
700 : 4125 : }
701 : :
702 : : void
703 : 546450 : WENOLimiting( const tk::Fields& U,
704 : : const std::vector< int >& esuel,
705 : : std::size_t e,
706 : : inciter::ncomp_t c,
707 : : std::size_t rdof,
708 : : inciter::ncomp_t offset,
709 : : tk::real cweight,
710 : : std::array< std::vector< tk::real >, 3 >& limU )
711 : : // *****************************************************************************
712 : : // WENO limiter function calculation for P1 dofs
713 : : //! \param[in] U High-order solution vector which is to be limited
714 : : //! \param[in] esuel Elements surrounding elements
715 : : //! \param[in] e Id of element whose solution is to be limited
716 : : //! \param[in] c Index of component which is to be limited
717 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
718 : : //! \param[in] offset Index for equation systems
719 : : //! \param[in] cweight Weight of the central stencil
720 : : //! \param[in,out] limU Limited gradients of component c
721 : : // *****************************************************************************
722 : : {
723 : : std::array< std::array< tk::real, 3 >, 5 > gradu;
724 : : std::array< tk::real, 5 > wtStencil, osc, wtDof;
725 : :
726 : 546450 : auto mark = c*rdof;
727 : :
728 : : // reset all stencil values to zero
729 [ + + ][ + - ]: 3278700 : for (auto& g : gradu) g.fill(0.0);
730 [ + - ]: 546450 : osc.fill(0);
731 [ + - ]: 546450 : wtDof.fill(0);
732 [ + - ]: 546450 : wtStencil.fill(0);
733 : :
734 : : // The WENO limiter uses solution data from the neighborhood in the form
735 : : // of stencils to enforce non-oscillatory conditions. The immediate
736 : : // (Von Neumann) neighborhood of a tetrahedral cell consists of the 4
737 : : // cells that share faces with it. These are the 4 neighborhood-stencils
738 : : // for the tetrahedron. The primary stencil is the tet itself. Weights are
739 : : // assigned to these stencils, with the primary stencil usually assigned
740 : : // the highest weight. The lower the primary/central weight, the more
741 : : // dissipative the limiting effect. This central weight is usually problem
742 : : // dependent. It is set higher for relatively weaker discontinuities, and
743 : : // lower for stronger discontinuities.
744 : :
745 : : // primary stencil
746 [ + - ]: 546450 : gradu[0][0] = U(e, mark+1, offset);
747 [ + - ]: 546450 : gradu[0][1] = U(e, mark+2, offset);
748 [ + - ]: 546450 : gradu[0][2] = U(e, mark+3, offset);
749 : 546450 : wtStencil[0] = cweight;
750 : :
751 : : // stencils from the neighborhood
752 [ + + ]: 2732250 : for (std::size_t is=1; is<5; ++is)
753 : : {
754 : 2185800 : auto nel = esuel[ 4*e+(is-1) ];
755 : :
756 : : // ignore physical domain ghosts
757 [ + + ]: 2185800 : if (nel == -1)
758 : : {
759 [ + - ]: 359700 : gradu[is].fill(0.0);
760 : 359700 : wtStencil[is] = 0.0;
761 : 359700 : continue;
762 : : }
763 : :
764 : 1826100 : std::size_t n = static_cast< std::size_t >( nel );
765 [ + - ]: 1826100 : gradu[is][0] = U(n, mark+1, offset);
766 [ + - ]: 1826100 : gradu[is][1] = U(n, mark+2, offset);
767 [ + - ]: 1826100 : gradu[is][2] = U(n, mark+3, offset);
768 : 1826100 : wtStencil[is] = 1.0;
769 : : }
770 : :
771 : : // From these stencils, an oscillation indicator is calculated, which
772 : : // determines the effective weights for the high-order solution DOFs.
773 : : // These effective weights determine the contribution of each of the
774 : : // stencils to the high-order solution DOFs of the current cell which are
775 : : // being limited. If this indicator detects a large oscillation in the
776 : : // solution of the current cell, it reduces the effective weight for the
777 : : // central stencil contribution to its high-order DOFs. This results in
778 : : // a more dissipative and well-behaved solution in the troubled cell.
779 : :
780 : : // oscillation indicators
781 [ + + ]: 3278700 : for (std::size_t is=0; is<5; ++is)
782 : 2732250 : osc[is] = std::sqrt( tk::dot(gradu[is], gradu[is]) );
783 : :
784 : 546450 : tk::real wtotal = 0;
785 : :
786 : : // effective weights for dofs
787 [ + + ]: 3278700 : for (std::size_t is=0; is<5; ++is)
788 : : {
789 : : // A small number (1.0e-8) is needed here to avoid dividing by a zero in
790 : : // the case of a constant solution, where osc would be zero. The number
791 : : // is not set to machine zero because it is squared, and a number
792 : : // between 1.0e-8 to 1.0e-6 is needed.
793 : 2732250 : wtDof[is] = wtStencil[is] * pow( (1.0e-8 + osc[is]), -2 );
794 : 2732250 : wtotal += wtDof[is];
795 : : }
796 : :
797 [ + + ]: 3278700 : for (std::size_t is=0; is<5; ++is)
798 : : {
799 : 2732250 : wtDof[is] = wtDof[is]/wtotal;
800 : : }
801 : :
802 : 546450 : limU[0][e] = 0.0;
803 : 546450 : limU[1][e] = 0.0;
804 : 546450 : limU[2][e] = 0.0;
805 : :
806 : : // limiter function
807 [ + + ]: 3278700 : for (std::size_t is=0; is<5; ++is)
808 : : {
809 : 2732250 : limU[0][e] += wtDof[is]*gradu[is][0];
810 : 2732250 : limU[1][e] += wtDof[is]*gradu[is][1];
811 : 2732250 : limU[2][e] += wtDof[is]*gradu[is][2];
812 : : }
813 : 546450 : }
814 : :
815 : : std::vector< tk::real >
816 : 4791603 : SuperbeeLimiting( const tk::Fields& U,
817 : : const std::vector< int >& esuel,
818 : : const std::vector< std::size_t >& inpoel,
819 : : const tk::UnsMesh::Coords& coord,
820 : : std::size_t e,
821 : : std::size_t ndof,
822 : : std::size_t rdof,
823 : : std::size_t dof_el,
824 : : inciter::ncomp_t offset,
825 : : inciter:: ncomp_t ncomp,
826 : : tk::real beta_lim )
827 : : // *****************************************************************************
828 : : // Superbee limiter function calculation for P1 dofs
829 : : //! \param[in] U High-order solution vector which is to be limited
830 : : //! \param[in] esuel Elements surrounding elements
831 : : //! \param[in] inpoel Element connectivity
832 : : //! \param[in] coord Array of nodal coordinates
833 : : //! \param[in] e Id of element whose solution is to be limited
834 : : //! \param[in] ndof Maximum number of degrees of freedom
835 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
836 : : //! \param[in] dof_el Local number of degrees of freedom
837 : : //! \param[in] offset Index for equation systems
838 : : //! \param[in] ncomp Number of scalar components in this PDE system
839 : : //! \param[in] beta_lim Parameter which is equal to 2 for Superbee and 1 for
840 : : //! minmod limiter
841 : : //! \return phi Limiter function for solution in element e
842 : : // *****************************************************************************
843 : : {
844 : : // Superbee is a TVD limiter, which uses min-max bounds that the
845 : : // high-order solution should satisfy, to ensure TVD properties. For a
846 : : // high-order method like DG, this involves the following steps:
847 : : // 1. Find min-max bounds in the immediate neighborhood of cell.
848 : : // 2. Calculate the Superbee function for all the points where solution
849 : : // needs to be reconstructed to (all quadrature points). From these,
850 : : // use the minimum value of the limiter function.
851 : :
852 [ + - ][ + - ]: 9583206 : std::vector< tk::real > uMin(ncomp, 0.0), uMax(ncomp, 0.0);
853 : :
854 [ + + ]: 13449018 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
855 : : {
856 : 8657415 : auto mark = c*rdof;
857 [ + - ]: 8657415 : uMin[c] = U(e, mark, offset);
858 [ + - ]: 8657415 : uMax[c] = U(e, mark, offset);
859 : : }
860 : :
861 : : // ----- Step-1: find min/max in the neighborhood
862 [ + + ]: 23958015 : for (std::size_t is=0; is<4; ++is)
863 : : {
864 : 19166412 : auto nel = esuel[ 4*e+is ];
865 : :
866 : : // ignore physical domain ghosts
867 [ + + ]: 19166412 : if (nel == -1) continue;
868 : :
869 : 16077936 : auto n = static_cast< std::size_t >( nel );
870 [ + + ]: 45336816 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
871 : : {
872 : 29258880 : auto mark = c*rdof;
873 [ + - ]: 29258880 : uMin[c] = std::min(uMin[c], U(n, mark, offset));
874 [ + - ]: 29258880 : uMax[c] = std::max(uMax[c], U(n, mark, offset));
875 : : }
876 : : }
877 : :
878 : : // ----- Step-2: loop over all quadrature points to get limiter function
879 : :
880 : : // to loop over all the quadrature points of all faces of element e,
881 : : // coordinates of the quadrature points are needed.
882 : : // Number of quadrature points for face integration
883 [ + - ]: 4791603 : auto ng = tk::NGfa(ndof);
884 : :
885 : : // arrays for quadrature points
886 : 9583206 : std::array< std::vector< tk::real >, 2 > coordgp;
887 : 9583206 : std::vector< tk::real > wgp;
888 : :
889 [ + - ]: 4791603 : coordgp[0].resize( ng );
890 [ + - ]: 4791603 : coordgp[1].resize( ng );
891 [ + - ]: 4791603 : wgp.resize( ng );
892 : :
893 : : // get quadrature point weights and coordinates for triangle
894 [ + - ]: 4791603 : tk::GaussQuadratureTri( ng, coordgp, wgp );
895 : :
896 : 4791603 : const auto& cx = coord[0];
897 : 4791603 : const auto& cy = coord[1];
898 : 4791603 : const auto& cz = coord[2];
899 : :
900 : : // Extract the element coordinates
901 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
902 : 14374809 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
903 : 14374809 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
904 : 14374809 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
905 : 43124427 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
906 : :
907 : : // Compute the determinant of Jacobian matrix
908 : : auto detT =
909 : 4791603 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
910 : :
911 : : // initialize limiter function
912 [ + - ]: 4791603 : std::vector< tk::real > phi(ncomp, 1.0);
913 [ + + ]: 23958015 : for (std::size_t lf=0; lf<4; ++lf)
914 : : {
915 : : // Extract the face coordinates
916 : 19166412 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
917 : 19166412 : inpoel[4*e+tk::lpofa[lf][1]],
918 : 38332824 : inpoel[4*e+tk::lpofa[lf][2]] }};
919 : :
920 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
921 : 57499236 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
922 : 57499236 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
923 : 114998472 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
924 : :
925 : : // Gaussian quadrature
926 [ + + ]: 60806364 : for (std::size_t igp=0; igp<ng; ++igp)
927 : : {
928 : : // Compute the coordinates of quadrature point at physical domain
929 [ + - ]: 41639952 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
930 : :
931 : : //Compute the basis functions
932 : : auto B_l = tk::eval_basis( rdof,
933 : 41639952 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
934 : 41639952 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
935 [ + - ]: 166559808 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
936 : :
937 [ + - ]: 83279904 : auto state = tk::eval_state( ncomp, offset, rdof, dof_el, e, U, B_l, {0, ncomp-1} );
938 : :
939 [ - + ][ - - ]: 41639952 : Assert( state.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
940 : :
941 : : // compute the limiter function
942 [ + + ]: 118691712 : for (inciter::ncomp_t c=0; c<ncomp; ++c)
943 : : {
944 : 77051760 : auto phi_gp = 1.0;
945 : 77051760 : auto mark = c*rdof;
946 [ + - ]: 77051760 : auto uNeg = state[c] - U(e, mark, offset);
947 [ + + ]: 77051760 : if (uNeg > 1.0e-14)
948 : : {
949 : 8507096 : uNeg = std::max(uNeg, 1.0e-08);
950 [ + - ]: 8507096 : phi_gp = std::min( 1.0, (uMax[c]-U(e, mark, offset))/(2.0*uNeg) );
951 : : }
952 [ + + ]: 68544664 : else if (uNeg < -1.0e-14)
953 : : {
954 : 8959703 : uNeg = std::min(uNeg, -1.0e-08);
955 [ + - ]: 8959703 : phi_gp = std::min( 1.0, (uMin[c]-U(e, mark, offset))/(2.0*uNeg) );
956 : : }
957 : : else
958 : : {
959 : 59584961 : phi_gp = 1.0;
960 : : }
961 : 154103520 : phi_gp = std::max( 0.0,
962 : 154103520 : std::max( std::min(beta_lim*phi_gp, 1.0),
963 : 77051760 : std::min(phi_gp, beta_lim) ) );
964 : 77051760 : phi[c] = std::min( phi[c], phi_gp );
965 : : }
966 : : }
967 : : }
968 : :
969 : 9583206 : return phi;
970 : : }
971 : :
972 : : void
973 : 2668014 : VertexBasedLimiting( const std::vector< std::vector< tk::real > >& unk,
974 : : const tk::Fields& U,
975 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
976 : : const std::vector< std::size_t >& inpoel,
977 : : const tk::UnsMesh::Coords& coord,
978 : : const tk::Fields& geoElem,
979 : : std::size_t e,
980 : : std::size_t rdof,
981 : : std::size_t dof_el,
982 : : std::size_t offset,
983 : : std::size_t ncomp,
984 : : std::vector< tk::real >& phi,
985 : : const std::array< std::size_t, 2 >& VarRange )
986 : : // *****************************************************************************
987 : : // Kuzmin's vertex-based limiter function calculation for P1 dofs
988 : : //! \param[in] U High-order solution vector which is to be limited
989 : : //! \param[in] esup Elements surrounding points
990 : : //! \param[in] inpoel Element connectivity
991 : : //! \param[in] coord Array of nodal coordinates
992 : : //! \param[in] geoElem Element geometry array
993 : : //! \param[in] e Id of element whose solution is to be limited
994 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
995 : : //! \param[in] dof_el Local number of degrees of freedom
996 : : //! \param[in] offset Index for equation systems
997 : : //! \param[in] ncomp Number of scalar components in this PDE system
998 : : //! \return phi Limiter function for solution in element e
999 : : // *****************************************************************************
1000 : : {
1001 : : // Kuzmin's vertex-based TVD limiter uses min-max bounds that the
1002 : : // high-order solution should satisfy, to ensure TVD properties. For a
1003 : : // high-order method like DG, this involves the following steps:
1004 : : // 1. Find min-max bounds in the nodal-neighborhood of cell.
1005 : : // 2. Calculate the limiter function (Superbee) for all the vertices of cell.
1006 : : // From these, use the minimum value of the limiter function.
1007 : :
1008 : : // Prepare for calculating Basis functions
1009 : 2668014 : const auto& cx = coord[0];
1010 : 2668014 : const auto& cy = coord[1];
1011 : 2668014 : const auto& cz = coord[2];
1012 : :
1013 : : // Extract the element coordinates
1014 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
1015 : 8004042 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
1016 : 8004042 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
1017 : 8004042 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
1018 : 24012126 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
1019 : :
1020 : : // Compute the determinant of Jacobian matrix
1021 : : auto detT =
1022 : 2668014 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
1023 : :
1024 [ + - ]: 5336028 : std::vector< tk::real > uMin(VarRange[1]-VarRange[0]+1, 0.0),
1025 [ + - ]: 5336028 : uMax(VarRange[1]-VarRange[0]+1, 0.0);
1026 : :
1027 : : // loop over all nodes of the element e
1028 [ + + ]: 13340070 : for (std::size_t lp=0; lp<4; ++lp)
1029 : : {
1030 : : // reset min/max
1031 [ + + ]: 45887816 : for (std::size_t c=VarRange[0]; c<=VarRange[1]; ++c)
1032 : : {
1033 : 35215760 : auto mark = c*rdof;
1034 : 35215760 : auto cmark = c-VarRange[0];
1035 [ + - ]: 35215760 : uMin[cmark] = U(e, mark, offset);
1036 [ + - ]: 35215760 : uMax[cmark] = U(e, mark, offset);
1037 : : }
1038 : 10672056 : auto p = inpoel[4*e+lp];
1039 [ + - ]: 10672056 : const auto& pesup = tk::cref_find(esup, p);
1040 : :
1041 : : // ----- Step-1: find min/max in the neighborhood of node p
1042 : : // loop over all the internal elements surrounding this node p
1043 [ + + ]: 176243396 : for (auto er : pesup)
1044 : : {
1045 [ + + ]: 712196460 : for (std::size_t c=VarRange[0]; c<=VarRange[1]; ++c)
1046 : : {
1047 : 546625120 : auto mark = c*rdof;
1048 : 546625120 : auto cmark = c-VarRange[0];
1049 [ + - ]: 546625120 : uMin[cmark] = std::min(uMin[cmark], U(er, mark, offset));
1050 [ + - ]: 546625120 : uMax[cmark] = std::max(uMax[cmark], U(er, mark, offset));
1051 : : }
1052 : : }
1053 : :
1054 : : // ----- Step-2: compute the limiter function at this node
1055 : : // find high-order solution
1056 [ + - ]: 21344112 : std::vector< tk::real > state( ncomp, 0.0 );
1057 [ + + ]: 10672056 : if(rdof == 4)
1058 : : {
1059 : : // If DG(P1), evaluate high order solution based on dubiner basis
1060 : 10657640 : std::array< tk::real, 3 > gp{cx[p], cy[p], cz[p]};
1061 : : auto B_p = tk::eval_basis( rdof,
1062 : 10657640 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
1063 : 10657640 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
1064 [ + - ]: 31972920 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
1065 [ + - ]: 10657640 : state = tk::eval_state( ncomp, offset, rdof, dof_el, e, U, B_p, VarRange );
1066 : : }
1067 : : else { // If DG(P2), evaluate high order solution based on Taylor basis
1068 : : // The nodal and central coordinates
1069 : 14416 : std::array< tk::real, 3 > node{cx[p], cy[p], cz[p]};
1070 : : std::array< tk::real, 3 > x_center
1071 [ + - ][ + - ]: 14416 : { geoElem(e,1,0), geoElem(e,2,0), geoElem(e,3,0) };
[ + - ]
1072 [ + - ]: 28832 : auto B_p = tk::eval_TaylorBasis( rdof, node, x_center, coordel );
1073 : :
1074 [ + + ]: 86496 : for (ncomp_t c=0; c<ncomp; ++c)
1075 [ + + ]: 360400 : for(std::size_t idof = 0; idof < 4; idof++)
1076 : 288320 : state[c] += unk[c][idof] * B_p[idof];
1077 : : }
1078 : :
1079 [ - + ][ - - ]: 10672056 : Assert( state.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
1080 : :
1081 : : // compute the limiter function
1082 [ + + ]: 45887816 : for (std::size_t c=VarRange[0]; c<=VarRange[1]; ++c)
1083 : : {
1084 : 35215760 : auto phi_gp = 1.0;
1085 : 35215760 : auto mark = c*rdof;
1086 [ + - ]: 35215760 : auto uNeg = state[c] - U(e, mark, offset);
1087 [ + - ]: 35215760 : auto uref = std::max(std::fabs(U(e,mark,offset)), 1e-14);
1088 : 35215760 : auto cmark = c - VarRange[0];
1089 [ + + ]: 35215760 : if (uNeg > 1.0e-06*uref)
1090 : : {
1091 [ + - ]: 5847521 : phi_gp = std::min( 1.0, (uMax[cmark]-U(e, mark, offset))/uNeg );
1092 : : }
1093 [ + + ]: 29368239 : else if (uNeg < -1.0e-06*uref)
1094 : : {
1095 [ + - ]: 5932582 : phi_gp = std::min( 1.0, (uMin[cmark]-U(e, mark, offset))/uNeg );
1096 : : }
1097 : : else
1098 : : {
1099 : 23435657 : phi_gp = 1.0;
1100 : : }
1101 : :
1102 : : // ----- Step-3: take the minimum of the nodal-limiter functions
1103 : 35215760 : phi[c] = std::min( phi[c], phi_gp );
1104 : : }
1105 : : }
1106 : 2668014 : }
1107 : :
1108 : : std::vector< tk::real >
1109 : 3604 : VertexBasedLimiting_P2( const std::vector< std::vector< tk::real > >& unk,
1110 : : const tk::Fields& U,
1111 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
1112 : : const std::vector< std::size_t >& inpoel,
1113 : : const tk::UnsMesh::Coords& coord,
1114 : : const tk::Fields& geoElem,
1115 : : std::size_t e,
1116 : : std::size_t rdof,
1117 : : [[maybe_unused]] std::size_t dof_el,
1118 : : std::size_t offset,
1119 : : std::size_t ncomp,
1120 : : const std::vector< std::size_t >& gid,
1121 : : const std::unordered_map< std::size_t, std::size_t >& bid,
1122 : : const std::vector< std::vector<tk::real> >& NodalExtrm )
1123 : : // *****************************************************************************
1124 : : // Kuzmin's vertex-based limiter function calculation for P2 dofs
1125 : : //! \param[in] U High-order solution vector which is to be limited
1126 : : //! \param[in] esup Elements surrounding points
1127 : : //! \param[in] inpoel Element connectivity
1128 : : //! \param[in] coord Array of nodal coordinates
1129 : : //! \param[in] geoElem Element geometry array
1130 : : //! \param[in] e Id of element whose solution is to be limited
1131 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
1132 : : //! \param[in] dof_el Local number of degrees of freedom
1133 : : //! \param[in] offset Index for equation systems
1134 : : //! \param[in] ncomp Number of scalar components in this PDE system
1135 : : //! \param[in] gid Local->global node id map
1136 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
1137 : : //! global node ids (key)
1138 : : //! \param[in] NodalExtrm Chare-boundary nodal extrema
1139 : : //! \return phi Limiter function for solution in element e
1140 : : //! \details This function limits the P2 dofs of P2 solution in a hierachical
1141 : : //! way to P1 dof limiting. Here we treat the first order derivatives the same
1142 : : //! way as cell average while second order derivatives represent the gradients
1143 : : //! to be limited in the P1 limiting procedure.
1144 : : // *****************************************************************************
1145 : : {
1146 : 3604 : const auto nelem = inpoel.size() / 4;
1147 : :
1148 : : // Prepare for calculating Basis functions
1149 : 3604 : const auto& cx = coord[0];
1150 : 3604 : const auto& cy = coord[1];
1151 : 3604 : const auto& cz = coord[2];
1152 : :
1153 [ + - ]: 3604 : std::vector< tk::real > phi(ncomp, 1.0);
1154 : 7208 : std::vector< std::vector< tk::real > > uMin, uMax;
1155 [ + - ][ + - ]: 3604 : uMin.resize( ncomp, std::vector<tk::real>(3, 0.0) );
1156 [ + - ][ + - ]: 3604 : uMax.resize( ncomp, std::vector<tk::real>(3, 0.0) );
1157 : :
1158 : : // The coordinates of centroid in the reference domain
1159 : 7208 : std::array< std::vector< tk::real >, 3 > center;
1160 [ + - ]: 3604 : center[0].resize(1, 0.25);
1161 [ + - ]: 3604 : center[1].resize(1, 0.25);
1162 [ + - ]: 3604 : center[2].resize(1, 0.25);
1163 : :
1164 : : // loop over all nodes of the element e
1165 [ + + ]: 18020 : for (std::size_t lp=0; lp<4; ++lp)
1166 : : {
1167 : : // Find the max/min first-order derivatives for internal element
1168 [ + + ]: 86496 : for (std::size_t c=0; c<ncomp; ++c)
1169 : : {
1170 [ + + ]: 288320 : for (std::size_t idir=1; idir < 4; ++idir)
1171 : : {
1172 : 216240 : uMin[c][idir-1] = unk[c][idir];
1173 : 216240 : uMax[c][idir-1] = unk[c][idir];
1174 : : }
1175 : : }
1176 : :
1177 : 14416 : auto p = inpoel[4*e+lp];
1178 [ + - ]: 14416 : const auto& pesup = tk::cref_find(esup, p);
1179 : :
1180 : : // Step-1: find min/max first order derivative at the centroid in the
1181 : : // neighborhood of node p
1182 [ + + ]: 194892 : for (auto er : pesup)
1183 : : {
1184 [ + - ]: 180476 : if(er < nelem) // If this is internal element
1185 : : {
1186 : : // Coordinates of the neighboring element
1187 : : std::array< std::array< tk::real, 3>, 4 > coorder {{
1188 : 541428 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
1189 : 541428 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
1190 : 541428 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
1191 : 1624284 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
1192 : :
1193 : : auto jacInv_er =
1194 : 180476 : tk::inverseJacobian( coorder[0], coorder[1], coorder[2], coorder[3] );
1195 : :
1196 : : // Compute the derivatives of basis function in the physical domain
1197 [ + - ]: 360952 : auto dBdx_er = tk::eval_dBdx_p1( rdof, jacInv_er );
1198 : :
1199 [ + - ]: 180476 : if(rdof > 4)
1200 [ + - ]: 180476 : tk::eval_dBdx_p2(0, center, jacInv_er, dBdx_er);
1201 : :
1202 [ + + ]: 1082856 : for (std::size_t c=0; c<ncomp; ++c)
1203 : : {
1204 : 902380 : auto mark = c*rdof;
1205 [ + + ]: 3609520 : for (std::size_t idir=0; idir < 3; ++idir)
1206 : : {
1207 : : // The first order derivative at the centroid of element er
1208 : 2707140 : tk::real slope_er(0.0);
1209 [ + + ]: 27071400 : for(std::size_t idof = 1; idof < rdof; idof++)
1210 [ + - ]: 24364260 : slope_er += U(er, mark+idof, offset) * dBdx_er[idir][idof];
1211 : :
1212 : 2707140 : uMin[c][idir] = std::min(uMin[c][idir], slope_er);
1213 : 2707140 : uMax[c][idir] = std::max(uMax[c][idir], slope_er);
1214 : :
1215 : : }
1216 : : }
1217 : : }
1218 : : }
1219 : : // If node p is the chare-boundary node, find min/max by comparing with
1220 : : // the chare-boundary nodal extrema from vector NodalExtrm
1221 [ + - ]: 14416 : auto gip = bid.find( gid[p] );
1222 [ - + ]: 14416 : if(gip != end(bid))
1223 : : {
1224 : 0 : auto ndof_NodalExtrm = NodalExtrm[0].size() / (ncomp * 2);
1225 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
1226 : : {
1227 [ - - ]: 0 : for (std::size_t idir = 0; idir < 3; idir++)
1228 : : {
1229 : 0 : auto max_mark = 2*c*ndof_NodalExtrm + 2*idir;
1230 : 0 : auto min_mark = max_mark + 1;
1231 : 0 : auto& ex = NodalExtrm[gip->second];
1232 : 0 : uMax[c][idir] = std::max(ex[max_mark], uMax[c][idir]);
1233 : 0 : uMin[c][idir] = std::min(ex[min_mark], uMin[c][idir]);
1234 : : }
1235 : : }
1236 : : }
1237 : :
1238 : : //Step-2: compute the limiter function at this node
1239 : 14416 : std::array< tk::real, 3 > node{cx[p], cy[p], cz[p]};
1240 : : std::array< tk::real, 3 >
1241 [ + - ][ + - ]: 14416 : centroid_physical{geoElem(e,1,0), geoElem(e,2,0), geoElem(e,3,0)};
[ + - ]
1242 : :
1243 : : // find high-order solution
1244 : 28832 : std::vector< std::vector< tk::real > > state;
1245 [ + - ][ + - ]: 14416 : state.resize( ncomp, std::vector<tk::real>(3, 0.0) );
1246 : :
1247 [ + + ]: 86496 : for (ncomp_t c=0; c<ncomp; ++c)
1248 : : {
1249 : 72080 : auto dx = node[0] - centroid_physical[0];
1250 : 72080 : auto dy = node[1] - centroid_physical[1];
1251 : 72080 : auto dz = node[2] - centroid_physical[2];
1252 : :
1253 : 72080 : state[c][0] = unk[c][1] + unk[c][4] * dx + unk[c][7] * dy + unk[c][8] * dz;
1254 : 72080 : state[c][1] = unk[c][2] + unk[c][5] * dy + unk[c][7] * dx + unk[c][9] * dz;
1255 : 72080 : state[c][2] = unk[c][3] + unk[c][6] * dz + unk[c][8] * dx + unk[c][9] * dy;
1256 : : }
1257 : :
1258 : : // compute the limiter function
1259 [ + + ]: 86496 : for (std::size_t c=0; c<ncomp; ++c)
1260 : : {
1261 : 72080 : tk::real phi_dir(1.0);
1262 [ + + ]: 216240 : for (std::size_t idir=1; idir < 3; ++idir)
1263 : : {
1264 : 144160 : phi_dir = 1.0;
1265 : 144160 : auto uNeg = state[c][idir-1] - unk[c][idir];
1266 : 144160 : auto uref = std::max(std::fabs(unk[c][idir]), 1e-14);
1267 [ + + ]: 144160 : if (uNeg > 1.0e-6*uref)
1268 : : {
1269 : 67360 : phi_dir =
1270 : 67360 : std::min( 1.0, ( uMax[c][idir-1] - unk[c][idir])/uNeg );
1271 : : }
1272 [ + - ]: 76800 : else if (uNeg < -1.0e-6*uref)
1273 : : {
1274 : 76800 : phi_dir =
1275 : 76800 : std::min( 1.0, ( uMin[c][idir-1] - unk[c][idir])/uNeg );
1276 : : }
1277 : : else
1278 : : {
1279 : 0 : phi_dir = 1.0;
1280 : : }
1281 : :
1282 : 144160 : phi[c] = std::min( phi[c], phi_dir );
1283 : : }
1284 : : }
1285 : : }
1286 : :
1287 : 7208 : return phi;
1288 : : }
1289 : :
1290 : :
1291 : :
1292 : 898870 : void consistentMultiMatLimiting_P1(
1293 : : std::size_t nmat,
1294 : : ncomp_t offset,
1295 : : std::size_t rdof,
1296 : : std::size_t e,
1297 : : tk::Fields& U,
1298 : : [[maybe_unused]] tk::Fields& P,
1299 : : std::vector< tk::real >& phic,
1300 : : [[maybe_unused]] std::vector< tk::real >& phip )
1301 : : // *****************************************************************************
1302 : : // Consistent limiter modifications for P1 dofs
1303 : : //! \param[in] nmat Number of materials in this PDE system
1304 : : //! \param[in] offset Index for equation system
1305 : : //! \param[in] rdof Total number of reconstructed dofs
1306 : : //! \param[in] e Element being checked for consistency
1307 : : //! \param[in,out] U Second-order solution vector which gets modified near
1308 : : //! material interfaces for consistency
1309 : : //! \param[in,out] P Second-order vector of primitive quantities which gets
1310 : : //! modified near material interfaces for consistency
1311 : : //! \param[in,out] phic Vector of limiter functions for the conserved quantities
1312 : : //! \param[in,out] phip Vector of limiter functions for the primitive quantities
1313 : : // *****************************************************************************
1314 : : {
1315 [ - + ][ - - ]: 898870 : Assert(phic.size() == U.nprop()/rdof, "Number of unknowns in vector of "
[ - - ][ - - ]
1316 : : "conserved quantities incorrect");
1317 [ - + ][ - - ]: 898870 : Assert(phip.size() == P.nprop()/rdof, "Number of unknowns in vector of "
[ - - ][ - - ]
1318 : : "primitive quantities incorrect");
1319 : :
1320 : : // find the limiter-function for volume-fractions
1321 : 898870 : auto phi_al(1.0), almax(0.0), dalmax(0.0);
1322 : : //std::size_t nmax(0);
1323 [ + + ]: 2696610 : for (std::size_t k=0; k<nmat; ++k)
1324 : : {
1325 : 1797740 : phi_al = std::min( phi_al, phic[volfracIdx(nmat, k)] );
1326 [ + - ][ + + ]: 1797740 : if (almax < U(e,volfracDofIdx(nmat, k, rdof, 0),offset))
1327 : : {
1328 : : //nmax = k;
1329 [ + - ]: 1253748 : almax = U(e,volfracDofIdx(nmat, k, rdof, 0),offset);
1330 : : }
1331 : 1797740 : auto dmax = std::max(
1332 : : std::max(
1333 [ + - ]: 1797740 : std::abs(U(e,volfracDofIdx(nmat, k, rdof, 1),offset)),
1334 [ + - ]: 1797740 : std::abs(U(e,volfracDofIdx(nmat, k, rdof, 2),offset)) ),
1335 [ + - ]: 3595480 : std::abs(U(e,volfracDofIdx(nmat, k, rdof, 3),offset)) );
1336 : 1797740 : dalmax = std::max( dalmax, dmax );
1337 : : }
1338 : :
1339 : 898870 : auto al_band = 1e-4;
1340 : :
1341 : : //phi_al = phic[nmax];
1342 : :
1343 : : // determine if cell is a material-interface cell based on ad-hoc tolerances.
1344 : : // if interface-cell, then modify high-order dofs of conserved unknowns
1345 : : // consistently and use same limiter for all equations.
1346 : : // Slopes of solution variables \alpha_k \rho_k and \alpha_k \rho_k E_k need
1347 : : // to be modified in interface cells, such that slopes in the \rho_k and
1348 : : // \rho_k E_k part are ignored and only slopes in \alpha_k are considered.
1349 : : // Ideally, we would like to not do this, but this is a necessity to avoid
1350 : : // limiter-limiter interactions in multiphase CFD (see "K.-M. Shyue, F. Xiao,
1351 : : // An Eulerian interface sharpening algorithm for compressible two-phase flow:
1352 : : // the algebraic THINC approach, Journal of Computational Physics 268, 2014,
1353 : : // 326–354. doi:10.1016/j.jcp.2014.03.010." and "A. Chiapolino, R. Saurel,
1354 : : // B. Nkonga, Sharpening diffuse interfaces with compressible fluids on
1355 : : // unstructured meshes, Journal of Computational Physics 340 (2017) 389–417.
1356 : : // doi:10.1016/j.jcp.2017.03.042."). This approximation should be applied in
1357 : : // as narrow a band of interface-cells as possible. The following if-test
1358 : : // defines this band of interface-cells. This tests checks the value of the
1359 : : // maximum volume-fraction in the cell (almax) and the maximum change in
1360 : : // volume-fraction in the cell (dalmax, calculated from second-order DOFs),
1361 : : // to determine the band of interface-cells where the aforementioned fix needs
1362 : : // to be applied. This if-test says that, the fix is applied when the change
1363 : : // in volume-fraction across a cell is greater than 0.1, *and* the
1364 : : // volume-fraction is between 0.1 and 0.9.
1365 [ + + ][ + - ]: 898870 : if ( dalmax > al_band &&
1366 [ + + ]: 71068 : (almax > al_band && almax < (1.0-al_band)) )
1367 : : {
1368 : : // 1. consistent high-order dofs
1369 [ + + ]: 81480 : for (std::size_t k=0; k<nmat; ++k)
1370 : : {
1371 [ + - ]: 54320 : auto alk = std::max( 1.0e-14, U(e,volfracDofIdx(nmat, k, rdof, 0),offset) );
1372 [ + - ]: 54320 : auto rhok = U(e,densityDofIdx(nmat, k, rdof, 0),offset)/alk;
1373 [ + + ]: 217280 : for (std::size_t idir=1; idir<=3; ++idir)
1374 : : {
1375 [ + - ]: 162960 : U(e,densityDofIdx(nmat, k, rdof, idir),offset) = rhok *
1376 [ + - ]: 162960 : U(e,volfracDofIdx(nmat, k, rdof, idir),offset);
1377 : : }
1378 : : }
1379 : :
1380 : : // 2. same limiter for all volume-fractions and densities
1381 [ + + ]: 81480 : for (std::size_t k=0; k<nmat; ++k)
1382 : : {
1383 : 54320 : phic[volfracIdx(nmat, k)] = phi_al;
1384 : 54320 : phic[densityIdx(nmat, k)] = phi_al;
1385 : 27160 : }
1386 : : }
1387 : : else
1388 : : {
1389 : : // same limiter for all volume-fractions
1390 [ + + ]: 2615130 : for (std::size_t k=0; k<nmat; ++k)
1391 : 1743420 : phic[volfracIdx(nmat, k)] = phi_al;
1392 : : }
1393 : 898870 : }
1394 : :
1395 : 227400 : void BoundPreservingLimiting( std::size_t nmat,
1396 : : ncomp_t offset,
1397 : : std::size_t ndof,
1398 : : std::size_t e,
1399 : : const std::vector< std::size_t >& inpoel,
1400 : : const tk::UnsMesh::Coords& coord,
1401 : : const tk::Fields& U,
1402 : : std::vector< tk::real >& phic )
1403 : : // *****************************************************************************
1404 : : // Bound preserving limiter for P1 dofs when MulMat scheme is selected
1405 : : //! \param[in] nmat Number of materials in this PDE system
1406 : : //! \param[in] offset Index for equation system
1407 : : //! \param[in] ndof Total number of reconstructed dofs
1408 : : //! \param[in] e Element being checked for consistency
1409 : : //! \param[in] inpoel Element connectivity
1410 : : //! \param[in] coord Array of nodal coordinates
1411 : : //! \param[in,out] U Second-order solution vector which gets modified near
1412 : : //! material interfaces for consistency
1413 : : //! \param[in,out] phic Vector of limiter functions for the conserved quantities
1414 : : //! \details This bound-preserving limiter is specifically meant to enforce
1415 : : //! bounds [0,1], but it does not suppress oscillations like the other 'TVD'
1416 : : //! limiters. TVD limiters on the other hand, do not preserve such bounds. A
1417 : : //! combination of oscillation-suppressing and bound-preserving limiters can
1418 : : //! obtain a non-oscillatory and bounded solution.
1419 : : // *****************************************************************************
1420 : : {
1421 : 227400 : const auto& cx = coord[0];
1422 : 227400 : const auto& cy = coord[1];
1423 : 227400 : const auto& cz = coord[2];
1424 : :
1425 : : // Extract the element coordinates
1426 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
1427 : 682200 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
1428 : 682200 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
1429 : 682200 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
1430 : 2046600 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
1431 : :
1432 : : // Compute the determinant of Jacobian matrix
1433 : : auto detT =
1434 : 227400 : tk::Jacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
1435 : :
1436 [ + - ]: 454800 : std::vector< tk::real > phi_bound(nmat, 1.0);
1437 : :
1438 : : // loop over all faces of the element e
1439 [ + + ]: 1137000 : for (std::size_t lf=0; lf<4; ++lf)
1440 : : {
1441 : : // Extract the face coordinates
1442 : 909600 : std::array< std::size_t, 3 > inpofa_l {{ inpoel[4*e+tk::lpofa[lf][0]],
1443 : 909600 : inpoel[4*e+tk::lpofa[lf][1]],
1444 : 1819200 : inpoel[4*e+tk::lpofa[lf][2]] }};
1445 : :
1446 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
1447 : 2728800 : {{ cx[ inpofa_l[0] ], cy[ inpofa_l[0] ], cz[ inpofa_l[0] ] }},
1448 : 2728800 : {{ cx[ inpofa_l[1] ], cy[ inpofa_l[1] ], cz[ inpofa_l[1] ] }},
1449 : 5457600 : {{ cx[ inpofa_l[2] ], cy[ inpofa_l[2] ], cz[ inpofa_l[2] ] }} }};
1450 : :
1451 [ + - ]: 909600 : auto ng = tk::NGfa(ndof);
1452 : :
1453 : : // arrays for quadrature points
1454 : 1819200 : std::array< std::vector< tk::real >, 2 > coordgp;
1455 : 1819200 : std::vector< tk::real > wgp;
1456 : :
1457 [ + - ]: 909600 : coordgp[0].resize( ng );
1458 [ + - ]: 909600 : coordgp[1].resize( ng );
1459 [ + - ]: 909600 : wgp.resize( ng );
1460 : :
1461 : : // get quadrature point weights and coordinates for triangle
1462 [ + - ]: 909600 : tk::GaussQuadratureTri( ng, coordgp, wgp );
1463 : :
1464 : : // Compute the upper and lower bound for volume fraction
1465 : 909600 : tk::real min = 1e-14;
1466 : 909600 : tk::real max = 1.0 - min;
1467 : :
1468 : : // Gaussian quadrature
1469 [ + + ]: 3638400 : for (std::size_t igp=0; igp<ng; ++igp)
1470 : : {
1471 : : // Compute the coordinates of quadrature point at physical domain
1472 [ + - ]: 2728800 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
1473 : :
1474 : : //Compute the basis functions
1475 : : auto B = tk::eval_basis( ndof,
1476 : 2728800 : tk::Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT,
1477 : 2728800 : tk::Jacobian( coordel[0], coordel[1], gp, coordel[3] ) / detT,
1478 [ + - ]: 10915200 : tk::Jacobian( coordel[0], coordel[1], coordel[2], gp ) / detT );
1479 : :
1480 : 2728800 : auto state = eval_state( U.nprop()/ndof, offset, ndof, ndof, e, U, B,
1481 [ + - ]: 5457600 : {0, U.nprop()/ndof-1} );
1482 : :
1483 [ + + ]: 8186400 : for(std::size_t imat = 0; imat < nmat; imat++)
1484 : : {
1485 : 5457600 : tk::real phi(1.0);
1486 : 5457600 : auto al = state[volfracIdx(nmat, imat)];
1487 [ + + ]: 5457600 : if(al > 1.0)
1488 : : {
1489 : 255707 : phi = std::fabs(
1490 [ + - ]: 255707 : (max - U(e,volfracDofIdx(nmat, imat, ndof, 0),offset))
1491 [ + - ]: 255707 : / (al - U(e,volfracDofIdx(nmat, imat, ndof, 0),offset)) );
1492 : : }
1493 [ + + ]: 5201893 : else if(al < 1e-14)
1494 : : {
1495 : 5374 : phi = std::fabs(
1496 [ + - ]: 5374 : (min - U(e,volfracDofIdx(nmat, imat, ndof, 0),offset))
1497 [ + - ]: 5374 : / (al - U(e,volfracDofIdx(nmat, imat, ndof, 0),offset)) );
1498 : : }
1499 : :
1500 : 5457600 : phi_bound[imat] = std::min( phi_bound[imat], phi );
1501 : : }
1502 : : }
1503 : : }
1504 : :
1505 [ + + ]: 682200 : for(std::size_t imat = 0; imat < nmat; imat++)
1506 : 454800 : phic[imat] = phi_bound[imat] * phic[imat];
1507 : 227400 : }
1508 : :
1509 : : bool
1510 : 8981850 : interfaceIndicator( std::size_t nmat,
1511 : : const std::vector< tk::real >& al,
1512 : : std::vector< std::size_t >& matInt )
1513 : : // *****************************************************************************
1514 : : // Interface indicator function, which checks element for material interface
1515 : : //! \param[in] nmat Number of materials in this PDE system
1516 : : //! \param[in] al Cell-averaged volume fractions
1517 : : //! \param[in] matInt Array indicating which material has an interface
1518 : : //! \return Boolean which indicates if the element contains a material interface
1519 : : // *****************************************************************************
1520 : : {
1521 : 8981850 : bool intInd = false;
1522 : :
1523 : : // limits under which compression is to be performed
1524 : 8981850 : auto al_eps = 1e-08;
1525 : 8981850 : auto loLim = 2.0 * al_eps;
1526 : 8981850 : auto hiLim = 1.0 - loLim;
1527 : :
1528 : 8981850 : auto almax = 0.0;
1529 [ + + ]: 25306200 : for (std::size_t k=0; k<nmat; ++k)
1530 : : {
1531 : 16324350 : almax = std::max(almax, al[k]);
1532 : 16324350 : matInt[k] = 0;
1533 [ + + ][ + + ]: 16324350 : if ((al[k] > loLim) && (al[k] < hiLim)) matInt[k] = 1;
[ + + ]
1534 : : }
1535 : :
1536 [ + + ][ + + ]: 8981850 : if ((almax > loLim) && (almax < hiLim)) intInd = true;
1537 : :
1538 : 8981850 : return intInd;
1539 : : }
1540 : :
1541 : 0 : void MarkShockCells ( const std::size_t nelem,
1542 : : const std::size_t nmat,
1543 : : const std::size_t system,
1544 : : const std::size_t offset,
1545 : : const std::size_t ndof,
1546 : : const std::size_t rdof,
1547 : : const std::vector< std::size_t >& ndofel,
1548 : : const std::vector< std::size_t >& inpoel,
1549 : : const tk::UnsMesh::Coords& coord,
1550 : : const inciter::FaceData& fd,
1551 : : const tk::Fields& geoFace,
1552 : : const tk::Fields& geoElem,
1553 : : const tk::Fields& U,
1554 : : const tk::Fields& P,
1555 : : std::vector< std::size_t >& shockmarker )
1556 : : // *****************************************************************************
1557 : : // Mark the cells that contain discontinuity according to the interface
1558 : : // condition
1559 : : //! \param[in] nelem Number of elements
1560 : : //! \param[in] nmat Number of materials in this PDE system
1561 : : //! \param[in] system Equation system index
1562 : : //! \param[in] offset Offset this PDE system operates from
1563 : : //! \param[in] ndof Maximum number of degrees of freedom
1564 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
1565 : : //! \param[in] ndofel Vector of local number of degrees of freedome
1566 : : //! \param[in] inpoel Element-node connectivity
1567 : : //! \param[in] coord Array of nodal coordinates
1568 : : //! \param[in] fd Face connectivity and boundary conditions object
1569 : : //! \param[in] geoFace Face geometry array
1570 : : //! \param[in] geoElem Element geometry array
1571 : : //! \param[in] U Solution vector at recent time step
1572 : : //! \param[in] P Vector of primitives at recent time step
1573 : : //! \param[in, out] shockmarker Vector of the shock indicator
1574 : : //! \details This function computes the discontinuity indicator based on
1575 : : //! interface conditon. It is based on the following paper:
1576 : : //! Hong L., Gianni A., Robert N. (2021) A moving discontinuous Galerkin
1577 : : //! finite element method with interface condition enforcement for
1578 : : //! compressible flows. Journal of Computational Physics,
1579 : : //! doi: https://doi.org/10.1016/j.jcp.2021.110618
1580 : : // *****************************************************************************
1581 : : {
1582 [ - - ]: 0 : std::vector< tk::real > IC(U.nunk(), 0.0);
1583 : 0 : const auto& esuf = fd.Esuf();
1584 : 0 : const auto& inpofa = fd.Inpofa();
1585 : :
1586 : 0 : const auto& cx = coord[0];
1587 : 0 : const auto& cy = coord[1];
1588 : 0 : const auto& cz = coord[2];
1589 : :
1590 : 0 : auto ncomp = U.nprop()/rdof;
1591 : 0 : auto nprim = P.nprop()/rdof;
1592 : :
1593 : : // Loop over faces
1594 [ - - ][ - - ]: 0 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f) {
1595 [ - - ][ - - ]: 0 : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
[ - - ][ - - ]
[ - - ]
1596 : : "as -1" );
1597 : :
1598 : 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
1599 : 0 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
1600 : :
1601 : : // When the number of gauss points for the left and right element are
1602 : : // different, choose the larger ng
1603 [ - - ]: 0 : auto ng_l = tk::NGfa(ndofel[el]);
1604 [ - - ]: 0 : auto ng_r = tk::NGfa(ndofel[er]);
1605 : :
1606 : 0 : auto ng = std::max( ng_l, ng_r );
1607 : :
1608 : : std::array< std::vector< tk::real >, 2 > coordgp
1609 [ - - ][ - - ]: 0 : { std::vector<tk::real>(ng), std::vector<tk::real>(ng) };
1610 [ - - ]: 0 : std::vector< tk::real > wgp( ng );
1611 : :
1612 [ - - ]: 0 : tk::GaussQuadratureTri( ng, coordgp, wgp );
1613 : :
1614 : : // Extract the element coordinates
1615 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
1616 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
1617 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
1618 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
1619 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
1620 : :
1621 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
1622 : 0 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
1623 : 0 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
1624 : 0 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
1625 : 0 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
1626 : :
1627 : : // Compute the determinant of Jacobian matrix
1628 : : auto detT_l =
1629 : 0 : tk::Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
1630 : : auto detT_r =
1631 : 0 : tk::Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
1632 : :
1633 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
1634 : 0 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
1635 : 0 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
1636 : 0 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
1637 : :
1638 : : std::array< tk::real, 3 >
1639 [ - - ][ - - ]: 0 : fn{{ geoFace(f,1,0), geoFace(f,2,0), geoFace(f,3,0) }};
[ - - ]
1640 : :
1641 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp) {
1642 [ - - ]: 0 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
1643 : : std::size_t dof_el, dof_er;
1644 [ - - ]: 0 : if (rdof > ndof)
1645 : : {
1646 : 0 : dof_el = rdof;
1647 : 0 : dof_er = rdof;
1648 : : }
1649 : : else
1650 : : {
1651 : 0 : dof_el = ndofel[el];
1652 : 0 : dof_er = ndofel[er];
1653 : : }
1654 : : std::array< tk::real, 3> ref_gp_l{
1655 : 0 : tk::Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
1656 : 0 : tk::Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
1657 : 0 : tk::Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
1658 : : std::array< tk::real, 3> ref_gp_r{
1659 : 0 : tk::Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
1660 : 0 : tk::Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
1661 : 0 : tk::Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
1662 [ - - ]: 0 : auto B_l = tk::eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
1663 [ - - ]: 0 : auto B_r = tk::eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
1664 : :
1665 [ - - ]: 0 : auto wt = wgp[igp] * geoFace(f,0,0);
1666 : :
1667 : 0 : std::array< std::vector< tk::real >, 2 > state;
1668 : :
1669 : : // Evaluate the high order solution at the qudrature point
1670 [ - - ]: 0 : state[0] = tk::evalPolynomialSol(system, offset, 0, ncomp, nprim, rdof,
1671 : 0 : nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
1672 [ - - ]: 0 : state[1] = tk::evalPolynomialSol(system, offset, 0, ncomp, nprim, rdof,
1673 : 0 : nmat, er, dof_er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P);
1674 : :
1675 [ - - ][ - - ]: 0 : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
1676 : : "appended boundary state vector" );
1677 [ - - ][ - - ]: 0 : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
1678 : : "appended boundary state vector" );
1679 : :
1680 : : // Evaluate the bulk density
1681 : 0 : tk::real rhol(0.0), rhor(0.0);
1682 [ - - ]: 0 : for(std::size_t k = 0; k < nmat; k++) {
1683 : 0 : rhol += state[0][densityIdx(nmat,k)];
1684 : 0 : rhor += state[1][densityIdx(nmat,k)];
1685 : : }
1686 : :
1687 : : // Evaluate the flux for the density
1688 : 0 : tk::real fl(0.0), fr(0.0);
1689 [ - - ]: 0 : for(std::size_t i = 0; i < 3; i++) {
1690 : 0 : fl += rhol * state[0][ncomp+velocityIdx(nmat,i)] * fn[i];
1691 : 0 : fr += rhor * state[1][ncomp+velocityIdx(nmat,i)] * fn[i];
1692 : : }
1693 : :
1694 : 0 : tk::real rhs = wt * fabs(fl - fr);
1695 : 0 : IC[el] += rhs;
1696 : 0 : IC[er] += rhs;
1697 : : }
1698 : : }
1699 : :
1700 : : // Loop over element to mark shock cell
1701 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e) {
1702 [ - - ]: 0 : if(fabs(IC[e]) > 1e-6)
1703 : 0 : shockmarker[e] = 1;
1704 : : else
1705 : 0 : shockmarker[e] = 0;
1706 : : }
1707 : 0 : }
1708 : :
1709 : : } // inciter::
|
