MFEM  v4.6.0
Finite element discretization library
ex6p.cpp
Go to the documentation of this file.
1 // MFEM Example 6 - Parallel Version
2 //
3 // Compile with: make ex6p
4 //
5 // Sample runs: mpirun -np 4 ex6p -m ../data/star-hilbert.mesh -o 2
6 // mpirun -np 4 ex6p -m ../data/square-disc.mesh -rm 1 -o 1
7 // mpirun -np 4 ex6p -m ../data/square-disc.mesh -rm 1 -o 2 -h1
8 // mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 2 -cs
9 // mpirun -np 4 ex6p -m ../data/square-disc-nurbs.mesh -o 2
10 // mpirun -np 4 ex6p -m ../data/fichera.mesh -o 2
11 // mpirun -np 4 ex6p -m ../data/escher.mesh -rm 2 -o 2
12 // mpirun -np 4 ex6p -m ../data/escher.mesh -o 2 -cs
13 // mpirun -np 4 ex6p -m ../data/disc-nurbs.mesh -o 2
14 // mpirun -np 4 ex6p -m ../data/ball-nurbs.mesh
15 // mpirun -np 4 ex6p -m ../data/pipe-nurbs.mesh
16 // mpirun -np 4 ex6p -m ../data/star-surf.mesh -o 2
17 // mpirun -np 4 ex6p -m ../data/square-disc-surf.mesh -rm 2 -o 2
18 // mpirun -np 4 ex6p -m ../data/inline-segment.mesh -o 1 -md 200
19 // mpirun -np 4 ex6p -m ../data/amr-quad.mesh
20 // mpirun -np 4 ex6p --restart
21 //
22 // Device sample runs:
23 // mpirun -np 4 ex6p -pa -d cuda
24 // mpirun -np 4 ex6p -pa -d occa-cuda
25 // mpirun -np 4 ex6p -pa -d raja-omp
26 // mpirun -np 4 ex6p -pa -d ceed-cpu
27 // * mpirun -np 4 ex6p -pa -d ceed-cuda
28 // mpirun -np 4 ex6p -pa -d ceed-cuda:/gpu/cuda/shared
29 //
30 // Description: This is a version of Example 1 with a simple adaptive mesh
31 // refinement loop. The problem being solved is again the Laplace
32 // equation -Delta u = 1 with homogeneous Dirichlet boundary
33 // conditions. The problem is solved on a sequence of meshes which
34 // are locally refined in a conforming (triangles, tetrahedrons)
35 // or non-conforming (quadrilaterals, hexahedra) manner according
36 // to a simple ZZ error estimator.
37 //
38 // The example demonstrates MFEM's capability to work with both
39 // conforming and nonconforming refinements, in 2D and 3D, on
40 // linear, curved and surface meshes. Interpolation of functions
41 // from coarse to fine meshes, restarting from a checkpoint, as
42 // well as persistent GLVis visualization are also illustrated.
43 //
44 // We recommend viewing Example 1 before viewing this example.
45 
46 #include "mfem.hpp"
47 #include <fstream>
48 #include <iostream>
49 
50 using namespace std;
51 using namespace mfem;
52 
53 int main(int argc, char *argv[])
54 {
55  // 1. Initialize MPI and HYPRE.
56  Mpi::Init(argc, argv);
57  int num_procs = Mpi::WorldSize();
58  int myid = Mpi::WorldRank();
59  Hypre::Init();
60 
61  // 2. Parse command-line options.
62  const char *mesh_file = "../data/star.mesh";
63  int order = 1;
64  bool pa = false;
65  const char *device_config = "cpu";
66  bool nc_simplices = true;
67  int reorder_mesh = 0;
68  int max_dofs = 100000;
69  bool smooth_rt = true;
70  bool restart = false;
71  bool visualization = true;
72 
73  OptionsParser args(argc, argv);
74  args.AddOption(&mesh_file, "-m", "--mesh",
75  "Mesh file to use.");
76  args.AddOption(&order, "-o", "--order",
77  "Finite element order (polynomial degree).");
78  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
79  "--no-partial-assembly", "Enable Partial Assembly.");
80  args.AddOption(&device_config, "-d", "--device",
81  "Device configuration string, see Device::Configure().");
82  args.AddOption(&reorder_mesh, "-rm", "--reorder-mesh",
83  "Reorder elements of the coarse mesh to improve "
84  "dynamic partitioning: 0=none, 1=hilbert, 2=gecko.");
85  args.AddOption(&nc_simplices, "-ns", "--nonconforming-simplices",
86  "-cs", "--conforming-simplices",
87  "For simplicial meshes, enable/disable nonconforming"
88  " refinement");
89  args.AddOption(&max_dofs, "-md", "--max-dofs",
90  "Stop after reaching this many degrees of freedom.");
91  args.AddOption(&smooth_rt, "-rt", "--smooth-rt", "-h1", "--smooth-h1",
92  "Represent the smooth flux in RT or vector H1 space.");
93  args.AddOption(&restart, "-res", "--restart", "-no-res", "--no-restart",
94  "Restart computation from the last checkpoint.");
95  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
96  "--no-visualization",
97  "Enable or disable GLVis visualization.");
98  args.Parse();
99  if (!args.Good())
100  {
101  if (myid == 0)
102  {
103  args.PrintUsage(cout);
104  }
105  return 1;
106  }
107  if (myid == 0)
108  {
109  args.PrintOptions(cout);
110  }
111 
112  // 3. Enable hardware devices such as GPUs, and programming models such as
113  // CUDA, OCCA, RAJA and OpenMP based on command line options.
114  Device device(device_config);
115  if (myid == 0) { device.Print(); }
116 
117  ParMesh *pmesh;
118  if (!restart)
119  {
120  // 4. Read the (serial) mesh from the given mesh file on all processors.
121  // We can handle triangular, quadrilateral, tetrahedral, hexahedral,
122  // surface and volume meshes with the same code.
123  Mesh mesh(mesh_file, 1, 1);
124 
125  // 5. A NURBS mesh cannot be refined locally so we refine it uniformly
126  // and project it to a standard curvilinear mesh of order 2.
127  if (mesh.NURBSext)
128  {
129  mesh.UniformRefinement();
130  mesh.SetCurvature(2);
131  }
132 
133  // 6. MFEM supports dynamic partitioning (load balancing) of parallel non-
134  // conforming meshes based on space-filling curve (SFC) partitioning.
135  // SFC partitioning is extremely fast and scales to hundreds of
136  // thousands of processors, but requires the coarse mesh to be ordered,
137  // ideally as a sequence of face-neighbors. The mesh may already be
138  // ordered (like star-hilbert.mesh) or we can order it here. Ordering
139  // type 1 is a fast spatial sort of the mesh, type 2 is a high quality
140  // optimization algorithm suitable for ordering general unstructured
141  // meshes.
142  if (reorder_mesh)
143  {
144  Array<int> ordering;
145  switch (reorder_mesh)
146  {
147  case 1: mesh.GetHilbertElementOrdering(ordering); break;
148  case 2: mesh.GetGeckoElementOrdering(ordering); break;
149  default: MFEM_ABORT("Unknown mesh reodering type " << reorder_mesh);
150  }
151  mesh.ReorderElements(ordering);
152  }
153 
154  // 7. Make sure the mesh is in the non-conforming mode to enable local
155  // refinement of quadrilaterals/hexahedra, and the above partitioning
156  // algorithm. Simplices can be refined either in conforming or in non-
157  // conforming mode. The conforming mode however does not support
158  // dynamic partitioning.
159  mesh.EnsureNCMesh(nc_simplices);
160 
161  // 8. Define a parallel mesh by partitioning the serial mesh.
162  // Once the parallel mesh is defined, the serial mesh can be deleted.
163  pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
164  }
165  else
166  {
167  // 9. We can also restart the computation by loading the mesh from a
168  // previously saved check-point.
169  string fname(MakeParFilename("ex6p-checkpoint.", myid));
170  ifstream ifs(fname);
171  MFEM_VERIFY(ifs.good(), "Checkpoint file " << fname << " not found.");
172  pmesh = new ParMesh(MPI_COMM_WORLD, ifs);
173  }
174 
175  int dim = pmesh->Dimension();
176  int sdim = pmesh->SpaceDimension();
177 
178  MFEM_VERIFY(pmesh->bdr_attributes.Size() > 0,
179  "Boundary attributes required in the mesh.");
180  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
181  ess_bdr = 1;
182 
183  // 10. Define a finite element space on the mesh. The polynomial order is
184  // one (linear) by default, but this can be changed on the command line.
185  H1_FECollection fec(order, dim);
186  ParFiniteElementSpace fespace(pmesh, &fec);
187 
188  // 11. As in Example 1p, we set up bilinear and linear forms corresponding to
189  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
190  // problem yet, this will be done in the main loop.
191  ParBilinearForm a(&fespace);
192  if (pa)
193  {
194  a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
195  a.SetDiagonalPolicy(Operator::DIAG_ONE);
196  }
197  ParLinearForm b(&fespace);
198 
199  ConstantCoefficient one(1.0);
200 
202  a.AddDomainIntegrator(integ);
203  b.AddDomainIntegrator(new DomainLFIntegrator(one));
204 
205  // 12. The solution vector x and the associated finite element grid function
206  // will be maintained over the AMR iterations. We initialize it to zero.
207  ParGridFunction x(&fespace);
208  x = 0;
209 
210  // 13. Connect to GLVis.
211  char vishost[] = "localhost";
212  int visport = 19916;
213 
214  socketstream sout;
215  if (visualization)
216  {
217  sout.open(vishost, visport);
218  if (!sout)
219  {
220  if (myid == 0)
221  {
222  cout << "Unable to connect to GLVis server at "
223  << vishost << ':' << visport << endl;
224  cout << "GLVis visualization disabled.\n";
225  }
226  visualization = false;
227  }
228 
229  sout.precision(8);
230  }
231 
232  // 14. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
233  // with L2 projection in the smoothing step to better handle hanging
234  // nodes and parallel partitioning. We need to supply a space for the
235  // discontinuous flux (L2) and a space for the smoothed flux.
236  L2_FECollection flux_fec(order, dim);
237  ParFiniteElementSpace flux_fes(pmesh, &flux_fec, sdim);
238  FiniteElementCollection *smooth_flux_fec = NULL;
239  ParFiniteElementSpace *smooth_flux_fes = NULL;
240  if (smooth_rt && dim > 1)
241  {
242  // Use an H(div) space for the smoothed flux (this is the default).
243  smooth_flux_fec = new RT_FECollection(order-1, dim);
244  smooth_flux_fes = new ParFiniteElementSpace(pmesh, smooth_flux_fec, 1);
245  }
246  else
247  {
248  // Another possible option for the smoothed flux space: H1^dim space
249  smooth_flux_fec = new H1_FECollection(order, dim);
250  smooth_flux_fes = new ParFiniteElementSpace(pmesh, smooth_flux_fec, dim);
251  }
252  L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fes, *smooth_flux_fes);
253 
254  // 15. A refiner selects and refines elements based on a refinement strategy.
255  // The strategy here is to refine elements with errors larger than a
256  // fraction of the maximum element error. Other strategies are possible.
257  // The refiner will call the given error estimator.
258  ThresholdRefiner refiner(estimator);
259  refiner.SetTotalErrorFraction(0.7);
260 
261  // 16. The main AMR loop. In each iteration we solve the problem on the
262  // current mesh, visualize the solution, and refine the mesh.
263  for (int it = 0; ; it++)
264  {
265  HYPRE_BigInt global_dofs = fespace.GlobalTrueVSize();
266  if (myid == 0)
267  {
268  cout << "\nAMR iteration " << it << endl;
269  cout << "Number of unknowns: " << global_dofs << endl;
270  }
271 
272  // 17. Assemble the right-hand side and determine the list of true
273  // (i.e. parallel conforming) essential boundary dofs.
274  Array<int> ess_tdof_list;
275  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
276  b.Assemble();
277 
278  // 18. Assemble the stiffness matrix. Note that MFEM doesn't care at this
279  // point that the mesh is nonconforming and parallel. The FE space is
280  // considered 'cut' along hanging edges/faces, and also across
281  // processor boundaries.
282  a.Assemble();
283 
284  // 19. Create the parallel linear system: eliminate boundary conditions.
285  // The system will be solved for true (unconstrained/unique) DOFs only.
286  OperatorPtr A;
287  Vector B, X;
288 
289  const int copy_interior = 1;
290  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
291 
292  // 20. Solve the linear system A X = B.
293  // * With full assembly, use the BoomerAMG preconditioner from hypre.
294  // * With partial assembly, use a diagonal preconditioner.
295  Solver *M = NULL;
296  if (pa)
297  {
298  M = new OperatorJacobiSmoother(a, ess_tdof_list);
299  }
300  else
301  {
302  HypreBoomerAMG *amg = new HypreBoomerAMG;
303  amg->SetPrintLevel(0);
304  M = amg;
305  }
306  CGSolver cg(MPI_COMM_WORLD);
307  cg.SetRelTol(1e-6);
308  cg.SetMaxIter(2000);
309  cg.SetPrintLevel(3); // print the first and the last iterations only
310  cg.SetPreconditioner(*M);
311  cg.SetOperator(*A);
312  cg.Mult(B, X);
313  delete M;
314 
315  // 21. Switch back to the host and extract the parallel grid function
316  // corresponding to the finite element approximation X. This is the
317  // local solution on each processor.
318  a.RecoverFEMSolution(X, b, x);
319 
320  // 22. Send the solution by socket to a GLVis server.
321  if (visualization)
322  {
323  sout << "parallel " << num_procs << " " << myid << "\n";
324  sout << "solution\n" << *pmesh << x << flush;
325  }
326 
327  if (global_dofs >= max_dofs)
328  {
329  if (myid == 0)
330  {
331  cout << "Reached the maximum number of dofs. Stop." << endl;
332  }
333  break;
334  }
335 
336  // 23. Call the refiner to modify the mesh. The refiner calls the error
337  // estimator to obtain element errors, then it selects elements to be
338  // refined and finally it modifies the mesh. The Stop() method can be
339  // used to determine if a stopping criterion was met.
340  refiner.Apply(*pmesh);
341  if (refiner.Stop())
342  {
343  if (myid == 0)
344  {
345  cout << "Stopping criterion satisfied. Stop." << endl;
346  }
347  break;
348  }
349 
350  // 24. Update the finite element space (recalculate the number of DOFs,
351  // etc.) and create a grid function update matrix. Apply the matrix
352  // to any GridFunctions over the space. In this case, the update
353  // matrix is an interpolation matrix so the updated GridFunction will
354  // still represent the same function as before refinement.
355  fespace.Update();
356  x.Update();
357 
358  // 25. Load balance the mesh, and update the space and solution. Currently
359  // available only for nonconforming meshes.
360  if (pmesh->Nonconforming())
361  {
362  pmesh->Rebalance();
363 
364  // Update the space and the GridFunction. This time the update matrix
365  // redistributes the GridFunction among the processors.
366  fespace.Update();
367  x.Update();
368  }
369 
370  // 26. Inform also the bilinear and linear forms that the space has
371  // changed.
372  a.Update();
373  b.Update();
374 
375  // 27. Save the current state of the mesh every 5 iterations. The
376  // computation can be restarted from this point. Note that unlike in
377  // visualization, we need to use the 'ParPrint' method to save all
378  // internal parallel data structures.
379  if ((it + 1) % 5 == 0)
380  {
381  ofstream ofs(MakeParFilename("ex6p-checkpoint.", myid));
382  ofs.precision(8);
383  pmesh->ParPrint(ofs);
384 
385  if (myid == 0)
386  {
387  cout << "\nCheckpoint saved." << endl;
388  }
389  }
390  }
391 
392  delete smooth_flux_fes;
393  delete smooth_flux_fec;
394  delete pmesh;
395 
396  return 0;
397 }
double GetGeckoElementOrdering(Array< int > &ordering, int iterations=4, int window=4, int period=2, int seed=0, bool verbose=false, double time_limit=0)
Definition: mesh.cpp:2050
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
int visport
Conjugate gradient method.
Definition: solvers.hpp:493
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
std::string MakeParFilename(const std::string &prefix, const int myid, const std::string suffix, const int width)
Construct a string of the form "<prefix><myid><suffix>" where the integer myid is padded with leading...
Definition: globals.cpp:43
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3323
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
bool Nonconforming() const
Definition: mesh.hpp:1969
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
The L2ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure where t...
Definition: estimators.hpp:328
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
bool Apply(Mesh &mesh)
Perform the mesh operation.
void Rebalance()
Definition: pmesh.cpp:4044
Class for parallel linear form.
Definition: plinearform.hpp:26
void ParPrint(std::ostream &out) const
Save the mesh in a parallel mesh format.
Definition: pmesh.cpp:6295
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9888
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
Mesh refinement operator using an error threshold.
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
void GetHilbertElementOrdering(Array< int > &ordering)
Definition: mesh.cpp:2217
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
HYPRE_Int HYPRE_BigInt
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
void ReorderElements(const Array< int > &ordering, bool reorder_vertices=true)
Definition: mesh.cpp:2269
int dim
Definition: ex24.cpp:53
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Vector data type.
Definition: vector.hpp:58
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
Class for parallel meshes.
Definition: pmesh.hpp:32
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
int main(int argc, char *argv[])
Definition: ex6p.cpp:53