MFEM  v3.4
Finite element discretization library
ex10p.cpp
Go to the documentation of this file.
1 // MFEM Example 10 - Parallel Version
2 // SUNDIALS Modification
3 //
4 // Compile with: make ex10p
5 //
6 // Sample runs:
7 // mpirun -np 4 ex10p -m ../../data/beam-quad.mesh -rp 1 -o 2 -s 5 -dt 0.15 -vs 10
8 // mpirun -np 4 ex10p -m ../../data/beam-tri.mesh -rp 1 -o 2 -s 7 -dt 0.25 -vs 10
9 // mpirun -np 4 ex10p -m ../../data/beam-hex.mesh -rp 0 -o 2 -s 5 -dt 0.15 -vs 10
10 // mpirun -np 4 ex10p -m ../../data/beam-tri.mesh -rp 1 -o 2 -s 2 -dt 3 -nls kinsol
11 // mpirun -np 4 ex10p -m ../../data/beam-quad.mesh -rp 1 -o 2 -s 2 -dt 3 -nls kinsol
12 // mpirun -np 4 ex10p -m ../../data/beam-hex.mesh -rs 1 -o 2 -s 2 -dt 3 -nls kinsol
13 // mpirun -np 4 ex10p -m ../../data/beam-quad.mesh -rp 1 -o 2 -s 15 -dt 3e-3 -vs 120
14 // mpirun -np 4 ex10p -m ../../data/beam-tri.mesh -rp 1 -o 2 -s 16 -dt 5e-3 -vs 60
15 // mpirun -np 4 ex10p -m ../../data/beam-hex.mesh -rp 0 -o 2 -s 15 -dt 5e-3 -vs 60
16 // mpirun -np 4 ex10p -m ../../data/beam-quad-amr.mesh -rp 1 -o 2 -s 5 -dt 0.15 -vs 10
17 //
18 // Description: This examples solves a time dependent nonlinear elasticity
19 // problem of the form dv/dt = H(x) + S v, dx/dt = v, where H is a
20 // hyperelastic model and S is a viscosity operator of Laplacian
21 // type. The geometry of the domain is assumed to be as follows:
22 //
23 // +---------------------+
24 // boundary --->| |
25 // attribute 1 | |
26 // (fixed) +---------------------+
27 //
28 // The example demonstrates the use of nonlinear operators (the
29 // class HyperelasticOperator defining H(x)), as well as their
30 // implicit time integration using a Newton method for solving an
31 // associated reduced backward-Euler type nonlinear equation
32 // (class ReducedSystemOperator). Each Newton step requires the
33 // inversion of a Jacobian matrix, which is done through a
34 // (preconditioned) inner solver. Note that implementing the
35 // method HyperelasticOperator::ImplicitSolve is the only
36 // requirement for high-order implicit (SDIRK) time integration.
37 //
38 // We recommend viewing examples 2 and 9 before viewing this
39 // example.
40 
41 #include "mfem.hpp"
42 #include <memory>
43 #include <iostream>
44 #include <fstream>
45 #include <string>
46 #include <map>
47 
48 #ifndef MFEM_USE_SUNDIALS
49 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
50 #endif
51 
52 using namespace std;
53 using namespace mfem;
54 
55 class ReducedSystemOperator;
56 class SundialsJacSolver;
57 
58 /** After spatial discretization, the hyperelastic model can be written as a
59  * system of ODEs:
60  * dv/dt = -M^{-1}*(H(x) + S*v)
61  * dx/dt = v,
62  * where x is the vector representing the deformation, v is the velocity field,
63  * M is the mass matrix, S is the viscosity matrix, and H(x) is the nonlinear
64  * hyperelastic operator.
65  *
66  * Class HyperelasticOperator represents the right-hand side of the above
67  * system of ODEs. */
68 class HyperelasticOperator : public TimeDependentOperator
69 {
70 protected:
71  ParFiniteElementSpace &fespace;
72  Array<int> ess_tdof_list;
73 
74  ParBilinearForm M, S;
76  double viscosity;
77  HyperelasticModel *model;
78 
79  HypreParMatrix *Mmat; // Mass matrix from ParallelAssemble()
80  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
81  HypreSmoother M_prec; // Preconditioner for the mass matrix M
82 
83  /** Nonlinear operator defining the reduced backward Euler equation for the
84  velocity. Used in the implementation of method ImplicitSolve. */
85  ReducedSystemOperator *reduced_oper;
86 
87  /// Newton solver for the reduced backward Euler equation
88  NewtonSolver *newton_solver;
89 
90  /// Solver for the Jacobian solve in the Newton method
91  Solver *J_solver;
92  /// Preconditioner for the Jacobian solve in the Newton method
93  Solver *J_prec;
94 
95  mutable Vector z; // auxiliary vector
96 
97 public:
98  /// Solver type to use in the ImplicitSolve() method, used by SDIRK methods.
99  enum NonlinearSolverType
100  {
101  NEWTON = 0, ///< Use MFEM's plain NewtonSolver
102  KINSOL = 1 ///< Use SUNDIALS' KINSOL (through MFEM's class KinSolver)
103  };
104 
105  HyperelasticOperator(ParFiniteElementSpace &f, Array<int> &ess_bdr,
106  double visc, double mu, double K,
107  NonlinearSolverType nls_type);
108 
109  /// Compute the right-hand side of the ODE system.
110  virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
111  /** Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown k.
112  This is the only requirement for high-order SDIRK implicit integration.*/
113  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
114 
115  /** Connect the Jacobian linear system solver (SundialsJacSolver) used by
116  SUNDIALS' CVODE and ARKODE time integrators to the internal objects
117  created by HyperelasticOperator. This method is called by the InitSystem
118  method of SundialsJacSolver. */
119  void InitSundialsJacSolver(SundialsJacSolver &sjsolv);
120 
121  double ElasticEnergy(const ParGridFunction &x) const;
122  double KineticEnergy(const ParGridFunction &v) const;
123  void GetElasticEnergyDensity(const ParGridFunction &x,
124  ParGridFunction &w) const;
125 
126  virtual ~HyperelasticOperator();
127 };
128 
129 /** Nonlinear operator of the form:
130  k --> (M + dt*S)*k + H(x + dt*v + dt^2*k) + S*v,
131  where M and S are given BilinearForms, H is a given NonlinearForm, v and x
132  are given vectors, and dt is a scalar. */
133 class ReducedSystemOperator : public Operator
134 {
135 private:
136  ParBilinearForm *M, *S;
137  ParNonlinearForm *H;
138  mutable HypreParMatrix *Jacobian;
139  double dt;
140  const Vector *v, *x;
141  mutable Vector w, z;
142  const Array<int> &ess_tdof_list;
143 
144 public:
145  ReducedSystemOperator(ParBilinearForm *M_, ParBilinearForm *S_,
146  ParNonlinearForm *H_, const Array<int> &ess_tdof_list);
147 
148  /// Set current dt, v, x values - needed to compute action and Jacobian.
149  void SetParameters(double dt_, const Vector *v_, const Vector *x_);
150 
151  /// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
152  virtual void Mult(const Vector &k, Vector &y) const;
153 
154  /// Compute J = M + dt S + dt^2 grad_H(x + dt (v + dt k)).
155  virtual Operator &GetGradient(const Vector &k) const;
156 
157  virtual ~ReducedSystemOperator();
158 };
159 
160 /// Custom Jacobian system solver for the SUNDIALS time integrators.
161 /** For the ODE system represented by HyperelasticOperator
162 
163  M dv/dt = -(H(x) + S*v)
164  dx/dt = v,
165 
166  this class facilitates the solution of linear systems of the form
167 
168  (M + γS) yv + γJ yx = M bv, J=(dH/dx)(x)
169  - γ yv + yx = bx
170 
171  for given bv, bx, x, and γ = GetTimeStep(). */
172 class SundialsJacSolver : public SundialsODELinearSolver
173 {
174 private:
175  ParBilinearForm *M, *S;
176  ParNonlinearForm *H;
177  const SparseMatrix *local_grad_H;
178  HypreParMatrix *Jacobian;
179  Solver *J_solver;
180  const Array<int> *ess_tdof_list;
181 
182 public:
183  SundialsJacSolver()
184  : M(), S(), H(), local_grad_H(), Jacobian(), J_solver() { }
185 
186  /// Connect the solver to the objects created inside HyperelasticOperator.
187  void SetOperators(ParBilinearForm &M_, ParBilinearForm &S_,
188  ParNonlinearForm &H_, Solver &solver,
189  const Array<int> &ess_tdof_list_)
190  {
191  M = &M_; S = &S_; H = &H_; J_solver = &solver;
192  ess_tdof_list = &ess_tdof_list_;
193  }
194 
195  /** Linear solve applicable to the SUNDIALS format.
196  Solves (Mass - dt J) y = Mass b, where in our case:
197  Mass = | M 0 | J = | -S -grad_H | y = | v_hat | b = | b_v |
198  | 0 I | | I 0 | | x_hat | | b_x |
199  The result replaces the rhs b.
200  We substitute x_hat = b_x + dt v_hat and solve
201  (M + dt S + dt^2 grad_H) v_hat = M b_v - dt grad_H b_x. */
202  int InitSystem(void *sundials_mem);
203  int SetupSystem(void *sundials_mem, int conv_fail,
204  const Vector &y_pred, const Vector &f_pred, int &jac_cur,
205  Vector &v_temp1, Vector &v_temp2, Vector &v_temp3);
206  int SolveSystem(void *sundials_mem, Vector &b, const Vector &weight,
207  const Vector &y_cur, const Vector &f_cur);
208  int FreeSystem(void *sundials_mem);
209 };
210 
211 
212 /** Function representing the elastic energy density for the given hyperelastic
213  model+deformation. Used in HyperelasticOperator::GetElasticEnergyDensity. */
214 class ElasticEnergyCoefficient : public Coefficient
215 {
216 private:
217  HyperelasticModel &model;
218  const ParGridFunction &x;
219  DenseMatrix J;
220 
221 public:
222  ElasticEnergyCoefficient(HyperelasticModel &m, const ParGridFunction &x_)
223  : model(m), x(x_) { }
224  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
225  virtual ~ElasticEnergyCoefficient() { }
226 };
227 
228 void InitialDeformation(const Vector &x, Vector &y);
229 
230 void InitialVelocity(const Vector &x, Vector &v);
231 
232 void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
233  ParGridFunction *field, const char *field_name = NULL,
234  bool init_vis = false);
235 
236 
237 int main(int argc, char *argv[])
238 {
239  // 1. Initialize MPI.
240  int num_procs, myid;
241  MPI_Init(&argc, &argv);
242  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
243  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
244 
245  // 2. Parse command-line options.
246  const char *mesh_file = "../../data/beam-quad.mesh";
247  int ser_ref_levels = 2;
248  int par_ref_levels = 0;
249  int order = 2;
250  int ode_solver_type = 3;
251  double t_final = 300.0;
252  double dt = 3.0;
253  double visc = 1e-2;
254  double mu = 0.25;
255  double K = 5.0;
256  bool visualization = true;
257  const char *nls = "newton";
258  int vis_steps = 1;
259 
260  // Relative and absolute tolerances for CVODE and ARKODE.
261  const double reltol = 1e-1, abstol = 1e-1;
262 
263  OptionsParser args(argc, argv);
264  args.AddOption(&mesh_file, "-m", "--mesh",
265  "Mesh file to use.");
266  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
267  "Number of times to refine the mesh uniformly in serial.");
268  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
269  "Number of times to refine the mesh uniformly in parallel.");
270  args.AddOption(&order, "-o", "--order",
271  "Order (degree) of the finite elements.");
272  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
273  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
274  " 4 - CVODE implicit, approximate Jacobian,\n\t"
275  " 5 - CVODE implicit, specified Jacobian,\n\t"
276  " 6 - ARKODE implicit, approximate Jacobian,\n\t"
277  " 7 - ARKODE implicit, specified Jacobian,\n\t"
278  " 11 - Forward Euler, 12 - RK2,\n\t"
279  " 13 - RK3 SSP, 14 - RK4,\n\t"
280  " 15 - CVODE (adaptive order) explicit,\n\t"
281  " 16 - ARKODE default (4th order) explicit.");
282  args.AddOption(&nls, "-nls", "--nonlinear-solver",
283  "Nonlinear systems solver: "
284  "\"newton\" (plain Newton) or \"kinsol\" (KINSOL).");
285  args.AddOption(&t_final, "-tf", "--t-final",
286  "Final time; start time is 0.");
287  args.AddOption(&dt, "-dt", "--time-step",
288  "Time step.");
289  args.AddOption(&visc, "-v", "--viscosity",
290  "Viscosity coefficient.");
291  args.AddOption(&mu, "-mu", "--shear-modulus",
292  "Shear modulus in the Neo-Hookean hyperelastic model.");
293  args.AddOption(&K, "-K", "--bulk-modulus",
294  "Bulk modulus in the Neo-Hookean hyperelastic model.");
295  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
296  "--no-visualization",
297  "Enable or disable GLVis visualization.");
298  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
299  "Visualize every n-th timestep.");
300  args.Parse();
301  if (!args.Good())
302  {
303  if (myid == 0)
304  {
305  args.PrintUsage(cout);
306  }
307  MPI_Finalize();
308  return 1;
309  }
310  if (myid == 0)
311  {
312  args.PrintOptions(cout);
313  }
314 
315  // 3. Read the serial mesh from the given mesh file on all processors. We can
316  // handle triangular, quadrilateral, tetrahedral and hexahedral meshes
317  // with the same code.
318  Mesh *mesh = new Mesh(mesh_file, 1, 1);
319  int dim = mesh->Dimension();
320 
321  // 4. Define the ODE solver used for time integration. Several implicit
322  // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
323  // explicit Runge-Kutta methods are available.
324  ODESolver *ode_solver;
325  CVODESolver *cvode = NULL;
326  ARKODESolver *arkode = NULL;
327  SundialsJacSolver *sjsolver = NULL;
328  switch (ode_solver_type)
329  {
330  // Implicit L-stable methods
331  case 1: ode_solver = new BackwardEulerSolver; break;
332  case 2: ode_solver = new SDIRK23Solver(2); break;
333  case 3: ode_solver = new SDIRK33Solver; break;
334  case 4:
335  case 5:
336  cvode = new CVODESolver(MPI_COMM_WORLD, CV_BDF, CV_NEWTON);
337  cvode->SetSStolerances(reltol, abstol);
338  cvode->SetMaxStep(dt);
339  if (ode_solver_type == 5)
340  {
341  sjsolver = new SundialsJacSolver;
342  cvode->SetLinearSolver(*sjsolver); // Custom Jacobian inversion.
343  }
344  ode_solver = cvode; break;
345  case 6:
346  case 7:
347  arkode = new ARKODESolver(MPI_COMM_WORLD, ARKODESolver::IMPLICIT);
348  arkode->SetSStolerances(reltol, abstol);
349  arkode->SetMaxStep(dt);
350  if (ode_solver_type == 7)
351  {
352  sjsolver = new SundialsJacSolver;
353  arkode->SetLinearSolver(*sjsolver); // Custom Jacobian inversion.
354  }
355  ode_solver = arkode; break;
356  // Explicit methods
357  case 11: ode_solver = new ForwardEulerSolver; break;
358  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
359  case 13: ode_solver = new RK3SSPSolver; break;
360  case 14: ode_solver = new RK4Solver; break;
361  case 15:
362  cvode = new CVODESolver(MPI_COMM_WORLD, CV_ADAMS, CV_FUNCTIONAL);
363  cvode->SetSStolerances(reltol, abstol);
364  cvode->SetMaxStep(dt);
365  ode_solver = cvode; break;
366  case 16:
367  arkode = new ARKODESolver(MPI_COMM_WORLD, ARKODESolver::EXPLICIT);
368  arkode->SetSStolerances(reltol, abstol);
369  arkode->SetMaxStep(dt);
370  ode_solver = arkode; break;
371  // Implicit A-stable methods (not L-stable)
372  case 22: ode_solver = new ImplicitMidpointSolver; break;
373  case 23: ode_solver = new SDIRK23Solver; break;
374  case 24: ode_solver = new SDIRK34Solver; break;
375  default:
376  if (myid == 0)
377  {
378  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
379  }
380  delete mesh;
381  MPI_Finalize();
382  return 3;
383  }
384 
385  map<string,HyperelasticOperator::NonlinearSolverType> nls_map;
386  nls_map["newton"] = HyperelasticOperator::NEWTON;
387  nls_map["kinsol"] = HyperelasticOperator::KINSOL;
388  if (nls_map.find(nls) == nls_map.end())
389  {
390  if (myid == 0)
391  {
392  cout << "Unknown type of nonlinear solver: " << nls << endl;
393  }
394  delete ode_solver;
395  delete mesh;
396  MPI_Finalize();
397  return 4;
398  }
399 
400  // 5. Refine the mesh in serial to increase the resolution. In this example
401  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
402  // a command-line parameter.
403  for (int lev = 0; lev < ser_ref_levels; lev++)
404  {
405  mesh->UniformRefinement();
406  }
407 
408  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
409  // this mesh further in parallel to increase the resolution. Once the
410  // parallel mesh is defined, the serial mesh can be deleted.
411  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
412  delete mesh;
413  for (int lev = 0; lev < par_ref_levels; lev++)
414  {
415  pmesh->UniformRefinement();
416  }
417 
418  // 7. Define the parallel vector finite element spaces representing the mesh
419  // deformation x_gf, the velocity v_gf, and the initial configuration,
420  // x_ref. Define also the elastic energy density, w_gf, which is in a
421  // discontinuous higher-order space. Since x and v are integrated in time
422  // as a system, we group them together in block vector vx, on the unique
423  // parallel degrees of freedom, with offsets given by array true_offset.
424  H1_FECollection fe_coll(order, dim);
425  ParFiniteElementSpace fespace(pmesh, &fe_coll, dim);
426 
427  HYPRE_Int glob_size = fespace.GlobalTrueVSize();
428  if (myid == 0)
429  {
430  cout << "Number of velocity/deformation unknowns: " << glob_size << endl;
431  }
432  int true_size = fespace.TrueVSize();
433  Array<int> true_offset(3);
434  true_offset[0] = 0;
435  true_offset[1] = true_size;
436  true_offset[2] = 2*true_size;
437 
438  BlockVector vx(true_offset);
439  ParGridFunction v_gf, x_gf;
440  v_gf.MakeTRef(&fespace, vx, true_offset[0]);
441  x_gf.MakeTRef(&fespace, vx, true_offset[1]);
442 
443  ParGridFunction x_ref(&fespace);
444  pmesh->GetNodes(x_ref);
445 
446  L2_FECollection w_fec(order + 1, dim);
447  ParFiniteElementSpace w_fespace(pmesh, &w_fec);
448  ParGridFunction w_gf(&w_fespace);
449 
450  // 8. Set the initial conditions for v_gf, x_gf and vx, and define the
451  // boundary conditions on a beam-like mesh (see description above).
453  v_gf.ProjectCoefficient(velo);
454  v_gf.SetTrueVector();
456  x_gf.ProjectCoefficient(deform);
457  x_gf.SetTrueVector();
458 
459  v_gf.SetFromTrueVector(); x_gf.SetFromTrueVector();
460 
461  Array<int> ess_bdr(fespace.GetMesh()->bdr_attributes.Max());
462  ess_bdr = 0;
463  ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed
464 
465  // 9. Initialize the hyperelastic operator, the GLVis visualization and print
466  // the initial energies.
467  HyperelasticOperator oper(fespace, ess_bdr, visc, mu, K, nls_map[nls]);
468 
469  socketstream vis_v, vis_w;
470  if (visualization)
471  {
472  char vishost[] = "localhost";
473  int visport = 19916;
474  vis_v.open(vishost, visport);
475  vis_v.precision(8);
476  visualize(vis_v, pmesh, &x_gf, &v_gf, "Velocity", true);
477  // Make sure all ranks have sent their 'v' solution before initiating
478  // another set of GLVis connections (one from each rank):
479  MPI_Barrier(pmesh->GetComm());
480  vis_w.open(vishost, visport);
481  if (vis_w)
482  {
483  oper.GetElasticEnergyDensity(x_gf, w_gf);
484  vis_w.precision(8);
485  visualize(vis_w, pmesh, &x_gf, &w_gf, "Elastic energy density", true);
486  }
487  }
488 
489  double ee0 = oper.ElasticEnergy(x_gf);
490  double ke0 = oper.KineticEnergy(v_gf);
491  if (myid == 0)
492  {
493  cout << "initial elastic energy (EE) = " << ee0 << endl;
494  cout << "initial kinetic energy (KE) = " << ke0 << endl;
495  cout << "initial total energy (TE) = " << (ee0 + ke0) << endl;
496  }
497 
498  double t = 0.0;
499  oper.SetTime(t);
500  ode_solver->Init(oper);
501 
502  // 10. Perform time-integration
503  // (looping over the time iterations, ti, with a time-step dt).
504  bool last_step = false;
505  for (int ti = 1; !last_step; ti++)
506  {
507  double dt_real = min(dt, t_final - t);
508 
509  ode_solver->Step(vx, t, dt_real);
510 
511  last_step = (t >= t_final - 1e-8*dt);
512 
513  if (last_step || (ti % vis_steps) == 0)
514  {
515  v_gf.SetFromTrueVector(); x_gf.SetFromTrueVector();
516 
517  double ee = oper.ElasticEnergy(x_gf);
518  double ke = oper.KineticEnergy(v_gf);
519 
520  if (myid == 0)
521  {
522  cout << "step " << ti << ", t = " << t << ", EE = " << ee
523  << ", KE = " << ke << ", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
524 
525  if (cvode) { cvode->PrintInfo(); }
526  else if (arkode) { arkode->PrintInfo(); }
527  }
528 
529  if (visualization)
530  {
531  visualize(vis_v, pmesh, &x_gf, &v_gf);
532  if (vis_w)
533  {
534  oper.GetElasticEnergyDensity(x_gf, w_gf);
535  visualize(vis_w, pmesh, &x_gf, &w_gf);
536  }
537  }
538  }
539  }
540 
541  // 11. Save the displaced mesh, the velocity and elastic energy.
542  {
543  v_gf.SetFromTrueVector(); x_gf.SetFromTrueVector();
544  GridFunction *nodes = &x_gf;
545  int owns_nodes = 0;
546  pmesh->SwapNodes(nodes, owns_nodes);
547 
548  ostringstream mesh_name, velo_name, ee_name;
549  mesh_name << "deformed." << setfill('0') << setw(6) << myid;
550  velo_name << "velocity." << setfill('0') << setw(6) << myid;
551  ee_name << "elastic_energy." << setfill('0') << setw(6) << myid;
552 
553  ofstream mesh_ofs(mesh_name.str().c_str());
554  mesh_ofs.precision(8);
555  pmesh->Print(mesh_ofs);
556  pmesh->SwapNodes(nodes, owns_nodes);
557  ofstream velo_ofs(velo_name.str().c_str());
558  velo_ofs.precision(8);
559  v_gf.Save(velo_ofs);
560  ofstream ee_ofs(ee_name.str().c_str());
561  ee_ofs.precision(8);
562  oper.GetElasticEnergyDensity(x_gf, w_gf);
563  w_gf.Save(ee_ofs);
564  }
565 
566  // 12. Free the used memory.
567  delete ode_solver;
568  delete sjsolver;
569  delete pmesh;
570 
571  MPI_Finalize();
572 
573  return 0;
574 }
575 
576 void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
577  ParGridFunction *field, const char *field_name, bool init_vis)
578 {
579  if (!out)
580  {
581  return;
582  }
583 
584  GridFunction *nodes = deformed_nodes;
585  int owns_nodes = 0;
586 
587  mesh->SwapNodes(nodes, owns_nodes);
588 
589  out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
590  out << "solution\n" << *mesh << *field;
591 
592  mesh->SwapNodes(nodes, owns_nodes);
593 
594  if (init_vis)
595  {
596  out << "window_size 800 800\n";
597  out << "window_title '" << field_name << "'\n";
598  if (mesh->SpaceDimension() == 2)
599  {
600  out << "view 0 0\n"; // view from top
601  out << "keys jl\n"; // turn off perspective and light
602  }
603  out << "keys cm\n"; // show colorbar and mesh
604  out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
605  out << "pause\n";
606  }
607  out << flush;
608 }
609 
610 
611 ReducedSystemOperator::ReducedSystemOperator(
613  const Array<int> &ess_tdof_list_)
614  : Operator(M_->ParFESpace()->TrueVSize()), M(M_), S(S_), H(H_),
615  Jacobian(NULL), dt(0.0), v(NULL), x(NULL), w(height), z(height),
616  ess_tdof_list(ess_tdof_list_)
617 { }
618 
619 void ReducedSystemOperator::SetParameters(double dt_, const Vector *v_,
620  const Vector *x_)
621 {
622  dt = dt_; v = v_; x = x_;
623 }
624 
625 void ReducedSystemOperator::Mult(const Vector &k, Vector &y) const
626 {
627  // compute: y = H(x + dt*(v + dt*k)) + M*k + S*(v + dt*k)
628  add(*v, dt, k, w);
629  add(*x, dt, w, z);
630  H->Mult(z, y);
631  M->TrueAddMult(k, y);
632  S->TrueAddMult(w, y);
633  y.SetSubVector(ess_tdof_list, 0.0);
634 }
635 
636 Operator &ReducedSystemOperator::GetGradient(const Vector &k) const
637 {
638  delete Jacobian;
639  SparseMatrix *localJ = Add(1.0, M->SpMat(), dt, S->SpMat());
640  add(*v, dt, k, w);
641  add(*x, dt, w, z);
642  localJ->Add(dt*dt, H->GetLocalGradient(z));
643  Jacobian = M->ParallelAssemble(localJ);
644  delete localJ;
645  HypreParMatrix *Je = Jacobian->EliminateRowsCols(ess_tdof_list);
646  delete Je;
647  return *Jacobian;
648 }
649 
650 ReducedSystemOperator::~ReducedSystemOperator()
651 {
652  delete Jacobian;
653 }
654 
655 
656 int SundialsJacSolver::InitSystem(void *sundials_mem)
657 {
658  TimeDependentOperator *td_oper = GetTimeDependentOperator(sundials_mem);
659  HyperelasticOperator *he_oper;
660 
661  // During development, we use dynamic_cast<> to ensure the setup is correct:
662  he_oper = dynamic_cast<HyperelasticOperator*>(td_oper);
663  MFEM_VERIFY(he_oper, "operator is not HyperelasticOperator");
664 
665  // When the implementation is finalized, we can switch to static_cast<>:
666  // he_oper = static_cast<HyperelasticOperator*>(td_oper);
667 
668  he_oper->InitSundialsJacSolver(*this);
669  return 0;
670 }
671 
672 int SundialsJacSolver::SetupSystem(void *sundials_mem, int conv_fail,
673  const Vector &y_pred, const Vector &f_pred,
674  int &jac_cur, Vector &v_temp1,
675  Vector &v_temp2, Vector &v_temp3)
676 {
677  int sc = y_pred.Size() / 2;
678  const Vector x(y_pred.GetData() + sc, sc);
679  double dt = GetTimeStep(sundials_mem);
680 
681  // J = M + dt*(S + dt*grad(H))
682  delete Jacobian;
683  SparseMatrix *localJ = Add(1.0, M->SpMat(), dt, S->SpMat());
684  local_grad_H = &H->GetLocalGradient(x);
685  localJ->Add(dt*dt, *local_grad_H);
686  Jacobian = M->ParallelAssemble(localJ);
687  delete localJ;
688  HypreParMatrix *Je = Jacobian->EliminateRowsCols(*ess_tdof_list);
689  delete Je;
690 
691  J_solver->SetOperator(*Jacobian);
692 
693  jac_cur = 1;
694  return 0;
695 }
696 
697 int SundialsJacSolver::SolveSystem(void *sundials_mem, Vector &b,
698  const Vector &weight, const Vector &y_cur,
699  const Vector &f_cur)
700 {
701  int sc = b.Size() / 2;
702  ParFiniteElementSpace *fes = H->ParFESpace();
703  // Vector x(y_cur.GetData() + sc, sc);
704  Vector b_v(b.GetData() + 0, sc);
705  Vector b_x(b.GetData() + sc, sc);
706  Vector rhs(sc);
707  double dt = GetTimeStep(sundials_mem);
708 
709  // We can assume that b_v and b_x have zeros at essential tdofs.
710 
711  // rhs = M b_v - dt*grad(H) b_x
712  ParGridFunction lb_x(fes), lrhs(fes);
713  lb_x.Distribute(b_x);
714  local_grad_H->Mult(lb_x, lrhs);
715  lrhs.ParallelAssemble(rhs);
716  rhs *= -dt;
717  M->TrueAddMult(b_v, rhs);
718  rhs.SetSubVector(*ess_tdof_list, 0.0);
719 
720  J_solver->iterative_mode = false;
721  J_solver->Mult(rhs, b_v);
722 
723  b_x.Add(dt, b_v);
724 
725  return 0;
726 }
727 
728 int SundialsJacSolver::FreeSystem(void *sundials_mem)
729 {
730  delete Jacobian;
731  return 0;
732 }
733 
734 
735 HyperelasticOperator::HyperelasticOperator(ParFiniteElementSpace &f,
736  Array<int> &ess_bdr, double visc,
737  double mu, double K,
738  NonlinearSolverType nls_type)
739  : TimeDependentOperator(2*f.TrueVSize(), 0.0), fespace(f),
740  M(&fespace), S(&fespace), H(&fespace),
741  viscosity(visc), M_solver(f.GetComm()), z(height/2)
742 {
743  const double rel_tol = 1e-8;
744  const int skip_zero_entries = 0;
745 
746  const double ref_density = 1.0; // density in the reference configuration
747  ConstantCoefficient rho0(ref_density);
748  M.AddDomainIntegrator(new VectorMassIntegrator(rho0));
749  M.Assemble(skip_zero_entries);
750  M.Finalize(skip_zero_entries);
751  Mmat = M.ParallelAssemble();
752  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
753  HypreParMatrix *Me = Mmat->EliminateRowsCols(ess_tdof_list);
754  delete Me;
755 
756  M_solver.iterative_mode = false;
757  M_solver.SetRelTol(rel_tol);
758  M_solver.SetAbsTol(0.0);
759  M_solver.SetMaxIter(30);
760  M_solver.SetPrintLevel(0);
761  M_prec.SetType(HypreSmoother::Jacobi);
762  M_solver.SetPreconditioner(M_prec);
763  M_solver.SetOperator(*Mmat);
764 
765  model = new NeoHookeanModel(mu, K);
766  H.AddDomainIntegrator(new HyperelasticNLFIntegrator(model));
767  H.SetEssentialTrueDofs(ess_tdof_list);
768 
769  ConstantCoefficient visc_coeff(viscosity);
770  S.AddDomainIntegrator(new VectorDiffusionIntegrator(visc_coeff));
771  S.Assemble(skip_zero_entries);
772  S.Finalize(skip_zero_entries);
773 
774  reduced_oper = new ReducedSystemOperator(&M, &S, &H, ess_tdof_list);
775 
776  HypreSmoother *J_hypreSmoother = new HypreSmoother;
777  J_hypreSmoother->SetType(HypreSmoother::l1Jacobi);
778  J_hypreSmoother->SetPositiveDiagonal(true);
779  J_prec = J_hypreSmoother;
780 
781  MINRESSolver *J_minres = new MINRESSolver(f.GetComm());
782  J_minres->SetRelTol(rel_tol);
783  J_minres->SetAbsTol(0.0);
784  J_minres->SetMaxIter(300);
785  J_minres->SetPrintLevel(-1);
786  J_minres->SetPreconditioner(*J_prec);
787  J_solver = J_minres;
788 
789  if (nls_type == KINSOL)
790  {
791  KinSolver *kinsolver = new KinSolver(f.GetComm(), KIN_NONE, true);
792  kinsolver->SetMaxSetupCalls(4);
793  newton_solver = kinsolver;
794  newton_solver->SetMaxIter(200);
795  newton_solver->SetRelTol(rel_tol);
796  newton_solver->SetPrintLevel(0);
797  }
798  else
799  {
800  newton_solver = new NewtonSolver(f.GetComm());
801  newton_solver->SetMaxIter(10);
802  newton_solver->SetRelTol(rel_tol);
803  newton_solver->SetPrintLevel(-1);
804  }
805  newton_solver->SetSolver(*J_solver);
806  newton_solver->iterative_mode = false;
807  newton_solver->SetOperator(*reduced_oper);
808 }
809 
810 void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
811 {
812  // Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt
813  int sc = height/2;
814  Vector v(vx.GetData() + 0, sc);
815  Vector x(vx.GetData() + sc, sc);
816  Vector dv_dt(dvx_dt.GetData() + 0, sc);
817  Vector dx_dt(dvx_dt.GetData() + sc, sc);
818 
819  H.Mult(x, z);
820  if (viscosity != 0.0)
821  {
822  S.TrueAddMult(v, z);
823  z.SetSubVector(ess_tdof_list, 0.0);
824  }
825  z.Neg(); // z = -z
826  M_solver.Mult(z, dv_dt);
827 
828  dx_dt = v;
829 }
830 
831 void HyperelasticOperator::ImplicitSolve(const double dt,
832  const Vector &vx, Vector &dvx_dt)
833 {
834  int sc = height/2;
835  Vector v(vx.GetData() + 0, sc);
836  Vector x(vx.GetData() + sc, sc);
837  Vector dv_dt(dvx_dt.GetData() + 0, sc);
838  Vector dx_dt(dvx_dt.GetData() + sc, sc);
839 
840  // By eliminating kx from the coupled system:
841  // kv = -M^{-1}*[H(x + dt*kx) + S*(v + dt*kv)]
842  // kx = v + dt*kv
843  // we reduce it to a nonlinear equation for kv, represented by the
844  // reduced_oper. This equation is solved with the newton_solver
845  // object (using J_solver and J_prec internally).
846  reduced_oper->SetParameters(dt, &v, &x);
847  Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
848  newton_solver->Mult(zero, dv_dt);
849  MFEM_VERIFY(newton_solver->GetConverged(),
850  "Nonlinear solver did not converge.");
851 #ifdef MFEM_DEBUG
852  if (fespace.GetMyRank() == 0)
853  {
854  cout << " num nonlin sol iters = " << newton_solver->GetNumIterations()
855  << ", final norm = " << newton_solver->GetFinalNorm() << '\n';
856  }
857 #endif
858  add(v, dt, dv_dt, dx_dt);
859 }
860 
861 void HyperelasticOperator::InitSundialsJacSolver(SundialsJacSolver &sjsolv)
862 {
863  sjsolv.SetOperators(M, S, H, *J_solver, ess_tdof_list);
864 }
865 
866 double HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const
867 {
868  return H.GetEnergy(x);
869 }
870 
871 double HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
872 {
873  double loc_energy = 0.5*M.InnerProduct(v, v);
874  double energy;
875  MPI_Allreduce(&loc_energy, &energy, 1, MPI_DOUBLE, MPI_SUM,
876  fespace.GetComm());
877  return energy;
878 }
879 
880 void HyperelasticOperator::GetElasticEnergyDensity(
881  const ParGridFunction &x, ParGridFunction &w) const
882 {
883  ElasticEnergyCoefficient w_coeff(*model, x);
884  w.ProjectCoefficient(w_coeff);
885 }
886 
887 HyperelasticOperator::~HyperelasticOperator()
888 {
889  delete newton_solver;
890  delete J_solver;
891  delete J_prec;
892  delete reduced_oper;
893  delete model;
894  delete Mmat;
895 }
896 
897 
899  const IntegrationPoint &ip)
900 {
901  model.SetTransformation(T);
902  x.GetVectorGradient(T, J);
903  // return model.EvalW(J); // in reference configuration
904  return model.EvalW(J)/J.Det(); // in deformed configuration
905 }
906 
907 
908 void InitialDeformation(const Vector &x, Vector &y)
909 {
910  // set the initial configuration to be the same as the reference, stress
911  // free, configuration
912  y = x;
913 }
914 
915 void InitialVelocity(const Vector &x, Vector &v)
916 {
917  const int dim = x.Size();
918  const double s = 0.1/64.;
919 
920  v = 0.0;
921  v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
922  v(0) = -s*x(0)*x(0);
923 }
void InitialDeformation(const Vector &x, Vector &y)
Definition: ex10p.cpp:661
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:1368
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:494
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Definition: coefficient.hpp:45
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:2035
Conjugate gradient method.
Definition: solvers.hpp:111
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:224
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
Definition: sundials.cpp:980
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:122
Subclass constant coefficient.
Definition: coefficient.hpp:57
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.hpp:37
int Dimension() const
Definition: mesh.hpp:645
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Definition: gridfunc.cpp:1040
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Base abstract class for time dependent operators.
Definition: operator.hpp:151
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
Definition: mesh.cpp:5460
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:109
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Parallel non-linear operator on the true dofs.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Definition: sundials.cpp:233
bool Good() const
Definition: optparser.hpp:120
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int main(int argc, char *argv[])
Definition: ex10p.cpp:161
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:290
STL namespace.
MINRES method.
Definition: solvers.hpp:221
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:3427
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:212
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:260
void SetSStolerances(double reltol, double abstol)
Specify the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:477
void SetPositiveDiagonal(bool pos=true)
After computing l1-norms, replace them with their absolute values.
Definition: hypre.hpp:605
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:116
int dim
Definition: ex3.cpp:47
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:234
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:72
Wrapper for SUNDIALS&#39; KINSOL library – Nonlinear solvers.
Definition: sundials.hpp:312
Data type sparse matrix.
Definition: sparsemat.hpp:38
Wrapper for SUNDIALS&#39; CVODE library – Multi-step time integration.
Definition: sundials.hpp:133
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:247
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6741
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition: gridfunc.cpp:199
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
void SetMaxStep(double dt_max)
Set the maximum time step of the linear multistep method.
Definition: sundials.hpp:183
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:259
Parallel smoothers in hypre.
Definition: hypre.hpp:517
MPI_Comm GetComm() const
Definition: pmesh.hpp:124
Wrapper for SUNDIALS&#39; ARKODE library – Runge-Kutta time integration.
Definition: sundials.hpp:220
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Set a custom Jacobian system solver for implicit methods.
Definition: sundials.cpp:486
int GetMyRank() const
Definition: pmesh.hpp:126
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:224
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
void SetAbsTol(double atol)
Definition: solvers.hpp:62
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:174
void SetRelTol(double rtol)
Definition: solvers.hpp:61
void PrintInfo() const
Print CVODE statistics.
Definition: sundials.cpp:400
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
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
void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes, ParGridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10p.cpp:434
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:134
void InitialVelocity(const Vector &x, Vector &v)
Definition: ex10p.cpp:668
Abstract base class, wrapping the custom linear solvers interface in SUNDIALS&#39; CVODE and ARKODE solve...
Definition: sundials.hpp:47
int SpaceDimension() const
Definition: mesh.hpp:646
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:398
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:225
Class for integration point with weight.
Definition: intrules.hpp:25
int GetNRanks() const
Definition: pmesh.hpp:125
Class for parallel bilinear form.
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
MPI_Comm GetComm() const
Definition: pfespace.hpp:235
Vector data type.
Definition: vector.hpp:48
void PrintInfo() const
Print ARKODE statistics.
Definition: sundials.cpp:691
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:79
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5422
Base class for solvers.
Definition: operator.hpp:268
Class for parallel grid function.
Definition: pgridfunc.hpp:32
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
The classical forward Euler method.
Definition: ode.hpp:101
void SetMaxStep(double dt_max)
Set the maximum time step of the Runge-Kutta method.
Definition: sundials.hpp:276
Abstract operator.
Definition: operator.hpp:21
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:1816
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:128