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