Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/PrefIndicator.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 Adaptive indicators for p-adaptive discontiunous Galerkin methods
9 : : \details This file contains functions that provide adaptive indicator
10 : : function calculations for marking the number of degree of freedom of each
11 : : element.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include "PrefIndicator.hpp"
16 : :
17 : : #include "Tags.hpp"
18 : : #include "Vector.hpp"
19 : : #include "Integrate/Basis.hpp"
20 : : #include "Integrate/Quadrature.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : 1636 : void spectral_decay( std::size_t nmat,
25 : : std::size_t nunk,
26 : : const std::vector< int >& esuel,
27 : : const tk::Fields& unk,
28 : : const tk::Fields& prim,
29 : : std::size_t ndof,
30 : : std::size_t ndofmax,
31 : : tk::real tolref,
32 : : std::vector< std::size_t >& ndofel )
33 : : // *****************************************************************************
34 : : //! Evaluate the spectral-decay indicator and mark the ndof for each element
35 : : //! \param[in] nmat Number of materials in this PDE system
36 : : //! \param[in] nunk Number of unknowns
37 : : //! \param[in] esuel Elements surrounding elements
38 : : //! \param[in] unk Array of unknowns
39 : : //! \param[in] prim Array of primitive quantities
40 : : //! \param[in] ndof Number of degrees of freedom in the solution
41 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
42 : : //! \param[in] tolref Tolerance for p-refinement
43 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
44 : : //! \details The spectral decay indicator, implemented in this functiopn,
45 : : //! calculates the difference between the projections of the numerical
46 : : //! solutions on finite element space of order p and p-1.
47 : : //! \see F. Naddei, et. al., "A comparison of refinement indicators for the
48 : : //! p-adaptive simulation of steady and unsteady flows with discontinuous
49 : : //! Galerkin methods" at https://doi.org/10.1016/j.jcp.2018.09.045, and G.
50 : : //! Gassner, et al., "Explicit discontinuous Galerkin schemes with adaptation
51 : : //! in space and time"
52 : : // *****************************************************************************
53 : : {
54 : 1636 : const auto ncomp = unk.nprop() / ndof;
55 : 1636 : const auto nprim = prim.nprop() / ndof;
56 : :
57 : : // The array storing the adaptive indicator for each elements
58 : 1636 : std::vector< tk::real > Ind(nunk, 0);
59 : :
60 [ + + ]: 379148 : for (std::size_t e=0; e<esuel.size()/4; ++e)
61 [ + + ]: 377512 : if(ndofel[e] > 1) {
62 [ + - ]: 125871 : if(nmat == 1)
63 : 125871 : Ind[e] =
64 [ + - ]: 125871 : evalDiscIndicator_CompFlow(e, ncomp, ndof, ndofel[e], unk);
65 : : else
66 [ - - ]: 0 : Ind[e] = evalDiscIndicator_MultiMat(e, nmat, ncomp, nprim, ndof,
67 : : ndofel[e], unk, prim);
68 : : }
69 : :
70 : : // As for spectral-decay indicator, rho_p - rho_(p-1) actually is the leading
71 : : // term of discretization error for the numerical solution of p-1. Therefore,
72 : : // this function represents the qualitative behavior of the discretization
73 : : // error. If the value is less than epsL which means the discretization error
74 : : // is already a really small number, then the element should be de-refined. On
75 : : // the other hand, if the value is larger than epsH which means the
76 : : // discretization error is relatively large, then it should be refined.
77 : :
78 : : // Note: Spectral-decay indicator is a measurement of the continuity of the
79 : : // numerical solution inside this element. So when this indicator appears
80 : : // to be relatively large, there might be a shock inside this element and a
81 : : // derefinement or h-refinement should be applied. This condition will be
82 : : // implemented later.
83 : :
84 : : // As for the discretiazation-error based indicator, like spectral-decay
85 : : // indicator, the choices for epsH and epsL are based on the data from
86 : : // numerical experiments. Empirically, we found that when the epsH belongs
87 : : // to [-4, -8] and epsL belongs to [-6, -14], decent results could be
88 : : // observed. And then a linear projection is performed to map epsL and espH
89 : : // to the range of [0, 1] so that it could be controlled by tolref.
90 : :
91 : 1636 : auto epsH = std::pow(10, -4 - tolref * 4.0);
92 : 1636 : auto epsL = std::pow(10, -6 - tolref * 8.0);
93 : :
94 : : // Marke the ndof according to the adaptive indicator
95 [ + + ]: 379148 : for (std::size_t e=0; e<esuel.size()/4; ++e)
96 : : {
97 [ + + ]: 377512 : if(Ind[e] < epsL) // Derefinement
98 : : {
99 [ + + ]: 329584 : if(ndofel[e] == 4)
100 : 43049 : ndofel[e] = 1;
101 [ + + ]: 286535 : else if(ndofel[e] == 10)
102 : 34894 : ndofel[e] = 4;
103 : : }
104 [ + + ]: 47928 : else if(Ind[e] > epsH) // Refinement
105 : : {
106 [ + + ][ - + ]: 18590 : if(ndofel[e] == 4 && ndofmax > 4)
107 : 0 : ndofel[e] = 10;
108 : : }
109 : : }
110 : 1636 : }
111 : :
112 : 0 : void non_conformity( std::size_t nunk,
113 : : std::size_t Nbfac,
114 : : const std::vector< std::size_t >& inpoel,
115 : : const tk::UnsMesh::Coords& coord,
116 : : const std::vector< int >& esuel,
117 : : const std::vector< int >& esuf,
118 : : const std::vector< std::size_t >& inpofa,
119 : : const tk::Fields& unk,
120 : : std::size_t ndof,
121 : : std::size_t ndofmax,
122 : : std::vector< std::size_t >& ndofel )
123 : : // *****************************************************************************
124 : : //! Evaluate the non-conformity indicator and mark the ndof for each element
125 : : //! \param[in] nunk Number of unknowns
126 : : //! \param[in] Nbfac Number of internal faces
127 : : //! \param[in] inpoel Element-node connectivity
128 : : //! \param[in] coord Array of nodal coordinates
129 : : //! \param[in] esuel Elements surrounding elements
130 : : //! \param[in] esuf Elements surrounding faces
131 : : //! \param[in] inpofa Face-node connectivity
132 : : //! \param[in] unk Array of unknowns
133 : : //! \param[in] ndof Number of degrees of freedom in the solution
134 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
135 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
136 : : //! \details The non-conformity indicator, this function implements, evaluates
137 : : //! the jump in the numerical solutions as a measure of the numerical error.
138 : : //! \see F. Naddei, et. al., "A comparison of refinement indicators for the
139 : : //! p-adaptive simulation of steady and unsteady flows with discontinuous
140 : : //! Galerkin methods at https://doi.org/10.1016/j.jcp.2018.09.045.
141 : : //! \warning This indicator can only be applied in serial, i.e., single CPU, for
142 : : //! now because the solution communication happens before eval_ndof() in DG,
143 : : //! which will lead to incorrect evaluation of the numerical solution at the
144 : : //! neighboring cells.
145 : : // *****************************************************************************
146 : : {
147 : 0 : const auto ncomp = unk.nprop() / ndof;
148 : :
149 : : const auto& cx = coord[0];
150 : : const auto& cy = coord[1];
151 : : const auto& cz = coord[2];
152 : :
153 : : // The array storing the adaptive indicator for each elements
154 : 0 : std::vector< tk::real > Ind(nunk, 0);
155 : :
156 : : // compute error indicator for each face
157 [ - - ]: 0 : for (auto f=Nbfac; f<esuf.size()/2; ++f)
158 : : {
159 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
160 : : "as -1" );
161 : :
162 [ - - ]: 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
163 : 0 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
164 : :
165 [ - - ]: 0 : auto ng_l = tk::NGfa(ndofel[el]);
166 [ - - ][ - - ]: 0 : auto ng_r = tk::NGfa(ndofel[er]);
167 : :
168 : : // When the number of gauss points for the left and right element are
169 : : // different, choose the larger ng
170 : 0 : auto ng = std::max( ng_l, ng_r );
171 : :
172 : : // arrays for quadrature points
173 : : std::array< std::vector< tk::real >, 2 > coordgp;
174 : : std::vector< tk::real > wgp;
175 : :
176 [ - - ]: 0 : coordgp[0].resize( ng );
177 [ - - ]: 0 : coordgp[1].resize( ng );
178 [ - - ]: 0 : wgp.resize( ng );
179 : :
180 : : // get quadrature point weights and coordinates for triangle
181 [ - - ]: 0 : tk::GaussQuadratureTri( ng, coordgp, wgp );
182 : :
183 : : // Extract the element coordinates
184 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
185 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
186 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
187 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
188 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
189 : :
190 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
191 : 0 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
192 : 0 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
193 : 0 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
194 : 0 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
195 : :
196 : : // Compute the determinant of Jacobian matrix
197 : : auto detT_l =
198 : 0 : tk::Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
199 : : auto detT_r =
200 : 0 : tk::Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
201 : :
202 : : // Extract the face coordinates
203 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
204 : 0 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
205 : 0 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
206 : 0 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
207 : :
208 : : // Gaussian quadrature
209 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
210 : : {
211 : : // Compute the coordinates of quadrature point at physical domain
212 [ - - ]: 0 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
213 : :
214 : : //Compute the basis functions
215 : 0 : auto B_l = tk::eval_basis( ndofel[el],
216 : 0 : tk::Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
217 : 0 : tk::Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
218 [ - - ]: 0 : tk::Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l );
219 : 0 : auto B_r = tk::eval_basis( ndofel[er],
220 : 0 : tk::Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
221 : 0 : tk::Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
222 [ - - ]: 0 : tk::Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r );
223 : :
224 : : std::array< std::vector< tk::real >, 2 > state;
225 : :
226 [ - - ]: 0 : state[0] = tk::eval_state( ncomp, ndof, ndofel[el], el, unk, B_l );
227 [ - - ]: 0 : state[1] = tk::eval_state( ncomp, ndof, ndofel[er], er, unk, B_r );
228 : :
229 : : Assert( unk[0].size() == ncomp, "Size mismatch" );
230 : : Assert( unk[1].size() == ncomp, "Size mismatch" );
231 : :
232 : 0 : auto rhoL = state[0][0];
233 : 0 : auto rhoR = state[1][0];
234 : :
235 : 0 : auto ind = fabs( rhoL - rhoR ) / 2.0 * ( rhoL + rhoR );
236 [ - - ][ - - ]: 0 : Ind[el] = std::max( ind, Ind[el] );
237 : 0 : Ind[er] = std::max( ind, Ind[er] );
238 : : }
239 : : }
240 : :
241 : : // By assuming a smooth solution, we use the non-conformity indicator to
242 : : // represent the error for the numerical solution qualitatively. If the value
243 : : // is less than epsL which means the error is already a really small number,
244 : : // then the element should be de-refined. On the other hand, if the value is
245 : : // larger than epsH which means the error is relatively large, then it should
246 : : // be refined.
247 : :
248 : : // Marke the ndof according to the adaptive indicator
249 [ - - ]: 0 : for (std::size_t e=0; e<esuel.size()/4; ++e)
250 : : {
251 [ - - ]: 0 : if(Ind[e] < 1e-4) // Derefinement
252 : : {
253 [ - - ]: 0 : if(ndofel[e] == 10)
254 : 0 : ndofel[e] = 4;
255 [ - - ]: 0 : else if(ndofel[e] == 4)
256 : 0 : ndofel[e] = 1;
257 : : }
258 [ - - ]: 0 : else if(Ind[e] > 1e-3) // Refinement
259 : : {
260 [ - - ][ - - ]: 0 : if(ndofel[e] == 4 && ndofmax > 4)
261 : 0 : ndofel[e] = 10;
262 [ - - ]: 0 : else if(ndofel[e] == 1)
263 : 0 : ndofel[e] = 4;
264 : : }
265 : : }
266 : 0 : }
267 : :
268 : 125871 : tk::real evalDiscIndicator_CompFlow( std::size_t e,
269 : : ncomp_t ncomp,
270 : : const std::size_t ndof,
271 : : const std::size_t ndofel,
272 : : const tk::Fields& unk )
273 : : // *****************************************************************************
274 : : //! Evaluate the spectral decay indicator
275 : : //! \param[in] e Index for the tetrahedron element
276 : : //! \param[in] ncomp Number of scalar components in this PDE system
277 : : //! \param[in] ndof Number of degrees of freedom in the solution
278 : : //! \param[in] ndofel Local number of degrees of freedom
279 : : //! \param[in] unk Array of unknowns
280 : : //! \return The value of spectral indicator for the element
281 : : // *****************************************************************************
282 : : {
283 : 125871 : auto ng = tk::NGvol(ndofel);
284 : :
285 : : // arrays for quadrature points
286 : : std::array< std::vector< tk::real >, 3 > coordgp;
287 [ + - ]: 125871 : std::vector< tk::real > wgp( ng );
288 : :
289 [ + - ]: 125871 : coordgp[0].resize( ng );
290 [ + - ]: 125871 : coordgp[1].resize( ng );
291 [ + - ]: 125871 : coordgp[2].resize( ng );
292 : :
293 [ + - ]: 125871 : tk::GaussQuadratureTet( ng, coordgp, wgp );
294 : :
295 : : tk::real dU(0.0), U(0.0), Ind(0.0);
296 : :
297 : : // Gaussian quadrature
298 [ + + ]: 1245876 : for (std::size_t igp=0; igp<ng; ++igp)
299 : : {
300 : : // Compute the basis function
301 [ + - ]: 1120005 : auto B = tk::eval_basis( ndofel, coordgp[0][igp], coordgp[1][igp],
302 [ + - ]: 1120005 : coordgp[2][igp] );
303 : :
304 [ + - ]: 1120005 : auto state = tk::eval_state( ncomp, ndof, ndofel, e, unk, B );
305 : :
306 [ + + ]: 1120005 : U += wgp[igp] * state[0] * state[0];
307 : :
308 [ + + ]: 1120005 : if(ndofel > 4)
309 : : {
310 : 899525 : auto dU_p2 = unk(e, 4) * B[4]
311 : 899525 : + unk(e, 5) * B[5]
312 : 899525 : + unk(e, 6) * B[6]
313 : 899525 : + unk(e, 7) * B[7]
314 : 899525 : + unk(e, 8) * B[8]
315 : 899525 : + unk(e, 9) * B[9];
316 : :
317 : 899525 : dU += wgp[igp] * dU_p2 * dU_p2;
318 : : }
319 : : else
320 : : {
321 : 220480 : auto dU_p1 = unk(e, 1) * B[1]
322 : 220480 : + unk(e, 2) * B[2]
323 : 220480 : + unk(e, 3) * B[3];
324 : :
325 : 220480 : dU += wgp[igp] * dU_p1 * dU_p1;
326 : : }
327 : : }
328 : :
329 [ + - ]: 125871 : Ind = dU / U;
330 : :
331 : 125871 : return Ind;
332 : : }
333 : :
334 : 0 : tk::real evalDiscIndicator_MultiMat( std::size_t e,
335 : : std::size_t nmat,
336 : : ncomp_t ncomp,
337 : : ncomp_t nprim,
338 : : const std::size_t ndof,
339 : : const std::size_t ndofel,
340 : : const tk::Fields& unk,
341 : : const tk::Fields& prim )
342 : : // *****************************************************************************
343 : : //! Evaluate the spectral decay indicator
344 : : //! \param[in] e Index for the tetrahedron element
345 : : //! \param[in] nmat Number of materials in this PDE system
346 : : //! \param[in] ncomp Number of scalar components in this PDE system
347 : : //! \param[in] nprim Number of primitive quantities stored for this PDE system
348 : : //! \param[in] ndof Number of degrees of freedom in the solution
349 : : //! \param[in] ndofel Local number of degrees of freedom
350 : : //! \param[in] unk Array of unknowns
351 : : //! \param[in] prim Array of primitive quantities
352 : : //! \return The value of spectral indicator for the element
353 : : // *****************************************************************************
354 : : {
355 : 0 : auto ng = tk::NGvol(ndof);
356 : :
357 : : // arrays for quadrature points
358 : : std::array< std::vector< tk::real >, 3 > coordgp;
359 [ - - ]: 0 : std::vector< tk::real > wgp( ng );
360 : :
361 [ - - ]: 0 : coordgp[0].resize( ng );
362 [ - - ]: 0 : coordgp[1].resize( ng );
363 [ - - ]: 0 : coordgp[2].resize( ng );
364 : :
365 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
366 : :
367 [ - - ][ - - ]: 0 : std::vector<tk::real> dU(nmat,0.0);
368 [ - - ][ - - ]: 0 : std::vector<tk::real> U(nmat,0.0);
369 [ - - ][ - - ]: 0 : std::vector<std::size_t> marker(nmat,0);
370 : 0 : std::array<tk::real, 3> V{0, 0, 0}, dV{0, 0, 0};
371 : :
372 : : // Gaussian quadrature
373 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
374 : : {
375 : : // Compute the basis function
376 [ - - ]: 0 : auto B = tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp],
377 [ - - ]: 0 : coordgp[2][igp] );
378 : :
379 [ - - ]: 0 : auto state = tk::eval_state( ncomp, ndof, ndofel, e, unk, B );
380 [ - - ]: 0 : auto pstate = tk::eval_state( nprim, ndof, ndofel, e, unk, B );
381 : :
382 : : // indicator based on material density
383 [ - - ]: 0 : for(std::size_t k = 0; k < nmat; k++) {
384 [ - - ]: 0 : if(unk(e, volfracDofIdx(nmat, k, ndof, 0)) > 1e-2) {
385 [ - - ]: 0 : marker[k] = 1;
386 [ - - ]: 0 : U[k] += wgp[igp] *
387 : 0 : state[densityIdx(nmat, k)] * state[densityIdx(nmat, k)];
388 : :
389 [ - - ]: 0 : if(ndofel > 4)
390 : : {
391 : 0 : auto dU_p2 = unk(e, densityDofIdx(nmat, k, ndof, 4)) * B[4]
392 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 5)) * B[5]
393 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 6)) * B[6]
394 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 7)) * B[7]
395 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 8)) * B[8]
396 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 9)) * B[9];
397 : :
398 : 0 : dU[k] += wgp[igp] * dU_p2 * dU_p2;
399 : : }
400 : : else
401 : : {
402 : 0 : auto dU_p1 = unk(e, densityDofIdx(nmat, k, ndof, 1)) * B[1]
403 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 2)) * B[2]
404 : 0 : + unk(e, densityDofIdx(nmat, k, ndof, 3)) * B[3];
405 : :
406 : 0 : dU[k] += wgp[igp] * dU_p1 * dU_p1;
407 : : }
408 : : }
409 : : }
410 : :
411 : : // indicator based on velocity
412 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
413 [ - - ]: 0 : V[i] += wgp[igp] * pstate[velocityIdx(nmat,i)] *
414 : : pstate[velocityIdx(nmat,i)];
415 : :
416 [ - - ]: 0 : if(ndofel > 4)
417 : : {
418 : 0 : auto dV_p2 = prim(e, velocityDofIdx(nmat, i, ndof, 4)) * B[4]
419 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 5)) * B[5]
420 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 6)) * B[6]
421 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 7)) * B[7]
422 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 8)) * B[8]
423 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 9)) * B[9];
424 : :
425 : 0 : dV[i] += wgp[igp] * dV_p2 * dV_p2;
426 : : }
427 : : else
428 : : {
429 : 0 : auto dV_p1 = prim(e, velocityDofIdx(nmat, i, ndof, 1)) * B[1]
430 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 2)) * B[2]
431 : 0 : + prim(e, velocityDofIdx(nmat, i, ndof, 3)) * B[3];
432 : :
433 : 0 : dV[i] += wgp[igp] * dV_p1 * dV_p1;
434 : : }
435 : : }
436 : : }
437 : :
438 : : // The max indicator value among all the materials
439 : 0 : tk::real Indmax = 0.0;
440 [ - - ]: 0 : for(std::size_t k = 0; k < nmat; k++)
441 [ - - ][ - - ]: 0 : if(marker[k]) Indmax = std::max(dU[k]/U[k], Indmax);
442 : :
443 : 0 : std::array<tk::real, 3> vavg{prim(e, velocityDofIdx(nmat,0,ndof,0)),
444 : 0 : prim(e, velocityDofIdx(nmat,1,ndof,0)),
445 : 0 : prim(e, velocityDofIdx(nmat,2,ndof,0))};
446 : 0 : auto vmag = std::sqrt(tk::dot(vavg, vavg));
447 : :
448 : : // The max indicator value among all the velocity components
449 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
450 [ - - ][ - - ]: 0 : if (std::abs(V[i]) > 1e-4*std::max(1.0,vmag))
451 [ - - ]: 0 : Indmax = std::max((dV[i])/(V[i]), Indmax);
452 : : }
453 : :
454 [ - - ]: 0 : return Indmax;
455 : : }
456 : :
457 : : }
458 : : // inciter::
|
