MFEM  v4.5.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);
95  K->AddDomainIntegrator(new DiffusionIntegrator(*c2));
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);
103  M->AddDomainIntegrator(new MassIntegrator());
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);
202  args.AddOption(&mesh_file, "-m", "--mesh",
203  "Mesh file to use.");
204  args.AddOption(&ref_levels, "-r", "--refine",
205  "Number of times to refine the mesh uniformly.");
206  args.AddOption(&order, "-o", "--order",
207  "Order (degree) of the finite elements.");
208  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
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");
212  args.AddOption(&t_final, "-tf", "--t-final",
213  "Final time; start time is 0.");
214  args.AddOption(&dt, "-dt", "--time-step",
215  "Time step.");
216  args.AddOption(&speed, "-c", "--speed",
217  "Wave speed.");
218  args.AddOption(&dirichlet, "-dir", "--dirichlet", "-neu",
219  "--neumann",
220  "BC switch.");
221  args.AddOption(&ref_dir, "-r", "--ref",
222  "Reference directory for checking final solution.");
223  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
224  "--no-visualization",
225  "Enable or disable GLVis visualization.");
226  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
227  "--no-visit-datafiles",
228  "Save data files for VisIt (visit.llnl.gov) visualization.");
229  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
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 
292  // 6. Set the initial conditions for u. All boundaries are considered
293  // natural.
295  u_gf.ProjectCoefficient(u_0);
296  Vector u;
297  u_gf.GetTrueDofs(u);
298 
300  dudt_gf.ProjectCoefficient(dudt_0);
301  Vector dudt;
302  dudt_gf.GetTrueDofs(dudt);
303 
304  // 7. Initialize the conduction operator and the visualization.
305  Array<int> ess_bdr;
306  if (mesh->bdr_attributes.Size())
307  {
308  ess_bdr.SetSize(mesh->bdr_attributes.Max());
309 
310  if (dirichlet)
311  {
312  ess_bdr = 1;
313  }
314  else
315  {
316  ess_bdr = 0;
317  }
318  }
319 
320  WaveOperator oper(fespace, ess_bdr, speed);
321 
322  u_gf.SetFromTrueDofs(u);
323  {
324  ofstream omesh("ex23.mesh");
325  omesh.precision(precision);
326  mesh->Print(omesh);
327  ofstream osol("ex23-init.gf");
328  osol.precision(precision);
329  u_gf.Save(osol);
330  dudt_gf.Save(osol);
331  }
332 
333  VisItDataCollection visit_dc("Example23", mesh);
334  visit_dc.RegisterField("solution", &u_gf);
335  visit_dc.RegisterField("rate", &dudt_gf);
336  if (visit)
337  {
338  visit_dc.SetCycle(0);
339  visit_dc.SetTime(0.0);
340  visit_dc.Save();
341  }
342 
343  socketstream sout;
344  if (visualization)
345  {
346  char vishost[] = "localhost";
347  int visport = 19916;
348  sout.open(vishost, visport);
349  if (!sout)
350  {
351  cout << "Unable to connect to GLVis server at "
352  << vishost << ':' << visport << endl;
353  visualization = false;
354  cout << "GLVis visualization disabled.\n";
355  }
356  else
357  {
358  sout.precision(precision);
359  sout << "solution\n" << *mesh << dudt_gf;
360  sout << "pause\n";
361  sout << flush;
362  cout << "GLVis visualization paused."
363  << " Press space (in the GLVis window) to resume it.\n";
364  }
365  }
366 
367  // 8. Perform time-integration (looping over the time iterations, ti, with a
368  // time-step dt).
369  ode_solver->Init(oper);
370  double t = 0.0;
371 
372  bool last_step = false;
373  for (int ti = 1; !last_step; ti++)
374  {
375 
376  if (t + dt >= t_final - dt/2)
377  {
378  last_step = true;
379  }
380 
381  ode_solver->Step(u, dudt, t, dt);
382 
383  if (last_step || (ti % vis_steps) == 0)
384  {
385  cout << "step " << ti << ", t = " << t << endl;
386 
387  u_gf.SetFromTrueDofs(u);
388  dudt_gf.SetFromTrueDofs(dudt);
389  if (visualization)
390  {
391  sout << "solution\n" << *mesh << 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 
404  // 9. Save the final solution. This output can be viewed later using GLVis:
405  // "glvis -m ex23.mesh -g ex23-final.gf".
406  {
407  ofstream osol("ex23-final.gf");
408  osol.precision(precision);
409  u_gf.Save(osol);
410  dudt_gf.Save(osol);
411  }
412 
413  // 10. Free the used memory.
414  delete ode_solver;
415  delete mesh;
416 
417  return 0;
418 }
Conjugate gradient method.
Definition: solvers.hpp:465
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
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
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void Save()
Save the collection and a VisIt root file.
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
STL namespace.
Abstract class for solving systems of ODEs: d2x/dt2 = f(x,dx/dt,t)
Definition: ode.hpp:607
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2268
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition: gridfunc.cpp:367
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]...
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
The classical midpoint method.
Definition: ode.hpp:786
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9802
constexpr int visport
Data collection with VisIt I/O routines.
int main(int argc, char *argv[])
Definition: ex23.cpp:182
virtual void Init(SecondOrderTimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:989
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:235
void SetTime(double t)
Set physical time (for time-dependent simulations)
double InitialSolution(const Vector &x)
Definition: ex23.cpp:171
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
Base abstract class for second order time dependent operators.
Definition: operator.hpp:607
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
double InitialRate(const Vector &x)
Definition: ex23.cpp:176
int dim
Definition: ex24.cpp:53
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2399
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;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
RefCoord t[3]
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1671
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3676
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:382
double f(const Vector &p)