MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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.
156  int num_procs, myid;
157  MPI_Init(&argc, &argv);
158  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
159  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
160 
161  // 2. Parse command-line options.
162  problem = 0;
163  const char *mesh_file = "../../data/periodic-hexagon.mesh";
164  int ser_ref_levels = 2;
165  int par_ref_levels = 0;
166  int order = 3;
167  bool pa = false;
168  bool ea = false;
169  bool fa = false;
170  const char *device_config = "cpu";
171  int ode_solver_type = 7;
172  double t_final = 10.0;
173  double dt = 0.01;
174  bool visualization = false;
175  bool visit = false;
176  bool paraview = false;
177  bool adios2 = false;
178  bool binary = false;
179  int vis_steps = 5;
180 
181  // Relative and absolute tolerances for CVODE and ARKODE.
182  const double reltol = 1e-2, abstol = 1e-2;
183 
184  int precision = 8;
185  cout.precision(precision);
186 
187  OptionsParser args(argc, argv);
188  args.AddOption(&mesh_file, "-m", "--mesh",
189  "Mesh file to use.");
190  args.AddOption(&problem, "-p", "--problem",
191  "Problem setup to use. See options in velocity_function().");
192  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
193  "Number of times to refine the mesh uniformly in serial.");
194  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
195  "Number of times to refine the mesh uniformly in parallel.");
196  args.AddOption(&order, "-o", "--order",
197  "Order (degree) of the finite elements.");
198  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
199  "--no-partial-assembly", "Enable Partial Assembly.");
200  args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
201  "--no-element-assembly", "Enable Element Assembly.");
202  args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
203  "--no-full-assembly", "Enable Full Assembly.");
204  args.AddOption(&device_config, "-d", "--device",
205  "Device configuration string, see Device::Configure().");
206  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
207  "ODE solver:\n\t"
208  "1 - Forward Euler,\n\t"
209  "2 - RK2 SSP,\n\t"
210  "3 - RK3 SSP,\n\t"
211  "4 - RK4,\n\t"
212  "6 - RK6,\n\t"
213  "7 - CVODE (adaptive order implicit Adams),\n\t"
214  "8 - ARKODE default (4th order) explicit,\n\t"
215  "9 - ARKODE RK8.");
216  args.AddOption(&t_final, "-tf", "--t-final",
217  "Final time; start time is 0.");
218  args.AddOption(&dt, "-dt", "--time-step",
219  "Time step.");
220  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
221  "--no-visualization",
222  "Enable or disable GLVis visualization.");
223  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
224  "--no-visit-datafiles",
225  "Save data files for VisIt (visit.llnl.gov) visualization.");
226  args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
227  "--no-paraview-datafiles",
228  "Save data files for ParaView (paraview.org) visualization.");
229  args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
230  "--no-adios2-streams",
231  "Save data using adios2 streams.");
232  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
233  "--ascii-datafiles",
234  "Use binary (Sidre) or ascii format for VisIt data files.");
235  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
236  "Visualize every n-th timestep.");
237  args.Parse();
238  if (!args.Good())
239  {
240  if (myid == 0)
241  {
242  args.PrintUsage(cout);
243  }
244  MPI_Finalize();
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  MPI_Finalize();
260  return 3;
261  }
262 
263  Device device(device_config);
264  if (myid == 0) { device.Print(); }
265 
266  // 3. Read the serial mesh from the given mesh file on all processors. We can
267  // handle geometrically periodic meshes in this code.
268  Mesh *mesh = new Mesh(mesh_file, 1, 1);
269  int dim = mesh->Dimension();
270 
271  // 4. Refine the mesh in serial to increase the resolution. In this example
272  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
273  // a command-line parameter. If the mesh is of NURBS type, we convert it
274  // to a (piecewise-polynomial) high-order mesh.
275  for (int lev = 0; lev < ser_ref_levels; lev++)
276  {
277  mesh->UniformRefinement();
278  }
279  if (mesh->NURBSext)
280  {
281  mesh->SetCurvature(max(order, 1));
282  }
283  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
284 
285  // 5. Define the parallel mesh by a partitioning of the serial mesh. Refine
286  // this mesh further in parallel to increase the resolution. Once the
287  // parallel mesh is defined, the serial mesh can be deleted.
288  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
289  delete mesh;
290  for (int lev = 0; lev < par_ref_levels; lev++)
291  {
292  pmesh->UniformRefinement();
293  }
294 
295  // 6. Define the parallel discontinuous DG finite element space on the
296  // parallel refined mesh of the given polynomial order.
297  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
298  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
299 
300  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
301  if (myid == 0)
302  {
303  cout << "Number of unknowns: " << global_vSize << endl;
304  }
305 
306  // 7. Set up and assemble the parallel bilinear and linear forms (and the
307  // parallel hypre matrices) corresponding to the DG discretization. The
308  // DGTraceIntegrator involves integrals over mesh interior faces.
312 
313  ParBilinearForm *m = new ParBilinearForm(fes);
314  ParBilinearForm *k = new ParBilinearForm(fes);
315  if (pa)
316  {
317  m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
318  k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
319  }
320  else if (ea)
321  {
322  m->SetAssemblyLevel(AssemblyLevel::ELEMENT);
323  k->SetAssemblyLevel(AssemblyLevel::ELEMENT);
324  }
325  else if (fa)
326  {
327  m->SetAssemblyLevel(AssemblyLevel::FULL);
328  k->SetAssemblyLevel(AssemblyLevel::FULL);
329  }
330 
332  k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
334  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
336  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
337 
338  ParLinearForm *b = new ParLinearForm(fes);
340  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
341 
342  int skip_zeros = 0;
343  m->Assemble();
344  k->Assemble(skip_zeros);
345  b->Assemble();
346  m->Finalize();
347  k->Finalize(skip_zeros);
348 
350 
351  // 8. Define the initial conditions, save the corresponding grid function to
352  // a file and (optionally) save data in the VisIt format and initialize
353  // GLVis visualization.
354  ParGridFunction *u = new ParGridFunction(fes);
355  u->ProjectCoefficient(u0);
356  HypreParVector *U = u->GetTrueDofs();
357 
358  {
359  ostringstream mesh_name, sol_name;
360  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
361  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
362  ofstream omesh(mesh_name.str().c_str());
363  omesh.precision(precision);
364  pmesh->Print(omesh);
365  ofstream osol(sol_name.str().c_str());
366  osol.precision(precision);
367  u->Save(osol);
368  }
369 
370  // Create data collection for solution output: either VisItDataCollection for
371  // ascii data files, or SidreDataCollection for binary data files.
372  DataCollection *dc = NULL;
373  if (visit)
374  {
375  if (binary)
376  {
377 #ifdef MFEM_USE_SIDRE
378  dc = new SidreDataCollection("Example9-Parallel", pmesh);
379 #else
380  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
381 #endif
382  }
383  else
384  {
385  dc = new VisItDataCollection("Example9-Parallel", pmesh);
386  dc->SetPrecision(precision);
387  // To save the mesh using MFEM's parallel mesh format:
388  // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
389  }
390  dc->RegisterField("solution", u);
391  dc->SetCycle(0);
392  dc->SetTime(0.0);
393  dc->Save();
394  }
395 
396  ParaViewDataCollection *pd = NULL;
397  if (paraview)
398  {
399  pd = new ParaViewDataCollection("Example9P", pmesh);
400  pd->SetPrefixPath("ParaView");
401  pd->RegisterField("solution", u);
402  pd->SetLevelsOfDetail(order);
403  pd->SetDataFormat(VTKFormat::BINARY);
404  pd->SetHighOrderOutput(true);
405  pd->SetCycle(0);
406  pd->SetTime(0.0);
407  pd->Save();
408  }
409 
410  // Optionally output a BP (binary pack) file using ADIOS2. This can be
411  // visualized with the ParaView VTX reader.
412 #ifdef MFEM_USE_ADIOS2
413  ADIOS2DataCollection *adios2_dc = NULL;
414  if (adios2)
415  {
416  std::string postfix(mesh_file);
417  postfix.erase(0, std::string("../data/").size() );
418  postfix += "_o" + std::to_string(order);
419  const std::string collection_name = "ex9-p-" + postfix + ".bp";
420 
421  adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh);
422  // output data substreams are half the number of mpi processes
423  adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) );
424  adios2_dc->RegisterField("solution", u);
425  adios2_dc->SetCycle(0);
426  adios2_dc->SetTime(0.0);
427  adios2_dc->Save();
428  }
429 #endif
430 
431  socketstream sout;
432  if (visualization)
433  {
434  char vishost[] = "localhost";
435  int visport = 19916;
436  sout.open(vishost, visport);
437  if (!sout)
438  {
439  if (myid == 0)
440  cout << "Unable to connect to GLVis server at "
441  << vishost << ':' << visport << endl;
442  visualization = false;
443  if (myid == 0)
444  {
445  cout << "GLVis visualization disabled.\n";
446  }
447  }
448  else
449  {
450  sout << "parallel " << num_procs << " " << myid << "\n";
451  sout.precision(precision);
452  sout << "solution\n" << *pmesh << *u;
453  sout << "pause\n";
454  sout << flush;
455  if (myid == 0)
456  cout << "GLVis visualization paused."
457  << " Press space (in the GLVis window) to resume it.\n";
458  }
459  }
460 
461  // 9. Define the time-dependent evolution operator describing the ODE
462  // right-hand side, and define the ODE solver used for time integration.
463  FE_Evolution adv(*m, *k, *B);
464 
465  double t = 0.0;
466  adv.SetTime(t);
467 
468  // Create the time integrator
469  ODESolver *ode_solver = NULL;
470  CVODESolver *cvode = NULL;
471  ARKStepSolver *arkode = NULL;
472  switch (ode_solver_type)
473  {
474  case 1: ode_solver = new ForwardEulerSolver; break;
475  case 2: ode_solver = new RK2Solver(1.0); break;
476  case 3: ode_solver = new RK3SSPSolver; break;
477  case 4: ode_solver = new RK4Solver; break;
478  case 6: ode_solver = new RK6Solver; break;
479  case 7:
480  cvode = new CVODESolver(MPI_COMM_WORLD, CV_ADAMS);
481  cvode->Init(adv);
482  cvode->SetSStolerances(reltol, abstol);
483  cvode->SetMaxStep(dt);
484  cvode->UseSundialsLinearSolver();
485  ode_solver = cvode; break;
486  case 8:
487  case 9:
488  arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::EXPLICIT);
489  arkode->Init(adv);
490  arkode->SetSStolerances(reltol, abstol);
491  arkode->SetMaxStep(dt);
492  if (ode_solver_type == 9) { arkode->SetERKTableNum(FEHLBERG_13_7_8); }
493  ode_solver = arkode; break;
494  }
495 
496  // Initialize MFEM integrators, SUNDIALS integrators are initialized above
497  if (ode_solver_type < 7) { ode_solver->Init(adv); }
498 
499  // 10. Perform time-integration (looping over the time iterations, ti,
500  // with a time-step dt).
501  bool done = false;
502  for (int ti = 0; !done; )
503  {
504  double dt_real = min(dt, t_final - t);
505  ode_solver->Step(*U, t, dt_real);
506  ti++;
507 
508  done = (t >= t_final - 1e-8*dt);
509 
510  if (done || ti % vis_steps == 0)
511  {
512  if (myid == 0)
513  {
514  cout << "time step: " << ti << ", time: " << t << endl;
515  if (cvode) { cvode->PrintInfo(); }
516  if (arkode) { arkode->PrintInfo(); }
517  }
518 
519  // 11. Extract the parallel grid function corresponding to the finite
520  // element approximation U (the local solution on each processor).
521  *u = *U;
522 
523  if (visualization)
524  {
525  sout << "parallel " << num_procs << " " << myid << "\n";
526  sout << "solution\n" << *pmesh << *u << flush;
527  }
528 
529  if (visit)
530  {
531  dc->SetCycle(ti);
532  dc->SetTime(t);
533  dc->Save();
534  }
535 
536  if (paraview)
537  {
538  pd->SetCycle(ti);
539  pd->SetTime(t);
540  pd->Save();
541  }
542 
543 #ifdef MFEM_USE_ADIOS2
544  // transient solutions can be visualized with ParaView
545  if (adios2)
546  {
547  adios2_dc->SetCycle(ti);
548  adios2_dc->SetTime(t);
549  adios2_dc->Save();
550  }
551 #endif
552  }
553  }
554 
555  // 12. Save the final solution in parallel. This output can be viewed later
556  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
557  {
558  *u = *U;
559  ostringstream sol_name;
560  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
561  ofstream osol(sol_name.str().c_str());
562  osol.precision(precision);
563  u->Save(osol);
564  }
565 
566  // 13. Free the used memory.
567  delete U;
568  delete u;
569  delete B;
570  delete b;
571  delete k;
572  delete m;
573  delete fes;
574  delete pmesh;
575  delete ode_solver;
576  delete pd;
577 #ifdef MFEM_USE_ADIOS2
578  if (adios2)
579  {
580  delete adios2_dc;
581  }
582 #endif
583  delete dc;
584 
585  MPI_Finalize();
586  return 0;
587 }
588 
589 
590 // Implementation of class FE_Evolution
592  const Vector &b_)
593  : TimeDependentOperator(M_.Height()),
594  b(b_),
595  M_solver(M_.ParFESpace()->GetComm()),
596  z(M_.Height())
597 {
598  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
599  {
600  M.Reset(M_.ParallelAssemble(), true);
601  K.Reset(K_.ParallelAssemble(), true);
602  }
603  else
604  {
605  M.Reset(&M_, false);
606  K.Reset(&K_, false);
607  }
608 
609  M_solver.SetOperator(*M);
610 
611  Array<int> ess_tdof_list;
612  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
613  {
614  HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
615  HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
616  HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
617  M_prec = hypre_prec;
618 
619  dg_solver = new DG_Solver(M_mat, K_mat, *M_.FESpace());
620  }
621  else
622  {
623  M_prec = new OperatorJacobiSmoother(M_, ess_tdof_list);
624  dg_solver = NULL;
625  }
626 
627  M_solver.SetPreconditioner(*M_prec);
628  M_solver.iterative_mode = false;
629  M_solver.SetRelTol(1e-9);
630  M_solver.SetAbsTol(0.0);
631  M_solver.SetMaxIter(100);
632  M_solver.SetPrintLevel(0);
633 }
634 
635 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
636 {
637  K->Mult(x, z);
638  z += b;
639  dg_solver->SetTimeStep(dt);
640  dg_solver->Mult(z, k);
641 }
642 
643 void FE_Evolution::Mult(const Vector &x, Vector &y) const
644 {
645  // y = M^{-1} (K x + b)
646  K->Mult(x, z);
647  z += b;
648  M_solver.Mult(z, y);
649 }
650 
652 {
653  delete M_prec;
654  delete dg_solver;
655 }
656 
657 
658 // Velocity coefficient
659 void velocity_function(const Vector &x, Vector &v)
660 {
661  int dim = x.Size();
662 
663  // map to the reference [-1,1] domain
664  Vector X(dim);
665  for (int i = 0; i < dim; i++)
666  {
667  double center = (bb_min[i] + bb_max[i]) * 0.5;
668  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
669  }
670 
671  switch (problem)
672  {
673  case 0:
674  {
675  // Translations in 1D, 2D, and 3D
676  switch (dim)
677  {
678  case 1: v(0) = 1.0; break;
679  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
680  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
681  break;
682  }
683  break;
684  }
685  case 1:
686  case 2:
687  {
688  // Clockwise rotation in 2D around the origin
689  const double w = M_PI/2;
690  switch (dim)
691  {
692  case 1: v(0) = 1.0; break;
693  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
694  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
695  }
696  break;
697  }
698  case 3:
699  {
700  // Clockwise twisting rotation in 2D around the origin
701  const double w = M_PI/2;
702  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
703  d = d*d;
704  switch (dim)
705  {
706  case 1: v(0) = 1.0; break;
707  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
708  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
709  }
710  break;
711  }
712  }
713 }
714 
715 // Initial condition
716 double u0_function(const Vector &x)
717 {
718  int dim = x.Size();
719 
720  // map to the reference [-1,1] domain
721  Vector X(dim);
722  for (int i = 0; i < dim; i++)
723  {
724  double center = (bb_min[i] + bb_max[i]) * 0.5;
725  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
726  }
727 
728  switch (problem)
729  {
730  case 0:
731  case 1:
732  {
733  switch (dim)
734  {
735  case 1:
736  return exp(-40.*pow(X(0)-0.5,2));
737  case 2:
738  case 3:
739  {
740  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
741  if (dim == 3)
742  {
743  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
744  rx *= s;
745  ry *= s;
746  }
747  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
748  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
749  }
750  }
751  }
752  case 2:
753  {
754  double x_ = X(0), y_ = X(1), rho, phi;
755  rho = hypot(x_, y_);
756  phi = atan2(y_, x_);
757  return pow(sin(M_PI*rho),2)*sin(3*phi);
758  }
759  case 3:
760  {
761  const double f = M_PI;
762  return sin(f*X(0))*sin(f*X(1));
763  }
764  }
765  return 0.0;
766 }
767 
768 // Inflow boundary condition (zero for the problems considered in this example)
769 double inflow_function(const Vector &x)
770 {
771  switch (problem)
772  {
773  case 0:
774  case 1:
775  case 2:
776  case 3: return 0.0;
777  }
778  return 0.0;
779 }
Vector bb_max
Definition: ex9.cpp:66
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
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:2577
Conjugate gradient method.
Definition: solvers.hpp:316
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
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:698
void SetDataFormat(VTKFormat fmt)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:616
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
Helper class for ParaView visualization data.
Base abstract class for first order time dependent operators.
Definition: operator.hpp:282
AssemblyLevel GetAssemblyLevel() const
Returns the assembly level.
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: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]...
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:1540
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:326
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:493
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
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.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1931
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.
void SetPrintLevel(int print_lvl)
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:150
Data type sparse matrix.
Definition: sparsemat.hpp:41
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:130
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4882
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:716
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
void SetParameter(const std::string key, const std::string value) noexcept
Parallel smoothers in hypre.
Definition: hypre.hpp:840
void SetHighOrderOutput(bool high_order_output_)
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:83
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:66
A general vector function coefficient.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:99
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
void SetAbsTol(double atol)
Definition: solvers.hpp:99
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:483
void SetRelTol(double rtol)
Definition: solvers.hpp:98
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:500
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:323
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:62
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
HYPRE_Int HYPRE_BigInt
GMRES method.
Definition: solvers.hpp:348
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
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
Definition: ex18.hpp:102
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:206
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:549
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.
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:2388
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
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;.
double u0_function(const Vector &x)
Definition: ex9.cpp:557
RefCoord t[3]
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:330
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:92
virtual void Save() override
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:742
RefCoord s[3]
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Base class for solvers.
Definition: operator.hpp:648
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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: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:277
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:1297
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:1534
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:87
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int main()
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
double inflow_function(const Vector &x)
Definition: ex9.cpp:610
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:285
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
alpha (q . grad u, v)