Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
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__LIB_DISC__OPERATOR__COMPOSITE_CONVERGENCE_CHECK_IMPL__
34 : #define __H__LIB_DISC__OPERATOR__COMPOSITE_CONVERGENCE_CHECK_IMPL__
35 :
36 : #include "composite_conv_check.h"
37 : #include "common/util/string_util.h"
38 : #include <boost/core/enable_if.hpp>
39 :
40 : namespace ug{
41 :
42 : ////////////////////////////////////////////////////////////////////////////////
43 : // Composite convergence check
44 : ////////////////////////////////////////////////////////////////////////////////
45 :
46 : template <class TVector, class TDomain>
47 0 : CompositeConvCheck<TVector, TDomain>::
48 : CompositeConvCheck(SmartPtr<ApproximationSpace<TDomain> > spApproxSpace)
49 : : m_spApprox(spApproxSpace),
50 0 : m_bCheckRest(true), m_restMinDefect(1-12), m_restRelReduction(1-10),
51 0 : m_currentStep(0), m_maxSteps(100),
52 0 : m_verbose(true), m_supress_unsuccessful(false),
53 0 : m_offset(0), m_symbol('%'), m_name("Iteration"), m_info(""),
54 0 : m_bTimeMeas(true), m_bAdaptive(false)
55 : {
56 0 : set_level(GridLevel::TOP);
57 0 : }
58 :
59 : template <class TVector, class TDomain>
60 0 : CompositeConvCheck<TVector, TDomain>::
61 : CompositeConvCheck(SmartPtr<ApproximationSpace<TDomain> > spApproxSpace,
62 : int maxSteps, number minDefect, number relReduction)
63 : : m_spApprox(spApproxSpace),
64 0 : m_bCheckRest(true), m_restMinDefect(minDefect), m_restRelReduction(relReduction),
65 0 : m_currentStep(0), m_maxSteps(maxSteps),
66 0 : m_verbose(true), m_supress_unsuccessful(false),
67 0 : m_offset(0), m_symbol('%'), m_name("Iteration"), m_info(""),
68 0 : m_bTimeMeas(true), m_bAdaptive(false)
69 : {
70 0 : set_level(GridLevel::TOP);
71 0 : }
72 :
73 : template <class TVector, class TDomain>
74 0 : void CompositeConvCheck<TVector, TDomain>::set_level(int level)
75 : {
76 0 : ConstSmartPtr<DoFDistribution> dd = m_spApprox->dof_distribution(GridLevel(level, GridLevel::SURFACE));
77 0 : extract_dof_indices(dd);
78 :
79 0 : update_rest_check();
80 0 : }
81 :
82 :
83 : template <class TVector, class TDomain>
84 : template <typename TBaseElem>
85 0 : void CompositeConvCheck<TVector, TDomain>::
86 : extract_dof_indices(ConstSmartPtr<DoFDistribution> dd)
87 : {
88 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
89 :
90 0 : const SurfaceView& sv = *dd->surface_view();
91 0 : const MultiGrid& mg = *dd->multi_grid();
92 :
93 : // We need to iterate over SurfaceView::ALL unknowns here since, in parallel,
94 : // it is possible for the shadowing elems of SHADOW_RIM_COPY elem to be located
95 : // on a different processor!
96 : // (cf. DoFDistribution::reinit() implementation comments)
97 :
98 : // iterate all elements (including SHADOW_RIM_COPY!)
99 0 : iter = dd->template begin<TBaseElem>(SurfaceView::ALL);
100 0 : iterEnd = dd->template end<TBaseElem>(SurfaceView::ALL);
101 :
102 : // loop over all elements
103 : std::vector<DoFIndex> vInd;
104 0 : for (; iter != iterEnd; ++iter)
105 : {
106 : TBaseElem* elem = *iter;
107 0 : if (sv.is_contained(elem, dd->grid_level(), SurfaceView::SHADOW_RIM_COPY))
108 : {
109 0 : if (mg.num_children<TBaseElem>(elem) > 0)
110 : {
111 : TBaseElem* child = mg.get_child<TBaseElem>(elem, 0);
112 0 : if (sv.is_contained(child, dd->grid_level(), SurfaceView::SURFACE_RIM))
113 0 : continue;
114 : }
115 : }
116 :
117 0 : for (size_t fi = 0; fi < dd->num_fct(); fi++)
118 : {
119 : // inner multi indices of the grid object for a function component
120 0 : dd->inner_dof_indices(elem, fi, vInd);
121 :
122 : // remember multi indices
123 0 : for (size_t dof = 0; dof < vInd.size(); dof++)
124 0 : m_vNativCmpInfo[fi].vMultiIndex.push_back(vInd[dof]);
125 : }
126 : }
127 : // note: no duplicate indices possible
128 0 : }
129 :
130 : template <class TVector, class TDomain>
131 0 : void CompositeConvCheck<TVector, TDomain>::
132 : extract_dof_indices(ConstSmartPtr<DoFDistribution> dd)
133 : {
134 : // compute indices for faster access later
135 : m_vNativCmpInfo.clear();
136 0 : m_vNativCmpInfo.resize(dd->num_fct());
137 0 : if(dd->max_dofs(VERTEX)) extract_dof_indices<Vertex>(dd);
138 0 : if(dd->max_dofs(EDGE)) extract_dof_indices<Edge>(dd);
139 0 : if(dd->max_dofs(FACE)) extract_dof_indices<Face>(dd);
140 0 : if(dd->max_dofs(VOLUME)) extract_dof_indices<Volume>(dd);
141 :
142 0 : for(size_t i = 0; i < m_vNativCmpInfo.size(); ++i)
143 0 : m_vNativCmpInfo[i].name = dd->name(i);
144 :
145 0 : m_numAllDoFs = 0;
146 0 : for(size_t i = 0; i < m_vNativCmpInfo.size(); ++i)
147 0 : m_numAllDoFs += m_vNativCmpInfo[i].vMultiIndex.size();
148 0 : }
149 :
150 :
151 : template <class TVector, class TDomain>
152 0 : number CompositeConvCheck<TVector, TDomain>::
153 : norm(const TVector& vec, const std::vector<DoFIndex>& vMultiIndex)
154 : {
155 : #ifdef UG_PARALLEL
156 :
157 : // make vector d additive unique
158 : if (!const_cast<TVector*>(&vec)->change_storage_type(PST_UNIQUE))
159 : UG_THROW("CompositeConvCheck::norm(): Cannot change ParallelStorageType to unique.");
160 : #endif
161 :
162 : double norm = 0.0;
163 : size_t sz = vMultiIndex.size();
164 0 : for (size_t dof = 0; dof < sz; ++dof)
165 : {
166 0 : const number val = DoFRef(vec, vMultiIndex[dof]);
167 0 : norm += (double) (val*val);
168 : }
169 :
170 : #ifdef UG_PARALLEL
171 : // sum squared local norms
172 : //if (!vec.layouts()->proc_comm().empty())
173 : //norm = vec.layouts()->proc_comm().allreduce(norm, PCL_RO_SUM);
174 :
175 : // In the lines above,
176 : // racing conditions occur in cases where a process has no elements,
177 : // since the defect would be 0 for them then and iteration_ended() would return true;
178 : // ergo: the empty processors would wait at the next global communication involving them
179 : // while non-empty processors might encounter a different communication event before.
180 : // This results in error messages like MPI ERROR: MPI_ERR_TRUNCATE: message truncated.
181 :
182 : // The lines below, however, result in Bi-CGSTAB going down in the first step with
183 : // "minOrthogonality failed" in settings with empty processes as they will compute
184 : // a zero (local) dot_product r*r0, but with (global) r=r0 != 0.
185 : // Therefore switching to old alternative here. MPI_ERR_TRUNCATE should be prevented
186 : // by not communicating globally, but only with the processes that really need to be
187 : // involved (i.e. NO empty processes), if possible.
188 :
189 : pcl::ProcessCommunicator commWorld;
190 : norm = commWorld.allreduce(norm, PCL_RO_SUM);
191 : #endif
192 :
193 : // return global norm
194 0 : return sqrt((number) norm);
195 : }
196 :
197 :
198 : template <class TVector, class TDomain>
199 0 : SmartPtr<IConvergenceCheck<TVector> > CompositeConvCheck<TVector, TDomain>::clone()
200 : {
201 0 : SmartPtr<CompositeConvCheck<TVector, TDomain> > newInst(
202 0 : new CompositeConvCheck<TVector, TDomain>(m_spApprox));
203 :
204 : // use std assignment (implicit member-wise is fine here)
205 0 : *newInst = *this;
206 0 : return newInst;
207 : }
208 :
209 : template <class TVector, class TDomain>
210 0 : void CompositeConvCheck<TVector, TDomain>::
211 : set_component_check(const std::vector<std::string>& vFctName,
212 : const std::vector<number>& vMinDefect,
213 : const std::vector<number>& vRelReduction)
214 : {
215 0 : if(vFctName.size() != vMinDefect.size() || vFctName.size() != vRelReduction.size())
216 0 : UG_THROW("CompositeConvCheck: Please specify one value for each cmp.");
217 :
218 0 : for(size_t cmp = 0; cmp < vFctName.size(); ++cmp)
219 0 : set_component_check(vFctName[cmp], vMinDefect[cmp], vRelReduction[cmp]);
220 :
221 0 : }
222 :
223 : template <class TVector, class TDomain>
224 0 : void CompositeConvCheck<TVector, TDomain>::
225 : set_component_check(const std::string& vFctName,
226 : const std::vector<number>& vMinDefect,
227 : const std::vector<number>& vRelReduction)
228 : {
229 0 : set_component_check(TokenizeTrimString(vFctName), vMinDefect, vRelReduction);
230 0 : }
231 :
232 :
233 : template <class TVector, class TDomain>
234 0 : void CompositeConvCheck<TVector, TDomain>::
235 : set_component_check(const std::vector<std::string>& vFctName,
236 : const number minDefect,
237 : const number relReduction)
238 : {
239 0 : for(size_t cmp = 0; cmp < vFctName.size(); ++cmp)
240 0 : set_component_check(vFctName[cmp], minDefect, relReduction);
241 0 : }
242 :
243 : template <class TVector, class TDomain>
244 0 : void CompositeConvCheck<TVector, TDomain>::
245 : set_all_component_check(const number minDefect,
246 : const number relReduction)
247 : {
248 0 : for(size_t fct = 0; fct < m_spApprox->num_fct(); ++fct){
249 : std::string name = m_spApprox->name(fct);
250 0 : set_component_check(name, minDefect, relReduction);
251 : }
252 0 : }
253 :
254 : template <class TVector, class TDomain>
255 0 : void CompositeConvCheck<TVector, TDomain>::
256 : set_component_check(const std::string& fctNames,
257 : const number minDefect,
258 : const number relReduction)
259 : {
260 0 : std::vector<std::string> vName = TokenizeTrimString(fctNames);
261 :
262 0 : for(size_t i = 0; i < vName.size(); ++i)
263 : {
264 0 : const int fct = m_spApprox->fct_id_by_name(vName[i].c_str());
265 :
266 : // search for single component
267 : size_t cmp = 0;
268 0 : for(; cmp < m_CmpInfo.size(); ++cmp){
269 :
270 : // not checking rest
271 0 : if(m_CmpInfo[cmp].isRest) continue;
272 :
273 : // only single valued checked
274 0 : if(m_CmpInfo[cmp].vFct.size() != 1) continue;
275 :
276 : // check if component found
277 0 : if(m_CmpInfo[cmp].vFct[0] == fct)
278 : break;
279 : }
280 :
281 : // if not found add new one
282 0 : if(cmp == m_CmpInfo.size())
283 0 : m_CmpInfo.push_back(CmpInfo());
284 :
285 : // set values
286 0 : m_CmpInfo[cmp].isRest = false;
287 : m_CmpInfo[cmp].vFct.clear();
288 0 : m_CmpInfo[cmp].vFct.push_back(fct);
289 0 : m_CmpInfo[cmp].name = vName[i];
290 0 : m_CmpInfo[cmp].minDefect = minDefect;
291 0 : m_CmpInfo[cmp].relReduction = relReduction;
292 : }
293 :
294 : // update rest values
295 0 : update_rest_check();
296 0 : }
297 :
298 : template <class TVector, class TDomain>
299 0 : void CompositeConvCheck<TVector, TDomain>::
300 : set_group_check(const std::vector<std::string>& vFctName,
301 : const number minDefect,
302 : const number relReduction)
303 : {
304 0 : std::vector<int> vFct(vFctName.size());
305 0 : for(size_t i = 0; i < vFctName.size(); ++i){
306 0 : std::string name = TrimString(vFctName[i]);
307 0 : vFct[i] = m_spApprox->fct_id_by_name(name.c_str());
308 : }
309 :
310 : // search for group
311 : size_t cmp = 0;
312 0 : for(; cmp < m_CmpInfo.size(); ++cmp){
313 :
314 : // not checking rest
315 0 : if(m_CmpInfo[cmp].isRest) continue;
316 :
317 : // only single valued checked
318 0 : if(m_CmpInfo[cmp].vFct.size() != vFct.size()) continue;
319 :
320 : // check if component found
321 : bool bFound = true;
322 0 : for(size_t i = 0; i < vFct.size(); ++i){
323 0 : if(m_CmpInfo[cmp].vFct[i] != vFct[i])
324 : bFound = false;
325 : }
326 :
327 0 : if(bFound) break;
328 : }
329 :
330 : // if not found add new one
331 0 : if(cmp == m_CmpInfo.size())
332 0 : m_CmpInfo.push_back(CmpInfo());
333 :
334 0 : std::string name("");
335 0 : for(size_t fct = 0; fct < vFct.size(); ++fct){
336 0 : if(!name.empty()) name.append(", ");
337 0 : name.append(m_vNativCmpInfo[vFct[fct]].name);
338 : }
339 :
340 : // set values
341 0 : m_CmpInfo[cmp].isRest = false;
342 0 : m_CmpInfo[cmp].vFct = vFct;
343 0 : m_CmpInfo[cmp].name = name;
344 0 : m_CmpInfo[cmp].minDefect = minDefect;
345 0 : m_CmpInfo[cmp].relReduction = relReduction;
346 :
347 : // update rest values
348 0 : update_rest_check();
349 0 : }
350 :
351 : template <class TVector, class TDomain>
352 0 : void CompositeConvCheck<TVector, TDomain>::
353 : set_group_check(const std::string& fctNames,
354 : const number minDefect,
355 : const number relReduction)
356 : {
357 0 : set_group_check(TokenizeTrimString(fctNames), minDefect, relReduction);
358 0 : }
359 :
360 : template <class TVector, class TDomain>
361 0 : void CompositeConvCheck<TVector, TDomain>::
362 : update_rest_check()
363 : {
364 : // remove old rest
365 0 : for(size_t cmp = 0; cmp < m_CmpInfo.size(); ++cmp)
366 0 : if(m_CmpInfo[cmp].isRest)
367 0 : m_CmpInfo.erase(m_CmpInfo.begin()+cmp);
368 :
369 : // check if rest required
370 0 : if(!m_bCheckRest) return;
371 :
372 : // detect used functions
373 0 : std::vector<bool> vUsed(m_vNativCmpInfo.size(), false);
374 0 : for(size_t cmp = 0; cmp < m_CmpInfo.size(); ++cmp)
375 : {
376 0 : for(size_t i = 0; i < m_CmpInfo[cmp].vFct.size(); ++i)
377 0 : vUsed[ m_CmpInfo[cmp].vFct[i] ] = true;
378 : }
379 :
380 : std::vector<int> vFct;
381 0 : std::string name("");
382 :
383 0 : for(size_t fct = 0; fct < vUsed.size(); ++fct){
384 0 : if(vUsed[fct]) continue;
385 :
386 0 : vFct.push_back(fct);
387 0 : if(!name.empty()) name.append(", ");
388 : name.append(m_vNativCmpInfo[fct].name);
389 : }
390 :
391 : // if no rest --> not added
392 0 : if(vFct.empty()) return;
393 :
394 : // add new rest
395 0 : m_CmpInfo.push_back(CmpInfo());
396 0 : m_CmpInfo.back().isRest = true;
397 0 : m_CmpInfo.back().vFct = vFct;
398 0 : m_CmpInfo.back().name = name;
399 0 : m_CmpInfo.back().minDefect = m_restMinDefect;
400 0 : m_CmpInfo.back().relReduction = m_restRelReduction;
401 0 : }
402 :
403 : template <class TVector, class TDomain>
404 0 : void CompositeConvCheck<TVector, TDomain>::start_defect(number initialDefect)
405 : {
406 0 : UG_THROW( "This method cannot be used to set defect values,\n"
407 : "since obviously this class is meant for an individual\n"
408 : "defect calculation of more than one function\n"
409 : "(use start(TVector& d) instead).");
410 : }
411 :
412 :
413 : template <typename TVector, typename Enable = void>
414 : struct MyVectorTraits
415 : {
416 : static const size_t block_size = 1;
417 : };
418 :
419 : // specialization for all blocked vector types
420 : template <typename TVector>
421 : struct MyVectorTraits<TVector, typename boost::enable_if_c<TVector::value_type::is_static>::type>
422 : {
423 : static const size_t block_size = TVector::value_type::static_size;
424 : };
425 :
426 :
427 : template <class TVector, class TDomain>
428 0 : void CompositeConvCheck<TVector, TDomain>::start(const TVector& vec)
429 : {
430 : // if meshing is adaptive, prepare convCheck for possibly new grid
431 0 : if (m_bAdaptive)
432 0 : set_level(GridLevel::TOP);
433 :
434 : // assert correct number of dofs
435 0 : if (vec.size() * MyVectorTraits<TVector>::block_size != m_numAllDoFs)
436 : {
437 0 : UG_THROW("Number of dofs in CompositeConvCheck does not match "
438 : "number of dofs given in vector from algorithm (" << m_numAllDoFs
439 : << ", but " << vec.size() << " given). \nMake sure that you set "
440 : "the right grid level via set_level().");
441 : }
442 :
443 : // start time measurement
444 0 : if (m_bTimeMeas) m_stopwatch.start();
445 :
446 : // update native defects
447 0 : for (size_t fct = 0; fct < m_vNativCmpInfo.size(); fct++){
448 0 : m_vNativCmpInfo[fct].initDefect = norm(vec, m_vNativCmpInfo[fct].vMultiIndex);
449 0 : m_vNativCmpInfo[fct].currDefect = m_vNativCmpInfo[fct].initDefect;
450 : }
451 :
452 : // update grouped defects
453 0 : for (size_t cmp = 0; cmp < m_CmpInfo.size(); cmp++){
454 : CmpInfo& cmpInfo = m_CmpInfo[cmp];
455 :
456 0 : cmpInfo.currDefect = 0.0;
457 0 : for(size_t i = 0; i < cmpInfo.vFct.size(); ++i)
458 0 : cmpInfo.currDefect += pow(m_vNativCmpInfo[cmpInfo.vFct[i]].currDefect, 2);
459 0 : cmpInfo.currDefect = sqrt(cmpInfo.currDefect);
460 :
461 0 : cmpInfo.initDefect = cmpInfo.currDefect;
462 : }
463 :
464 0 : m_currentStep = 0;
465 :
466 0 : if (m_verbose)
467 : {
468 : UG_LOG("\n");
469 :
470 : // number of symbols to print before name and info
471 : int num_sym = 18;
472 : int num_line_length = 80;
473 :
474 0 : int max_length = std::max(m_name.length(), m_info.length());
475 0 : int space_left = std::max(num_line_length - max_length - num_sym, 0);
476 :
477 : // print name line
478 0 : print_offset();
479 0 : UG_LOG(repeat(m_symbol, num_sym));
480 0 : int pre_space = (int)(max_length -(int)m_name.length()) / 2;
481 0 : UG_LOG(repeat(' ', pre_space));
482 : UG_LOG(" "<< m_name << " ");
483 0 : UG_LOG(repeat(' ', max_length - pre_space -m_name.length()));
484 0 : UG_LOG(repeat(m_symbol, space_left));
485 : UG_LOG("\n");
486 : // print info line
487 0 : print_offset();
488 0 : if (m_info.length() > 0)
489 : {
490 0 : UG_LOG(repeat(m_symbol, num_sym));
491 : UG_LOG(" "<< m_info << " ");
492 0 : UG_LOG(repeat(' ', max_length-m_info.length()));
493 0 : UG_LOG(repeat(m_symbol, space_left));
494 : UG_LOG("\n");
495 : }
496 : else
497 : {
498 : UG_LOG("\n");
499 : }
500 :
501 : // start iteration output
502 0 : print_offset();
503 : UG_LOG(" Iter Defect Required Rate "
504 : "Reduction Required Component(s)\n");
505 :
506 0 : for (size_t cmp = 0; cmp < m_CmpInfo.size(); cmp++)
507 : {
508 : CmpInfo& cmpInfo = m_CmpInfo[cmp];
509 :
510 0 : print_offset();
511 0 : if(cmp != 0) {UG_LOG(" " );}
512 0 : else {UG_LOG(std::right << std::setw(5) << step() << " ");}
513 :
514 0 : UG_LOG(std::scientific << cmpInfo.currDefect << " ");
515 0 : UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.minDefect << " " << std::setprecision(6) );
516 : UG_LOG(std::scientific << " ----- " << " ");
517 : UG_LOG(std::scientific << " -------- " << " ");
518 0 : UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.relReduction << " " << std::setprecision(6));
519 : UG_LOG(std::scientific << cmpInfo.name << "\n");
520 : }
521 : }
522 0 : }
523 :
524 :
525 : template <class TVector, class TDomain>
526 0 : void CompositeConvCheck<TVector, TDomain>::update_defect(number newDefect)
527 : {
528 0 : UG_THROW( "This method cannot be used to update defect values,\n"
529 : "since obviously this class is meant for an individual\n"
530 : "defect calculation of more than one function\n"
531 : "(use update(TVector& d) instead).");
532 : }
533 :
534 :
535 : template <class TVector, class TDomain>
536 0 : void CompositeConvCheck<TVector, TDomain>::update(const TVector& vec)
537 : {
538 : // assert correct number of dofs
539 0 : if (vec.size() * MyVectorTraits<TVector>::block_size != m_numAllDoFs)
540 : {
541 0 : UG_THROW("Number of dofs in CompositeConvCheck does not match"
542 : "number of dofs given in vector from algorithm (" << m_numAllDoFs <<
543 : ", but " << vec.size() << " given). \nMake sure that you set the "
544 : "right grid level via set_level().");
545 : }
546 :
547 : // update native defects
548 0 : for (size_t fct = 0; fct < m_vNativCmpInfo.size(); fct++){
549 0 : m_vNativCmpInfo[fct].lastDefect = m_vNativCmpInfo[fct].currDefect;
550 0 : m_vNativCmpInfo[fct].currDefect = norm(vec, m_vNativCmpInfo[fct].vMultiIndex);
551 : }
552 :
553 : // update grouped defects
554 0 : for (size_t cmp = 0; cmp < m_CmpInfo.size(); cmp++){
555 : CmpInfo& cmpInfo = m_CmpInfo[cmp];
556 :
557 0 : cmpInfo.lastDefect = cmpInfo.currDefect;
558 :
559 0 : cmpInfo.currDefect = 0.0;
560 0 : for(size_t i = 0; i < cmpInfo.vFct.size(); ++i)
561 0 : cmpInfo.currDefect += pow(m_vNativCmpInfo[cmpInfo.vFct[i]].currDefect, 2);
562 0 : cmpInfo.currDefect = sqrt(cmpInfo.currDefect);
563 : }
564 :
565 0 : m_currentStep++;
566 :
567 0 : if (m_verbose)
568 : {
569 0 : for (size_t cmp = 0; cmp < m_CmpInfo.size(); cmp++)
570 : {
571 : CmpInfo& cmpInfo = m_CmpInfo[cmp];
572 :
573 0 : print_offset();
574 0 : if(cmp != 0) {UG_LOG(" " );}
575 0 : else {UG_LOG(std::right << std::setw(5) << step() << " ");}
576 :
577 0 : UG_LOG(std::scientific << cmpInfo.currDefect << " ");
578 0 : UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.minDefect << " " << std::setprecision(6));
579 0 : if(cmpInfo.lastDefect != 0.0){
580 0 : UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.currDefect / cmpInfo.lastDefect << " "<< std::setprecision(6));
581 : } else {
582 : UG_LOG(std::scientific << " ----- " << " ");
583 : }
584 0 : if(cmpInfo.initDefect != 0.0){
585 0 : UG_LOG(std::scientific << cmpInfo.currDefect / cmpInfo.initDefect << " ");
586 : } else {
587 : UG_LOG(std::scientific << " -------- " << " ");
588 : }
589 0 : UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.relReduction << " " << std::setprecision(6));
590 : UG_LOG(std::scientific << cmpInfo.name << "\n");
591 : }
592 : }
593 0 : }
594 :
595 :
596 : template <class TVector, class TDomain>
597 0 : bool CompositeConvCheck<TVector, TDomain>::iteration_ended()
598 : {
599 0 : if (step() >= m_maxSteps) return true;
600 :
601 : bool ended = true;
602 0 : for (size_t cmp = 0; cmp < m_CmpInfo.size(); cmp++)
603 : {
604 : CmpInfo& cmpInfo = m_CmpInfo[cmp];
605 :
606 0 : if (!is_valid_number(cmpInfo.currDefect)) return true;
607 :
608 : bool cmpFinished = false;
609 0 : if(cmpInfo.currDefect < cmpInfo.minDefect) cmpFinished = true;
610 0 : if(cmpInfo.initDefect != 0.0)
611 0 : if((cmpInfo.currDefect/cmpInfo.initDefect) < cmpInfo.relReduction)
612 : cmpFinished = true;
613 :
614 0 : ended = ended && cmpFinished;
615 : }
616 :
617 : return ended;
618 : }
619 :
620 :
621 : template <class TVector, class TDomain>
622 0 : bool CompositeConvCheck<TVector, TDomain>::post()
623 : {
624 0 : if (m_bTimeMeas) m_stopwatch.stop();
625 :
626 : bool success = true;
627 :
628 0 : if (step() > m_maxSteps){
629 0 : print_offset();
630 0 : UG_LOG("Maximum numbers of "<< m_maxSteps <<
631 : " iterations reached without convergence.\n");
632 : success = false;
633 : }
634 :
635 0 : for (size_t cmp = 0; cmp < m_CmpInfo.size(); cmp++)
636 : {
637 : CmpInfo& cmpInfo = m_CmpInfo[cmp];
638 :
639 0 : if (!is_valid_number(cmpInfo.currDefect))
640 : {
641 : success = false;
642 0 : if (m_verbose)
643 : {
644 0 : print_offset();
645 : UG_LOG("Current defect for '" << cmpInfo.name <<
646 : "' is not a valid number.\n");
647 : }
648 : }
649 :
650 : bool cmpFinished = false;
651 0 : if (cmpInfo.currDefect < cmpInfo.minDefect)
652 : {
653 0 : if (m_verbose)
654 : {
655 0 : print_offset();
656 0 : UG_LOG("Absolute defect of " << cmpInfo.minDefect << " for '"
657 : << cmpInfo.name << "' reached after " << step() << " steps.\n");
658 : }
659 : cmpFinished = true;
660 : }
661 :
662 0 : if (cmpInfo.initDefect != 0.0)
663 : {
664 0 : if (cmpInfo.currDefect/cmpInfo.initDefect < cmpInfo.relReduction)
665 : {
666 0 : if (m_verbose)
667 : {
668 0 : print_offset();
669 0 : UG_LOG("Relative reduction of " << cmpInfo.relReduction << " for '"
670 : << cmpInfo.name << "' reached after " << step() << " steps.\n");
671 : }
672 : cmpFinished = true;
673 : }
674 : }
675 :
676 0 : if (!cmpFinished) success = false;
677 : }
678 :
679 0 : if (m_verbose)
680 : {
681 0 : print_offset();
682 0 : UG_LOG(repeat(m_symbol, 5));
683 0 : std::stringstream tmsg;
684 0 : if (m_bTimeMeas)
685 : {
686 0 : number time = m_stopwatch.ms()/1000.0;
687 : tmsg << " (t: " << std::setprecision(3) << time << "s; t/it: "
688 0 : << time / step() << "s)";
689 : }
690 0 : if (success) {UG_LOG(" Iteration converged" << tmsg.str() << " ");}
691 0 : else {UG_LOG(" Iteration not successful" << tmsg.str() << " ");}
692 0 : UG_LOG(repeat(m_symbol, 5));
693 : UG_LOG("\n\n");
694 0 : }
695 :
696 0 : return success || m_supress_unsuccessful;
697 : }
698 :
699 :
700 : template <class TVector, class TDomain>
701 0 : void CompositeConvCheck<TVector, TDomain>::print_offset()
702 : {
703 : // step 1: whitespace
704 0 : UG_LOG(repeat(' ', m_offset));
705 :
706 : // step 2: print style character
707 0 : UG_LOG(m_symbol << " ");
708 0 : }
709 :
710 : template <class TVector, class TDomain>
711 0 : void CompositeConvCheck<TVector, TDomain>::print_line(std::string line)
712 : {
713 0 : print_offset();
714 : UG_LOG(line << "\n");
715 0 : }
716 :
717 :
718 : template <class TVector, class TDomain>
719 : bool CompositeConvCheck<TVector, TDomain>::is_valid_number(number value)
720 : {
721 0 : if (value == 0.0) return true;
722 : else return (value >= std::numeric_limits<number>::min()
723 0 : && value <= std::numeric_limits<number>::max()
724 0 : && value == value && value >= 0.0);
725 : }
726 :
727 : } // end namespace ug
728 :
729 :
730 : #endif /* __H__LIB_DISC__OPERATOR__COMPOSITE_CONVERGENCE_CHECK_IMPL__ */
731 :
|
