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