MFEM  v4.6.0
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 // Type of preconditioner for implicit time integrator
67 enum class PrecType : int
68 {
69  ILU = 0,
70  AIR = 1
71 };
72 
73 #if MFEM_HYPRE_VERSION >= 21800
74 // Algebraic multigrid preconditioner for advective problems based on
75 // approximate ideal restriction (AIR). Most effective when matrix is
76 // first scaled by DG block inverse, and AIR applied to scaled matrix.
77 // See https://doi.org/10.1137/17M1144350.
78 class AIR_prec : public Solver
79 {
80 private:
81  const HypreParMatrix *A;
82  // Copy of A scaled by block-diagonal inverse
83  HypreParMatrix A_s;
84 
85  HypreBoomerAMG *AIR_solver;
86  int blocksize;
87 
88 public:
89  AIR_prec(int blocksize_) : AIR_solver(NULL), blocksize(blocksize_) { }
90 
91  void SetOperator(const Operator &op)
92  {
93  width = op.Width();
94  height = op.Height();
95 
96  A = dynamic_cast<const HypreParMatrix *>(&op);
97  MFEM_VERIFY(A != NULL, "AIR_prec requires a HypreParMatrix.")
98 
99  // Scale A by block-diagonal inverse
100  BlockInverseScale(A, &A_s, NULL, NULL, blocksize,
101  BlockInverseScaleJob::MATRIX_ONLY);
102  delete AIR_solver;
103  AIR_solver = new HypreBoomerAMG(A_s);
104  AIR_solver->SetAdvectiveOptions(1, "", "FA");
105  AIR_solver->SetPrintLevel(0);
106  AIR_solver->SetMaxLevels(50);
107  }
108 
109  virtual void Mult(const Vector &x, Vector &y) const
110  {
111  // Scale the rhs by block inverse and solve system
112  HypreParVector z_s;
113  BlockInverseScale(A, NULL, &x, &z_s, blocksize,
114  BlockInverseScaleJob::RHS_ONLY);
115  AIR_solver->Mult(z_s, y);
116  }
117 
118  ~AIR_prec()
119  {
120  delete AIR_solver;
121  }
122 };
123 #endif
124 
125 
126 class DG_Solver : public Solver
127 {
128 private:
129  HypreParMatrix &M, &K;
130  SparseMatrix M_diag;
131  HypreParMatrix *A;
132  GMRESSolver linear_solver;
133  Solver *prec;
134  double dt;
135 public:
136  DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
137  PrecType prec_type)
138  : M(M_),
139  K(K_),
140  A(NULL),
141  linear_solver(M.GetComm()),
142  dt(-1.0)
143  {
144  int block_size = fes.GetFE(0)->GetDof();
145  if (prec_type == PrecType::ILU)
146  {
147  prec = new BlockILU(block_size,
148  BlockILU::Reordering::MINIMUM_DISCARDED_FILL);
149  }
150  else if (prec_type == PrecType::AIR)
151  {
152 #if MFEM_HYPRE_VERSION >= 21800
153  prec = new AIR_prec(block_size);
154 #else
155  MFEM_ABORT("Must have MFEM_HYPRE_VERSION >= 21800 to use AIR.\n");
156 #endif
157  }
158  linear_solver.iterative_mode = false;
159  linear_solver.SetRelTol(1e-9);
160  linear_solver.SetAbsTol(0.0);
161  linear_solver.SetMaxIter(100);
162  linear_solver.SetPrintLevel(0);
163  linear_solver.SetPreconditioner(*prec);
164 
165  M.GetDiag(M_diag);
166  }
167 
168  void SetTimeStep(double dt_)
169  {
170  if (dt_ != dt)
171  {
172  dt = dt_;
173  // Form operator A = M - dt*K
174  delete A;
175  A = Add(-dt, K, 0.0, K);
176  SparseMatrix A_diag;
177  A->GetDiag(A_diag);
178  A_diag.Add(1.0, M_diag);
179  // this will also call SetOperator on the preconditioner
180  linear_solver.SetOperator(*A);
181  }
182  }
183 
184  void SetOperator(const Operator &op)
185  {
186  linear_solver.SetOperator(op);
187  }
188 
189  virtual void Mult(const Vector &x, Vector &y) const
190  {
191  linear_solver.Mult(x, y);
192  }
193 
194  ~DG_Solver()
195  {
196  delete prec;
197  delete A;
198  }
199 };
200 
201 
202 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
203  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
204  and advection matrices, and b describes the flow on the boundary. This can
205  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
206  used to evaluate the right-hand side. */
207 class FE_Evolution : public TimeDependentOperator
208 {
209 private:
210  OperatorHandle M, K;
211  const Vector &b;
212  Solver *M_prec;
213  CGSolver M_solver;
214  DG_Solver *dg_solver;
215 
216  mutable Vector z;
217 
218 public:
220  PrecType prec_type);
221 
222  virtual void Mult(const Vector &x, Vector &y) const;
223  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
224 
225  virtual ~FE_Evolution();
226 };
227 
228 
229 int main(int argc, char *argv[])
230 {
231  // 1. Initialize MPI, HYPRE, and SUNDIALS.
232  Mpi::Init(argc, argv);
233  int num_procs = Mpi::WorldSize();
234  int myid = Mpi::WorldRank();
235  Hypre::Init();
236  Sundials::Init();
237 
238  // 2. Parse command-line options.
239  problem = 0;
240  const char *mesh_file = "../../data/periodic-hexagon.mesh";
241  int ser_ref_levels = 2;
242  int par_ref_levels = 0;
243  int order = 3;
244  bool pa = false;
245  bool ea = false;
246  bool fa = false;
247  const char *device_config = "cpu";
248  int ode_solver_type = 7;
249  double t_final = 10.0;
250  double dt = 0.01;
251  bool visualization = false;
252  bool visit = false;
253  bool paraview = false;
254  bool adios2 = false;
255  bool binary = false;
256  int vis_steps = 5;
257 #if MFEM_HYPRE_VERSION >= 21800
258  PrecType prec_type = PrecType::AIR;
259 #else
260  PrecType prec_type = PrecType::ILU;
261 #endif
262 
263  // Relative and absolute tolerances for CVODE and ARKODE.
264  const double reltol = 1e-2, abstol = 1e-2;
265 
266  int precision = 8;
267  cout.precision(precision);
268 
269  OptionsParser args(argc, argv);
270  args.AddOption(&mesh_file, "-m", "--mesh",
271  "Mesh file to use.");
272  args.AddOption(&problem, "-p", "--problem",
273  "Problem setup to use. See options in velocity_function().");
274  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
275  "Number of times to refine the mesh uniformly in serial.");
276  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
277  "Number of times to refine the mesh uniformly in parallel.");
278  args.AddOption(&order, "-o", "--order",
279  "Order (degree) of the finite elements.");
280  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
281  "--no-partial-assembly", "Enable Partial Assembly.");
282  args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
283  "--no-element-assembly", "Enable Element Assembly.");
284  args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
285  "--no-full-assembly", "Enable Full Assembly.");
286  args.AddOption(&device_config, "-d", "--device",
287  "Device configuration string, see Device::Configure().");
288  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
289  "ODE solver:\n\t"
290  "1 - Forward Euler,\n\t"
291  "2 - RK2 SSP,\n\t"
292  "3 - RK3 SSP,\n\t"
293  "4 - RK4,\n\t"
294  "6 - RK6,\n\t"
295  "7 - CVODE (adaptive order implicit Adams),\n\t"
296  "8 - ARKODE default (4th order) explicit,\n\t"
297  "9 - ARKODE RK8.");
298  args.AddOption(&t_final, "-tf", "--t-final",
299  "Final time; start time is 0.");
300  args.AddOption(&dt, "-dt", "--time-step",
301  "Time step.");
302  args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
303  "implicit solves. 0 for ILU, 1 for pAIR-AMG.");
304  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
305  "--no-visualization",
306  "Enable or disable GLVis visualization.");
307  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
308  "--no-visit-datafiles",
309  "Save data files for VisIt (visit.llnl.gov) visualization.");
310  args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
311  "--no-paraview-datafiles",
312  "Save data files for ParaView (paraview.org) visualization.");
313  args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
314  "--no-adios2-streams",
315  "Save data using adios2 streams.");
316  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
317  "--ascii-datafiles",
318  "Use binary (Sidre) or ascii format for VisIt data files.");
319  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
320  "Visualize every n-th timestep.");
321  args.Parse();
322  if (!args.Good())
323  {
324  if (Mpi::Root())
325  {
326  args.PrintUsage(cout);
327  }
328  return 1;
329  }
330  if (Mpi::Root())
331  {
332  args.PrintOptions(cout);
333  }
334 
335  // check for valid ODE solver option
336  if (ode_solver_type < 1 || ode_solver_type > 9)
337  {
338  if (Mpi::Root())
339  {
340  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
341  }
342  return 3;
343  }
344 
345  Device device(device_config);
346  if (Mpi::Root()) { device.Print(); }
347 
348  // 3. Read the serial mesh from the given mesh file on all processors. We can
349  // handle geometrically periodic meshes in this code.
350  Mesh *mesh = new Mesh(mesh_file, 1, 1);
351  int dim = mesh->Dimension();
352 
353  // 4. Refine the mesh in serial to increase the resolution. In this example
354  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
355  // a command-line parameter. If the mesh is of NURBS type, we convert it
356  // to a (piecewise-polynomial) high-order mesh.
357  for (int lev = 0; lev < ser_ref_levels; lev++)
358  {
359  mesh->UniformRefinement();
360  }
361  if (mesh->NURBSext)
362  {
363  mesh->SetCurvature(max(order, 1));
364  }
365  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
366 
367  // 5. Define the parallel mesh by a partitioning of the serial mesh. Refine
368  // this mesh further in parallel to increase the resolution. Once the
369  // parallel mesh is defined, the serial mesh can be deleted.
370  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
371  delete mesh;
372  for (int lev = 0; lev < par_ref_levels; lev++)
373  {
374  pmesh->UniformRefinement();
375  }
376 
377  // 6. Define the parallel discontinuous DG finite element space on the
378  // parallel refined mesh of the given polynomial order.
379  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
380  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
381 
382  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
383  if (Mpi::Root())
384  {
385  cout << "Number of unknowns: " << global_vSize << endl;
386  }
387 
388  // 7. Set up and assemble the parallel bilinear and linear forms (and the
389  // parallel hypre matrices) corresponding to the DG discretization. The
390  // DGTraceIntegrator involves integrals over mesh interior faces.
394 
395  ParBilinearForm *m = new ParBilinearForm(fes);
396  ParBilinearForm *k = new ParBilinearForm(fes);
397  if (pa)
398  {
399  m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
400  k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
401  }
402  else if (ea)
403  {
404  m->SetAssemblyLevel(AssemblyLevel::ELEMENT);
405  k->SetAssemblyLevel(AssemblyLevel::ELEMENT);
406  }
407  else if (fa)
408  {
409  m->SetAssemblyLevel(AssemblyLevel::FULL);
410  k->SetAssemblyLevel(AssemblyLevel::FULL);
411  }
412 
414  constexpr double alpha = -1.0;
417  new NonconservativeDGTraceIntegrator(velocity, alpha));
419  new NonconservativeDGTraceIntegrator(velocity, alpha));
420 
421  ParLinearForm *b = new ParLinearForm(fes);
422  b->AddBdrFaceIntegrator(
423  new BoundaryFlowIntegrator(inflow, velocity, alpha));
424 
425  int skip_zeros = 0;
426  m->Assemble();
427  k->Assemble(skip_zeros);
428  b->Assemble();
429  m->Finalize();
430  k->Finalize(skip_zeros);
431 
432  HypreParVector *B = b->ParallelAssemble();
433 
434  // 8. Define the initial conditions, save the corresponding grid function to
435  // a file and (optionally) save data in the VisIt format and initialize
436  // GLVis visualization.
437  ParGridFunction *u = new ParGridFunction(fes);
438  u->ProjectCoefficient(u0);
439  HypreParVector *U = u->GetTrueDofs();
440 
441  {
442  ostringstream mesh_name, sol_name;
443  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
444  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
445  ofstream omesh(mesh_name.str().c_str());
446  omesh.precision(precision);
447  pmesh->Print(omesh);
448  ofstream osol(sol_name.str().c_str());
449  osol.precision(precision);
450  u->Save(osol);
451  }
452 
453  // Create data collection for solution output: either VisItDataCollection for
454  // ascii data files, or SidreDataCollection for binary data files.
455  DataCollection *dc = NULL;
456  if (visit)
457  {
458  if (binary)
459  {
460 #ifdef MFEM_USE_SIDRE
461  dc = new SidreDataCollection("Example9-Parallel", pmesh);
462 #else
463  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
464 #endif
465  }
466  else
467  {
468  dc = new VisItDataCollection("Example9-Parallel", pmesh);
469  dc->SetPrecision(precision);
470  // To save the mesh using MFEM's parallel mesh format:
471  // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
472  }
473  dc->RegisterField("solution", u);
474  dc->SetCycle(0);
475  dc->SetTime(0.0);
476  dc->Save();
477  }
478 
479  ParaViewDataCollection *pd = NULL;
480  if (paraview)
481  {
482  pd = new ParaViewDataCollection("Example9P", pmesh);
483  pd->SetPrefixPath("ParaView");
484  pd->RegisterField("solution", u);
485  pd->SetLevelsOfDetail(order);
486  pd->SetDataFormat(VTKFormat::BINARY);
487  pd->SetHighOrderOutput(true);
488  pd->SetCycle(0);
489  pd->SetTime(0.0);
490  pd->Save();
491  }
492 
493  // Optionally output a BP (binary pack) file using ADIOS2. This can be
494  // visualized with the ParaView VTX reader.
495 #ifdef MFEM_USE_ADIOS2
496  ADIOS2DataCollection *adios2_dc = NULL;
497  if (adios2)
498  {
499  std::string postfix(mesh_file);
500  postfix.erase(0, std::string("../data/").size() );
501  postfix += "_o" + std::to_string(order);
502  const std::string collection_name = "ex9-p-" + postfix + ".bp";
503 
504  adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh);
505  // output data substreams are half the number of mpi processes
506  adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) );
507  adios2_dc->RegisterField("solution", u);
508  adios2_dc->SetCycle(0);
509  adios2_dc->SetTime(0.0);
510  adios2_dc->Save();
511  }
512 #endif
513 
514  socketstream sout;
515  if (visualization)
516  {
517  char vishost[] = "localhost";
518  int visport = 19916;
519  sout.open(vishost, visport);
520  if (!sout)
521  {
522  if (Mpi::Root())
523  {
524  cout << "Unable to connect to GLVis server at "
525  << vishost << ':' << visport << endl;
526  }
527  visualization = false;
528  if (Mpi::Root())
529  {
530  cout << "GLVis visualization disabled.\n";
531  }
532  }
533  else
534  {
535  sout << "parallel " << num_procs << " " << myid << "\n";
536  sout.precision(precision);
537  sout << "solution\n" << *pmesh << *u;
538  sout << "pause\n";
539  sout << flush;
540  if (Mpi::Root())
541  {
542  cout << "GLVis visualization paused."
543  << " Press space (in the GLVis window) to resume it.\n";
544  }
545  }
546  }
547 
548  // 9. Define the time-dependent evolution operator describing the ODE
549  // right-hand side, and define the ODE solver used for time integration.
550  FE_Evolution adv(*m, *k, *B, prec_type);
551 
552  double t = 0.0;
553  adv.SetTime(t);
554 
555  // Create the time integrator
556  ODESolver *ode_solver = NULL;
557  CVODESolver *cvode = NULL;
558  ARKStepSolver *arkode = NULL;
559  switch (ode_solver_type)
560  {
561  case 1: ode_solver = new ForwardEulerSolver; break;
562  case 2: ode_solver = new RK2Solver(1.0); break;
563  case 3: ode_solver = new RK3SSPSolver; break;
564  case 4: ode_solver = new RK4Solver; break;
565  case 6: ode_solver = new RK6Solver; break;
566  case 7:
567  cvode = new CVODESolver(MPI_COMM_WORLD, CV_ADAMS);
568  cvode->Init(adv);
569  cvode->SetSStolerances(reltol, abstol);
570  cvode->SetMaxStep(dt);
571  cvode->UseSundialsLinearSolver();
572  ode_solver = cvode; break;
573  case 8:
574  case 9:
575  arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::EXPLICIT);
576  arkode->Init(adv);
577  arkode->SetSStolerances(reltol, abstol);
578  arkode->SetMaxStep(dt);
579  if (ode_solver_type == 9)
580  {
582  }
583  ode_solver = arkode; break;
584  }
585 
586  // Initialize MFEM integrators, SUNDIALS integrators are initialized above
587  if (ode_solver_type < 7) { ode_solver->Init(adv); }
588 
589  // 10. Perform time-integration (looping over the time iterations, ti,
590  // with a time-step dt).
591  bool done = false;
592  for (int ti = 0; !done; )
593  {
594  double dt_real = min(dt, t_final - t);
595  ode_solver->Step(*U, t, dt_real);
596  ti++;
597 
598  done = (t >= t_final - 1e-8*dt);
599 
600  if (done || ti % vis_steps == 0)
601  {
602  if (Mpi::Root())
603  {
604  cout << "time step: " << ti << ", time: " << t << endl;
605  if (cvode) { cvode->PrintInfo(); }
606  if (arkode) { arkode->PrintInfo(); }
607  }
608 
609  // 11. Extract the parallel grid function corresponding to the finite
610  // element approximation U (the local solution on each processor).
611  *u = *U;
612 
613  if (visualization)
614  {
615  sout << "parallel " << num_procs << " " << myid << "\n";
616  sout << "solution\n" << *pmesh << *u << flush;
617  }
618 
619  if (visit)
620  {
621  dc->SetCycle(ti);
622  dc->SetTime(t);
623  dc->Save();
624  }
625 
626  if (paraview)
627  {
628  pd->SetCycle(ti);
629  pd->SetTime(t);
630  pd->Save();
631  }
632 
633 #ifdef MFEM_USE_ADIOS2
634  // transient solutions can be visualized with ParaView
635  if (adios2)
636  {
637  adios2_dc->SetCycle(ti);
638  adios2_dc->SetTime(t);
639  adios2_dc->Save();
640  }
641 #endif
642  }
643  }
644 
645  // 12. Save the final solution in parallel. This output can be viewed later
646  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
647  {
648  *u = *U;
649  ostringstream sol_name;
650  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
651  ofstream osol(sol_name.str().c_str());
652  osol.precision(precision);
653  u->Save(osol);
654  }
655 
656  // 13. Free the used memory.
657  delete U;
658  delete u;
659  delete B;
660  delete b;
661  delete k;
662  delete m;
663  delete fes;
664  delete pmesh;
665  delete ode_solver;
666  delete pd;
667 #ifdef MFEM_USE_ADIOS2
668  if (adios2)
669  {
670  delete adios2_dc;
671  }
672 #endif
673  delete dc;
674 
675  return 0;
676 }
677 
678 
679 // Implementation of class FE_Evolution
681  const Vector &b_, PrecType prec_type)
682  : TimeDependentOperator(M_.Height()),
683  b(b_),
684  M_solver(M_.ParFESpace()->GetComm()),
685  z(M_.Height())
686 {
687  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
688  {
689  M.Reset(M_.ParallelAssemble(), true);
690  K.Reset(K_.ParallelAssemble(), true);
691  }
692  else
693  {
694  M.Reset(&M_, false);
695  K.Reset(&K_, false);
696  }
697 
698  M_solver.SetOperator(*M);
699 
700  Array<int> ess_tdof_list;
701  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
702  {
703  HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
704  HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
705  HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
706  M_prec = hypre_prec;
707 
708  dg_solver = new DG_Solver(M_mat, K_mat, *M_.FESpace(), prec_type);
709  }
710  else
711  {
712  M_prec = new OperatorJacobiSmoother(M_, ess_tdof_list);
713  dg_solver = NULL;
714  }
715 
716  M_solver.SetPreconditioner(*M_prec);
717  M_solver.iterative_mode = false;
718  M_solver.SetRelTol(1e-9);
719  M_solver.SetAbsTol(0.0);
720  M_solver.SetMaxIter(100);
721  M_solver.SetPrintLevel(0);
722 }
723 
724 // Solve the equation:
725 // u_t = M^{-1}(Ku + b),
726 // by solving associated linear system
727 // (M - dt*K) d = K*u + b
728 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
729 {
730  K->Mult(x, z);
731  z += b;
732  dg_solver->SetTimeStep(dt);
733  dg_solver->Mult(z, k);
734 }
735 
736 void FE_Evolution::Mult(const Vector &x, Vector &y) const
737 {
738  // y = M^{-1} (K x + b)
739  K->Mult(x, z);
740  z += b;
741  M_solver.Mult(z, y);
742 }
743 
745 {
746  delete M_prec;
747  delete dg_solver;
748 }
749 
750 
751 // Velocity coefficient
752 void velocity_function(const Vector &x, Vector &v)
753 {
754  int dim = x.Size();
755 
756  // map to the reference [-1,1] domain
757  Vector X(dim);
758  for (int i = 0; i < dim; i++)
759  {
760  double center = (bb_min[i] + bb_max[i]) * 0.5;
761  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
762  }
763 
764  switch (problem)
765  {
766  case 0:
767  {
768  // Translations in 1D, 2D, and 3D
769  switch (dim)
770  {
771  case 1: v(0) = 1.0; break;
772  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
773  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
774  break;
775  }
776  break;
777  }
778  case 1:
779  case 2:
780  {
781  // Clockwise rotation in 2D around the origin
782  const double w = M_PI/2;
783  switch (dim)
784  {
785  case 1: v(0) = 1.0; break;
786  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
787  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
788  }
789  break;
790  }
791  case 3:
792  {
793  // Clockwise twisting rotation in 2D around the origin
794  const double w = M_PI/2;
795  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
796  d = d*d;
797  switch (dim)
798  {
799  case 1: v(0) = 1.0; break;
800  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
801  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
802  }
803  break;
804  }
805  }
806 }
807 
808 // Initial condition
809 double u0_function(const Vector &x)
810 {
811  int dim = x.Size();
812 
813  // map to the reference [-1,1] domain
814  Vector X(dim);
815  for (int i = 0; i < dim; i++)
816  {
817  double center = (bb_min[i] + bb_max[i]) * 0.5;
818  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
819  }
820 
821  switch (problem)
822  {
823  case 0:
824  case 1:
825  {
826  switch (dim)
827  {
828  case 1:
829  return exp(-40.*pow(X(0)-0.5,2));
830  case 2:
831  case 3:
832  {
833  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
834  if (dim == 3)
835  {
836  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
837  rx *= s;
838  ry *= s;
839  }
840  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
841  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
842  }
843  }
844  }
845  case 2:
846  {
847  double x_ = X(0), y_ = X(1), rho, phi;
848  rho = hypot(x_, y_);
849  phi = atan2(y_, x_);
850  return pow(sin(M_PI*rho),2)*sin(3*phi);
851  }
852  case 3:
853  {
854  const double f = M_PI;
855  return sin(f*X(0))*sin(f*X(1));
856  }
857  }
858  return 0.0;
859 }
860 
861 // Inflow boundary condition (zero for the problems considered in this example)
862 double inflow_function(const Vector &x)
863 {
864  switch (problem)
865  {
866  case 0:
867  case 1:
868  case 2:
869  case 3: return 0.0;
870  }
871  return 0.0;
872 }
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:709
int visport
Conjugate gradient method.
Definition: solvers.hpp:493
double u0_function(const Vector &x)
Definition: ex9p.cpp:788
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:875
void SetDataFormat(VTKFormat fmt)
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
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:1731
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
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:136
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
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:718
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:980
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:23
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:1719
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:159
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.
void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C, const Vector *b, HypreParVector *d, int blocksize, BlockInverseScaleJob job)
Definition: hypre.cpp:2766
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:2841
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:672
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:2320
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
Class for parallel linear form.
Definition: plinearform.hpp:26
PrecType
Definition: ex9p.cpp:71
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
char vishost[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
BlockInverseScaleJob
Definition: hypre.hpp:915
virtual void Save()
Save the collection to disk.
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:385
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:10232
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:5635
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:893
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:1756
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:281
void Add(const int i, const int j, const double val)
Definition: sparsemat.cpp:2787
A general vector function coefficient.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:163
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:919
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
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:150
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:323
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
constexpr ARKODE_ERKTableID ARKODE_FEHLBERG_13_7_8
Definition: sundials.hpp:65
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:841
void SetLevelsOfDetail(int levels_of_detail_)
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
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:731
virtual void Save() override
RefCoord s[3]
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4825
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
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:117
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:1474
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:1713
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1736
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:855
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
double f(const Vector &p)
alpha (q . grad u, v)