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