MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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 -mf
11// ex18 -p 2 -r 0 -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.
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 // 1. Parse command-line options.
63 int problem = 1;
64 const real_t specific_heat_ratio = 1.4;
65 const real_t gas_constant = 1.0;
66
67 string mesh_file = "";
68 int IntOrderOffset = 1;
69 int ref_levels = 1;
70 int order = 3;
71 int ode_solver_type = 4;
72 real_t t_final = 2.0;
73 real_t dt = -0.01;
74 real_t cfl = 0.3;
75 bool visualization = true;
76 bool preassembleWeakDiv = 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. If not provided, then a periodic square"
85 " mesh will be used.");
86 args.AddOption(&problem, "-p", "--problem",
87 "Problem setup to use. See EulerInitialCondition().");
88 args.AddOption(&ref_levels, "-r", "--refine",
89 "Number of times to refine the mesh uniformly.");
90 args.AddOption(&order, "-o", "--order",
91 "Order (degree) of the finite elements.");
92 args.AddOption(&ode_solver_type, "-s", "--ode-solver",
93 "ODE solver: 1 - Forward Euler,\n\t"
94 " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
95 args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
96 args.AddOption(&dt, "-dt", "--time-step",
97 "Time step. Positive number skips CFL timestep calculation.");
98 args.AddOption(&cfl, "-c", "--cfl-number",
99 "CFL number for timestep calculation.");
100 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
101 "--no-visualization",
102 "Enable or disable GLVis visualization.");
103 args.AddOption(&preassembleWeakDiv, "-ea", "--element-assembly-divergence",
104 "-mf", "--matrix-free-divergence",
105 "Weak divergence assembly level\n"
106 " ea - Element assembly with interpolated F\n"
107 " mf - Nonlinear assembly in matrix-free manner");
108 args.AddOption(&vis_steps, "-vs", "--visualization-steps",
109 "Visualize every n-th timestep.");
110 args.ParseCheck();
111
112 // 2. Read the mesh from the given mesh file. When the user does not provide
113 // mesh file, use the default mesh file for the problem.
114 Mesh mesh = mesh_file.empty() ? EulerMesh(problem) : Mesh(mesh_file);
115 const int dim = mesh.Dimension();
116 const int num_equations = dim + 2;
117
118 // Refine the mesh to increase the resolution. In this example we do
119 // 'ref_levels' of uniform refinement, where 'ref_levels' is a command-line
120 // parameter.
121 for (int lev = 0; lev < ref_levels; lev++)
122 {
123 mesh.UniformRefinement();
124 }
125
126 // 3. 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 cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
138 return 3;
139 }
140
141 // 4. Define the discontinuous DG finite element space of the given
142 // polynomial order on the refined mesh.
143 DG_FECollection fec(order, dim);
144 // Finite element space for a scalar (thermodynamic quantity)
145 FiniteElementSpace fes(&mesh, &fec);
146 // Finite element space for a mesh-dim vector quantity (momentum)
147 FiniteElementSpace dfes(&mesh, &fec, dim, Ordering::byNODES);
148 // Finite element space for all variables together (total thermodynamic state)
149 FiniteElementSpace vfes(&mesh, &fec, num_equations, Ordering::byNODES);
150
151 // This example depends on this ordering of the space.
152 MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
153
154 cout << "Number of unknowns: " << vfes.GetVSize() << endl;
155
156 // 5. Define the initial conditions, save the corresponding mesh and grid
157 // functions to files. These can be opened with GLVis using:
158 // "glvis -m euler-mesh.mesh -g euler-1-init.gf" (for x-momentum).
159
160 // Initialize the state.
162 specific_heat_ratio,
163 gas_constant);
164 GridFunction sol(&vfes);
165 sol.ProjectCoefficient(u0);
166 GridFunction mom(&dfes, sol.GetData() + fes.GetNDofs());
167 // Output the initial solution.
168 {
169 ostringstream mesh_name;
170 mesh_name << "euler-mesh.mesh";
171 ofstream mesh_ofs(mesh_name.str().c_str());
172 mesh_ofs.precision(precision);
173 mesh_ofs << mesh;
174
175 for (int k = 0; k < num_equations; k++)
176 {
177 GridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
178 ostringstream sol_name;
179 sol_name << "euler-" << k << "-init.gf";
180 ofstream sol_ofs(sol_name.str().c_str());
181 sol_ofs.precision(precision);
182 sol_ofs << uk;
183 }
184 }
185
186 // 6. Set up the nonlinear form with euler flux and numerical flux
187 EulerFlux flux(dim, specific_heat_ratio);
188 RusanovFlux numericalFlux(flux);
190 vfes, std::unique_ptr<HyperbolicFormIntegrator>(
191 new HyperbolicFormIntegrator(numericalFlux, IntOrderOffset)),
192 preassembleWeakDiv);
193
194 // 7. Visualize momentum with its magnitude
195 socketstream sout;
196 if (visualization)
197 {
198 char vishost[] = "localhost";
199 int visport = 19916;
200
201 sout.open(vishost, visport);
202 if (!sout)
203 {
204 visualization = false;
205 cout << "Unable to connect to GLVis server at " << vishost << ':'
206 << visport << endl;
207 cout << "GLVis visualization disabled.\n";
208 }
209 else
210 {
211 sout.precision(precision);
212 // Plot magnitude of vector-valued momentum
213 sout << "solution\n" << mesh << mom;
214 sout << "window_title 'momentum, t = 0'\n";
215 sout << "view 0 0\n"; // view from top
216 sout << "keys jlm\n"; // turn off perspective and light, show mesh
217 sout << "pause\n";
218 sout << flush;
219 cout << "GLVis visualization paused."
220 << " Press space (in the GLVis window) to resume it.\n";
221 }
222 }
223
224 // 8. Time integration
225
226 // When dt is not specified, use CFL condition.
227 // Compute h_min and initial maximum characteristic speed
228 real_t hmin = infinity();
229 if (cfl > 0)
230 {
231 for (int i = 0; i < mesh.GetNE(); i++)
232 {
233 hmin = min(mesh.GetElementSize(i, 1), hmin);
234 }
235 // Find a safe dt, using a temporary vector. Calling Mult() computes the
236 // maximum char speed at all quadrature points on all faces (and all
237 // elements with -mf).
238 Vector z(sol.Size());
239 euler.Mult(sol, z);
240
241 real_t max_char_speed = euler.GetMaxCharSpeed();
242 dt = cfl * hmin / max_char_speed / (2 * order + 1);
243 }
244
245 // Start the timer.
246 tic_toc.Clear();
247 tic_toc.Start();
248
249 // Init time integration
250 real_t t = 0.0;
251 euler.SetTime(t);
252 ode_solver->Init(euler);
253
254 // Integrate in time.
255 bool done = false;
256 for (int ti = 0; !done;)
257 {
258 real_t dt_real = min(dt, t_final - t);
259
260 ode_solver->Step(sol, t, dt_real);
261 if (cfl > 0) // update time step size with CFL
262 {
263 real_t max_char_speed = euler.GetMaxCharSpeed();
264 dt = cfl * hmin / max_char_speed / (2 * order + 1);
265 }
266 ti++;
267
268 done = (t >= t_final - 1e-8 * dt);
269 if (done || ti % vis_steps == 0)
270 {
271 cout << "time step: " << ti << ", time: " << t << endl;
272 if (visualization)
273 {
274 sout << "window_title 'momentum, t = " << t << "'\n";
275 sout << "solution\n" << mesh << mom << flush;
276 }
277 }
278 }
279
280 tic_toc.Stop();
281 cout << " done, " << tic_toc.RealTime() << "s." << endl;
282
283 // 9. Save the final solution. This output can be viewed later using GLVis:
284 // "glvis -m euler-mesh-final.mesh -g euler-1-final.gf" (for x-momentum).
285 {
286 ostringstream mesh_name;
287 mesh_name << "euler-mesh-final.mesh";
288 ofstream mesh_ofs(mesh_name.str().c_str());
289 mesh_ofs.precision(precision);
290 mesh_ofs << mesh;
291
292 for (int k = 0; k < num_equations; k++)
293 {
294 GridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
295 ostringstream sol_name;
296 sol_name << "euler-" << k << "-final.gf";
297 ofstream sol_ofs(sol_name.str().c_str());
298 sol_ofs.precision(precision);
299 sol_ofs << uk;
300 }
301 }
302
303 // 10. Compute the L2 solution error summed for all components.
304 const real_t error = sol.ComputeLpError(2, u0);
305 cout << "Solution error: " << error << endl;
306
307 // Free the used memory.
308 delete ode_solver;
309
310 return 0;
311}
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
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
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:713
The classical forward Euler method.
Definition: ode.hpp:118
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
virtual real_t ComputeLpError(const real_t p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:3298
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2346
Abstract hyperbolic form integrator, (F(u, x), ∇v) and (F̂(u±, x, n))
Definition: hyperbolic.hpp:162
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Mesh data type.
Definition: mesh.hpp:56
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
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)
Definition: optparser.cpp:255
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
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 λ...
Definition: hyperbolic.hpp:249
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 visport
char vishost[]
int main()
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
RefCoord t[3]