MFEM  v4.5.2
Finite element discretization library
ex9p.cpp
Go to the documentation of this file.
1 // MFEM Example 9 - Parallel Version
2 // SUNDIALS Modification
3 //
4 // Compile with: make ex9p
5 //
6 // Sample runs:
7 // mpirun -np 4 ex9p -m ../../data/periodic-segment.mesh -p 1 -rp 1 -s 7 -dt 0.0025
8 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 1 -rp 1 -s 8 -dt 0.0025 -tf 9
9 // mpirun -np 4 ex9p -m ../../data/periodic-hexagon.mesh -p 0 -rp 1 -s 7 -dt 0.0009 -vs 25
10 // mpirun -np 4 ex9p -m ../../data/periodic-hexagon.mesh -p 0 -rp 1 -s 9 -dt 0.005 -vs 15
11 // mpirun -np 4 ex9p -m ../../data/amr-quad.mesh -p 1 -rp 1 -s 9 -dt 0.001 -tf 9
12 // mpirun -np 4 ex9p -m ../../data/star-q3.mesh -p 1 -rp 1 -s 9 -dt 0.0025 -tf 9
13 // mpirun -np 4 ex9p -m ../../data/disc-nurbs.mesh -p 1 -rp 2 -s 7 -dt 0.0025 -tf 9
14 // mpirun -np 4 ex9p -m ../../data/periodic-cube.mesh -p 0 -rp 1 -s 8 -dt 0.01 -tf 8 -o 2
15 //
16 // Device sample runs:
17 // mpirun -np 4 ex9p -pa
18 // mpirun -np 4 ex9p -ea
19 // mpirun -np 4 ex9p -fa
20 // mpirun -np 4 ex9p -pa -m ../../data/periodic-cube.mesh
21 // mpirun -np 4 ex9p -pa -m ../../data/periodic-cube.mesh -d cuda
22 // mpirun -np 4 ex9p -ea -m ../../data/periodic-cube.mesh -d cuda
23 // mpirun -np 4 ex9p -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. Saving of time-dependent data files for visualization
35 // with VisIt (visit.llnl.gov) and ParaView (paraview.org), as
36 // well as the optional saving with ADIOS2 (adios2.readthedocs.io)
37 // are also illustrated.
38 
39 #include "mfem.hpp"
40 #include <fstream>
41 #include <iostream>
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  HypreParMatrix &M, &K;
70  SparseMatrix M_diag;
71  HypreParMatrix *A;
72  GMRESSolver linear_solver;
73  BlockILU prec;
74  double dt;
75 public:
76  DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
77  : M(M_),
78  K(K_),
79  A(NULL),
80  linear_solver(M.GetComm()),
81  prec(fes.GetFE(0)->GetDof(),
82  BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
83  dt(-1.0)
84  {
85  linear_solver.iterative_mode = false;
86  linear_solver.SetRelTol(1e-9);
87  linear_solver.SetAbsTol(0.0);
88  linear_solver.SetMaxIter(100);
89  linear_solver.SetPrintLevel(0);
90  linear_solver.SetPreconditioner(prec);
91 
92  M.GetDiag(M_diag);
93  }
94 
95  void SetTimeStep(double dt_)
96  {
97  if (dt_ != dt)
98  {
99  dt = dt_;
100  // Form operator A = M - dt*K
101  delete A;
102  A = Add(-dt, K, 0.0, K);
103  SparseMatrix A_diag;
104  A->GetDiag(A_diag);
105  A_diag.Add(1.0, M_diag);
106  // this will also call SetOperator on the preconditioner
107  linear_solver.SetOperator(*A);
108  }
109  }
110 
111  void SetOperator(const Operator &op)
112  {
113  linear_solver.SetOperator(op);
114  }
115 
116  virtual void Mult(const Vector &x, Vector &y) const
117  {
118  linear_solver.Mult(x, y);
119  }
120 
121  ~DG_Solver()
122  {
123  delete A;
124  }
125 };
126 
127 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
128  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
129  and advection matrices, and b describes the flow on the boundary. This can
130  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
131  used to evaluate the right-hand side. */
132 class FE_Evolution : public TimeDependentOperator
133 {
134 private:
135  OperatorHandle M, K;
136  const Vector &b;
137  Solver *M_prec;
138  CGSolver M_solver;
139  DG_Solver *dg_solver;
140 
141  mutable Vector z;
142 
143 public:
144  FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, const Vector &b_);
145 
146  virtual void Mult(const Vector &x, Vector &y) const;
147  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
148 
149  virtual ~FE_Evolution();
150 };
151 
152 
153 int main(int argc, char *argv[])
154 {
155  // 1. Initialize MPI, HYPRE, and SUNDIALS.
156  Mpi::Init(argc, argv);
157  int num_procs = Mpi::WorldSize();
158  int myid = Mpi::WorldRank();
159  Hypre::Init();
160  Sundials::Init();
161 
162  // 2. Parse command-line options.
163  problem = 0;
164  const char *mesh_file = "../../data/periodic-hexagon.mesh";
165  int ser_ref_levels = 2;
166  int par_ref_levels = 0;
167  int order = 3;
168  bool pa = false;
169  bool ea = false;
170  bool fa = false;
171  const char *device_config = "cpu";
172  int ode_solver_type = 7;
173  double t_final = 10.0;
174  double dt = 0.01;
175  bool visualization = false;
176  bool visit = false;
177  bool paraview = false;
178  bool adios2 = false;
179  bool binary = false;
180  int vis_steps = 5;
181 
182  // Relative and absolute tolerances for CVODE and ARKODE.
183  const double reltol = 1e-2, abstol = 1e-2;
184 
185  int precision = 8;
186  cout.precision(precision);
187 
188  OptionsParser args(argc, argv);
189  args.AddOption(&mesh_file, "-m", "--mesh",
190  "Mesh file to use.");
191  args.AddOption(&problem, "-p", "--problem",
192  "Problem setup to use. See options in velocity_function().");
193  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
194  "Number of times to refine the mesh uniformly in serial.");
195  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
196  "Number of times to refine the mesh uniformly in parallel.");
197  args.AddOption(&order, "-o", "--order",
198  "Order (degree) of the finite elements.");
199  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
200  "--no-partial-assembly", "Enable Partial Assembly.");
201  args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
202  "--no-element-assembly", "Enable Element Assembly.");
203  args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
204  "--no-full-assembly", "Enable Full Assembly.");
205  args.AddOption(&device_config, "-d", "--device",
206  "Device configuration string, see Device::Configure().");
207  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
208  "ODE solver:\n\t"
209  "1 - Forward Euler,\n\t"
210  "2 - RK2 SSP,\n\t"
211  "3 - RK3 SSP,\n\t"
212  "4 - RK4,\n\t"
213  "6 - RK6,\n\t"
214  "7 - CVODE (adaptive order implicit Adams),\n\t"
215  "8 - ARKODE default (4th order) explicit,\n\t"
216  "9 - ARKODE RK8.");
217  args.AddOption(&t_final, "-tf", "--t-final",
218  "Final time; start time is 0.");
219  args.AddOption(&dt, "-dt", "--time-step",
220  "Time step.");
221  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
222  "--no-visualization",
223  "Enable or disable GLVis visualization.");
224  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
225  "--no-visit-datafiles",
226  "Save data files for VisIt (visit.llnl.gov) visualization.");
227  args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
228  "--no-paraview-datafiles",
229  "Save data files for ParaView (paraview.org) visualization.");
230  args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
231  "--no-adios2-streams",
232  "Save data using adios2 streams.");
233  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
234  "--ascii-datafiles",
235  "Use binary (Sidre) or ascii format for VisIt data files.");
236  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
237  "Visualize every n-th timestep.");
238  args.Parse();
239  if (!args.Good())
240  {
241  if (myid == 0)
242  {
243  args.PrintUsage(cout);
244  }
245  return 1;
246  }
247  if (myid == 0)
248  {
249  args.PrintOptions(cout);
250  }
251 
252  // check for valid ODE solver option
253  if (ode_solver_type < 1 || ode_solver_type > 9)
254  {
255  if (myid == 0)
256  {
257  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
258  }
259  return 3;
260  }
261 
262  Device device(device_config);
263  if (myid == 0) { device.Print(); }
264 
265  // 3. Read the serial mesh from the given mesh file on all processors. We can
266  // handle geometrically periodic meshes in this code.
267  Mesh *mesh = new Mesh(mesh_file, 1, 1);
268  int dim = mesh->Dimension();
269 
270  // 4. Refine the mesh in serial to increase the resolution. In this example
271  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
272  // a command-line parameter. If the mesh is of NURBS type, we convert it
273  // to a (piecewise-polynomial) high-order mesh.
274  for (int lev = 0; lev < ser_ref_levels; lev++)
275  {
276  mesh->UniformRefinement();
277  }
278  if (mesh->NURBSext)
279  {
280  mesh->SetCurvature(max(order, 1));
281  }
282  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
283 
284  // 5. Define the parallel mesh by a partitioning of the serial mesh. Refine
285  // this mesh further in parallel to increase the resolution. Once the
286  // parallel mesh is defined, the serial mesh can be deleted.
287  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
288  delete mesh;
289  for (int lev = 0; lev < par_ref_levels; lev++)
290  {
291  pmesh->UniformRefinement();
292  }
293 
294  // 6. Define the parallel discontinuous DG finite element space on the
295  // parallel refined mesh of the given polynomial order.
296  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
297  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
298 
299  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
300  if (myid == 0)
301  {
302  cout << "Number of unknowns: " << global_vSize << endl;
303  }
304 
305  // 7. Set up and assemble the parallel bilinear and linear forms (and the
306  // parallel hypre matrices) corresponding to the DG discretization. The
307  // DGTraceIntegrator involves integrals over mesh interior faces.
311 
312  ParBilinearForm *m = new ParBilinearForm(fes);
313  ParBilinearForm *k = new ParBilinearForm(fes);
314  if (pa)
315  {
316  m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
317  k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
318  }
319  else if (ea)
320  {
321  m->SetAssemblyLevel(AssemblyLevel::ELEMENT);
322  k->SetAssemblyLevel(AssemblyLevel::ELEMENT);
323  }
324  else if (fa)
325  {
326  m->SetAssemblyLevel(AssemblyLevel::FULL);
327  k->SetAssemblyLevel(AssemblyLevel::FULL);
328  }
329 
331  k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
333  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
335  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
336 
337  ParLinearForm *b = new ParLinearForm(fes);
338  b->AddBdrFaceIntegrator(
339  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
340 
341  int skip_zeros = 0;
342  m->Assemble();
343  k->Assemble(skip_zeros);
344  b->Assemble();
345  m->Finalize();
346  k->Finalize(skip_zeros);
347 
348  HypreParVector *B = b->ParallelAssemble();
349 
350  // 8. Define the initial conditions, save the corresponding grid function to
351  // a file and (optionally) save data in the VisIt format and initialize
352  // GLVis visualization.
353  ParGridFunction *u = new ParGridFunction(fes);
354  u->ProjectCoefficient(u0);
355  HypreParVector *U = u->GetTrueDofs();
356 
357  {
358  ostringstream mesh_name, sol_name;
359  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
360  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
361  ofstream omesh(mesh_name.str().c_str());
362  omesh.precision(precision);
363  pmesh->Print(omesh);
364  ofstream osol(sol_name.str().c_str());
365  osol.precision(precision);
366  u->Save(osol);
367  }
368 
369  // Create data collection for solution output: either VisItDataCollection for
370  // ascii data files, or SidreDataCollection for binary data files.
371  DataCollection *dc = NULL;
372  if (visit)
373  {
374  if (binary)
375  {
376 #ifdef MFEM_USE_SIDRE
377  dc = new SidreDataCollection("Example9-Parallel", pmesh);
378 #else
379  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
380 #endif
381  }
382  else
383  {
384  dc = new VisItDataCollection("Example9-Parallel", pmesh);
385  dc->SetPrecision(precision);
386  // To save the mesh using MFEM's parallel mesh format:
387  // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
388  }
389  dc->RegisterField("solution", u);
390  dc->SetCycle(0);
391  dc->SetTime(0.0);
392  dc->Save();
393  }
394 
395  ParaViewDataCollection *pd = NULL;
396  if (paraview)
397  {
398  pd = new ParaViewDataCollection("Example9P", pmesh);
399  pd->SetPrefixPath("ParaView");
400  pd->RegisterField("solution", u);
401  pd->SetLevelsOfDetail(order);
402  pd->SetDataFormat(VTKFormat::BINARY);
403  pd->SetHighOrderOutput(true);
404  pd->SetCycle(0);
405  pd->SetTime(0.0);
406  pd->Save();
407  }
408 
409  // Optionally output a BP (binary pack) file using ADIOS2. This can be
410  // visualized with the ParaView VTX reader.
411 #ifdef MFEM_USE_ADIOS2
412  ADIOS2DataCollection *adios2_dc = NULL;
413  if (adios2)
414  {
415  std::string postfix(mesh_file);
416  postfix.erase(0, std::string("../data/").size() );
417  postfix += "_o" + std::to_string(order);
418  const std::string collection_name = "ex9-p-" + postfix + ".bp";
419 
420  adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh);
421  // output data substreams are half the number of mpi processes
422  adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) );
423  adios2_dc->RegisterField("solution", u);
424  adios2_dc->SetCycle(0);
425  adios2_dc->SetTime(0.0);
426  adios2_dc->Save();
427  }
428 #endif
429 
430  socketstream sout;
431  if (visualization)
432  {
433  char vishost[] = "localhost";
434  int visport = 19916;
435  sout.open(vishost, visport);
436  if (!sout)
437  {
438  if (myid == 0)
439  cout << "Unable to connect to GLVis server at "
440  << vishost << ':' << visport << endl;
441  visualization = false;
442  if (myid == 0)
443  {
444  cout << "GLVis visualization disabled.\n";
445  }
446  }
447  else
448  {
449  sout << "parallel " << num_procs << " " << myid << "\n";
450  sout.precision(precision);
451  sout << "solution\n" << *pmesh << *u;
452  sout << "pause\n";
453  sout << flush;
454  if (myid == 0)
455  cout << "GLVis visualization paused."
456  << " Press space (in the GLVis window) to resume it.\n";
457  }
458  }
459 
460  // 9. Define the time-dependent evolution operator describing the ODE
461  // right-hand side, and define the ODE solver used for time integration.
462  FE_Evolution adv(*m, *k, *B);
463 
464  double t = 0.0;
465  adv.SetTime(t);
466 
467  // Create the time integrator
468  ODESolver *ode_solver = NULL;
469  CVODESolver *cvode = NULL;
470  ARKStepSolver *arkode = NULL;
471  switch (ode_solver_type)
472  {
473  case 1: ode_solver = new ForwardEulerSolver; break;
474  case 2: ode_solver = new RK2Solver(1.0); break;
475  case 3: ode_solver = new RK3SSPSolver; break;
476  case 4: ode_solver = new RK4Solver; break;
477  case 6: ode_solver = new RK6Solver; break;
478  case 7:
479  cvode = new CVODESolver(MPI_COMM_WORLD, CV_ADAMS);
480  cvode->Init(adv);
481  cvode->SetSStolerances(reltol, abstol);
482  cvode->SetMaxStep(dt);
483  cvode->UseSundialsLinearSolver();
484  ode_solver = cvode; break;
485  case 8:
486  case 9:
487  arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::EXPLICIT);
488  arkode->Init(adv);
489  arkode->SetSStolerances(reltol, abstol);
490  arkode->SetMaxStep(dt);
491  if (ode_solver_type == 9)
492  {
494  }
495  ode_solver = arkode; break;
496  }
497 
498  // Initialize MFEM integrators, SUNDIALS integrators are initialized above
499  if (ode_solver_type < 7) { ode_solver->Init(adv); }
500 
501  // 10. Perform time-integration (looping over the time iterations, ti,
502  // with a time-step dt).
503  bool done = false;
504  for (int ti = 0; !done; )
505  {
506  double dt_real = min(dt, t_final - t);
507  ode_solver->Step(*U, t, dt_real);
508  ti++;
509 
510  done = (t >= t_final - 1e-8*dt);
511 
512  if (done || ti % vis_steps == 0)
513  {
514  if (myid == 0)
515  {
516  cout << "time step: " << ti << ", time: " << t << endl;
517  if (cvode) { cvode->PrintInfo(); }
518  if (arkode) { arkode->PrintInfo(); }
519  }
520 
521  // 11. Extract the parallel grid function corresponding to the finite
522  // element approximation U (the local solution on each processor).
523  *u = *U;
524 
525  if (visualization)
526  {
527  sout << "parallel " << num_procs << " " << myid << "\n";
528  sout << "solution\n" << *pmesh << *u << flush;
529  }
530 
531  if (visit)
532  {
533  dc->SetCycle(ti);
534  dc->SetTime(t);
535  dc->Save();
536  }
537 
538  if (paraview)
539  {
540  pd->SetCycle(ti);
541  pd->SetTime(t);
542  pd->Save();
543  }
544 
545 #ifdef MFEM_USE_ADIOS2
546  // transient solutions can be visualized with ParaView
547  if (adios2)
548  {
549  adios2_dc->SetCycle(ti);
550  adios2_dc->SetTime(t);
551  adios2_dc->Save();
552  }
553 #endif
554  }
555  }
556 
557  // 12. Save the final solution in parallel. This output can be viewed later
558  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
559  {
560  *u = *U;
561  ostringstream sol_name;
562  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
563  ofstream osol(sol_name.str().c_str());
564  osol.precision(precision);
565  u->Save(osol);
566  }
567 
568  // 13. Free the used memory.
569  delete U;
570  delete u;
571  delete B;
572  delete b;
573  delete k;
574  delete m;
575  delete fes;
576  delete pmesh;
577  delete ode_solver;
578  delete pd;
579 #ifdef MFEM_USE_ADIOS2
580  if (adios2)
581  {
582  delete adios2_dc;
583  }
584 #endif
585  delete dc;
586 
587  return 0;
588 }
589 
590 
591 // Implementation of class FE_Evolution
593  const Vector &b_)
594  : TimeDependentOperator(M_.Height()),
595  b(b_),
596  M_solver(M_.ParFESpace()->GetComm()),
597  z(M_.Height())
598 {
599  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
600  {
601  M.Reset(M_.ParallelAssemble(), true);
602  K.Reset(K_.ParallelAssemble(), true);
603  }
604  else
605  {
606  M.Reset(&M_, false);
607  K.Reset(&K_, false);
608  }
609 
610  M_solver.SetOperator(*M);
611 
612  Array<int> ess_tdof_list;
613  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
614  {
615  HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
616  HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
617  HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
618  M_prec = hypre_prec;
619 
620  dg_solver = new DG_Solver(M_mat, K_mat, *M_.FESpace());
621  }
622  else
623  {
624  M_prec = new OperatorJacobiSmoother(M_, ess_tdof_list);
625  dg_solver = NULL;
626  }
627 
628  M_solver.SetPreconditioner(*M_prec);
629  M_solver.iterative_mode = false;
630  M_solver.SetRelTol(1e-9);
631  M_solver.SetAbsTol(0.0);
632  M_solver.SetMaxIter(100);
633  M_solver.SetPrintLevel(0);
634 }
635 
636 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
637 {
638  K->Mult(x, z);
639  z += b;
640  dg_solver->SetTimeStep(dt);
641  dg_solver->Mult(z, k);
642 }
643 
644 void FE_Evolution::Mult(const Vector &x, Vector &y) const
645 {
646  // y = M^{-1} (K x + b)
647  K->Mult(x, z);
648  z += b;
649  M_solver.Mult(z, y);
650 }
651 
653 {
654  delete M_prec;
655  delete dg_solver;
656 }
657 
658 
659 // Velocity coefficient
660 void velocity_function(const Vector &x, Vector &v)
661 {
662  int dim = x.Size();
663 
664  // map to the reference [-1,1] domain
665  Vector X(dim);
666  for (int i = 0; i < dim; i++)
667  {
668  double center = (bb_min[i] + bb_max[i]) * 0.5;
669  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
670  }
671 
672  switch (problem)
673  {
674  case 0:
675  {
676  // Translations in 1D, 2D, and 3D
677  switch (dim)
678  {
679  case 1: v(0) = 1.0; break;
680  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
681  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
682  break;
683  }
684  break;
685  }
686  case 1:
687  case 2:
688  {
689  // Clockwise rotation in 2D around the origin
690  const double w = M_PI/2;
691  switch (dim)
692  {
693  case 1: v(0) = 1.0; break;
694  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
695  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
696  }
697  break;
698  }
699  case 3:
700  {
701  // Clockwise twisting rotation in 2D around the origin
702  const double w = M_PI/2;
703  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
704  d = d*d;
705  switch (dim)
706  {
707  case 1: v(0) = 1.0; break;
708  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
709  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
710  }
711  break;
712  }
713  }
714 }
715 
716 // Initial condition
717 double u0_function(const Vector &x)
718 {
719  int dim = x.Size();
720 
721  // map to the reference [-1,1] domain
722  Vector X(dim);
723  for (int i = 0; i < dim; i++)
724  {
725  double center = (bb_min[i] + bb_max[i]) * 0.5;
726  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
727  }
728 
729  switch (problem)
730  {
731  case 0:
732  case 1:
733  {
734  switch (dim)
735  {
736  case 1:
737  return exp(-40.*pow(X(0)-0.5,2));
738  case 2:
739  case 3:
740  {
741  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
742  if (dim == 3)
743  {
744  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
745  rx *= s;
746  ry *= s;
747  }
748  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
749  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
750  }
751  }
752  }
753  case 2:
754  {
755  double x_ = X(0), y_ = X(1), rho, phi;
756  rho = hypot(x_, y_);
757  phi = atan2(y_, x_);
758  return pow(sin(M_PI*rho),2)*sin(3*phi);
759  }
760  case 3:
761  {
762  const double f = M_PI;
763  return sin(f*X(0))*sin(f*X(1));
764  }
765  }
766  return 0.0;
767 }
768 
769 // Inflow boundary condition (zero for the problems considered in this example)
770 double inflow_function(const Vector &x)
771 {
772  switch (problem)
773  {
774  case 0:
775  case 1:
776  case 2:
777  case 3: return 0.0;
778  }
779  return 0.0;
780 }
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:696
Conjugate gradient method.
Definition: solvers.hpp:493
double u0_function(const Vector &x)
Definition: ex9p.cpp:784
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:862
void SetDataFormat(VTKFormat fmt)
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
int Dimension() const
Definition: mesh.hpp:1047
Helper class for ParaView visualization data.
void SetERKTableNum(ARKODE_ERKTableID table_id)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:1718
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
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 Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:712
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:974
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:1706
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:360
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
STL namespace.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
AssemblyLevel GetAssemblyLevel() const
Returns the assembly level.
Vector bb_max
Definition: ex9p.cpp:68
virtual 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:2783
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:661
Data collection with Sidre routines following the Conduit mesh blueprint specification.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2321
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:50
constexpr char vishost[]
virtual void Save()
Save the collection to disk.
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:374
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
constexpr int visport
Vector bb_min
Definition: ex9p.cpp:68
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5356
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:880
void Assemble(int skip_zeros=1)
Assemble the local matrix.
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:107
void SetParameter(const std::string key, const std::string value) noexcept
Parallel smoothers in hypre.
Definition: hypre.hpp:971
void SetHighOrderOutput(bool high_order_output_)
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:1743
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1481
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
void Add(const int i, const int j, const double val)
Definition: sparsemat.cpp:2787
A general vector function coefficient.
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
void SetAbsTol(double atol)
Definition: solvers.hpp:200
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:484
void SetRelTol(double rtol)
Definition: solvers.hpp:199
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:906
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
HYPRE_Int HYPRE_BigInt
GMRES method.
Definition: solvers.hpp:525
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
constexpr ARKODE_ERKTableID ARKODE_FEHLBERG_13_7_8
Definition: sundials.hpp:54
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
Definition: ex18.hpp:81
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
int main(int argc, char *argv[])
Definition: ex9p.cpp:233
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
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;.
RefCoord t[3]
double inflow_function(const Vector &x)
Definition: ex9p.cpp:837
void SetLevelsOfDetail(int levels_of_detail_)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
void velocity_function(const Vector &x, Vector &v)
Definition: ex9p.cpp:727
virtual void Save() override
RefCoord s[3]
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4839
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
The classical forward Euler method.
Definition: ode.hpp:116
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:1461
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:1700
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int problem
Definition: ex9p.cpp:56
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:842
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
double f(const Vector &p)
alpha (q . grad u, v)