MFEM  v4.1.0 Finite element discretization library
ex23.cpp
Go to the documentation of this file.
1 // MFEM Example 23
2 //
3 // Compile with: make ex23
4 //
5 // Sample runs: ex23
6 // ex23 -o 4 -tf 5
7 // ex23 -m ../data/square-disc.mesh -o 2 -tf 2 --neumann
8 // ex23 -m ../data/disc-nurbs.mesh -r 3 -o 4 -tf 2
9 // ex23 -m ../data/inline-hex.mesh -o 1 -tf 2 --neumann
10 // ex23 -m ../data/inline-tet.mesh -o 1 -tf 2 --neumann
11 //
12 // Description: This example solves the wave equation problem of the form:
13 //
14 // d^2u/dt^2 = c^2 \Delta u.
15 //
16 // The example demonstrates the use of time dependent operators,
17 // implicit solvers and second order time integration.
18 //
19 // We recommend viewing examples 9 and 10 before viewing this
20 // example.
21
22 #include "mfem.hpp"
23 #include <fstream>
24 #include <iostream>
25
26 using namespace std;
27 using namespace mfem;
28
29 /** After spatial discretization, the conduction model can be written as:
30  *
31  * d^2u/dt^2 = M^{-1}(-Ku)
32  *
33  * where u is the vector representing the temperature, M is the mass matrix,
34  * and K is the diffusion operator with diffusivity depending on u:
35  * (\kappa + \alpha u).
36  *
37  * Class WaveOperator represents the right-hand side of the above ODE.
38  */
39 class WaveOperator : public SecondOrderTimeDependentOperator
40 {
41 protected:
42  FiniteElementSpace &fespace;
43  Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
44
45  BilinearForm *M;
46  BilinearForm *K;
47
48  SparseMatrix Mmat, Kmat, Kmat0;
49  SparseMatrix *T; // T = M + dt K
50  double current_dt;
51
52  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
53  DSmoother M_prec; // Preconditioner for the mass matrix M
54
55  CGSolver T_solver; // Implicit solver for T = M + fac0*K
56  DSmoother T_prec; // Preconditioner for the implicit solver
57
58  Coefficient *c2;
59  mutable Vector z; // auxiliary vector
60
61 public:
62  WaveOperator(FiniteElementSpace &f, Array<int> &ess_bdr,double speed);
63
65  virtual void Mult(const Vector &u, const Vector &du_dt,
66  Vector &d2udt2) const;
67
68  /** Solve the Backward-Euler equation:
69  d2udt2 = f(u + fac0*d2udt2,dudt + fac1*d2udt2, t),
70  for the unknown d2udt2. */
71  using SecondOrderTimeDependentOperator::ImplicitSolve;
72  virtual void ImplicitSolve(const double fac0, const double fac1,
73  const Vector &u, const Vector &dudt, Vector &d2udt2);
74
75  ///
76  void SetParameters(const Vector &u);
77
78  virtual ~WaveOperator();
79 };
80
81
82 WaveOperator::WaveOperator(FiniteElementSpace &f,
83  Array<int> &ess_bdr, double speed)
84  : SecondOrderTimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL),
85  K(NULL),
86  T(NULL), current_dt(0.0), z(height)
87 {
88  const double rel_tol = 1e-8;
89
90  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
91
92  c2 = new ConstantCoefficient(speed*speed);
93
94  K = new BilinearForm(&fespace);
96  K->Assemble();
97
98  Array<int> dummy;
99  K->FormSystemMatrix(dummy, Kmat0);
100  K->FormSystemMatrix(ess_tdof_list, Kmat);
101
102  M = new BilinearForm(&fespace);
104  M->Assemble();
105  M->FormSystemMatrix(ess_tdof_list, Mmat);
106
107  M_solver.iterative_mode = false;
108  M_solver.SetRelTol(rel_tol);
109  M_solver.SetAbsTol(0.0);
110  M_solver.SetMaxIter(30);
111  M_solver.SetPrintLevel(0);
112  M_solver.SetPreconditioner(M_prec);
113  M_solver.SetOperator(Mmat);
114
115  T_solver.iterative_mode = false;
116  T_solver.SetRelTol(rel_tol);
117  T_solver.SetAbsTol(0.0);
118  T_solver.SetMaxIter(100);
119  T_solver.SetPrintLevel(0);
120  T_solver.SetPreconditioner(T_prec);
121
122  T = NULL;
123 }
124
125 void WaveOperator::Mult(const Vector &u, const Vector &du_dt,
126  Vector &d2udt2) const
127 {
128  // Compute:
129  // d2udt2 = M^{-1}*-K(u)
130  // for d2udt2
131  Kmat.Mult(u, z);
132  z.Neg(); // z = -z
133  M_solver.Mult(z, d2udt2);
134 }
135
136 void WaveOperator::ImplicitSolve(const double fac0, const double fac1,
137  const Vector &u, const Vector &dudt, Vector &d2udt2)
138 {
139  // Solve the equation:
140  // d2udt2 = M^{-1}*[-K(u + fac0*d2udt2)]
141  // for d2udt2
142  if (!T)
143  {
144  T = Add(1.0, Mmat, fac0, Kmat);
145  T_solver.SetOperator(*T);
146  }
147  Kmat0.Mult(u, z);
148  z.Neg();
149
150  for (int i = 0; i < ess_tdof_list.Size(); i++)
151  {
152  z[ess_tdof_list[i]] = 0.0;
153  }
154  T_solver.Mult(z, d2udt2);
155 }
156
157 void WaveOperator::SetParameters(const Vector &u)
158 {
159  delete T;
160  T = NULL; // re-compute T on the next ImplicitSolve
161 }
162
163 WaveOperator::~WaveOperator()
164 {
165  delete T;
166  delete M;
167  delete K;
168  delete c2;
169 }
170
171 double InitialSolution(const Vector &x)
172 {
173  return exp(-x.Norml2()*x.Norml2()*30);
174 }
175
176 double InitialRate(const Vector &x)
177 {
178  return 0.0;
179 }
180
181
182 int main(int argc, char *argv[])
183 {
184  // 1. Parse command-line options.
185  const char *mesh_file = "../data/star.mesh";
186  const char *ref_dir = "";
187  int ref_levels = 2;
188  int order = 2;
189  int ode_solver_type = 10;
190  double t_final = 0.5;
191  double dt = 1.0e-2;
192  double speed = 1.0;
193  bool visualization = true;
194  bool visit = true;
195  bool dirichlet = true;
196  int vis_steps = 5;
197
198  int precision = 8;
199  cout.precision(precision);
200
201  OptionsParser args(argc, argv);
203  "Mesh file to use.");
205  "Number of times to refine the mesh uniformly.");
207  "Order (degree) of the finite elements.");
209  "ODE solver: [0--10] - GeneralizedAlpha(0.1 * s),\n\t"
210  "\t 11 - Average Acceleration, 12 - Linear Acceleration\n"
211  "\t 13 - CentralDifference, 14 - FoxGoodwin");
213  "Final time; start time is 0.");
215  "Time step.");
217  "Wave speed.");
219  "--neumann",
220  "BC switch.");
222  "Reference directory for checking final solution.");
224  "--no-visualization",
225  "Enable or disable GLVis visualization.");
227  "--no-visit-datafiles",
228  "Save data files for VisIt (visit.llnl.gov) visualization.");
230  "Visualize every n-th timestep.");
231  args.Parse();
232  if (!args.Good())
233  {
234  args.PrintUsage(cout);
235  return 1;
236  }
237  args.PrintOptions(cout);
238
239  // 2. Read the mesh from the given mesh file. We can handle triangular,
240  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
241  Mesh *mesh = new Mesh(mesh_file, 1, 1);
242  int dim = mesh->Dimension();
243
244  // 3. Define the ODE solver used for time integration. Several second order
245  // time integrators are available.
246  SecondOrderODESolver *ode_solver;
247  switch (ode_solver_type)
248  {
249  // Implicit methods
250  case 0: ode_solver = new GeneralizedAlpha2Solver(0.0); break;
251  case 1: ode_solver = new GeneralizedAlpha2Solver(0.1); break;
252  case 2: ode_solver = new GeneralizedAlpha2Solver(0.2); break;
253  case 3: ode_solver = new GeneralizedAlpha2Solver(0.3); break;
254  case 4: ode_solver = new GeneralizedAlpha2Solver(0.4); break;
255  case 5: ode_solver = new GeneralizedAlpha2Solver(0.5); break;
256  case 6: ode_solver = new GeneralizedAlpha2Solver(0.6); break;
257  case 7: ode_solver = new GeneralizedAlpha2Solver(0.7); break;
258  case 8: ode_solver = new GeneralizedAlpha2Solver(0.8); break;
259  case 9: ode_solver = new GeneralizedAlpha2Solver(0.9); break;
260  case 10: ode_solver = new GeneralizedAlpha2Solver(1.0); break;
261
262  case 11: ode_solver = new AverageAccelerationSolver(); break;
263  case 12: ode_solver = new LinearAccelerationSolver(); break;
264  case 13: ode_solver = new CentralDifferenceSolver(); break;
265  case 14: ode_solver = new FoxGoodwinSolver(); break;
266
267  default:
268  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
269  delete mesh;
270  return 3;
271  }
272
273  // 4. Refine the mesh to increase the resolution. In this example we do
274  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
275  // command-line parameter.
276  for (int lev = 0; lev < ref_levels; lev++)
277  {
278  mesh->UniformRefinement();
279  }
280
281  // 5. Define the vector finite element space representing the current and the
282  // initial temperature, u_ref.
283  H1_FECollection fe_coll(order, dim);
284  FiniteElementSpace fespace(mesh, &fe_coll);
285
286  int fe_size = fespace.GetTrueVSize();
287  cout << "Number of temperature unknowns: " << fe_size << endl;
288
289  GridFunction u_gf(&fespace);
290  GridFunction dudt_gf(&fespace);
291  // 6. Set the initial conditions for u. All boundaries are considered
292  // natural.
294  u_gf.ProjectCoefficient(u_0);
295  Vector u;
296  u_gf.GetTrueDofs(u);
297
299  dudt_gf.ProjectCoefficient(dudt_0);
300  Vector dudt;
301  dudt_gf.GetTrueDofs(dudt);
302
303  // 7. Initialize the conduction operator and the visualization.
304  Array<int> ess_bdr;
305  if (mesh->bdr_attributes.Size())
306  {
307  ess_bdr.SetSize(mesh->bdr_attributes.Max());
308
309  if (dirichlet)
310  {
311  ess_bdr = 1;
312  }
313  else
314  {
315  ess_bdr = 0;
316  }
317  }
318
319  WaveOperator oper(fespace, ess_bdr, speed);
320
321  u_gf.SetFromTrueDofs(u);
322  {
323  ofstream omesh("ex23.mesh");
324  omesh.precision(precision);
325  mesh->Print(omesh);
326  ofstream osol("ex23-init.gf");
327  osol.precision(precision);
328  u_gf.Save(osol);
329  dudt_gf.Save(osol);
330  }
331
332  VisItDataCollection visit_dc("Example23", mesh);
333  visit_dc.RegisterField("solution", &u_gf);
334  visit_dc.RegisterField("rate", &dudt_gf);
335  if (visit)
336  {
337  visit_dc.SetCycle(0);
338  visit_dc.SetTime(0.0);
339  visit_dc.Save();
340  }
341
342  socketstream sout;
343  if (visualization)
344  {
345  char vishost[] = "localhost";
346  int visport = 19916;
347  sout.open(vishost, visport);
348  if (!sout)
349  {
350  cout << "Unable to connect to GLVis server at "
351  << vishost << ':' << visport << endl;
352  visualization = false;
353  cout << "GLVis visualization disabled.\n";
354  }
355  else
356  {
357  sout.precision(precision);
358  sout << "solution\n" << *mesh << dudt_gf;
359  sout << "pause\n";
360  sout << flush;
361  cout << "GLVis visualization paused."
362  << " Press space (in the GLVis window) to resume it.\n";
363  }
364  }
365
366  // 8. Perform time-integration (looping over the time iterations, ti, with a
367  // time-step dt).
368  ode_solver->Init(oper);
369  double t = 0.0;
370
371  bool last_step = false;
372  for (int ti = 1; !last_step; ti++)
373  {
374
375  if (t + dt >= t_final - dt/2)
376  {
377  last_step = true;
378  }
379
380  ode_solver->Step(u, dudt, t, dt);
381
382  if (last_step || (ti % vis_steps) == 0)
383  {
384  cout << "step " << ti << ", t = " << t << endl;
385
386  u_gf.SetFromTrueDofs(u);
387  dudt_gf.SetFromTrueDofs(dudt);
388  if (visualization)
389  {
390  sout << "solution\n" << *mesh << 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
403  // 9. Save the final solution. This output can be viewed later using GLVis:
404  // "glvis -m ex23.mesh -g ex23-final.gf".
405  {
406  ofstream osol("ex23-final.gf");
407  osol.precision(precision);
408  u_gf.Save(osol);
409  dudt_gf.Save(osol);
410  }
411
412  // 10. Free the used memory.
413  delete ode_solver;
414  delete mesh;
415
416  return 0;
417 }
int Size() const
Logical size of the array.
Definition: array.hpp:124
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
Definition: solvers.hpp:152
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
Subclass constant coefficient.
Definition: coefficient.hpp:67
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition: operator.cpp:230
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:711
virtual void Save()
Save the collection and a VisIt root file.
Abstract class for solving systems of ODEs: d2x/dt2 = f(x,dx/dt,t)
Definition: ode.hpp:528
int main(int argc, char *argv[])
Definition: ex1.cpp:62
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2731
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1926
virtual void Step(Vector &x, Vector &dxdt, 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]...
Data type sparse matrix.
Definition: sparsemat.hpp:40
The classical midpoint method.
Definition: ode.hpp:686
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
Data collection with VisIt I/O routines.
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:817
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction. If all dofs are true, then tv will be set to point to th...
Definition: gridfunc.cpp:323
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:384
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
double InitialSolution(const Vector &x)
Definition: ex23.cpp:171
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
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:76
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
Base abstract class for second order time dependent operators.
Definition: operator.hpp:451
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
double InitialRate(const Vector &x)
Definition: ex23.cpp:176
int dim
Definition: ex24.cpp:43
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1685
int open(const char hostname[], int port)
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:338
bool Good() const
Definition: optparser.hpp:125