61 int main(
int argc,
char *argv[])
65 const char *mesh_file =
"../data/periodic-square.mesh";
68 int ode_solver_type = 4;
72 bool visualization =
true;
76 cout.precision(precision);
78 OptionsParser args(argc, argv);
79 args.AddOption(&mesh_file,
"-m",
"--mesh",
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",
98 "Enable or disable GLVis visualization.");
99 args.AddOption(&vis_steps,
"-vs",
"--visualization-steps",
100 "Visualize every n-th timestep.");
105 args.PrintUsage(cout);
108 args.PrintOptions(cout);
112 Mesh mesh(mesh_file, 1, 1);
113 const int dim = mesh.Dimension();
115 MFEM_ASSERT(dim == 2,
"Need a two-dimensional mesh for the problem definition");
119 ODESolver *ode_solver = NULL;
120 switch (ode_solver_type)
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;
128 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
135 for (
int lev = 0; lev < ref_levels; lev++)
137 mesh.UniformRefinement();
144 FiniteElementSpace fes(&mesh, &fec);
146 FiniteElementSpace dfes(&mesh, &fec, dim, Ordering::byNODES);
148 FiniteElementSpace vfes(&mesh, &fec,
num_equation, Ordering::byNODES);
151 MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES,
"");
153 cout <<
"Number of unknowns: " << vfes.GetVSize() << endl;
161 for (
int k = 0; k <=
num_equation; k++) { offsets[k] = k * vfes.GetNDofs(); }
162 BlockVector u_block(offsets);
165 GridFunction mom(&dfes, u_block.GetData() + offsets[1]);
169 GridFunction sol(&vfes, u_block.GetData());
170 sol.ProjectCoefficient(u0);
174 ofstream mesh_ofs(
"vortex.mesh");
175 mesh_ofs.precision(precision);
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);
191 MixedBilinearForm Aflux(&dfes, &fes);
195 NonlinearForm A(&vfes);
211 sout.open(vishost, visport);
214 cout <<
"Unable to connect to GLVis server at "
215 << vishost <<
':' << visport << endl;
216 visualization =
false;
217 cout <<
"GLVis visualization disabled.\n";
221 sout.precision(precision);
222 sout <<
"solution\n" << mesh << mom;
225 cout <<
"GLVis visualization paused."
226 <<
" Press space (in the GLVis window) to resume it.\n";
234 hmin = mesh.GetElementSize(0, 1);
235 for (
int i = 1; i < mesh.GetNE(); i++)
237 hmin = min(mesh.GetElementSize(i, 1), hmin);
247 ode_solver->Init(euler);
256 dt = cfl * hmin / max_char_speed / (2*order+1);
261 for (
int ti = 0; !done; )
263 double dt_real = min(dt, t_final - t);
265 ode_solver->Step(sol, t, dt_real);
272 done = (t >= t_final - 1e-8*dt);
273 if (done || ti % vis_steps == 0)
275 cout <<
"time step: " << ti <<
", time: " << t << endl;
278 sout <<
"solution\n" << mesh << mom << flush;
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);
301 const double error = sol.ComputeLpError(2, u0);
302 cout <<
"Solution error: " << error << endl;
void Stop()
Stop the stopwatch.
const double specific_heat_ratio
void Start()
Start the stopwatch. The elapsed time is not cleared.
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
void InitialCondition(const Vector &x, Vector &y)
const double gas_constant
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.