Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG_DISC__ERROR_INDICATOR_UTIL__
34 : #define __H__UG_DISC__ERROR_INDICATOR_UTIL__
35 :
36 : #include "lib_grid/multi_grid.h"
37 : #include "lib_grid/refinement/refiner_interface.h"
38 : #include "lib_disc/dof_manager/dof_distribution.h"
39 :
40 : namespace ug{
41 :
42 :
43 : /// helper function that computes min/max and total of error indicators
44 : /**
45 : * @param[out] min minimal eta_i
46 : * @param[out] max maximal eta_i
47 : * @param[out] sum sum of eta_i
48 : * @param[out] errSq sum of eta_i^2
49 : * @param[out] numElem number of elements
50 : * @return average eta
51 : */
52 : template<typename TElem>
53 0 : number ComputeAvg
54 : ( MultiGrid::AttachmentAccessor<TElem, ug::Attachment<number> >& aaError2,
55 : ConstSmartPtr<DoFDistribution> dd,
56 : number& min, number& max, number& sum, number& errSq, size_t& numElem,
57 : number& minLocal, number& maxLocal, number& sumLocal, size_t& numElemLocal
58 : )
59 : {
60 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
61 :
62 : // reset maximum of error
63 0 : max = 0.0, min = std::numeric_limits<number>::max();
64 0 : sum = 0.0;
65 0 : errSq = 0.0;
66 0 : numElem = 0;
67 :
68 : // get element iterator
69 0 : const_iterator iter = dd->template begin<TElem>();
70 0 : const_iterator iterEnd = dd->template end<TElem>();
71 :
72 : // loop all elements to find the maximum of the error
73 0 : for (; iter != iterEnd; ++iter)
74 : {
75 : // get element
76 : TElem* elem = *iter;
77 :
78 0 : const number elemEta2 = aaError2[elem];
79 :
80 : // if no error value exists: ignore (might be newly added by refinement);
81 : // newly added elements are supposed to have a negative error estimator
82 0 : if (aaError2[elem] < 0) UG_THROW("error value invalid!");//continue;
83 :
84 :
85 0 : const number elemEta = sqrt(aaError2[elem]);
86 :
87 : // search for maximum and minimum
88 0 : if (elemEta > max) max = elemEta;
89 0 : if (elemEta < min) min = elemEta;
90 :
91 : // sum up total error
92 0 : sum += elemEta;
93 0 : errSq += elemEta2;
94 0 : ++numElem;
95 : }
96 :
97 : // set local variables
98 0 : maxLocal = max;
99 0 : minLocal = min;
100 0 : sumLocal = sum;
101 : #ifdef UG_PARALLEL
102 : number errLocal = errSq;
103 : #endif
104 0 : numElemLocal = numElem;
105 :
106 : #ifdef UG_PARALLEL
107 : if (pcl::NumProcs() > 1)
108 : {
109 : pcl::ProcessCommunicator com;
110 : max = com.allreduce(maxLocal, PCL_RO_MAX);
111 : min = com.allreduce(minLocal, PCL_RO_MIN);
112 : sum = com.allreduce(sumLocal, PCL_RO_SUM);
113 : errSq = com.allreduce(errLocal, PCL_RO_SUM);
114 : numElem = com.allreduce(numElemLocal, PCL_RO_SUM);
115 : }
116 : #endif
117 0 : UG_LOG(" +++++ Error indicator on " << numElem << " elements +++++\n");
118 0 : UG_LOG(" +++ Element errors: maxEta=" << max << ", minEta="
119 : << min << ", sumEta=" << sum << ", sumEtaSq=" << errSq << ".\n");
120 :
121 0 : return (sum/numElem);
122 : }
123 :
124 :
125 : /// helper function that computes min/max and total of error indicators
126 : /**
127 : * @param[out] min minimal eta_i^2
128 : * @param[out] max maximal eta_i^2
129 : * @param[out] totalErr sum of eta_i^2
130 : * @param[out] numElem number of elements
131 : *
132 : */
133 : template<typename TElem>
134 0 : void ComputeMinMax
135 : ( MultiGrid::AttachmentAccessor<TElem, ug::Attachment<number> >& aaError,
136 : ConstSmartPtr<DoFDistribution> dd,
137 : number& min, number& max, number& totalErr, size_t& numElem,
138 : number& minLocal, number& maxLocal, number& totalErrLocal, size_t& numElemLocal
139 : )
140 : {
141 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
142 :
143 : // reset maximum of error
144 0 : max = 0.0, min = std::numeric_limits<number>::max();
145 0 : totalErr = 0.0;
146 0 : numElem = 0;
147 :
148 : // get element iterator
149 0 : const_iterator iter = dd->template begin<TElem>();
150 0 : const_iterator iterEnd = dd->template end<TElem>();
151 :
152 : // loop all elements to find the maximum of the error
153 0 : for (; iter != iterEnd; ++iter)
154 : {
155 : // get element
156 : TElem* elem = *iter;
157 :
158 : // if no error value exists: ignore (might be newly added by refinement);
159 : // newly added elements are supposed to have a negative error estimator
160 0 : if (aaError[elem] < 0) continue;
161 :
162 : // search for maximum and minimum
163 0 : if (aaError[elem] > max) max = aaError[elem];
164 0 : if (aaError[elem] < min) min = aaError[elem];
165 :
166 : // sum up total error
167 0 : totalErr += aaError[elem];
168 0 : ++numElem;
169 : }
170 :
171 : // set local variables
172 0 : maxLocal = max;
173 0 : minLocal = min;
174 0 : totalErrLocal = totalErr;
175 0 : numElemLocal = numElem;
176 :
177 : #ifdef UG_PARALLEL
178 : if (pcl::NumProcs() > 1)
179 : {
180 : pcl::ProcessCommunicator com;
181 : max = com.allreduce(maxLocal, PCL_RO_MAX);
182 : min = com.allreduce(minLocal, PCL_RO_MIN);
183 : totalErr = com.allreduce(totalErrLocal, PCL_RO_SUM);
184 : numElem = com.allreduce(numElemLocal, PCL_RO_SUM);
185 : }
186 : #endif
187 0 : UG_LOG(" +++++ Error indicator on " << numElem << " elements +++++\n");
188 0 : UG_LOG(" +++ Element errors: maxEtaSq=" << max << ", minEtaSq="
189 : << min << ", sumEtaSq=" << totalErr << ".\n");
190 0 : }
191 :
192 :
193 : /// helper function that computes min/max and total of error indicators
194 : template<typename TElem>
195 0 : void ComputeMinMaxTotal
196 : ( MultiGrid::AttachmentAccessor<TElem, ug::Attachment<number> >& aaError2,
197 : ConstSmartPtr<DoFDistribution> dd,
198 : number& min, number& max, number& totalErr, size_t& numElem
199 :
200 : )
201 : {
202 : number minLocal, maxLocal, totalErrLocal;
203 : size_t numElemLocal;
204 0 : ComputeMinMax(aaError2, dd, min, max, totalErr, numElem, minLocal, maxLocal, totalErrLocal, numElemLocal);
205 0 : }
206 : /// marks elements according to an attached error value field
207 : /**
208 : * This function marks elements for refinement and coarsening. The passed error attachment
209 : * is used as a weight for the amount of the error an each element. All elements
210 : * that have an indicated error with s* max <= err <= max are marked for refinement.
211 : * Here, max is the maximum error measured, s is a scaling quantity chosen by
212 : * the user. In addition, all elements with an error smaller than TOL
213 : * (user defined) are not refined.
214 : *
215 : * \param[in, out] refiner Refiner, elements marked on exit
216 : * \param[in] dd dof distribution
217 : * \param[in] TOL Minimum error, such that an element is marked
218 : * \param[in] scale scaling factor indicating lower bound for marking
219 : * \param[in] aaError Error value attachment to elements (\f$ \eta_i^2 \f$)
220 : */
221 : template<typename TElem>
222 0 : void MarkElements(MultiGrid::AttachmentAccessor<TElem, ug::Attachment<number> >& aaError,
223 : IRefiner& refiner,
224 : ConstSmartPtr<DoFDistribution> dd,
225 : number TOL,
226 : number refineFrac, number coarseFrac,
227 : int maxLevel)
228 : {
229 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
230 :
231 : // compute minimal/maximal/ total error and number of elements
232 : number min, max, totalErr;
233 : size_t numElem;
234 0 : ComputeMinMaxTotal(aaError, dd, min, max, totalErr, numElem);
235 :
236 : // check if total error is smaller than tolerance. If that is the case we're done
237 0 : if(totalErr < TOL)
238 : {
239 : UG_LOG(" +++ Total error "<<totalErr<<" smaller than TOL ("<<TOL<<"). done.");
240 0 : return;
241 : }
242 :
243 : // Compute minimum
244 0 : number minErrToRefine = max * refineFrac;
245 : UG_LOG(" +++ Refining elements if error greater " << refineFrac << "*" << max <<
246 : " = " << minErrToRefine << ".\n");
247 0 : number maxErrToCoarse = min * (1 + coarseFrac);
248 0 : if(maxErrToCoarse < TOL/numElem) maxErrToCoarse = TOL/numElem;
249 : UG_LOG(" +++ Coarsening elements if error smaller " << maxErrToCoarse << ".\n");
250 :
251 : // reset counter
252 : int numMarkedRefine = 0, numMarkedCoarse = 0;
253 :
254 0 : const_iterator iter = dd->template begin<TElem>();
255 0 : const_iterator iterEnd = dd->template end<TElem>();
256 :
257 : // loop elements for marking
258 0 : for(; iter != iterEnd; ++iter)
259 : {
260 : // get element
261 : TElem* elem = *iter;
262 :
263 : // marks for refinement
264 0 : if(aaError[elem] >= minErrToRefine)
265 0 : if(dd->multi_grid()->get_level(elem) <= maxLevel)
266 : {
267 0 : refiner.mark(elem, RM_REFINE);
268 0 : numMarkedRefine++;
269 : }
270 :
271 : // marks for coarsening
272 0 : if(aaError[elem] <= maxErrToCoarse)
273 : {
274 0 : refiner.mark(elem, RM_COARSEN);
275 0 : numMarkedCoarse++;
276 : }
277 : }
278 :
279 : #ifdef UG_PARALLEL
280 : if(pcl::NumProcs() > 1)
281 : {
282 : UG_LOG(" +++ Marked for refinement on Proc "<<pcl::ProcRank()<<": " << numMarkedRefine << " Elements.\n");
283 : UG_LOG(" +++ Marked for coarsening on Proc "<<pcl::ProcRank()<<": " << numMarkedCoarse << " Elements.\n");
284 : pcl::ProcessCommunicator com;
285 : int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
286 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
287 : numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
288 : }
289 : #endif
290 :
291 0 : UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " Elements.\n");
292 0 : UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
293 : }
294 :
295 : /// marks elements according for refinement to an attached error value field
296 : /**
297 : * This function marks elements for refinement. The passed error attachment
298 : * is used as an indicator for the the error on each element.
299 : * Elements are refined if the error sum taken over all elements is greater
300 : * than the tolerance value tol. In that case, elements with an indicated
301 : * error of err >= tol / #elems are marked for refinement if and only if their
302 : * multigrid level is below the tolerated maximum of maxLevel.
303 : *
304 : * \param[in] aaError error value attachment to elements (\f$ \eta_i^2 \f$)
305 : * \param[in, out] refiner refiner, elements marked on exit
306 : * \param[in] dd dof distribution
307 : * \param[in] tol tolerated global error (no refinement if error below)
308 : * \param[in] maxLevel maximal refinement level in multigrid
309 : */
310 :
311 : template<typename TElem>
312 0 : void MarkElementsForRefinement
313 : ( MultiGrid::AttachmentAccessor<TElem, ug::Attachment<number> >& aaError,
314 : IRefiner& refiner,
315 : ConstSmartPtr<DoFDistribution> dd,
316 : number tol,
317 : int maxLevel
318 : )
319 : {
320 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
321 :
322 : // compute minimal/maximal/ total error and number of elements
323 : number min, max, totalErr;
324 : size_t numElem;
325 0 : ComputeMinMaxTotal(aaError, dd, min, max, totalErr, numElem);
326 :
327 : // check if total error is smaller than tolerance; if that is the case we're done
328 0 : if (totalErr < tol)
329 : {
330 : UG_LOG(" +++ Total error "<<totalErr<<" smaller than TOL (" << tol << "). "
331 : "No refinement necessary." << std::endl);
332 0 : return;
333 : }
334 :
335 : // compute minimum
336 : //number minErrToRefine = max * refineFrac;
337 0 : number minErrToRefine = tol / numElem;
338 : UG_LOG(" +++ Refining elements if error greater " << tol << "/" << numElem <<
339 : " = " << minErrToRefine << ".\n");
340 :
341 : // reset counter
342 : size_t numMarkedRefine = 0;
343 :
344 0 : const_iterator iter = dd->template begin<TElem>();
345 0 : const_iterator iterEnd = dd->template end<TElem>();
346 :
347 : // loop elements for marking
348 0 : for (; iter != iterEnd; ++iter)
349 : {
350 : // get element
351 : TElem* elem = *iter;
352 :
353 : // if no error value exists: ignore (might be newly added by refinement);
354 : // newly added elements are supposed to have a negative error estimator
355 0 : if (aaError[elem] < 0) continue;
356 :
357 : // marks for refinement
358 0 : if (aaError[elem] >= minErrToRefine)
359 0 : if (dd->multi_grid()->get_level(elem) < maxLevel)
360 : {
361 0 : refiner.mark(elem, RM_REFINE);
362 0 : numMarkedRefine++;
363 : }
364 : }
365 :
366 : #ifdef UG_PARALLEL
367 : if (pcl::NumProcs() > 1)
368 : {
369 : pcl::ProcessCommunicator com;
370 : size_t numMarkedRefineLocal = numMarkedRefine;
371 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
372 : }
373 : #endif
374 :
375 : UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " elements.\n");
376 : }
377 :
378 : /// marks elements for coarsening according to an attached error value field
379 : /**
380 : * This function marks elements for coarsening. The passed error attachment
381 : * is used as an indicator for the the error on each element.
382 : * Elements one level below surface are marked if the average error of its
383 : * children is lower than tol / #elems / 4 / safety. The safety factor is
384 : * supposed to ensure that elements are not refined and coarsened back and
385 : * forth in a dynamic adaptive simulation.
386 : *
387 : * \param[in] aaError error value attachment to elements (\f$ \eta_i^2 \f$)
388 : * \param[in, out] refiner refiner, elements marked on exit
389 : * \param[in] dd dof distribution
390 : * \param[in] tol tolerated global error
391 : * \param[in] safety safety factor
392 : */
393 : template<typename TElem>
394 0 : void MarkElementsForCoarsening
395 : ( MultiGrid::AttachmentAccessor<TElem,
396 : ug::Attachment<number> >& aaError,
397 : IRefiner& refiner,
398 : ConstSmartPtr<DoFDistribution> dd,
399 : number TOL,
400 : number safety,
401 : int minLevel = 0
402 : )
403 : {
404 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
405 :
406 : // compute minimal/maximal/ total error and number of elements
407 : number min, max, totalErr;
408 : size_t numElem;
409 0 : ComputeMinMaxTotal(aaError, dd, min, max, totalErr, numElem);
410 :
411 : // compute maximum
412 : //number maxErrToCoarse = min * (1+coarseFrac);
413 : //if (maxErrToCoarse < TOL/numElem) maxErrToCoarse = TOL/numElem;
414 0 : number maxErrToCoarse = TOL / numElem / 4.0 / safety; // = tol_per_elem / (2^2*safety_factor)
415 : UG_LOG(" +++ Coarsening elements if avg child error smaller than "<< maxErrToCoarse << ".\n");
416 :
417 : // reset counter
418 : size_t numMarkedCoarse = 0;
419 :
420 0 : const_iterator iter = dd->template begin<TElem>();
421 0 : const_iterator iterEnd = dd->template end<TElem>();
422 :
423 : // loop elements for marking
424 0 : for (; iter != iterEnd; ++iter)
425 : {
426 : // get element
427 : TElem* elem = *iter;
428 :
429 : // ignore if already marked for coarsening
430 0 : if (refiner.get_mark(elem) & RM_COARSEN) continue;
431 :
432 : // ignore if level too low
433 0 : if (dd->multi_grid()->get_level(elem) <= minLevel)
434 0 : continue;
435 :
436 : // get parent
437 0 : TElem* parent = dynamic_cast<TElem*>(dd->multi_grid()->get_parent(elem));
438 0 : if (!parent) continue; // either dynamic casting failed or parent does not exist
439 :
440 : // sum up error values over all children of parent
441 : number sum = 0.0;
442 : size_t cnt = 0;
443 0 : size_t nCh = dd->multi_grid()->num_children<TElem>(parent);
444 0 : for (size_t ch = 0; ch < nCh; ch++)
445 : {
446 0 : TElem* child = dd->multi_grid()->get_child<TElem>(parent, ch);
447 :
448 : // if no error value exists: ignore (might be newly added by refinement);
449 : // newly added elements are supposed to have a negative error estimator
450 0 : if (aaError[child] < 0) continue;
451 :
452 0 : sum += aaError[child];
453 0 : cnt++;
454 : }
455 :
456 : // marks for coarsening
457 0 : if (sum/cnt <= maxErrToCoarse)
458 : {
459 0 : for (size_t ch = 0; ch < nCh; ch++)
460 : {
461 0 : TElem* child = dd->multi_grid()->get_child<TElem>(parent, ch);
462 0 : if (refiner.get_mark(child) & RM_COARSEN) continue;
463 0 : refiner.mark(child, RM_COARSEN);
464 0 : numMarkedCoarse++;
465 : }
466 : }
467 : }
468 :
469 : #ifdef UG_PARALLEL
470 : if (pcl::NumProcs() > 1)
471 : {
472 : pcl::ProcessCommunicator com;
473 : size_t numMarkedCoarseLocal = numMarkedCoarse;
474 : numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
475 : }
476 : #endif
477 :
478 : UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
479 0 : }
480 :
481 : /// marks elements according to an attached error value field
482 : /**
483 : * This function marks elements for refinement. The passed error attachment
484 : * is used as a weight for the amount of the error an each element. All elements
485 : * that have an indicated error > refineTol are marked for refinement and
486 : * elements with an error < coarsenTol are marked for coarsening
487 : *
488 : * \param[in, out] refiner Refiner, elements marked on exit
489 : * \param[in] dd dof distribution
490 : * \param[in] refTol all elements with error > refTol are marked for refinement.
491 : * If refTol is negative, no element will be marked for refinement.
492 : * \param[in] coarsenTol all elements with error < coarsenTol are marked for coarsening.
493 : * If coarsenTol is negative, no element will be marked for coarsening.
494 : * \param[in] aaError Error value attachment to elements
495 : */
496 : template<typename TElem>
497 0 : void MarkElementsAbsolute(MultiGrid::AttachmentAccessor<TElem, ug::Attachment<number> >& aaError,
498 : IRefiner& refiner,
499 : ConstSmartPtr<DoFDistribution> dd,
500 : number refTol,
501 : number coarsenTol,
502 : int minLevel,
503 : int maxLevel,
504 : bool refTopLvlOnly = false)
505 : {
506 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
507 :
508 : int numMarkedRefine = 0, numMarkedCoarse = 0;
509 0 : const_iterator iter = dd->template begin<TElem>();
510 0 : const_iterator iterEnd = dd->template end<TElem>();
511 0 : const MultiGrid* mg = dd->multi_grid().get();
512 : int topLvl = 0;
513 0 : if(mg)
514 0 : topLvl = (int)mg->top_level();
515 : else
516 : refTopLvlOnly = false;
517 :
518 : // loop elements for marking
519 0 : for(; iter != iterEnd; ++iter)
520 : {
521 : TElem* elem = *iter;
522 :
523 : // marks for refinement
524 0 : if((refTol >= 0)
525 0 : && (aaError[elem] > refTol)
526 0 : && (dd->multi_grid()->get_level(elem) < maxLevel)
527 0 : && ((!refTopLvlOnly) || (mg->get_level(elem) == topLvl)))
528 : {
529 0 : refiner.mark(elem, RM_REFINE);
530 0 : numMarkedRefine++;
531 : }
532 :
533 : // marks for coarsening
534 0 : if((coarsenTol >= 0)
535 0 : && (aaError[elem] < coarsenTol)
536 0 : && (dd->multi_grid()->get_level(elem) > minLevel))
537 : {
538 0 : refiner.mark(elem, RM_COARSEN);
539 0 : numMarkedCoarse++;
540 : }
541 : }
542 :
543 : #ifdef UG_PARALLEL
544 : if(pcl::NumProcs() > 1)
545 : {
546 : UG_LOG(" +++ Marked for refinement on Proc "<<pcl::ProcRank()<<": " << numMarkedRefine << " Elements.\n");
547 : UG_LOG(" +++ Marked for coarsening on Proc "<<pcl::ProcRank()<<": " << numMarkedCoarse << " Elements.\n");
548 : pcl::ProcessCommunicator com;
549 : int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
550 : numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
551 : numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
552 : }
553 : #endif
554 :
555 0 : UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " Elements.\n");
556 0 : UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
557 0 : }
558 :
559 : }// end of namespace
560 :
561 : #endif
|
