MFEM  v4.6.0
Finite element discretization library
par_example.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // MFEM AD Example - Parallel Version
13 //
14 // Compile with: make par_example
15 //
16 // Sample runs: mpirun -np 2 par_example -m ../data/beam-quad.mesh -pp 3.8
17 // mpirun -np 2 par_example -m ../data/beam-tri.mesh -pp 7.2
18 // mpirun -np 2 par_example -m ../data/beam-hex.mesh
19 // mpirun -np 2 par_example -m ../data/beam-tet.mesh
20 // mpirun -np 2 par_example -m ../data/beam-wedge.mesh
21 //
22 // Description: This examples solves a quasi-static nonlinear p-Laplacian
23 // problem with zero Dirichlet boundary conditions applied on all
24 // defined boundaries
25 //
26 // The example demonstrates the use of nonlinear operators
27 // combined with automatic differentiation (AD). The integrators
28 // are defined in example.hpp. Selecting integrator = 0 will use
29 // the manually implemented integrator. Selecting integrator = 1
30 // or 2 will utilize one of the AD integrators.
31 //
32 // We recommend viewing examples 1 and 19, before viewing this
33 // example.
34 
35 #include "example.hpp"
36 
37 using namespace mfem;
38 
40 {
44 };
45 
46 /// Non-linear solver for the p-Laplacian problem.
47 class ParNLSolverPLaplacian
48 {
49 public:
50  /// Constructor Input: imesh - FE mesh, finite element space, power for the
51  /// p-Laplacian, external load (source, input), regularization parameter
52  ParNLSolverPLaplacian(MPI_Comm comm, ParMesh& imesh,
53  ParFiniteElementSpace& ifespace,
54  double powerp=2,
55  Coefficient* load=nullptr,
56  double regularizationp=1e-7)
57  {
58  lcomm = comm;
59 
60  // default parameters for the Newton solver
61  newton_rtol = 1e-4;
62  newton_atol = 1e-8;
63  newton_iter = 10;
64 
65  // linear solver
66  linear_rtol = 1e-7;
67  linear_atol = 1e-15;
68  linear_iter = 500;
69 
70  print_level = 0;
71 
72  // set the mesh
73  mesh=&imesh;
74 
75  // set the fespace
76  fespace=&ifespace;
77 
78  // set the parameters
79  plap_epsilon=new ConstantCoefficient(regularizationp);
80  plap_power=new ConstantCoefficient(powerp);
81  if (load==nullptr)
82  {
83  plap_input=new ConstantCoefficient(1.0);
84  input_ownership=true;
85  }
86  else
87  {
88  plap_input=load;
89  input_ownership=false;
90  }
91 
92  nlform=nullptr;
93  nsolver=nullptr;
94  gmres=nullptr;
95  prec=nullptr;
96 
97  // set the default integrator
99  }
100 
101  ~ParNLSolverPLaplacian()
102  {
103  delete nlform;
104  delete nsolver;
105  delete prec;
106  delete gmres;
107  if (input_ownership) { delete plap_input;}
108  delete plap_epsilon;
109  delete plap_power;
110  }
111 
112  /// Set the integrator.
113  /// 0 - hand coded, 1 - AD based (compute only Hessian by AD),
114  /// 2 - AD based (compute residual and Hessian by AD)
115  void SetIntegrator(IntegratorType intr)
116  {
117  integ=intr;
118  }
119 
120  // set relative tolerance for the Newton solver
121  void SetNRRTol(double rtol)
122  {
123  newton_rtol=rtol;
124  }
125 
126  // set absolute tolerance for the Newton solver
127  void SetNRATol(double atol)
128  {
129  newton_atol=atol;
130  }
131 
132  // set max iterations for the NR solver
133  void SetMaxNRIter(int miter)
134  {
135  newton_iter=miter;
136  }
137 
138  void SetLSRTol(double rtol)
139  {
140  linear_rtol=rtol;
141  }
142 
143  void SetLSATol(double atol)
144  {
145  linear_atol=atol;
146  }
147 
148  // set max iterations for the linear solver
149  void SetMaxLSIter(int miter)
150  {
151  linear_iter=miter;
152  }
153 
154  // set the print level
155  void SetPrintLevel(int plev)
156  {
157  print_level=plev;
158  }
159 
160  /// The state vector is used as initial condition for the NR solver. On
161  /// return the statev holds the solution to the problem.
162  void Solve(Vector& statev)
163  {
164  if (nlform==nullptr)
165  {
166  AllocSolvers();
167  }
168  Vector b; // RHS is zero
169  nsolver->Mult(b, statev);
170  }
171 
172  /// Compute the energy
173  double GetEnergy(Vector& statev)
174  {
175  if (nlform==nullptr)
176  {
177  // allocate the solvers
178  AllocSolvers();
179  }
180  return nlform->GetEnergy(statev);
181  }
182 
183 private:
184  void AllocSolvers()
185  {
186  if (nlform!=nullptr) { delete nlform;}
187  if (nsolver!=nullptr) { delete nsolver;}
188  if (gmres!=nullptr) { delete gmres;}
189  if (prec!=nullptr) { delete prec;}
190 
191  // Define the essential boundary attributes
192  Array<int> ess_bdr(mesh->bdr_attributes.Max());
193  ess_bdr = 1;
194 
195  nlform = new ParNonlinearForm(fespace);
197  {
198  nlform->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,
199  *plap_input));
200  }
201  else if (integ==IntegratorType::ADJacobianIntegrator)
202  {
203  // The template integrator is based on automatic differentiation. For
204  // ADJacobianIntegrator the residual (vector function) at an
205  // integration point is implemented as a functor by MyResidualFunctor.
206  // The vector function has a return size of four(4), four state
207  // arguments, and three(3) parameters. MyResidualFunctor is a template
208  // argument to the actual template class performing the differentiation
209  // - in this case, QVectorFuncAutoDiff. The derivatives are used in the
210  // integration loop in the integrator pLaplaceAD.
211  nlform->AddDomainIntegrator(new
213  *plap_epsilon,*plap_input));
214  }
215  else if (integ==IntegratorType::ADHessianIntegrator)
216  {
217  // The main difference from the previous case is that the user has to
218  // implement only a functional evaluation at an integration point. The
219  // implementation is in MyEnergyFunctor, which takes four state
220  // arguments and three parameters. The residual vector is the first
221  // derivative of the energy/functional with respect to the state
222  // variables, and the Hessian is the second derivative. Automatic
223  // differentiation is used for evaluating both of them.
224  nlform->AddDomainIntegrator(new
226  *plap_epsilon,*plap_input));
227  }
228 
229  nlform->SetEssentialBC(ess_bdr);
230 
231  prec = new HypreBoomerAMG();
232  prec->SetPrintLevel(print_level);
233 
234  gmres = new GMRESSolver(lcomm);
235  gmres->SetAbsTol(linear_atol);
236  gmres->SetRelTol(linear_rtol);
237  gmres->SetMaxIter(linear_iter);
238  gmres->SetPrintLevel(print_level);
239  gmres->SetPreconditioner(*prec);
240 
241  nsolver = new NewtonSolver(lcomm);
242 
243  nsolver->iterative_mode = true;
244  nsolver->SetSolver(*gmres);
245  nsolver->SetOperator(*nlform);
246  nsolver->SetPrintLevel(print_level);
247  nsolver->SetRelTol(newton_rtol);
248  nsolver->SetAbsTol(newton_atol);
249  nsolver->SetMaxIter(newton_iter);
250  }
251 
252  double newton_rtol;
253  double newton_atol;
254  int newton_iter;
255 
256  double linear_rtol;
257  double linear_atol;
258  int linear_iter;
259 
260  int print_level;
261 
262  // power of the p-laplacian
263  Coefficient* plap_power;
264  // regularization parameter
265  Coefficient* plap_epsilon;
266  // load(input) parameter
267  Coefficient* plap_input;
268  // flag indicating the ownership of plap_input
269  bool input_ownership;
270 
271  MPI_Comm lcomm;
272 
273  ParMesh *mesh;
274  ParFiniteElementSpace *fespace;
275 
276  ParNonlinearForm *nlform;
277 
278  HypreBoomerAMG *prec;
279  GMRESSolver *gmres;
280  NewtonSolver *nsolver;
281  IntegratorType integ;
282 
283 };
284 
285 
286 int main(int argc, char *argv[])
287 {
288  // 1. Initialize MPI and HYPRE.
289  Mpi::Init(argc, argv);
290  int myrank = Mpi::WorldRank();
291  Hypre::Init();
292  // Define Caliper ConfigManager
293 #ifdef MFEM_USE_CALIPER
294  cali::ConfigManager mgr;
295 #endif
296  // Caliper instrumentation
297  MFEM_PERF_FUNCTION;
298 
299  // 2. Parse command-line options
300  const char *mesh_file = "../../data/beam-tet.mesh";
301  int ser_ref_levels = 3;
302  int par_ref_levels = 1;
303  int order = 1;
304  bool visualization = true;
305  double newton_rel_tol = 1e-4;
306  double newton_abs_tol = 1e-6;
307  int newton_iter = 10;
308  int print_level = 0;
309 
310  double pp = 2.0; // p-Laplacian power
311 
313  int int_integrator = integrator;
314  // HandCodedIntegrator = 0 - do not use AD (hand coded)
315  // ADJacobianIntegrator = 1 - use AD for Hessian only
316  // ADHessianIntegrator = 2 - use AD for Residual and Hessian
317 
318  const char* cali_config = "runtime-report";
319 
320  OptionsParser args(argc, argv);
321  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
322  args.AddOption(&ser_ref_levels,
323  "-rs",
324  "--refine-serial",
325  "Number of times to refine the mesh uniformly in serial.");
326  args.AddOption(&par_ref_levels,
327  "-rp",
328  "--refine-parallel",
329  "Number of times to refine the mesh uniformly in parallel.");
330  args.AddOption(&order,
331  "-o",
332  "--order",
333  "Order (degree) of the finite elements.");
334  args.AddOption(&visualization,
335  "-vis",
336  "--visualization",
337  "-no-vis",
338  "--no-visualization",
339  "Enable or disable GLVis visualization.");
340  args.AddOption(&newton_rel_tol,
341  "-rel",
342  "--relative-tolerance",
343  "Relative tolerance for the Newton solve.");
344  args.AddOption(&newton_abs_tol,
345  "-abs",
346  "--absolute-tolerance",
347  "Absolute tolerance for the Newton solve.");
348  args.AddOption(&newton_iter,
349  "-it",
350  "--newton-iterations",
351  "Maximum iterations for the Newton solve.");
352  args.AddOption(&pp,
353  "-pp",
354  "--power-parameter",
355  "Power parameter (>=2.0) for the p-Laplacian.");
356  args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
357  args.AddOption(&int_integrator,
358  "-int",
359  "--integrator",
360  "Integrator 0: standard; 1: AD for Hessian; 2: AD for residual and Hessian");
361  args.AddOption(&cali_config, "-p", "--caliper",
362  "Caliper configuration string.");
363  args.Parse();
364  if (!args.Good())
365  {
366  if (myrank == 0)
367  {
368  args.PrintUsage(std::cout);
369  }
370  return 1;
371  }
372  if (myrank == 0)
373  {
374  args.PrintOptions(std::cout);
375  }
376  integrator = static_cast<IntegratorType>(int_integrator);
377 
378  StopWatch *timer = new StopWatch();
379 
380  // Caliper configuration
381 #ifdef MFEM_USE_CALIPER
382  mgr.add(cali_config);
383  mgr.start();
384 #endif
385  // 3. Read the (serial) mesh from the given mesh file on all processors. We
386  // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
387  // with the same code.
388  Mesh *mesh = new Mesh(mesh_file, 1, 1);
389  int dim = mesh->Dimension();
390 
391  // 4. Refine the mesh in serial to increase the resolution. In this example
392  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
393  // a command-line parameter.
394  for (int lev = 0; lev < ser_ref_levels; lev++)
395  {
396  mesh->UniformRefinement();
397  }
398 
399  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
400  // this mesh further in parallel to increase the resolution. Once the
401  // parallel mesh is defined, the serial mesh can be deleted.
402  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
403  delete mesh;
404  for (int lev = 0; lev < par_ref_levels; lev++)
405  {
406  pmesh->UniformRefinement();
407  }
408 
409  // 6. Define the load for the p-Laplacian
410  ConstantCoefficient load(1.00);
411 
412  // 7. Define the finite element spaces for the solution
413  H1_FECollection fec(order, dim);
414  ParFiniteElementSpace fespace(pmesh, &fec, 1, Ordering::byVDIM);
415  HYPRE_Int glob_size = fespace.GlobalTrueVSize();
416  if (myrank == 0)
417  {
418  std::cout << "Number of finite element unknowns: " << glob_size
419  << std::endl;
420  }
421 
422  // 8. Define the solution vector x as a parallel finite element grid function
423  // corresponding to fespace. Initialize x with initial guess of zero,
424  // which satisfies the boundary conditions.
425  ParGridFunction x(&fespace);
426  x = 0.0;
427  HypreParVector *sv = x.GetTrueDofs();
428 
429  // 9. Define ParaView DataCollection
430  ParaViewDataCollection *dacol = new ParaViewDataCollection("Example",
431  pmesh);
432  dacol->SetLevelsOfDetail(order);
433  dacol->RegisterField("sol", &x);
434 
435  // 10. Define the NR solver
436  ParNLSolverPLaplacian* nr;
437 
438  // 11. Start with linear diffusion - solvable for any initial guess
439  nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, 2.0, &load);
440  nr->SetIntegrator(integrator);
441  nr->SetMaxNRIter(newton_iter);
442  nr->SetNRATol(newton_abs_tol);
443  nr->SetNRRTol(newton_rel_tol);
444  nr->SetPrintLevel(print_level);
445  timer->Clear();
446  timer->Start();
447  nr->Solve(*sv);
448  timer->Stop();
449  if (myrank==0)
450  {
451  std::cout << "[pp=2] The solution time is: " << timer->RealTime()
452  << std::endl;
453  }
454  // Compute the energy
455  double energy = nr->GetEnergy(*sv);
456  if (myrank==0)
457  {
458  std::cout << "[pp=2] The total energy of the system is E=" << energy
459  << std::endl;
460  }
461  delete nr;
462  x.SetFromTrueDofs(*sv);
463  dacol->SetTime(2.0);
464  dacol->SetCycle(2);
465  dacol->Save();
466 
467  // 12. Continue with powers higher than 2
468  for (int i = 3; i < pp; i++)
469  {
470  nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, (double)i, &load);
471  nr->SetIntegrator(integrator);
472  nr->SetMaxNRIter(newton_iter);
473  nr->SetNRATol(newton_abs_tol);
474  nr->SetNRRTol(newton_rel_tol);
475  nr->SetPrintLevel(print_level);
476  timer->Clear();
477  timer->Start();
478  nr->Solve(*sv);
479  timer->Stop();
480  if (myrank==0)
481  {
482  std::cout << "[pp="<<i<<"] The solution time is: " << timer->RealTime()
483  << std::endl;
484  }
485  // Compute the energy
486  energy = nr->GetEnergy(*sv);
487  if (myrank==0)
488  {
489  std::cout << "[pp="<<i<<"] The total energy of the system is E=" << energy
490  << std::endl;
491  }
492  delete nr;
493  x.SetFromTrueDofs(*sv);
494  dacol->SetTime((double)i);
495  dacol->SetCycle(i);
496  dacol->Save();
497  }
498 
499  // 13. Continue with the final power
500  if (std::abs(pp - 2.0) > std::numeric_limits<double>::epsilon())
501  {
502  nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, pp, &load);
503  nr->SetIntegrator(integrator);
504  nr->SetMaxNRIter(newton_iter);
505  nr->SetNRATol(newton_abs_tol);
506  nr->SetNRRTol(newton_rel_tol);
507  nr->SetPrintLevel(print_level);
508  timer->Clear();
509  timer->Start();
510  nr->Solve(*sv);
511  timer->Stop();
512  if (myrank==0)
513  {
514  std::cout << "[pp="<<pp<<"] The solution time is: " << timer->RealTime()
515  << std::endl;
516  }
517  // Compute the energy
518  energy = nr->GetEnergy(*sv);
519  if (myrank==0)
520  {
521  std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
522  << std::endl;
523  }
524  delete nr;
525  x.SetFromTrueDofs(*sv);
526  dacol->SetTime(pp);
527  if (pp < 2.0)
528  {
529  dacol->SetCycle(static_cast<int>(std::floor(pp)));
530  }
531  else
532  {
533  dacol->SetCycle(static_cast<int>(std::ceil(pp)));
534  }
535  dacol->Save();
536  }
537 
538  // 14. Free the used memory
539  delete dacol;
540  delete sv;
541  delete pmesh;
542  delete timer;
543 
544  // Flush output before MPI_finalize
545 #ifdef MFEM_USE_CALIPER
546  mgr.flush();
547 #endif
548 
549  return 0;
550 }
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:80
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
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
double epsilon
Definition: ex25.cpp:141
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Parallel non-linear operator on the true dofs.
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition: tic_toc.cpp:429
int main(int argc, char *argv[])
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
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:419
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
Timing object.
Definition: tic_toc.hpp:35
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:641
static void Init()
Singleton creation with Mpi::Init();.
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:408
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
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
GMRES method.
Definition: solvers.hpp:525
int dim
Definition: ex24.cpp:53
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
void SetLevelsOfDetail(int levels_of_detail_)
IntegratorType
Definition: par_example.cpp:39
Vector data type.
Definition: vector.hpp:58
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:152
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual void Save() override
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:403