MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex20p.cpp
Go to the documentation of this file.
1 // MFEM Example 20 - Parallel Version
2 //
3 // Compile with: make ex20p
4 //
5 // Sample runs: mpirun -np 4 ex20p
6 //
7 // Description: This example demonstrates the use of the variable order,
8 // symplectic ODE integration algorithm. Symplectic integration
9 // algorithms are designed to conserve energy when integrating, in
10 // time, systems of ODEs which are derived from Hamiltonian
11 // systems.
12 //
13 // Hamiltonian systems define the energy of a system as a function
14 // of time (t), a set of generalized coordinates (q), and their
15 // corresponding generalized momenta (p).
16 //
17 // H(q,p,t) = T(p) + V(q,t)
18 //
19 // Hamilton's equations then specify how q and p evolve in time:
20 //
21 // dq/dt = dH/dp
22 // dp/dt = -dH/dq
23 //
24 // To use the symplectic integration classes we need to define an
25 // mfem::Operator P which evaluates the action of dH/dp, and an
26 // mfem::TimeDependentOperator F which computes -dH/dq.
27 //
28 // This example offers five simple 1D Hamiltonians:
29 // 0) Simple Harmonic Oscillator (mass on a spring)
30 // H = ( p^2 / m + q^2 / k ) / 2
31 // 1) Pendulum
32 // H = ( p^2 / m - k ( 1 - cos(q) ) ) / 2
33 // 2) Gaussian Potential Well
34 // H = ( p^2 / m ) / 2 - k exp(-q^2 / 2)
35 // 3) Quartic Potential
36 // H = ( p^2 / m + k ( 1 + q^2 ) q^2 ) / 2
37 // 4) Negative Quartic Potential
38 // H = ( p^2 / m + k ( 1 - q^2 /8 ) q^2 ) / 2
39 //
40 // In all cases these Hamiltonians are shifted by constant values
41 // so that the energy will remain positive. The mean and standard
42 // deviation of the computed energies at each time step are
43 // displayed upon completion. When run in parallel the same
44 // Hamiltonian system is evolved on each processor but starting
45 // from different initial conditions.
46 //
47 // We then use GLVis to visualize the results in a non-standard way
48 // by defining the axes to be q, p, and t rather than x, y, and z.
49 // In this space we build a ribbon-like mesh on each processor with
50 // nodes at (0,0,t) and (q,p,t). When these ribbons are bonded
51 // together on the t-axis they resemble a Rotini pasta. Finally we
52 // plot the energy as a function of time as a scalar field on this
53 // Rotini-like mesh.
54 //
55 // For a more traditional plot of the results, including q, p, and
56 // H from each processor, can be obtained by selecting the "-gp"
57 // option. This creates a collection of data files and an input
58 // deck for the GnuPlot application (not included with MFEM). To
59 // visualize these results on most linux systems type the command
60 // "gnuplot gnuplot_ex20p.inp". The data files, named
61 // "ex20p_?????.dat", should be simple enough to display with other
62 // plotting programs as well.
63 
64 #include "mfem.hpp"
65 #include <fstream>
66 #include <iostream>
67 
68 using namespace std;
69 using namespace mfem;
70 
71 // Constants used in the Hamiltonian
72 static int prob_ = 0;
73 static double m_ = 1.0;
74 static double k_ = 1.0;
75 
76 // Hamiltonian functional, see below for implementation
77 double hamiltonian(double q, double p, double t);
78 
79 class GradT : public Operator
80 {
81 public:
82  GradT() : Operator(1) {}
83  void Mult(const Vector &x, Vector &y) const { y.Set(1.0/m_, x); }
84 };
85 
86 class NegGradV : public TimeDependentOperator
87 {
88 public:
89  NegGradV() : TimeDependentOperator(1) {}
90  void Mult(const Vector &x, Vector &y) const;
91 };
92 
93 int main(int argc, char *argv[])
94 {
95  // 1. Initialize MPI.
96  int num_procs, myid;
97  MPI_Comm comm = MPI_COMM_WORLD;
98  MPI_Init(&argc, &argv);
99  MPI_Comm_size(comm, &num_procs);
100  MPI_Comm_rank(comm, &myid);
101 
102  // 2. Parse command-line options.
103  int order = 1;
104  int nsteps = 100;
105  double dt = 0.1;
106  bool visualization = true;
107  bool gnuplot = false;
108 
109  OptionsParser args(argc, argv);
110  args.AddOption(&order, "-o", "--order",
111  "Time integration order.");
112  args.AddOption(&prob_, "-p", "--problem-type",
113  "Problem Type:\n"
114  "\t 0 - Simple Harmonic Oscillator\n"
115  "\t 1 - Pendulum\n"
116  "\t 2 - Gaussian Potential Well\n"
117  "\t 3 - Quartic Potential\n"
118  "\t 4 - Negative Quartic Potential");
119  args.AddOption(&nsteps, "-n", "--number-of-steps",
120  "Number of time steps.");
121  args.AddOption(&dt, "-dt", "--time-step",
122  "Time step size.");
123  args.AddOption(&m_, "-m", "--mass",
124  "Mass.");
125  args.AddOption(&k_, "-k", "--spring-const",
126  "Spring constant.");
127  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
128  "--no-visualization",
129  "Enable or disable GLVis visualization.");
130  args.AddOption(&gnuplot, "-gp", "--gnuplot", "-no-gp", "--no-gnuplot",
131  "Enable or disable GnuPlot visualization.");
132  args.Parse();
133  if (!args.Good())
134  {
135  if (myid == 0)
136  {
137  args.PrintUsage(cout);
138  }
139  MPI_Finalize();
140  return 1;
141  }
142  if (myid == 0)
143  {
144  args.PrintOptions(cout);
145  }
146 
147  // 3. Create and Initialize the Symplectic Integration Solver
148  SIAVSolver siaSolver(order);
149  GradT P;
150  NegGradV F;
151  siaSolver.Init(P,F);
152 
153  // 4. Set the initial conditions
154  double t = 0.0;
155  Vector q(1), p(1);
156  Vector e(nsteps+1);
157  q(0) = sin(2.0*M_PI*(double)myid/num_procs);
158  p(0) = cos(2.0*M_PI*(double)myid/num_procs);
159 
160  // 5. Prepare GnuPlot output file if needed
161  ostringstream oss;
162  ofstream ofs;
163  if (gnuplot)
164  {
165  oss << "ex20p_" << setfill('0') << setw(5) << myid << ".dat";
166  ofs.open(oss.str().c_str());
167  ofs << t << "\t" << q(0) << "\t" << p(0) << endl;
168  }
169 
170  // 6. Create a Mesh for visualization in phase space
171  int nverts = (visualization) ? (num_procs+1)*(nsteps+1) : 0;
172  int nelems = (visualization) ? (nsteps * num_procs) : 0;
173  Mesh mesh(2, nverts, nelems, 0, 3);
174 
175  int *part = (visualization) ? (new int[nelems]) : NULL;
176  int v[4];
177  Vector x0(3); x0 = 0.0;
178  Vector x1(3); x1 = 0.0;
179 
180  // 7. Perform time-stepping
181  double e_mean = 0.0;
182 
183  for (int i = 0; i < nsteps; i++)
184  {
185  // 7a. Record initial state
186  if (i == 0)
187  {
188  e[0] = hamiltonian(q(0),p(0),t);
189  e_mean += e[0];
190 
191  if (visualization)
192  {
193  mesh.AddVertex(x0);
194  for (int j = 0; j < num_procs; j++)
195  {
196  x1[0] = q(0);
197  x1[1] = p(0);
198  x1[2] = 0.0;
199  mesh.AddVertex(x1);
200  }
201  }
202  }
203 
204  // 7b. Advance the state of the system
205  siaSolver.Step(q,p,t,dt);
206  e[i+1] = hamiltonian(q(0),p(0),t);
207  e_mean += e[i+1];
208 
209  // 7c. Record the state of the system
210  if (gnuplot)
211  {
212  ofs << t << "\t" << q(0) << "\t" << p(0) << "\t" << e[i+1] << endl;
213  }
214 
215  // 7d. Add results to GLVis visualization
216  if (visualization)
217  {
218  x0[2] = t;
219  mesh.AddVertex(x0);
220  for (int j = 0; j < num_procs; j++)
221  {
222  x1[0] = q(0);
223  x1[1] = p(0);
224  x1[2] = t;
225  mesh.AddVertex(x1);
226  v[0] = (num_procs + 1) * i;
227  v[1] = (num_procs + 1) * (i + 1);
228  v[2] = (num_procs + 1) * (i + 1) + j + 1;
229  v[3] = (num_procs + 1) * i + j + 1;
230  mesh.AddQuad(v);
231  part[num_procs * i + j] = j;
232  }
233  }
234  }
235 
236  // 8. Compute and display mean and standard deviation of the energy
237  e_mean /= (nsteps + 1);
238  double e_var = 0.0;
239  for (int i = 0; i <= nsteps; i++)
240  {
241  e_var += pow(e[i] - e_mean, 2);
242  }
243  e_var /= (nsteps + 1);
244  double e_sd = sqrt(e_var);
245 
246  if (myid == 0)
247  {
248  cout << endl << "Mean and standard deviation of the energy" << endl;
249  }
250  for (int i = 0; i < num_procs; i++)
251  {
252  if (myid == i)
253  {
254  cout << myid << ": " << e_mean << "\t" << e_sd << endl;
255  }
256  MPI_Barrier(comm);
257  }
258 
259  // 9. Finalize the GnuPlot output
260  if (gnuplot)
261  {
262  ofs.close();
263  if (myid == 0)
264  {
265  ofs.open("gnuplot_ex20p.inp");
266  for (int i = 0; i < num_procs; i++)
267  {
268  ostringstream ossi;
269  ossi << "ex20p_" << setfill('0') << setw(5) << i << ".dat";
270  if (i == 0)
271  {
272  ofs << "plot";
273  }
274  ofs << " '" << ossi.str() << "' using 1:2 w l t 'q" << i << "',"
275  << " '" << ossi.str() << "' using 1:3 w l t 'p" << i << "',"
276  << " '" << ossi.str() << "' using 1:4 w l t 'H" << i << "'";
277  if (i < num_procs-1)
278  {
279  ofs << ",";
280  }
281  else
282  {
283  ofs << ";" << endl;
284  }
285  }
286  ofs.close();
287  }
288  }
289 
290  // 10. Finalize the GLVis output
291  if (visualization)
292  {
293  mesh.FinalizeQuadMesh(1);
294  ParMesh pmesh(comm, mesh, part);
295  delete [] part;
296 
297  H1_FECollection fec(order = 1, 2);
298  ParFiniteElementSpace fespace(&pmesh, &fec);
299  ParGridFunction energy(&fespace);
300  energy = 0.0;
301  for (int i = 0; i <= nsteps; i++)
302  {
303  energy[2*i+0] = e[i];
304  energy[2*i+1] = e[i];
305  }
306 
307  char vishost[] = "localhost";
308  int visport = 19916;
309  socketstream sock(vishost, visport);
310  sock.precision(8);
311  sock << "parallel " << num_procs << " " << myid << "\n"
312  << "solution\n" << pmesh << energy
313  << "window_title 'Energy in Phase Space'\n"
314  << "keys\n maac\n" << "axis_labels 'q' 'p' 't'\n"<< flush;
315  }
316 
317  MPI_Finalize();
318 }
319 
320 double hamiltonian(double q, double p, double t)
321 {
322  double h = 1.0 - 0.5 / m_ + 0.5 * p * p / m_;
323  switch (prob_)
324  {
325  case 1:
326  h += k_ * (1.0 - cos(q));
327  break;
328  case 2:
329  h += k_ * (1.0 - exp(-0.5 * q * q));
330  break;
331  case 3:
332  h += 0.5 * k_ * (1.0 + q * q) * q * q;
333  break;
334  case 4:
335  h += 0.5 * k_ * (1.0 - 0.125 * q * q) * q * q;
336  break;
337  default:
338  h += 0.5 * k_ * q * q;
339  break;
340  }
341  return h;
342 }
343 
344 void NegGradV::Mult(const Vector &x, Vector &y) const
345 {
346  switch (prob_)
347  {
348  case 1:
349  y(0) = - k_* sin(x(0));
350  break;
351  case 2:
352  y(0) = - k_ * x(0) * exp(-0.5 * x(0) * x(0));
353  break;
354  case 3:
355  y(0) = - k_ * (1.0 + 2.0 * x(0) * x(0)) * x(0);
356  break;
357  case 4:
358  y(0) = - k_ * (1.0 - 0.25 * x(0) * x(0)) * x(0);
359  break;
360  default:
361  y(0) = - k_ * x(0);
362  break;
363  };
364 }
double hamiltonian(double q, double p, double t)
Definition: ex20.cpp:254
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:713
Base abstract class for first order time dependent operators.
Definition: operator.hpp:259
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int main(int argc, char *argv[])
Definition: ex1.cpp:62
void AddVertex(const double *)
Definition: mesh.cpp:1155
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Definition: optparser.hpp:76
void AddQuad(const int *vi, int attr=1)
Definition: mesh.cpp:1181
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:205
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:792
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1357
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
Vector data type.
Definition: vector.hpp:48
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Abstract operator.
Definition: operator.hpp:24
Class for parallel meshes.
Definition: pmesh.hpp:32
bool Good() const
Definition: optparser.hpp:125