MFEM  v4.5.2
Finite element discretization library
seq_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 - Serial Version
13 //
14 // Compile with: make seq_example
15 //
16 // Sample runs: seq_example -m ../data/beam-quad.mesh -pp 3.5
17 // seq_example -m ../data/beam-tri.mesh -pp 4.6
18 // seq_example -m ../data/beam-hex.mesh
19 // seq_example -m ../data/beam-tet.mesh
20 // seq_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 NLSolverPLaplacian
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  NLSolverPLaplacian(Mesh& imesh, FiniteElementSpace& ifespace,
53  double powerp=2,
54  Coefficient* load=nullptr,
55  double regularizationp=1e-7)
56  {
57  // default parameters for the Newton solver
58  newton_rtol = 1e-4;
59  newton_atol = 1e-8;
60  newton_iter = 10;
61 
62  // linear solver
63  linear_rtol = 1e-7;
64  linear_atol = 1e-15;
65  linear_iter = 500;
66 
67  print_level = 0;
68 
69  // set the mesh
70  mesh=&imesh;
71 
72  // set the fespace
73  fespace=&ifespace;
74 
75  // set the parameters
76  plap_epsilon=new ConstantCoefficient(regularizationp);
77  plap_power=new ConstantCoefficient(powerp);
78  if (load==nullptr)
79  {
80  plap_input=new ConstantCoefficient(1.0);
81  input_ownership=true;
82  }
83  else
84  {
85  plap_input=load;
86  input_ownership=false;
87  }
88 
89  // set the nonlinear form
90  nlform=nullptr;
91  lsolver=nullptr;
92  prec=nullptr;
93  nsolver=nullptr;
94 
95  // set the default integrator
96  integ=IntegratorType::HandCodedIntegrator; // hand coded
97  }
98 
99  ~NLSolverPLaplacian()
100  {
101  if (nlform!=nullptr) { delete nlform;}
102  if (nsolver!=nullptr) { delete nsolver;}
103  if (prec!=nullptr) { delete prec;}
104  if (lsolver!=nullptr) { delete lsolver;}
105  if (input_ownership) { delete plap_input;}
106  delete plap_epsilon;
107  delete plap_power;
108  }
109 
110  /// Set the integrator.
111  /// 0 - hand coded, 1 - AD based (compute only Hessian by AD),
112  /// 2 - AD based (compute residual and Hessian by AD)
113  void SetIntegrator(IntegratorType intr)
114  {
115  integ=intr;
116  }
117 
118 
119  // set relative tolerance for the Newton solver
120  void SetNRRTol(double rtol)
121  {
122  newton_rtol=rtol;
123  }
124 
125  // set absolute tolerance for the Newton solver
126  void SetNRATol(double atol)
127  {
128  newton_atol=atol;
129  }
130 
131  // set max iterations for the NR solver
132  void SetMaxNRIter(int miter)
133  {
134  newton_iter=miter;
135  }
136 
137  void SetLSRTol(double rtol)
138  {
139  linear_rtol=rtol;
140  }
141 
142  void SetLSATol(double atol)
143  {
144  linear_atol=atol;
145  }
146 
147  // set max iterations for the linear solver
148  void SetMaxLSIter(int miter)
149  {
150  linear_iter=miter;
151  }
152 
153  // set the print level
154  void SetPrintLevel(int plev)
155  {
156  print_level=plev;
157  }
158 
159  /// The state vector is used as initial condition for the NR solver. On
160  /// return the statev holds the solution to the problem.
161  void Solve(Vector& statev)
162  {
163  if (nlform==nullptr)
164  {
165  AllocSolvers();
166  }
167  Vector b; // RHS is zero
168  nsolver->Mult(b, statev);
169  }
170 
171  /// Compute the energy
172  double GetEnergy(Vector& statev)
173  {
174  if (nlform==nullptr)
175  {
176  // allocate the solvers
177  AllocSolvers();
178  }
179  return nlform->GetEnergy(statev);
180  }
181 
182 
183 private:
184 
185  void AllocSolvers()
186  {
187  if (nlform!=nullptr) { delete nlform;}
188  if (nsolver!=nullptr) {delete nsolver;}
189  if (prec!=nullptr) {delete prec;}
190  if (lsolver!=nullptr) { delete lsolver;}
191 
192  // Define the essential boundary attributes
193  Array<int> ess_bdr(mesh->bdr_attributes.Max());
194  ess_bdr = 1;
195 
196  nlform = new NonlinearForm(fespace);
197 
199  {
200  // standard hand coded integrator
201  nlform->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,
202  *plap_input));
203  }
204  else if (integ==IntegratorType::ADJacobianIntegrator)
205  {
206  // The template integrator is based on automatic differentiation. For
207  // ADJacobianIntegrator the residual (vector function) at an
208  // integration point is implemented as a functor by MyVFunctor. The
209  // vector function has a return size of four(4), four state arguments,
210  // and three(3) parameters. MyVFunctor is a template argument to the
211  // actual template class performing the differentiation - in this case,
212  // QVectorFuncAutoDiff. The derivatives are used in the integration
213  // loop in the integrator pLaplaceAD.
214  nlform->AddDomainIntegrator(new
216  *plap_epsilon,*plap_input));
217  }
218  else // IntegratorType::ADHessianIntegrator
219  {
220  // The main difference from the previous case is that the user has to
221  // implement only a functional evaluation at an integration point. The
222  // implementation is in MyQFunctor, which takes four state arguments
223  // and three parameters. The residual vector is the first derivative of
224  // the energy/functional with respect to the state variables, and the
225  // Hessian is the second derivative. Automatic differentiation is used
226  // for evaluating both of them.
227  nlform->AddDomainIntegrator(new
229  *plap_epsilon,*plap_input));
230  }
231 
232  nlform->SetEssentialBC(ess_bdr);
233 
234 #ifdef MFEM_USE_SUITESPARSE
235  prec = new UMFPackSolver();
236 #else
237  prec = new GSSmoother();
238 #endif
239 
240  // allocate the linear solver
241  lsolver=new CGSolver();
242  lsolver->SetRelTol(linear_rtol);
243  lsolver->SetAbsTol(linear_atol);
244  lsolver->SetMaxIter(linear_iter);
245  lsolver->SetPrintLevel(print_level);
246  lsolver->SetPreconditioner(*prec);
247 
248  // allocate the NR solver
249  nsolver = new NewtonSolver();
250  nsolver->iterative_mode = true;
251  nsolver->SetSolver(*lsolver);
252  nsolver->SetOperator(*nlform);
253  nsolver->SetPrintLevel(print_level);
254  nsolver->SetRelTol(newton_rtol);
255  nsolver->SetAbsTol(newton_atol);
256  nsolver->SetMaxIter(newton_iter);
257  }
258 
259  double newton_rtol;
260  double newton_atol;
261  int newton_iter;
262 
263  double linear_rtol;
264  double linear_atol;
265  int linear_iter;
266 
267  int print_level;
268 
269  // reference to the mesh
270  Mesh* mesh;
271  // reference to the fespace
272  FiniteElementSpace *fespace;
273 
274  // nonlinear form for the p-laplacian
275  NonlinearForm *nlform;
276  CGSolver *lsolver; // linear solver
277  Solver *prec; // preconditioner for the linear solver
278  NewtonSolver *nsolver; // NR solver
279  IntegratorType integ;
280 
281  // power of the p-laplacian
282  Coefficient* plap_power;
283  // regularization parameter
284  Coefficient* plap_epsilon;
285  // load(input) parameter
286  Coefficient* plap_input;
287  // flag indicating the ownership of plap_input
288  bool input_ownership;
289 };
290 
291 int main(int argc, char *argv[])
292 {
293  // 1. Parse command-line options
294  const char *mesh_file = "../../data/beam-tet.mesh";
295  int ser_ref_levels = 3;
296  int order = 1;
297  bool visualization = true;
298  double newton_rel_tol = 1e-4;
299  double newton_abs_tol = 1e-6;
300  int newton_iter = 10;
301  int print_level = 0;
302 
303  double pp = 2.0; // p-Laplacian power
304 
306  int int_integrator = integrator;
307  // HandCodedIntegrator = 0 - do not use AD (hand coded)
308  // ADJacobianIntegrator = 1 - use AD for Hessian only
309  // ADHessianIntegrator = 2 - use AD for Residual and Hessian
310  StopWatch *timer = new StopWatch();
311  OptionsParser args(argc, argv);
312  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
313  args.AddOption(&ser_ref_levels,
314  "-rs",
315  "--refine-serial",
316  "Number of times to refine the mesh uniformly in serial.");
317  args.AddOption(&order,
318  "-o",
319  "--order",
320  "Order (degree) of the finite elements.");
321  args.AddOption(&visualization,
322  "-vis",
323  "--visualization",
324  "-no-vis",
325  "--no-visualization",
326  "Enable or disable GLVis visualization.");
327  args.AddOption(&newton_rel_tol,
328  "-rel",
329  "--relative-tolerance",
330  "Relative tolerance for the Newton solve.");
331  args.AddOption(&newton_abs_tol,
332  "-abs",
333  "--absolute-tolerance",
334  "Absolute tolerance for the Newton solve.");
335  args.AddOption(&newton_iter,
336  "-it",
337  "--newton-iterations",
338  "Maximum iterations for the Newton solve.");
339  args.AddOption(&pp,
340  "-pp",
341  "--power-parameter",
342  "Power parameter (>=2.0) for the p-Laplacian.");
343  args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
344  args.AddOption(&int_integrator,
345  "-int",
346  "--integrator",
347  "Integrator 0: standard; 1: AD for Hessian; 2: AD for residual and Hessian");
348  args.Parse();
349  if (!args.Good())
350  {
351  args.PrintUsage(std::cout);
352  return 1;
353  }
354  args.PrintOptions(std::cout);
355  integrator = static_cast<IntegratorType>(int_integrator);
356 
357  // 2. Read the (serial) mesh from the given mesh file.
358  Mesh *mesh = new Mesh(mesh_file, 1, 1);
359  int dim = mesh->Dimension();
360 
361  // 3. Refine the mesh in serial to increase the resolution. In this example
362  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
363  // a command-line parameter.
364  for (int lev = 0; lev < ser_ref_levels; lev++)
365  {
366  mesh->UniformRefinement();
367  }
368 
369  // 4. Define the load parameter for the p-Laplacian
370  ConstantCoefficient load(1.00);
371 
372  // 5. Define the finite element spaces for the solution
373  H1_FECollection fec(order, dim);
374  FiniteElementSpace fespace(mesh, &fec, 1, Ordering::byVDIM);
375  int glob_size = fespace.GetTrueVSize();
376 
377  std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
378 
379  // 6. Define the solution grid function
380  GridFunction x(&fespace);
381  x = 0.0;
382 
383  // 7. Define the solution true vector
384  Vector sv(fespace.GetTrueVSize());
385  sv = 0.0;
386 
387  // 8. Define ParaView DataCollection
388  ParaViewDataCollection *dacol = new ParaViewDataCollection("Example", mesh);
389  dacol->SetLevelsOfDetail(order);
390  dacol->RegisterField("sol", &x);
391 
392  // 9. Define the nonlinear p-Laplacian solver
393  NLSolverPLaplacian* nr;
394 
395  // 10. Start with linear diffusion - solvable for any initial guess
396  nr=new NLSolverPLaplacian(*mesh, fespace, 2.0, &load);
397  nr->SetIntegrator(integrator);
398  nr->SetMaxNRIter(newton_iter);
399  nr->SetNRATol(newton_abs_tol);
400  nr->SetNRRTol(newton_rel_tol);
401  timer->Clear();
402  timer->Start();
403  nr->Solve(sv);
404  timer->Stop();
405  std::cout << "[pp=2] The solution time is: " << timer->RealTime()
406  << std::endl;
407  // Compute the energy
408  double energy = nr->GetEnergy(sv);
409  std::cout << "[pp=2] The total energy of the system is E=" << energy
410  << std::endl;
411  delete nr;
412  x.SetFromTrueDofs(sv);
413  dacol->SetTime(2.0);
414  dacol->SetCycle(2);
415  dacol->Save();
416 
417 
418  // 11. Continue with powers higher than 2
419  for (int i = 3; i < pp; i++)
420  {
421  nr=new NLSolverPLaplacian(*mesh, fespace, (double)i, &load);
422  nr->SetIntegrator(integrator);
423  nr->SetMaxNRIter(newton_iter);
424  nr->SetNRATol(newton_abs_tol);
425  nr->SetNRRTol(newton_rel_tol);
426  timer->Clear();
427  timer->Start();
428  nr->Solve(sv);
429  timer->Stop();
430  std::cout << "[pp=" << i
431  << "] The solution time is: " << timer->RealTime() << std::endl;
432  energy = nr->GetEnergy(sv);
433  std::cout << "[pp="<< i<<"] The total energy of the system is E=" << energy
434  << std::endl;
435  delete nr;
436  x.SetFromTrueDofs(sv);
437  dacol->SetTime(i);
438  dacol->SetCycle(i);
439  dacol->Save();
440  }
441 
442  // 12. Continue with the final power
443  if (std::abs(pp - 2.0) > std::numeric_limits<double>::epsilon())
444  {
445  nr=new NLSolverPLaplacian(*mesh, fespace, pp, &load);
446  nr->SetIntegrator(integrator);
447  nr->SetMaxNRIter(newton_iter);
448  nr->SetNRATol(newton_abs_tol);
449  nr->SetNRRTol(newton_rel_tol);
450  timer->Clear();
451  timer->Start();
452  nr->Solve(sv);
453  timer->Stop();
454  std::cout << "[pp=" << pp
455  << "] The solution time is: " << timer->RealTime() << std::endl;
456  energy = nr->GetEnergy(sv);
457  std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
458  << std::endl;
459  delete nr;
460  x.SetFromTrueDofs(sv);
461  dacol->SetTime(pp);
462  if (pp < 2.0)
463  {
464  dacol->SetCycle(static_cast<int>(std::floor(pp)));
465  }
466  else
467  {
468  dacol->SetCycle(static_cast<int>(std::ceil(pp)));
469  }
470  dacol->Save();
471  }
472 
473  // 13. Free the memory
474  delete dacol;
475  delete mesh;
476  delete timer;
477  return 0;
478 }
Conjugate gradient method.
Definition: solvers.hpp:493
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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:324
double epsilon
Definition: ex25.cpp:140
int Dimension() const
Definition: mesh.hpp:1047
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
double RealTime()
Definition: tic_toc.cpp:426
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Data type for Gauss-Seidel smoother of sparse matrix.
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:416
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1053
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
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:9878
Timing object.
Definition: tic_toc.hpp:34
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:641
void SetTime(double t)
Set physical time (for time-dependent simulations)
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:411
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
int dim
Definition: ex24.cpp:53
void SetLevelsOfDetail(int levels_of_detail_)
IntegratorType
Definition: par_example.cpp:39
Vector data type.
Definition: vector.hpp:60
int main(int argc, char *argv[])
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
virtual void Save() override
Base class for solvers.
Definition: operator.hpp:682
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:382
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:406