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