MFEM  v4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex20.cpp
Go to the documentation of this file.
1 // MFEM Example 20
2 //
3 // Compile with: make ex20
4 //
5 // Sample runs: ex20
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.
44 //
45 // We then use GLVis to visualize the results in a non-standard way
46 // by defining the axes to be q, p, and t rather than x, y, and z.
47 // In this space we build a ribbon-like mesh with nodes at (0,0,t)
48 // and (q,p,t). Finally we plot the energy as a function of time
49 // as a scalar field on this ribbon-like mesh.
50 //
51 // For a more traditional plot of the results, including q, p, and
52 // H, can be obtained by selecting the "-gp" option. This creates
53 // a data file and input deck for the GnuPlot application (not
54 // included with MFEM). To visualize these results on most Linux
55 // systems type the command "gnuplot gnuplot_ex20.inp". The data
56 // file, named "ex20.dat", should be simple enough to display with
57 // other plotting programs as well.
58 
59 #include "mfem.hpp"
60 #include <fstream>
61 #include <iostream>
62 
63 using namespace std;
64 using namespace mfem;
65 
66 // Constants used in the Hamiltonian
67 static int prob_ = 0;
68 static double m_ = 1.0;
69 static double k_ = 1.0;
70 
71 // Hamiltonian functional, see below for implementation
72 double hamiltonian(double q, double p, double t);
73 
74 class GradT : public Operator
75 {
76 public:
77  GradT() : Operator(1) {}
78  void Mult(const Vector &x, Vector &y) const { y.Set(1.0/m_, x); }
79 };
80 
81 class NegGradV : public TimeDependentOperator
82 {
83 public:
84  NegGradV() : TimeDependentOperator(1) {}
85  void Mult(const Vector &x, Vector &y) const;
86 };
87 
88 int main(int argc, char *argv[])
89 {
90  // 1. Parse command-line options.
91  int order = 1;
92  int nsteps = 100;
93  double dt = 0.1;
94  bool visualization = true;
95  bool gnuplot = false;
96 
97  OptionsParser args(argc, argv);
98  args.AddOption(&order, "-o", "--order",
99  "Time integration order.");
100  args.AddOption(&prob_, "-p", "--problem-type",
101  "Problem Type:\n"
102  "\t 0 - Simple Harmonic Oscillator\n"
103  "\t 1 - Pendulum\n"
104  "\t 2 - Gaussian Potential Well\n"
105  "\t 3 - Quartic Potential\n"
106  "\t 4 - Negative Quartic Potential");
107  args.AddOption(&nsteps, "-n", "--number-of-steps",
108  "Number of time steps.");
109  args.AddOption(&dt, "-dt", "--time-step",
110  "Time step size.");
111  args.AddOption(&m_, "-m", "--mass",
112  "Mass.");
113  args.AddOption(&k_, "-k", "--spring-const",
114  "Spring constant.");
115  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
116  "--no-visualization",
117  "Enable or disable GLVis visualization.");
118  args.AddOption(&gnuplot, "-gp", "--gnuplot", "-no-gp", "--no-gnuplot",
119  "Enable or disable GnuPlot visualization.");
120  args.Parse();
121  if (!args.Good())
122  {
123  args.PrintUsage(cout);
124  return 1;
125  }
126  args.PrintOptions(cout);
127 
128  // 2. Create and Initialize the Symplectic Integration Solver
129  SIAVSolver siaSolver(order);
130  GradT P;
131  NegGradV F;
132  siaSolver.Init(P,F);
133 
134  // 3. Set the initial conditions
135  double t = 0.0;
136  Vector q(1), p(1);
137  Vector e(nsteps+1);
138  q(0) = 0.0;
139  p(0) = 1.0;
140 
141  // 4. Prepare GnuPlot output file if needed
142  ofstream ofs;
143  if (gnuplot)
144  {
145  ofs.open("ex20.dat");
146  ofs << t << "\t" << q(0) << "\t" << p(0) << endl;
147  }
148 
149  // 5. Create a Mesh for visualization in phase space
150  int nverts = (visualization) ? 2*(nsteps+1) : 0;
151  int nelems = (visualization) ? nsteps : 0;
152  Mesh mesh(2, nverts, nelems, 0, 3);
153 
154  int v[4];
155  Vector x0(3); x0 = 0.0;
156  Vector x1(3); x1 = 0.0;
157 
158  // 6. Perform time-stepping
159  double e_mean = 0.0;
160 
161  for (int i = 0; i < nsteps; i++)
162  {
163  // 6a. Record initial state
164  if (i == 0)
165  {
166  e[0] = hamiltonian(q(0),p(0),t);
167  e_mean += e[0];
168 
169  if (visualization)
170  {
171  x1[0] = q(0);
172  x1[1] = p(0);
173  x1[2] = 0.0;
174  mesh.AddVertex(x0);
175  mesh.AddVertex(x1);
176  }
177  }
178 
179  // 6b. Advance the state of the system
180  siaSolver.Step(q,p,t,dt);
181  e[i+1] = hamiltonian(q(0),p(0),t);
182  e_mean += e[i+1];
183 
184  // 6c. Record the state of the system
185  if (gnuplot)
186  {
187  ofs << t << "\t" << q(0) << "\t" << p(0) << "\t" << e[i+1] << endl;
188  }
189 
190  // 6d. Add results to GLVis visualization
191  if (visualization)
192  {
193  x0[2] = t;
194  x1[0] = q(0);
195  x1[1] = p(0);
196  x1[2] = t;
197  mesh.AddVertex(x0);
198  mesh.AddVertex(x1);
199  v[0] = 2*i;
200  v[1] = 2*(i+1);
201  v[2] = 2*(i+1)+1;
202  v[3] = 2*i+1;
203  mesh.AddQuad(v);
204  }
205  }
206 
207  // 7. Compute and display mean and standard deviation of the energy
208  e_mean /= (nsteps + 1);
209  double e_var = 0.0;
210  for (int i=0; i<=nsteps; i++)
211  {
212  e_var += pow(e[i] - e_mean, 2);
213  }
214  e_var /= (nsteps + 1);
215  double e_sd = sqrt(e_var);
216  cout << endl << "Mean and standard deviation of the energy" << endl;
217  cout << e_mean << "\t" << e_sd << endl;
218 
219  // 8. Finalize the GnuPlot output
220  if (gnuplot)
221  {
222  ofs.close();
223 
224  ofs.open("gnuplot_ex20.inp");
225  ofs << "plot 'ex20.dat' using 1:2 w l t 'q', "
226  << "'ex20.dat' using 1:3 w l t 'p', "
227  << "'ex20.dat' using 1:4 w l t 'H'" << endl;
228  ofs.close();
229  }
230 
231  // 9. Finalize the GLVis output
232  if (visualization)
233  {
234  H1_FECollection fec(order = 1, 2);
235  FiniteElementSpace fespace(&mesh, &fec);
236  GridFunction energy(&fespace);
237  energy = 0.0;
238  for (int i = 0; i <= nsteps; i++)
239  {
240  energy[2*i+0] = e[i];
241  energy[2*i+1] = e[i];
242  }
243 
244  char vishost[] = "localhost";
245  int visport = 19916;
246  socketstream sock(vishost, visport);
247  sock.precision(8);
248  sock << "solution\n" << mesh << energy
249  << "window_title 'Energy in Phase Space'\n"
250  << "keys\n maac\n" << "axis_labels 'q' 'p' 't'\n"<< flush;
251  }
252 }
253 
254 double hamiltonian(double q, double p, double t)
255 {
256  double h = 1.0 - 0.5 / m_ + 0.5 * p * p / m_;
257  switch (prob_)
258  {
259  case 1:
260  h += k_ * (1.0 - cos(q));
261  break;
262  case 2:
263  h += k_ * (1.0 - exp(-0.5 * q * q));
264  break;
265  case 3:
266  h += 0.5 * k_ * (1.0 + q * q) * q * q;
267  break;
268  case 4:
269  h += 0.5 * k_ * (1.0 - 0.125 * q * q) * q * q;
270  break;
271  default:
272  h += 0.5 * k_ * q * q;
273  break;
274  }
275  return h;
276 }
277 
278 void NegGradV::Mult(const Vector &x, Vector &y) const
279 {
280  switch (prob_)
281  {
282  case 1:
283  y(0) = - k_* sin(x(0));
284  break;
285  case 2:
286  y(0) = - k_ * x(0) * exp(-0.5 * x(0) * x(0));
287  break;
288  case 3:
289  y(0) = - k_ * (1.0 + 2.0 * x(0) * x(0)) * x(0);
290  break;
291  case 4:
292  y(0) = - k_ * (1.0 - 0.25 * x(0) * x(0)) * x(0);
293  break;
294  default:
295  y(0) = - k_ * x(0);
296  break;
297  };
298 }
double hamiltonian(double q, double p, double t)
Definition: ex20.cpp:254
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:571
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
Base abstract class for time dependent operators.
Definition: operator.hpp:162
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int main(int argc, char *argv[])
Definition: ex1.cpp:58
void AddVertex(const double *)
Definition: mesh.cpp:1117
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:85
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:1138
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:650
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
Abstract operator.
Definition: operator.hpp:21
bool Good() const
Definition: optparser.hpp:122