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