MFEM  v3.1 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/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
12 // ex9 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
13 // ex9 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
14 //
15 // Description: This example code solves the time-dependent advection equation
16 // du/dt = v.grad(u), where v is a given fluid velocity, and
17 // u0(x)=u(0,x) is a given initial condition.
18 //
19 // The example demonstrates the use of Discontinuous Galerkin (DG)
20 // bilinear forms in MFEM (face integrators), the use of explicit
21 // ODE time integrators, the definition of periodic boundary
22 // conditions through periodic meshes, as well as the use of GLVis
23 // for persistent visualization of a time-evolving solution. The
24 // saving of time-dependent data files for external visualization
25 // with VisIt (visit.llnl.gov) is also illustrated.
26
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30 #include <algorithm>
31
32 using namespace std;
33 using namespace mfem;
34
35 // Choice for the problem setup. The fluid velocity, initial condition and
36 // inflow boundary condition are chosen based on this parameter.
37 int problem;
38
39 // Velocity coefficient
40 void velocity_function(const Vector &x, Vector &v);
41
42 // Initial condition
43 double u0_function(const Vector &x);
44
45 // Inflow boundary condition
46 double inflow_function(const Vector &x);
47
48
54 class FE_Evolution : public TimeDependentOperator
55 {
56 private:
57  SparseMatrix &M, &K;
58  const Vector &b;
59  DSmoother M_prec;
60  CGSolver M_solver;
61
62  mutable Vector z;
63
64 public:
65  FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b);
66
67  virtual void Mult(const Vector &x, Vector &y) const;
68
69  virtual ~FE_Evolution() { }
70 };
71
72
73 int main(int argc, char *argv[])
74 {
75  // 1. Parse command-line options.
76  problem = 0;
77  const char *mesh_file = "../data/periodic-hexagon.mesh";
78  int ref_levels = 2;
79  int order = 3;
80  int ode_solver_type = 4;
81  double t_final = 10.0;
82  double dt = 0.01;
83  bool visualization = true;
84  bool visit = false;
85  int vis_steps = 5;
86
87  int precision = 8;
88  cout.precision(precision);
89
90  OptionsParser args(argc, argv);
92  "Mesh file to use.");
94  "Problem setup to use. See options in velocity_function().");
96  "Number of times to refine the mesh uniformly.");
98  "Order (degree) of the finite elements.");
100  "ODE solver: 1 - Forward Euler, 2 - RK2 SSP, 3 - RK3 SSP,"
101  " 4 - RK4, 6 - RK6.");
103  "Final time; start time is 0.");
105  "Time step.");
107  "--no-visualization",
108  "Enable or disable GLVis visualization.");
110  "--no-visit-datafiles",
111  "Save data files for VisIt (visit.llnl.gov) visualization.");
113  "Visualize every n-th timestep.");
114  args.Parse();
115  if (!args.Good())
116  {
117  args.PrintUsage(cout);
118  return 1;
119  }
120  args.PrintOptions(cout);
121
122  // 2. Read the mesh from the given mesh file. We can handle geometrically
123  // periodic meshes in this code.
124  Mesh *mesh;
125  ifstream imesh(mesh_file);
126  if (!imesh)
127  {
128  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
129  return 2;
130  }
131  mesh = new Mesh(imesh, 1, 1);
132  imesh.close();
133  int dim = mesh->Dimension();
134
135  // 3. Define the ODE solver used for time integration. Several explicit
136  // Runge-Kutta methods are available.
137  ODESolver *ode_solver = NULL;
138  switch (ode_solver_type)
139  {
140  case 1: ode_solver = new ForwardEulerSolver; break;
141  case 2: ode_solver = new RK2Solver(1.0); break;
142  case 3: ode_solver = new RK3SSPSolver; break;
143  case 4: ode_solver = new RK4Solver; break;
144  case 6: ode_solver = new RK6Solver; break;
145  default:
146  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
147  return 3;
148  }
149
150  // 4. Refine the mesh to increase the resolution. In this example we do
151  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
152  // command-line parameter. If the mesh is of NURBS type, we convert it to
153  // a (piecewise-polynomial) high-order mesh.
154  for (int lev = 0; lev < ref_levels; lev++)
155  {
156  mesh->UniformRefinement();
157  }
158
159  if (mesh->NURBSext)
160  {
161  int mesh_order = std::max(order, 1);
162  FiniteElementCollection *mfec = new H1_FECollection(mesh_order, dim);
163  FiniteElementSpace *mfes = new FiniteElementSpace(mesh, mfec, dim);
164  mesh->SetNodalFESpace(mfes);
165  mesh->GetNodes()->MakeOwner(mfec);
166  }
167
168  // 5. Define the discontinuous DG finite element space of the given
169  // polynomial order on the refined mesh.
170  DG_FECollection fec(order, dim);
171  FiniteElementSpace fes(mesh, &fec);
172
173  cout << "Number of unknowns: " << fes.GetVSize() << endl;
174
175  // 6. Set up and assemble the bilinear and linear forms corresponding to the
176  // DG discretization. The DGTraceIntegrator involves integrals over mesh
177  // interior faces.
181
182  BilinearForm m(&fes);
184  BilinearForm k(&fes);
187  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
189  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
190
191  LinearForm b(&fes);
193  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
194
195  m.Assemble();
196  m.Finalize();
197  int skip_zeros = 0;
198  k.Assemble(skip_zeros);
199  k.Finalize(skip_zeros);
200  b.Assemble();
201
202  // 7. Define the initial conditions, save the corresponding grid function to
203  // a file and (optionally) save data in the VisIt format and initialize
204  // GLVis visualization.
205  GridFunction u(&fes);
206  u.ProjectCoefficient(u0);
207
208  {
209  ofstream omesh("ex9.mesh");
210  omesh.precision(precision);
211  mesh->Print(omesh);
212  ofstream osol("ex9-init.gf");
213  osol.precision(precision);
214  u.Save(osol);
215  }
216
217  VisItDataCollection visit_dc("Example9", mesh);
218  visit_dc.RegisterField("solution", &u);
219  if (visit)
220  {
221  visit_dc.SetCycle(0);
222  visit_dc.SetTime(0.0);
223  visit_dc.Save();
224  }
225
226  socketstream sout;
227  if (visualization)
228  {
229  char vishost[] = "localhost";
230  int visport = 19916;
231  sout.open(vishost, visport);
232  if (!sout)
233  {
234  cout << "Unable to connect to GLVis server at "
235  << vishost << ':' << visport << endl;
236  visualization = false;
237  cout << "GLVis visualization disabled.\n";
238  }
239  else
240  {
241  sout.precision(precision);
242  sout << "solution\n" << *mesh << u;
243  sout << "pause\n";
244  sout << flush;
245  cout << "GLVis visualization paused."
246  << " Press space (in the GLVis window) to resume it.\n";
247  }
248  }
249
250  // 8. Define the time-dependent evolution operator describing the ODE
251  // right-hand side, and perform time-integration (looping over the time
252  // iterations, ti, with a time-step dt).
255
256  double t = 0.0;
257  for (int ti = 0; true; )
258  {
259  if (t >= t_final - dt/2)
260  {
261  break;
262  }
263
264  ode_solver->Step(u, t, dt);
265  ti++;
266
267  if (ti % vis_steps == 0)
268  {
269  cout << "time step: " << ti << ", time: " << t << endl;
270
271  if (visualization)
272  {
273  sout << "solution\n" << *mesh << u << flush;
274  }
275
276  if (visit)
277  {
278  visit_dc.SetCycle(ti);
279  visit_dc.SetTime(t);
280  visit_dc.Save();
281  }
282  }
283  }
284
285  // 9. Save the final solution. This output can be viewed later using GLVis:
286  // "glvis -m ex9.mesh -g ex9-final.gf".
287  {
288  ofstream osol("ex9-final.gf");
289  osol.precision(precision);
290  u.Save(osol);
291  }
292
293  // 10. Free the used memory.
294  delete ode_solver;
295  delete mesh;
296
297  return 0;
298 }
299
300
301 // Implementation of class FE_Evolution
302 FE_Evolution::FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b)
303  : TimeDependentOperator(_M.Size()), M(_M), K(_K), b(_b), z(_M.Size())
304 {
305  M_solver.SetPreconditioner(M_prec);
306  M_solver.SetOperator(M);
307
308  M_solver.iterative_mode = false;
309  M_solver.SetRelTol(1e-9);
310  M_solver.SetAbsTol(0.0);
311  M_solver.SetMaxIter(100);
312  M_solver.SetPrintLevel(0);
313 }
314
315 void FE_Evolution::Mult(const Vector &x, Vector &y) const
316 {
317  // y = M^{-1} (K x + b)
318  K.Mult(x, z);
319  z += b;
320  M_solver.Mult(z, y);
321 }
322
323
324 // Velocity coefficient
325 void velocity_function(const Vector &x, Vector &v)
326 {
327  int dim = x.Size();
328
329  switch (problem)
330  {
331  case 0:
332  {
333  // Translations in 1D, 2D, and 3D
334  switch (dim)
335  {
336  case 1: v(0) = 1.0; break;
337  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
338  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
339  break;
340  }
341  break;
342  }
343  case 1:
344  case 2:
345  {
346  // Clockwise rotation in 2D around the origin
347  const double w = M_PI/2;
348  switch (dim)
349  {
350  case 1: v(0) = 1.0; break;
351  case 2: v(0) = w*x(1); v(1) = -w*x(0); break;
352  case 3: v(0) = w*x(1); v(1) = -w*x(0); v(2) = 0.0; break;
353  }
354  break;
355  }
356  case 3:
357  {
358  // Clockwise twisting rotation in 2D around the origin
359  const double w = M_PI/2;
360  double d = max((x(0)+1.)*(1.-x(0)),0.) * max((x(1)+1.)*(1.-x(1)),0.);
361  d = d*d;
362  switch (dim)
363  {
364  case 1: v(0) = 1.0; break;
365  case 2: v(0) = d*w*x(1); v(1) = -d*w*x(0); break;
366  case 3: v(0) = d*w*x(1); v(1) = -d*w*x(0); v(2) = 0.0; break;
367  }
368  break;
369  }
370  }
371 }
372
373 // Initial condition
374 double u0_function(const Vector &x)
375 {
376  int dim = x.Size();
377
378  switch (problem)
379  {
380  case 0:
381  case 1:
382  {
383  switch (dim)
384  {
385  case 1:
386  return exp(-40.*pow(x(0)-0.5,2));
387  case 2:
388  case 3:
389  {
390  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
391  if (dim == 3)
392  {
393  const double s = (1. + 0.25*cos(2*M_PI*x(2)));
394  rx *= s;
395  ry *= s;
396  }
397  return ( erfc(w*(x(0)-cx-rx))*erfc(-w*(x(0)-cx+rx)) *
398  erfc(w*(x(1)-cy-ry))*erfc(-w*(x(1)-cy+ry)) )/16;
399  }
400  }
401  }
402  case 2:
403  {
404  const double r = sqrt(8.);
405  double x_ = x(0), y_ = x(1), rho, phi;
406  rho = hypot(x_, y_) / r;
407  phi = atan2(y_, x_);
408  return pow(sin(M_PI*rho),2)*sin(3*phi);
409  }
410  case 3:
411  {
412  const double f = M_PI;
413  return sin(f*x(0))*sin(f*x(1));
414  }
415  }
416  return 0.0;
417 }
418
419 // Inflow boundary condition (zero for the problems considered in this example)
420 double inflow_function(const Vector &x)
421 {
422  switch (problem)
423  {
424  case 0:
425  case 1:
426  case 2:
427  case 3: return 0.0;
428  }
429  return 0.0;
430 }
int GetVSize() const
Definition: fespace.hpp:164
Definition: solvers.hpp:109
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:445
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:84
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
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2138
int dim
Definition: ex3.cpp:48
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
Data collection with VisIt I/O routines.
Definition: linearform.cpp:29
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3719
int problem
Definition: ex9.cpp:37
int Dimension() const
Definition: mesh.hpp:475
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
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:8332
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:85
int main(int argc, char *argv[])
Definition: ex1.cpp:45
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:325
Abstract finite element space.
Definition: fespace.hpp:62
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:72
NURBSExtension * NURBSext
Definition: mesh.hpp:142
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:1140
int open(const char hostname[], int port)
double u0_function(const Vector &x)
Definition: ex9.cpp:374
class for C-function coefficient
Vector data type.
Definition: vector.hpp:33
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5844
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:54
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:420