MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
seq_example.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
37using namespace mfem;
38
45
46/// Non-linear solver for the p-Laplacian problem.
47class NLSolverPLaplacian
48{
49public:
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 real_t powerp=2,
54 Coefficient* load=nullptr,
55 real_t 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(real_t rtol)
121 {
122 newton_rtol=rtol;
123 }
124
125 // set absolute tolerance for the Newton solver
126 void SetNRATol(real_t 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(real_t rtol)
138 {
139 linear_rtol=rtol;
140 }
141
142 void SetLSATol(real_t 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 real_t GetEnergy(Vector& statev)
173 {
174 if (nlform==nullptr)
175 {
176 // allocate the solvers
177 AllocSolvers();
178 }
179 return nlform->GetEnergy(statev);
180 }
181
182
183private:
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 }
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 real_t newton_rtol;
260 real_t newton_atol;
261 int newton_iter;
262
263 real_t linear_rtol;
264 real_t 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
291int 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 real_t newton_rel_tol = 1e-4;
299 real_t newton_abs_tol = 1e-6;
300 int newton_iter = 10;
301 int print_level = 0;
302
303 real_t 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 real_t 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, (real_t)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<real_t>::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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
Conjugate gradient method.
Definition solvers.hpp:513
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
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)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition gridfunc.cpp:381
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Newton's method for solving F(x)=b for a given operator F.
Definition solvers.hpp:667
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:1810
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition solvers.hpp:714
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition solvers.cpp:1821
virtual real_t GetEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
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.
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
Base class for solvers.
Definition operator.hpp:683
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Timing object.
Definition tic_toc.hpp:36
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Stop()
Stop the stopwatch.
Definition tic_toc.cpp:422
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.
Definition tic_toc.cpp:406
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
int main()
real_t b
Definition lissajous.cpp:42
float real_t
Definition config.hpp:43
IntegratorType
@ ADJacobianIntegrator
@ HandCodedIntegrator
@ ADHessianIntegrator
@ ADJacobianIntegrator
@ HandCodedIntegrator
@ ADHessianIntegrator