MFEM  v4.1.0 Finite element discretization library
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 7 -dt 0.005
8 // ex9 -m ../../data/periodic-square.mesh -p 1 -r 2 -s 8 -dt 0.005 -tf 9
9 // ex9 -m ../../data/periodic-hexagon.mesh -p 0 -r 2 -s 7 -dt 0.0018 -vs 25
10 // ex9 -m ../../data/periodic-hexagon.mesh -p 0 -r 2 -s 9 -dt 0.01 -vs 15
11 // ex9 -m ../../data/amr-quad.mesh -p 1 -r 2 -s 9 -dt 0.002 -tf 9
12 // ex9 -m ../../data/star-q3.mesh -p 1 -r 2 -s 9 -dt 0.005 -tf 9
13 // ex9 -m ../../data/disc-nurbs.mesh -p 1 -r 3 -s 7 -dt 0.005 -tf 9
14 // ex9 -m ../../data/periodic-cube.mesh -p 0 -r 2 -s 8 -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);
104  "Mesh file to use.");
106  "Problem setup to use. See options in velocity_function().");
108  "Number of times to refine the mesh uniformly.");
110  "Order (degree) of the finite elements.");
112  "ODE solver:\n\t"
113  "1 - Forward Euler,\n\t"
114  "2 - RK2 SSP,\n\t"
115  "3 - RK3 SSP,\n\t"
116  "4 - RK4,\n\t"
117  "6 - RK6,\n\t"
119  "8 - ARKODE default (4th order) explicit,\n\t"
120  "9 - ARKODE RK8.");
122  "Final time; start time is 0.");
124  "Time step.");
126  "--no-visualization",
127  "Enable or disable GLVis visualization.");
129  "--no-visit-datafiles",
130  "Save data files for VisIt (visit.llnl.gov) visualization.");
132  "--ascii-datafiles",
133  "Use binary (Sidre) or ascii format for VisIt data files.");
135  "Visualize every n-th timestep.");
136  args.Parse();
137  if (!args.Good())
138  {
139  args.PrintUsage(cout);
140  return 1;
141  }
142  // check for vaild ODE solver option
143  if (ode_solver_type < 1 || ode_solver_type > 9)
144  {
145  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
146  return 3;
147  }
148  args.PrintOptions(cout);
149
150  // 2. Read the mesh from the given mesh file. We can handle geometrically
151  // periodic meshes in this code.
152  Mesh mesh(mesh_file, 1, 1);
153  int dim = mesh.Dimension();
154
155  // 3. Refine the mesh to increase the resolution. In this example we do
156  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
157  // command-line parameter. If the mesh is of NURBS type, we convert it to
158  // a (piecewise-polynomial) high-order mesh.
159  for (int lev = 0; lev < ref_levels; lev++)
160  {
161  mesh.UniformRefinement();
162  }
163  if (mesh.NURBSext)
164  {
165  mesh.SetCurvature(max(order, 1));
166  }
167  mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
168
169  // 4. Define the discontinuous DG finite element space of the given
170  // polynomial order on the refined mesh.
171  DG_FECollection fec(order, dim);
172  FiniteElementSpace fes(&mesh, &fec);
173
174  cout << "Number of unknowns: " << fes.GetVSize() << endl;
175
176  // 5. Set up and assemble the bilinear and linear forms corresponding to the
177  // DG discretization. The DGTraceIntegrator involves integrals over mesh
178  // interior faces.
182
183  BilinearForm m(&fes);
185  BilinearForm k(&fes);
188  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
190  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
191
192  LinearForm b(&fes);
194  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
195
196  m.Assemble();
197  m.Finalize();
198  int skip_zeros = 0;
199  k.Assemble(skip_zeros);
200  k.Finalize(skip_zeros);
201  b.Assemble();
202
203  // 6. Define the initial conditions, save the corresponding grid function to
204  // a file and (optionally) save data in the VisIt format and initialize
205  // GLVis visualization.
206  GridFunction u(&fes);
207  u.ProjectCoefficient(u0);
208
209  {
210  ofstream omesh("ex9.mesh");
211  omesh.precision(precision);
212  mesh.Print(omesh);
213  ofstream osol("ex9-init.gf");
214  osol.precision(precision);
215  u.Save(osol);
216  }
217
218  // Create data collection for solution output: either VisItDataCollection for
219  // ascii data files, or SidreDataCollection for binary data files.
220  DataCollection *dc = NULL;
221  if (visit)
222  {
223  if (binary)
224  {
225 #ifdef MFEM_USE_SIDRE
226  dc = new SidreDataCollection("Example9", &mesh);
227 #else
228  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
229 #endif
230  }
231  else
232  {
233  dc = new VisItDataCollection("Example9", &mesh);
234  dc->SetPrecision(precision);
235  }
236  dc->RegisterField("solution", &u);
237  dc->SetCycle(0);
238  dc->SetTime(0.0);
239  dc->Save();
240  }
241
242  socketstream sout;
243  if (visualization)
244  {
245  char vishost[] = "localhost";
246  int visport = 19916;
247  sout.open(vishost, visport);
248  if (!sout)
249  {
250  cout << "Unable to connect to GLVis server at "
251  << vishost << ':' << visport << endl;
252  visualization = false;
253  cout << "GLVis visualization disabled.\n";
254  }
255  else
256  {
257  sout.precision(precision);
258  sout << "solution\n" << mesh << u;
259  sout << "pause\n";
260  sout << flush;
261  cout << "GLVis visualization paused."
262  << " Press space (in the GLVis window) to resume it.\n";
263  }
264  }
265
266  // 7. Define the time-dependent evolution operator describing the ODE
267  // right-hand side, and define the ODE solver used for time integration.
269
270  double t = 0.0;
272
273  // Create the time integrator
274  ODESolver *ode_solver = NULL;
275  CVODESolver *cvode = NULL;
276  ARKStepSolver *arkode = NULL;
277  switch (ode_solver_type)
278  {
279  case 1: ode_solver = new ForwardEulerSolver; break;
280  case 2: ode_solver = new RK2Solver(1.0); break;
281  case 3: ode_solver = new RK3SSPSolver; break;
282  case 4: ode_solver = new RK4Solver; break;
283  case 6: ode_solver = new RK6Solver; break;
284  case 7:
287  cvode->SetSStolerances(reltol, abstol);
288  cvode->SetMaxStep(dt);
289  cvode->UseSundialsLinearSolver();
290  ode_solver = cvode; break;
291  case 8:
292  case 9:
293  arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
295  arkode->SetSStolerances(reltol, abstol);
296  arkode->SetMaxStep(dt);
297  if (ode_solver_type == 9) { arkode->SetERKTableNum(FEHLBERG_13_7_8); }
298  ode_solver = arkode; break;
299  }
300
301  // Initialize MFEM integrators, SUNDIALS integrators are initialized above
302  if (ode_solver_type < 7) { ode_solver->Init(adv); }
303
304  // 8. Perform time-integration (looping over the time iterations, ti,
305  // with a time-step dt).
306  bool done = false;
307  for (int ti = 0; !done; )
308  {
309  double dt_real = min(dt, t_final - t);
310  ode_solver->Step(u, t, dt_real);
311  ti++;
312
313  done = (t >= t_final - 1e-8*dt);
314
315  if (done || ti % vis_steps == 0)
316  {
317  cout << "time step: " << ti << ", time: " << t << endl;
318  if (cvode) { cvode->PrintInfo(); }
319  if (arkode) { arkode->PrintInfo(); }
320
321  if (visualization)
322  {
323  sout << "solution\n" << mesh << u << flush;
324  }
325
326  if (visit)
327  {
328  dc->SetCycle(ti);
329  dc->SetTime(t);
330  dc->Save();
331  }
332  }
333  }
334
335  // 9. Save the final solution. This output can be viewed later using GLVis:
336  // "glvis -m ex9.mesh -g ex9-final.gf".
337  {
338  ofstream osol("ex9-final.gf");
339  osol.precision(precision);
340  u.Save(osol);
341  }
342
343  // 10. Free the used memory.
344  delete ode_solver;
345  delete dc;
346
347  return 0;
348 }
349
350
351 // Implementation of class FE_Evolution
353  : TimeDependentOperator(_M.Size()), M(_M), K(_K), b(_b), z(_M.Size())
354 {
355  M_solver.SetPreconditioner(M_prec);
356  M_solver.SetOperator(M);
357
358  M_solver.iterative_mode = false;
359  M_solver.SetRelTol(1e-9);
360  M_solver.SetAbsTol(0.0);
361  M_solver.SetMaxIter(100);
362  M_solver.SetPrintLevel(0);
363 }
364
365 void FE_Evolution::Mult(const Vector &x, Vector &y) const
366 {
367  // y = M^{-1} (K x + b)
368  K.Mult(x, z);
369  z += b;
370  M_solver.Mult(z, y);
371 }
372
373
374 // Velocity coefficient
375 void velocity_function(const Vector &x, Vector &v)
376 {
377  int dim = x.Size();
378
379  // map to the reference [-1,1] domain
380  Vector X(dim);
381  for (int i = 0; i < dim; i++)
382  {
383  double center = (bb_min[i] + bb_max[i]) * 0.5;
384  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
385  }
386
387  switch (problem)
388  {
389  case 0:
390  {
391  // Translations in 1D, 2D, and 3D
392  switch (dim)
393  {
394  case 1: v(0) = 1.0; break;
395  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
396  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
397  break;
398  }
399  break;
400  }
401  case 1:
402  case 2:
403  {
404  // Clockwise rotation in 2D around the origin
405  const double w = M_PI/2;
406  switch (dim)
407  {
408  case 1: v(0) = 1.0; break;
409  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
410  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
411  }
412  break;
413  }
414  case 3:
415  {
416  // Clockwise twisting rotation in 2D around the origin
417  const double w = M_PI/2;
418  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
419  d = d*d;
420  switch (dim)
421  {
422  case 1: v(0) = 1.0; break;
423  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
424  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
425  }
426  break;
427  }
428  }
429 }
430
431 // Initial condition
432 double u0_function(const Vector &x)
433 {
434  int dim = x.Size();
435
436  // map to the reference [-1,1] domain
437  Vector X(dim);
438  for (int i = 0; i < dim; i++)
439  {
440  double center = (bb_min[i] + bb_max[i]) * 0.5;
441  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
442  }
443
444  switch (problem)
445  {
446  case 0:
447  case 1:
448  {
449  switch (dim)
450  {
451  case 1:
452  return exp(-40.*pow(X(0)-0.5,2));
453  case 2:
454  case 3:
455  {
456  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
457  if (dim == 3)
458  {
459  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
460  rx *= s;
461  ry *= s;
462  }
463  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
464  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
465  }
466  }
467  }
468  case 2:
469  {
470  double x_ = X(0), y_ = X(1), rho, phi;
471  rho = hypot(x_, y_);
472  phi = atan2(y_, x_);
473  return pow(sin(M_PI*rho),2)*sin(3*phi);
474  }
475  case 3:
476  {
477  const double f = M_PI;
478  return sin(f*X(0))*sin(f*X(1));
479  }
480  }
481  return 0.0;
482 }
483
484 // Inflow boundary condition (zero for the problems considered in this example)
485 double inflow_function(const Vector &x)
486 {
487  switch (problem)
488  {
489  case 0:
490  case 1:
491  case 2:
492  case 3: return 0.0;
493  }
494  return 0.0;
495 }
Vector bb_max
Definition: ex9.cpp:60
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
void Init(TimeDependentOperator &f_)
Initialize CVODE: calls CVodeCreate() to create the CVODE memory and set some defaults.
Definition: sundials.cpp:151
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
Definition: solvers.hpp:152
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
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:342
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:358
Base abstract class for first order time dependent operators.
Definition: operator.hpp:259
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:120
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:303
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:835
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:504
int main(int argc, char *argv[])
Definition: ex1.cpp:62
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:104
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:197
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:2731
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
Data type sparse matrix.
Definition: sparsemat.hpp:40
virtual void Save()
Save the collection to disk.
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:92
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:65
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:348
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
int Dimension() const
Definition: mesh.hpp:744
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:60
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:145
void SetAbsTol(double atol)
Definition: solvers.hpp:64
void SetRelTol(double rtol)
Definition: solvers.hpp:63
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:478
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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:54
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:191
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:871
int dim
Definition: ex24.cpp:43
Adds new interior Face Integrator. Assumes ownership of 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:1685
int open(const char hostname[], int port)
double u0_function(const Vector &x)
Definition: ex9.cpp:535
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:166
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:91
Class for linear form - Vector with associated FE space and LFIntegrators.
Definition: linearform.hpp:23
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:360
const SparseMatrix & SpMat() const
Returns a reference to the sparse matrix.
void SetERKTableNum(int table_num)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:847
The classical forward Euler method.
Definition: ode.hpp:99
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:565
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:829
double inflow_function(const Vector &x)
Definition: ex9.cpp:588
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:322