MFEM  v3.3.2
Finite element discretization library
sundials/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  return 3;
248  }
249 
250  // Since we want to update the diffusion coefficient after every time step,
251  // we need to use the "one-step" mode of the SUNDIALS solvers.
252  if (cvode) { cvode->SetStepMode(CV_ONE_STEP); }
253  if (arkode) { arkode->SetStepMode(ARK_ONE_STEP); }
254 
255  // 5. Refine the mesh in serial to increase the resolution. In this example
256  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
257  // a command-line parameter.
258  for (int lev = 0; lev < ser_ref_levels; lev++)
259  {
260  mesh->UniformRefinement();
261  }
262 
263  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
264  // this mesh further in parallel to increase the resolution. Once the
265  // parallel mesh is defined, the serial mesh can be deleted.
266  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
267  delete mesh;
268  for (int lev = 0; lev < par_ref_levels; lev++)
269  {
270  pmesh->UniformRefinement();
271  }
272 
273  // 7. Define the vector finite element space representing the current and the
274  // initial temperature, u_ref.
275  H1_FECollection fe_coll(order, dim);
276  ParFiniteElementSpace fespace(pmesh, &fe_coll);
277 
278  int fe_size = fespace.GlobalTrueVSize();
279  if (myid == 0)
280  {
281  cout << "Number of temperature unknowns: " << fe_size << endl;
282  }
283 
284  ParGridFunction u_gf(&fespace);
285 
286  // 8. Set the initial conditions for u. All boundaries are considered
287  // natural.
289  u_gf.ProjectCoefficient(u_0);
290  Vector u;
291  u_gf.GetTrueDofs(u);
292 
293  // 9. Initialize the conduction operator and the VisIt visualization.
294  ConductionOperator oper(fespace, alpha, kappa, u);
295 
296  u_gf.SetFromTrueDofs(u);
297  {
298  ostringstream mesh_name, sol_name;
299  mesh_name << "ex16-mesh." << setfill('0') << setw(6) << myid;
300  sol_name << "ex16-init." << setfill('0') << setw(6) << myid;
301  ofstream omesh(mesh_name.str().c_str());
302  omesh.precision(precision);
303  pmesh->Print(omesh);
304  ofstream osol(sol_name.str().c_str());
305  osol.precision(precision);
306  u_gf.Save(osol);
307  }
308 
309  VisItDataCollection visit_dc("Example16-Parallel", pmesh);
310  visit_dc.RegisterField("temperature", &u_gf);
311  if (visit)
312  {
313  visit_dc.SetCycle(0);
314  visit_dc.SetTime(0.0);
315  visit_dc.Save();
316  }
317 
318  socketstream sout;
319  if (visualization)
320  {
321  char vishost[] = "localhost";
322  int visport = 19916;
323  sout.open(vishost, visport);
324  sout << "parallel " << num_procs << " " << myid << endl;
325  int good = sout.good(), all_good;
326  MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm());
327  if (!all_good)
328  {
329  sout.close();
330  visualization = false;
331  if (myid == 0)
332  {
333  cout << "Unable to connect to GLVis server at "
334  << vishost << ':' << visport << endl;
335  cout << "GLVis visualization disabled.\n";
336  }
337  }
338  else
339  {
340  sout.precision(precision);
341  sout << "solution\n" << *pmesh << u_gf;
342  sout << "pause\n";
343  sout << flush;
344  if (myid == 0)
345  {
346  cout << "GLVis visualization paused."
347  << " Press space (in the GLVis window) to resume it.\n";
348  }
349  }
350  }
351 
352  // 10. Perform time-integration (looping over the time iterations, ti, with a
353  // time-step dt).
354  if (myid == 0)
355  {
356  cout << "Integrating the ODE ..." << endl;
357  }
358  tic_toc.Clear();
359  tic_toc.Start();
360  ode_solver->Init(oper);
361  double t = 0.0;
362 
363  bool last_step = false;
364  for (int ti = 1; !last_step; ti++)
365  {
366  double dt_real = min(dt, t_final - t);
367 
368  // Note that since we are using the "one-step" mode of the SUNDIALS
369  // solvers, they will, generally, step over the final time and will not
370  // explicitly perform the interpolation to t_final as they do in the
371  // "normal" step mode.
372 
373  ode_solver->Step(u, t, dt_real);
374 
375  last_step = (t >= t_final - 1e-8*dt);
376 
377  if (last_step || (ti % vis_steps) == 0)
378  {
379  if (myid == 0)
380  {
381  cout << "step " << ti << ", t = " << t << endl;
382  if (cvode) { cvode->PrintInfo(); }
383  if (arkode) { arkode->PrintInfo(); }
384  }
385 
386  u_gf.SetFromTrueDofs(u);
387  if (visualization)
388  {
389  sout << "parallel " << num_procs << " " << myid << "\n";
390  sout << "solution\n" << *pmesh << u_gf << flush;
391  }
392 
393  if (visit)
394  {
395  visit_dc.SetCycle(ti);
396  visit_dc.SetTime(t);
397  visit_dc.Save();
398  }
399  }
400  oper.SetParameters(u);
401  }
402  tic_toc.Stop();
403  if (myid == 0)
404  {
405  cout << "Done, " << tic_toc.RealTime() << "s." << endl;
406  }
407 
408  // 11. Save the final solution in parallel. This output can be viewed later
409  // using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
410  {
411  ostringstream sol_name;
412  sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
413  ofstream osol(sol_name.str().c_str());
414  osol.precision(precision);
415  u_gf.Save(osol);
416  }
417 
418  // 12. Free the used memory.
419  delete ode_solver;
420  delete pmesh;
421 
422  MPI_Finalize();
423 
424  return 0;
425 }
426 
427 ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
428  double kap, const Vector &u)
429  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
430  T(NULL), current_dt(0.0),
431  M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
432 {
433  const double rel_tol = 1e-8;
434 
435  M = new ParBilinearForm(&fespace);
436  M->AddDomainIntegrator(new MassIntegrator());
437  M->Assemble(0); // keep sparsity pattern of M and K the same
438  M->FormSystemMatrix(ess_tdof_list, Mmat);
439 
440  M_solver.iterative_mode = false;
441  M_solver.SetRelTol(rel_tol);
442  M_solver.SetAbsTol(0.0);
443  M_solver.SetMaxIter(100);
444  M_solver.SetPrintLevel(0);
445  M_prec.SetType(HypreSmoother::Jacobi);
446  M_solver.SetPreconditioner(M_prec);
447  M_solver.SetOperator(Mmat);
448 
449  alpha = al;
450  kappa = kap;
451 
452  T_solver.iterative_mode = false;
453  T_solver.SetRelTol(rel_tol);
454  T_solver.SetAbsTol(0.0);
455  T_solver.SetMaxIter(100);
456  T_solver.SetPrintLevel(0);
457  T_solver.SetPreconditioner(T_prec);
458 
459  SetParameters(u);
460 }
461 
462 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
463 {
464  // Compute:
465  // du_dt = M^{-1}*-K(u)
466  // for du_dt
467  Kmat.Mult(u, z);
468  z.Neg(); // z = -z
469  M_solver.Mult(z, du_dt);
470 }
471 
472 void ConductionOperator::ImplicitSolve(const double dt,
473  const Vector &u, Vector &du_dt)
474 {
475  // Solve the equation:
476  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
477  // for du_dt
478  if (!T)
479  {
480  T = Add(1.0, Mmat, dt, Kmat);
481  current_dt = dt;
482  T_solver.SetOperator(*T);
483  }
484  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
485  Kmat.Mult(u, z);
486  z.Neg();
487  T_solver.Mult(z, du_dt);
488 }
489 
490 void ConductionOperator::SundialsSolve(const double dt, Vector &b)
491 {
492  // Solve the system (M + dt K) y = M b. The result y replaces the input b.
493  if (!T || dt != current_dt)
494  {
495  delete T;
496  T = Add(1.0, Mmat, dt, Kmat);
497  current_dt = dt;
498  T_solver.SetOperator(*T);
499  }
500  Mmat.Mult(b, z);
501  T_solver.Mult(z, b);
502 }
503 
504 void ConductionOperator::SetParameters(const Vector &u)
505 {
506  ParGridFunction u_alpha_gf(&fespace);
507  u_alpha_gf.SetFromTrueDofs(u);
508  for (int i = 0; i < u_alpha_gf.Size(); i++)
509  {
510  u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
511  }
512 
513  delete K;
514  K = new ParBilinearForm(&fespace);
515 
516  GridFunctionCoefficient u_coeff(&u_alpha_gf);
517 
518  K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
519  K->Assemble(0); // keep sparsity pattern of M and K the same
520  K->FormSystemMatrix(ess_tdof_list, Kmat);
521  delete T;
522  T = NULL; // re-compute T on the next ImplicitSolve or SundialsSolve
523 }
524 
525 ConductionOperator::~ConductionOperator()
526 {
527  delete T;
528  delete M;
529  delete K;
530 }
531 
532 
533 int SundialsJacSolver::InitSystem(void *sundials_mem)
534 {
535  TimeDependentOperator *td_oper = GetTimeDependentOperator(sundials_mem);
536 
537  // During development, we use dynamic_cast<> to ensure the setup is correct:
538  oper = dynamic_cast<ConductionOperator*>(td_oper);
539  MFEM_VERIFY(oper, "operator is not ConductionOperator");
540 
541  // When the implementation is finalized, we can switch to static_cast<>:
542  // oper = static_cast<ConductionOperator*>(td_oper);
543 
544  return 0;
545 }
546 
547 int SundialsJacSolver::SetupSystem(void *sundials_mem, int conv_fail,
548  const Vector &y_pred, const Vector &f_pred,
549  int &jac_cur, Vector &v_temp1,
550  Vector &v_temp2, Vector &v_temp3)
551 {
552  jac_cur = 1;
553 
554  return 0;
555 }
556 
557 int SundialsJacSolver::SolveSystem(void *sundials_mem, Vector &b,
558  const Vector &weight, const Vector &y_cur,
559  const Vector &f_cur)
560 {
561  oper->SundialsSolve(GetTimeStep(sundials_mem), b);
562 
563  return 0;
564 }
565 
566 int SundialsJacSolver::FreeSystem(void *sundials_mem)
567 {
568  return 0;
569 }
570 
571 
572 double InitialTemperature(const Vector &x)
573 {
574  if (x.Norml2() < 0.5)
575  {
576  return 2.0;
577  }
578  else
579  {
580  return 1.0;
581  }
582 }
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
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:142
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:670
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.
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:361
StopWatch tic_toc
Definition: tic_toc.cpp:447
double RealTime()
Definition: tic_toc.cpp:426
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Definition: sundials.cpp:233
Abstract parallel finite element space.
Definition: pfespace.hpp:31
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:275
STL namespace.
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:130
int main(int argc, char *argv[])
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:6718
Data collection with VisIt I/O routines.
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:178
void SetMaxStep(double dt_max)
Set the maximum time step of the linear multistep method.
Definition: sundials.hpp:183
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:3399
Parallel smoothers in hypre.
Definition: hypre.hpp:504
int Dimension() const
Definition: mesh.hpp:641
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
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:486
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
MPI_Comm GetComm() const
Definition: pmesh.hpp:117
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
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
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:126
Class for parallel bilinear form.
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:41
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
void PrintInfo() const
Print CVODE statistics.
Definition: sundials.cpp:400
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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's ParCSR matrix class.
Definition: hypre.hpp:175
double InitialTemperature(const Vector &x)
Class for parallel meshes.
Definition: pmesh.hpp:29
void PrintInfo() const
Print ARKODE statistics.
Definition: sundials.cpp:691
void SetStepMode(int itask)
CVode supports two modes, specified by itask: CV_NORMAL (default) and CV_ONE_STEP.
Definition: sundials.cpp:254
bool Good() const
Definition: optparser.hpp:120