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