MFEM  v4.2.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 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 // Device sample runs:
17 // ex9 -pa
18 // ex9 -ea
19 // ex9 -fa
20 // ex9 -pa -m ../../data/periodic-cube.mesh
21 // ex9 -pa -m ../../data/periodic-cube.mesh -d cuda
22 // ex9 -ea -m ../../data/periodic-cube.mesh -d cuda
23 // ex9 -fa -m ../../data/periodic-cube.mesh -d cuda
24 //
25 // Description: This example code solves the time-dependent advection equation
26 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
27 // u0(x)=u(0,x) is a given initial condition.
28 //
29 // The example demonstrates the use of Discontinuous Galerkin (DG)
30 // bilinear forms in MFEM (face integrators), the use of implicit
31 // and explicit ODE time integrators, the definition of periodic
32 // boundary conditions through periodic meshes, as well as the use
33 // of GLVis for persistent visualization of a time-evolving
34 // solution. The saving of time-dependent data files for external
35 // visualization with VisIt (visit.llnl.gov) and ParaView
36 // (paraview.org) is also illustrated.
37 
38 #include "mfem.hpp"
39 #include <fstream>
40 #include <iostream>
41 #include <algorithm>
42 
43 #ifndef MFEM_USE_SUNDIALS
44 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
45 #endif
46 
47 using namespace std;
48 using namespace mfem;
49 
50 // Choice for the problem setup. The fluid velocity, initial condition and
51 // inflow boundary condition are chosen based on this parameter.
52 int problem;
53 
54 // Velocity coefficient
55 void velocity_function(const Vector &x, Vector &v);
56 
57 // Initial condition
58 double u0_function(const Vector &x);
59 
60 // Inflow boundary condition
61 double inflow_function(const Vector &x);
62 
63 // Mesh bounding box
65 
66 class DG_Solver : public Solver
67 {
68 private:
69  SparseMatrix &M, &K, A;
70  GMRESSolver linear_solver;
71  BlockILU prec;
72  double dt;
73 public:
74  DG_Solver(SparseMatrix &M_, SparseMatrix &K_, const FiniteElementSpace &fes)
75  : M(M_),
76  K(K_),
77  prec(fes.GetFE(0)->GetDof(),
78  BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
79  dt(-1.0)
80  {
81  linear_solver.iterative_mode = false;
82  linear_solver.SetRelTol(1e-9);
83  linear_solver.SetAbsTol(0.0);
84  linear_solver.SetMaxIter(100);
85  linear_solver.SetPrintLevel(0);
86  linear_solver.SetPreconditioner(prec);
87  }
88 
89  void SetTimeStep(double dt_)
90  {
91  if (dt_ != dt)
92  {
93  dt = dt_;
94  // Form operator A = M - dt*K
95  A = K;
96  A *= -dt;
97  A += M;
98 
99  // this will also call SetOperator on the preconditioner
100  linear_solver.SetOperator(A);
101  }
102  }
103 
104  void SetOperator(const Operator &op)
105  {
106  linear_solver.SetOperator(op);
107  }
108 
109  virtual void Mult(const Vector &x, Vector &y) const
110  {
111  linear_solver.Mult(x, y);
112  }
113 };
114 
115 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
116  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
117  and advection matrices, and b describes the flow on the boundary. This can
118  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
119  used to evaluate the right-hand side. */
120 class FE_Evolution : public TimeDependentOperator
121 {
122 private:
123  BilinearForm &M, &K;
124  const Vector &b;
125  Solver *M_prec;
126  CGSolver M_solver;
127  DG_Solver *dg_solver;
128 
129  mutable Vector z;
130 
131 public:
132  FE_Evolution(BilinearForm &_M, BilinearForm &_K, const Vector &_b);
133 
134  virtual void Mult(const Vector &x, Vector &y) const;
135  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
136 
137  virtual ~FE_Evolution();
138 };
139 
140 
141 int main(int argc, char *argv[])
142 {
143  // 1. Parse command-line options.
144  problem = 0;
145  const char *mesh_file = "../../data/periodic-hexagon.mesh";
146  int ref_levels = 2;
147  int order = 3;
148  bool pa = false;
149  bool ea = false;
150  bool fa = false;
151  const char *device_config = "cpu";
152  int ode_solver_type = 7;
153  double t_final = 10.0;
154  double dt = 0.01;
155  bool visualization = false;
156  bool visit = false;
157  bool paraview = false;
158  bool binary = false;
159  int vis_steps = 5;
160 
161  // Relative and absolute tolerances for CVODE and ARKODE.
162  const double reltol = 1e-2, abstol = 1e-2;
163 
164  int precision = 8;
165  cout.precision(precision);
166 
167  OptionsParser args(argc, argv);
168  args.AddOption(&mesh_file, "-m", "--mesh",
169  "Mesh file to use.");
170  args.AddOption(&problem, "-p", "--problem",
171  "Problem setup to use. See options in velocity_function().");
172  args.AddOption(&ref_levels, "-r", "--refine",
173  "Number of times to refine the mesh uniformly.");
174  args.AddOption(&order, "-o", "--order",
175  "Order (degree) of the finite elements.");
176  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
177  "--no-partial-assembly", "Enable Partial Assembly.");
178  args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
179  "--no-element-assembly", "Enable Element Assembly.");
180  args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
181  "--no-full-assembly", "Enable Full Assembly.");
182  args.AddOption(&device_config, "-d", "--device",
183  "Device configuration string, see Device::Configure().");
184  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
185  "ODE solver:\n\t"
186  "1 - Forward Euler,\n\t"
187  "2 - RK2 SSP,\n\t"
188  "3 - RK3 SSP,\n\t"
189  "4 - RK4,\n\t"
190  "6 - RK6,\n\t"
191  "7 - CVODE (adaptive order implicit Adams),\n\t"
192  "8 - ARKODE default (4th order) explicit,\n\t"
193  "9 - ARKODE RK8.");
194  args.AddOption(&t_final, "-tf", "--t-final",
195  "Final time; start time is 0.");
196  args.AddOption(&dt, "-dt", "--time-step",
197  "Time step.");
198  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
199  "--no-visualization",
200  "Enable or disable GLVis visualization.");
201  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
202  "--no-visit-datafiles",
203  "Save data files for VisIt (visit.llnl.gov) visualization.");
204  args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
205  "--no-paraview-datafiles",
206  "Save data files for ParaView (paraview.org) visualization.");
207  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
208  "--ascii-datafiles",
209  "Use binary (Sidre) or ascii format for VisIt data files.");
210  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
211  "Visualize every n-th timestep.");
212  args.Parse();
213  if (!args.Good())
214  {
215  args.PrintUsage(cout);
216  return 1;
217  }
218  // check for valid ODE solver option
219  if (ode_solver_type < 1 || ode_solver_type > 9)
220  {
221  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
222  return 3;
223  }
224  args.PrintOptions(cout);
225 
226  Device device(device_config);
227  device.Print();
228 
229  // 2. Read the mesh from the given mesh file. We can handle geometrically
230  // periodic meshes in this code.
231  Mesh mesh(mesh_file, 1, 1);
232  int dim = mesh.Dimension();
233 
234  // 3. Refine the mesh to increase the resolution. In this example we do
235  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
236  // command-line parameter. If the mesh is of NURBS type, we convert it to
237  // a (piecewise-polynomial) high-order mesh.
238  for (int lev = 0; lev < ref_levels; lev++)
239  {
240  mesh.UniformRefinement();
241  }
242  if (mesh.NURBSext)
243  {
244  mesh.SetCurvature(max(order, 1));
245  }
246  mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
247 
248  // 4. Define the discontinuous DG finite element space of the given
249  // polynomial order on the refined mesh.
250  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
251  FiniteElementSpace fes(&mesh, &fec);
252 
253  cout << "Number of unknowns: " << fes.GetVSize() << endl;
254 
255  // 5. Set up and assemble the bilinear and linear forms corresponding to the
256  // DG discretization. The DGTraceIntegrator involves integrals over mesh
257  // interior faces.
261 
262  BilinearForm m(&fes);
263  BilinearForm k(&fes);
264  if (pa)
265  {
266  m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
267  k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
268  }
269  else if (ea)
270  {
271  m.SetAssemblyLevel(AssemblyLevel::ELEMENT);
272  k.SetAssemblyLevel(AssemblyLevel::ELEMENT);
273  }
274  else if (fa)
275  {
276  m.SetAssemblyLevel(AssemblyLevel::FULL);
277  k.SetAssemblyLevel(AssemblyLevel::FULL);
278  }
280  k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
282  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
284  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
285 
286  LinearForm b(&fes);
288  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
289 
290  m.Assemble();
291  int skip_zeros = 0;
292  k.Assemble(skip_zeros);
293  b.Assemble();
294  m.Finalize();
295  k.Finalize(skip_zeros);
296 
297  // 6. Define the initial conditions, save the corresponding grid function to
298  // a file and (optionally) save data in the VisIt format and initialize
299  // GLVis visualization.
300  GridFunction u(&fes);
301  u.ProjectCoefficient(u0);
302 
303  {
304  ofstream omesh("ex9.mesh");
305  omesh.precision(precision);
306  mesh.Print(omesh);
307  ofstream osol("ex9-init.gf");
308  osol.precision(precision);
309  u.Save(osol);
310  }
311 
312  // Create data collection for solution output: either VisItDataCollection for
313  // ascii data files, or SidreDataCollection for binary data files.
314  DataCollection *dc = NULL;
315  if (visit)
316  {
317  if (binary)
318  {
319 #ifdef MFEM_USE_SIDRE
320  dc = new SidreDataCollection("Example9", &mesh);
321 #else
322  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
323 #endif
324  }
325  else
326  {
327  dc = new VisItDataCollection("Example9", &mesh);
328  dc->SetPrecision(precision);
329  }
330  dc->RegisterField("solution", &u);
331  dc->SetCycle(0);
332  dc->SetTime(0.0);
333  dc->Save();
334  }
335 
336  ParaViewDataCollection *pd = NULL;
337  if (paraview)
338  {
339  pd = new ParaViewDataCollection("Example9", &mesh);
340  pd->SetPrefixPath("ParaView");
341  pd->RegisterField("solution", &u);
342  pd->SetLevelsOfDetail(order);
343  pd->SetDataFormat(VTKFormat::BINARY);
344  pd->SetHighOrderOutput(true);
345  pd->SetCycle(0);
346  pd->SetTime(0.0);
347  pd->Save();
348  }
349 
350  socketstream sout;
351  if (visualization)
352  {
353  char vishost[] = "localhost";
354  int visport = 19916;
355  sout.open(vishost, visport);
356  if (!sout)
357  {
358  cout << "Unable to connect to GLVis server at "
359  << vishost << ':' << visport << endl;
360  visualization = false;
361  cout << "GLVis visualization disabled.\n";
362  }
363  else
364  {
365  sout.precision(precision);
366  sout << "solution\n" << mesh << u;
367  sout << "pause\n";
368  sout << flush;
369  cout << "GLVis visualization paused."
370  << " Press space (in the GLVis window) to resume it.\n";
371  }
372  }
373 
374  // 7. Define the time-dependent evolution operator describing the ODE
375  // right-hand side, and define the ODE solver used for time integration.
376  FE_Evolution adv(m, k, b);
377 
378  double t = 0.0;
379  adv.SetTime(t);
380 
381  // Create the time integrator
382  ODESolver *ode_solver = NULL;
383  CVODESolver *cvode = NULL;
384  ARKStepSolver *arkode = NULL;
385  switch (ode_solver_type)
386  {
387  case 1: ode_solver = new ForwardEulerSolver; break;
388  case 2: ode_solver = new RK2Solver(1.0); break;
389  case 3: ode_solver = new RK3SSPSolver; break;
390  case 4: ode_solver = new RK4Solver; break;
391  case 6: ode_solver = new RK6Solver; break;
392  case 7:
393  cvode = new CVODESolver(CV_ADAMS);
394  cvode->Init(adv);
395  cvode->SetSStolerances(reltol, abstol);
396  cvode->SetMaxStep(dt);
397  cvode->UseSundialsLinearSolver();
398  ode_solver = cvode; break;
399  case 8:
400  arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
401  arkode->Init(adv);
402  arkode->SetSStolerances(reltol, abstol);
403  arkode->SetMaxStep(dt);
404  arkode->SetOrder(4);
405  ode_solver = arkode; break;
406  case 9:
407  arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
408  arkode->Init(adv);
409  arkode->SetSStolerances(reltol, abstol);
410  arkode->SetMaxStep(dt);
411  arkode->SetERKTableNum(FEHLBERG_13_7_8);
412  ode_solver = arkode; break;
413  }
414 
415  // Initialize MFEM integrators, SUNDIALS integrators are initialized above
416  if (ode_solver_type < 7) { ode_solver->Init(adv); }
417 
418  // 8. Perform time-integration (looping over the time iterations, ti,
419  // with a time-step dt).
420  bool done = false;
421  for (int ti = 0; !done; )
422  {
423  double dt_real = min(dt, t_final - t);
424  ode_solver->Step(u, t, dt_real);
425  ti++;
426 
427  done = (t >= t_final - 1e-8*dt);
428 
429  if (done || ti % vis_steps == 0)
430  {
431  cout << "time step: " << ti << ", time: " << t << endl;
432  if (cvode) { cvode->PrintInfo(); }
433  if (arkode) { arkode->PrintInfo(); }
434 
435  if (visualization)
436  {
437  sout << "solution\n" << mesh << u << flush;
438  }
439 
440  if (visit)
441  {
442  dc->SetCycle(ti);
443  dc->SetTime(t);
444  dc->Save();
445  }
446 
447  if (paraview)
448  {
449  pd->SetCycle(ti);
450  pd->SetTime(t);
451  pd->Save();
452  }
453  }
454  }
455 
456  // 9. Save the final solution. This output can be viewed later using GLVis:
457  // "glvis -m ex9.mesh -g ex9-final.gf".
458  {
459  ofstream osol("ex9-final.gf");
460  osol.precision(precision);
461  u.Save(osol);
462  }
463 
464  // 10. Free the used memory.
465  delete ode_solver;
466  delete pd;
467  delete dc;
468 
469  return 0;
470 }
471 
472 
473 // Implementation of class FE_Evolution
475  : TimeDependentOperator(_M.Height()), M(_M), K(_K), b(_b), z(_M.Height())
476 {
477  Array<int> ess_tdof_list;
478  if (M.GetAssemblyLevel() == AssemblyLevel::LEGACYFULL)
479  {
480  M_prec = new DSmoother(M.SpMat());
481  M_solver.SetOperator(M.SpMat());
482  dg_solver = new DG_Solver(M.SpMat(), K.SpMat(), *M.FESpace());
483  }
484  else
485  {
486  M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
487  M_solver.SetOperator(M);
488  dg_solver = NULL;
489  }
490  M_solver.SetPreconditioner(*M_prec);
491  M_solver.iterative_mode = false;
492  M_solver.SetRelTol(1e-9);
493  M_solver.SetAbsTol(0.0);
494  M_solver.SetMaxIter(100);
495  M_solver.SetPrintLevel(0);
496 }
497 
498 void FE_Evolution::Mult(const Vector &x, Vector &y) const
499 {
500  // y = M^{-1} (K x + b)
501  K.Mult(x, z);
502  z += b;
503  M_solver.Mult(z, y);
504 }
505 
506 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
507 {
508  MFEM_VERIFY(dg_solver != NULL,
509  "Implicit time integration is not supported with partial assembly");
510  K.Mult(x, z);
511  z += b;
512  dg_solver->SetTimeStep(dt);
513  dg_solver->Mult(z, k);
514 }
515 
517 {
518  delete M_prec;
519  delete dg_solver;
520 }
521 
522 // Velocity coefficient
523 void velocity_function(const Vector &x, Vector &v)
524 {
525  int dim = x.Size();
526 
527  // map to the reference [-1,1] domain
528  Vector X(dim);
529  for (int i = 0; i < dim; i++)
530  {
531  double center = (bb_min[i] + bb_max[i]) * 0.5;
532  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
533  }
534 
535  switch (problem)
536  {
537  case 0:
538  {
539  // Translations in 1D, 2D, and 3D
540  switch (dim)
541  {
542  case 1: v(0) = 1.0; break;
543  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
544  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
545  break;
546  }
547  break;
548  }
549  case 1:
550  case 2:
551  {
552  // Clockwise rotation in 2D around the origin
553  const double w = M_PI/2;
554  switch (dim)
555  {
556  case 1: v(0) = 1.0; break;
557  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
558  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
559  }
560  break;
561  }
562  case 3:
563  {
564  // Clockwise twisting rotation in 2D around the origin
565  const double w = M_PI/2;
566  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
567  d = d*d;
568  switch (dim)
569  {
570  case 1: v(0) = 1.0; break;
571  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
572  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
573  }
574  break;
575  }
576  }
577 }
578 
579 // Initial condition
580 double u0_function(const Vector &x)
581 {
582  int dim = x.Size();
583 
584  // map to the reference [-1,1] domain
585  Vector X(dim);
586  for (int i = 0; i < dim; i++)
587  {
588  double center = (bb_min[i] + bb_max[i]) * 0.5;
589  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
590  }
591 
592  switch (problem)
593  {
594  case 0:
595  case 1:
596  {
597  switch (dim)
598  {
599  case 1:
600  return exp(-40.*pow(X(0)-0.5,2));
601  case 2:
602  case 3:
603  {
604  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
605  if (dim == 3)
606  {
607  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
608  rx *= s;
609  ry *= s;
610  }
611  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
612  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
613  }
614  }
615  }
616  case 2:
617  {
618  double x_ = X(0), y_ = X(1), rho, phi;
619  rho = hypot(x_, y_);
620  phi = atan2(y_, x_);
621  return pow(sin(M_PI*rho),2)*sin(3*phi);
622  }
623  case 3:
624  {
625  const double f = M_PI;
626  return sin(f*X(0))*sin(f*X(1));
627  }
628  }
629  return 0.0;
630 }
631 
632 // Inflow boundary condition (zero for the problems considered in this example)
633 double inflow_function(const Vector &x)
634 {
635  switch (problem)
636  {
637  case 0:
638  case 1:
639  case 2:
640  case 3: return 0.0;
641  }
642  return 0.0;
643 }
Vector bb_max
Definition: ex9.cpp:64
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:532
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:400
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1234
Conjugate gradient method.
Definition: solvers.hpp:258
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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:698
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void SetDataFormat(VTKFormat fmt)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:535
Helper class for ParaView visualization data.
Base abstract class for first order time dependent operators.
Definition: operator.hpp:267
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 SetOperator(const Operator &a)
Set/update the solver for the given operator.
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:311
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
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:1540
void SetOrder(int order)
Chooses integration order for all explicit / implicit / IMEX methods.
Definition: sundials.cpp:1546
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:261
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
int main(int argc, char *argv[])
Definition: ex1.cpp:66
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:102
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:539
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:3417
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
virtual void Save()
Save the collection to disk.
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:252
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:128
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4165
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:716
void SetHighOrderOutput(bool high_order_output_)
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:64
A general vector function coefficient.
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition: ex9.cpp:480
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:497
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:315
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int problem
Definition: ex15.cpp:54
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
GMRES method.
Definition: solvers.hpp:290
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:128
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:203
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:1576
int dim
Definition: ex24.cpp:53
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.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:1798
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:2252
virtual ~FE_Evolution()
Definition: ex18.hpp:43
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
double u0_function(const Vector &x)
Definition: ex9.cpp:554
void SetLevelsOfDetail(int levels_of_detail_)
A general function coefficient.
Vector data type.
Definition: vector.hpp:51
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
virtual void Save() override
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:742
Base class for solvers.
Definition: operator.hpp:634
void SetERKTableNum(int table_num)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:1552
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:118
The classical forward Euler method.
Definition: ode.hpp:116
Abstract operator.
Definition: operator.hpp:24
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:1297
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:1534
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
double inflow_function(const Vector &x)
Definition: ex9.cpp:607
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:678
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:221
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145
alpha (q . grad u, v)