MFEM  v3.3.2
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 0 -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 0 -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, DomainIntegrator 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.
64  MPI_Session mpi(argc, argv);
65 
66  // 2. Parse command-line options.
67  problem = 1;
68  const char *mesh_file = "../data/periodic-square.mesh";
69  int ser_ref_levels = 0;
70  int par_ref_levels = 1;
71  int order = 3;
72  int ode_solver_type = 4;
73  double t_final = 2.0;
74  double dt = -0.01;
75  double cfl = 0.3;
76  bool visualization = true;
77  int vis_steps = 50;
78 
79  int precision = 8;
80  cout.precision(precision);
81 
82  OptionsParser args(argc, argv);
83  args.AddOption(&mesh_file, "-m", "--mesh",
84  "Mesh file to use.");
85  args.AddOption(&problem, "-p", "--problem",
86  "Problem setup to use. See options in velocity_function().");
87  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
88  "Number of times to refine the mesh uniformly before parallel"
89  " partitioning, -1 for auto.");
90  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
91  "Number of times to refine the mesh uniformly after parallel"
92  " partitioning.");
93  args.AddOption(&order, "-o", "--order",
94  "Order (degree) of the finite elements.");
95  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
96  "ODE solver: 1 - Forward Euler,\n\t"
97  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
98  args.AddOption(&t_final, "-tf", "--t-final",
99  "Final time; start time is 0.");
100  args.AddOption(&dt, "-dt", "--time-step",
101  "Time step. Positive number skips CFL timestep calculation.");
102  args.AddOption(&cfl, "-c", "--cfl-number",
103  "CFL number for timestep calculation.");
104  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
105  "--no-visualization",
106  "Enable or disable GLVis visualization.");
107  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
108  "Visualize every n-th timestep.");
109 
110  args.Parse();
111  if (!args.Good())
112  {
113  if (mpi.Root()) { args.PrintUsage(cout); }
114  return 1;
115  }
116  if (mpi.Root()) { args.PrintOptions(cout); }
117 
118  // 3. Read the mesh from the given mesh file. This example requires a 2D
119  // periodic mesh, such as ../data/periodic-square.mesh.
120  Mesh mesh(mesh_file, 1, 1);
121  const int dim = mesh.Dimension();
122 
123  MFEM_ASSERT(dim == 2, "Need a two-dimensional mesh for the problem definition");
124 
125  // 4. Define the ODE solver used for time integration. Several explicit
126  // Runge-Kutta methods are available.
127  ODESolver *ode_solver = NULL;
128  switch (ode_solver_type)
129  {
130  case 1: ode_solver = new ForwardEulerSolver; break;
131  case 2: ode_solver = new RK2Solver(1.0); break;
132  case 3: ode_solver = new RK3SSPSolver; break;
133  case 4: ode_solver = new RK4Solver; break;
134  case 6: ode_solver = new RK6Solver; break;
135  default:
136  if (mpi.Root())
137  {
138  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
139  }
140  return 3;
141  }
142 
143  // 5. Refine the mesh in serial to increase the resolution. In this example
144  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
145  // a command-line parameter.
146  for (int lev = 0; lev < ser_ref_levels; lev++)
147  {
148  mesh.UniformRefinement();
149  }
150 
151  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
152  // this mesh further in parallel to increase the resolution. Once the
153  // parallel mesh is defined, the serial mesh can be deleted.
154  ParMesh pmesh(MPI_COMM_WORLD, mesh);
155  mesh.Clear();
156  for (int lev = 0; lev < par_ref_levels; lev++)
157  {
158  pmesh.UniformRefinement();
159  }
160 
161  // 7. Define the discontinuous DG finite element space of the given
162  // polynomial order on the refined mesh.
163  DG_FECollection fec(order, dim);
164  // Finite element space for a scalar (thermodynamic quantity)
165  ParFiniteElementSpace fes(&pmesh, &fec);
166  // Finite element space for a mesh-dim vector quantity (momentum)
167  ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
168  // Finite element space for all variables together (total thermodynamic state)
169  ParFiniteElementSpace vfes(&pmesh, &fec, num_equation, Ordering::byNODES);
170 
171  // This example depends on this ordering of the space.
172  MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
173 
174  HYPRE_Int glob_size = vfes.GlobalTrueVSize();
175  if (mpi.Root()) { cout << "Number of unknowns: " << glob_size << endl; }
176 
177  // 8. Define the initial conditions, save the corresponding mesh and grid
178  // functions to a file. This can be opened with GLVis with the -gc option.
179 
180  // The solution u has components {density, x-momentum, y-momentum, energy}.
181  // These are stored contiguously in the BlockVector u_block.
182  Array<int> offsets(num_equation + 1);
183  for (int k = 0; k <= num_equation; k++) { offsets[k] = k * vfes.GetNDofs(); }
184  BlockVector u_block(offsets);
185 
186  // Momentum grid function on dfes for visualization.
187  ParGridFunction mom(&dfes, u_block.GetData() + offsets[1]);
188 
189  // Initialize the state.
190  VectorFunctionCoefficient u0(num_equation, InitialCondition);
191  ParGridFunction sol(&vfes, u_block.GetData());
192  sol.ProjectCoefficient(u0);
193 
194  // Output the initial solution.
195  {
196  ostringstream mesh_name;
197  mesh_name << "vortex-mesh." << setfill('0') << setw(6) << mpi.WorldRank();
198  ofstream mesh_ofs(mesh_name.str().c_str());
199  mesh_ofs.precision(precision);
200  mesh_ofs << pmesh;
201 
202  for (int k = 0; k < num_equation; k++)
203  {
204  ParGridFunction uk(&fes, u_block.GetBlock(k));
205  ostringstream sol_name;
206  sol_name << "vortex-" << k << "-init."
207  << setfill('0') << setw(6) << mpi.WorldRank();
208  ofstream sol_ofs(sol_name.str().c_str());
209  sol_ofs.precision(precision);
210  sol_ofs << uk;
211  }
212  }
213 
214  // 9. Set up the nonlinear form corresponding to the DG discretization of the
215  // flux divergence, and assemble the corresponding mass matrix.
216  MixedBilinearForm Aflux(&dfes, &fes);
217  Aflux.AddDomainIntegrator(new DomainIntegrator(dim));
218  Aflux.Assemble();
219 
220  ParNonlinearForm A(&vfes);
221  RiemannSolver rsolver;
222  A.AddInteriorFaceIntegrator(new FaceIntegrator(rsolver, dim));
223 
224  // 10. Define the time-dependent evolution operator describing the ODE
225  // right-hand side, and perform time-integration (looping over the time
226  // iterations, ti, with a time-step dt).
227  FE_Evolution euler(vfes, A, Aflux.SpMat());
228 
229  // Visualize the density
230  socketstream sout;
231  if (visualization)
232  {
233  char vishost[] = "localhost";
234  int visport = 19916;
235 
236  MPI_Barrier(pmesh.GetComm());
237  sout.open(vishost, visport);
238  if (!sout)
239  {
240  if (mpi.Root())
241  {
242  cout << "Unable to connect to GLVis server at "
243  << vishost << ':' << visport << endl;
244  }
245  visualization = false;
246  if (mpi.Root()) { cout << "GLVis visualization disabled.\n"; }
247  }
248  else
249  {
250  sout << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank() << "\n";
251  sout.precision(precision);
252  sout << "solution\n" << pmesh << mom;
253  sout << "pause\n";
254  sout << flush;
255  if (mpi.Root())
256  {
257  cout << "GLVis visualization paused."
258  << " Press space (in the GLVis window) to resume it.\n";
259  }
260  }
261  }
262 
263  // Determine the minimum element size.
264  double hmin;
265  if (cfl > 0)
266  {
267  double my_hmin = pmesh.GetElementSize(0, 1);
268  for (int i = 1; i < pmesh.GetNE(); i++)
269  {
270  my_hmin = min(pmesh.GetElementSize(i, 1), my_hmin);
271  }
272  // Reduce to find the global minimum element size
273  MPI_Allreduce(&my_hmin, &hmin, 1, MPI_DOUBLE, MPI_MIN, pmesh.GetComm());
274  }
275 
276  // Start the timer.
277  tic_toc.Clear();
278  tic_toc.Start();
279 
280  double t = 0.0;
281  euler.SetTime(t);
282  ode_solver->Init(euler);
283 
284  if (cfl > 0)
285  {
286  // Find a safe dt, using a temporary vector. Calling Mult() computes the
287  // maximum char speed at all quadrature points on all faces.
288  max_char_speed = 0.;
289  Vector z(sol.Size());
290  A.Mult(sol, z);
291  // Reduce to find the global maximum wave speed
292  {
293  double all_max_char_speed;
294  MPI_Allreduce(&max_char_speed, &all_max_char_speed,
295  1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
296  max_char_speed = all_max_char_speed;
297  }
298  dt = cfl * hmin / max_char_speed / (2*order+1);
299  }
300 
301  // Integrate in time.
302  bool done = false;
303  for (int ti = 0; !done; )
304  {
305  double dt_real = min(dt, t_final - t);
306 
307  ode_solver->Step(sol, t, dt_real);
308  if (cfl > 0)
309  {
310  // Reduce to find the global maximum wave speed
311  {
312  double all_max_char_speed;
313  MPI_Allreduce(&max_char_speed, &all_max_char_speed,
314  1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
315  max_char_speed = all_max_char_speed;
316  }
317  dt = cfl * hmin / max_char_speed / (2*order+1);
318  }
319  ti++;
320 
321  done = (t >= t_final - 1e-8*dt);
322  if (done || ti % vis_steps == 0)
323  {
324  if (mpi.Root())
325  {
326  cout << "time step: " << ti << ", time: " << t << endl;
327  }
328  if (visualization)
329  {
330  MPI_Barrier(pmesh.GetComm());
331  sout << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank() << "\n";
332  sout << "solution\n" << pmesh << mom << flush;
333  }
334  }
335  }
336 
337  tic_toc.Stop();
338  if (mpi.Root()) { cout << " done, " << tic_toc.RealTime() << "s." << endl; }
339 
340  // 11. Save the final solution. This output can be viewed later using GLVis:
341  // "glvis -np 4 -m vortex-mesh -g vortex-1-final".
342  for (int k = 0; k < num_equation; k++)
343  {
344  ParGridFunction uk(&fes, u_block.GetBlock(k));
345  ostringstream sol_name;
346  sol_name << "vortex-" << k << "-final."
347  << setfill('0') << setw(6) << mpi.WorldRank();
348  ofstream sol_ofs(sol_name.str().c_str());
349  sol_ofs.precision(precision);
350  sol_ofs << uk;
351  }
352 
353  // 12. Compute the L2 solution error summed for all components.
354  if (t_final == 2.0)
355  {
356  const double error = sol.ComputeLpError(2, u0);
357  if (mpi.Root()) { cout << "Solution error: " << error << endl; }
358  }
359 
360  // Free the used memory.
361  delete ode_solver;
362 
363  return 0;
364 }
StopWatch tic_toc
Definition: tic_toc.cpp:447
double RealTime()
Definition: tic_toc.cpp:426
const int num_equation
Definition: ex18p.cpp:54
int dim
Definition: ex3.cpp:47
double max_char_speed
Definition: ex18p.cpp:59
const double specific_heat_ratio
Definition: ex18p.cpp:55
L2_FECollection DG_FECollection
Definition: fe_coll.hpp:236
const double gas_constant
Definition: ex18p.cpp:56
int problem
Definition: ex18p.cpp:51
void InitialCondition(const Vector &x, Vector &y)
Definition: ex18.hpp:510
int main(int argc, char *argv[])
Definition: ex18p.cpp:61