MFEM  v3.2 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
60 class FE_Evolution : public TimeDependentOperator
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  int vis_steps = 5;
92
93  int precision = 8;
94  cout.precision(precision);
95
96  OptionsParser args(argc, argv);
98  "Mesh file to use.");
100  "Problem setup to use. See options in velocity_function().");
102  "Number of times to refine the mesh uniformly.");
104  "Order (degree) of the finite elements.");
106  "ODE solver: 1 - Forward Euler, 2 - RK2 SSP, 3 - RK3 SSP,"
107  " 4 - RK4, 6 - RK6.");
109  "Final time; start time is 0.");
111  "Time step.");
113  "--no-visualization",
114  "Enable or disable GLVis visualization.");
116  "--no-visit-datafiles",
117  "Save data files for VisIt (visit.llnl.gov) visualization.");
119  "Visualize every n-th timestep.");
120  args.Parse();
121  if (!args.Good())
122  {
123  args.PrintUsage(cout);
124  return 1;
125  }
126  args.PrintOptions(cout);
127
128  // 2. Read the mesh from the given mesh file. We can handle geometrically
129  // periodic meshes in this code.
130  Mesh *mesh = new Mesh(mesh_file, 1, 1);
131  int dim = mesh->Dimension();
132
133  // 3. Define the ODE solver used for time integration. Several explicit
134  // Runge-Kutta methods are available.
135  ODESolver *ode_solver = NULL;
136  switch (ode_solver_type)
137  {
138  case 1: ode_solver = new ForwardEulerSolver; break;
139  case 2: ode_solver = new RK2Solver(1.0); break;
140  case 3: ode_solver = new RK3SSPSolver; break;
141  case 4: ode_solver = new RK4Solver; break;
142  case 6: ode_solver = new RK6Solver; break;
143  default:
144  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
145  return 3;
146  }
147
148  // 4. Refine the mesh to increase the resolution. In this example we do
149  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
150  // command-line parameter. If the mesh is of NURBS type, we convert it to
151  // a (piecewise-polynomial) high-order mesh.
152  for (int lev = 0; lev < ref_levels; lev++)
153  {
154  mesh->UniformRefinement();
155  }
156  if (mesh->NURBSext)
157  {
158  mesh->SetCurvature(max(order, 1));
159  }
160  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
161
162  // 5. Define the discontinuous DG finite element space of the given
163  // polynomial order on the refined mesh.
164  DG_FECollection fec(order, dim);
165  FiniteElementSpace fes(mesh, &fec);
166
167  cout << "Number of unknowns: " << fes.GetVSize() << endl;
168
169  // 6. Set up and assemble the bilinear and linear forms corresponding to the
170  // DG discretization. The DGTraceIntegrator involves integrals over mesh
171  // interior faces.
175
176  BilinearForm m(&fes);
178  BilinearForm k(&fes);
181  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
183  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
184
185  LinearForm b(&fes);
187  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
188
189  m.Assemble();
190  m.Finalize();
191  int skip_zeros = 0;
192  k.Assemble(skip_zeros);
193  k.Finalize(skip_zeros);
194  b.Assemble();
195
196  // 7. Define the initial conditions, save the corresponding grid function to
197  // a file and (optionally) save data in the VisIt format and initialize
198  // GLVis visualization.
199  GridFunction u(&fes);
200  u.ProjectCoefficient(u0);
201
202  {
203  ofstream omesh("ex9.mesh");
204  omesh.precision(precision);
205  mesh->Print(omesh);
206  ofstream osol("ex9-init.gf");
207  osol.precision(precision);
208  u.Save(osol);
209  }
210
211  VisItDataCollection visit_dc("Example9", mesh);
212  visit_dc.RegisterField("solution", &u);
213  if (visit)
214  {
215  visit_dc.SetCycle(0);
216  visit_dc.SetTime(0.0);
217  visit_dc.Save();
218  }
219
220  socketstream sout;
221  if (visualization)
222  {
223  char vishost[] = "localhost";
224  int visport = 19916;
225  sout.open(vishost, visport);
226  if (!sout)
227  {
228  cout << "Unable to connect to GLVis server at "
229  << vishost << ':' << visport << endl;
230  visualization = false;
231  cout << "GLVis visualization disabled.\n";
232  }
233  else
234  {
235  sout.precision(precision);
236  sout << "solution\n" << *mesh << u;
237  sout << "pause\n";
238  sout << flush;
239  cout << "GLVis visualization paused."
240  << " Press space (in the GLVis window) to resume it.\n";
241  }
242  }
243
244  // 8. Define the time-dependent evolution operator describing the ODE
245  // right-hand side, and perform time-integration (looping over the time
246  // iterations, ti, with a time-step dt).
249
250  double t = 0.0;
251  for (int ti = 0; true; )
252  {
253  if (t >= t_final - dt/2)
254  {
255  break;
256  }
257
258  ode_solver->Step(u, t, dt);
259  ti++;
260
261  if (ti % vis_steps == 0)
262  {
263  cout << "time step: " << ti << ", time: " << t << endl;
264
265  if (visualization)
266  {
267  sout << "solution\n" << *mesh << u << flush;
268  }
269
270  if (visit)
271  {
272  visit_dc.SetCycle(ti);
273  visit_dc.SetTime(t);
274  visit_dc.Save();
275  }
276  }
277  }
278
279  // 9. Save the final solution. This output can be viewed later using GLVis:
280  // "glvis -m ex9.mesh -g ex9-final.gf".
281  {
282  ofstream osol("ex9-final.gf");
283  osol.precision(precision);
284  u.Save(osol);
285  }
286
287  // 10. Free the used memory.
288  delete ode_solver;
289  delete mesh;
290
291  return 0;
292 }
293
294
295 // Implementation of class FE_Evolution
296 FE_Evolution::FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b)
297  : TimeDependentOperator(_M.Size()), M(_M), K(_K), b(_b), z(_M.Size())
298 {
299  M_solver.SetPreconditioner(M_prec);
300  M_solver.SetOperator(M);
301
302  M_solver.iterative_mode = false;
303  M_solver.SetRelTol(1e-9);
304  M_solver.SetAbsTol(0.0);
305  M_solver.SetMaxIter(100);
306  M_solver.SetPrintLevel(0);
307 }
308
309 void FE_Evolution::Mult(const Vector &x, Vector &y) const
310 {
311  // y = M^{-1} (K x + b)
312  K.Mult(x, z);
313  z += b;
314  M_solver.Mult(z, y);
315 }
316
317
318 // Velocity coefficient
319 void velocity_function(const Vector &x, Vector &v)
320 {
321  int dim = x.Size();
322
323  // map to the reference [-1,1] domain
324  Vector X(dim);
325  for (int i = 0; i < dim; i++)
326  {
327  double center = (bb_min[i] + bb_max[i]) * 0.5;
328  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
329  }
330  switch (problem)
331  {
332  case 0:
333  {
334  // Translations in 1D, 2D, and 3D
335  switch (dim)
336  {
337  case 1: v(0) = 1.0; break;
338  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
339  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
340  break;
341  }
342  break;
343  }
344  case 1:
345  case 2:
346  {
347  // Clockwise rotation in 2D around the origin
348  const double w = M_PI/2;
349  switch (dim)
350  {
351  case 1: v(0) = 1.0; break;
352  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
353  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
354  }
355  break;
356  }
357  case 3:
358  {
359  // Clockwise twisting rotation in 2D around the origin
360  const double w = M_PI/2;
361  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
362  d = d*d;
363  switch (dim)
364  {
365  case 1: v(0) = 1.0; break;
366  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
367  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
368  }
369  break;
370  }
371  }
372 }
373
374 // Initial condition
375 double u0_function(const Vector &x)
376 {
377  int dim = x.Size();
378  // map to the reference [-1,1] domain
379  Vector X(dim);
380  for (int i = 0; i < dim; i++)
381  {
382  double center = (bb_min[i] + bb_max[i]) * 0.5;
383  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
384  }
385  switch (problem)
386  {
387  case 0:
388  case 1:
389  {
390  switch (dim)
391  {
392  case 1:
393  return exp(-40.*pow(X(0)-0.5,2));
394  case 2:
395  case 3:
396  {
397  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
398  if (dim == 3)
399  {
400  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
401  rx *= s;
402  ry *= s;
403  }
404  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
405  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
406  }
407  }
408  }
409  case 2:
410  {
411  double x_ = X(0), y_ = X(1), rho, phi;
412  rho = hypot(x_, y_) ;
413  phi = atan2(y_, x_);
414  return pow(sin(M_PI*rho),2)*sin(3*phi);
415  }
416  case 3:
417  {
418  const double f = M_PI;
419  return sin(f*X(0))*sin(f*X(1));
420  }
421  }
422  return 0.0;
423 }
424
425 // Inflow boundary condition (zero for the problems considered in this example)
426 double inflow_function(const Vector &x)
427 {
428  switch (problem)
429  {
430  case 0:
431  case 1:
432  case 2:
433  case 3: return 0.0;
434  }
435  return 0.0;
436 }
Vector bb_max
Definition: ex9.cpp:52
int GetVSize() const
Definition: fespace.hpp:161
Definition: solvers.hpp:110
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.
virtual void RegisterField(const char *field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
Base abstract class for time dependent operators: (x,t) -&gt; f(x,t)
Definition: operator.hpp:68
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:451
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Definition: mesh.cpp:87
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
virtual void Step(Vector &x, double &t, double &dt)=0
int Size() const
Returns the size of the vector.
Definition: vector.hpp:86
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
virtual void Save()
Save the collection and a VisIt root file.
virtual void Init(TimeDependentOperator &_f)
Definition: ode.hpp:30
int main(int argc, char *argv[])
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2116
int dim
Definition: ex3.cpp:47
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6402
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3020
Definition: linearform.cpp:29
int Dimension() const
Definition: mesh.hpp:523
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:52
virtual void Print(std::ostream &out=std::cout) const
Print the mesh to the given stream using the default MFEM mesh format.
Definition: mesh.cpp:6736
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:85
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:319
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
int problem
Definition: ex15.cpp:57
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:72
NURBSExtension * NURBSext
Definition: mesh.hpp:143
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1164
int open(const char hostname[], int port)
double u0_function(const Vector &x)
Definition: ex9.cpp:375
class for C-function coefficient
Vector data type.
Definition: vector.hpp:33
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
const SparseMatrix & SpMat() const
Returns a reference to the sparse matrix.
The classical forward Euler method.
Definition: ode.hpp:39
double inflow_function(const Vector &x)
Definition: ex9.cpp:426