oomph_utilities.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1295 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-07-27 08:03:47 +0100 (Thu, 27 Jul 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifdef OOMPH_HAS_MPI
31 #include "mpi.h"
32 #endif
33 
34 #include <algorithm>
35 #include <limits.h>
36 #include <cstring>
37 
38 #ifdef OOMPH_HAS_UNISTDH
39 #include <unistd.h> // for getpid()
40 #endif
41 
42 #include "oomph_utilities.h"
43 #include "Vector.h"
44 #include "matrices.h"
45 
46 namespace oomph
47 {
48 
49 //=====================================================================
50 /// Helper namespace containing function that computes second invariant
51 /// of tensor
52 //=====================================================================
53  namespace SecondInvariantHelper
54  {
55 
56  /// Compute second invariant of tensor
57  double second_invariant(const DenseMatrix<double>& tensor)
58  {
59  // get size of tensor
60  unsigned n=tensor.nrow();
61  double trace_of_tensor=0.0;
62  double trace_of_tensor_squared=0.0;
63  for(unsigned i=0;i<n;i++)
64  {
65  // Calculate the trace of the tensor: T_ii
66  trace_of_tensor += tensor(i,i);
67  // Calculate the trace of the tensor squared: T_ij T_ji
68  for(unsigned j=0;j<n;j++)
69  {
70  trace_of_tensor_squared += tensor(i,j)*tensor(j,i);
71  }
72  }
73 
74  // return the second invariant: 1.0/2.0*[(T_ii)^2 - T_ij T_ji]
75  return 0.5*(trace_of_tensor*trace_of_tensor-trace_of_tensor_squared);
76 
77  }
78 
79  }
80 
81 
82 //==============================================
83 /// Namespace for error messages for broken
84 /// copy constructors and assignment operators
85 //==============================================
86 namespace BrokenCopy
87 {
88 
89  /// Issue error message and terminate execution
90  void broken_assign(const std::string& class_name)
91  {
92  //Write the error message into a string
93  std::string error_message = "Assignment operator for class\n\n";
94  error_message += class_name;
95  error_message += "\n\n";
96  error_message += "is deliberately broken to avoid the accidental \n";
97  error_message += "use of the inappropriate C++ default.\n";
98  error_message += "If you really need an assignment operator\n";
99  error_message += "for this class, write it yourself...\n";
100 
101  throw OomphLibError(error_message,"broken_assign()",
102  OOMPH_EXCEPTION_LOCATION);
103  }
104 
105 
106  /// Issue error message and terminate execution
107  void broken_copy(const std::string& class_name)
108  {
109  //Write the error message into a string
110  std::string error_message = "Copy constructor for class\n\n";
111  error_message += class_name;
112  error_message += "\n\n";
113  error_message += "is deliberately broken to avoid the accidental\n";
114  error_message += "use of the inappropriate C++ default.\n";
115  error_message +=
116  "All function arguments should be passed by reference or\n";
117  error_message +=
118  "constant reference. If you really need a copy constructor\n";
119  error_message += "for this class, write it yourself...\n";
120 
121  throw OomphLibError(error_message,
122  OOMPH_CURRENT_FUNCTION,
123  OOMPH_EXCEPTION_LOCATION);
124  }
125 }
126 
127 
128 //////////////////////////////////////////////////////////////////
129 //////////////////////////////////////////////////////////////////
130 //////////////////////////////////////////////////////////////////
131 
132 
133 
134 //====================================================================
135 /// Namespace for global (cumulative) timings
136 //====================================================================
137 namespace CumulativeTimings
138 {
139 
140  /// (Re-)start i-th timer
141  void start(const unsigned& i)
142  {
143  Start_time[i]=clock();
144  }
145 
146  /// Halt i-th timer
147  void halt(const unsigned& i)
148  {
149  Timing[i]+=clock()-Start_time[i];
150  }
151 
152  /// Report time accumulated by i-th timer
153  double cumulative_time(const unsigned& i)
154  {
155  return double(Timing[i])/CLOCKS_PER_SEC;
156  }
157 
158  /// Reset i-th timer
159  void reset(const unsigned& i)
160  {
161  Timing[i]=clock_t(0.0);
162  }
163 
164  /// Reset all timers
165  void reset()
166  {
167  unsigned n=Timing.size();
168  for (unsigned i=0;i<n;i++)
169  {
170  Timing[i]=clock_t(0.0);
171  }
172  }
173 
174  /// Set number of timings that can be recorded in parallel
175  void set_ntimers(const unsigned& ntimers)
176  {
177  Timing.resize(ntimers,clock_t(0.0));
178  Start_time.resize(ntimers,clock_t(0.0));
179  }
180 
181  /// Cumulative timings
183 
184  /// Start times of active timers
186 
187 }
188 
189 
190 //======================================================================
191 /// Namespace for black-box FD Newton solver.
192 //======================================================================
193 namespace BlackBoxFDNewtonSolver
194 {
195 
196 
197 /// Function pointer for residual function: Parameters, unknowns, residuals
198 typedef void (*ResidualFctPt)(const Vector<double>&,
199  const Vector<double>&,
200  Vector<double>&);
201 
202  /// Max. # of Newton iterations
203  unsigned Max_iter=20;
204 
205  /// Number of Newton iterations taken in most recent invocation
206  unsigned N_iter_taken=0;
207 
208  /// \short Flag to indicate if progress of Newton iteration is to be
209  /// documented (defaults to false)
210  bool Doc_Progress=false;
211 
212  /// FD step
213  double FD_step=1.0e-8;
214 
215  /// Tolerance
216  double Tol=1.0e-8;
217 
218  /// Use steplength control do make globally convergent (default false)
220 
221 /// \short Black-box FD Newton solver:
222 /// Calling sequence for residual function is
223 /// \code residual_fct(parameters,unknowns,residuals) \endcode
224 /// where all arguments are double Vectors.
225 /// unknowns.size() = residuals.size()
227  const Vector<double>& params,
228  Vector<double>& unknowns)
229 {
230  // Jacobian, current and advanced residual Vectors
231  unsigned ndof=unknowns.size();
232  DenseDoubleMatrix jacobian(ndof);
233  Vector<double> residuals(ndof);
234  Vector<double> residuals_pls(ndof);
235  Vector<double> dx(ndof);
236  Vector<double> gradient(ndof);
237  Vector<double> newton_direction(ndof);
238 
239  double half_residual_squared=0.0;
240  double max_step=0.0;
241 
242  /// Reset number of Newton iterations taken in most recent invocation
243  N_iter_taken=0;
244 
245  // Newton iterations
246  for (unsigned iloop=0;iloop<Max_iter;iloop++)
247  {
248  // Evaluate current residuals
249  residual_fct(params,unknowns,residuals);
250 
251  // Get half of squared residual and find maximum step length
252  // for step length control
254  {
255  half_residual_squared=0.0;
256  double sum=0.0;
257  for (unsigned i=0;i<ndof;i++)
258  {
259  sum += unknowns[i]*unknowns[i];
260  half_residual_squared+=residuals[i]*residuals[i];
261  }
262  half_residual_squared*=0.5;
263  max_step=100.0*std::max(sqrt(sum),double(ndof));
264  }
265 
266 
267  // Check max. residuals
268  double max_res = std::fabs(*std::max_element(residuals.begin(),
269  residuals.end(),
270  AbsCmp<double>()));
271 
272 
273  // Doc progress?
274  if (Doc_Progress)
275  {
276  oomph_info << "\nNewton iteration iter=" << iloop
277  << "\ni residual[i] unknown[i] " << std::endl;
278  for (unsigned i=0;i<ndof;i++)
279  {
280  oomph_info << i << " " << residuals[i]
281  << " " << unknowns[i] << std::endl;
282  }
283  }
284 
285 
286  // Converged?
287  if (max_res<Tol)
288  {
289  return;
290  }
291 
292  // Next iteration...
293  N_iter_taken++;
294 
295  // FD loop for Jacobian
296  for (unsigned i=0;i<ndof;i++)
297  {
298  double backup=unknowns[i];
299  unknowns[i]+=FD_step;
300 
301  // Evaluate advanced residuals
302  residual_fct(params,unknowns,residuals_pls);
303 
304  // Do FD
305  for (unsigned j=0;j<ndof;j++)
306  {
307  jacobian(j,i)=(residuals_pls[j]-residuals[j])/FD_step;
308  }
309 
310  // Reset fd step
311  unknowns[i]=backup;
312  }
313 
314  // oomph_info << "\n\nJacobian: "<< std::endl;
315  // jacobian.sparse_indexed_output(std::cout);
316  // oomph_info << std::endl;
317 
318  // Get gradient
320  {
321  for (unsigned i=0;i<ndof;i++)
322  {
323  double sum=0.0;
324  for (unsigned j=0;j<ndof;j++)
325  {
326  sum += jacobian(j,i)*residuals[j];
327  }
328  gradient[i]=sum;
329  }
330  }
331 
332  // Solve
333  jacobian.solve(residuals,newton_direction);
334 
335  // Update
337  {
338  for (unsigned i=0;i<ndof;i++)
339  {
340  newton_direction[i]*=-1.0;
341  }
342  // Update with steplength control
343  Vector<double> unknowns_old(unknowns);
344  double half_residual_squared_old=half_residual_squared;
345  line_search(unknowns_old,
346  half_residual_squared_old,
347  gradient,
348  residual_fct,
349  params,
350  newton_direction,
351  unknowns,
352  half_residual_squared,
353  max_step);
354  }
355  else
356  {
357  // Direct Newton update:
358  for (unsigned i=0;i<ndof;i++)
359  {
360  unknowns[i]-=newton_direction[i];
361  }
362 
363  }
364 
365 
366  }
367 
368 
369  // Failed to converge
370  std::ostringstream error_stream;
371  error_stream<< "Newton solver did not converge in "
372  << Max_iter << " steps " << std::endl;
373 
374  throw OomphLibError(error_stream.str(),
375  OOMPH_CURRENT_FUNCTION,
376  OOMPH_EXCEPTION_LOCATION);
377 }
378 
379 
380 
381 
382 
383 
384 
385 
386 
387 //=======================================================================
388 /// Line search helper for globally convergent Newton method
389 //=======================================================================
390  void line_search(const Vector<double>& x_old,
391  const double half_residual_squared_old,
392  const Vector<double>& gradient,
393  ResidualFctPt residual_fct,
394  const Vector<double>& params,
395  Vector<double>& newton_dir,
396  Vector<double>& x,
397  double& half_residual_squared,
398  const double& stpmax)
399  {
400  const double min_fct_decrease=1.0e-4;
401  double convergence_tol_on_x=1.0e-16;
402  double f_aux=0.0;
403  double lambda_aux=0.0;
404  double proposed_lambda=0.0;
405  unsigned n=x_old.size();
406  double sum=0.0;
407  for (unsigned i=0;i<n;i++)
408  {
409  sum += newton_dir[i]*newton_dir[i];
410  }
411  sum=sqrt(sum);
412  if (sum > stpmax)
413  {
414  for (unsigned i=0;i<n;i++)
415  {
416  newton_dir[i] *= stpmax/sum;
417  }
418  }
419  double slope=0.0;
420  for (unsigned i=0;i<n;i++)
421  {
422  slope += gradient[i]*newton_dir[i];
423  }
424  if (slope >= 0.0)
425  {
426  std::ostringstream error_stream;
427  error_stream << "Roundoff problem in lnsrch: slope=" << slope << "\n";
428  throw OomphLibError(error_stream.str(),
429  OOMPH_CURRENT_FUNCTION,
430  OOMPH_EXCEPTION_LOCATION);
431  }
432  double test=0.0;
433  for (unsigned i=0;i<n;i++)
434  {
435  double temp=std::fabs(newton_dir[i])/std::max(std::fabs(x_old[i]),1.0);
436  if (temp > test) test=temp;
437  }
438  double lambda_min=convergence_tol_on_x/test;
439  double lambda=1.0;
440  while (true)
441  {
442  for (unsigned i=0;i<n;i++)
443  {
444  x[i]=x_old[i]+lambda*newton_dir[i];
445  }
446 
447  // Evaluate current residuals
448  Vector<double> residuals(n);
449  residual_fct(params,x,residuals);
450  half_residual_squared=0.0;
451  for (unsigned i=0;i<n;i++)
452  {
453  half_residual_squared+=residuals[i]*residuals[i];
454  }
455  half_residual_squared*=0.5;
456 
457  if (lambda < lambda_min)
458  {
459  for (unsigned i=0;i<n;i++) x[i]=x_old[i];
460 
461  //Create an Oomph Lib warning
462  OomphLibWarning("Warning: Line search converged on x only!",
463  "BlackBoxFDNewtonSolver::line_search()",
464  OOMPH_EXCEPTION_LOCATION);
465  return;
466  }
467  else if (half_residual_squared <= half_residual_squared_old+
468  min_fct_decrease*lambda*slope)
469  {
470  return;
471  }
472  else
473  {
474  if (lambda == 1.0)
475  {
476  proposed_lambda = -slope/(2.0*(half_residual_squared-
477  half_residual_squared_old-slope));
478  }
479  else
480  {
481  double r1=half_residual_squared-
482  half_residual_squared_old-lambda*slope;
483  double r2=f_aux-half_residual_squared_old-lambda_aux*slope;
484  double a_poly=(r1/(lambda*lambda)-r2/(lambda_aux*lambda_aux))/
485  (lambda-lambda_aux);
486  double b_poly=(-lambda_aux*r1/(lambda*lambda)+lambda*r2/
487  (lambda_aux*lambda_aux))/(lambda-lambda_aux);
488  if (a_poly == 0.0)
489  {
490  proposed_lambda = -slope/(2.0*b_poly);
491  }
492  else
493  {
494  double discriminant=b_poly*b_poly-3.0*a_poly*slope;
495  if (discriminant < 0.0)
496  {
497  proposed_lambda=0.5*lambda;
498  }
499  else if (b_poly <= 0.0)
500  {
501  proposed_lambda=(-b_poly+sqrt(discriminant))/(3.0*a_poly);
502  }
503  else
504  {
505  proposed_lambda=-slope/(b_poly+sqrt(discriminant));
506  }
507  }
508  if (proposed_lambda>0.5*lambda)
509  {
510  proposed_lambda=0.5*lambda;
511  }
512  }
513  }
514  lambda_aux=lambda;
515  f_aux = half_residual_squared;
516  lambda=std::max(proposed_lambda,0.1*lambda);
517  }
518  }
519 }
520 
521 //======================================================================
522 /// \short Set output directory (we try to open a file in it
523 /// to see if the directory exists -- if it doesn't we'll
524 /// issue a warning -- or, if directory_must_exist()==true,
525 /// die by throwing and OomphLibError
526 //======================================================================
527 void DocInfo::set_directory(const std::string& directory_)
528 {
529  // Try to open a file in output directory
530  std::ostringstream filename;
531  filename << directory_ << "/.dummy_check.dat";
532  std::ofstream some_file;
533  some_file.open(filename.str().c_str());
534  if (!some_file.is_open())
535  {
536  //Construct the error message
537  std::string error_message = "Problem opening output file.\n";
538  error_message += "I suspect you haven't created the output directory ";
539  error_message += directory_;
540  error_message += "\n";
541 
542  //Issue a warning if the directory does not have to exist
544  {
545  //Create an Oomph Lib warning
546  OomphLibWarning(error_message,"set_directory()",
547  OOMPH_EXCEPTION_LOCATION);
548  }
549  //Otherwise throw an error
550  else
551  {
552  error_message += "and the Directory_must_exist flag is true.\n";
553  throw OomphLibError(error_message,
554  OOMPH_CURRENT_FUNCTION,
555  OOMPH_EXCEPTION_LOCATION);
556  }
557  }
558  //Write to the dummy file
559  some_file << "Dummy file, opened to check if output directory " << std::endl;
560  some_file << "exists. Can be deleted...." << std::endl;
561  some_file.close();
562  // Set directory
563  Directory=directory_;
564 }
565 
566 
567 // =================================================================
568 /// Conversion functions for easily making strings (e.g. for filenames - to
569 /// avoid stack smashing problems with cstrings and long filenames).
570 // =================================================================
571 namespace StringConversion
572 {
573  /// \short Convert a string to lower case (outputs a copy).
575  {
576  std::string output(input);
577  std::string::iterator it;
578  for(it = output.begin(); it!= output.end(); ++it)
579  { ::tolower(*it); }
580 
581  return output;
582  }
583 
584  /// \short Convert a string to upper case (outputs a copy).
586  {
587  std::string output(input);
588  std::string::iterator it;
589  for(it = output.begin(); it!= output.end(); ++it)
590  { ::toupper(*it); }
591  return output;
592  }
593 
594  /// \short Split a string, s, into a vector of strings where ever there is
595  /// an instance of delimiter (i.e. is delimiter is " " will give a list of
596  /// words). Note that mutliple delimiters in a row will give empty
597  /// strings.
598  void split_string(const std::string &s, char delim, Vector<std::string> &elems)
599  {
600  // From http://stackoverflow.com/questions/236129/splitting-a-string-in-c
601  std::stringstream ss(s);
602  std::string item;
603  while (std::getline(ss, item, delim))
604  {
605  elems.push_back(item);
606  }
607  }
608 
609  /// \short Split a string, s, into a vector of strings where ever there is
610  /// an instance of delimiter (i.e. is delimiter is " " will give a list of
611  /// words). Note that mutliple delimiters in a row will give empty
612  /// strings. Return by value.
614  {
615  // From http://stackoverflow.com/questions/236129/splitting-a-string-in-c
616  Vector<std::string> elems;
617  split_string(s, delim, elems);
618  return elems;
619  }
620 
621 }
622 
623 
624 
625 //====================================================================
626 /// Namespace for command line arguments
627 //====================================================================
628 namespace CommandLineArgs
629 {
630 
631  /// Number of arguments + 1
632  int Argc;
633 
634  /// Arguments themselves
635  char** Argv;
636 
637  /// Map to indicate an input flag as having been set
638  std::map<std::string, ArgInfo<bool> > Specified_command_line_flag;
639 
640  /// Map to associate an input flag with a double -- specified via pointer
641  std::map<std::string, ArgInfo<double> >
643 
644  /// Map to associate an input flag with an int -- specified via pointer
645  std::map<std::string, ArgInfo<int> >
647 
648  /// Map to associate an input flag with an unsigned -- specified via pointer
649  std::map<std::string, ArgInfo<unsigned> >
651 
652  /// Map to associate an input flag with a string -- specified via pointer
653  std::map<std::string, ArgInfo<std::string> >
655 
656  /// Set values
657  void setup(int argc, char** argv)
658  {
659  Argc=argc;
660  Argv=argv;
661  }
662 
663  /// Doc the command line arguments
664  void output()
665  {
666  oomph_info << "You are running the program: "
667  << CommandLineArgs::Argv[0] << std::endl;
668  oomph_info << "with the following command line args: " << std::endl;
669  std::stringstream str;
670  for (int i=1;i<CommandLineArgs::Argc;i++)
671  {
672  str << CommandLineArgs::Argv[i] << " ";
673  }
674  oomph_info << str.str() << std::endl;
675  }
676 
677 
678  /// Specify possible argument-free command line flag
679  void specify_command_line_flag(const std::string& command_line_flag,
680  const std::string& doc)
681  {
682  Specified_command_line_flag[command_line_flag] =
683  ArgInfo<bool>(false, 0, doc);
684  }
685 
686  /// \short Specify possible command line flag that specifies a double, accessed
687  /// via pointer
688  void specify_command_line_flag(const std::string& command_line_flag,
689  double* arg_pt,
690  const std::string& doc)
691  {
692  Specified_command_line_double_pt[command_line_flag] =
693  ArgInfo<double>(false, arg_pt, doc);
694  }
695 
696  /// \short Specify possible command line flag that specifies an int, accessed
697  /// via pointer
698  void specify_command_line_flag(const std::string& command_line_flag,
699  int* arg_pt,
700  const std::string& doc)
701  {
702  Specified_command_line_int_pt[command_line_flag] =
703  ArgInfo<int>(false, arg_pt, doc);
704  }
705 
706  /// \short Specify possible command line flag that specifies an unsigned,
707  /// accessed via pointer
708  void specify_command_line_flag(const std::string& command_line_flag,
709  unsigned* arg_pt,
710  const std::string& doc)
711  {
712  Specified_command_line_unsigned_pt[command_line_flag] =
713  ArgInfo<unsigned>(false, arg_pt, doc);
714 }
715 
716  /// \short Specify possible command line flag that specifies a string,
717  /// accessed via pointer
718  void specify_command_line_flag(const std::string& command_line_flag,
719  std::string* arg_pt,
720  const std::string& doc)
721  {
722  Specified_command_line_string_pt[command_line_flag] =
723  ArgInfo<std::string>(false, arg_pt, doc);
724  }
725 
726 
727  /// \short Check if command line flag has been set (value will have been
728  /// assigned directly).
730  {
731  for (std::map<std::string, ArgInfo<bool> >::iterator it=
733  it!=Specified_command_line_flag.end();it++)
734  {
735  if (it->first==flag)
736  {
737  return it->second.is_set;
738  }
739  }
740 
741  for (std::map<std::string, ArgInfo<double> >::iterator it=
743  it!=Specified_command_line_double_pt.end();it++)
744  {
745  if (it->first==flag)
746  {
747  return (it->second).is_set;
748  }
749  }
750 
751  for (std::map<std::string, ArgInfo<int> >::iterator it=
753  it!=Specified_command_line_int_pt.end();it++)
754  {
755  if (it->first==flag)
756  {
757  return (it->second).is_set;
758  }
759  }
760 
761  for (std::map<std::string, ArgInfo<unsigned> >::iterator it=
763  it!=Specified_command_line_unsigned_pt.end();it++)
764  {
765  if (it->first==flag)
766  {
767  return (it->second).is_set;
768  }
769  }
770 
771  for (std::map<std::string, ArgInfo<std::string> >::iterator it=
773  it!=Specified_command_line_string_pt.end();it++)
774  {
775  if (it->first==flag)
776  {
777  return (it->second).is_set;
778  }
779  }
780 
781  return false;
782  }
783 
784  /// Document the values of all flags (specified or not).
785  void doc_all_flags(std::ostream& outstream)
786  {
787  for (std::map<std::string, ArgInfo<bool> >::iterator it=
789  it!=Specified_command_line_flag.end();it++)
790  {
791  outstream << it->first << " " << it->second.is_set << std::endl;
792  }
793  for (std::map<std::string, ArgInfo<double> >::iterator it=
795  it!=Specified_command_line_double_pt.end();it++)
796  {
797  outstream << it->first << " " << *(it->second.arg_pt)
798  << std::endl;
799  }
800  for (std::map<std::string, ArgInfo<int> >::iterator it=
802  it!=Specified_command_line_int_pt.end();it++)
803  {
804  outstream << it->first << " " << *(it->second.arg_pt)
805  << std::endl;
806  }
807  for (std::map<std::string, ArgInfo<unsigned> >::iterator it=
809  it!=Specified_command_line_unsigned_pt.end();it++)
810  {
811  outstream << it->first << " " << *(it->second.arg_pt)
812  << std::endl;
813  }
814  for (std::map<std::string, ArgInfo<std::string> >::iterator it=
816  it!=Specified_command_line_string_pt.end();it++)
817  {
818  // Quote blank strings, otherwise trying to parse the output in any way
819  // will go wrong.
820  std::string arg_string = *(it->second.arg_pt);
821  if(arg_string == "")
822  {
823  arg_string = "\"\"";
824  }
825 
826  outstream << it->first << " " << arg_string
827  << std::endl;
828  }
829  }
830 
831  /// Document specified command line flags
833  {
834  oomph_info << std::endl;
835  oomph_info << "Specified (and recognised) command line flags:\n" ;
836  oomph_info << "----------------------------------------------\n";
837 
838  for (std::map<std::string, ArgInfo<bool> >::iterator it=
840  it!=Specified_command_line_flag.end();it++)
841  {
842  if (it->second.is_set)
843  {
844  oomph_info << it->first << std::endl;
845  }
846  }
847 
848  for (std::map<std::string, ArgInfo<double> >::iterator it=
850  it!=Specified_command_line_double_pt.end();it++)
851  {
852  if (it->second.is_set)
853  {
854  oomph_info << it->first << " "
855  << *(it->second.arg_pt)
856  << std::endl;
857  }
858  }
859 
860  for (std::map<std::string, ArgInfo<int> >::iterator it=
862  it!=Specified_command_line_int_pt.end();it++)
863  {
864  if (it->second.is_set)
865  {
866  oomph_info << it->first << " "
867  << *(it->second.arg_pt)
868  << std::endl;
869  }
870  }
871 
872  for (std::map<std::string, ArgInfo<unsigned> >::iterator it=
874  it!=Specified_command_line_unsigned_pt.end();it++)
875  {
876  if (it->second.is_set)
877  {
878  oomph_info << it->first << " "
879  << *(it->second.arg_pt)
880  << std::endl;
881  }
882  }
883 
884  for (std::map<std::string, ArgInfo<std::string> >::iterator it=
886  it!=Specified_command_line_string_pt.end();it++)
887  {
888  if (it->second.is_set)
889  {
890  oomph_info << it->first << " "
891  << *(it->second.arg_pt)
892  << std::endl;
893  }
894  }
895 
896  oomph_info << std::endl;
897  }
898 
899 
900  /// Document available command line flags
902  {
903  oomph_info << std::endl;
904  oomph_info << "Available command line flags:\n" ;
905  oomph_info << "-----------------------------\n";
906 
907  for (std::map<std::string, ArgInfo<bool> >::iterator it=
909  it!=Specified_command_line_flag.end();it++)
910  {
911  oomph_info << it->first << std::endl
912  << it->second.doc << std::endl
913  << std::endl;
914  }
915 
916  for (std::map<std::string, ArgInfo<double> >::iterator it=
918  it!=Specified_command_line_double_pt.end();it++)
919  {
920  oomph_info << it->first << " <double> " << std::endl
921  << it->second.doc << std::endl
922  << std::endl;
923  }
924 
925  for (std::map<std::string, ArgInfo<int> >::iterator it=
927  it!=Specified_command_line_int_pt.end();it++)
928  {
929  oomph_info << it->first << " <int> " << std::endl
930  << it->second.doc << std::endl
931  << std::endl;
932  }
933 
934  for (std::map<std::string, ArgInfo<unsigned> >::iterator it=
936  it!=Specified_command_line_unsigned_pt.end();it++)
937  {
938  oomph_info << it->first << " <unsigned> " << std::endl
939  << it->second.doc << std::endl
940  << std::endl;
941  }
942 
943  for (std::map<std::string, ArgInfo<std::string> >::iterator it=
945  it!=Specified_command_line_string_pt.end();it++)
946  {
947  oomph_info << it->first << " <string> " << std::endl
948  << it->second.doc << std::endl
949  << std::endl;
950  }
951 
952  oomph_info << std::endl;
953  }
954 
955 
956 
957  /// Helper function to check if command line index is legal
958  void check_arg_index(const int& argc,const int& arg_index)
959  {
960  if (arg_index>=argc)
961  {
962  output();
964  std::stringstream error_stream;
965  error_stream << "Tried to read more command line arguments than\n"
966  << "specified. This tends to happen if a required argument\n"
967  << "to a command line flag was omitted, e.g. by running \n\n"
968  << " ./a.out -some_double \n\n rather than\n\n"
969  << " ./a.out -some_double 1.23 \n\n"
970  << "To aid the debugging I've output the available\n"
971  << "command line arguments above.\n";
972  throw OomphLibError(
973  error_stream.str(),
974  OOMPH_CURRENT_FUNCTION,
975  OOMPH_EXCEPTION_LOCATION);
976  }
977  }
978 
979 
980 
981  /// \short Parse command line, check for recognised flags and assign
982  /// associated values
983  void parse_and_assign(int argc, char *argv[], const bool& throw_on_unrecognised_args)
984  {
985 
986  //Keep looping over command line arguments
987  int arg_index = 1;
988  while (arg_index < argc)
989  {
990 
991  bool found_match=false;
992 
993  if ( (strcmp(argv[arg_index], "-help") == 0 ) ||
994  (strcmp(argv[arg_index], "--help")== 0 ) )
995  {
996  oomph_info << "NOTE: You entered --help or -help on the command line\n";
997  oomph_info << "so I'm going to tell you about this code's\n";
998  oomph_info << "available command line flags and then return.\n";
1000 
1001 #ifdef OOMPH_HAS_MPI
1002  int flag;
1003  MPI_Initialized(&flag);
1004  if (flag!=0) MPI_Helpers::finalize();
1005 #endif
1006  oomph_info << "Shutting down...\n";
1007  exit(0);
1008 
1009  }
1010 
1011 
1012  //Check if the flag has been previously specified as a simple argument free
1013  //command line argument
1014  for (std::map<std::string, ArgInfo<bool> >::iterator it=
1016  it!=Specified_command_line_flag.end();it++)
1017  {
1018  if (it->first==argv[arg_index])
1019  {
1020  Specified_command_line_flag[argv[arg_index]].is_set = true;
1021  found_match=true;
1022  break;
1023  }
1024  }
1025 
1026  if (!found_match)
1027  {
1028  //Check if the flag has been previously specified as a
1029  //command line argument that specifies (and is followed by) a double
1030  for (std::map<std::string, ArgInfo<double> >::iterator it=
1032  it!=Specified_command_line_double_pt.end();it++)
1033  {
1034  if (it->first==argv[arg_index])
1035  {
1036  //Next command line argument specifies the double. Set it via
1037  //the pointer
1038  arg_index++;
1039  check_arg_index(argc,arg_index);
1040  it->second.is_set=true;
1041  *(it->second.arg_pt)=atof(argv[arg_index]);
1042  found_match=true;
1043  break;
1044  }
1045  }
1046  }
1047 
1048 
1049  if (!found_match)
1050  {
1051  //Check if the flag has been previously specified as a
1052  //command line argument that specifies (and is followed by) an int
1053  for (std::map<std::string, ArgInfo<int> >::iterator it=
1055  it!=Specified_command_line_int_pt.end();it++)
1056  {
1057  if (it->first==argv[arg_index])
1058  {
1059  //Next command line argument specifies the int. Set it via
1060  //the pointer
1061  arg_index++;
1062  check_arg_index(argc,arg_index);
1063  it->second.is_set=true;
1064  *(it->second.arg_pt)=atoi(argv[arg_index]);
1065  found_match=true;
1066  break;
1067  }
1068  }
1069  }
1070 
1071 
1072  if (!found_match)
1073  {
1074  //Check if the flag has been previously specified as a
1075  //command line argument that specifies (and is followed by) an unsigned
1076  for (std::map<std::string, ArgInfo<unsigned> >::iterator it=
1078  it!=Specified_command_line_unsigned_pt.end();it++)
1079  {
1080  if (it->first==argv[arg_index])
1081  {
1082  //Next command line argument specifies the unsigned. Set it via
1083  //the pointer
1084  arg_index++;
1085  check_arg_index(argc,arg_index);
1086  it->second.is_set=true;
1087  *(it->second.arg_pt)=unsigned(atoi(argv[arg_index]));
1088  found_match=true;
1089  break;
1090  }
1091  }
1092  }
1093 
1094 
1095  if (!found_match)
1096  {
1097  //Check if the flag has been previously specified as a
1098  //command line argument that specifies (and is followed by) a string
1099  for (std::map<std::string, ArgInfo<std::string> >::iterator it=
1101  it!=Specified_command_line_string_pt.end();it++)
1102  {
1103  if (it->first==argv[arg_index])
1104  {
1105  //Next command line argument specifies the string. Set it via
1106  //the pointer
1107  arg_index++;
1108  check_arg_index(argc,arg_index);
1109  it->second.is_set=true;
1110  *(it->second.arg_pt)=argv[arg_index];
1111  found_match=true;
1112  break;
1113  }
1114  }
1115  }
1116 
1117 
1118  // Oh dear, we still haven't found the argument in the list.
1119  // Maybe it was specified wrongly -- issue warning.
1120  if (!found_match)
1121  {
1122 
1123  //Construct the error message
1124  std::string error_message = "Command line argument\n\n";
1125  error_message += argv[arg_index];
1126  error_message += "\n\nwas not recognised. This may be legal\n";
1127  error_message += "but seems sufficiently suspicious to flag up.\n";
1128  error_message += "If it should have been recognised, make sure you\n";
1129  error_message += "used the right number of \"-\" signs...\n";
1130 
1131  if(throw_on_unrecognised_args)
1132  {
1133  error_message += "Throwing an error because you requested it with";
1134  error_message += " throw_on_unrecognised_args option.";
1135  throw OomphLibError(error_message, OOMPH_CURRENT_FUNCTION,
1136  OOMPH_EXCEPTION_LOCATION);
1137  }
1138  else
1139  {
1140  //Create an Oomph Lib warning
1141  OomphLibWarning(error_message,OOMPH_CURRENT_FUNCTION,
1142  OOMPH_EXCEPTION_LOCATION);
1143  }
1144  }
1145 
1146 
1147 
1148  arg_index++;
1149  }
1150  }
1151 
1152 
1153 
1154  /// \short Parse previously specified command line, check for
1155  /// recognised flags and assign associated values
1156  void parse_and_assign(const bool& throw_on_unrecognised_args)
1157  {
1158  parse_and_assign(Argc, Argv, throw_on_unrecognised_args);
1159  }
1160 
1161 }
1162 
1163 #ifdef OOMPH_HAS_MPI
1164 //========================================================================
1165 /// Single (global) instantiation of the mpi output modifier
1166 //========================================================================
1168 
1169 //========================================================================
1170 /// Precede the output by the processor ID but output everything
1171 //========================================================================
1172  bool MPIOutputModifier::operator()(std::ostream &stream)
1173  {
1174  int my_rank = Communicator_pt->my_rank();
1175 
1177  {
1178  stream << "Processor " << my_rank << ": ";
1179  // Continue processing
1180  return true;
1181  }
1182  else
1183  {
1184  if (unsigned(my_rank)==Output_rank)
1185  {
1186  stream << "Processor " << my_rank << ": ";
1187  // Continue processing
1188  return true;
1189  }
1190  else
1191  {
1192  return false;
1193  }
1194  }
1195  }
1196 
1197 #endif
1198 
1199 //=============================================================================
1200 /// Initialize mpi. If optional boolean flag is set to false, we use
1201 /// MPI_COMM_WORLD itself as oomph-lib's communicator. Defaults to true.
1202 //=============================================================================
1203  void MPI_Helpers::init(int argc, char **argv,
1204  const bool& make_duplicate_of_mpi_comm_world)
1205 {
1206 #ifdef OOMPH_HAS_MPI
1207  // call mpi int
1208  MPI_Init(&argc,&argv);
1209 
1210 
1211  // By default, create the oomph-lib communicator using MPI_Comm_dup so that
1212  // the communicator has the same group of processes but a new context
1213  MPI_Comm oomph_comm_world=MPI_COMM_WORLD;
1214  if (make_duplicate_of_mpi_comm_world)
1215  {
1216  MPI_Comm_dup(MPI_COMM_WORLD,&oomph_comm_world);
1217  }
1218 
1219  if (MPI_COMM_WORLD!=oomph_comm_world)
1220  {
1221  oomph_info << "Oomph-lib communicator is a duplicate of MPI_COMM_WORLD\n";
1222  }
1223  else
1224  {
1225  oomph_info << "Oomph-lib communicator is MPI_COMM_WORLD\n";
1226  }
1227 
1228  // create the oomph-lib communicator
1229  // note: oomph_comm_world is deleted when the destructor of
1230  // Communicator_pt is called
1231  Communicator_pt = new OomphCommunicator(oomph_comm_world,true);
1232 
1233  // Change MPI error handler so that error will return
1234  // rather than aborting
1235  MPI_Errhandler_set(oomph_comm_world, MPI_ERRORS_RETURN);
1236 
1237  // Use MPI output modifier: Each processor precedes its output
1238  // by its rank
1241 #else
1242  // create a serial communicator
1244 #endif
1246 }
1247 
1248 //=============================================================================
1249 /// finalize mpi
1250 //=============================================================================
1252 {
1253  // delete the communicator
1254  delete Communicator_pt;
1255 
1256  // and call MPI_Finalize
1257 #ifdef OOMPH_HAS_MPI
1258  MPI_Finalize();
1259 #endif
1260 }
1261 
1262 //=============================================================================
1263 /// access to the global oomph-lib communicator
1264 //=============================================================================
1266 {
1267 #ifdef MPI
1268 #ifdef PARANOID
1270  {
1271  std::ostringstream error_message_stream;
1272  error_message_stream
1273  << "MPI has not been initialised.\n Call MPI_Helpers::init(...)";
1274  throw OomphLibError(error_message_stream.str(),
1275  OOMPH_CURRENT_FUNCTION,
1276  OOMPH_EXCEPTION_LOCATION);
1277  }
1278 #endif
1279  return Communicator_pt;
1280 
1281 #else // ifndef MPI
1282 
1283  // If MPI is not enabled then we are unlikely to have called init, so to
1284  // be nice to users create a new serial comm_pt if none exits.
1285 
1287 
1288 // #ifdef PARANOID
1289 // if(!Communicator_pt->serial_communicator())
1290 // {
1291 // std::string error_msg =
1292 // "MPI_Helpers has somehow ended up with a non-serial\n"
1293 // + "communicator pointer even though MPI is disabled!";
1294 // throw OomphLibError(error_msg.str(),
1295 // OOMPH_CURRENT_FUNCTION,
1296 // OOMPH_EXCEPTION_LOCATION);
1297 // }
1298 // #endif
1299 
1300  return Communicator_pt;
1301 
1302 #endif // end ifdef MPI
1303 }
1304 
1307 
1308 
1309 //====================================================================
1310 /// Namespace for flagging up obsolete parts of the code
1311 //====================================================================
1312 namespace ObsoleteCode
1313 {
1314 
1315  /// Flag up obsolete parts of the code
1316  bool FlagObsoleteCode=true;
1317 
1318  /// Output warning message
1319  void obsolete()
1320  {
1321  if (FlagObsoleteCode)
1322  {
1323  std::string junk;
1324  oomph_info << "\n\n--------------------------------------------\n";
1325  oomph_info << "You are using obsolete code " << std::endl;
1326  oomph_info << "--------------------------------------------\n\n";
1327  oomph_info << "Enter: \"s\" to suppress further messages" << std::endl;
1328  oomph_info <<
1329  " \"k\" to crash the code to allow a trace back in the debugger"
1330  << std::endl;
1331 
1332  oomph_info << " any other key to continue\n \n";
1333  oomph_info <<
1334  " [Note: Insert \n \n ";
1335  oomph_info <<
1336  " ObsoleteCode::FlagObsoleteCode=false;\n \n";
1337  oomph_info <<
1338  " into your code to suppress these messages \n";
1339  oomph_info <<
1340  " altogether.] \n";
1341 
1342  std::cin >> junk;
1343  if (junk=="s")
1344  {
1345  FlagObsoleteCode=false;
1346  }
1347  if (junk=="k")
1348  {
1349  throw OomphLibError("Killed",
1350  OOMPH_CURRENT_FUNCTION,
1351  OOMPH_EXCEPTION_LOCATION);
1352  }
1353  }
1354  }
1355 
1356 
1357  ///Ouput a warning message with a string argument
1358  void obsolete(const std::string &message)
1359  {
1360  if(FlagObsoleteCode)
1361  {
1362  oomph_info << "\n\n------------------------------------" << std::endl;
1363  oomph_info << message << std::endl;
1364  oomph_info << "----------------------------------------" << std::endl;
1365 
1366  obsolete();
1367  }
1368  }
1369 
1370 }
1371 
1372 
1373 
1374 
1375 //====================================================================
1376 /// Namespace for tecplot stuff
1377 //====================================================================
1378 namespace TecplotNames
1379 {
1380 
1381  /// Tecplot colours
1383 
1384 
1385  /// Setup namespace
1386  void setup()
1387  {
1388  colour.resize(5);
1389  colour[0]="RED";
1390  colour[1]="GREEN";
1391  colour[2]="BLUE";
1392  colour[3]="CYAN";
1393  colour[4]="BLACK";
1394  }
1395 
1396 
1397 }
1398 
1399 
1400 
1401 #ifdef LEAK_CHECK
1402 
1403 //====================================================================
1404 /// Namespace for leak check: Keep a running count of all instantiated
1405 /// objects -- add your own if you want to...
1406 //====================================================================
1407 namespace LeakCheckNames
1408 {
1409 
1421 
1422  void reset()
1423  {
1424  QuadTree_build=0;
1425  OcTree_build=0;
1429  RefineableQElement<3>_build=0;
1431  HangInfo_build=0;
1432  Node_build=0;
1435  }
1436 
1437  void doc()
1438  {
1439  oomph_info <<
1440  "\n Leak check: # of builds - # of deletes for the following objects: \n\n";
1441  oomph_info << "LeakCheckNames::QuadTree_build "
1442  << LeakCheckNames::QuadTree_build << std::endl;
1443  oomph_info << "LeakCheckNames::QuadTreeForest_build "
1444  << LeakCheckNames::QuadTreeForest_build << std::endl;
1445  oomph_info << "LeakCheckNames::OcTree_build "
1446  << LeakCheckNames::OcTree_build << std::endl;
1447  oomph_info << "LeakCheckNames::OcTreeForest_build "
1448  << LeakCheckNames::OcTreeForest_build << std::endl;
1449  oomph_info << "LeakCheckNames::RefineableQElement<2>_build "
1450  << LeakCheckNames::RefineableQElement<2>_build << std::endl;
1451  oomph_info << "LeakCheckNames::RefineableQElement<3>_build "
1452  << LeakCheckNames::RefineableQElement<3>_build << std::endl;
1453  oomph_info << "LeakCheckNames::MacroElement_build "
1454  << LeakCheckNames::MacroElement_build<< std::endl;
1455  oomph_info << "LeakCheckNames::HangInfo_build "
1456  << LeakCheckNames::HangInfo_build<< std::endl;
1457  oomph_info << "LeakCheckNames::Node_build "
1458  << LeakCheckNames::Node_build<< std::endl;
1459  oomph_info << "LeakCheckNames::GeomReference_build "
1460  << LeakCheckNames::GeomReference_build<< std::endl;
1461  oomph_info << "LeakCheckNames::AlgebraicNode_build "
1462  << LeakCheckNames::AlgebraicNode_build<< std::endl;
1463  oomph_info << std::endl;
1464  }
1465 
1466 
1467 
1468 }
1469 
1470 
1471 #endif
1472 
1473 
1474 
1475 ///////////////////////////////////////////////////////////////////////////
1476 /////////////////////////////////////////////////////////////////////////////
1477 
1478 
1479 //====================================================================
1480 /// Namespace for pause() command
1481 //====================================================================
1482 namespace PauseFlags
1483 {
1484 
1485  /// Flag to enable pausing code -- pause the code by default
1486  bool PauseFlag=true;
1487 
1488 }
1489 
1490 //======================================================================
1491 /// Pause and display message
1492 //======================================================================
1493 void pause(std::string message)
1494 {
1495  std::string junk;
1497  {
1498  oomph_info << message
1499  << "\n hit any key to continue [hit \"S\" "
1500  << "to suppress further interruptions]\n";
1501  std::cin >> junk;
1502  if (junk=="S")
1503  {
1504  PauseFlags::PauseFlag=false;
1505  }
1506  }
1507  else
1508  {
1509  oomph_info << "\n[Suppressed pause message] \n";
1510  }
1511  }
1512 
1513 
1514 
1515 
1516 
1517 ///////////////////////////////////////////////////////////////////////////
1518 ///////////////////////////////////////////////////////////////////////////
1519 ///////////////////////////////////////////////////////////////////////////
1520 
1521 //=============================================================================
1522 /// Helper for recording execution time.
1523 //=============================================================================
1524 namespace TimingHelpers
1525 {
1526 
1527  /// returns the time in seconds after some point in past
1528  double timer()
1529  {
1530 #ifdef OOMPH_HAS_MPI
1532  {
1533  return MPI_Wtime();
1534  }
1535  else
1536 #endif
1537  {
1538  time_t t = clock();
1539  return double(t) / double(CLOCKS_PER_SEC);
1540  }
1541  }
1542 }//end of namespace TimingHelpers
1543 
1544 
1545 ////////////////////////////////////////////////////////////////////
1546 ////////////////////////////////////////////////////////////////////
1547 ////////////////////////////////////////////////////////////////////
1548 
1549 
1550 
1551 //===============================================================
1552 /// Namespace with helper functions to assess total memory usage
1553 /// on the fly using system() -- details are very machine specific! This just
1554 /// provides the overall machinery with default settings for
1555 /// our own (linux machines). Uses the system command
1556 /// to spawn a command that computes the total memory usage
1557 /// on the machine where this is called. [Disclaimer: works
1558 /// on my machine(s) -- no guarantees for any other platform;
1559 /// linux or not. MH]
1560 //===============================================================
1561 namespace MemoryUsage
1562 {
1563 
1564  /// \short Boolean to suppress synchronisation of doc memory
1565  /// usage on processors (via mpi barriers). True (i.e. sync is
1566  /// suppressed) by default because not all processors may
1567  /// reach the relevant doc memory usage statements
1568  /// causing the code to hang).
1570 
1571  /// \short String containing system command that obtains memory usage
1572  /// of all processes.
1573  /// Default assignment for linux. [Disclaimer: works on my machine(s) --
1574  /// no guarantees for any other platform; linux or not. MH]
1576 
1577  /// \short Bool allowing quick bypassing of ALL operations related
1578  /// to memory usage monitoring -- this allows the code to remain
1579  /// "instrumented" without incurring the heavy penalties associated
1580  /// with the system calls and i/o. Default setting: false.
1582 
1583  /// \short String containing name of file in which we document
1584  /// my memory usage -- you may want to change this to allow different
1585  /// processors to write to separate files (especially in mpi
1586  /// context). Note that file is appended to
1587  /// so it ought to be emptied (either manually or by calling
1588  /// helper function empty_memory_usage_file()
1589  std::string My_memory_usage_filename="my_memory_usage.dat";
1590 
1591  /// \short Function to empty file that records my memory usage in
1592  /// file whose name is specified by My_memory_usage_filename.
1594  {
1595  // bail out straight away?
1597 
1598  // Open without appending and write header
1599  std::ofstream the_file;
1600  the_file.open(My_memory_usage_filename.c_str());
1601  the_file << "# My memory usage: \n";
1602  the_file.close();
1603  }
1604 
1605 
1606 #ifdef OOMPH_HAS_UNISTDH
1607 
1608  /// \short Doc my memory usage, prepended by string (which allows
1609  /// identification from where the function is called, say) that records
1610  /// memory usage in file whose name is specified by My_memory_usage_filename.
1611  /// Data is appended to that file; wipe it with empty_my_memory_usage_file().
1612  /// Note: This requires getpid() which is defined in unistd.h so if you
1613  /// don't have that we won't build that function!
1614  void doc_my_memory_usage(const std::string& prefix_string)
1615  {
1616  // bail out straight away?
1618 
1619  // Write prefix
1620  std::ofstream the_file;
1621  the_file.open(My_memory_usage_filename.c_str(),std::ios::app);
1622  the_file << prefix_string << " ";
1623  the_file.close();
1624 
1625  // Sync all processors if in parallel (unless suppressed)
1626 #ifdef OOMPH_HAS_MPI
1629  {
1630  MPI_Barrier(MPI_Helpers::communicator_pt()->mpi_comm());
1631  }
1632 #endif
1633 
1634  // Process memory usage command and write to file
1635  std::stringstream tmp;
1637  << " | awk '{if ($2==" << getpid() << "){print $4}}' >> "
1639  int success=system(tmp.str().c_str());
1640 
1641  // Dummy command to shut up compiler warnings
1642  success+=1;
1643  }
1644 
1645 #endif
1646 
1647  /// \short String containing system command that obtains total memory usage.
1648  /// Default assignment for linux. [Disclaimer: works on my machine(s) --
1649  /// no guarantees for any other platform; linux or not. MH]
1651  "ps aux | awk 'BEGIN{sum=0}{sum+=$4}END{print sum}'";
1652 
1653  /// \short String containing name of file in which we document total
1654  /// memory usage -- you may want to change this to allow different
1655  /// processors to write to separate files (especially in mpi
1656  /// context). Note that file is appended to
1657  /// so it ought to be emptied (either manually or by calling
1658  /// helper function empty_memory_usage_file()
1660 
1661  /// \short Function to empty file that records total memory usage in
1662  /// file whose name is specified by Total_memory_usage_filename.
1664  {
1665  // bail out straight away?
1667 
1668  // Open without appending and write header
1669  std::ofstream the_file;
1670  the_file.open(Total_memory_usage_filename.c_str());
1671  the_file << "# Total memory usage: \n";
1672  the_file.close();
1673  }
1674 
1675  /// \short Doc total memory usage, prepended by string (which allows
1676  /// identification from where the function is called, say) that records
1677  /// memory usage in file whose name is specified by Total_memory_usage_filename.
1678  /// Data is appended to that file; wipe it with empty_memory_usage_file().
1679  void doc_total_memory_usage(const std::string& prefix_string)
1680  {
1681  // bail out straight away?
1683 
1684  // Write prefix
1685  std::ofstream the_file;
1686  the_file.open(Total_memory_usage_filename.c_str(),std::ios::app);
1687  the_file << prefix_string << " ";
1688  the_file.close();
1689 
1690  // Sync all processors if in parallel
1691 #ifdef OOMPH_HAS_MPI
1694  {
1695  MPI_Barrier(MPI_Helpers::communicator_pt()->mpi_comm());
1696  }
1697 #endif
1698 
1699  // Process memory usage command and write to file
1700  std::stringstream tmp;
1701  tmp << Total_memory_usage_system_string << " >> "
1703  int success=system(tmp.str().c_str());
1704 
1705  // Dummy command to shut up compiler warnings
1706  success+=1;
1707  }
1708 
1709 
1710 
1711  /// \short Function to empty file that records total and local memory usage
1712  /// in appropriate files
1714  {
1715  // bail out straight away?
1717 
1720  empty_top_file();
1721  }
1722 
1723  /// \short Doc total and local memory usage, prepended by string (which allows
1724  /// identification from where the function is called, say). NOTE: Local
1725  /// memory usage only works if we have unistd.h header
1726  void doc_memory_usage(const std::string& prefix_string)
1727  {
1728  // bail out straight away?
1730 
1731 #ifdef OOMPH_HAS_UNISTDH
1732  doc_my_memory_usage(prefix_string);
1733 #endif
1734 
1735  doc_total_memory_usage(prefix_string);
1736  }
1737 
1738 
1739 
1740 
1741  /// \short String containing system command that runs "top" (or equivalent)
1742  /// "indefinitely" and writes to file specified in Top_output_filename.
1743  /// Default assignment for linux. [Disclaimer: works on my machine(s) --
1744  /// no guarantees for any other platform; linux or not. MH]
1745  std::string Top_system_string="while true; do top -b -n 2 ; done ";
1746 
1747  /// \short String containing name of file in which we document "continuous"
1748  /// output from "top" (or equivalent)-- you may want to change this to
1749  /// allow different processors to write to separate files (especially in mpi
1750  /// context). Note that file is appended to
1751  /// so it ought to be emptied (either manually or by calling
1752  /// helper function empty_top_file()
1754 
1755  /// \short Function to empty file that records continuous output from top in
1756  /// file whose name is specified by Top_output_filename
1758  {
1759  // bail out straight away?
1761 
1762  // Open without appending and write header
1763  std::ofstream the_file;
1764  the_file.open(Top_output_filename.c_str());
1765  the_file << "# Continuous output from top obtained with: \n";
1766  the_file << "# " << Top_system_string << "\n";
1767  the_file.close();
1768  }
1769 
1770 
1771  /// \short Start running top continuously and output (append) into
1772  /// file specified by Top_output_filename. Wipe that file with
1773  /// empty_top_file() if you wish. Note that this is again quite Linux specific
1774  /// and unlikely to work on other operating systems.
1775  /// Insert optional comment into output file before starting.
1776  void run_continous_top(const std::string& comment)
1777  {
1778  // bail out straight away?
1780 
1781  // Sync all processors if in parallel
1782  std::string modifier="";
1783 
1784 #ifdef OOMPH_HAS_MPI
1786  {
1788  {
1789  MPI_Barrier(MPI_Helpers::communicator_pt()->mpi_comm());
1790  }
1791  std::stringstream tmp;
1792  tmp << "_proc" << MPI_Helpers::communicator_pt()->my_rank();
1793  modifier=tmp.str();
1794  }
1795 #endif
1796 
1797  // Process memory usage command and write to file
1798  std::stringstream tmp;
1799 
1800  // String stream seems unhappy about direct insertions of these
1801  std::string backslash="\\";
1802  std::string dollar="$";
1803 
1804  // Create files that spawn and kill continuous top
1805  tmp << "echo \"#/bin/bash\" > run_continuous_top" << modifier
1806  << ".bash; "
1807  << "echo \" echo " << backslash << "\" kill "
1808  << backslash << dollar
1809  << backslash << dollar << " " << backslash << "\" > kill_continuous_top"
1810  << modifier << ".bash; chmod a+x kill_continuous_top"
1811  << modifier << ".bash; " << Top_system_string
1812  << " \" >> run_continuous_top" << modifier
1813  << ".bash; chmod a+x run_continuous_top" << modifier << ".bash ";
1814  int success=system(tmp.str().c_str());
1815 
1816  // Add comment to annotate?
1817  if (comment!="")
1818  {
1820  }
1821 
1822  // Start spawning
1823  std::stringstream tmp2;
1824  tmp2 << "./run_continuous_top" << modifier << ".bash >> "
1825  << Top_output_filename << " & ";
1826  success=system(tmp2.str().c_str());
1827 
1828  // Dummy command to shut up compiler warnings
1829  success+=1;
1830  }
1831 
1832 
1833 
1834  /// \short Stop running top continuously. Note that this is
1835  /// again quite Linux specific and unlikely to work on other operating
1836  /// systems. Insert optional comment into output file before stopping.
1837  void stop_continous_top(const std::string& comment)
1838  {
1839  // bail out straight away?
1841 
1842  // Sync all processors if in parallel
1843  std::string modifier="";
1844 
1845 #ifdef OOMPH_HAS_MPI
1847  {
1849  {
1850  MPI_Barrier(MPI_Helpers::communicator_pt()->mpi_comm());
1851  }
1852  std::stringstream tmp;
1853  tmp << "_proc" << MPI_Helpers::communicator_pt()->my_rank();
1854  modifier=tmp.str();
1855  }
1856 #endif
1857 
1858  // Add comment to annotate?
1859  if (comment!="")
1860  {
1862  }
1863 
1864  // Kill
1865  std::stringstream tmp2;
1866  tmp2 << "./kill_continuous_top" << modifier << ".bash >> "
1867  << Top_output_filename << " & ";
1868  int success=system(tmp2.str().c_str());
1869 
1870  // Dummy command to shut up compiler warnings
1871  success+=1;
1872  }
1873 
1874 
1875  /// \short Insert comment into running continuous top output
1877  {
1878  // bail out straight away?
1880 
1881  std::stringstream tmp;
1882  tmp << " echo \"OOMPH-LIB EVENT: " << comment << "\" >> "
1884  int success=system(tmp.str().c_str());
1885 
1886  // Dummy command to shut up compiler warnings
1887  success+=1;
1888  }
1889 
1890 
1891 
1892 } // end of namespace MemoryUsage
1893 
1894 
1895 }
void doc_total_memory_usage(const std::string &prefix_string)
Doc total memory usage, prepended by string (which allows identification from where the function is c...
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void specify_command_line_flag(const std::string &command_line_flag, const std::string &doc)
Specify possible argument-free command line flag.
char ** Argv
Arguments themselves.
std::string to_lower(const std::string &input)
Convert a string to lower case (outputs a copy).
Structure to store information on a command line argument.
Vector< clock_t > Timing
Cumulative timings.
bool Doc_Progress
Flag to indicate if progress of Newton iteration is to be documented (defaults to false) ...
std::string My_memory_usage_system_string
String containing system command that obtains memory usage of all processes. Default assignment for l...
cstr elem_len * i
Definition: cfortran.h:607
std::string Top_system_string
String containing system command that runs "top" (or equivalent) "indefinitely" and writes to file sp...
void output()
Doc the command line arguments.
std::map< std::string, ArgInfo< int > > Specified_command_line_int_pt
Map to associate an input flag with an int – specified via pointer.
void check_arg_index(const int &argc, const int &arg_index)
Helper function to check if command line index is legal.
void obsolete()
Output warning message.
std::string Directory
Directory name.
static void finalize()
finalize mpi
void run_continous_top(const std::string &comment)
Start running top continuously and output (append) into file specified by Top_output_filename. Wipe that file with empty_top_file() if you wish. Note that this is again quite Linux specific and unlikely to work on other operating systems. Insert optional comment into output file before starting.
char t
Definition: cfortran.h:572
std::map< std::string, ArgInfo< std::string > > Specified_command_line_string_pt
Map to associate an input flag with a string – specified via pointer.
unsigned N_iter_taken
Number of Newton iterations taken in most recent invocation.
void empty_total_memory_usage_file()
Function to empty file that records total memory usage in file whose name is specified by Total_memor...
static OomphCommunicator * communicator_pt()
access to the global oomph-lib communicator
OomphInfo oomph_info
unsigned Output_rank
Rank of single processor that produces output (only used if Output_from_single_processor=true.
static bool mpi_has_been_initialised()
return true if MPI has been initialised
OomphCommunicator *& communicator_pt()
Return pointer to communicator.
std::map< std::string, ArgInfo< bool > > Specified_command_line_flag
Map to indicate an input flag as having been set.
void set_ntimers(const unsigned &ntimers)
Set number of timings that can be recorded in parallel.
bool Use_step_length_control
Use steplength control do make globally convergent (default false)
void empty_memory_usage_files()
Function to empty file that records total and local memory usage in appropriate files.
std::map< std::string, ArgInfo< unsigned > > Specified_command_line_unsigned_pt
Map to associate an input flag with an unsigned – specified via pointer.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
OomphCommunicator * Communicator_pt
Communicator.
std::string to_upper(const std::string &input)
Convert a string to upper case (outputs a copy).
void insert_comment_to_continous_top(const std::string &comment)
Insert comment into running continuous top output.
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)...
Definition: matrices.cc:55
static OomphCommunicator * Communicator_pt
the global communicator
std::map< std::string, ArgInfo< double > > Specified_command_line_double_pt
Map to associate an input flag with a double – specified via pointer.
void reset(const unsigned &i)
Reset i-th timer.
virtual bool operator()(std::ostream &stream)
Precede the output by the processor ID but output everything.
static bool MPI_has_been_initialised
Bool set to true if MPI has been initialised.
void start(const unsigned &i)
(Re-)start i-th timer
void setup()
Setup namespace.
void split_string(const std::string &s, char delim, Vector< std::string > &elems)
Split a string, s, into a vector of strings where ever there is an instance of delimiter (i...
static char t char * s
Definition: cfortran.h:572
std::string Top_output_filename
String containing name of file in which we document "continuous" output from "top" (or equivalent)– y...
void parse_and_assign(int argc, char *argv[], const bool &throw_on_unrecognised_args)
Parse command line, check for recognised flags and assign associated values.
unsigned Max_iter
Max. # of Newton iterations.
std::string My_memory_usage_filename
String containing name of file in which we document my memory usage – you may want to change this to ...
void setup(int argc, char **argv)
Set values.
MPIOutputModifier oomph_mpi_output
Single (global) instantiation of the mpi output modifier.
double timer()
returns the time in seconds after some point in past
void output(std::ostream &outfile)
Output with default number of plot points.
bool command_line_flag_has_been_set(const std::string &flag)
Check if command line flag has been set (value will have been assigned directly). ...
void empty_my_memory_usage_file()
Function to empty file that records my memory usage in file whose name is specified by My_memory_usag...
void doc_my_memory_usage(const std::string &prefix_string)
Doc my memory usage, prepended by string (which allows identification from where the function is call...
void black_box_fd_newton_solve(ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &unknowns)
Black-box FD Newton solver: Calling sequence for residual function is.
static void init(int argc, char **argv, const bool &make_duplicate_of_mpi_comm_world=true)
initialise mpi (oomph-libs equivalent of MPI_Init(...)) Initialises MPI and creates the global oomph-...
int Argc
Number of arguments + 1.
OutputModifier *& output_modifier_pt()
Access function for the output modifier pointer.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
std::string Total_memory_usage_filename
String containing name of file in which we document total memory usage – you may want to change this ...
bool Suppress_mpi_synchronisation
Boolean to suppress synchronisation of doc memory usage on processors (via mpi barriers). True (i.e. sync is suppressed) by default because not all processors may reach the relevant doc memory usage statements causing the code to hang).
bool PauseFlag
Flag to enable pausing code – pause the code by default.
Vector< clock_t > Start_time
Start times of active timers.
bool Bypass_all_memory_usage_monitoring
Bool allowing quick bypassing of ALL operations related to memory usage monitoring – this allows the ...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void set_directory(const std::string &directory)
Set output directory (we try to open a file in it to see if the directory exists – if it doesn't we'l...
std::string Total_memory_usage_system_string
String containing system command that obtains total memory usage. Default assignment for linux...
void pause(std::string message)
Pause and display message.
void empty_top_file()
Function to empty file that records continuous output from top in file whose name is specified by Top...
Vector< std::string > colour
Tecplot colours.
void doc_available_flags()
Document available command line flags.
long RefineableQElement< 2 > _build
void stop_continous_top(const std::string &comment)
Stop running top continuously. Note that this is again quite Linux specific and unlikely to work on o...
void line_search(const Vector< double > &x_old, const double half_residual_squared_old, const Vector< double > &gradient, ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &newton_dir, Vector< double > &x, double &half_residual_squared, const double &stpmax)
Line search helper for globally convergent Newton method.
double cumulative_time(const unsigned &i)
Report time accumulated by i-th timer.
void doc_memory_usage(const std::string &prefix_string)
Doc total and local memory usage, prepended by string (which allows identification from where the fun...
void doc_specified_flags()
Document specified command line flags.
void doc_all_flags(std::ostream &outstream)
Document the values of all flags (specified or not).
void(* ResidualFctPt)(const Vector< double > &, const Vector< double > &, Vector< double > &)
Function pointer for residual function: Parameters, unknowns, residuals.
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
bool FlagObsoleteCode
Flag up obsolete parts of the code.
void halt(const unsigned &i)
Halt i-th timer.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57