MFEM  v4.6.0
Finite element discretization library
ex18p.cpp
Go to the documentation of this file.
1 // MFEM Example 18 - Parallel Version
2 //
3 // Compile with: make ex18
4 //
5 // Sample runs:
6 //
7 // mpirun -np 4 ex18p -p 1 -rs 2 -rp 1 -o 1 -s 3
8 // mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 3 -s 4
9 // mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 5 -s 6
10 // mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 1 -s 3
11 // mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 3 -s 3
12 //
13 // Description: This example code solves the compressible Euler system of
14 // equations, a model nonlinear hyperbolic PDE, with a
15 // discontinuous Galerkin (DG) formulation.
16 //
17 // Specifically, it solves for an exact solution of the equations
18 // whereby a vortex is transported by a uniform flow. Since all
19 // boundaries are periodic here, the method's accuracy can be
20 // assessed by measuring the difference between the solution and
21 // the initial condition at a later time when the vortex returns
22 // to its initial location.
23 //
24 // Note that as the order of the spatial discretization increases,
25 // the timestep must become smaller. This example currently uses a
26 // simple estimate derived by Cockburn and Shu for the 1D RKDG
27 // method. An additional factor can be tuned by passing the --cfl
28 // (or -c shorter) flag.
29 //
30 // The example demonstrates user-defined bilinear and nonlinear
31 // form integrators for systems of equations that are defined with
32 // block vectors, and how these are used with an operator for
33 // explicit time integrators. In this case the system also
34 // involves an external approximate Riemann solver for the DG
35 // interface flux. It also demonstrates how to use GLVis for
36 // in-situ visualization of vector grid functions.
37 //
38 // We recommend viewing examples 9, 14 and 17 before viewing this
39 // example.
40 
41 #include "mfem.hpp"
42 #include <fstream>
43 #include <sstream>
44 #include <iostream>
45 
46 // Classes FE_Evolution, RiemannSolver, and FaceIntegrator
47 // shared between the serial and parallel version of the example.
48 #include "ex18.hpp"
49 
50 // Choice for the problem setup. See InitialCondition in ex18.hpp.
51 int problem;
52 
53 // Equation constant parameters.
54 const int num_equation = 4;
55 const double specific_heat_ratio = 1.4;
56 const double gas_constant = 1.0;
57 
58 // Maximum characteristic speed (updated by integrators)
60 
61 int main(int argc, char *argv[])
62 {
63  // 1. Initialize MPI and HYPRE.
64  Mpi::Init(argc, argv);
65  Hypre::Init();
66 
67  // 2. Parse command-line options.
68  problem = 1;
69  const char *mesh_file = "../data/periodic-square.mesh";
70  int ser_ref_levels = 0;
71  int par_ref_levels = 1;
72  int order = 3;
73  int ode_solver_type = 4;
74  double t_final = 2.0;
75  double dt = -0.01;
76  double cfl = 0.3;
77  bool visualization = true;
78  int vis_steps = 50;
79 
80  int precision = 8;
81  cout.precision(precision);
82 
83  OptionsParser args(argc, argv);
84  args.AddOption(&mesh_file, "-m", "--mesh",
85  "Mesh file to use.");
86  args.AddOption(&problem, "-p", "--problem",
87  "Problem setup to use. See options in velocity_function().");
88  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
89  "Number of times to refine the mesh uniformly before parallel"
90  " partitioning, -1 for auto.");
91  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
92  "Number of times to refine the mesh uniformly after parallel"
93  " partitioning.");
94  args.AddOption(&order, "-o", "--order",
95  "Order (degree) of the finite elements.");
96  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
97  "ODE solver: 1 - Forward Euler,\n\t"
98  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
99  args.AddOption(&t_final, "-tf", "--t-final",
100  "Final time; start time is 0.");
101  args.AddOption(&dt, "-dt", "--time-step",
102  "Time step. Positive number skips CFL timestep calculation.");
103  args.AddOption(&cfl, "-c", "--cfl-number",
104  "CFL number for timestep calculation.");
105  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
106  "--no-visualization",
107  "Enable or disable GLVis visualization.");
108  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
109  "Visualize every n-th timestep.");
110 
111  args.Parse();
112  if (!args.Good())
113  {
114  if (Mpi::Root()) { args.PrintUsage(cout); }
115  return 1;
116  }
117  if (Mpi::Root()) { args.PrintOptions(cout); }
118 
119  // 3. Read the mesh from the given mesh file. This example requires a 2D
120  // periodic mesh, such as ../data/periodic-square.mesh.
121  Mesh mesh(mesh_file, 1, 1);
122  const int dim = mesh.Dimension();
123 
124  MFEM_ASSERT(dim == 2, "Need a two-dimensional mesh for the problem definition");
125 
126  // 4. Define the ODE solver used for time integration. Several explicit
127  // Runge-Kutta methods are available.
128  ODESolver *ode_solver = NULL;
129  switch (ode_solver_type)
130  {
131  case 1: ode_solver = new ForwardEulerSolver; break;
132  case 2: ode_solver = new RK2Solver(1.0); break;
133  case 3: ode_solver = new RK3SSPSolver; break;
134  case 4: ode_solver = new RK4Solver; break;
135  case 6: ode_solver = new RK6Solver; break;
136  default:
137  if (Mpi::Root())
138  {
139  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
140  }
141  return 3;
142  }
143 
144  // 5. Refine the mesh in serial to increase the resolution. In this example
145  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
146  // a command-line parameter.
147  for (int lev = 0; lev < ser_ref_levels; lev++)
148  {
149  mesh.UniformRefinement();
150  }
151 
152  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
153  // this mesh further in parallel to increase the resolution. Once the
154  // parallel mesh is defined, the serial mesh can be deleted.
155  ParMesh pmesh(MPI_COMM_WORLD, mesh);
156  mesh.Clear();
157  for (int lev = 0; lev < par_ref_levels; lev++)
158  {
159  pmesh.UniformRefinement();
160  }
161 
162  // 7. Define the discontinuous DG finite element space of the given
163  // polynomial order on the refined mesh.
164  DG_FECollection fec(order, dim);
165  // Finite element space for a scalar (thermodynamic quantity)
166  ParFiniteElementSpace fes(&pmesh, &fec);
167  // Finite element space for a mesh-dim vector quantity (momentum)
168  ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
169  // Finite element space for all variables together (total thermodynamic state)
170  ParFiniteElementSpace vfes(&pmesh, &fec, num_equation, Ordering::byNODES);
171 
172  // This example depends on this ordering of the space.
173  MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
174 
175  HYPRE_BigInt glob_size = vfes.GlobalTrueVSize();
176  if (Mpi::Root())
177  {
178  cout << "Number of unknowns: " << glob_size << endl;
179  }
180 
181  // 8. Define the initial conditions, save the corresponding mesh and grid
182  // functions to a file. This can be opened with GLVis with the -gc option.
183 
184  // The solution u has components {density, x-momentum, y-momentum, energy}.
185  // These are stored contiguously in the BlockVector u_block.
186  Array<int> offsets(num_equation + 1);
187  for (int k = 0; k <= num_equation; k++) { offsets[k] = k * vfes.GetNDofs(); }
188  BlockVector u_block(offsets);
189 
190  // Momentum grid function on dfes for visualization.
191  ParGridFunction mom(&dfes, u_block.GetData() + offsets[1]);
192 
193  // Initialize the state.
194  VectorFunctionCoefficient u0(num_equation, InitialCondition);
195  ParGridFunction sol(&vfes, u_block.GetData());
196  sol.ProjectCoefficient(u0);
197 
198  // Output the initial solution.
199  {
200  ostringstream mesh_name;
201  mesh_name << "vortex-mesh." << setfill('0')
202  << setw(6) << Mpi::WorldRank();
203  ofstream mesh_ofs(mesh_name.str().c_str());
204  mesh_ofs.precision(precision);
205  mesh_ofs << pmesh;
206 
207  for (int k = 0; k < num_equation; k++)
208  {
209  ParGridFunction uk(&fes, u_block.GetBlock(k));
210  ostringstream sol_name;
211  sol_name << "vortex-" << k << "-init."
212  << setfill('0') << setw(6) << Mpi::WorldRank();
213  ofstream sol_ofs(sol_name.str().c_str());
214  sol_ofs.precision(precision);
215  sol_ofs << uk;
216  }
217  }
218 
219  // 9. Set up the nonlinear form corresponding to the DG discretization of the
220  // flux divergence, and assemble the corresponding mass matrix.
221  MixedBilinearForm Aflux(&dfes, &fes);
222  Aflux.AddDomainIntegrator(new TransposeIntegrator(new GradientIntegrator()));
223  Aflux.Assemble();
224 
225  ParNonlinearForm A(&vfes);
226  RiemannSolver rsolver;
227  A.AddInteriorFaceIntegrator(new FaceIntegrator(rsolver, dim));
228 
229  // 10. Define the time-dependent evolution operator describing the ODE
230  // right-hand side, and perform time-integration (looping over the time
231  // iterations, ti, with a time-step dt).
232  FE_Evolution euler(vfes, A, Aflux.SpMat());
233 
234  // Visualize the density
235  socketstream sout;
236  if (visualization)
237  {
238  char vishost[] = "localhost";
239  int visport = 19916;
240 
241  MPI_Barrier(pmesh.GetComm());
242  sout.open(vishost, visport);
243  if (!sout)
244  {
245  if (Mpi::Root())
246  {
247  cout << "Unable to connect to GLVis server at "
248  << vishost << ':' << visport << endl;
249  }
250  visualization = false;
251  if (Mpi::Root())
252  {
253  cout << "GLVis visualization disabled.\n";
254  }
255  }
256  else
257  {
258  sout << "parallel " << Mpi::WorldSize()
259  << " " << Mpi::WorldRank() << "\n";
260  sout.precision(precision);
261  sout << "solution\n" << pmesh << mom;
262  sout << "pause\n";
263  sout << flush;
264  if (Mpi::Root())
265  {
266  cout << "GLVis visualization paused."
267  << " Press space (in the GLVis window) to resume it.\n";
268  }
269  }
270  }
271 
272  // Determine the minimum element size.
273  double hmin;
274  if (cfl > 0)
275  {
276  double my_hmin = pmesh.GetElementSize(0, 1);
277  for (int i = 1; i < pmesh.GetNE(); i++)
278  {
279  my_hmin = min(pmesh.GetElementSize(i, 1), my_hmin);
280  }
281  // Reduce to find the global minimum element size
282  MPI_Allreduce(&my_hmin, &hmin, 1, MPI_DOUBLE, MPI_MIN, pmesh.GetComm());
283  }
284 
285  // Start the timer.
286  tic_toc.Clear();
287  tic_toc.Start();
288 
289  double t = 0.0;
290  euler.SetTime(t);
291  ode_solver->Init(euler);
292 
293  if (cfl > 0)
294  {
295  // Find a safe dt, using a temporary vector. Calling Mult() computes the
296  // maximum char speed at all quadrature points on all faces.
297  max_char_speed = 0.;
298  Vector z(sol.Size());
299  A.Mult(sol, z);
300  // Reduce to find the global maximum wave speed
301  {
302  double all_max_char_speed;
303  MPI_Allreduce(&max_char_speed, &all_max_char_speed,
304  1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
305  max_char_speed = all_max_char_speed;
306  }
307  dt = cfl * hmin / max_char_speed / (2*order+1);
308  }
309 
310  // Integrate in time.
311  bool done = false;
312  for (int ti = 0; !done; )
313  {
314  double dt_real = min(dt, t_final - t);
315 
316  ode_solver->Step(sol, t, dt_real);
317  if (cfl > 0)
318  {
319  // Reduce to find the global maximum wave speed
320  {
321  double all_max_char_speed;
322  MPI_Allreduce(&max_char_speed, &all_max_char_speed,
323  1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
324  max_char_speed = all_max_char_speed;
325  }
326  dt = cfl * hmin / max_char_speed / (2*order+1);
327  }
328  ti++;
329 
330  done = (t >= t_final - 1e-8*dt);
331  if (done || ti % vis_steps == 0)
332  {
333  if (Mpi::Root())
334  {
335  cout << "time step: " << ti << ", time: " << t << endl;
336  }
337  if (visualization)
338  {
339  MPI_Barrier(pmesh.GetComm());
340  sout << "parallel " << Mpi::WorldSize()
341  << " " << Mpi::WorldRank() << "\n";
342  sout << "solution\n" << pmesh << mom << flush;
343  }
344  }
345  }
346 
347  tic_toc.Stop();
348  if (Mpi::Root())
349  {
350  cout << " done, " << tic_toc.RealTime() << "s." << endl;
351  }
352 
353  // 11. Save the final solution. This output can be viewed later using GLVis:
354  // "glvis -np 4 -m vortex-mesh -g vortex-1-final".
355  for (int k = 0; k < num_equation; k++)
356  {
357  ParGridFunction uk(&fes, u_block.GetBlock(k));
358  ostringstream sol_name;
359  sol_name << "vortex-" << k << "-final."
360  << setfill('0') << setw(6) << Mpi::WorldRank();
361  ofstream sol_ofs(sol_name.str().c_str());
362  sol_ofs.precision(precision);
363  sol_ofs << uk;
364  }
365 
366  // 12. Compute the L2 solution error summed for all components.
367  if (t_final == 2.0)
368  {
369  const double error = sol.ComputeLpError(2, u0);
370  if (Mpi::Root())
371  {
372  cout << "Solution error: " << error << endl;
373  }
374  }
375 
376  // Free the used memory.
377  delete ode_solver;
378 
379  return 0;
380 }
int visport
StopWatch tic_toc
Definition: tic_toc.cpp:447
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition: tic_toc.cpp:429
const int num_equation
Definition: ex18p.cpp:54
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:419
char vishost[]
double max_char_speed
Definition: ex18p.cpp:59
const double specific_heat_ratio
Definition: ex18p.cpp:55
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:408
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
Definition: fe_coll.hpp:377
const double gas_constant
Definition: ex18p.cpp:56
HYPRE_Int HYPRE_BigInt
int problem
Definition: ex18p.cpp:51
int dim
Definition: ex24.cpp:53
void InitialCondition(const Vector &x, Vector &y)
Definition: ex18.hpp:430
RefCoord t[3]
int main(int argc, char *argv[])
Definition: ex18p.cpp:61
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:403