MFEM  v4.5.0
Finite element discretization library
ex16p.cpp
Go to the documentation of this file.
1 // MFEM Example 16 - Parallel Version
2 //
3 // Compile with: make ex16p
4 //
5 // Sample runs: mpirun -np 4 ex16p
6 // mpirun -np 4 ex16p -m ../data/inline-tri.mesh
7 // mpirun -np 4 ex16p -m ../data/disc-nurbs.mesh -tf 2
8 // mpirun -np 4 ex16p -s 1 -a 0.0 -k 1.0
9 // mpirun -np 4 ex16p -s 2 -a 1.0 -k 0.0
10 // mpirun -np 8 ex16p -s 3 -a 0.5 -k 0.5 -o 4
11 // mpirun -np 4 ex16p -s 14 -dt 1.0e-4 -tf 4.0e-2 -vs 40
12 // mpirun -np 16 ex16p -m ../data/fichera-q2.mesh
13 // mpirun -np 16 ex16p -m ../data/fichera-mixed.mesh
14 // mpirun -np 16 ex16p -m ../data/escher-p2.mesh
15 // mpirun -np 8 ex16p -m ../data/beam-tet.mesh -tf 10 -dt 0.1
16 // mpirun -np 4 ex16p -m ../data/amr-quad.mesh -o 4 -rs 0 -rp 0
17 // mpirun -np 4 ex16p -m ../data/amr-hex.mesh -o 2 -rs 0 -rp 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. In this example,
28 // the diffusion operator is linearized by evaluating with the
29 // lagged solution from the previous timestep, so there is only
30 // a linear solve. Optional saving with ADIOS2
31 // (adios2.readthedocs.io) is also illustrated.
32 //
33 // We recommend viewing examples 2, 9 and 10 before viewing this
34 // example.
35 
36 #include "mfem.hpp"
37 #include <fstream>
38 #include <iostream>
39 
40 using namespace std;
41 using namespace mfem;
42 
43 /** After spatial discretization, the conduction model can be written as:
44  *
45  * du/dt = M^{-1}(-Ku)
46  *
47  * where u is the vector representing the temperature, M is the mass matrix,
48  * and K is the diffusion operator with diffusivity depending on u:
49  * (\kappa + \alpha u).
50  *
51  * Class ConductionOperator represents the right-hand side of the above ODE.
52  */
53 class ConductionOperator : public TimeDependentOperator
54 {
55 protected:
56  ParFiniteElementSpace &fespace;
57  Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
58 
59  ParBilinearForm *M;
60  ParBilinearForm *K;
61 
62  HypreParMatrix Mmat;
63  HypreParMatrix Kmat;
64  HypreParMatrix *T; // T = M + dt K
65  double current_dt;
66 
67  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
68  HypreSmoother M_prec; // Preconditioner for the mass matrix M
69 
70  CGSolver T_solver; // Implicit solver for T = M + dt K
71  HypreSmoother T_prec; // Preconditioner for the implicit solver
72 
73  double alpha, kappa;
74 
75  mutable Vector z; // auxiliary vector
76 
77 public:
78  ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
79  const Vector &u);
80 
81  virtual void Mult(const Vector &u, Vector &du_dt) const;
82  /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
83  This is the only requirement for high-order SDIRK implicit integration.*/
84  virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
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 double InitialTemperature(const Vector &x);
93 
94 int main(int argc, char *argv[])
95 {
96  // 1. Initialize MPI and HYPRE.
97  Mpi::Init(argc, argv);
98  int num_procs = Mpi::WorldSize();
99  int myid = Mpi::WorldRank();
100  Hypre::Init();
101 
102  // 2. Parse command-line options.
103  const char *mesh_file = "../data/star.mesh";
104  int ser_ref_levels = 2;
105  int par_ref_levels = 1;
106  int order = 2;
107  int ode_solver_type = 3;
108  double t_final = 0.5;
109  double dt = 1.0e-2;
110  double alpha = 1.0e-2;
111  double kappa = 0.5;
112  bool visualization = true;
113  bool visit = false;
114  int vis_steps = 5;
115  bool adios2 = false;
116 
117  int precision = 8;
118  cout.precision(precision);
119 
120  OptionsParser args(argc, argv);
121  args.AddOption(&mesh_file, "-m", "--mesh",
122  "Mesh file to use.");
123  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
124  "Number of times to refine the mesh uniformly in serial.");
125  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
126  "Number of times to refine the mesh uniformly in parallel.");
127  args.AddOption(&order, "-o", "--order",
128  "Order (degree) of the finite elements.");
129  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
130  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
131  "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
132  args.AddOption(&t_final, "-tf", "--t-final",
133  "Final time; start time is 0.");
134  args.AddOption(&dt, "-dt", "--time-step",
135  "Time step.");
136  args.AddOption(&alpha, "-a", "--alpha",
137  "Alpha coefficient.");
138  args.AddOption(&kappa, "-k", "--kappa",
139  "Kappa coefficient offset.");
140  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
141  "--no-visualization",
142  "Enable or disable GLVis visualization.");
143  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
144  "--no-visit-datafiles",
145  "Save data files for VisIt (visit.llnl.gov) visualization.");
146  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
147  "Visualize every n-th timestep.");
148  args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
149  "--no-adios2-streams",
150  "Save data using adios2 streams.");
151  args.Parse();
152  if (!args.Good())
153  {
154  args.PrintUsage(cout);
155  return 1;
156  }
157 
158  if (myid == 0)
159  {
160  args.PrintOptions(cout);
161  }
162 
163  // 3. Read the serial mesh from the given mesh file on all processors. We can
164  // handle triangular, quadrilateral, tetrahedral and hexahedral meshes
165  // with the same code.
166  Mesh *mesh = new Mesh(mesh_file, 1, 1);
167  int dim = mesh->Dimension();
168 
169  // 4. Define the ODE solver used for time integration. Several implicit
170  // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
171  // explicit Runge-Kutta methods are available.
172  ODESolver *ode_solver;
173  switch (ode_solver_type)
174  {
175  // Implicit L-stable methods
176  case 1: ode_solver = new BackwardEulerSolver; break;
177  case 2: ode_solver = new SDIRK23Solver(2); break;
178  case 3: ode_solver = new SDIRK33Solver; break;
179  // Explicit methods
180  case 11: ode_solver = new ForwardEulerSolver; break;
181  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
182  case 13: ode_solver = new RK3SSPSolver; break;
183  case 14: ode_solver = new RK4Solver; break;
184  case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
185  // Implicit A-stable methods (not L-stable)
186  case 22: ode_solver = new ImplicitMidpointSolver; break;
187  case 23: ode_solver = new SDIRK23Solver; break;
188  case 24: ode_solver = new SDIRK34Solver; break;
189  default:
190  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
191  delete mesh;
192  return 3;
193  }
194 
195  // 5. Refine the mesh in serial to increase the resolution. In this example
196  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
197  // a command-line parameter.
198  for (int lev = 0; lev < ser_ref_levels; lev++)
199  {
200  mesh->UniformRefinement();
201  }
202 
203  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
204  // this mesh further in parallel to increase the resolution. Once the
205  // parallel mesh is defined, the serial mesh can be deleted.
206  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
207  delete mesh;
208  for (int lev = 0; lev < par_ref_levels; lev++)
209  {
210  pmesh->UniformRefinement();
211  }
212 
213  // 7. Define the vector finite element space representing the current and the
214  // initial temperature, u_ref.
215  H1_FECollection fe_coll(order, dim);
216  ParFiniteElementSpace fespace(pmesh, &fe_coll);
217 
218  HYPRE_BigInt fe_size = fespace.GlobalTrueVSize();
219  if (myid == 0)
220  {
221  cout << "Number of temperature unknowns: " << fe_size << endl;
222  }
223 
224  ParGridFunction u_gf(&fespace);
225 
226  // 8. Set the initial conditions for u. All boundaries are considered
227  // natural.
229  u_gf.ProjectCoefficient(u_0);
230  Vector u;
231  u_gf.GetTrueDofs(u);
232 
233  // 9. Initialize the conduction operator and the VisIt visualization.
234  ConductionOperator oper(fespace, alpha, kappa, u);
235 
236  u_gf.SetFromTrueDofs(u);
237  {
238  ostringstream mesh_name, sol_name;
239  mesh_name << "ex16-mesh." << setfill('0') << setw(6) << myid;
240  sol_name << "ex16-init." << setfill('0') << setw(6) << myid;
241  ofstream omesh(mesh_name.str().c_str());
242  omesh.precision(precision);
243  pmesh->Print(omesh);
244  ofstream osol(sol_name.str().c_str());
245  osol.precision(precision);
246  u_gf.Save(osol);
247  }
248 
249  VisItDataCollection visit_dc("Example16-Parallel", pmesh);
250  visit_dc.RegisterField("temperature", &u_gf);
251  if (visit)
252  {
253  visit_dc.SetCycle(0);
254  visit_dc.SetTime(0.0);
255  visit_dc.Save();
256  }
257 
258  // Optionally output a BP (binary pack) file using ADIOS2. This can be
259  // visualized with the ParaView VTX reader.
260 #ifdef MFEM_USE_ADIOS2
261  ADIOS2DataCollection* adios2_dc = NULL;
262  if (adios2)
263  {
264  std::string postfix(mesh_file);
265  postfix.erase(0, std::string("../data/").size() );
266  postfix += "_o" + std::to_string(order);
267  postfix += "_solver" + std::to_string(ode_solver_type);
268  const std::string collection_name = "ex16-p-" + postfix + ".bp";
269 
270  adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh);
271  adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) );
272  adios2_dc->RegisterField("temperature", &u_gf);
273  adios2_dc->SetCycle(0);
274  adios2_dc->SetTime(0.0);
275  adios2_dc->Save();
276  }
277 #endif
278 
279  socketstream sout;
280  if (visualization)
281  {
282  char vishost[] = "localhost";
283  int visport = 19916;
284  sout.open(vishost, visport);
285  sout << "parallel " << num_procs << " " << myid << endl;
286  int good = sout.good(), all_good;
287  MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm());
288  if (!all_good)
289  {
290  sout.close();
291  visualization = false;
292  if (myid == 0)
293  {
294  cout << "Unable to connect to GLVis server at "
295  << vishost << ':' << visport << endl;
296  cout << "GLVis visualization disabled.\n";
297  }
298  }
299  else
300  {
301  sout.precision(precision);
302  sout << "solution\n" << *pmesh << u_gf;
303  sout << "pause\n";
304  sout << flush;
305  if (myid == 0)
306  {
307  cout << "GLVis visualization paused."
308  << " Press space (in the GLVis window) to resume it.\n";
309  }
310  }
311  }
312 
313  // 10. Perform time-integration (looping over the time iterations, ti, with a
314  // time-step dt).
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  if (t + dt >= t_final - dt/2)
322  {
323  last_step = true;
324  }
325 
326  ode_solver->Step(u, t, dt);
327 
328  if (last_step || (ti % vis_steps) == 0)
329  {
330  if (myid == 0)
331  {
332  cout << "step " << ti << ", t = " << t << endl;
333  }
334 
335  u_gf.SetFromTrueDofs(u);
336  if (visualization)
337  {
338  sout << "parallel " << num_procs << " " << myid << "\n";
339  sout << "solution\n" << *pmesh << u_gf << flush;
340  }
341 
342  if (visit)
343  {
344  visit_dc.SetCycle(ti);
345  visit_dc.SetTime(t);
346  visit_dc.Save();
347  }
348 
349 #ifdef MFEM_USE_ADIOS2
350  if (adios2)
351  {
352  adios2_dc->SetCycle(ti);
353  adios2_dc->SetTime(t);
354  adios2_dc->Save();
355  }
356 #endif
357  }
358  oper.SetParameters(u);
359  }
360 
361 #ifdef MFEM_USE_ADIOS2
362  if (adios2)
363  {
364  delete adios2_dc;
365  }
366 #endif
367 
368  // 11. Save the final solution in parallel. This output can be viewed later
369  // using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
370  {
371  ostringstream sol_name;
372  sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
373  ofstream osol(sol_name.str().c_str());
374  osol.precision(precision);
375  u_gf.Save(osol);
376  }
377 
378  // 12. Free the used memory.
379  delete ode_solver;
380  delete pmesh;
381 
382  return 0;
383 }
384 
385 ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
386  double kap, const Vector &u)
387  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
388  T(NULL), current_dt(0.0),
389  M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
390 {
391  const double rel_tol = 1e-8;
392 
393  M = new ParBilinearForm(&fespace);
394  M->AddDomainIntegrator(new MassIntegrator());
395  M->Assemble(0); // keep sparsity pattern of M and K the same
396  M->FormSystemMatrix(ess_tdof_list, Mmat);
397 
398  M_solver.iterative_mode = false;
399  M_solver.SetRelTol(rel_tol);
400  M_solver.SetAbsTol(0.0);
401  M_solver.SetMaxIter(100);
402  M_solver.SetPrintLevel(0);
403  M_prec.SetType(HypreSmoother::Jacobi);
404  M_solver.SetPreconditioner(M_prec);
405  M_solver.SetOperator(Mmat);
406 
407  alpha = al;
408  kappa = kap;
409 
410  T_solver.iterative_mode = false;
411  T_solver.SetRelTol(rel_tol);
412  T_solver.SetAbsTol(0.0);
413  T_solver.SetMaxIter(100);
414  T_solver.SetPrintLevel(0);
415  T_solver.SetPreconditioner(T_prec);
416 
417  SetParameters(u);
418 }
419 
420 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
421 {
422  // Compute:
423  // du_dt = M^{-1}*-Ku
424  // for du_dt, where K is linearized by using u from the previous timestep
425  Kmat.Mult(u, z);
426  z.Neg(); // z = -z
427  M_solver.Mult(z, du_dt);
428 }
429 
430 void ConductionOperator::ImplicitSolve(const double dt,
431  const Vector &u, Vector &du_dt)
432 {
433  // Solve the equation:
434  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
435  // for du_dt, where K is linearized by using u from the previous timestep
436  if (!T)
437  {
438  T = Add(1.0, Mmat, dt, Kmat);
439  current_dt = dt;
440  T_solver.SetOperator(*T);
441  }
442  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
443  Kmat.Mult(u, z);
444  z.Neg();
445  T_solver.Mult(z, du_dt);
446 }
447 
448 void ConductionOperator::SetParameters(const Vector &u)
449 {
450  ParGridFunction u_alpha_gf(&fespace);
451  u_alpha_gf.SetFromTrueDofs(u);
452  for (int i = 0; i < u_alpha_gf.Size(); i++)
453  {
454  u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
455  }
456 
457  delete K;
458  K = new ParBilinearForm(&fespace);
459 
460  GridFunctionCoefficient u_coeff(&u_alpha_gf);
461 
462  K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
463  K->Assemble(0); // keep sparsity pattern of M and K the same
464  K->FormSystemMatrix(ess_tdof_list, Kmat);
465  delete T;
466  T = NULL; // re-compute T on the next ImplicitSolve
467 }
468 
469 ConductionOperator::~ConductionOperator()
470 {
471  delete T;
472  delete M;
473  delete K;
474 }
475 
476 double InitialTemperature(const Vector &x)
477 {
478  if (x.Norml2() < 0.5)
479  {
480  return 2.0;
481  }
482  else
483  {
484  return 1.0;
485  }
486 }
Conjugate gradient method.
Definition: solvers.hpp:465
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
Base abstract class for first order time dependent operators.
Definition: operator.hpp:289
int main(int argc, char *argv[])
Definition: ex16p.cpp:94
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]...
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.
double kappa
Definition: ex24.cpp:54
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
STL namespace.
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:391
int close()
Close the socketstream.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2268
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9802
constexpr int visport
Data collection with VisIt I/O routines.
void SetParameter(const std::string key, const std::string value) noexcept
Parallel smoothers in hypre.
Definition: hypre.hpp:962
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
int precision
Precision (number of digits) used for the text output of doubles.
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:149
HYPRE_Int HYPRE_BigInt
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:404
int dim
Definition: ex24.cpp:53
Class for parallel bilinear form.
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:812
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:60
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
double InitialTemperature(const Vector &x)
Definition: ex16p.cpp:476
The classical forward Euler method.
Definition: ode.hpp:116
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
double f(const Vector &p)