MFEM  v3.3 Finite element discretization library
ex15p.cpp
Go to the documentation of this file.
1 // MFEM Example 15 - Parallel Version
2 //
3 // Compile with: make ex15p
4 //
5 // Sample runs: mpirun -np 4 ex15p
6 // mpirun -np 4 ex15p -o 1 -y 0.2
7 // mpirun -np 4 ex15p -o 4 -y 0.1
8 // mpirun -np 4 ex15p -n 5
9 // mpirun -np 4 ex15p -p 1 -n 3
10 //
11 // Other meshes:
12 //
13 // mpirun -np 4 ex15p -m ../data/square-disc-nurbs.mesh
14 // mpirun -np 4 ex15p -m ../data/disc-nurbs.mesh
15 // mpirun -np 4 ex15p -m ../data/fichera.mesh -tf 0.5
16 // mpirun -np 4 ex15p -m ../data/ball-nurbs.mesh -tf 0.5
17 // mpirun -np 4 ex15p -m ../data/mobius-strip.mesh
18 // mpirun -np 4 ex15p -m ../data/amr-quad.mesh
19 //
20 // Conforming meshes (no load balancing and derefinement):
21 //
22 // mpirun -np 4 ex15p -m ../data/square-disc.mesh
23 // mpirun -np 4 ex15p -m ../data/escher.mesh -r 2 -tf 0.3
24 //
25 // Description: Building on Example 6, this example demonstrates dynamic AMR.
26 // The mesh is adapted to a time-dependent solution by refinement
27 // as well as by derefinement. For simplicity, the solution is
28 // prescribed and no time integration is done. However, the error
29 // estimation and refinement/derefinement decisions are realistic.
30 //
31 // At each outer iteration the right hand side function is changed
32 // to mimic a time dependent problem. Within each inner iteration
33 // the problem is solved on a sequence of meshes which are locally
34 // refined according to a simple ZZ error estimator. At the end
35 // of the inner iteration the error estimates are also used to
36 // identify any elements which may be over-refined and a single
37 // derefinement step is performed. After each refinement or
38 // derefinement step a rebalance operation is performed to keep
39 // the mesh evenly distributed among the available processors.
40 //
41 // The example demonstrates MFEM's capability to refine, derefine
42 // and load balance nonconforming meshes, in 2D and 3D, and on
43 // linear, curved and surface meshes. Interpolation of functions
44 // between coarse and fine meshes, persistent GLVis visualization,
45 // and saving of time-dependent fields for external visualization
46 // with VisIt (visit.llnl.gov) are also illustrated.
47 //
48 // We recommend viewing Examples 1, 6 and 9 before viewing this
49 // example.
50
51 #include "mfem.hpp"
52 #include <fstream>
53 #include <iostream>
54
55 using namespace std;
56 using namespace mfem;
57
58 // Choices for the problem setup. Affect bdr_func and rhs_func.
59 int problem;
61
62 // Prescribed time-dependent boundary and right-hand side functions.
63 double bdr_func(const Vector &pt, double t);
64 double rhs_func(const Vector &pt, double t);
65
66 // Update the finite element space, interpolate the solution and perform
70  ParLinearForm &b);
71
72
73 int main(int argc, char *argv[])
74 {
75  // 1. Initialize MPI.
76  int num_procs, myid;
77  MPI_Init(&argc, &argv);
78  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
79  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
80
81  // 2. Parse command-line options.
82  problem = 0;
83  nfeatures = 1;
84  const char *mesh_file = "../data/star-hilbert.mesh";
85  int order = 2;
86  double t_final = 1.0;
87  double max_elem_error = 1.0e-4;
88  double hysteresis = 0.25; // derefinement safety coefficient
89  int ref_levels = 0;
90  int nc_limit = 3; // maximum level of hanging nodes
91  bool visualization = true;
92  bool visit = false;
93
94  OptionsParser args(argc, argv);
96  "Mesh file to use.");
98  "Problem setup to use: 0 = spherical front, 1 = ball.");
100  "Number of solution features (fronts/balls).");
102  "Finite element order (polynomial degree).");
104  "Maximum element error");
106  "Derefinement safety coefficient.");
108  "Number of initial uniform refinement levels.");
110  "Maximum level of hanging nodes.");
112  "Final time; start time is 0.");
114  "--no-visualization",
115  "Enable or disable GLVis visualization.");
117  "--no-visit-datafiles",
118  "Save data files for VisIt (visit.llnl.gov) visualization.");
119  args.Parse();
120  if (!args.Good())
121  {
122  if (myid == 0)
123  {
124  args.PrintUsage(cout);
125  }
126  MPI_Finalize();
127  return 1;
128  }
129  if (myid == 0)
130  {
131  args.PrintOptions(cout);
132  }
133
134  // 3. Read the (serial) mesh from the given mesh file on all processors. We
135  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
136  // and volume meshes with the same code.
137  Mesh *mesh = new Mesh(mesh_file, 1, 1);
138  int dim = mesh->Dimension();
139  int sdim = mesh->SpaceDimension();
140
141  // 4. Project a NURBS mesh to a piecewise-quadratic curved mesh. Make sure
142  // that the mesh is non-conforming if it has quads or hexes and refine it.
143  if (mesh->NURBSext)
144  {
145  mesh->UniformRefinement();
146  if (ref_levels > 0) { ref_levels--; }
147  mesh->SetCurvature(2);
148  }
149  mesh->EnsureNCMesh();
150  for (int l = 0; l < ref_levels; l++)
151  {
152  mesh->UniformRefinement();
153  }
154
155  // 5. Define a parallel mesh by partitioning the serial mesh. Once the
156  // parallel mesh is defined, the serial mesh can be deleted.
157  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
158  delete mesh;
159
160  MFEM_VERIFY(pmesh.bdr_attributes.Size() > 0,
161  "Boundary attributes required in the mesh.");
162  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
163  ess_bdr = 1;
164
165  // 6. Define a finite element space on the mesh. The polynomial order is one
166  // (linear) by default, but this can be changed on the command line.
167  H1_FECollection fec(order, dim);
168  ParFiniteElementSpace fespace(&pmesh, &fec);
169
170  // 7. As in Example 1p, we set up bilinear and linear forms corresponding to
171  // the Laplace problem -\Delta u = 1. We don't assemble the discrete
172  // problem yet, this will be done in the inner loop.
173  ParBilinearForm a(&fespace);
174  ParLinearForm b(&fespace);
175
176  ConstantCoefficient one(1.0);
179
183
184  // 8. The solution vector x and the associated finite element grid function
185  // will be maintained over the AMR iterations.
186  ParGridFunction x(&fespace);
187
188  // 9. Connect to GLVis. Prepare for VisIt output.
189  char vishost[] = "localhost";
190  int visport = 19916;
191
192  socketstream sout;
193  if (visualization)
194  {
195  sout.open(vishost, visport);
196  if (!sout)
197  {
198  if (myid == 0)
199  {
200  cout << "Unable to connect to GLVis server at "
201  << vishost << ':' << visport << endl;
202  cout << "GLVis visualization disabled.\n";
203  }
204  visualization = false;
205  }
206  sout.precision(8);
207  }
208
209  VisItDataCollection visit_dc("Example15-Parallel", &pmesh);
210  visit_dc.RegisterField("solution", &x);
211  int vis_cycle = 0;
212
213  // 10. As in Example 6p, we set up a Zienkiewicz-Zhu estimator that will be
214  // used to obtain element error indicators. The integrator needs to
215  // provide the method ComputeElementFlux. We supply an L2 space for the
216  // discontinuous flux and an H(div) space for the smoothed flux.
217  L2_FECollection flux_fec(order, dim);
218  ParFiniteElementSpace flux_fes(&pmesh, &flux_fec, sdim);
219  RT_FECollection smooth_flux_fec(order-1, dim);
220  ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec);
221  L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fes, smooth_flux_fes);
222
223  // 11. As in Example 6p, we also need a refiner. This time the refinement
224  // strategy is based on a fixed threshold that is applied locally to each
225  // element. The global threshold is turned off by setting the total error
226  // fraction to zero. We also enforce a maximum refinement ratio between
228  ThresholdRefiner refiner(estimator);
229  refiner.SetTotalErrorFraction(0.0); // use purely local threshold
230  refiner.SetLocalErrorGoal(max_elem_error);
231  refiner.PreferConformingRefinement();
232  refiner.SetNCLimit(nc_limit);
233
234  // 12. A derefiner selects groups of elements that can be coarsened to form
235  // a larger element. A conservative enough threshold needs to be set to
236  // prevent derefining elements that would immediately be refined again.
237  ThresholdDerefiner derefiner(estimator);
238  derefiner.SetThreshold(hysteresis * max_elem_error);
239  derefiner.SetNCLimit(nc_limit);
240
241  // 13. The outer time loop. In each iteration we update the right hand side,
242  // solve the problem on the current mesh, visualize the solution and
243  // refine the mesh as many times as necessary. Then we derefine any
244  // elements which have very small errors.
245  for (double time = 0.0; time < t_final + 1e-10; time += 0.01)
246  {
247  if (myid == 0)
248  {
249  cout << "\nTime " << time << "\n\nRefinement:" << endl;
250  }
251
252  // Set the current time in the coefficients
253  bdr.SetTime(time);
254  rhs.SetTime(time);
255
256  // Make sure errors will be recomputed in the following.
257  refiner.Reset();
258  derefiner.Reset();
259
260  // 14. The inner refinement loop. At the end we want to have the current
261  // time step resolved to the prescribed tolerance in each element.
262  for (int ref_it = 1; ; ref_it++)
263  {
264  HYPRE_Int global_dofs = fespace.GlobalTrueVSize();
265  if (myid == 0)
266  {
267  cout << "Iteration: " << ref_it << ", number of unknowns: "
268  << global_dofs << flush;
269  }
270
271  // 15. Recompute the field on the current mesh: assemble the stiffness
272  // matrix and the right-hand side.
273  a.Assemble();
274  b.Assemble();
275
276  // 16. Project the exact solution to the essential DOFs.
277  x.ProjectBdrCoefficient(bdr, ess_bdr);
278
279  // 17. Create and solve the parallel linear system.
280  Array<int> ess_tdof_list;
281  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
282
283  HypreParMatrix A;
284  Vector B, X;
285  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
286
287  HypreBoomerAMG amg(A);
288  amg.SetPrintLevel(0);
289  HyprePCG pcg(A);
290  pcg.SetTol(1e-12);
291  pcg.SetMaxIter(200);
292  pcg.SetPrintLevel(0);
293  pcg.SetPreconditioner(amg);
294  pcg.Mult(B, X);
295
296  // 18. Extract the local solution on each processor.
297  a.RecoverFEMSolution(X, b, x);
298
299  // 19. Send the solution by socket to a GLVis server and optionally
300  // save it in VisIt format.
301  if (visualization)
302  {
303  sout << "parallel " << num_procs << " " << myid << "\n";
304  sout << "solution\n" << pmesh << x << flush;
305  }
306  if (visit)
307  {
308  visit_dc.SetCycle(vis_cycle++);
309  visit_dc.SetTime(time);
310  visit_dc.Save();
311  }
312
313  // 20. Apply the refiner on the mesh. The refiner calls the error
314  // estimator to obtain element errors, then it selects elements to
315  // be refined and finally it modifies the mesh. The Stop() method
316  // determines if all elements satisfy the local threshold.
317  refiner.Apply(pmesh);
318  if (myid == 0)
319  {
320  cout << ", total error: " << estimator.GetTotalError() << endl;
321  }
322
323  // 21. Quit the AMR loop if the termination criterion has been met
324  if (refiner.Stop())
325  {
326  a.Update(); // Free the assembled data
327  break;
328  }
329
330  // 22. Update the space, interpolate the solution, rebalance the mesh.
331  UpdateAndRebalance(pmesh, fespace, x, a, b);
332  }
333
334  // 23. Use error estimates from the last inner iteration to check for
335  // possible derefinements. The derefiner works similarly as the
336  // refiner. The errors are not recomputed because the mesh did not
337  // change (and also the estimator was not Reset() at this time).
338  if (derefiner.Apply(pmesh))
339  {
340  if (myid == 0)
341  {
342  cout << "\nDerefined elements." << endl;
343  }
344
345  // 24. Update the space and the solution, rebalance the mesh.
346  UpdateAndRebalance(pmesh, fespace, x, a, b);
347  }
348  }
349
350  // 25. Exit
351  MPI_Finalize();
352  return 0;
353 }
354
355
358  ParLinearForm &b)
359 {
360  // Update the space: recalculate the number of DOFs and construct a matrix
361  // that will adjust any GridFunctions to the new mesh state.
362  fespace.Update();
363
364  // Interpolate the solution on the new mesh by applying the transformation
365  // matrix computed in the finite element space. Multiple GridFunctions could
366  // be updated here.
367  x.Update();
368
369  if (pmesh.Nonconforming())
370  {
371  // Load balance the mesh.
372  pmesh.Rebalance();
373
374  // Update the space again, this time a GridFunction redistribution matrix
375  // is created. Apply it to the solution.
376  fespace.Update();
377  x.Update();
378  }
379
380  // Inform the linear and bilinear forms that the space has changed.
381  a.Update();
382  b.Update();
383
384  // Free any transformation matrices to save memory.
386 }
387
388
389 const double alpha = 0.02;
390
391 // Spherical front with a Gaussian cross section and radius t
392 double front(double x, double y, double z, double t, int)
393 {
394  double r = sqrt(x*x + y*y + z*z);
395  return exp(-0.5*pow((r - t)/alpha, 2));
396 }
397
398 double front_laplace(double x, double y, double z, double t, int dim)
399 {
400  double x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
401  double r = sqrt(x2 + y2 + z2);
402  double a2 = alpha*alpha, a4 = a2*a2;
403  return -exp(-0.5*pow((r - t)/alpha, 2)) / a4 *
404  (-2*t*(x2 + y2 + z2 - (dim-1)*a2/2)/r + x2 + y2 + z2 + t2 - dim*a2);
405 }
406
407 // Smooth spherical step function with radius t
408 double ball(double x, double y, double z, double t, int)
409 {
410  double r = sqrt(x*x + y*y + z*z);
411  return -atan(2*(r - t)/alpha);
412 }
413
414 double ball_laplace(double x, double y, double z, double t, int dim)
415 {
416  double x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
417  double r = sqrt(x2 + y2 + z2);
418  double a2 = alpha*alpha;
419  double den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
420  return (dim == 2) ? 2*alpha*(a2 + t2 - 4*x2 - 4*y2)/r/den
421  /* */ : 4*alpha*(a2 + t2 - 4*r*t)/r/den;
422 }
423
424 // Composes several features into one function
425 template<typename F0, typename F1>
426 double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
427 {
428  int dim = pt.Size();
429  double x = pt(0), y = pt(1), z = 0.0;
430  if (dim == 3) { z = pt(2); }
431
432  if (problem == 0)
433  {
434  if (nfeatures <= 1)
435  {
436  return f0(x, y, z, t, dim);
437  }
438  else
439  {
440  double sum = 0.0;
441  for (int i = 0; i < nfeatures; i++)
442  {
443  double x0 = 0.5*cos(2*M_PI * i / nfeatures);
444  double y0 = 0.5*sin(2*M_PI * i / nfeatures);
445  sum += f0(x - x0, y - y0, z, t, dim);
446  }
447  return sum;
448  }
449  }
450  else
451  {
452  double sum = 0.0;
453  for (int i = 0; i < nfeatures; i++)
454  {
455  double x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
456  double y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
457  sum += f1(x - x0, y - y0, z, 0.25, dim);
458  }
459  return sum;
460  }
461 }
462
463 // Exact solution, used for the Dirichlet BC.
464 double bdr_func(const Vector &pt, double t)
465 {
466  return composite_func(pt, t, front, ball);
467 }
468
469 // Laplace of the exact solution, used for the right hand side.
470 double rhs_func(const Vector &pt, double t)
471 {
472  return composite_func(pt, t, front_laplace, ball_laplace);
473 }
void SetTol(double tol)
Definition: hypre.cpp:1993
double front_laplace(double x, double y, double z, double t, int dim)
Definition: ex15.cpp:346
int Size() const
Logical size of the array.
Definition: array.hpp:109
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
double bdr_func(const Vector &pt, double t)
Definition: ex15.cpp:412
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Subclass constant coefficient.
Definition: coefficient.hpp:57
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:2198
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:42
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
int Size() const
Returns the size of the vector.
Definition: vector.hpp:106
virtual void Save()
Save the collection and a VisIt root file.
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int main(int argc, char *argv[])
double rhs_func(const Vector &pt, double t)
Definition: ex15.cpp:418
The L2ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure where t...
Definition: estimators.hpp:183
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2008
void Update(ParFiniteElementSpace *pf=NULL)
Definition: plinearform.cpp:21
virtual void Update(FiniteElementSpace *nfes=NULL)
The BoomerAMG solver in hypre.
Definition: hypre.hpp:770
bool Apply(Mesh &mesh)
Perform the mesh operation.
void Rebalance()
Load balance the mesh. NC meshes only.
Definition: pmesh.cpp:2713
bool Nonconforming() const
Definition: mesh.hpp:968
Class for parallel linear form.
Definition: plinearform.hpp:26
int dim
Definition: ex3.cpp:47
virtual void Reset()
Reset the associated estimator.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
virtual void Reset()
Reset the associated estimator.
void SetPrintLevel(int print_level)
Definition: hypre.hpp:811
Mesh refinement operator using an error threshold.
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3278
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:108
int nfeatures
Definition: ex15.cpp:58
void Assemble(int skip_zeros=1)
Assemble the local matrix.
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
void SetTime(double t)
Definition: coefficient.hpp:39
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1998
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:239
void UpdateAndRebalance(ParMesh &pmesh, ParFiniteElementSpace &fespace, ParGridFunction &x, ParBilinearForm &a, ParLinearForm &b)
Definition: ex15p.cpp:356
int SpaceDimension() const
Definition: mesh.hpp:612
Definition: linearform.cpp:19
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:22
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:142
double ball_laplace(double x, double y, double z, double t, int dim)
Definition: ex15.cpp:362
PCG solver in hypre.
Definition: hypre.hpp:630
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: pfespace.cpp:521
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:74
int problem
Definition: ex15.cpp:57
void PreferConformingRefinement()
Use conforming refinement, if possible (triangles, tetrahedra) – this is the default.
double GetTotalError() const
Return the total error from the last error estimate.
Definition: estimators.hpp:267
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:144
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2013
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
Class for parallel bilinear form.
void EnsureNCMesh(bool triangles_nonconforming=false)
Definition: mesh.cpp:6350
double ball(double x, double y, double z, double t, int)
Definition: ex15.cpp:356
int open(const char hostname[], int port)
void SetThreshold(double thresh)
Set the de-refinement threshold. The default value is zero.
void SetLocalErrorGoal(double err_goal)
Set the local stopping criterion: stop when local_err_i &lt;= local_err_goal. The default value is zero...
const double alpha
Definition: ex15.cpp:337
class for C-function coefficient
Vector data type.
Definition: vector.hpp:36
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
Class for parallel grid function.
Definition: pgridfunc.hpp:31
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:174
Free ParGridFunction transformation matrix (if any), to save memory.
Definition: pfespace.hpp:268
double front(double x, double y, double z, double t, int)
Definition: ex15.cpp:340
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:2034
Class for parallel meshes.
Definition: pmesh.hpp:28
double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
Definition: ex15.cpp:374
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Definition: gridfunc.hpp:197
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:195
De-refinement operator using an error threshold.
bool Good() const
Definition: optparser.hpp:120