MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
41using namespace std;
42using 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.
47
48// Velocity coefficient
49void velocity_function(const Vector &x, Vector &v);
50
51// Initial condition
52real_t u0_function(const Vector &x);
53
54// Inflow boundary condition
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. */
67class FE_Evolution : public TimeDependentOperator
68{
69private:
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
81public:
82 FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, const Vector &b_,
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 real_t 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
96class UserMonitor : public PetscSolverMonitor
97{
98private:
99 socketstream& sout;
100 ParMesh* pmesh;
102 int vt;
103 bool pause;
104
105public:
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
137int 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 real_t t_final = 10.0;
157 real_t 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.
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);
338 b->AddBdrFaceIntegrator(
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
349 HypreParVector *B = b->ParallelAssemble();
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.
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 real_t 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 real_t 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
521FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
522 const Vector &b_,bool M_in_lhs)
523 : TimeDependentOperator(M_.ParFESpace()->GetTrueVSize(), 0.0,
524 M_in_lhs ? TimeDependentOperator::IMPLICIT
525 : TimeDependentOperator::EXPLICIT),
526 b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(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
566void 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
584void 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
597void 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
606Operator& FE_Evolution::GetExplicitGradient(const Vector &x) const
607{
608 delete rJacobian;
609 Operator::Type otype = (KAlev == AssemblyLevel::LEGACY ?
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
623Operator& FE_Evolution::GetImplicitGradient(const Vector &x, const Vector &xp,
624 real_t shift) const
625{
626 Operator::Type otype = (MAlev == AssemblyLevel::LEGACY ?
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
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 real_t 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 real_t 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 real_t w = M_PI/2;
686 real_t 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
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 real_t 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 real_t rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
725 if (dim == 3)
726 {
727 const real_t 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 real_t 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 real_t 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)
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}
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
AssemblyLevel GetAssemblyLevel() const
Returns the assembly level.
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
virtual void Save()
Save the collection to disk.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
The classical forward Euler method.
Definition ode.hpp:118
A general function coefficient.
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:206
Parallel smoothers in hypre.
Definition hypre.hpp:1020
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition mesh.cpp:138
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition ode.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
virtual void Run(Vector &x, real_t &t, real_t &dt, real_t tf)
Perform time integration from time t [in] to time tf [in].
Definition ode.hpp:90
Pointer to an Operator of a specified type.
Definition handle.hpp:34
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
Abstract operator.
Definition operator.hpp:25
Type
Enumeration defining IDs for some classes derived from Operator.
Definition operator.hpp:284
@ ANY_TYPE
ID for the base class Operator, i.e. any type.
Definition operator.hpp:285
@ PETSC_MATAIJ
ID for class PetscParMatrix, MATAIJ format.
Definition operator.hpp:288
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Abstract class for PETSc's ODE solvers.
Definition petsc.hpp:950
virtual void Init(TimeDependentOperator &f_, enum PetscODESolver::Type type)
Initialize the ODE solver.
Definition petsc.cpp:4173
Wrapper for PETSc's matrix class.
Definition petsc.hpp:319
Abstract class for monitoring PETSc's solvers.
Definition petsc.hpp:985
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
Definition petsc.cpp:2538
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition ode.hpp:151
The classical explicit forth-order Runge-Kutta method, RK4.
Definition ode.hpp:164
Data collection with Sidre routines following the Conduit mesh blueprint specification.
Base class for solvers.
Definition operator.hpp:683
Base abstract class for first order time dependent operators.
Definition operator.hpp:317
bool isExplicit() const
True if type is EXPLICIT.
Definition operator.hpp:363
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition operator.hpp:365
virtual void SetTime(const real_t t_)
Set the current time.
Definition operator.hpp:360
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
Data collection with VisIt I/O routines.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void mfem_error(const char *msg)
Definition error.cpp:154
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
Definition petsc.cpp:201
void MFEMFinalizePetsc()
Definition petsc.cpp:241
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]
RefCoord t[3]
RefCoord s[3]
void velocity_function(const Vector &x, Vector &v)
Definition ex9p.cpp:643
real_t inflow_function(const Vector &x)
Definition ex9p.cpp:753
int problem
Definition ex9p.cpp:46
real_t u0_function(const Vector &x)
Definition ex9p.cpp:700
Vector bb_min
Definition ex9p.cpp:58
Vector bb_max
Definition ex9p.cpp:58
HYPRE_Int PetscInt
Definition petsc.hpp:50
real_t PetscReal
Definition petsc.hpp:52