ug4
Loading...
Searching...
No Matches
algebra_conv_check_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 * Author: Arne Nägel
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/*
34 * algebraic version of composite_conv_check.h by M. Breit
35 */
36
37#ifndef __H__LIB_ALGEBRA__OPERATOR__ALGEBRA_CONV_CHECK_IMPL__
38#define __H__LIB_ALGEBRA__OPERATOR__ALGEBRA_CONV_CHECK_IMPL__
39
40#include "algebra_conv_check.h"
43
44namespace ug{
45
47// Algebraic (composite) convergence check
49
50template <class TVector>
52:
53 m_vCmpInfo(ncmp, CmpInfo(1e-12, 1e-10)),
54 m_maxSteps(100), m_minDefect(1e-12), m_relReduction(1e-10), m_verbose(true),
55 m_currentStep(0), m_offset(0), m_symbol('%'), m_name("Iteration"), m_info(""),
56 m_bTimeMeas(true)
57{}
58
59
60template <class TVector>
62AlgebraicConvCheck(size_t ncmp, int maxSteps, number minDefect, number relReduction)
63: m_maxSteps(maxSteps), m_minDefect(minDefect), m_relReduction(relReduction), m_verbose(true),
64 m_currentStep(0), m_offset(0), m_symbol('%'), m_name("Iteration"), m_info(""),
65 m_bTimeMeas(true), m_vCmpInfo(ncmp, CmpInfo(minDefect, relReduction))
66{
67 //const CmpInfo(minDefect, relReduction);
68
69}
70
71template <class TVector>
73AlgebraicConvCheck(size_t ncmp, int maxSteps, number minDefect, number relReduction, bool verbose)
74: m_maxSteps(maxSteps), m_minDefect(minDefect), m_relReduction(relReduction), m_verbose(verbose),
75 m_currentStep(0), m_offset(0), m_symbol('%'), m_name("Iteration"), m_info(""),
76 m_bTimeMeas(true)
77{
78 //const CmpInfo(minDefect, relReduction);
79 //m_vCmpInfo(ncmp, a);
80}
81
82
83
84
85
86template <class TVector>
88{
89 SmartPtr<AlgebraicConvCheck<TVector> > newInst(new AlgebraicConvCheck<TVector>(this->m_vCmpInfo.size()));
90
91 // use std assignment (implicit member-wise is fine here)
92 *newInst = *this;
93 return newInst;
94}
95
96
97
98template <class TVector>
100{
101 UG_THROW( "This method cannot be used to set defect values,\n"
102 "since obviously this class is meant for an individual\n"
103 "defect calculation of more than one function\n"
104 "(use start(TVector& d) instead).");
105}
106
107
108template <class TVector>
110{
111 // start time measurement
112 if (m_bTimeMeas) m_stopwatch.start();
113
114 // update defects
115 for (size_t fct = 0; fct < m_vCmpInfo.size(); fct++){
116 m_vCmpInfo[fct].currDefect = m_vCmpInfo[fct].initDefect = norm(vec, fct);
117 }
118
119 m_currentStep = 0;
120
121 if (m_verbose)
122 {
123 UG_LOG("\n");
124
125 // number of symbols to print before name and info
126 int num_sym = 18;
127 int num_line_length = 80;
128
129 int max_length = std::max(m_name.length(), m_info.length());
130 int space_left = std::max(num_line_length - max_length - num_sym, 0);
131
132 // print name line
133 print_offset();
134 UG_LOG(repeat(m_symbol, num_sym));
135 int pre_space = (int)(max_length -(int)m_name.length()) / 2;
136 UG_LOG(repeat(' ', pre_space));
137 UG_LOG(" "<< m_name << " ");
138 UG_LOG(repeat(' ', max_length - pre_space -m_name.length()));
139 UG_LOG(repeat(m_symbol, space_left));
140 UG_LOG("\n");
141 // print info line
142 print_offset();
143 if (m_info.length() > 0)
144 {
145 UG_LOG(repeat(m_symbol, num_sym));
146 UG_LOG(" "<< m_info << " ");
147 UG_LOG(repeat(' ', max_length-m_info.length()));
148 UG_LOG(repeat(m_symbol, space_left))
149 UG_LOG("\n");
150 }
151 else
152 {
153 UG_LOG("\n");
154 }
155
156 // start iteration output
157 print_offset();
158 UG_LOG(" Iter Defect Rate \n");
159
160 for (size_t cmp = 0; cmp < m_vCmpInfo.size(); cmp++)
161 {
162 CmpInfo& cmpInfo = m_vCmpInfo[cmp];
163
164 print_offset();
165 if(cmp != 0) {UG_LOG(" " )}
166 else {UG_LOG(std::setw(4) << step() << ": ");}
167
168
169 UG_LOG(std::scientific << cmpInfo.currDefect << " ")
170 UG_LOG(std::scientific << " ----- " << " c=" << cmp << std::endl)
171
172 //UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.minDefect << " " << std::setprecision(6) )
173
174 //UG_LOG(std::scientific << " -------- " << " ")
175 //UG_LOG(std::scientific << std::setprecision(3) << cmpInfo.relReduction << " " << std::setprecision(6))
176 //UG_LOG(std::scientific << cmpInfo.name << "\n");
177 }
178 }
179}
180
181
182template <class TVector>
184{
185 UG_THROW( "This method cannot be used to update defect values,\n"
186 "since obviously this class is meant for an individual\n"
187 "defect calculation of more than one function\n"
188 "(use update(TVector& d) instead).");
189}
190
191
192template <class TVector>
194{
195 // update defects
196 for (size_t fct = 0; fct < m_vCmpInfo.size(); fct++){
197 m_vCmpInfo[fct].lastDefect = m_vCmpInfo[fct].currDefect;
198 m_vCmpInfo[fct].currDefect = norm(vec, fct);
199 }
200
201 //vec.print(); // debug only!
202
203 m_currentStep++;
204
205 if (m_verbose)
206 {
207 for (size_t cmp = 0; cmp < m_vCmpInfo.size(); cmp++)
208 {
209 CmpInfo& cmpInfo = m_vCmpInfo[cmp];
210
211 print_offset();
212 if(cmp != 0) {UG_LOG(" " )}
213 else {UG_LOG(std::setw(4) << step() << ": ");}
214
215 UG_LOG(std::scientific << cmpInfo.currDefect << " ")
216 if(cmpInfo.lastDefect != 0.0){
217 UG_LOG(std::scientific << std::setprecision(6) << cmpInfo.currDefect / cmpInfo.lastDefect << " c="<< cmp<< std::setprecision(6) << std::endl)
218 } else {
219 UG_LOG(std::scientific << " ----- " << " c="<< cmp << std::endl)
220 }
221 }
222 }
223}
224
225
226template <class TVector>
228{
229 if (step() >= m_maxSteps) return true;
230
231 bool ended = true;
232 for (size_t cmp = 0; cmp < m_vCmpInfo.size(); cmp++)
233 {
234 CmpInfo& cmpInfo = m_vCmpInfo[cmp];
235
236 if (!is_valid_number(cmpInfo.currDefect)) return true;
237
238 bool cmpFinished = false;
239 if(cmpInfo.currDefect < cmpInfo.minDefect) cmpFinished = true;
240 if(cmpInfo.initDefect != 0.0)
241 if((cmpInfo.currDefect/cmpInfo.initDefect) < cmpInfo.relReduction)
242 cmpFinished = true;
243
244 ended = ended && cmpFinished;
245 }
246
247 return ended;
248}
249
250
251template <class TVector>
253{
254 if (m_bTimeMeas) m_stopwatch.stop();
255
256 bool success = true;
257
258 if (step() >= m_maxSteps){
259 print_offset();
260 UG_LOG("Maximum numbers of "<< m_maxSteps <<
261 " iterations reached without convergence.\n");
262 success = false;
263 }
264
265 for (size_t cmp = 0; cmp < m_vCmpInfo.size(); cmp++)
266 {
267 CmpInfo& cmpInfo = m_vCmpInfo[cmp];
268
269 if (!is_valid_number(cmpInfo.currDefect)){
270 success = false;
271 if(m_verbose){
272 print_offset();
273 UG_LOG("Current defect for '" << cmpInfo.name <<
274 "' is not a valid number.\n");
275 }
276 }
277
278 bool cmpFinished = false;
279 if (cmpInfo.currDefect < cmpInfo.minDefect)
280 {
281 print_offset();
282 UG_LOG("Absolute defect of " << cmpInfo.minDefect << " for '"
283 << cmpInfo.name << "' reached after " << step() << " steps.\n");
284 cmpFinished = true;
285 }
286
287 if(cmpInfo.initDefect != 0.0){
288 if (cmpInfo.currDefect/cmpInfo.initDefect < cmpInfo.relReduction)
289 {
290 print_offset();
291 UG_LOG("Relative reduction of " << cmpInfo.relReduction << " for '"
292 << cmpInfo.name << "' reached after " << step() << " steps.\n");
293 cmpFinished = true;
294 }
295 }
296
297 if(!cmpFinished) success = false;
298 }
299
300 if(m_verbose){
301 print_offset();
302 UG_LOG(repeat(m_symbol, 5));
303 std::stringstream tmsg;
304 if (m_bTimeMeas)
305 {
306 number time = m_stopwatch.ms()/1000.0;
307 tmsg << " (t: " << std::setprecision(3) << time << "s; t/it: "
308 << time / step() << "s)";
309 }
310 if (success) {UG_LOG(" Iteration converged" << tmsg.str() << " ");}
311 else {UG_LOG(" Iteration not successful" << tmsg.str() << " ");}
312 UG_LOG(repeat(m_symbol, 5));
313 UG_LOG("\n\n");
314 }
315
316 return success;
317}
318
319
320template <class TVector>
322{
323 // step 1: whitespace
324 UG_LOG(repeat(' ', m_offset));
325
326 // step 2: print style character
327 UG_LOG(m_symbol << " ");
328}
329
330template <class TVector>
332{
333 print_offset();
334 UG_LOG(line << "\n");
335}
336
337
338template <class TVector>
340{
341 if (value == 0.0) return true;
342 else return (value >= std::numeric_limits<number>::min()
343 && value <= std::numeric_limits<number>::max()
344 && value == value && value >= 0.0);
345}
346
347
348
349
350
351template <class TVector>
353norm(const TVector& vec, size_t cmp)
354{
355#ifdef UG_PARALLEL
356
357 // make vector d additive unique
358 if (!const_cast<TVector*>(&vec)->change_storage_type(PST_UNIQUE))
359 UG_THROW("AlgebraicConvCheck::norm(): Cannot change ParallelStorageType to unique.");
360#endif
361
362
364 double norm = VecNormSquared(dummy);
365/* for (size_t dof = 0; dof < vMultiIndex.size(); ++dof)
366 {
367 const number val = BlockRef(vec, vMultiIndex[dof]);
368 norm += (double) (val*val);
369 }
370*/
371#ifdef UG_PARALLEL
372 // sum squared local norms
373 norm = vec.layouts()->proc_comm().allreduce(norm, PCL_RO_SUM);
374#endif
375
376 // return global norm
377 return sqrt((number) norm);
378}
379
380} // end namespace ug
381
382
383#endif /* __H__LIB_DISC__OPERATOR__COMPOSITE_CONVERGENCE_CHECK_IMPL__ */
384
location verbose
Definition checkpoint_util.lua:128
Definition smart_pointer.h:108
Definition algebra_conv_check.h:70
virtual SmartPtr< IConvergenceCheck< TVector > > clone()
clones this instance
Definition algebra_conv_check_impl.h:87
void update(const TVector &d)
computes the defect and sets it a the next defect value
Definition algebra_conv_check_impl.h:193
bool is_valid_number(number value)
Definition algebra_conv_check_impl.h:339
number norm(const TVector &vec, size_t cmp)
calculates the 2-norm of the entries of the vector vec specified by index
Definition algebra_conv_check_impl.h:353
void update_defect(number newDefect)
sets the update for the current defect
Definition algebra_conv_check_impl.h:183
bool iteration_ended()
Definition algebra_conv_check_impl.h:227
void start_defect(number initialDefect)
defect control
Definition algebra_conv_check_impl.h:99
AlgebraicConvCheck(size_t ncmp)
Definition algebra_conv_check_impl.h:51
void print_line(std::string line)
prints a line using prefixes
Definition algebra_conv_check_impl.h:331
void start(const TVector &d)
computes the start defect and set it
Definition algebra_conv_check_impl.h:109
void print_offset()
Definition algebra_conv_check_impl.h:321
bool post()
Definition algebra_conv_check_impl.h:252
Definition scalar_subvector_adapter.hh:53
@ PST_UNIQUE
Definition parallel_storage_type.h:70
#define PCL_RO_SUM
Definition pcl_methods.h:63
string repeat(char c, int nr)
Builds a string with specified repetitions of given character.
Definition string_util.cpp:346
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
double VecNormSquared(const double &a)
returns norm_2^2(a)
Definition operations_vec.h:100
Definition algebra_conv_check.h:156
number minDefect
Minimal required Defect of component.
Definition algebra_conv_check.h:164
number initDefect
Initial Defect of component.
Definition algebra_conv_check.h:160
number lastDefect
Last Defect if component.
Definition algebra_conv_check.h:162
number currDefect
Current Defect of component.
Definition algebra_conv_check.h:161
number relReduction
Relative reduction required for component.
Definition algebra_conv_check.h:165
std::string name
Name of components.
Definition algebra_conv_check.h:158