MFEM  v3.4
Finite element discretization library
ex10.cpp
Go to the documentation of this file.
1 // MFEM Example 10
2 // SUNDIALS Modification
3 //
4 // Compile with: make ex10
5 //
6 // Sample runs:
7 // ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 5 -dt 0.15 -vs 10
8 // ex10 -m ../../data/beam-tri.mesh -r 2 -o 2 -s 7 -dt 0.3 -vs 5
9 // ex10 -m ../../data/beam-hex.mesh -r 1 -o 2 -s 5 -dt 0.2 -vs 5
10 // ex10 -m ../../data/beam-tri.mesh -r 2 -o 2 -s 2 -dt 3 -nls kinsol
11 // ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 2 -dt 3 -nls kinsol
12 // ex10 -m ../../data/beam-hex.mesh -r 1 -o 2 -s 2 -dt 3 -nls kinsol
13 // ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 15 -dt 5e-3 -vs 60
14 // ex10 -m ../../data/beam-tri.mesh -r 2 -o 2 -s 16 -dt 0.01 -vs 30
15 // ex10 -m ../../data/beam-hex.mesh -r 1 -o 2 -s 15 -dt 0.01 -vs 30
16 // ex10 -m ../../data/beam-quad-amr.mesh -r 2 -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  FiniteElementSpace &fespace;
72 
73  BilinearForm M, S;
74  NonlinearForm H;
75  double viscosity;
76  HyperelasticModel *model;
77 
78  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
79  DSmoother 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(FiniteElementSpace &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(const Vector &x) const;
120  double KineticEnergy(const Vector &v) const;
121  void GetElasticEnergyDensity(const GridFunction &x, GridFunction &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  BilinearForm *M, *S;
134  NonlinearForm *H;
135  mutable SparseMatrix *Jacobian;
136  double dt;
137  const Vector *v, *x;
138  mutable Vector w, z;
139 
140 public:
141  ReducedSystemOperator(BilinearForm *M_, BilinearForm *S_, NonlinearForm *H_);
142 
143  /// Set current dt, v, x values - needed to compute action and Jacobian.
144  void SetParameters(double dt_, const Vector *v_, const Vector *x_);
145 
146  /// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
147  virtual void Mult(const Vector &k, Vector &y) const;
148 
149  /// Compute J = M + dt S + dt^2 grad_H(x + dt (v + dt k)).
150  virtual Operator &GetGradient(const Vector &k) const;
151 
152  virtual ~ReducedSystemOperator();
153 };
154 
155 /// Custom Jacobian system solver for the SUNDIALS time integrators.
156 /** For the ODE system represented by HyperelasticOperator
157 
158  M dv/dt = -(H(x) + S*v)
159  dx/dt = v,
160 
161  this class facilitates the solution of linear systems of the form
162 
163  (M + γS) yv + γJ yx = M bv, J=(dH/dx)(x)
164  - γ yv + yx = bx
165 
166  for given bv, bx, x, and γ = GetTimeStep(). */
167 class SundialsJacSolver : public SundialsODELinearSolver
168 {
169 private:
170  BilinearForm *M, *S;
171  NonlinearForm *H;
172  SparseMatrix *grad_H, *Jacobian;
173  Solver *J_solver;
174 
175 public:
176  SundialsJacSolver()
177  : M(), S(), H(), grad_H(), Jacobian(), J_solver() { }
178 
179  /// Connect the solver to the objects created inside HyperelasticOperator.
180  void SetOperators(BilinearForm &M_, BilinearForm &S_,
181  NonlinearForm &H_, Solver &solver)
182  {
183  M = &M_; S = &S_; H = &H_; J_solver = &solver;
184  }
185 
186  /** Linear solve applicable to the SUNDIALS format.
187  Solves (Mass - dt J) y = Mass b, where in our case:
188  Mass = | M 0 | J = | -S -grad_H | y = | v_hat | b = | b_v |
189  | 0 I | | I 0 | | x_hat | | b_x |
190  The result replaces the rhs b.
191  We substitute x_hat = b_x + dt v_hat and solve
192  (M + dt S + dt^2 grad_H) v_hat = M b_v - dt grad_H b_x. */
193  int InitSystem(void *sundials_mem);
194  int SetupSystem(void *sundials_mem, int conv_fail,
195  const Vector &y_pred, const Vector &f_pred, int &jac_cur,
196  Vector &v_temp1, Vector &v_temp2, Vector &v_temp3);
197  int SolveSystem(void *sundials_mem, Vector &b, const Vector &weight,
198  const Vector &y_cur, const Vector &f_cur);
199  int FreeSystem(void *sundials_mem);
200 };
201 
202 
203 /** Function representing the elastic energy density for the given hyperelastic
204  model+deformation. Used in HyperelasticOperator::GetElasticEnergyDensity. */
205 class ElasticEnergyCoefficient : public Coefficient
206 {
207 private:
208  HyperelasticModel &model;
209  const GridFunction &x;
210  DenseMatrix J;
211 
212 public:
213  ElasticEnergyCoefficient(HyperelasticModel &m, const GridFunction &x_)
214  : model(m), x(x_) { }
215  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
216  virtual ~ElasticEnergyCoefficient() { }
217 };
218 
219 void InitialDeformation(const Vector &x, Vector &y);
220 
221 void InitialVelocity(const Vector &x, Vector &v);
222 
223 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
224  GridFunction *field, const char *field_name = NULL,
225  bool init_vis = false);
226 
227 
228 int main(int argc, char *argv[])
229 {
230  // 1. Parse command-line options.
231  const char *mesh_file = "../../data/beam-quad.mesh";
232  int ref_levels = 2;
233  int order = 2;
234  int ode_solver_type = 3;
235  double t_final = 300.0;
236  double dt = 3.0;
237  double visc = 1e-2;
238  double mu = 0.25;
239  double K = 5.0;
240  bool visualization = true;
241  const char *nls = "newton";
242  int vis_steps = 1;
243 
244  // Relative and absolute tolerances for CVODE and ARKODE.
245  const double reltol = 1e-1, abstol = 1e-1;
246 
247  OptionsParser args(argc, argv);
248  args.AddOption(&mesh_file, "-m", "--mesh",
249  "Mesh file to use.");
250  args.AddOption(&ref_levels, "-r", "--refine",
251  "Number of times to refine the mesh uniformly.");
252  args.AddOption(&order, "-o", "--order",
253  "Order (degree) of the finite elements.");
254  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
255  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
256  " 4 - CVODE implicit, approximate Jacobian,\n\t"
257  " 5 - CVODE implicit, specified Jacobian,\n\t"
258  " 6 - ARKODE implicit, approximate Jacobian,\n\t"
259  " 7 - ARKODE implicit, specified Jacobian,\n\t"
260  " 11 - Forward Euler, 12 - RK2,\n\t"
261  " 13 - RK3 SSP, 14 - RK4,\n\t"
262  " 15 - CVODE (adaptive order) explicit,\n\t"
263  " 16 - ARKODE default (4th order) explicit.");
264  args.AddOption(&nls, "-nls", "--nonlinear-solver",
265  "Nonlinear systems solver: "
266  "\"newton\" (plain Newton) or \"kinsol\" (KINSOL).");
267  args.AddOption(&t_final, "-tf", "--t-final",
268  "Final time; start time is 0.");
269  args.AddOption(&dt, "-dt", "--time-step",
270  "Time step.");
271  args.AddOption(&visc, "-v", "--viscosity",
272  "Viscosity coefficient.");
273  args.AddOption(&mu, "-mu", "--shear-modulus",
274  "Shear modulus in the Neo-Hookean hyperelastic model.");
275  args.AddOption(&K, "-K", "--bulk-modulus",
276  "Bulk modulus in the Neo-Hookean hyperelastic model.");
277  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
278  "--no-visualization",
279  "Enable or disable GLVis visualization.");
280  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
281  "Visualize every n-th timestep.");
282  args.Parse();
283  if (!args.Good())
284  {
285  args.PrintUsage(cout);
286  return 1;
287  }
288  args.PrintOptions(cout);
289 
290  // 2. Read the mesh from the given mesh file. We can handle triangular,
291  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
292  Mesh *mesh = new Mesh(mesh_file, 1, 1);
293  int dim = mesh->Dimension();
294 
295  // 3. Define the ODE solver used for time integration. Several implicit
296  // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
297  // explicit Runge-Kutta methods are available.
298  ODESolver *ode_solver;
299  CVODESolver *cvode = NULL;
300  ARKODESolver *arkode = NULL;
301  SundialsJacSolver *sjsolver = NULL;
302  switch (ode_solver_type)
303  {
304  // Implicit L-stable methods
305  case 1: ode_solver = new BackwardEulerSolver; break;
306  case 2: ode_solver = new SDIRK23Solver(2); break;
307  case 3: ode_solver = new SDIRK33Solver; break;
308  case 4:
309  case 5:
310  cvode = new CVODESolver(CV_BDF, CV_NEWTON);
311  cvode->SetSStolerances(reltol, abstol);
312  cvode->SetMaxStep(dt);
313  if (ode_solver_type == 5)
314  {
315  sjsolver = new SundialsJacSolver;
316  cvode->SetLinearSolver(*sjsolver);
317  }
318  ode_solver = cvode; break;
319  case 6:
320  case 7:
321  arkode = new ARKODESolver(ARKODESolver::IMPLICIT);
322  arkode->SetSStolerances(reltol, abstol);
323  arkode->SetMaxStep(dt);
324  if (ode_solver_type == 7)
325  {
326  // Custom Jacobian inversion.
327  sjsolver = new SundialsJacSolver;
328  arkode->SetLinearSolver(*sjsolver);
329  }
330  ode_solver = arkode; break;
331  // Explicit methods
332  case 11: ode_solver = new ForwardEulerSolver; break;
333  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
334  case 13: ode_solver = new RK3SSPSolver; break;
335  case 14: ode_solver = new RK4Solver; break;
336  case 15:
337  cvode = new CVODESolver(CV_ADAMS, CV_FUNCTIONAL);
338  cvode->SetSStolerances(reltol, abstol);
339  cvode->SetMaxStep(dt);
340  ode_solver = cvode; break;
341  case 16:
342  arkode = new ARKODESolver(ARKODESolver::IMPLICIT);
343  arkode->SetSStolerances(reltol, abstol);
344  arkode->SetMaxStep(dt);
345  ode_solver = arkode; break;
346  // Implicit A-stable methods (not L-stable)
347  case 22: ode_solver = new ImplicitMidpointSolver; break;
348  case 23: ode_solver = new SDIRK23Solver; break;
349  case 24: ode_solver = new SDIRK34Solver; break;
350  default:
351  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
352  delete mesh;
353  return 3;
354  }
355 
356  map<string,HyperelasticOperator::NonlinearSolverType> nls_map;
357  nls_map["newton"] = HyperelasticOperator::NEWTON;
358  nls_map["kinsol"] = HyperelasticOperator::KINSOL;
359  if (nls_map.find(nls) == nls_map.end())
360  {
361  cout << "Unknown type of nonlinear solver: " << nls << endl;
362  return 4;
363  }
364 
365  // 4. Refine the mesh to increase the resolution. In this example we do
366  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
367  // command-line parameter.
368  for (int lev = 0; lev < ref_levels; lev++)
369  {
370  mesh->UniformRefinement();
371  }
372 
373  // 5. Define the vector finite element spaces representing the mesh
374  // deformation x, the velocity v, and the initial configuration, x_ref.
375  // Define also the elastic energy density, w, which is in a discontinuous
376  // higher-order space. Since x and v are integrated in time as a system,
377  // we group them together in block vector vx, with offsets given by the
378  // fe_offset array.
379  H1_FECollection fe_coll(order, dim);
380  FiniteElementSpace fespace(mesh, &fe_coll, dim);
381 
382  int fe_size = fespace.GetTrueVSize();
383  cout << "Number of velocity/deformation unknowns: " << fe_size << endl;
384  Array<int> fe_offset(3);
385  fe_offset[0] = 0;
386  fe_offset[1] = fe_size;
387  fe_offset[2] = 2*fe_size;
388 
389  BlockVector vx(fe_offset);
390  GridFunction v, x;
391  v.MakeTRef(&fespace, vx.GetBlock(0), 0);
392  x.MakeTRef(&fespace, vx.GetBlock(1), 0);
393 
394  GridFunction x_ref(&fespace);
395  mesh->GetNodes(x_ref);
396 
397  L2_FECollection w_fec(order + 1, dim);
398  FiniteElementSpace w_fespace(mesh, &w_fec);
399  GridFunction w(&w_fespace);
400 
401  // 6. Set the initial conditions for v and x, and the boundary conditions on
402  // a beam-like mesh (see description above).
404  v.ProjectCoefficient(velo);
405  v.SetTrueVector();
407  x.ProjectCoefficient(deform);
408  x.SetTrueVector();
409 
410  Array<int> ess_bdr(fespace.GetMesh()->bdr_attributes.Max());
411  ess_bdr = 0;
412  ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed
413 
414  // 7. Initialize the hyperelastic operator, the GLVis visualization and print
415  // the initial energies.
416  HyperelasticOperator oper(fespace, ess_bdr, visc, mu, K, nls_map[nls]);
417 
418  socketstream vis_v, vis_w;
419  if (visualization)
420  {
421  char vishost[] = "localhost";
422  int visport = 19916;
423  vis_v.open(vishost, visport);
424  vis_v.precision(8);
426  visualize(vis_v, mesh, &x, &v, "Velocity", true);
427  vis_w.open(vishost, visport);
428  if (vis_w)
429  {
430  oper.GetElasticEnergyDensity(x, w);
431  vis_w.precision(8);
432  visualize(vis_w, mesh, &x, &w, "Elastic energy density", true);
433  }
434  }
435 
436  double ee0 = oper.ElasticEnergy(x.GetTrueVector());
437  double ke0 = oper.KineticEnergy(v.GetTrueVector());
438  cout << "initial elastic energy (EE) = " << ee0 << endl;
439  cout << "initial kinetic energy (KE) = " << ke0 << endl;
440  cout << "initial total energy (TE) = " << (ee0 + ke0) << endl;
441 
442  double t = 0.0;
443  oper.SetTime(t);
444  ode_solver->Init(oper);
445 
446  // 8. Perform time-integration (looping over the time iterations, ti, with a
447  // time-step dt).
448  bool last_step = false;
449  for (int ti = 1; !last_step; ti++)
450  {
451  double dt_real = min(dt, t_final - t);
452 
453  ode_solver->Step(vx, t, dt_real);
454 
455  last_step = (t >= t_final - 1e-8*dt);
456 
457  if (last_step || (ti % vis_steps) == 0)
458  {
459  double ee = oper.ElasticEnergy(x.GetTrueVector());
460  double ke = oper.KineticEnergy(v.GetTrueVector());
461 
462  cout << "step " << ti << ", t = " << t << ", EE = " << ee << ", KE = "
463  << ke << ", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
464 
465  if (cvode) { cvode->PrintInfo(); }
466  else if (arkode) { arkode->PrintInfo(); }
467 
468  if (visualization)
469  {
471  visualize(vis_v, mesh, &x, &v);
472  if (vis_w)
473  {
474  oper.GetElasticEnergyDensity(x, w);
475  visualize(vis_w, mesh, &x, &w);
476  }
477  }
478  }
479  }
480 
481  // 9. Save the displaced mesh, the velocity and elastic energy.
482  {
484  GridFunction *nodes = &x;
485  int owns_nodes = 0;
486  mesh->SwapNodes(nodes, owns_nodes);
487  ofstream mesh_ofs("deformed.mesh");
488  mesh_ofs.precision(8);
489  mesh->Print(mesh_ofs);
490  mesh->SwapNodes(nodes, owns_nodes);
491  ofstream velo_ofs("velocity.sol");
492  velo_ofs.precision(8);
493  v.Save(velo_ofs);
494  ofstream ee_ofs("elastic_energy.sol");
495  ee_ofs.precision(8);
496  oper.GetElasticEnergyDensity(x, w);
497  w.Save(ee_ofs);
498  }
499 
500  // 10. Free the used memory.
501  delete ode_solver;
502  delete sjsolver;
503  delete mesh;
504 
505  return 0;
506 }
507 
508 
509 void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
510  GridFunction *field, const char *field_name, bool init_vis)
511 {
512  if (!out)
513  {
514  return;
515  }
516 
517  GridFunction *nodes = deformed_nodes;
518  int owns_nodes = 0;
519 
520  mesh->SwapNodes(nodes, owns_nodes);
521 
522  out << "solution\n" << *mesh << *field;
523 
524  mesh->SwapNodes(nodes, owns_nodes);
525 
526  if (init_vis)
527  {
528  out << "window_size 800 800\n";
529  out << "window_title '" << field_name << "'\n";
530  if (mesh->SpaceDimension() == 2)
531  {
532  out << "view 0 0\n"; // view from top
533  out << "keys jl\n"; // turn off perspective and light
534  }
535  out << "keys cm\n"; // show colorbar and mesh
536  out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
537  out << "pause\n";
538  }
539  out << flush;
540 }
541 
542 
543 ReducedSystemOperator::ReducedSystemOperator(
545  : Operator(M_->Height()), M(M_), S(S_), H(H_), Jacobian(NULL),
546  dt(0.0), v(NULL), x(NULL), w(height), z(height)
547 { }
548 
549 void ReducedSystemOperator::SetParameters(double dt_, const Vector *v_,
550  const Vector *x_)
551 {
552  dt = dt_; v = v_; x = x_;
553 }
554 
555 void ReducedSystemOperator::Mult(const Vector &k, Vector &y) const
556 {
557  // compute: y = H(x + dt*(v + dt*k)) + M*k + S*(v + dt*k)
558  add(*v, dt, k, w);
559  add(*x, dt, w, z);
560  H->Mult(z, y);
561  M->AddMult(k, y);
562  S->AddMult(w, y);
563 }
564 
565 Operator &ReducedSystemOperator::GetGradient(const Vector &k) const
566 {
567  delete Jacobian;
568  Jacobian = Add(1.0, M->SpMat(), dt, S->SpMat());
569  add(*v, dt, k, w);
570  add(*x, dt, w, z);
571  SparseMatrix *grad_H = dynamic_cast<SparseMatrix *>(&H->GetGradient(z));
572  Jacobian->Add(dt*dt, *grad_H);
573  return *Jacobian;
574 }
575 
576 ReducedSystemOperator::~ReducedSystemOperator()
577 {
578  delete Jacobian;
579 }
580 
581 
582 int SundialsJacSolver::InitSystem(void *sundials_mem)
583 {
584  TimeDependentOperator *td_oper = GetTimeDependentOperator(sundials_mem);
585  HyperelasticOperator *he_oper;
586 
587  // During development, we use dynamic_cast<> to ensure the setup is correct:
588  he_oper = dynamic_cast<HyperelasticOperator*>(td_oper);
589  MFEM_VERIFY(he_oper, "operator is not HyperelasticOperator");
590 
591  // When the implementation is finalized, we can switch to static_cast<>:
592  // he_oper = static_cast<HyperelasticOperator*>(td_oper);
593 
594  he_oper->InitSundialsJacSolver(*this);
595  return 0;
596 }
597 
598 int SundialsJacSolver::SetupSystem(void *sundials_mem, int conv_fail,
599  const Vector &y_pred, const Vector &f_pred,
600  int &jac_cur, Vector &v_temp1,
601  Vector &v_temp2, Vector &v_temp3)
602 {
603  int sc = y_pred.Size() / 2;
604  const Vector x(y_pred.GetData() + sc, sc);
605  double dt = GetTimeStep(sundials_mem);
606 
607  // J = M + dt*(S + dt*grad(H))
608  delete Jacobian;
609  Jacobian = Add(1.0, M->SpMat(), dt, S->SpMat());
610  grad_H = dynamic_cast<SparseMatrix *>(&H->GetGradient(x));
611  Jacobian->Add(dt * dt, *grad_H);
612 
613  J_solver->SetOperator(*Jacobian);
614 
615  jac_cur = 1;
616  return 0;
617 }
618 
619 int SundialsJacSolver::SolveSystem(void *sundials_mem, Vector &b,
620  const Vector &weight, const Vector &y_cur,
621  const Vector &f_cur)
622 {
623  int sc = b.Size() / 2;
624  // Vector x(y_cur.GetData() + sc, sc);
625  Vector b_v(b.GetData() + 0, sc);
626  Vector b_x(b.GetData() + sc, sc);
627  Vector rhs(sc);
628  double dt = GetTimeStep(sundials_mem);
629 
630  // rhs = M b_v - dt*grad(H) b_x
631  grad_H->Mult(b_x, rhs);
632  rhs *= -dt;
633  M->AddMult(b_v, rhs);
634 
635  J_solver->iterative_mode = false;
636  J_solver->Mult(rhs, b_v);
637 
638  b_x.Add(dt, b_v);
639 
640  return 0;
641 }
642 
643 int SundialsJacSolver::FreeSystem(void *sundials_mem)
644 {
645  delete Jacobian;
646  return 0;
647 }
648 
649 
650 HyperelasticOperator::HyperelasticOperator(FiniteElementSpace &f,
651  Array<int> &ess_bdr, double visc,
652  double mu, double K,
653  NonlinearSolverType nls_type)
654  : TimeDependentOperator(2*f.GetTrueVSize(), 0.0), fespace(f),
655  M(&fespace), S(&fespace), H(&fespace),
656  viscosity(visc), z(height/2)
657 {
658  const double rel_tol = 1e-8;
659  const int skip_zero_entries = 0;
660 
661  const double ref_density = 1.0; // density in the reference configuration
662  ConstantCoefficient rho0(ref_density);
663  M.AddDomainIntegrator(new VectorMassIntegrator(rho0));
664  M.Assemble(skip_zero_entries);
665  Array<int> ess_tdof_list;
666  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
667  SparseMatrix tmp;
668  M.FormSystemMatrix(ess_tdof_list, tmp);
669 
670  M_solver.iterative_mode = false;
671  M_solver.SetRelTol(rel_tol);
672  M_solver.SetAbsTol(0.0);
673  M_solver.SetMaxIter(30);
674  M_solver.SetPrintLevel(0);
675  M_solver.SetPreconditioner(M_prec);
676  M_solver.SetOperator(M.SpMat());
677 
678  model = new NeoHookeanModel(mu, K);
679  H.AddDomainIntegrator(new HyperelasticNLFIntegrator(model));
680  H.SetEssentialTrueDofs(ess_tdof_list);
681 
682  ConstantCoefficient visc_coeff(viscosity);
683  S.AddDomainIntegrator(new VectorDiffusionIntegrator(visc_coeff));
684  S.Assemble(skip_zero_entries);
685  S.FormSystemMatrix(ess_tdof_list, tmp);
686 
687  reduced_oper = new ReducedSystemOperator(&M, &S, &H);
688 
689 #ifndef MFEM_USE_SUITESPARSE
690  J_prec = new DSmoother(1);
691  MINRESSolver *J_minres = new MINRESSolver;
692  J_minres->SetRelTol(rel_tol);
693  J_minres->SetAbsTol(0.0);
694  J_minres->SetMaxIter(300);
695  J_minres->SetPrintLevel(-1);
696  J_minres->SetPreconditioner(*J_prec);
697  J_solver = J_minres;
698 #else
699  J_solver = new UMFPackSolver;
700  J_prec = NULL;
701 #endif
702 
703  if (nls_type == KINSOL)
704  {
705  KinSolver *kinsolver = new KinSolver(KIN_NONE, true);
706  kinsolver->SetMaxSetupCalls(4);
707  newton_solver = kinsolver;
708  newton_solver->SetMaxIter(200);
709  newton_solver->SetRelTol(rel_tol);
710  newton_solver->SetPrintLevel(0);
711  }
712  else
713  {
714  newton_solver = new NewtonSolver();
715  newton_solver->SetMaxIter(10);
716  newton_solver->SetRelTol(rel_tol);
717  newton_solver->SetPrintLevel(-1);
718  }
719  newton_solver->SetSolver(*J_solver);
720  newton_solver->iterative_mode = false;
721  newton_solver->SetOperator(*reduced_oper);
722 }
723 
724 void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
725 {
726  // Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt
727  int sc = height/2;
728  Vector v(vx.GetData() + 0, sc);
729  Vector x(vx.GetData() + sc, sc);
730  Vector dv_dt(dvx_dt.GetData() + 0, sc);
731  Vector dx_dt(dvx_dt.GetData() + sc, sc);
732 
733  H.Mult(x, z);
734  if (viscosity != 0.0)
735  {
736  S.AddMult(v, z);
737  }
738  z.Neg(); // z = -z
739  M_solver.Mult(z, dv_dt);
740 
741  dx_dt = v;
742 }
743 
744 void HyperelasticOperator::ImplicitSolve(const double dt,
745  const Vector &vx, Vector &dvx_dt)
746 {
747  int sc = height/2;
748  Vector v(vx.GetData() + 0, sc);
749  Vector x(vx.GetData() + sc, sc);
750  Vector dv_dt(dvx_dt.GetData() + 0, sc);
751  Vector dx_dt(dvx_dt.GetData() + sc, sc);
752 
753  // By eliminating kx from the coupled system:
754  // kv = -M^{-1}*[H(x + dt*kx) + S*(v + dt*kv)]
755  // kx = v + dt*kv
756  // we reduce it to a nonlinear equation for kv, represented by the
757  // reduced_oper. This equation is solved with the newton_solver
758  // object (using J_solver and J_prec internally).
759  reduced_oper->SetParameters(dt, &v, &x);
760  Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
761  newton_solver->Mult(zero, dv_dt);
762  MFEM_VERIFY(newton_solver->GetConverged(),
763  "Nonlinear solver did not converge.");
764 #ifdef MFEM_DEBUG
765  cout << " num nonlin sol iters = " << newton_solver->GetNumIterations()
766  << ", final norm = " << newton_solver->GetFinalNorm() << '\n';
767 #endif
768  add(v, dt, dv_dt, dx_dt);
769 }
770 
771 void HyperelasticOperator::InitSundialsJacSolver(SundialsJacSolver &sjsolv)
772 {
773  sjsolv.SetOperators(M, S, H, *J_solver);
774 }
775 
776 double HyperelasticOperator::ElasticEnergy(const Vector &x) const
777 {
778  return H.GetEnergy(x);
779 }
780 
781 double HyperelasticOperator::KineticEnergy(const Vector &v) const
782 {
783  return 0.5*M.InnerProduct(v, v);
784 }
785 
786 void HyperelasticOperator::GetElasticEnergyDensity(
787  const GridFunction &x, GridFunction &w) const
788 {
789  ElasticEnergyCoefficient w_coeff(*model, x);
790  w.ProjectCoefficient(w_coeff);
791 }
792 
793 HyperelasticOperator::~HyperelasticOperator()
794 {
795  delete newton_solver;
796  delete J_solver;
797  delete J_prec;
798  delete reduced_oper;
799  delete model;
800 }
801 
802 
804  const IntegrationPoint &ip)
805 {
806  model.SetTransformation(T);
807  x.GetVectorGradient(T, J);
808  // return model.EvalW(J); // in reference configuration
809  return model.EvalW(J)/J.Det(); // in deformed configuration
810 }
811 
812 
813 void InitialDeformation(const Vector &x, Vector &y)
814 {
815  // set the initial configuration to be the same as the reference, stress
816  // free, configuration
817  y = x;
818 }
819 
820 void InitialVelocity(const Vector &x, Vector &v)
821 {
822  const int dim = x.Size();
823  const double s = 0.1/64.;
824 
825  v = 0.0;
826  v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
827  v(0) = -s*x(0)*x(0);
828 }
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:376
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
void InitialDeformation(const Vector &x, Vector &y)
Definition: ex10.cpp:588
Conjugate gradient method.
Definition: solvers.hpp:111
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
Data type for scaled Jacobi-type smoother of sparse matrix.
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
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
STL namespace.
MINRES method.
Definition: solvers.hpp:221
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 InitialVelocity(const Vector &x, Vector &v)
Definition: ex10.cpp:595
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:344
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
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
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
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1034
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:256
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:66
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
int main(int argc, char *argv[])
Definition: ex10.cpp:156
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
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:105
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:225
Class for integration point with weight.
Definition: intrules.hpp:25
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1377
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
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
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2349
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
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:77
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:128
void Neg()
(*this) = -(*this)
Definition: vector.cpp:252