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