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 : : #include "Inciter/InputDeck/InputDeck.hpp"
22 : :
23 : : namespace inciter {
24 : :
25 : : extern ctr::InputDeck g_inputdeck;
26 : :
27 : 1636 : void spectral_decay( std::size_t nmat,
28 : : std::size_t nunk,
29 : : const std::vector< int >& esuel,
30 : : const tk::Fields& unk,
31 : : std::size_t ndof,
32 : : std::size_t ndofmax,
33 : : tk::real tolref,
34 : : std::vector< std::size_t >& ndofel )
35 : : // *****************************************************************************
36 : : //! Evaluate the spectral-decay indicator and mark the ndof for each element
37 : : //! \param[in] nmat Number of materials in this PDE system
38 : : //! \param[in] nunk Number of unknowns
39 : : //! \param[in] esuel Elements surrounding elements
40 : : //! \param[in] unk Array of unknowns
41 : : //! \param[in] ndof Number of degrees of freedom in the solution
42 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
43 : : //! \param[in] tolref Tolerance for p-refinement
44 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
45 : : //! \details The spectral decay indicator, implemented in this functiopn,
46 : : //! calculates the difference between the projections of the numerical
47 : : //! solutions on finite element space of order p and p-1.
48 : : //! \see F. Naddei, et. al., "A comparison of refinement indicators for the
49 : : //! p-adaptive simulation of steady and unsteady flows with discontinuous
50 : : //! Galerkin methods" at https://doi.org/10.1016/j.jcp.2018.09.045, and G.
51 : : //! Gassner, et al., "Explicit discontinuous Galerkin schemes with adaptation
52 : : //! in space and time"
53 : : // *****************************************************************************
54 : : {
55 : 1636 : const auto ncomp = unk.nprop() / ndof;
56 : :
57 : : // The array storing the adaptive indicator for each elements
58 [ + - ]: 3272 : 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 [ + - ]: 100262 : if(nmat == 1)
63 : 100262 : Ind[e] =
64 [ + - ]: 100262 : evalDiscIndicator_CompFlow(e, ncomp, ndof, ndofel[e], unk);
65 : : else
66 : 0 : Ind[e] =
67 [ - - ]: 0 : evalDiscIndicator_MultiMat(e, nmat, ncomp, ndof, ndofel[e], unk);
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 [ + + ]: 326120 : if(ndofel[e] == 4)
100 : 33160 : ndofel[e] = 1;
101 [ + + ]: 292960 : else if(ndofel[e] == 10)
102 : 15710 : ndofel[e] = 4;
103 : : }
104 [ + + ]: 51392 : else if(Ind[e] > epsH) // Refinement
105 : : {
106 [ + + ][ - + ]: 21086 : 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 : 0 : const auto& cx = coord[0];
150 : 0 : const auto& cy = coord[1];
151 : 0 : 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 [ - - ][ - - ]: 0 : 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 : 0 : std::array< std::vector< tk::real >, 2 > coordgp;
174 : 0 : 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 : 0 : std::array< std::vector< tk::real >, 2 > state;
225 : :
226 [ - - ]: 0 : state[0] = tk::eval_state( ncomp, 0, ndof, ndofel[el], el, unk, B_l,
227 : 0 : {0, ncomp-1} );
228 [ - - ]: 0 : state[1] = tk::eval_state( ncomp, 0, ndof, ndofel[er], er, unk, B_r,
229 : 0 : {0, ncomp-1} );
230 : :
231 [ - - ][ - - ]: 0 : Assert( unk[0].size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
[ - - ]
232 [ - - ][ - - ]: 0 : Assert( unk[1].size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
[ - - ]
233 : :
234 : 0 : auto rhoL = state[0][0];
235 : 0 : auto rhoR = state[1][0];
236 : :
237 : 0 : auto ind = fabs( rhoL - rhoR ) / 2.0 * ( rhoL + rhoR );
238 : 0 : Ind[el] = std::max( ind, Ind[el] );
239 : 0 : Ind[er] = std::max( ind, Ind[er] );
240 : : }
241 : : }
242 : :
243 : : // By assuming a smooth solution, we use the non-conformity indicator to
244 : : // represent the error for the numerical solution qualitatively. If the value
245 : : // is less than epsL which means the error is already a really small number,
246 : : // then the element should be de-refined. On the other hand, if the value is
247 : : // larger than epsH which means the error is relatively large, then it should
248 : : // be refined.
249 : :
250 : : // Marke the ndof according to the adaptive indicator
251 [ - - ]: 0 : for (std::size_t e=0; e<esuel.size()/4; ++e)
252 : : {
253 [ - - ]: 0 : if(Ind[e] < 1e-4) // Derefinement
254 : : {
255 [ - - ]: 0 : if(ndofel[e] == 10)
256 : 0 : ndofel[e] = 4;
257 [ - - ]: 0 : else if(ndofel[e] == 4)
258 : 0 : ndofel[e] = 1;
259 : : }
260 [ - - ]: 0 : else if(Ind[e] > 1e-3) // Refinement
261 : : {
262 [ - - ][ - - ]: 0 : if(ndofel[e] == 4 && ndofmax > 4)
[ - - ]
263 : 0 : ndofel[e] = 10;
264 [ - - ]: 0 : else if(ndofel[e] == 1)
265 : 0 : ndofel[e] = 4;
266 : : }
267 : : }
268 : 0 : }
269 : :
270 : 537422 : tk::real evalDiscIndicator_CompFlow( std::size_t e,
271 : : ncomp_t ncomp,
272 : : const std::size_t ndof,
273 : : const std::size_t ndofel,
274 : : const tk::Fields& unk )
275 : : // *****************************************************************************
276 : : //! Evaluate the spectral decay indicator
277 : : //! \param[in] e Index for the tetrahedron element
278 : : //! \param[in] ncomp Number of scalar components in this PDE system
279 : : //! \param[in] ndof Number of degrees of freedom in the solution
280 : : //! \param[in] ndofel Local number of degrees of freedom
281 : : //! \param[in] unk Array of unknowns
282 : : //! \return The value of spectral indicator for the element
283 : : // *****************************************************************************
284 : : {
285 [ + - ]: 537422 : auto ng = tk::NGvol(ndofel);
286 : :
287 : : // arrays for quadrature points
288 : 1074844 : std::array< std::vector< tk::real >, 3 > coordgp;
289 [ + - ]: 537422 : std::vector< tk::real > wgp( ng );
290 : :
291 [ + - ]: 537422 : coordgp[0].resize( ng );
292 [ + - ]: 537422 : coordgp[1].resize( ng );
293 [ + - ]: 537422 : coordgp[2].resize( ng );
294 : :
295 [ + - ]: 537422 : tk::GaussQuadratureTet( ng, coordgp, wgp );
296 : :
297 : 537422 : tk::real dU(0.0), U(0.0), Ind(0.0);
298 : :
299 : : // Gaussian quadrature
300 [ + + ]: 6202302 : for (std::size_t igp=0; igp<ng; ++igp)
301 : : {
302 : : // Compute the basis function
303 : 11329760 : auto B = tk::eval_basis( ndofel, coordgp[0][igp], coordgp[1][igp],
304 [ + - ]: 22659520 : coordgp[2][igp] );
305 : :
306 [ + - ]: 11329760 : auto state = tk::eval_state( ncomp, 0, ndof, ndofel, e, unk, B, {0, ncomp-1} );
307 : :
308 : 5664880 : U += wgp[igp] * state[0] * state[0];
309 : :
310 [ + + ]: 5664880 : if(ndofel > 4)
311 : : {
312 [ + - ]: 5459245 : auto dU_p2 = unk(e, 4, 0) * B[4]
313 [ + - ]: 5459245 : + unk(e, 5, 0) * B[5]
314 [ + - ]: 5459245 : + unk(e, 6, 0) * B[6]
315 [ + - ]: 5459245 : + unk(e, 7, 0) * B[7]
316 [ + - ]: 5459245 : + unk(e, 8, 0) * B[8]
317 [ + - ]: 5459245 : + unk(e, 9, 0) * B[9];
318 : :
319 : 5459245 : dU += wgp[igp] * dU_p2 * dU_p2;
320 : : }
321 : : else
322 : : {
323 [ + - ]: 205635 : auto dU_p1 = unk(e, 1, 0) * B[1]
324 [ + - ]: 205635 : + unk(e, 2, 0) * B[2]
325 [ + - ]: 205635 : + unk(e, 3, 0) * B[3];
326 : :
327 : 205635 : dU += wgp[igp] * dU_p1 * dU_p1;
328 : : }
329 : : }
330 : :
331 : 537422 : Ind = dU / U;
332 : :
333 : 1074844 : return Ind;
334 : : }
335 : :
336 : 454800 : tk::real evalDiscIndicator_MultiMat( std::size_t e,
337 : : std::size_t nmat,
338 : : ncomp_t ncomp,
339 : : const std::size_t ndof,
340 : : const std::size_t ndofel,
341 : : const tk::Fields& unk )
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] ndof Number of degrees of freedom in the solution
348 : : //! \param[in] ndofel Local number of degrees of freedom
349 : : //! \param[in] unk Array of unknowns
350 : : //! \return The value of spectral indicator for the element
351 : : // *****************************************************************************
352 : : {
353 [ + - ]: 454800 : auto ng = tk::NGvol(ndof);
354 : :
355 : : // arrays for quadrature points
356 : 909600 : std::array< std::vector< tk::real >, 3 > coordgp;
357 [ + - ]: 909600 : std::vector< tk::real > wgp( ng );
358 : :
359 [ + - ]: 454800 : coordgp[0].resize( ng );
360 [ + - ]: 454800 : coordgp[1].resize( ng );
361 [ + - ]: 454800 : coordgp[2].resize( ng );
362 : :
363 [ + - ]: 454800 : tk::GaussQuadratureTet( ng, coordgp, wgp );
364 : :
365 [ + - ]: 909600 : std::vector<tk::real> dU(nmat,0.0);
366 [ + - ]: 909600 : std::vector<tk::real> U(nmat,0.0);
367 [ + - ]: 454800 : std::vector<std::size_t> marker(nmat,0);
368 : :
369 : : // Gaussian quadrature
370 [ + + ]: 2728800 : for (std::size_t igp=0; igp<ng; ++igp)
371 : : {
372 : : // Compute the basis function
373 : 4548000 : auto B = tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp],
374 [ + - ]: 9096000 : coordgp[2][igp] );
375 : :
376 [ + - ]: 4548000 : auto state = tk::eval_state( ncomp, 0, ndof, ndofel, e, unk, B, {0, ncomp-1} );
377 : :
378 [ + + ]: 6822000 : for(std::size_t k = 0; k < nmat; k++) {
379 [ + - ][ + + ]: 4548000 : if(unk(e, volfracDofIdx(nmat, k, ndof, 0), 0) > 1e-2) {
380 : 2336410 : marker[k] = 1;
381 : 2336410 : U[k] += wgp[igp] *
382 : 2336410 : state[densityIdx(nmat, k)] * state[densityIdx(nmat, k)];
383 : :
384 [ - + ]: 2336410 : if(ndofel > 4)
385 : : {
386 [ - - ]: 0 : auto dU_p2 = unk(e, densityDofIdx(nmat, k, ndof, 4), 0) * B[4]
387 [ - - ]: 0 : + unk(e, densityDofIdx(nmat, k, ndof, 5), 0) * B[5]
388 [ - - ]: 0 : + unk(e, densityDofIdx(nmat, k, ndof, 6), 0) * B[6]
389 [ - - ]: 0 : + unk(e, densityDofIdx(nmat, k, ndof, 7), 0) * B[7]
390 [ - - ]: 0 : + unk(e, densityDofIdx(nmat, k, ndof, 8), 0) * B[8]
391 [ - - ]: 0 : + unk(e, densityDofIdx(nmat, k, ndof, 9), 0) * B[9];
392 : :
393 : 0 : dU[k] += wgp[igp] * dU_p2 * dU_p2;
394 : : }
395 : : else
396 : : {
397 [ + - ]: 2336410 : auto dU_p1 = unk(e, densityDofIdx(nmat, k, ndof, 1), 0) * B[1]
398 [ + - ]: 2336410 : + unk(e, densityDofIdx(nmat, k, ndof, 2), 0) * B[2]
399 [ + - ]: 2336410 : + unk(e, densityDofIdx(nmat, k, ndof, 3), 0) * B[3];
400 : :
401 : 2336410 : dU[k] += wgp[igp] * dU_p1 * dU_p1;
402 : : }
403 : : }
404 : : }
405 : : }
406 : :
407 : : // Return the max indicator value among all the materials
408 : 454800 : tk::real Indmax = 0.0;
409 [ + + ]: 1364400 : for(std::size_t k = 0; k < nmat; k++)
410 [ + + ]: 909600 : if(marker[k]) Indmax = std::max(dU[k]/U[k], Indmax);
411 : 909600 : return Indmax;
412 : : }
413 : :
414 : : }
415 : : // inciter::
|
