ElmerFoamFSI  2.0
ElmerFoamFSI is fluid-solid interaction simulation application built up from OpenFOAM CFD and Elmer CSM coupled through the IMPACT multiphysics software integration infrastructure.
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Groups Pages
src/ExampleParallelProgram.C
Go to the documentation of this file.
1 
9 #include "ExampleProgram.H"
10 #include "ExampleHeader.H"
11 #include <cmath>
12 #include <iomanip>
13 
14 namespace ElmerFoamFSI {
15 
16  namespace ExampleProgram {
17 
18  int ParallelProgram::Run()
19  {
20  // FunctionEntry("NAME"): Updates the user-defined stack and
21  // the program profiles with timing information. The placement
22  // of FunctionEntry and FunctionExit calls is at the developer's
23  // discretion.
24  FunctionEntry("Run");
25 
26  int myid = _communicator.Rank();
27  int nproc = _communicator.Size();
28  bool use_file = !output_name.empty();
29 
30  // Put out a quick blurb about number of procs just
31  // to give the user confidence we are actually running
32  // in parallel. Note that when using the "StdOut" method
33  // that it automatically does output only on rank 0. If
34  // the developer wants asyncrhonous output, then they
35  // have to revert to using standard streams.
36  std::ostringstream RepStr;
37  if(verblevel > 1)
38  RepStr << "Running on " << nproc << " processors." << std::endl;
39  StdOut(RepStr.str());
40  _communicator.Barrier();
41  if(verblevel > 1)
42  StdOut("All procesors ready.\n");
43 
44  // Open the specified output file for writing on rank 0
45  if(use_file && !myid){
46  Ouf.open(output_name.c_str(),std::ios::app);
47  if(!Ouf){
48  // If the output file failed to open, notify
49  // to error stream and return non-zero
50  std::ostringstream Ostr;
51  Ostr << "Error: Unable to open output file, " << output_name << ".";
52  ErrOut(Ostr.str());
53  // In parallel, we don't return right away since
54  // this part of the code is only done on proc 0.
55  // Instead, the error value is set in the communicator
56  // which will indicate to all processors that there
57  // has been some error.
58  _communicator.SetExit(1);
59  }
60  }
61  // Check to see if an error condition was set
62  // in the file open block above. If so, then
63  // return with an error code.
64  if(_communicator.Check()){
65  // don't forget to tell the profiler/stacker the
66  // function is exiting.
67  FunctionExit("Run");
68  return(1);
69  }
70 
71  // Correct value of PI to several places
72  double PIVAL = 3.141592653589793238462643383;
73 
74  // All processors should already have this as it
75  // came from the command line - but just in case,
76  // we broadcast it here to make sure.
77  FunctionEntry("Broadcast");
78  _communicator.BroadCast(ndiv,0);
79  FunctionExit("Broadcast");
80 
81  // This block partitions the
82  // domain [0,1] with ndiv
83  // intervals among nproc processors
84  int nper = ndiv/nproc;
85  int nvol = nper*nproc;
86  int leftover = ndiv - nvol;
87  int myn = nper;
88  if(myid < leftover)
89  myn++;
90  int mystart = 0;
91  for(int i = 0; i < myid;i++){
92  mystart += nper;
93  if(i < leftover)
94  mystart++;
95  }
96  // If there are less divisions than processors, then
97  // just forget about domain decomposition and do everything
98  // on processor 0.
99  if(nper == 0){
100  if(myid == 0) myn = ndiv;
101  else myn = 0;
102  }
103 
104  // Use the results from domain decomposition to
105  // set the domain [a,b] for this processor and get the
106  // stepsize.
107  double h = 1.0/(static_cast<double>(ndiv));
108  double a = h*mystart;
109  double b = a + h*myn;
110 
111  // Integrate f on this processor's subdomain using trapezoid quadrature
112  FunctionEntry("TrapezoidMethod");
113  double my_tpi = 0.0;
114  if(myn > 0){ // if there are points in this processor's domain
115  try {
116  FunctionEntry("TrapezoidQuadrature");
117  my_tpi = ElmerFoamFSI::TrapezoidQuadrature(f,a,b,myn);
118  FunctionExit("TrapezoidQuadrature");
119  } catch (...){
120  StdOut("Numerical limits of ElmerFoamFSI::TrapezoidQuadrature exceeded.");
121  my_tpi = 0.0;
122  }
123  }
124  double tpi = 0.0;
125  // Sum up the contributions from each processor and store the results on
126  // processor 0 in "tpi".
127  FunctionEntry("Reduce"); // time the communication
128  _communicator.Reduce(my_tpi, tpi,IRAD::Comm::DTDOUBLE, IRAD::Comm::SUMOP, 0);
129  FunctionExit("Reduce"); // end of the communication
130  FunctionExit("TrapezoidMethod");
131 
132  // Integrate f on this processor's subdomain using midpoint quadrature
133  FunctionEntry("MidPointMethod");
134  double my_mppi = 0.0;
135  if(myn > 0) { // if there are points in this processor's domain
136  try{
137  FunctionEntry("MidPointQuadrature");
138  my_mppi = ElmerFoamFSI::MidPointQuadrature(f,a,b,myn);
139  FunctionExit("MidPointQuadrature");
140  } catch (...){
141  StdOut("Numerical limits of ElmerFoamFSI::MidPointQuadrature exceeded.");
142  my_mppi = 0.0;
143  }
144  }
145  double mppi = 0.0;
146  // Sum up the contributions from each processor and store the results on
147  // processor 0 in "mppi".
148  FunctionEntry("Reduce"); // time the communication
149  _communicator.Reduce(my_mppi, mppi,IRAD::Comm::DTDOUBLE, IRAD::Comm::SUMOP, 0);
150  FunctionExit("Reduce"); // end of the communication
151  FunctionExit("MidPointMethod");
152 
153  // Report the results to stdout or to file if one was specified.
154  FunctionEntry("IO");
155  if (!myid) {
156  if(use_file){
157  Ouf << ndiv << " " << std::setprecision(16) << tpi << " "
158  << std::setprecision(4) << std::fabs(tpi-PIVAL) << " "
159  << std::setprecision(16) << mppi << " "
160  << std::setprecision(4) << std::fabs(mppi-PIVAL) << " "
161  << std::endl;
162  Ouf.close();
163  }
164  else if(verblevel) {
165  std::ostringstream Ostr;
166  Ostr << "With " << ndiv << " divisions, PI was calculated:" << std::endl
167  << "MidPointQuadrature: "
168  << std::setprecision(16) << mppi << "\t\t"
169  << std::setprecision(4) << std::fabs(mppi - PIVAL) << std::endl
170  << "TrapezoidQuadrature: "
171  << std::setprecision(16) << tpi << "\t\t"
172  << std::setprecision(4) << std::fabs(tpi - PIVAL) << std::endl;
173  StdOut(Ostr.str());
174  }
175  }
176  FunctionExit("IO");
177  //
178  // ---------- Program End ----------------
179 
180  // Update the stacker/profiler that we are exiting
181  // this function.
182  FunctionExit("Run");
183  // return 0 for success
184  return(0);
185  };
186  };
187 };
188 
double MidPointQuadrature(double(*f)(double), double x0, double xn, int n)
Integrates f with composite midpoint rule.
double TrapezoidQuadrature(double(*f)(double), double x0, double xn, int n)
Integrates f with composite trapezoid rule.