MFEM  v4.6.0 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 8 -a 1.0 -k 0.0 -dt 1e-4 -tf 5e-2 -vs 25
11 // ex16 -s 9 -a 0.5 -k 0.5 -o 4 -dt 1e-4 -tf 2e-2 -vs 25
12 // ex16 -s 10 -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
62  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
63  DSmoother M_prec; // Preconditioner for the mass matrix M
64
65  CGSolver T_solver; // Implicit solver for T = M + dt K
66  DSmoother T_prec; // Preconditioner for the implicit solver
67
68  double alpha, kappa;
69
70  mutable Vector z; // auxiliary vector
71
72 public:
73  ConductionOperator(FiniteElementSpace &f, double alpha, double kappa,
74  const Vector &u);
75
76  virtual void Mult(const Vector &u, Vector &du_dt) const;
77
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  /// Custom Jacobian system solver for the SUNDIALS time integrators.
83  /** For the ODE system represented by ConductionOperator
84
85  M du/dt = -K(u),
86
87  this class facilitates the solution of linear systems of the form
88
89  (M + γK) y = M b,
90
91  for given b, u (not used), and γ = GetTimeStep(). */
92
93  /** Setup the system (M + dt K) x = M b. This method is used by the implicit
94  SUNDIALS solvers. */
95  virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
96  int jok, int *jcur, double gamma);
97
98  /** Solve the system (M + dt K) x = M b. This method is used by the implicit
99  SUNDIALS solvers. */
100  virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
101
102  /// Update the diffusion BilinearForm K using the given true-dof vector u.
103  void SetParameters(const Vector &u);
104
105  virtual ~ConductionOperator();
106 };
107
108 double InitialTemperature(const Vector &x);
109
110 int main(int argc, char *argv[])
111 {
112  // 0. Initialize SUNDIALS.
113  Sundials::Init();
114
115  // 1. Parse command-line options.
116  const char *mesh_file = "../../data/star.mesh";
117  int ref_levels = 2;
118  int order = 2;
119  int ode_solver_type = 9; // CVODE implicit BDF
120  double t_final = 0.5;
121  double dt = 1.0e-2;
122  double alpha = 1.0e-2;
123  double kappa = 0.5;
124  bool visualization = true;
125  bool visit = false;
126  int vis_steps = 5;
127
128  // Relative and absolute tolerances for CVODE and ARKODE.
129  const double reltol = 1e-4, abstol = 1e-4;
130
131  int precision = 8;
132  cout.precision(precision);
133
134  OptionsParser args(argc, argv);
135  args.AddOption(&mesh_file, "-m", "--mesh",
136  "Mesh file to use.");
137  args.AddOption(&ref_levels, "-r", "--refine",
138  "Number of times to refine the mesh uniformly.");
139  args.AddOption(&order, "-o", "--order",
140  "Order (degree) of the finite elements.");
141  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
142  "ODE solver:\n\t"
143  "1 - Forward Euler,\n\t"
144  "2 - RK2,\n\t"
145  "3 - RK3 SSP,\n\t"
146  "4 - RK4,\n\t"
147  "5 - Backward Euler,\n\t"
148  "6 - SDIRK 2,\n\t"
149  "7 - SDIRK 3,\n\t"
150  "8 - CVODE (implicit Adams),\n\t"
151  "9 - CVODE (implicit BDF),\n\t"
152  "10 - ARKODE (default explicit),\n\t"
153  "11 - ARKODE (explicit Fehlberg-6-4-5),\n\t"
154  "12 - ARKODE (default impicit).");
155  args.AddOption(&t_final, "-tf", "--t-final",
156  "Final time; start time is 0.");
157  args.AddOption(&dt, "-dt", "--time-step",
158  "Time step.");
159  args.AddOption(&alpha, "-a", "--alpha",
160  "Alpha coefficient.");
161  args.AddOption(&kappa, "-k", "--kappa",
162  "Kappa coefficient offset.");
163  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
164  "--no-visualization",
165  "Enable or disable GLVis visualization.");
166  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
167  "--no-visit-datafiles",
168  "Save data files for VisIt (visit.llnl.gov) visualization.");
169  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
170  "Visualize every n-th timestep.");
171  args.Parse();
172  if (!args.Good())
173  {
174  args.PrintUsage(cout);
175  return 1;
176  }
177  if (ode_solver_type < 1 || ode_solver_type > 12)
178  {
179  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
180  return 3;
181  }
182  args.PrintOptions(cout);
183
184  // 2. Read the mesh from the given mesh file. We can handle triangular,
185  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
186  Mesh *mesh = new Mesh(mesh_file, 1, 1);
187  int dim = mesh->Dimension();
188
189  // 3. Refine the mesh to increase the resolution. In this example we do
190  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
191  // command-line parameter.
192  for (int lev = 0; lev < ref_levels; lev++)
193  {
194  mesh->UniformRefinement();
195  }
196
197  // 4. Define the vector finite element space representing the current and the
198  // initial temperature, u_ref.
199  H1_FECollection fe_coll(order, dim);
200  FiniteElementSpace fespace(mesh, &fe_coll);
201
202  int fe_size = fespace.GetTrueVSize();
203  cout << "Number of temperature unknowns: " << fe_size << endl;
204
205  GridFunction u_gf(&fespace);
206
207  // 5. Set the initial conditions for u. All boundaries are considered
208  // natural.
210  u_gf.ProjectCoefficient(u_0);
211  Vector u;
212  u_gf.GetTrueDofs(u);
213
214  // 6. Initialize the conduction operator and the visualization.
215  ConductionOperator oper(fespace, alpha, kappa, u);
216
217  u_gf.SetFromTrueDofs(u);
218  {
219  ofstream omesh("ex16.mesh");
220  omesh.precision(precision);
221  mesh->Print(omesh);
222  ofstream osol("ex16-init.gf");
223  osol.precision(precision);
224  u_gf.Save(osol);
225  }
226
227  VisItDataCollection visit_dc("Example16", mesh);
228  visit_dc.RegisterField("temperature", &u_gf);
229  if (visit)
230  {
231  visit_dc.SetCycle(0);
232  visit_dc.SetTime(0.0);
233  visit_dc.Save();
234  }
235
236  socketstream sout;
237  if (visualization)
238  {
239  char vishost[] = "localhost";
240  int visport = 19916;
241  sout.open(vishost, visport);
242  if (!sout)
243  {
244  cout << "Unable to connect to GLVis server at "
245  << vishost << ':' << visport << endl;
246  visualization = false;
247  cout << "GLVis visualization disabled.\n";
248  }
249  else
250  {
251  sout.precision(precision);
252  sout << "solution\n" << *mesh << u_gf;
253  sout << "pause\n";
254  sout << flush;
255  cout << "GLVis visualization paused."
256  << " Press space (in the GLVis window) to resume it.\n";
257  }
258  }
259
260  // 7. Define the ODE solver used for time integration.
261  double t = 0.0;
262  ODESolver *ode_solver = NULL;
263  CVODESolver *cvode = NULL;
264  ARKStepSolver *arkode = NULL;
265  switch (ode_solver_type)
266  {
267  // MFEM explicit methods
268  case 1: ode_solver = new ForwardEulerSolver; break;
269  case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
270  case 3: ode_solver = new RK3SSPSolver; break;
271  case 4: ode_solver = new RK4Solver; break;
272  // MFEM implicit L-stable methods
273  case 5: ode_solver = new BackwardEulerSolver; break;
274  case 6: ode_solver = new SDIRK23Solver(2); break;
275  case 7: ode_solver = new SDIRK33Solver; break;
276  // CVODE
277  case 8:
278  cvode = new CVODESolver(CV_ADAMS);
279  cvode->Init(oper);
280  cvode->SetSStolerances(reltol, abstol);
281  cvode->SetMaxStep(dt);
282  ode_solver = cvode; break;
283  case 9:
284  cvode = new CVODESolver(CV_BDF);
285  cvode->Init(oper);
286  cvode->SetSStolerances(reltol, abstol);
287  cvode->SetMaxStep(dt);
288  ode_solver = cvode; break;
289  // ARKODE
290  case 10:
291  case 11:
292  arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
293  arkode->Init(oper);
294  arkode->SetSStolerances(reltol, abstol);
295  arkode->SetMaxStep(dt);
296  if (ode_solver_type == 11)
297  {
299  }
300  ode_solver = arkode; break;
301  case 12:
302  arkode = new ARKStepSolver(ARKStepSolver::IMPLICIT);
303  arkode->Init(oper);
304  arkode->SetSStolerances(reltol, abstol);
305  arkode->SetMaxStep(dt);
306  ode_solver = arkode; break;
307  }
308
309  // Initialize MFEM integrators, SUNDIALS integrators are initialized above
310  if (ode_solver_type < 8) { ode_solver->Init(oper); }
311
312  // Since we want to update the diffusion coefficient after every time step,
313  // we need to use the "one-step" mode of the SUNDIALS solvers.
314  if (cvode) { cvode->SetStepMode(CV_ONE_STEP); }
315  if (arkode) { arkode->SetStepMode(ARK_ONE_STEP); }
316
317  // 8. Perform time-integration (looping over the time iterations, ti, with a
318  // time-step dt).
319  cout << "Integrating the ODE ..." << endl;
320  tic_toc.Clear();
321  tic_toc.Start();
322
323  bool last_step = false;
324  for (int ti = 1; !last_step; ti++)
325  {
326  double dt_real = min(dt, t_final - t);
327
328  // Note that since we are using the "one-step" mode of the SUNDIALS
329  // solvers, they will, generally, step over the final time and will not
330  // explicitly perform the interpolation to t_final as they do in the
331  // "normal" step mode.
332
333  ode_solver->Step(u, t, dt_real);
334
335  last_step = (t >= t_final - 1e-8*dt);
336
337  if (last_step || (ti % vis_steps) == 0)
338  {
339  cout << "step " << ti << ", t = " << t << endl;
340  if (cvode) { cvode->PrintInfo(); }
341  if (arkode) { arkode->PrintInfo(); }
342
343  u_gf.SetFromTrueDofs(u);
344  if (visualization)
345  {
346  sout << "solution\n" << *mesh << u_gf << flush;
347  }
348
349  if (visit)
350  {
351  visit_dc.SetCycle(ti);
352  visit_dc.SetTime(t);
353  visit_dc.Save();
354  }
355  }
356  oper.SetParameters(u);
357  }
358  tic_toc.Stop();
359  cout << "Done, " << tic_toc.RealTime() << "s." << endl;
360
361  // 9. Save the final solution. This output can be viewed later using GLVis:
362  // "glvis -m ex16.mesh -g ex16-final.gf".
363  {
364  ofstream osol("ex16-final.gf");
365  osol.precision(precision);
366  u_gf.Save(osol);
367  }
368
369  // 10. Free the used memory.
370  delete ode_solver;
371  delete mesh;
372
373  return 0;
374 }
375
376 ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
377  double kap, const Vector &u)
378  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
379  T(NULL), z(height)
380 {
381  const double rel_tol = 1e-8;
382
383  M = new BilinearForm(&fespace);
385  M->Assemble();
386  M->FormSystemMatrix(ess_tdof_list, Mmat);
387
388  M_solver.iterative_mode = false;
389  M_solver.SetRelTol(rel_tol);
390  M_solver.SetAbsTol(0.0);
391  M_solver.SetMaxIter(50);
392  M_solver.SetPrintLevel(0);
393  M_solver.SetPreconditioner(M_prec);
394  M_solver.SetOperator(Mmat);
395
396  alpha = al;
397  kappa = kap;
398
399  T_solver.iterative_mode = false;
400  T_solver.SetRelTol(rel_tol);
401  T_solver.SetAbsTol(0.0);
402  T_solver.SetMaxIter(100);
403  T_solver.SetPrintLevel(0);
404  T_solver.SetPreconditioner(T_prec);
405
406  SetParameters(u);
407 }
408
409 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
410 {
411  // Compute:
412  // du_dt = M^{-1}*-K(u)
413  // for du_dt
414  Kmat.Mult(u, z);
415  z.Neg(); // z = -z
416  M_solver.Mult(z, du_dt);
417 }
418
419 void ConductionOperator::ImplicitSolve(const double dt,
420  const Vector &u, Vector &du_dt)
421 {
422  // Solve the equation:
423  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
424  // for du_dt
425  if (T) { delete T; }
426  T = Add(1.0, Mmat, dt, Kmat);
427  T_solver.SetOperator(*T);
428  Kmat.Mult(u, z);
429  z.Neg();
430  T_solver.Mult(z, du_dt);
431 }
432
433 void ConductionOperator::SetParameters(const Vector &u)
434 {
435  GridFunction u_alpha_gf(&fespace);
436  u_alpha_gf.SetFromTrueDofs(u);
437  for (int i = 0; i < u_alpha_gf.Size(); i++)
438  {
439  u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
440  }
441
442  delete K;
443  K = new BilinearForm(&fespace);
444
445  GridFunctionCoefficient u_coeff(&u_alpha_gf);
446
448  K->Assemble();
449  K->FormSystemMatrix(ess_tdof_list, Kmat);
450 }
451
452 int ConductionOperator::SUNImplicitSetup(const Vector &x,
453  const Vector &fx, int jok, int *jcur,
454  double gamma)
455 {
456  // Setup the ODE Jacobian T = M + gamma K.
457  if (T) { delete T; }
458  T = Add(1.0, Mmat, gamma, Kmat);
459  T_solver.SetOperator(*T);
460  *jcur = 1;
461  return (0);
462 }
463
464 int ConductionOperator::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
465 {
466  // Solve the system A x = z => (M - gamma K) x = M b.
467  Mmat.Mult(b, z);
468  T_solver.Mult(z, x);
469  return (0);
470 }
471
472 ConductionOperator::~ConductionOperator()
473 {
474  delete T;
475  delete M;
476  delete K;
477 }
478
479 double InitialTemperature(const Vector &x)
480 {
481  if (x.Norml2() < 0.5)
482  {
483  return 2.0;
484  }
485  else
486  {
487  return 1.0;
488  }
489 }
void Init(TimeDependentOperator &f_)
Initialize CVODE: calls CVodeCreate() to create the CVODE memory and set some defaults.
Definition: sundials.cpp:709
int visport
Definition: solvers.hpp:493
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
double InitialTemperature(const Vector &x)
Definition: ex16.cpp:385
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:875
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
Select the ARKode step mode: ARK_NORMAL (default) or ARK_ONE_STEP.
Definition: sundials.cpp:1708
void SetERKTableNum(ARKODE_ERKTableID table_id)
Choose a specific Butcher table for an explicit RK method.
Definition: sundials.cpp:1731
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
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]...
StopWatch tic_toc
Definition: tic_toc.cpp:447
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:23
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:1719
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition: tic_toc.cpp:429
double kappa
Definition: ex24.cpp:54
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:411
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:419
Interface to ARKode&#39;s ARKStep module – additive Runge-Kutta methods.
Definition: sundials.hpp:672
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2320
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition: gridfunc.cpp:366
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:385
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
Data collection with VisIt I/O routines.
void SetMaxStep(double dt_max)
Set the maximum time step.
Definition: sundials.cpp:893
void PrintInfo() const
Print various ARKStep statistics.
Definition: sundials.cpp:1756
void SetTime(double t)
Set physical time (for time-dependent simulations)
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:408
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:163
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:919
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
virtual void Save() override
Save the collection and a VisIt root file.
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:150
constexpr ARKODE_ERKTableID ARKODE_FEHLBERG_13_7_8
Definition: sundials.hpp:65
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2415
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:835
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
RefCoord t[3]
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
The classical forward Euler method.
Definition: ode.hpp:117
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:381
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
Definition: sundials.cpp:1474
int main(int argc, char *argv[])
Definition: ex16.cpp:92
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:1713
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:403
double f(const Vector &p)
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.