MFEM  v3.4
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/escher-p2.mesh
14 // mpirun -np 8 ex16p -m ../data/beam-tet.mesh -tf 10 -dt 0.1
15 // mpirun -np 4 ex16p -m ../data/amr-quad.mesh -o 4 -rs 0 -rp 0
16 // mpirun -np 4 ex16p -m ../data/amr-hex.mesh -o 2 -rs 0 -rp 0
17 //
18 // Description: This example solves a time dependent nonlinear heat equation
19 // problem of the form du/dt = C(u), with a non-linear diffusion
20 // operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
21 //
22 // The example demonstrates the use of nonlinear operators (the
23 // class ConductionOperator defining C(u)), as well as their
24 // implicit time integration. Note that implementing the method
25 // ConductionOperator::ImplicitSolve is the only requirement for
26 // high-order implicit (SDIRK) time integration.
27 //
28 // We recommend viewing examples 2, 9 and 10 before viewing this
29 // example.
30 
31 #include "mfem.hpp"
32 #include <fstream>
33 #include <iostream>
34 
35 using namespace std;
36 using namespace mfem;
37 
38 /** After spatial discretization, the conduction model can be written as:
39  *
40  * du/dt = M^{-1}(-Ku)
41  *
42  * where u is the vector representing the temperature, M is the mass matrix,
43  * and K is the diffusion operator with diffusivity depending on u:
44  * (\kappa + \alpha u).
45  *
46  * Class ConductionOperator represents the right-hand side of the above ODE.
47  */
48 class ConductionOperator : public TimeDependentOperator
49 {
50 protected:
51  ParFiniteElementSpace &fespace;
52  Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
53 
54  ParBilinearForm *M;
55  ParBilinearForm *K;
56 
57  HypreParMatrix Mmat;
58  HypreParMatrix Kmat;
59  HypreParMatrix *T; // T = M + dt K
60  double current_dt;
61 
62  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
63  HypreSmoother M_prec; // Preconditioner for the mass matrix M
64 
65  CGSolver T_solver; // Implicit solver for T = M + dt K
66  HypreSmoother T_prec; // Preconditioner for the implicit solver
67 
68  double alpha, kappa;
69 
70  mutable Vector z; // auxiliary vector
71 
72 public:
73  ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
74  const Vector &u);
75 
76  virtual void Mult(const Vector &u, Vector &du_dt) const;
77  /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
78  This is the only requirement for high-order SDIRK implicit integration.*/
79  virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
80 
81  /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
82  void SetParameters(const Vector &u);
83 
84  virtual ~ConductionOperator();
85 };
86 
87 double InitialTemperature(const Vector &x);
88 
89 int main(int argc, char *argv[])
90 {
91  // 1. Initialize MPI.
92  int num_procs, myid;
93  MPI_Init(&argc, &argv);
94  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
95  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
96 
97  // 2. Parse command-line options.
98  const char *mesh_file = "../data/star.mesh";
99  int ser_ref_levels = 2;
100  int par_ref_levels = 1;
101  int order = 2;
102  int ode_solver_type = 3;
103  double t_final = 0.5;
104  double dt = 1.0e-2;
105  double alpha = 1.0e-2;
106  double kappa = 0.5;
107  bool visualization = true;
108  bool visit = false;
109  int vis_steps = 5;
110 
111  int precision = 8;
112  cout.precision(precision);
113 
114  OptionsParser args(argc, argv);
115  args.AddOption(&mesh_file, "-m", "--mesh",
116  "Mesh file to use.");
117  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
118  "Number of times to refine the mesh uniformly in serial.");
119  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
120  "Number of times to refine the mesh uniformly in parallel.");
121  args.AddOption(&order, "-o", "--order",
122  "Order (degree) of the finite elements.");
123  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
124  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
125  "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
126  args.AddOption(&t_final, "-tf", "--t-final",
127  "Final time; start time is 0.");
128  args.AddOption(&dt, "-dt", "--time-step",
129  "Time step.");
130  args.AddOption(&alpha, "-a", "--alpha",
131  "Alpha coefficient.");
132  args.AddOption(&kappa, "-k", "--kappa",
133  "Kappa coefficient offset.");
134  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
135  "--no-visualization",
136  "Enable or disable GLVis visualization.");
137  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
138  "--no-visit-datafiles",
139  "Save data files for VisIt (visit.llnl.gov) visualization.");
140  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
141  "Visualize every n-th timestep.");
142  args.Parse();
143  if (!args.Good())
144  {
145  args.PrintUsage(cout);
146  MPI_Finalize();
147  return 1;
148  }
149 
150  if (myid == 0)
151  {
152  args.PrintOptions(cout);
153  }
154 
155  // 3. Read the serial mesh from the given mesh file on all processors. We can
156  // handle triangular, quadrilateral, tetrahedral and hexahedral meshes
157  // with the same code.
158  Mesh *mesh = new Mesh(mesh_file, 1, 1);
159  int dim = mesh->Dimension();
160 
161  // 4. Define the ODE solver used for time integration. Several implicit
162  // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
163  // explicit Runge-Kutta methods are available.
164  ODESolver *ode_solver;
165  switch (ode_solver_type)
166  {
167  // Implicit L-stable methods
168  case 1: ode_solver = new BackwardEulerSolver; break;
169  case 2: ode_solver = new SDIRK23Solver(2); break;
170  case 3: ode_solver = new SDIRK33Solver; break;
171  // Explicit methods
172  case 11: ode_solver = new ForwardEulerSolver; break;
173  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
174  case 13: ode_solver = new RK3SSPSolver; break;
175  case 14: ode_solver = new RK4Solver; break;
176  case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
177  // Implicit A-stable methods (not L-stable)
178  case 22: ode_solver = new ImplicitMidpointSolver; break;
179  case 23: ode_solver = new SDIRK23Solver; break;
180  case 24: ode_solver = new SDIRK34Solver; break;
181  default:
182  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
183  delete mesh;
184  return 3;
185  }
186 
187  // 5. Refine the mesh in serial to increase the resolution. In this example
188  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
189  // a command-line parameter.
190  for (int lev = 0; lev < ser_ref_levels; lev++)
191  {
192  mesh->UniformRefinement();
193  }
194 
195  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
196  // this mesh further in parallel to increase the resolution. Once the
197  // parallel mesh is defined, the serial mesh can be deleted.
198  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
199  delete mesh;
200  for (int lev = 0; lev < par_ref_levels; lev++)
201  {
202  pmesh->UniformRefinement();
203  }
204 
205  // 7. Define the vector finite element space representing the current and the
206  // initial temperature, u_ref.
207  H1_FECollection fe_coll(order, dim);
208  ParFiniteElementSpace fespace(pmesh, &fe_coll);
209 
210  int fe_size = fespace.GlobalTrueVSize();
211  if (myid == 0)
212  {
213  cout << "Number of temperature unknowns: " << fe_size << endl;
214  }
215 
216  ParGridFunction u_gf(&fespace);
217 
218  // 8. Set the initial conditions for u. All boundaries are considered
219  // natural.
221  u_gf.ProjectCoefficient(u_0);
222  Vector u;
223  u_gf.GetTrueDofs(u);
224 
225  // 9. Initialize the conduction operator and the VisIt visualization.
226  ConductionOperator oper(fespace, alpha, kappa, u);
227 
228  u_gf.SetFromTrueDofs(u);
229  {
230  ostringstream mesh_name, sol_name;
231  mesh_name << "ex16-mesh." << setfill('0') << setw(6) << myid;
232  sol_name << "ex16-init." << setfill('0') << setw(6) << myid;
233  ofstream omesh(mesh_name.str().c_str());
234  omesh.precision(precision);
235  pmesh->Print(omesh);
236  ofstream osol(sol_name.str().c_str());
237  osol.precision(precision);
238  u_gf.Save(osol);
239  }
240 
241  VisItDataCollection visit_dc("Example16-Parallel", pmesh);
242  visit_dc.RegisterField("temperature", &u_gf);
243  if (visit)
244  {
245  visit_dc.SetCycle(0);
246  visit_dc.SetTime(0.0);
247  visit_dc.Save();
248  }
249 
250  socketstream sout;
251  if (visualization)
252  {
253  char vishost[] = "localhost";
254  int visport = 19916;
255  sout.open(vishost, visport);
256  sout << "parallel " << num_procs << " " << myid << endl;
257  int good = sout.good(), all_good;
258  MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm());
259  if (!all_good)
260  {
261  sout.close();
262  visualization = false;
263  if (myid == 0)
264  {
265  cout << "Unable to connect to GLVis server at "
266  << vishost << ':' << visport << endl;
267  cout << "GLVis visualization disabled.\n";
268  }
269  }
270  else
271  {
272  sout.precision(precision);
273  sout << "solution\n" << *pmesh << u_gf;
274  sout << "pause\n";
275  sout << flush;
276  if (myid == 0)
277  {
278  cout << "GLVis visualization paused."
279  << " Press space (in the GLVis window) to resume it.\n";
280  }
281  }
282  }
283 
284  // 10. Perform time-integration (looping over the time iterations, ti, with a
285  // time-step dt).
286  ode_solver->Init(oper);
287  double t = 0.0;
288 
289  bool last_step = false;
290  for (int ti = 1; !last_step; ti++)
291  {
292  if (t + dt >= t_final - dt/2)
293  {
294  last_step = true;
295  }
296 
297  ode_solver->Step(u, t, dt);
298 
299  if (last_step || (ti % vis_steps) == 0)
300  {
301  if (myid == 0)
302  {
303  cout << "step " << ti << ", t = " << t << endl;
304  }
305 
306  u_gf.SetFromTrueDofs(u);
307  if (visualization)
308  {
309  sout << "parallel " << num_procs << " " << myid << "\n";
310  sout << "solution\n" << *pmesh << u_gf << flush;
311  }
312 
313  if (visit)
314  {
315  visit_dc.SetCycle(ti);
316  visit_dc.SetTime(t);
317  visit_dc.Save();
318  }
319  }
320  oper.SetParameters(u);
321  }
322 
323  // 11. Save the final solution in parallel. This output can be viewed later
324  // using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
325  {
326  ostringstream sol_name;
327  sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
328  ofstream osol(sol_name.str().c_str());
329  osol.precision(precision);
330  u_gf.Save(osol);
331  }
332 
333  // 12. Free the used memory.
334  delete ode_solver;
335  delete pmesh;
336 
337  MPI_Finalize();
338 
339  return 0;
340 }
341 
342 ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
343  double kap, const Vector &u)
344  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
345  T(NULL), current_dt(0.0),
346  M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
347 {
348  const double rel_tol = 1e-8;
349 
350  M = new ParBilinearForm(&fespace);
351  M->AddDomainIntegrator(new MassIntegrator());
352  M->Assemble(0); // keep sparsity pattern of M and K the same
353  M->FormSystemMatrix(ess_tdof_list, Mmat);
354 
355  M_solver.iterative_mode = false;
356  M_solver.SetRelTol(rel_tol);
357  M_solver.SetAbsTol(0.0);
358  M_solver.SetMaxIter(100);
359  M_solver.SetPrintLevel(0);
360  M_prec.SetType(HypreSmoother::Jacobi);
361  M_solver.SetPreconditioner(M_prec);
362  M_solver.SetOperator(Mmat);
363 
364  alpha = al;
365  kappa = kap;
366 
367  T_solver.iterative_mode = false;
368  T_solver.SetRelTol(rel_tol);
369  T_solver.SetAbsTol(0.0);
370  T_solver.SetMaxIter(100);
371  T_solver.SetPrintLevel(0);
372  T_solver.SetPreconditioner(T_prec);
373 
374  SetParameters(u);
375 }
376 
377 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
378 {
379  // Compute:
380  // du_dt = M^{-1}*-K(u)
381  // for du_dt
382  Kmat.Mult(u, z);
383  z.Neg(); // z = -z
384  M_solver.Mult(z, du_dt);
385 }
386 
387 void ConductionOperator::ImplicitSolve(const double dt,
388  const Vector &u, Vector &du_dt)
389 {
390  // Solve the equation:
391  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
392  // for du_dt
393  if (!T)
394  {
395  T = Add(1.0, Mmat, dt, Kmat);
396  current_dt = dt;
397  T_solver.SetOperator(*T);
398  }
399  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
400  Kmat.Mult(u, z);
401  z.Neg();
402  T_solver.Mult(z, du_dt);
403 }
404 
405 void ConductionOperator::SetParameters(const Vector &u)
406 {
407  ParGridFunction u_alpha_gf(&fespace);
408  u_alpha_gf.SetFromTrueDofs(u);
409  for (int i = 0; i < u_alpha_gf.Size(); i++)
410  {
411  u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
412  }
413 
414  delete K;
415  K = new ParBilinearForm(&fespace);
416 
417  GridFunctionCoefficient u_coeff(&u_alpha_gf);
418 
419  K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
420  K->Assemble(0); // keep sparsity pattern of M and K the same
421  K->FormSystemMatrix(ess_tdof_list, Kmat);
422  delete T;
423  T = NULL; // re-compute T on the next ImplicitSolve
424 }
425 
426 ConductionOperator::~ConductionOperator()
427 {
428  delete T;
429  delete M;
430  delete K;
431 }
432 
433 double InitialTemperature(const Vector &x)
434 {
435  if (x.Norml2() < 0.5)
436  {
437  return 2.0;
438  }
439  else
440  {
441  return 1.0;
442  }
443 }
Conjugate gradient method.
Definition: solvers.hpp:111
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
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]...
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.
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 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
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.
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
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
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
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.
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:225
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
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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
Class for parallel meshes.
Definition: pmesh.hpp:32