MFEM  v4.5.1
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 // PETSc Modification
3 //
4 // Compile with: make ex9p
5 //
6 // Sample runs:
7 // mpirun -np 4 ex9p -m ../../data/periodic-hexagon.mesh --petscopts rc_ex9p_expl
8 // mpirun -np 4 ex9p -m ../../data/periodic-hexagon.mesh --petscopts rc_ex9p_impl -implicit
9 //
10 // Description: This example code solves the time-dependent advection equation
11 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
12 // u0(x)=u(0,x) is a given initial condition.
13 //
14 // The example demonstrates the use of Discontinuous Galerkin (DG)
15 // bilinear forms in MFEM (face integrators), the use of explicit
16 // ODE time integrators, the definition of periodic boundary
17 // conditions through periodic meshes, as well as the use of GLVis
18 // for persistent visualization of a time-evolving solution. The
19 // saving of time-dependent data files for external visualization
20 // with VisIt (visit.llnl.gov) is also illustrated.
21 //
22 // The example also demonstrates how to use PETSc ODE solvers and
23 // customize them by command line (see rc_ex9p_expl and
24 // rc_ex9p_impl). The split in left-hand side and right-hand side
25 // of the TimeDependentOperator is amenable for IMEX methods.
26 // When using fully implicit methods, just the left-hand side of
27 // the operator should be provided for efficiency reasons when
28 // assembling the Jacobians. Here, we provide two Jacobian
29 // routines just to illustrate the capabilities of the
30 // PetscODESolver class. We also show how to monitor the time
31 // dependent solution inside a call to PetscODESolver:Mult.
32 
33 #include "mfem.hpp"
34 #include <fstream>
35 #include <iostream>
36 
37 #ifndef MFEM_USE_PETSC
38 #error This example requires that MFEM is built with MFEM_USE_PETSC=YES
39 #endif
40 
41 using namespace std;
42 using namespace mfem;
43 
44 // Choice for the problem setup. The fluid velocity, initial condition and
45 // inflow boundary condition are chosen based on this parameter.
46 int problem;
47 
48 // Velocity coefficient
49 void velocity_function(const Vector &x, Vector &v);
50 
51 // Initial condition
52 double u0_function(const Vector &x);
53 
54 // Inflow boundary condition
55 double inflow_function(const Vector &x);
56 
57 // Mesh bounding box
59 
60 
61 /** A time-dependent operator for the ODE as F(u,du/dt,t) = G(u,t)
62  The DG weak form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
63  and advection matrices, and b describes the flow on the boundary. This can
64  be also written as a general ODE with the right-hand side only as
65  du/dt = M^{-1} (K u + b).
66  This class is used to evaluate the right-hand side and the left-hand side. */
68 {
69 private:
70  OperatorHandle M, K;
71  const Vector &b;
72  MPI_Comm comm;
73  Solver *M_prec;
74  CGSolver M_solver;
75  AssemblyLevel MAlev,KAlev;
76 
77  mutable Vector z;
78  mutable PetscParMatrix* iJacobian;
79  mutable PetscParMatrix* rJacobian;
80 
81 public:
83  bool implicit);
84 
85  virtual void ExplicitMult(const Vector &x, Vector &y) const;
86  virtual void ImplicitMult(const Vector &x, const Vector &xp, Vector &y) const;
87  virtual void Mult(const Vector &x, Vector &y) const;
88  virtual Operator& GetExplicitGradient(const Vector &x) const;
89  virtual Operator& GetImplicitGradient(const Vector &x, const Vector &xp,
90  double shift) const;
91  virtual ~FE_Evolution() { delete iJacobian; delete rJacobian; }
92 };
93 
94 
95 // Monitor the solution at time step "step", explicitly in the time loop
96 class UserMonitor : public PetscSolverMonitor
97 {
98 private:
99  socketstream& sout;
100  ParMesh* pmesh;
102  int vt;
103  bool pause;
104 
105 public:
106  UserMonitor(socketstream& s_, ParMesh* m_, ParGridFunction* u_, int vt_) :
107  PetscSolverMonitor(true,false), sout(s_), pmesh(m_), u(u_), vt(vt_),
108  pause(true) {}
109 
110  void MonitorSolution(PetscInt step, PetscReal norm, const Vector &X)
111  {
112  if (step % vt == 0)
113  {
114  int num_procs, myid;
115 
116  *u = X;
117  MPI_Comm_size(pmesh->GetComm(),&num_procs);
118  MPI_Comm_rank(pmesh->GetComm(),&myid);
119  sout << "parallel " << num_procs << " " << myid << "\n";
120  sout << "solution\n" << *pmesh << *u;
121  if (pause) { sout << "pause\n"; }
122  sout << flush;
123  if (pause)
124  {
125  pause = false;
126  if (myid == 0)
127  {
128  cout << "GLVis visualization paused."
129  << " Press space (in the GLVis window) to resume it.\n";
130  }
131  }
132  }
133  }
134 };
135 
136 
137 int main(int argc, char *argv[])
138 {
139  // 1. Initialize MPI and HYPRE.
140  Mpi::Init(argc, argv);
141  int num_procs = Mpi::WorldSize();
142  int myid = Mpi::WorldRank();
143  Hypre::Init();
144 
145  // 2. Parse command-line options.
146  problem = 0;
147  const char *mesh_file = "../../data/periodic-hexagon.mesh";
148  int ser_ref_levels = 2;
149  int par_ref_levels = 0;
150  int order = 3;
151  bool pa = false;
152  bool ea = false;
153  bool fa = false;
154  const char *device_config = "cpu";
155  int ode_solver_type = 4;
156  double t_final = 10.0;
157  double dt = 0.01;
158  bool visualization = true;
159  bool visit = false;
160  bool binary = false;
161  int vis_steps = 5;
162  bool use_petsc = true;
163  bool implicit = false;
164  bool use_step = true;
165  const char *petscrc_file = "";
166 
167  int precision = 8;
168  cout.precision(precision);
169 
170  OptionsParser args(argc, argv);
171  args.AddOption(&mesh_file, "-m", "--mesh",
172  "Mesh file to use.");
173  args.AddOption(&problem, "-p", "--problem",
174  "Problem setup to use. See options in velocity_function().");
175  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
176  "Number of times to refine the mesh uniformly in serial.");
177  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
178  "Number of times to refine the mesh uniformly in parallel.");
179  args.AddOption(&order, "-o", "--order",
180  "Order (degree) of the finite elements.");
181  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
182  "--no-partial-assembly", "Enable Partial Assembly.");
183  args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
184  "--no-element-assembly", "Enable Element Assembly.");
185  args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
186  "--no-full-assembly", "Enable Full Assembly.");
187  args.AddOption(&device_config, "-d", "--device",
188  "Device configuration string, see Device::Configure().");
189  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
190  "ODE solver: 1 - Forward Euler,\n\t"
191  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
192  args.AddOption(&t_final, "-tf", "--t-final",
193  "Final time; start time is 0.");
194  args.AddOption(&dt, "-dt", "--time-step",
195  "Time step.");
196  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
197  "--no-visualization",
198  "Enable or disable GLVis visualization.");
199  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
200  "--no-visit-datafiles",
201  "Save data files for VisIt (visit.llnl.gov) visualization.");
202  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
203  "--ascii-datafiles",
204  "Use binary (Sidre) or ascii format for VisIt data files.");
205  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
206  "Visualize every n-th timestep.");
207  args.AddOption(&use_petsc, "-usepetsc", "--usepetsc", "-no-petsc",
208  "--no-petsc",
209  "Use or not PETSc to solve the ODE system.");
210  args.AddOption(&petscrc_file, "-petscopts", "--petscopts",
211  "PetscOptions file to use.");
212  args.AddOption(&use_step, "-usestep", "--usestep", "-no-step",
213  "--no-step",
214  "Use the Step() or Run() method to solve the ODE system.");
215  args.AddOption(&implicit, "-implicit", "--implicit", "-no-implicit",
216  "--no-implicit",
217  "Use or not an implicit method in PETSc to solve the ODE system.");
218  args.Parse();
219  if (!args.Good())
220  {
221  if (myid == 0)
222  {
223  args.PrintUsage(cout);
224  }
225  return 1;
226  }
227  if (myid == 0)
228  {
229  args.PrintOptions(cout);
230  }
231 
232  Device device(device_config);
233  if (myid == 0) { device.Print(); }
234 
235  // 3. Read the serial mesh from the given mesh file on all processors. We can
236  // handle geometrically periodic meshes in this code.
237  Mesh *mesh = new Mesh(mesh_file, 1, 1);
238  int dim = mesh->Dimension();
239 
240  // 4. Define the ODE solver used for time integration. Several explicit
241  // Runge-Kutta methods are available.
242  ODESolver *ode_solver = NULL;
243  PetscODESolver *pode_solver = NULL;
244  UserMonitor *pmon = NULL;
245  if (!use_petsc)
246  {
247  switch (ode_solver_type)
248  {
249  case 1: ode_solver = new ForwardEulerSolver; break;
250  case 2: ode_solver = new RK2Solver(1.0); break;
251  case 3: ode_solver = new RK3SSPSolver; break;
252  case 4: ode_solver = new RK4Solver; break;
253  case 6: ode_solver = new RK6Solver; break;
254  default:
255  if (myid == 0)
256  {
257  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
258  }
259  return 3;
260  }
261  }
262  else
263  {
264  // When using PETSc, we just create the ODE solver. We use command line
265  // customization to select a specific solver.
266  MFEMInitializePetsc(NULL, NULL, petscrc_file, NULL);
267  ode_solver = pode_solver = new PetscODESolver(MPI_COMM_WORLD);
268  }
269 
270  // 5. Refine the mesh in serial to increase the resolution. In this example
271  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
272  // a command-line parameter. If the mesh is of NURBS type, we convert it
273  // to a (piecewise-polynomial) high-order mesh.
274  for (int lev = 0; lev < ser_ref_levels; lev++)
275  {
276  mesh->UniformRefinement();
277  }
278  if (mesh->NURBSext)
279  {
280  mesh->SetCurvature(max(order, 1));
281  }
282  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
283 
284  // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
285  // this mesh further in parallel to increase the resolution. Once the
286  // parallel mesh is defined, the serial mesh can be deleted.
287  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
288  delete mesh;
289  for (int lev = 0; lev < par_ref_levels; lev++)
290  {
291  pmesh->UniformRefinement();
292  }
293 
294  // 7. Define the parallel discontinuous DG finite element space on the
295  // parallel refined mesh of the given polynomial order.
296  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
297  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
298 
299  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
300  if (myid == 0)
301  {
302  cout << "Number of unknowns: " << global_vSize << endl;
303  }
304 
305  // 8. Set up and assemble the parallel bilinear and linear forms (and the
306  // parallel hypre matrices) corresponding to the DG discretization. The
307  // DGTraceIntegrator involves integrals over mesh interior faces.
311 
312  ParBilinearForm *m = new ParBilinearForm(fes);
313  ParBilinearForm *k = new ParBilinearForm(fes);
314  if (pa)
315  {
316  m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
317  k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
318  }
319  else if (ea)
320  {
321  m->SetAssemblyLevel(AssemblyLevel::ELEMENT);
322  k->SetAssemblyLevel(AssemblyLevel::ELEMENT);
323  }
324  else if (fa)
325  {
326  m->SetAssemblyLevel(AssemblyLevel::FULL);
327  k->SetAssemblyLevel(AssemblyLevel::FULL);
328  }
329 
331  k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
333  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
335  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
336 
337  ParLinearForm *b = new ParLinearForm(fes);
339  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
340 
341  int skip_zeros = 0;
342  m->Assemble();
343  k->Assemble(skip_zeros);
344  b->Assemble();
345  m->Finalize();
346  k->Finalize(skip_zeros);
347 
348 
350 
351  // 9. Define the initial conditions, save the corresponding grid function to
352  // a file and (optionally) save data in the VisIt format and initialize
353  // GLVis visualization.
354  ParGridFunction *u = new ParGridFunction(fes);
355  u->ProjectCoefficient(u0);
356  HypreParVector *U = u->GetTrueDofs();
357 
358  {
359  ostringstream mesh_name, sol_name;
360  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
361  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
362  ofstream omesh(mesh_name.str().c_str());
363  omesh.precision(precision);
364  pmesh->Print(omesh);
365  ofstream osol(sol_name.str().c_str());
366  osol.precision(precision);
367  u->Save(osol);
368  }
369 
370  // Create data collection for solution output: either VisItDataCollection for
371  // ascii data files, or SidreDataCollection for binary data files.
372  DataCollection *dc = NULL;
373  if (visit)
374  {
375  if (binary)
376  {
377 #ifdef MFEM_USE_SIDRE
378  dc = new SidreDataCollection("Example9-Parallel", pmesh);
379 #else
380  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
381 #endif
382  }
383  else
384  {
385  dc = new VisItDataCollection("Example9-Parallel", pmesh);
386  dc->SetPrecision(precision);
387  }
388  dc->RegisterField("solution", u);
389  dc->SetCycle(0);
390  dc->SetTime(0.0);
391  dc->Save();
392  }
393 
394  socketstream sout;
395  if (visualization)
396  {
397  char vishost[] = "localhost";
398  int visport = 19916;
399  sout.open(vishost, visport);
400  if (!sout)
401  {
402  if (myid == 0)
403  cout << "Unable to connect to GLVis server at "
404  << vishost << ':' << visport << endl;
405  visualization = false;
406  if (myid == 0)
407  {
408  cout << "GLVis visualization disabled.\n";
409  }
410  }
411  else if (use_step)
412  {
413  sout << "parallel " << num_procs << " " << myid << "\n";
414  sout.precision(precision);
415  sout << "solution\n" << *pmesh << *u;
416  sout << "pause\n";
417  sout << flush;
418  if (myid == 0)
419  cout << "GLVis visualization paused."
420  << " Press space (in the GLVis window) to resume it.\n";
421  }
422  else if (use_petsc)
423  {
424  // Set the monitoring routine for the PetscODESolver.
425  sout.precision(precision);
426  pmon = new UserMonitor(sout,pmesh,u,vis_steps);
427  pode_solver->SetMonitor(pmon);
428  }
429  }
430 
431  // 10. Define the time-dependent evolution operator describing the ODE
432  FE_Evolution *adv = new FE_Evolution(*m, *k, *B, implicit);
433 
434  double t = 0.0;
435  adv->SetTime(t);
436  if (use_petsc)
437  {
438  pode_solver->Init(*adv,PetscODESolver::ODE_SOLVER_LINEAR);
439  }
440  else
441  {
442  ode_solver->Init(*adv);
443  }
444 
445  // Explicitly perform time-integration (looping over the time iterations, ti,
446  // with a time-step dt), or use the Run method of the ODE solver class.
447  if (use_step)
448  {
449  bool done = false;
450  for (int ti = 0; !done; )
451  {
452  // We cannot match exactly the time history of the Run method
453  // since we are explicitly telling PETSc to use a time step
454  double dt_real = min(dt, t_final - t);
455  ode_solver->Step(*U, t, dt_real);
456  ti++;
457 
458  done = (t >= t_final - 1e-8*dt);
459 
460  if (done || ti % vis_steps == 0)
461  {
462  if (myid == 0)
463  {
464  cout << "time step: " << ti << ", time: " << t << endl;
465  }
466  // 11. Extract the parallel grid function corresponding to the finite
467  // element approximation U (the local solution on each processor).
468  *u = *U;
469 
470  if (visualization)
471  {
472  sout << "parallel " << num_procs << " " << myid << "\n";
473  sout << "solution\n" << *pmesh << *u << flush;
474  }
475 
476  if (visit)
477  {
478  dc->SetCycle(ti);
479  dc->SetTime(t);
480  dc->Save();
481  }
482  }
483  }
484  }
485  else { ode_solver->Run(*U, t, dt, t_final); }
486 
487  // 12. Save the final solution in parallel. This output can be viewed later
488  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
489  {
490  *u = *U;
491  ostringstream sol_name;
492  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
493  ofstream osol(sol_name.str().c_str());
494  osol.precision(precision);
495  u->Save(osol);
496  }
497 
498  // 13. Free the used memory.
499  delete U;
500  delete u;
501  delete B;
502  delete b;
503  delete k;
504  delete m;
505  delete fes;
506  delete pmesh;
507  delete ode_solver;
508  delete dc;
509  delete adv;
510 
511  delete pmon;
512 
513  // We finalize PETSc
514  if (use_petsc) { MFEMFinalizePetsc(); }
515 
516  return 0;
517 }
518 
519 
520 // Implementation of class FE_Evolution
522  const Vector &b_,bool M_in_lhs)
523  : TimeDependentOperator(M_.Height(), 0.0,
524  M_in_lhs ? TimeDependentOperator::IMPLICIT
525  : TimeDependentOperator::EXPLICIT),
526  b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(M_.Height()),
527  iJacobian(NULL), rJacobian(NULL)
528 {
529  MAlev = M_.GetAssemblyLevel();
530  KAlev = K_.GetAssemblyLevel();
531  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
532  {
533  M.Reset(M_.ParallelAssemble(), true);
534  K.Reset(K_.ParallelAssemble(), true);
535  }
536  else
537  {
538  M.Reset(&M_, false);
539  K.Reset(&K_, false);
540  }
541 
542  M_solver.SetOperator(*M);
543 
544  Array<int> ess_tdof_list;
545  if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
546  {
547  HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
548 
549  HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
550  M_prec = hypre_prec;
551  }
552  else
553  {
554  M_prec = new OperatorJacobiSmoother(M_, ess_tdof_list);
555  }
556 
557  M_solver.SetPreconditioner(*M_prec);
558  M_solver.iterative_mode = false;
559  M_solver.SetRelTol(1e-9);
560  M_solver.SetAbsTol(0.0);
561  M_solver.SetMaxIter(100);
562  M_solver.SetPrintLevel(0);
563 }
564 
565 // RHS evaluation
566 void FE_Evolution::ExplicitMult(const Vector &x, Vector &y) const
567 {
568  if (isExplicit())
569  {
570  // y = M^{-1} (K x + b)
571  K->Mult(x, z);
572  z += b;
573  M_solver.Mult(z, y);
574  }
575  else
576  {
577  // y = K x + b
578  K->Mult(x, y);
579  y += b;
580  }
581 }
582 
583 // LHS evaluation
584 void FE_Evolution::ImplicitMult(const Vector &x, const Vector &xp,
585  Vector &y) const
586 {
587  if (isImplicit())
588  {
589  M->Mult(xp, y);
590  }
591  else
592  {
593  y = xp;
594  }
595 }
596 
597 void FE_Evolution::Mult(const Vector &x, Vector &y) const
598 {
599  // y = M^{-1} (K x + b)
600  K->Mult(x, z);
601  z += b;
602  M_solver.Mult(z, y);
603 }
604 
605 // RHS Jacobian
607 {
608  delete rJacobian;
609  Operator::Type otype = (KAlev == AssemblyLevel::LEGACY ?
610  Operator::PETSC_MATAIJ : Operator::ANY_TYPE);
611  if (isImplicit())
612  {
613  rJacobian = new PetscParMatrix(comm, K.Ptr(), otype);
614  }
615  else
616  {
617  mfem_error("FE_Evolution::GetExplicitGradient(x): Capability not coded!");
618  }
619  return *rJacobian;
620 }
621 
622 // LHS Jacobian, evaluated as shift*F_du/dt + F_u
624  double shift) const
625 {
626  Operator::Type otype = (MAlev == AssemblyLevel::LEGACY ?
627  Operator::PETSC_MATAIJ : Operator::ANY_TYPE);
628  delete iJacobian;
629  if (isImplicit())
630  {
631  iJacobian = new PetscParMatrix(comm, M.Ptr(), otype);
632  *iJacobian *= shift;
633  }
634  else
635  {
636  mfem_error("FE_Evolution::GetImplicitGradient(x,xp,shift):"
637  " Capability not coded!");
638  }
639  return *iJacobian;
640 }
641 
642 // Velocity coefficient
643 void velocity_function(const Vector &x, Vector &v)
644 {
645  int dim = x.Size();
646 
647  // map to the reference [-1,1] domain
648  Vector X(dim);
649  for (int i = 0; i < dim; i++)
650  {
651  double center = (bb_min[i] + bb_max[i]) * 0.5;
652  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
653  }
654 
655  switch (problem)
656  {
657  case 0:
658  {
659  // Translations in 1D, 2D, and 3D
660  switch (dim)
661  {
662  case 1: v(0) = 1.0; break;
663  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
664  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
665  break;
666  }
667  break;
668  }
669  case 1:
670  case 2:
671  {
672  // Clockwise rotation in 2D around the origin
673  const double w = M_PI/2;
674  switch (dim)
675  {
676  case 1: v(0) = 1.0; break;
677  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
678  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
679  }
680  break;
681  }
682  case 3:
683  {
684  // Clockwise twisting rotation in 2D around the origin
685  const double w = M_PI/2;
686  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
687  d = d*d;
688  switch (dim)
689  {
690  case 1: v(0) = 1.0; break;
691  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
692  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
693  }
694  break;
695  }
696  }
697 }
698 
699 // Initial condition
700 double u0_function(const Vector &x)
701 {
702  int dim = x.Size();
703 
704  // map to the reference [-1,1] domain
705  Vector X(dim);
706  for (int i = 0; i < dim; i++)
707  {
708  double center = (bb_min[i] + bb_max[i]) * 0.5;
709  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
710  }
711 
712  switch (problem)
713  {
714  case 0:
715  case 1:
716  {
717  switch (dim)
718  {
719  case 1:
720  return exp(-40.*pow(X(0)-0.5,2));
721  case 2:
722  case 3:
723  {
724  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
725  if (dim == 3)
726  {
727  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
728  rx *= s;
729  ry *= s;
730  }
731  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
732  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
733  }
734  }
735  }
736  case 2:
737  {
738  double x_ = X(0), y_ = X(1), rho, phi;
739  rho = hypot(x_, y_);
740  phi = atan2(y_, x_);
741  return pow(sin(M_PI*rho),2)*sin(3*phi);
742  }
743  case 3:
744  {
745  const double f = M_PI;
746  return sin(f*X(0))*sin(f*X(1));
747  }
748  }
749  return 0.0;
750 }
751 
752 // Inflow boundary condition (zero for the problems considered in this example)
753 double inflow_function(const Vector &x)
754 {
755  switch (problem)
756  {
757  case 0:
758  case 1:
759  case 2:
760  case 3: return 0.0;
761  }
762  return 0.0;
763 }
virtual void ImplicitMult(const Vector &x, const Vector &xp, Vector &y) const
Perform the action of the implicit part of the operator, F: y = F(x, k, t) where t is the current tim...
Definition: ex9p.cpp:584
Vector bb_max
Definition: ex9.cpp:67
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
Conjugate gradient method.
Definition: solvers.hpp:465
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition: operator.hpp:338
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
HYPRE_Int PetscInt
Definition: petsc.hpp:47
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
Wrapper for PETSc&#39;s matrix class.
Definition: petsc.hpp:315
Base abstract class for first order time dependent operators.
Definition: operator.hpp:289
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:475
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
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:873
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:333
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:525
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:659
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
Definition: petsc.cpp:2471
Abstract class for PETSc&#39;s ODE solvers.
Definition: petsc.hpp:938
double PetscReal
Definition: petsc.hpp:49
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.
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void Init(TimeDependentOperator &f_, enum PetscODESolver::Type type)
Initialize the ODE solver.
Definition: petsc.cpp:4103
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
Abstract class for monitoring PETSc&#39;s solvers.
Definition: petsc.hpp:973
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
virtual void Save()
Save the collection to disk.
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
constexpr int visport
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5343
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Parallel smoothers in hypre.
Definition: hypre.hpp:962
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:85
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:336
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:67
A general vector function coefficient.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:256
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
void SetAbsTol(double atol)
Definition: solvers.hpp:199
void SetRelTol(double rtol)
Definition: solvers.hpp:198
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:501
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
virtual Operator & GetExplicitGradient(const Vector &x) const
Return an Operator representing dG/dx at the given point x and the currently set time.
Definition: ex9p.cpp:606
virtual ~FE_Evolution()
Definition: ex9p.cpp:91
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
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
Definition: ex18.hpp:81
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
void MFEMFinalizePetsc()
Definition: petsc.cpp:192
virtual Operator & GetImplicitGradient(const Vector &x, const Vector &xp, double shift) const
Return an Operator representing (dF/dk shift + dF/dx) at the given x, k, and the currently set time...
Definition: ex9p.cpp:623
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
Definition: petsc.cpp:168
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.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
Class for parallel bilinear form.
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:558
RefCoord t[3]
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
virtual void Run(Vector &x, double &t, double &dt, double tf)
Perform time integration from time t [in] to time tf [in].
Definition: ode.hpp:89
RefCoord s[3]
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Base class for solvers.
Definition: operator.hpp:655
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
The classical forward Euler method.
Definition: ode.hpp:116
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void ExplicitMult(const Vector &x, Vector &y) const
Perform the action of the explicit part of the operator, G: y = G(x, t) where t is the current time...
Definition: ex9p.cpp:566
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:95
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int main()
double inflow_function(const Vector &x)
Definition: ex9.cpp:611
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:288
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)