MFEM  v3.4 Finite element discretization library
ex16.cpp
Go to the documentation of this file.
1 // MFEM Example 16
2 // SUNDIALS Modification
3 //
4 // Compile with: make ex16
5 //
6 // Sample runs: ex16
7 // ex16 -m ../../data/inline-tri.mesh
8 // ex16 -m ../../data/disc-nurbs.mesh -tf 2
9 // ex16 -s 12 -a 0.0 -k 1.0
10 // ex16 -s 1 -a 1.0 -k 0.0 -dt 1e-4 -tf 5e-2 -vs 25
11 // ex16 -s 2 -a 0.5 -k 0.5 -o 4 -dt 1e-4 -tf 2e-2 -vs 25
12 // ex16 -s 3 -dt 1.0e-4 -tf 4.0e-2 -vs 40
13 // ex16 -m ../../data/fichera-q2.mesh
14 // ex16 -m ../../data/escher.mesh
15 // ex16 -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
16 // ex16 -m ../../data/amr-quad.mesh -o 4 -r 0
17 // ex16 -m ../../data/amr-hex.mesh -o 2 -r 0
18 //
19 // Description: This example solves a time dependent nonlinear heat equation
20 // problem of the form du/dt = C(u), with a non-linear diffusion
21 // operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
22 //
23 // The example demonstrates the use of nonlinear operators (the
24 // class ConductionOperator defining C(u)), as well as their
25 // implicit time integration. Note that implementing the method
26 // ConductionOperator::ImplicitSolve is the only requirement for
27 // high-order implicit (SDIRK) time integration. By default, this
28 // example uses the SUNDIALS ODE solvers from CVODE and ARKODE.
29 //
30 // We recommend viewing examples 2, 9 and 10 before viewing this
31 // example.
32
33 #include "mfem.hpp"
34 #include <fstream>
35 #include <iostream>
36
37 using namespace std;
38 using namespace mfem;
39
40 /** After spatial discretization, the conduction model can be written as:
41  *
42  * du/dt = M^{-1}(-Ku)
43  *
44  * where u is the vector representing the temperature, M is the mass matrix,
45  * and K is the diffusion operator with diffusivity depending on u:
46  * (\kappa + \alpha u).
47  *
48  * Class ConductionOperator represents the right-hand side of the above ODE.
49  */
50 class ConductionOperator : public TimeDependentOperator
51 {
52 protected:
53  FiniteElementSpace &fespace;
54  Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
55
56  BilinearForm *M;
57  BilinearForm *K;
58
59  SparseMatrix Mmat, Kmat;
60  SparseMatrix *T; // T = M + dt K
61  double current_dt;
62
63  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
64  DSmoother M_prec; // Preconditioner for the mass matrix M
65
66  CGSolver T_solver; // Implicit solver for T = M + dt K
67  DSmoother T_prec; // Preconditioner for the implicit solver
68
69  double alpha, kappa;
70
71  mutable Vector z; // auxiliary vector
72
73 public:
74  ConductionOperator(FiniteElementSpace &f, double alpha, double kappa,
75  const Vector &u);
76
77  virtual void Mult(const Vector &u, Vector &du_dt) const;
78  /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
79  This is the only requirement for high-order SDIRK implicit integration.*/
80  virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
81
82  /** Solve the system (M + dt K) y = M b. The result y replaces the input b.
83  This method is used by the implicit SUNDIALS solvers. */
84  void SundialsSolve(const double dt, Vector &b);
85
86  /// Update the diffusion BilinearForm K using the given true-dof vector u.
87  void SetParameters(const Vector &u);
88
89  virtual ~ConductionOperator();
90 };
91
92 /// Custom Jacobian system solver for the SUNDIALS time integrators.
93 /** For the ODE system represented by ConductionOperator
94
95  M du/dt = -K(u),
96
97  this class facilitates the solution of linear systems of the form
98
99  (M + γK) y = M b,
100
101  for given b, u (not used), and γ = GetTimeStep(). */
102 class SundialsJacSolver : public SundialsODELinearSolver
103 {
104 private:
105  ConductionOperator *oper;
106
107 public:
108  SundialsJacSolver() : oper(NULL) { }
109
110  int InitSystem(void *sundials_mem);
111  int SetupSystem(void *sundials_mem, int conv_fail,
112  const Vector &y_pred, const Vector &f_pred, int &jac_cur,
113  Vector &v_temp1, Vector &v_temp2, Vector &v_temp3);
114  int SolveSystem(void *sundials_mem, Vector &b, const Vector &weight,
115  const Vector &y_cur, const Vector &f_cur);
116  int FreeSystem(void *sundials_mem);
117 };
118
119 double InitialTemperature(const Vector &x);
120
121 int main(int argc, char *argv[])
122 {
123  // 1. Parse command-line options.
124  const char *mesh_file = "../../data/star.mesh";
125  int ref_levels = 2;
126  int order = 2;
127  int ode_solver_type = 11; // 11 = CVODE implicit
128  double t_final = 0.5;
129  double dt = 1.0e-2;
130  double alpha = 1.0e-2;
131  double kappa = 0.5;
132  bool visualization = true;
133  bool visit = false;
134  int vis_steps = 5;
135
136  // Relative and absolute tolerances for CVODE and ARKODE.
137  const double reltol = 1e-4, abstol = 1e-4;
138
139  int precision = 8;
140  cout.precision(precision);
141
142  OptionsParser args(argc, argv);
144  "Mesh file to use.");
146  "Number of times to refine the mesh uniformly.");
148  "Order (degree) of the finite elements.");
150  "ODE solver:\n"
151  "\t 1/11 - CVODE (explicit/implicit),\n"
152  "\t 2/12 - ARKODE (default explicit/implicit),\n"
153  "\t 3 - ARKODE (Fehlberg-6-4-5)\n"
154  "\t 4 - Forward Euler, 5 - RK2, 6 - RK3 SSP, 7 - RK4,\n"
155  "\t 8 - Backward Euler, 9 - SDIRK23, 10 - SDIRK33.");
157  "Final time; start time is 0.");
159  "Time step.");
161  "Alpha coefficient.");
163  "Kappa coefficient offset.");
165  "--no-visualization",
166  "Enable or disable GLVis visualization.");
168  "--no-visit-datafiles",
169  "Save data files for VisIt (visit.llnl.gov) visualization.");
171  "Visualize every n-th timestep.");
172  args.Parse();
173  if (!args.Good())
174  {
175  args.PrintUsage(cout);
176  return 1;
177  }
178  args.PrintOptions(cout);
179
180  // 2. Read the mesh from the given mesh file. We can handle triangular,
181  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
182  Mesh *mesh = new Mesh(mesh_file, 1, 1);
183  int dim = mesh->Dimension();
184
185  // 3. Define the ODE solver used for time integration. Several
186  // SUNDIALS solvers are available, as well as included both
187  // explicit and implicit MFEM ODE solvers.
188  ODESolver *ode_solver = NULL;
189  CVODESolver *cvode = NULL;
190  ARKODESolver *arkode = NULL;
191  SundialsJacSolver sun_solver; // Used by the implicit SUNDIALS ode solvers.
192  switch (ode_solver_type)
193  {
194  // SUNDIALS solvers
195  case 1:
196  cvode = new CVODESolver(CV_ADAMS, CV_FUNCTIONAL);
197  cvode->SetSStolerances(reltol, abstol);
198  cvode->SetMaxStep(dt);
199  ode_solver = cvode; break;
200  case 11:
201  cvode = new CVODESolver(CV_BDF, CV_NEWTON);
202  cvode->SetLinearSolver(sun_solver);
203  cvode->SetSStolerances(reltol, abstol);
204  cvode->SetMaxStep(dt);
205  ode_solver = cvode; break;
206  case 2:
207  case 3:
208  arkode = new ARKODESolver(ARKODESolver::EXPLICIT);
209  arkode->SetSStolerances(reltol, abstol);
210  arkode->SetMaxStep(dt);
211  if (ode_solver_type == 3) { arkode->SetERKTableNum(FEHLBERG_13_7_8); }
212  ode_solver = arkode; break;
213  case 12:
214  arkode = new ARKODESolver(ARKODESolver::IMPLICIT);
215  arkode->SetLinearSolver(sun_solver);
216  arkode->SetSStolerances(reltol, abstol);
217  arkode->SetMaxStep(dt);
218  ode_solver = arkode; break;
219  // Other MFEM explicit methods
220  case 4: ode_solver = new ForwardEulerSolver; break;
221  case 5: ode_solver = new RK2Solver(0.5); break; // midpoint method
222  case 6: ode_solver = new RK3SSPSolver; break;
223  case 7: ode_solver = new RK4Solver; break;
224  // MFEM implicit L-stable methods
225  case 8: ode_solver = new BackwardEulerSolver; break;
226  case 9: ode_solver = new SDIRK23Solver(2); break;
227  case 10: ode_solver = new SDIRK33Solver; break;
228  default:
229  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
230  delete mesh;
231  return 3;
232  }
233
234  // Since we want to update the diffusion coefficient after every time step,
235  // we need to use the "one-step" mode of the SUNDIALS solvers.
236  if (cvode) { cvode->SetStepMode(CV_ONE_STEP); }
237  if (arkode) { arkode->SetStepMode(ARK_ONE_STEP); }
238
239  // 4. Refine the mesh to increase the resolution. In this example we do
240  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
241  // command-line parameter.
242  for (int lev = 0; lev < ref_levels; lev++)
243  {
244  mesh->UniformRefinement();
245  }
246
247  // 5. Define the vector finite element space representing the current and the
248  // initial temperature, u_ref.
249  H1_FECollection fe_coll(order, dim);
250  FiniteElementSpace fespace(mesh, &fe_coll);
251
252  int fe_size = fespace.GetTrueVSize();
253  cout << "Number of temperature unknowns: " << fe_size << endl;
254
255  GridFunction u_gf(&fespace);
256
257  // 6. Set the initial conditions for u. All boundaries are considered
258  // natural.
260  u_gf.ProjectCoefficient(u_0);
261  Vector u;
262  u_gf.GetTrueDofs(u);
263
264  // 7. Initialize the conduction operator and the visualization.
265  ConductionOperator oper(fespace, alpha, kappa, u);
266
267  u_gf.SetFromTrueDofs(u);
268  {
269  ofstream omesh("ex16.mesh");
270  omesh.precision(precision);
271  mesh->Print(omesh);
272  ofstream osol("ex16-init.gf");
273  osol.precision(precision);
274  u_gf.Save(osol);
275  }
276
277  VisItDataCollection visit_dc("Example16", mesh);
278  visit_dc.RegisterField("temperature", &u_gf);
279  if (visit)
280  {
281  visit_dc.SetCycle(0);
282  visit_dc.SetTime(0.0);
283  visit_dc.Save();
284  }
285
286  socketstream sout;
287  if (visualization)
288  {
289  char vishost[] = "localhost";
290  int visport = 19916;
291  sout.open(vishost, visport);
292  if (!sout)
293  {
294  cout << "Unable to connect to GLVis server at "
295  << vishost << ':' << visport << endl;
296  visualization = false;
297  cout << "GLVis visualization disabled.\n";
298  }
299  else
300  {
301  sout.precision(precision);
302  sout << "solution\n" << *mesh << u_gf;
303  sout << "pause\n";
304  sout << flush;
305  cout << "GLVis visualization paused."
306  << " Press space (in the GLVis window) to resume it.\n";
307  }
308  }
309
310  // 8. Perform time-integration (looping over the time iterations, ti, with a
311  // time-step dt).
312  cout << "Integrating the ODE ..." << endl;
313  tic_toc.Clear();
314  tic_toc.Start();
315  ode_solver->Init(oper);
316  double t = 0.0;
317
318  bool last_step = false;
319  for (int ti = 1; !last_step; ti++)
320  {
321  double dt_real = min(dt, t_final - t);
322
323  // Note that since we are using the "one-step" mode of the SUNDIALS
324  // solvers, they will, generally, step over the final time and will not
325  // explicitly perform the interpolation to t_final as they do in the
326  // "normal" step mode.
327
328  ode_solver->Step(u, t, dt_real);
329
330  last_step = (t >= t_final - 1e-8*dt);
331
332  if (last_step || (ti % vis_steps) == 0)
333  {
334  cout << "step " << ti << ", t = " << t << endl;
335  if (cvode) { cvode->PrintInfo(); }
336  if (arkode) { arkode->PrintInfo(); }
337
338  u_gf.SetFromTrueDofs(u);
339  if (visualization)
340  {
341  sout << "solution\n" << *mesh << u_gf << flush;
342  }
343
344  if (visit)
345  {
346  visit_dc.SetCycle(ti);
347  visit_dc.SetTime(t);
348  visit_dc.Save();
349  }
350  }
351  oper.SetParameters(u);
352  }
353  tic_toc.Stop();
354  cout << "Done, " << tic_toc.RealTime() << "s." << endl;
355
356  // 9. Save the final solution. This output can be viewed later using GLVis:
357  // "glvis -m ex16.mesh -g ex16-final.gf".
358  {
359  ofstream osol("ex16-final.gf");
360  osol.precision(precision);
361  u_gf.Save(osol);
362  }
363
364  // 10. Free the used memory.
365  delete ode_solver;
366  delete mesh;
367
368  return 0;
369 }
370
371 ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
372  double kap, const Vector &u)
373  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
374  T(NULL), current_dt(0.0), z(height)
375 {
376  const double rel_tol = 1e-8;
377
378  M = new BilinearForm(&fespace);
380  M->Assemble();
381  M->FormSystemMatrix(ess_tdof_list, Mmat);
382
383  M_solver.iterative_mode = false;
384  M_solver.SetRelTol(rel_tol);
385  M_solver.SetAbsTol(0.0);
386  M_solver.SetMaxIter(50);
387  M_solver.SetPrintLevel(0);
388  M_solver.SetPreconditioner(M_prec);
389  M_solver.SetOperator(Mmat);
390
391  alpha = al;
392  kappa = kap;
393
394  T_solver.iterative_mode = false;
395  T_solver.SetRelTol(rel_tol);
396  T_solver.SetAbsTol(0.0);
397  T_solver.SetMaxIter(100);
398  T_solver.SetPrintLevel(0);
399  T_solver.SetPreconditioner(T_prec);
400
401  SetParameters(u);
402 }
403
404 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
405 {
406  // Compute:
407  // du_dt = M^{-1}*-K(u)
408  // for du_dt
409  Kmat.Mult(u, z);
410  z.Neg(); // z = -z
411  M_solver.Mult(z, du_dt);
412 }
413
414 void ConductionOperator::ImplicitSolve(const double dt,
415  const Vector &u, Vector &du_dt)
416 {
417  // Solve the equation:
418  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
419  // for du_dt
420  if (!T)
421  {
422  T = Add(1.0, Mmat, dt, Kmat);
423  current_dt = dt;
424  T_solver.SetOperator(*T);
425  }
426  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
427  Kmat.Mult(u, z);
428  z.Neg();
429  T_solver.Mult(z, du_dt);
430 }
431
432 void ConductionOperator::SundialsSolve(const double dt, Vector &b)
433 {
434  // Solve the system (M + dt K) y = M b. The result y replaces the input b.
435  if (!T || dt != current_dt)
436  {
437  delete T;
438  T = Add(1.0, Mmat, dt, Kmat);
439  current_dt = dt;
440  T_solver.SetOperator(*T);
441  }
442  Mmat.Mult(b, z);
443  T_solver.Mult(z, b);
444 }
445
446 void ConductionOperator::SetParameters(const Vector &u)
447 {
448  GridFunction u_alpha_gf(&fespace);
449  u_alpha_gf.SetFromTrueDofs(u);
450  for (int i = 0; i < u_alpha_gf.Size(); i++)
451  {
452  u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
453  }
454
455  delete K;
456  K = new BilinearForm(&fespace);
457
458  GridFunctionCoefficient u_coeff(&u_alpha_gf);
459
461  K->Assemble();
462  K->FormSystemMatrix(ess_tdof_list, Kmat);
463  delete T;
464  T = NULL; // re-compute T on the next ImplicitSolve or SundialsSolve
465 }
466
467 ConductionOperator::~ConductionOperator()
468 {
469  delete T;
470  delete M;
471  delete K;
472 }
473
474
475 int SundialsJacSolver::InitSystem(void *sundials_mem)
476 {
477  TimeDependentOperator *td_oper = GetTimeDependentOperator(sundials_mem);
478
479  // During development, we use dynamic_cast<> to ensure the setup is correct:
480  oper = dynamic_cast<ConductionOperator*>(td_oper);
481  MFEM_VERIFY(oper, "operator is not ConductionOperator");
482
483  // When the implementation is finalized, we can switch to static_cast<>:
484  // oper = static_cast<ConductionOperator*>(td_oper);
485
486  return 0;
487 }
488
489 int SundialsJacSolver::SetupSystem(void *sundials_mem, int conv_fail,
490  const Vector &y_pred, const Vector &f_pred,
491  int &jac_cur, Vector &v_temp1,
492  Vector &v_temp2, Vector &v_temp3)
493 {
494  jac_cur = 1;
495
496  return 0;
497 }
498
499 int SundialsJacSolver::SolveSystem(void *sundials_mem, Vector &b,
500  const Vector &weight, const Vector &y_cur,
501  const Vector &f_cur)
502 {
503  oper->SundialsSolve(GetTimeStep(sundials_mem), b);
504
505  return 0;
506 }
507
508 int SundialsJacSolver::FreeSystem(void *sundials_mem)
509 {
510  return 0;
511 }
512
513
514 double InitialTemperature(const Vector &x)
515 {
516  if (x.Norml2() < 0.5)
517  {
518  return 2.0;
519  }
520  else
521  {
522  return 1.0;
523  }
524 }
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1034
Definition: solvers.hpp:111
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
double InitialTemperature(const Vector &x)
Definition: ex16.cpp:381
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
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
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:151
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:666
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]...
ARKode supports two modes, specified by itask: ARK_NORMAL (default) and ARK_ONE_STEP.
Definition: sundials.cpp:509
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
virtual void Save()
Save the collection and a VisIt root file.
StopWatch tic_toc
Definition: tic_toc.cpp:447
double RealTime()
Definition: tic_toc.cpp:426
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Definition: sundials.cpp:233
int main(int argc, char *argv[])
Definition: ex1.cpp:45
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:212
void SetSStolerances(double reltol, double abstol)
Specify the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:477
void SetERKTableNum(int table_num)
Choose a specific Butcher table for explicit RK method.
Definition: sundials.cpp:528
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2349
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
int dim
Definition: ex3.cpp:47
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
Data collection with VisIt I/O routines.
void SetMaxStep(double dt_max)
Set the maximum time step of the linear multistep method.
Definition: sundials.hpp:183
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction. If all dofs are true, then tv will be set to point to th...
Definition: gridfunc.cpp:312
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:256
int Dimension() const
Definition: mesh.hpp:645
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
Wrapper for SUNDIALS&#39; 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:486
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:66
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&#39; CVODE and ARKODE solve...
Definition: sundials.hpp:47
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1377
int open(const char hostname[], int port)
const double alpha
Definition: ex15.cpp:337
class for C-function coefficient
double kappa
Definition: ex3.cpp:46
Vector data type.
Definition: vector.hpp:48
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:79
void PrintInfo() const
Print CVODE statistics.
Definition: sundials.cpp:400
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
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:327
void PrintInfo() const
Print ARKODE statistics.
Definition: sundials.cpp:691