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