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