MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
72using namespace std;
73using namespace mfem;
74
75// Constants used in the Hamiltonian
76static int prob_ = 0;
77static real_t m_ = 1.0;
78static real_t k_ = 1.0;
79
80// Hamiltonian functional, see below for implementation
82
84{
85public:
87 void Mult(const Vector &x, Vector &y) const { y.Set(1.0/m_, x); }
88};
89
91{
92public:
94 void Mult(const Vector &x, Vector &y) const;
95};
96
97int 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 real_t dt = 0.1;
110 bool visualization = true;
111 bool gnuplot = false;
112
113 OptionsParser args(argc, argv);
115 "Time integration order.");
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");
124 "Number of time steps.");
126 "Time step size.");
128 "Mass.");
130 "Spring constant.");
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);
154 siaSolver.Init(P,F);
155
156 // 4. Set the initial conditions
157 real_t t = 0.0;
158 Vector q(1), p(1);
159 Vector e(nsteps+1);
160 q(0) = sin(2.0*M_PI*(real_t)myid/num_procs);
161 p(0) = cos(2.0*M_PI*(real_t)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 real_t 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 {
199 x1[0] = q(0);
200 x1[1] = p(0);
201 x1[2] = 0.0;
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 {
225 x1[0] = q(0);
226 x1[1] = p(0);
227 x1[2] = t;
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;
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 real_t 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 real_t e_sd = sqrt(e_var);
248
249 real_t e_loc_stats[2];
250 real_t *e_stats = (myid == 0) ? new real_t[2 * num_procs] : (real_t*)NULL;
251
252 e_loc_stats[0] = e_mean;
253 e_loc_stats[1] = e_sd;
254 MPI_Gather(e_loc_stats, 2, MPITypeMap<real_t>::mpi_type, e_stats, 2,
256
257 if (myid == 0)
258 {
259 cout << endl << "Mean and standard deviation of the energy "
260 << "for different initial conditions" << endl;
261 for (int i = 0; i < num_procs; i++)
262 {
263 cout << i << ": " << e_stats[2 * i + 0]
264 << "\t" << e_stats[2 * i + 1] << endl;
265 }
266 delete [] e_stats;
267 }
268
269 // 9. Finalize the GnuPlot output
270 if (gnuplot)
271 {
272 ofs.close();
273 if (myid == 0)
274 {
275 ofs.open("gnuplot_ex20p.inp");
276 for (int i = 0; i < num_procs; i++)
277 {
278 ostringstream ossi;
279 ossi << "ex20p_" << setfill('0') << setw(5) << i << ".dat";
280 if (i == 0)
281 {
282 ofs << "plot";
283 }
284 ofs << " '" << ossi.str() << "' using 1:2 w l t 'q" << i << "',"
285 << " '" << ossi.str() << "' using 1:3 w l t 'p" << i << "',"
286 << " '" << ossi.str() << "' using 1:4 w l t 'H" << i << "'";
287 if (i < num_procs-1)
288 {
289 ofs << ",";
290 }
291 else
292 {
293 ofs << ";" << endl;
294 }
295 }
296 ofs.close();
297 }
298 }
299
300 // 10. Finalize the GLVis output
301 if (visualization)
302 {
304 ParMesh pmesh(comm, mesh, part);
305 delete [] part;
306
307 H1_FECollection fec(order = 1, 2);
308 ParFiniteElementSpace fespace(&pmesh, &fec);
309 ParGridFunction energy(&fespace);
310 energy = 0.0;
311 for (int i = 0; i <= nsteps; i++)
312 {
313 energy[2*i+0] = e[i];
314 energy[2*i+1] = e[i];
315 }
316
317 char vishost[] = "localhost";
318 int visport = 19916;
320 sock.precision(8);
321 sock << "parallel " << num_procs << " " << myid << "\n"
322 << "solution\n" << pmesh << energy
323 << "window_title 'Energy in Phase Space'\n"
324 << "keys\n maac\n" << "axis_labels 'q' 'p' 't'\n"<< flush;
325 }
326}
327
329{
330 real_t h = 1.0 - 0.5 / m_ + 0.5 * p * p / m_;
331 switch (prob_)
332 {
333 case 1:
334 h += k_ * (1.0 - cos(q));
335 break;
336 case 2:
337 h += k_ * (1.0 - exp(-0.5 * q * q));
338 break;
339 case 3:
340 h += 0.5 * k_ * (1.0 + q * q) * q * q;
341 break;
342 case 4:
343 h += 0.5 * k_ * (1.0 - 0.125 * q * q) * q * q;
344 break;
345 default:
346 h += 0.5 * k_ * q * q;
347 break;
348 }
349 return h;
350}
351
352void NegGradV::Mult(const Vector &x, Vector &y) const
353{
354 switch (prob_)
355 {
356 case 1:
357 y(0) = - k_* sin(x(0));
358 break;
359 case 2:
360 y(0) = - k_ * x(0) * exp(-0.5 * x(0) * x(0));
361 break;
362 case 3:
363 y(0) = - k_ * (1.0 + 2.0 * x(0) * x(0)) * x(0);
364 break;
365 case 4:
366 y(0) = - k_ * (1.0 - 0.25 * x(0) * x(0)) * x(0);
367 break;
368 default:
369 y(0) = - k_ * x(0);
370 break;
371 };
372}
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Mesh data type.
Definition mesh.hpp:56
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition mesh.cpp:2177
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Abstract operator.
Definition operator.hpp:25
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
Class for parallel meshes.
Definition pmesh.hpp:34
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition ode.cpp:916
Variable order Symplectic Integration Algorithm (orders 1-4)
Definition ode.hpp:612
void Step(Vector &q, Vector &p, real_t &t, real_t &dt) override
Definition ode.cpp:995
Base abstract class for first order time dependent operators.
Definition operator.hpp:317
Vector data type.
Definition vector.hpp:80
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:262
real_t hamiltonian(real_t q, real_t p, real_t t)
Definition ex20p.cpp:328
int main()
const int visport
float real_t
Definition config.hpp:43
const char vishost[]
real_t p(const Vector &x, real_t t)
RefCoord t[3]
Helper struct to convert a C++ type to an MPI type.