MFEM  v4.1.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 //
20 // Device sample runs:
21 // mpirun -np 4 ex9p -pa
22 // mpirun -np 4 ex9p -pa -m ../data/periodic-cube.mesh
23 // mpirun -np 4 ex9p -pa -m ../data/periodic-cube.mesh -d cuda
24 //
25 // Description: This example code solves the time-dependent advection equation
26 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
27 // u0(x)=u(0,x) is a given initial condition.
28 //
29 // The example demonstrates the use of Discontinuous Galerkin (DG)
30 // bilinear forms in MFEM (face integrators), the use of implicit
31 // and explicit ODE time integrators, the definition of periodic
32 // boundary conditions through periodic meshes, as well as the use
33 // of GLVis for persistent visualization of a time-evolving
34 // solution. The saving of time-dependent data files for external
35 // visualization with VisIt (visit.llnl.gov) and ParaView
36 // (paraview.org) is also illustrated.
37 
38 #include "mfem.hpp"
39 #include <fstream>
40 #include <iostream>
41 
42 using namespace std;
43 using namespace mfem;
44 
45 // Choice for the problem setup. The fluid velocity, initial condition and
46 // inflow boundary condition are chosen based on this parameter.
47 int problem;
48 
49 // Velocity coefficient
50 void velocity_function(const Vector &x, Vector &v);
51 
52 // Initial condition
53 double u0_function(const Vector &x);
54 
55 // Inflow boundary condition
56 double inflow_function(const Vector &x);
57 
58 // Mesh bounding box
60 
61 class DG_Solver : public Solver
62 {
63 private:
64  HypreParMatrix &M, &K;
65  SparseMatrix M_diag;
66  HypreParMatrix *A;
67  GMRESSolver linear_solver;
68  BlockILU prec;
69  double dt;
70 public:
71  DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
72  : M(M_),
73  K(K_),
74  A(NULL),
75  linear_solver(M.GetComm()),
76  prec(fes.GetFE(0)->GetDof(),
77  BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
78  dt(-1.0)
79  {
80  linear_solver.iterative_mode = false;
81  linear_solver.SetRelTol(1e-9);
82  linear_solver.SetAbsTol(0.0);
83  linear_solver.SetMaxIter(100);
84  linear_solver.SetPrintLevel(0);
85  linear_solver.SetPreconditioner(prec);
86 
87  M.GetDiag(M_diag);
88  }
89 
90  void SetTimeStep(double dt_)
91  {
92  if (dt_ != dt)
93  {
94  dt = dt_;
95  // Form operator A = M - dt*K
96  delete A;
97  A = Add(-dt, K, 0.0, K);
98  SparseMatrix A_diag;
99  A->GetDiag(A_diag);
100  A_diag.Add(1.0, M_diag);
101  // this will also call SetOperator on the preconditioner
102  linear_solver.SetOperator(*A);
103  }
104  }
105 
106  void SetOperator(const Operator &op)
107  {
108  linear_solver.SetOperator(op);
109  }
110 
111  virtual void Mult(const Vector &x, Vector &y) const
112  {
113  linear_solver.Mult(x, y);
114  }
115 
116  ~DG_Solver()
117  {
118  delete A;
119  }
120 };
121 
122 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
123  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
124  and advection matrices, and b describes the flow on the boundary. This can
125  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
126  used to evaluate the right-hand side. */
127 class FE_Evolution : public TimeDependentOperator
128 {
129 private:
130  OperatorHandle M, K;
131  const Vector &b;
132  Solver *M_prec;
133  CGSolver M_solver;
134  DG_Solver *dg_solver;
135 
136  mutable Vector z;
137 
138 public:
139  FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b);
140 
141  virtual void Mult(const Vector &x, Vector &y) const;
142  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
143 
144  virtual ~FE_Evolution();
145 };
146 
147 
148 int main(int argc, char *argv[])
149 {
150  // 1. Initialize MPI.
151  int num_procs, myid;
152  MPI_Init(&argc, &argv);
153  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
154  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
155 
156  // 2. Parse command-line options.
157  problem = 0;
158  const char *mesh_file = "../data/periodic-hexagon.mesh";
159  int ser_ref_levels = 2;
160  int par_ref_levels = 0;
161  int order = 3;
162  bool pa = false;
163  const char *device_config = "cpu";
164  int ode_solver_type = 4;
165  double t_final = 10.0;
166  double dt = 0.01;
167  bool visualization = true;
168  bool visit = false;
169  bool paraview = false;
170  bool binary = false;
171  int vis_steps = 5;
172 
173  int precision = 8;
174  cout.precision(precision);
175 
176  OptionsParser args(argc, argv);
177  args.AddOption(&mesh_file, "-m", "--mesh",
178  "Mesh file to use.");
179  args.AddOption(&problem, "-p", "--problem",
180  "Problem setup to use. See options in velocity_function().");
181  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
182  "Number of times to refine the mesh uniformly in serial.");
183  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
184  "Number of times to refine the mesh uniformly in parallel.");
185  args.AddOption(&order, "-o", "--order",
186  "Order (degree) of the finite elements.");
187  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
188  "--no-partial-assembly", "Enable Partial Assembly.");
189  args.AddOption(&device_config, "-d", "--device",
190  "Device configuration string, see Device::Configure().");
191  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
192  "ODE solver: 1 - Forward Euler,\n\t"
193  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6,\n\t"
194  " 11 - Backward Euler,\n\t"
195  " 12 - SDIRK23 (L-stable), 13 - SDIRK33,\n\t"
196  " 22 - Implicit Midpoint Method,\n\t"
197  " 23 - SDIRK23 (A-stable), 24 - SDIRK34");
198  args.AddOption(&t_final, "-tf", "--t-final",
199  "Final time; start time is 0.");
200  args.AddOption(&dt, "-dt", "--time-step",
201  "Time step.");
202  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
203  "--no-visualization",
204  "Enable or disable GLVis visualization.");
205  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
206  "--no-visit-datafiles",
207  "Save data files for VisIt (visit.llnl.gov) visualization.");
208  args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
209  "--no-paraview-datafiles",
210  "Save data files for ParaView (paraview.org) visualization.");
211  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
212  "--ascii-datafiles",
213  "Use binary (Sidre) or ascii format for VisIt data files.");
214  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
215  "Visualize every n-th timestep.");
216  args.Parse();
217  if (!args.Good())
218  {
219  if (myid == 0)
220  {
221  args.PrintUsage(cout);
222  }
223  MPI_Finalize();
224  return 1;
225  }
226  if (myid == 0)
227  {
228  args.PrintOptions(cout);
229  }
230 
231  Device device(device_config);
232  if (myid == 0) { device.Print(); }
233 
234  // 3. Read the serial mesh from the given mesh file on all processors. We can
235  // handle geometrically periodic meshes in this code.
236  Mesh *mesh = new Mesh(mesh_file, 1, 1);
237  int dim = mesh->Dimension();
238 
239  // 4. Define the ODE solver used for time integration. Several explicit
240  // Runge-Kutta methods are available.
241  ODESolver *ode_solver = NULL;
242  switch (ode_solver_type)
243  {
244  // Explicit methods
245  case 1: ode_solver = new ForwardEulerSolver; break;
246  case 2: ode_solver = new RK2Solver(1.0); break;
247  case 3: ode_solver = new RK3SSPSolver; break;
248  case 4: ode_solver = new RK4Solver; break;
249  case 6: ode_solver = new RK6Solver; break;
250  // Implicit (L-stable) methods
251  case 11: ode_solver = new BackwardEulerSolver; break;
252  case 12: ode_solver = new SDIRK23Solver(2); break;
253  case 13: ode_solver = new SDIRK33Solver; break;
254  // Implicit A-stable methods (not L-stable)
255  case 22: ode_solver = new ImplicitMidpointSolver; break;
256  case 23: ode_solver = new SDIRK23Solver; break;
257  case 24: ode_solver = new SDIRK34Solver; break;
258  default:
259  if (myid == 0)
260  {
261  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
262  }
263  delete mesh;
264  MPI_Finalize();
265  return 3;
266  }
267 
268  // 5. Refine the mesh in serial to increase the resolution. In this example
269  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
270  // a command-line parameter. If the mesh is of NURBS type, we convert it
271  // to a (piecewise-polynomial) high-order mesh.
272  for (int lev = 0; lev < ser_ref_levels; lev++)
273  {
274  mesh->UniformRefinement();
275  }
276  if (mesh->NURBSext)
277  {
278  mesh->SetCurvature(max(order, 1));
279  }
280  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
281 
282  // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
283  // this mesh further in parallel to increase the resolution. Once the
284  // parallel mesh is defined, the serial mesh can be deleted.
285  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
286  delete mesh;
287  for (int lev = 0; lev < par_ref_levels; lev++)
288  {
289  pmesh->UniformRefinement();
290  }
291 
292  // 7. Define the parallel discontinuous DG finite element space on the
293  // parallel refined mesh of the given polynomial order.
294  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
295  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
296 
297  HYPRE_Int global_vSize = fes->GlobalTrueVSize();
298  if (myid == 0)
299  {
300  cout << "Number of unknowns: " << global_vSize << endl;
301  }
302 
303  // 8. Set up and assemble the parallel bilinear and linear forms (and the
304  // parallel hypre matrices) corresponding to the DG discretization. The
305  // DGTraceIntegrator involves integrals over mesh interior faces.
309 
310  ParBilinearForm *m = new ParBilinearForm(fes);
311  ParBilinearForm *k = new ParBilinearForm(fes);
312  if (pa)
313  {
314  m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
315  k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
316  }
318  k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
320  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
322  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
323 
324  ParLinearForm *b = new ParLinearForm(fes);
326  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
327 
328  int skip_zeros = 0;
329  m->Assemble();
330  k->Assemble(skip_zeros);
331  b->Assemble();
332  m->Finalize();
333  k->Finalize(skip_zeros);
334 
335 
337 
338  // 9. Define the initial conditions, save the corresponding grid function to
339  // a file and (optionally) save data in the VisIt format and initialize
340  // GLVis visualization.
341  ParGridFunction *u = new ParGridFunction(fes);
342  u->ProjectCoefficient(u0);
343  HypreParVector *U = u->GetTrueDofs();
344 
345  {
346  ostringstream mesh_name, sol_name;
347  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
348  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
349  ofstream omesh(mesh_name.str().c_str());
350  omesh.precision(precision);
351  pmesh->Print(omesh);
352  ofstream osol(sol_name.str().c_str());
353  osol.precision(precision);
354  u->Save(osol);
355  }
356 
357  // Create data collection for solution output: either VisItDataCollection for
358  // ascii data files, or SidreDataCollection for binary data files.
359  DataCollection *dc = NULL;
360  if (visit)
361  {
362  if (binary)
363  {
364 #ifdef MFEM_USE_SIDRE
365  dc = new SidreDataCollection("Example9-Parallel", pmesh);
366 #else
367  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
368 #endif
369  }
370  else
371  {
372  dc = new VisItDataCollection("Example9-Parallel", pmesh);
373  dc->SetPrecision(precision);
374  // To save the mesh using MFEM's parallel mesh format:
375  // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
376  }
377  dc->RegisterField("solution", u);
378  dc->SetCycle(0);
379  dc->SetTime(0.0);
380  dc->Save();
381  }
382 
383  ParaViewDataCollection *pd = NULL;
384  if (paraview)
385  {
386  pd = new ParaViewDataCollection("Example9P", pmesh);
387  pd->SetPrefixPath("ParaView");
388  pd->RegisterField("solution", u);
389  pd->SetLevelsOfDetail(order);
390  pd->SetDataFormat(VTKFormat::BINARY);
391  pd->SetHighOrderOutput(true);
392  pd->SetCycle(0);
393  pd->SetTime(0.0);
394  pd->Save();
395  }
396 
397  socketstream sout;
398  if (visualization)
399  {
400  char vishost[] = "localhost";
401  int visport = 19916;
402  sout.open(vishost, visport);
403  if (!sout)
404  {
405  if (myid == 0)
406  cout << "Unable to connect to GLVis server at "
407  << vishost << ':' << visport << endl;
408  visualization = false;
409  if (myid == 0)
410  {
411  cout << "GLVis visualization disabled.\n";
412  }
413  }
414  else
415  {
416  sout << "parallel " << num_procs << " " << myid << "\n";
417  sout.precision(precision);
418  sout << "solution\n" << *pmesh << *u;
419  sout << "pause\n";
420  sout << flush;
421  if (myid == 0)
422  cout << "GLVis visualization paused."
423  << " Press space (in the GLVis window) to resume it.\n";
424  }
425  }
426 
427  // 10. Define the time-dependent evolution operator describing the ODE
428  // right-hand side, and perform time-integration (looping over the time
429  // iterations, ti, with a time-step dt).
430  FE_Evolution adv(*m, *k, *B);
431 
432  double t = 0.0;
433  adv.SetTime(t);
434  ode_solver->Init(adv);
435 
436  bool done = false;
437  for (int ti = 0; !done; )
438  {
439  double dt_real = min(dt, t_final - t);
440  ode_solver->Step(*U, t, dt_real);
441  ti++;
442 
443  done = (t >= t_final - 1e-8*dt);
444 
445  if (done || ti % vis_steps == 0)
446  {
447  if (myid == 0)
448  {
449  cout << "time step: " << ti << ", time: " << t << endl;
450  }
451 
452  // 11. Extract the parallel grid function corresponding to the finite
453  // element approximation U (the local solution on each processor).
454  *u = *U;
455 
456  if (visualization)
457  {
458  sout << "parallel " << num_procs << " " << myid << "\n";
459  sout << "solution\n" << *pmesh << *u << flush;
460  }
461 
462  if (visit)
463  {
464  dc->SetCycle(ti);
465  dc->SetTime(t);
466  dc->Save();
467  }
468 
469  if (paraview)
470  {
471  pd->SetCycle(ti);
472  pd->SetTime(t);
473  pd->Save();
474  }
475  }
476  }
477 
478  // 12. Save the final solution in parallel. This output can be viewed later
479  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
480  {
481  *u = *U;
482  ostringstream sol_name;
483  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
484  ofstream osol(sol_name.str().c_str());
485  osol.precision(precision);
486  u->Save(osol);
487  }
488 
489  // 13. Free the used memory.
490  delete U;
491  delete u;
492  delete B;
493  delete b;
494  delete k;
495  delete m;
496  delete fes;
497  delete pmesh;
498  delete ode_solver;
499  delete pd;
500  delete dc;
501 
502  MPI_Finalize();
503  return 0;
504 }
505 
506 
507 // Implementation of class FE_Evolution
509  const Vector &_b)
510  : TimeDependentOperator(_M.Height()),
511  b(_b),
512  M_solver(_M.ParFESpace()->GetComm()),
513  z(_M.Height())
514 {
515  bool pa = _M.GetAssemblyLevel()==AssemblyLevel::PARTIAL;
516 
517  if (pa)
518  {
519  M.Reset(&_M, false);
520  K.Reset(&_K, false);
521  }
522  else
523  {
524  M.Reset(_M.ParallelAssemble(), true);
525  K.Reset(_K.ParallelAssemble(), true);
526  }
527 
528  M_solver.SetOperator(*M);
529 
530  Array<int> ess_tdof_list;
531  if (pa)
532  {
533  M_prec = new OperatorJacobiSmoother(_M, ess_tdof_list);
534  dg_solver = NULL;
535  }
536  else
537  {
538  HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
539  HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
540  HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
541  M_prec = hypre_prec;
542 
543  dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace());
544  }
545 
546  M_solver.SetPreconditioner(*M_prec);
547  M_solver.iterative_mode = false;
548  M_solver.SetRelTol(1e-9);
549  M_solver.SetAbsTol(0.0);
550  M_solver.SetMaxIter(100);
551  M_solver.SetPrintLevel(0);
552 }
553 
554 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
555 {
556  K->Mult(x, z);
557  z += b;
558  dg_solver->SetTimeStep(dt);
559  dg_solver->Mult(z, k);
560 }
561 
562 void FE_Evolution::Mult(const Vector &x, Vector &y) const
563 {
564  // y = M^{-1} (K x + b)
565  K->Mult(x, z);
566  z += b;
567  M_solver.Mult(z, y);
568 }
569 
571 {
572  delete M_prec;
573  delete dg_solver;
574 }
575 
576 
577 // Velocity coefficient
578 void velocity_function(const Vector &x, Vector &v)
579 {
580  int dim = x.Size();
581 
582  // map to the reference [-1,1] domain
583  Vector X(dim);
584  for (int i = 0; i < dim; i++)
585  {
586  double center = (bb_min[i] + bb_max[i]) * 0.5;
587  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
588  }
589 
590  switch (problem)
591  {
592  case 0:
593  {
594  // Translations in 1D, 2D, and 3D
595  switch (dim)
596  {
597  case 1: v(0) = 1.0; break;
598  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
599  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
600  break;
601  }
602  break;
603  }
604  case 1:
605  case 2:
606  {
607  // Clockwise rotation in 2D around the origin
608  const double w = M_PI/2;
609  switch (dim)
610  {
611  case 1: v(0) = 1.0; break;
612  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
613  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
614  }
615  break;
616  }
617  case 3:
618  {
619  // Clockwise twisting rotation in 2D around the origin
620  const double w = M_PI/2;
621  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
622  d = d*d;
623  switch (dim)
624  {
625  case 1: v(0) = 1.0; break;
626  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
627  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
628  }
629  break;
630  }
631  }
632 }
633 
634 // Initial condition
635 double u0_function(const Vector &x)
636 {
637  int dim = x.Size();
638 
639  // map to the reference [-1,1] domain
640  Vector X(dim);
641  for (int i = 0; i < dim; i++)
642  {
643  double center = (bb_min[i] + bb_max[i]) * 0.5;
644  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
645  }
646 
647  switch (problem)
648  {
649  case 0:
650  case 1:
651  {
652  switch (dim)
653  {
654  case 1:
655  return exp(-40.*pow(X(0)-0.5,2));
656  case 2:
657  case 3:
658  {
659  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
660  if (dim == 3)
661  {
662  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
663  rx *= s;
664  ry *= s;
665  }
666  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
667  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
668  }
669  }
670  }
671  case 2:
672  {
673  double x_ = X(0), y_ = X(1), rho, phi;
674  rho = hypot(x_, y_);
675  phi = atan2(y_, x_);
676  return pow(sin(M_PI*rho),2)*sin(3*phi);
677  }
678  case 3:
679  {
680  const double f = M_PI;
681  return sin(f*X(0))*sin(f*X(1));
682  }
683  }
684  return 0.0;
685 }
686 
687 // Inflow boundary condition (zero for the problems considered in this example)
688 double inflow_function(const Vector &x)
689 {
690  switch (problem)
691  {
692  case 0:
693  case 1:
694  case 2:
695  case 3: return 0.0;
696  }
697  return 0.0;
698 }
Vector bb_max
Definition: ex9.cpp:60
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:2281
Conjugate gradient method.
Definition: solvers.hpp:152
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:358
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:259
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:120
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
virtual void RegisterField(const std::string &field_name, mfem::GridFunction *gf) override
Add a grid function to the collection.
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:303
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
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:485
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:236
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:301
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:504
int main(int argc, char *argv[])
Definition: ex1.cpp:62
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:360
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:104
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:1926
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:70
Data type sparse matrix.
Definition: sparsemat.hpp:40
virtual void Save()
Save the collection to disk.
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:84
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:65
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4102
Parallel smoothers in hypre.
Definition: hypre.hpp:581
void SetHighOrderOutput(bool high_order_output_)
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:60
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:73
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:145
void SetAbsTol(double atol)
Definition: solvers.hpp:64
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:461
AssemblyLevel GetAssemblyLevel()
Get the assembly level.
void SetRelTol(double rtol)
Definition: solvers.hpp:63
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:478
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:314
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)
Definition: optparser.hpp:76
int problem
Definition: ex15.cpp:54
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:132
GMRES method.
Definition: solvers.hpp:184
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:130
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:373
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:489
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:43
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
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)
double u0_function(const Vector &x)
Definition: ex9.cpp:535
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:166
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
virtual void Save() override
Base class for solvers.
Definition: operator.hpp:500
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:114
The classical forward Euler method.
Definition: ode.hpp:99
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:187
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:34
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
double inflow_function(const Vector &x)
Definition: ex9.cpp:588
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:143
bool Good() const
Definition: optparser.hpp:125
alpha (q . grad u, v)