MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex18p.cpp
Go to the documentation of this file.
1// MFEM Example 18 - Parallel Version
2//
3// Compile with: make ex18p
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 -mf
11// mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 3 -s 3 -mf
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 in parallel.
16//
17// (u_t, v)_T - (F(u), ∇ v)_T + <F̂(u,n), [[v]]>_F = 0
18//
19// where (⋅,⋅)_T is volume integration, and <⋅,⋅>_F is face
20// integration, F is the Euler flux function, and F̂ is the
21// numerical flux.
22//
23// Specifically, it solves for an exact solution of the equations
24// whereby a vortex is transported by a uniform flow. Since all
25// boundaries are periodic here, the method's accuracy can be
26// assessed by measuring the difference between the solution and
27// the initial condition at a later time when the vortex returns
28// to its initial location.
29//
30// Note that as the order of the spatial discretization increases,
31// the timestep must become smaller. This example currently uses a
32// simple estimate derived by Cockburn and Shu for the 1D RKDG
33// method. An additional factor can be tuned by passing the --cfl
34// (or -c shorter) flag.
35//
36// The example demonstrates usage of DGHyperbolicConservationLaws
37// that wraps NonlinearFormIntegrators containing element and face
38// integration schemes. In this case the system also involves an
39// external approximate Riemann solver for the DG interface flux.
40// By default, weak-divergence is pre-assembled in element-wise
41// manner, which corresponds to (I_h(F(u_h)), ∇ v). This yields
42// better performance and similar accuracy for the included test
43// problems. This can be turned off and use nonlinear assembly
44// similar to matrix-free assembly when -mf flag is provided.
45// It also demonstrates how to use GLVis for in-situ visualization
46// of vector grid function and how to set top-view.
47//
48// We recommend viewing examples 9, 14 and 17 before viewing this
49// example.
50
51#include "mfem.hpp"
52#include <fstream>
53#include <iostream>
54#include <sstream>
55#include "ex18.hpp"
56
57using namespace std;
58using namespace mfem;
59
60int main(int argc, char *argv[])
61{
62 // 0. Parallel setup
63 Mpi::Init(argc, argv);
64 const int numProcs = Mpi::WorldSize();
65 const int myRank = Mpi::WorldRank();
67
68 // 1. Parse command-line options.
69 int problem = 1;
70 const real_t specific_heat_ratio = 1.4;
71 const real_t gas_constant = 1.0;
72
73 string mesh_file = "";
74 int IntOrderOffset = 1;
75 int ser_ref_levels = 0;
76 int par_ref_levels = 1;
77 int order = 3;
78 int ode_solver_type = 4;
79 real_t t_final = 2.0;
80 real_t dt = -0.01;
81 real_t cfl = 0.3;
82 bool visualization = true;
83 bool preassembleWeakDiv = true;
84 int vis_steps = 50;
85
86 int precision = 8;
87 cout.precision(precision);
88
89 OptionsParser args(argc, argv);
91 "Mesh file to use. If not provided, then a periodic square"
92 " mesh will be used.");
94 "Problem setup to use. See EulerInitialCondition().");
96 "Number of times to refine the serial mesh uniformly.");
98 "Number of times to refine the parallel mesh uniformly.");
100 "Order (degree) of the finite elements.");
102 "ODE solver: 1 - Forward Euler,\n\t"
103 " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
104 args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
106 "Time step. Positive number skips CFL timestep calculation.");
108 "CFL number for timestep calculation.");
110 "--no-visualization",
111 "Enable or disable GLVis visualization.");
113 "-mf", "--matrix-free-divergence",
114 "Weak divergence assembly level\n"
115 " ea - Element assembly with interpolated F\n"
116 " mf - Nonlinear assembly in matrix-free manner");
118 "Visualize every n-th timestep.");
119 args.ParseCheck();
120
121 // 2. Read the mesh from the given mesh file. When the user does not provide
122 // mesh file, use the default mesh file for the problem.
123 Mesh mesh = mesh_file.empty() ? EulerMesh(problem) : Mesh(mesh_file);
124 const int dim = mesh.Dimension();
125 const int num_equations = dim + 2;
126
127 // Refine the mesh to increase the resolution. In this example we do
128 // 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is a
129 // command-line parameter.
130 for (int lev = 0; lev < ser_ref_levels; lev++)
131 {
132 mesh.UniformRefinement();
133 }
134
135 // Define a parallel mesh by a partitioning of the serial mesh. Refine this
136 // mesh further in parallel to increase the resolution. Once the parallel
137 // mesh is defined, the serial mesh can be deleted.
138 ParMesh pmesh = ParMesh(MPI_COMM_WORLD, mesh);
139 mesh.Clear();
140
141 // Refine the mesh to increase the resolution. In this example we do
142 // 'par_ref_levels' of uniform refinement, where 'par_ref_levels' is a
143 // command-line parameter.
144 for (int lev = 0; lev < par_ref_levels; lev++)
145 {
146 pmesh.UniformRefinement();
147 }
148
149 // 3. Define the ODE solver used for time integration. Several explicit
150 // Runge-Kutta methods are available.
151 ODESolver *ode_solver = NULL;
152 switch (ode_solver_type)
153 {
154 case 1: ode_solver = new ForwardEulerSolver; break;
155 case 2: ode_solver = new RK2Solver(1.0); break;
156 case 3: ode_solver = new RK3SSPSolver; break;
157 case 4: ode_solver = new RK4Solver; break;
158 case 6: ode_solver = new RK6Solver; break;
159 default:
160 cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
161 return 3;
162 }
163
164 // 4. Define the discontinuous DG finite element space of the given
165 // polynomial order on the refined mesh.
166 DG_FECollection fec(order, dim);
167 // Finite element space for a scalar (thermodynamic quantity)
168 ParFiniteElementSpace fes(&pmesh, &fec);
169 // Finite element space for a mesh-dim vector quantity (momentum)
170 ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
171 // Finite element space for all variables together (total thermodynamic state)
172 ParFiniteElementSpace vfes(&pmesh, &fec, num_equations, Ordering::byNODES);
173
174 // This example depends on this ordering of the space.
175 MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
176
177 HYPRE_BigInt glob_size = vfes.GlobalTrueVSize();
178 if (Mpi::Root())
179 {
180 cout << "Number of unknowns: " << glob_size << endl;
181 }
182
183 // 5. Define the initial conditions, save the corresponding mesh and grid
184 // functions to files. These can be opened with GLVis using:
185 // "glvis -np 4 -m euler-mesh -g euler-1-init" (for x-momentum).
186
187 // Initialize the state.
189 specific_heat_ratio,
190 gas_constant);
191 ParGridFunction sol(&vfes);
192 sol.ProjectCoefficient(u0);
193 ParGridFunction mom(&dfes, sol.GetData() + fes.GetNDofs());
194 // Output the initial solution.
195 {
196 ostringstream mesh_name;
197 mesh_name << "euler-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_equations; k++)
203 {
204 ParGridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
205 ostringstream sol_name;
206 sol_name << "euler-" << k << "-init." << setfill('0') << setw(6)
207 << Mpi::WorldRank();
208 ofstream sol_ofs(sol_name.str().c_str());
209 sol_ofs.precision(precision);
210 sol_ofs << uk;
211 }
212 }
213
214 // 6. Set up the nonlinear form with euler flux and numerical flux
215 EulerFlux flux(dim, specific_heat_ratio);
216 RusanovFlux numericalFlux(flux);
218 vfes, std::unique_ptr<HyperbolicFormIntegrator>(
219 new HyperbolicFormIntegrator(numericalFlux, IntOrderOffset)),
220 preassembleWeakDiv);
221
222 // 7. Visualize momentum with its magnitude
223 socketstream sout;
224 if (visualization)
225 {
226 char vishost[] = "localhost";
227 int visport = 19916;
228
229 sout.open(vishost, visport);
230 if (!sout)
231 {
232 visualization = false;
233 if (Mpi::Root())
234 {
235 cout << "Unable to connect to GLVis server at " << vishost << ':'
236 << visport << endl;
237 cout << "GLVis visualization disabled.\n";
238 }
239 }
240 else
241 {
242 sout.precision(precision);
243 // Plot magnitude of vector-valued momentum
244 sout << "parallel " << numProcs << " " << myRank << "\n";
245 sout << "solution\n" << pmesh << mom;
246 sout << "window_title 'momentum, t = 0'\n";
247 sout << "view 0 0\n"; // view from top
248 sout << "keys jlm\n"; // turn off perspective and light, show mesh
249 sout << "pause\n";
250 sout << flush;
251 if (Mpi::Root())
252 {
253 cout << "GLVis visualization paused."
254 << " Press space (in the GLVis window) to resume it.\n";
255 }
256 MPI_Barrier(pmesh.GetComm());
257 }
258 }
259
260 // 8. Time integration
261
262 // When dt is not specified, use CFL condition.
263 // Compute h_min and initial maximum characteristic speed
264 real_t hmin = infinity();
265 if (cfl > 0)
266 {
267 for (int i = 0; i < pmesh.GetNE(); i++)
268 {
269 hmin = min(pmesh.GetElementSize(i, 1), hmin);
270 }
271 MPI_Allreduce(MPI_IN_PLACE, &hmin, 1, MPITypeMap<real_t>::mpi_type, MPI_MIN,
272 pmesh.GetComm());
273 // Find a safe dt, using a temporary vector. Calling Mult() computes the
274 // maximum char speed at all quadrature points on all faces (and all
275 // elements with -mf).
276 Vector z(sol.Size());
277 euler.Mult(sol, z);
278
279 real_t max_char_speed = euler.GetMaxCharSpeed();
280 MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1, MPITypeMap<real_t>::mpi_type,
281 MPI_MAX,
282 pmesh.GetComm());
283 dt = cfl * hmin / max_char_speed / (2 * order + 1);
284 }
285
286 // Start the timer.
287 tic_toc.Clear();
288 tic_toc.Start();
289
290 // Init time integration
291 real_t t = 0.0;
292 euler.SetTime(t);
293 ode_solver->Init(euler);
294
295 // Integrate in time.
296 bool done = false;
297 for (int ti = 0; !done;)
298 {
299 real_t dt_real = min(dt, t_final - t);
300
301 ode_solver->Step(sol, t, dt_real);
302 if (cfl > 0) // update time step size with CFL
303 {
304 real_t max_char_speed = euler.GetMaxCharSpeed();
305 MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1, MPITypeMap<real_t>::mpi_type,
306 MPI_MAX,
307 pmesh.GetComm());
308 dt = cfl * hmin / max_char_speed / (2 * order + 1);
309 }
310 ti++;
311
312 done = (t >= t_final - 1e-8 * dt);
313 if (done || ti % vis_steps == 0)
314 {
315 if (Mpi::Root())
316 {
317 cout << "time step: " << ti << ", time: " << t << endl;
318 }
319 if (visualization)
320 {
321 sout << "window_title 'momentum, t = " << t << "'\n";
322 sout << "parallel " << numProcs << " " << myRank << "\n";
323 sout << "solution\n" << pmesh << mom << flush;
324 }
325 }
326 }
327
328 tic_toc.Stop();
329 if (Mpi::Root())
330 {
331 cout << " done, " << tic_toc.RealTime() << "s." << endl;
332 }
333
334 // 9. Save the final solution. This output can be viewed later using GLVis:
335 // "glvis -np 4 -m euler-mesh-final -g euler-1-final" (for x-momentum).
336 {
337 ostringstream mesh_name;
338 mesh_name << "euler-mesh-final." << setfill('0') << setw(6)
339 << Mpi::WorldRank();
340 ofstream mesh_ofs(mesh_name.str().c_str());
341 mesh_ofs.precision(precision);
342 mesh_ofs << pmesh;
343
344 for (int k = 0; k < num_equations; k++)
345 {
346 ParGridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
347 ostringstream sol_name;
348 sol_name << "euler-" << k << "-final." << setfill('0') << setw(6)
349 << Mpi::WorldRank();
350 ofstream sol_ofs(sol_name.str().c_str());
351 sol_ofs.precision(precision);
352 sol_ofs << uk;
353 }
354 }
355
356 // 10. Compute the L2 solution error summed for all components.
357 const real_t error = sol.ComputeLpError(2, u0);
358 if (Mpi::Root())
359 {
360 cout << "Solution error: " << error << endl;
361 }
362
363 // Free the used memory.
364 delete ode_solver;
365
366 return 0;
367}
Time dependent DG operator for hyperbolic conservation laws.
Definition ex18.hpp:32
void Mult(const Vector &x, Vector &y) const override
Apply nonlinear form to obtain M⁻¹(DIVF + JUMP HAT(F))
Definition ex18.hpp:168
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:710
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:725
The classical forward Euler method.
Definition ode.hpp:118
Abstract hyperbolic form integrator, (F(u, x), ∇v) and (F̂(u±, x, n))
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Definition mesh.cpp:107
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
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 class for solving systems of ODEs: dx/dt = f(x,t)
Definition ode.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
void ParseCheck(std::ostream &out=mfem::out)
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
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeLpError(const real_t p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition ode.hpp:151
The classical explicit forth-order Runge-Kutta method, RK4.
Definition ode.hpp:164
Rusanov flux, also known as local Lax-Friedrichs, F̂ n = ½(F(u⁺,x)n + F(u⁻,x)n) - ½λ(u⁺ - u⁻) where λ...
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Stop()
Stop the stopwatch.
Definition tic_toc.cpp:422
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.
Definition tic_toc.cpp:406
virtual void SetTime(const real_t t_)
Set the current time.
Definition operator.hpp:360
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int problem
Definition ex15.cpp:62
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
const int visport
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:45
Mesh EulerMesh(const int problem)
Definition ex18.hpp:299
VectorFunctionCoefficient EulerInitialCondition(const int problem, const real_t specific_heat_ratio, const real_t gas_constant)
Definition ex18.hpp:317
StopWatch tic_toc
Definition tic_toc.cpp:450
float real_t
Definition config.hpp:43
const char vishost[]
RefCoord t[3]
Helper struct to convert a C++ type to an MPI type.