MFEM  v3.4
Finite element discretization library
ex9.cpp
Go to the documentation of this file.
1 // MFEM Example 9
2 //
3 // Compile with: make ex9
4 //
5 // Sample runs:
6 // ex9 -m ../data/periodic-segment.mesh -p 0 -r 2 -dt 0.005
7 // ex9 -m ../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10
8 // ex9 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.01 -tf 10
9 // ex9 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
10 // ex9 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
11 // ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
12 // ex9 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.005 -tf 9
13 // ex9 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
14 // ex9 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
15 // ex9 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
16 // ex9 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
17 //
18 // Description: This example code solves the time-dependent advection equation
19 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
20 // u0(x)=u(0,x) is a given initial condition.
21 //
22 // The example demonstrates the use of Discontinuous Galerkin (DG)
23 // bilinear forms in MFEM (face integrators), the use of explicit
24 // ODE time integrators, the definition of periodic boundary
25 // conditions through periodic meshes, as well as the use of GLVis
26 // for persistent visualization of a time-evolving solution. The
27 // saving of time-dependent data files for external visualization
28 // with VisIt (visit.llnl.gov) is also illustrated.
29 
30 #include "mfem.hpp"
31 #include <fstream>
32 #include <iostream>
33 #include <algorithm>
34 
35 using namespace std;
36 using namespace mfem;
37 
38 // Choice for the problem setup. The fluid velocity, initial condition and
39 // inflow boundary condition are chosen based on this parameter.
40 int problem;
41 
42 // Velocity coefficient
43 void velocity_function(const Vector &x, Vector &v);
44 
45 // Initial condition
46 double u0_function(const Vector &x);
47 
48 // Inflow boundary condition
49 double inflow_function(const Vector &x);
50 
51 // Mesh bounding box
53 
54 
55 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
56  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
57  and advection matrices, and b describes the flow on the boundary. This can
58  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
59  used to evaluate the right-hand side. */
61 {
62 private:
63  SparseMatrix &M, &K;
64  const Vector &b;
65  DSmoother M_prec;
66  CGSolver M_solver;
67 
68  mutable Vector z;
69 
70 public:
71  FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b);
72 
73  virtual void Mult(const Vector &x, Vector &y) const;
74 
75  virtual ~FE_Evolution() { }
76 };
77 
78 
79 int main(int argc, char *argv[])
80 {
81  // 1. Parse command-line options.
82  problem = 0;
83  const char *mesh_file = "../data/periodic-hexagon.mesh";
84  int ref_levels = 2;
85  int order = 3;
86  int ode_solver_type = 4;
87  double t_final = 10.0;
88  double dt = 0.01;
89  bool visualization = true;
90  bool visit = false;
91  bool binary = false;
92  int vis_steps = 5;
93 
94  int precision = 8;
95  cout.precision(precision);
96 
97  OptionsParser args(argc, argv);
98  args.AddOption(&mesh_file, "-m", "--mesh",
99  "Mesh file to use.");
100  args.AddOption(&problem, "-p", "--problem",
101  "Problem setup to use. See options in velocity_function().");
102  args.AddOption(&ref_levels, "-r", "--refine",
103  "Number of times to refine the mesh uniformly.");
104  args.AddOption(&order, "-o", "--order",
105  "Order (degree) of the finite elements.");
106  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
107  "ODE solver: 1 - Forward Euler,\n\t"
108  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
109  args.AddOption(&t_final, "-tf", "--t-final",
110  "Final time; start time is 0.");
111  args.AddOption(&dt, "-dt", "--time-step",
112  "Time step.");
113  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
114  "--no-visualization",
115  "Enable or disable GLVis visualization.");
116  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
117  "--no-visit-datafiles",
118  "Save data files for VisIt (visit.llnl.gov) visualization.");
119  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
120  "--ascii-datafiles",
121  "Use binary (Sidre) or ascii format for VisIt data files.");
122  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
123  "Visualize every n-th timestep.");
124  args.Parse();
125  if (!args.Good())
126  {
127  args.PrintUsage(cout);
128  return 1;
129  }
130  args.PrintOptions(cout);
131 
132  // 2. Read the mesh from the given mesh file. We can handle geometrically
133  // periodic meshes in this code.
134  Mesh *mesh = new Mesh(mesh_file, 1, 1);
135  int dim = mesh->Dimension();
136 
137  // 3. Define the ODE solver used for time integration. Several explicit
138  // Runge-Kutta methods are available.
139  ODESolver *ode_solver = NULL;
140  switch (ode_solver_type)
141  {
142  case 1: ode_solver = new ForwardEulerSolver; break;
143  case 2: ode_solver = new RK2Solver(1.0); break;
144  case 3: ode_solver = new RK3SSPSolver; break;
145  case 4: ode_solver = new RK4Solver; break;
146  case 6: ode_solver = new RK6Solver; break;
147  default:
148  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
149  delete mesh;
150  return 3;
151  }
152 
153  // 4. Refine the mesh to increase the resolution. In this example we do
154  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
155  // command-line parameter. If the mesh is of NURBS type, we convert it to
156  // a (piecewise-polynomial) high-order mesh.
157  for (int lev = 0; lev < ref_levels; lev++)
158  {
159  mesh->UniformRefinement();
160  }
161  if (mesh->NURBSext)
162  {
163  mesh->SetCurvature(max(order, 1));
164  }
165  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
166 
167  // 5. Define the discontinuous DG finite element space of the given
168  // polynomial order on the refined mesh.
169  DG_FECollection fec(order, dim);
170  FiniteElementSpace fes(mesh, &fec);
171 
172  cout << "Number of unknowns: " << fes.GetVSize() << endl;
173 
174  // 6. Set up and assemble the bilinear and linear forms corresponding to the
175  // DG discretization. The DGTraceIntegrator involves integrals over mesh
176  // interior faces.
180 
181  BilinearForm m(&fes);
183  BilinearForm k(&fes);
184  k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
186  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
188  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
189 
190  LinearForm b(&fes);
192  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
193 
194  m.Assemble();
195  m.Finalize();
196  int skip_zeros = 0;
197  k.Assemble(skip_zeros);
198  k.Finalize(skip_zeros);
199  b.Assemble();
200 
201  // 7. Define the initial conditions, save the corresponding grid function to
202  // a file and (optionally) save data in the VisIt format and initialize
203  // GLVis visualization.
204  GridFunction u(&fes);
205  u.ProjectCoefficient(u0);
206 
207  {
208  ofstream omesh("ex9.mesh");
209  omesh.precision(precision);
210  mesh->Print(omesh);
211  ofstream osol("ex9-init.gf");
212  osol.precision(precision);
213  u.Save(osol);
214  }
215 
216  // Create data collection for solution output: either VisItDataCollection for
217  // ascii data files, or SidreDataCollection for binary data files.
218  DataCollection *dc = NULL;
219  if (visit)
220  {
221  if (binary)
222  {
223 #ifdef MFEM_USE_SIDRE
224  dc = new SidreDataCollection("Example9", mesh);
225 #else
226  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
227 #endif
228  }
229  else
230  {
231  dc = new VisItDataCollection("Example9", mesh);
232  dc->SetPrecision(precision);
233  }
234  dc->RegisterField("solution", &u);
235  dc->SetCycle(0);
236  dc->SetTime(0.0);
237  dc->Save();
238  }
239 
240  socketstream sout;
241  if (visualization)
242  {
243  char vishost[] = "localhost";
244  int visport = 19916;
245  sout.open(vishost, visport);
246  if (!sout)
247  {
248  cout << "Unable to connect to GLVis server at "
249  << vishost << ':' << visport << endl;
250  visualization = false;
251  cout << "GLVis visualization disabled.\n";
252  }
253  else
254  {
255  sout.precision(precision);
256  sout << "solution\n" << *mesh << u;
257  sout << "pause\n";
258  sout << flush;
259  cout << "GLVis visualization paused."
260  << " Press space (in the GLVis window) to resume it.\n";
261  }
262  }
263 
264  // 8. Define the time-dependent evolution operator describing the ODE
265  // right-hand side, and perform time-integration (looping over the time
266  // iterations, ti, with a time-step dt).
267  FE_Evolution adv(m.SpMat(), k.SpMat(), b);
268 
269  double t = 0.0;
270  adv.SetTime(t);
271  ode_solver->Init(adv);
272 
273  bool done = false;
274  for (int ti = 0; !done; )
275  {
276  double dt_real = min(dt, t_final - t);
277  ode_solver->Step(u, t, dt_real);
278  ti++;
279 
280  done = (t >= t_final - 1e-8*dt);
281 
282  if (done || ti % vis_steps == 0)
283  {
284  cout << "time step: " << ti << ", time: " << t << endl;
285 
286  if (visualization)
287  {
288  sout << "solution\n" << *mesh << u << flush;
289  }
290 
291  if (visit)
292  {
293  dc->SetCycle(ti);
294  dc->SetTime(t);
295  dc->Save();
296  }
297  }
298  }
299 
300  // 9. Save the final solution. This output can be viewed later using GLVis:
301  // "glvis -m ex9.mesh -g ex9-final.gf".
302  {
303  ofstream osol("ex9-final.gf");
304  osol.precision(precision);
305  u.Save(osol);
306  }
307 
308  // 10. Free the used memory.
309  delete ode_solver;
310  delete dc;
311 
312  return 0;
313 }
314 
315 
316 // Implementation of class FE_Evolution
318  : TimeDependentOperator(_M.Size()), M(_M), K(_K), b(_b), z(_M.Size())
319 {
320  M_solver.SetPreconditioner(M_prec);
321  M_solver.SetOperator(M);
322 
323  M_solver.iterative_mode = false;
324  M_solver.SetRelTol(1e-9);
325  M_solver.SetAbsTol(0.0);
326  M_solver.SetMaxIter(100);
327  M_solver.SetPrintLevel(0);
328 }
329 
330 void FE_Evolution::Mult(const Vector &x, Vector &y) const
331 {
332  // y = M^{-1} (K x + b)
333  K.Mult(x, z);
334  z += b;
335  M_solver.Mult(z, y);
336 }
337 
338 
339 // Velocity coefficient
340 void velocity_function(const Vector &x, Vector &v)
341 {
342  int dim = x.Size();
343 
344  // map to the reference [-1,1] domain
345  Vector X(dim);
346  for (int i = 0; i < dim; i++)
347  {
348  double center = (bb_min[i] + bb_max[i]) * 0.5;
349  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
350  }
351 
352  switch (problem)
353  {
354  case 0:
355  {
356  // Translations in 1D, 2D, and 3D
357  switch (dim)
358  {
359  case 1: v(0) = 1.0; break;
360  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
361  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
362  break;
363  }
364  break;
365  }
366  case 1:
367  case 2:
368  {
369  // Clockwise rotation in 2D around the origin
370  const double w = M_PI/2;
371  switch (dim)
372  {
373  case 1: v(0) = 1.0; break;
374  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
375  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
376  }
377  break;
378  }
379  case 3:
380  {
381  // Clockwise twisting rotation in 2D around the origin
382  const double w = M_PI/2;
383  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
384  d = d*d;
385  switch (dim)
386  {
387  case 1: v(0) = 1.0; break;
388  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
389  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
390  }
391  break;
392  }
393  }
394 }
395 
396 // Initial condition
397 double u0_function(const Vector &x)
398 {
399  int dim = x.Size();
400 
401  // map to the reference [-1,1] domain
402  Vector X(dim);
403  for (int i = 0; i < dim; i++)
404  {
405  double center = (bb_min[i] + bb_max[i]) * 0.5;
406  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
407  }
408 
409  switch (problem)
410  {
411  case 0:
412  case 1:
413  {
414  switch (dim)
415  {
416  case 1:
417  return exp(-40.*pow(X(0)-0.5,2));
418  case 2:
419  case 3:
420  {
421  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
422  if (dim == 3)
423  {
424  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
425  rx *= s;
426  ry *= s;
427  }
428  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
429  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
430  }
431  }
432  }
433  case 2:
434  {
435  double x_ = X(0), y_ = X(1), rho, phi;
436  rho = hypot(x_, y_);
437  phi = atan2(y_, x_);
438  return pow(sin(M_PI*rho),2)*sin(3*phi);
439  }
440  case 3:
441  {
442  const double f = M_PI;
443  return sin(f*X(0))*sin(f*X(1));
444  }
445  }
446  return 0.0;
447 }
448 
449 // Inflow boundary condition (zero for the problems considered in this example)
450 double inflow_function(const Vector &x)
451 {
452  switch (problem)
453  {
454  case 0:
455  case 1:
456  case 2:
457  case 3: return 0.0;
458  }
459  return 0.0;
460 }
Vector bb_max
Definition: ex9.cpp:52
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:486
Conjugate gradient method.
Definition: solvers.hpp:111
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
const SparseMatrix & SpMat() const
Returns a reference to the sparse matrix.
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.hpp:37
int Dimension() const
Definition: mesh.hpp:645
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Base abstract class for time dependent operators.
Definition: operator.hpp:151
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Definition: mesh.cpp:108
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:51
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:295
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:181
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
bool Good() const
Definition: optparser.hpp:120
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:272
STL namespace.
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:104
Data collection with Sidre routines following the Conduit mesh blueprint specification.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
int dim
Definition: ex3.cpp:47
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:72
Data type sparse matrix.
Definition: sparsemat.hpp:38
virtual void Save()
Save the collection to disk.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6741
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3355
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: ex18.hpp:130
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator.
Definition: linearform.cpp:38
int problem
Definition: ex9.cpp:40
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:52
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1034
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
void SetAbsTol(double atol)
Definition: solvers.hpp:62
void SetRelTol(double rtol)
Definition: solvers.hpp:61
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:340
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:66
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)
Definition: optparser.hpp:74
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:134
virtual ~FE_Evolution()
Definition: ex9.cpp:75
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:176
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
int main(int argc, char *argv[])
Definition: ex9.cpp:79
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:253
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1377
int open(const char hostname[], int port)
double u0_function(const Vector &x)
Definition: ex9.cpp:397
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:125
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:93
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2349
The classical forward Euler method.
Definition: ode.hpp:101
double inflow_function(const Vector &x)
Definition: ex9.cpp:450
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:128
alpha (q . grad u, v)