MFEM  v3.3
Finite element discretization library
ex16.cpp
Go to the documentation of this file.
1 // MFEM Example 16
2 //
3 // Compile with: make ex16
4 //
5 // Sample runs: ex16
6 // ex16 -m ../data/inline-tri.mesh
7 // ex16 -m ../data/disc-nurbs.mesh -tf 2
8 // ex16 -s 1 -a 0.0 -k 1.0
9 // ex16 -s 2 -a 1.0 -k 0.0
10 // ex16 -s 3 -a 0.5 -k 0.5 -o 4
11 // ex16 -s 14 -dt 1.0e-4 -tf 4.0e-2 -vs 40
12 // ex16 -m ../data/fichera-q2.mesh
13 // ex16 -m ../data/escher.mesh
14 // ex16 -m ../data/beam-tet.mesh -tf 10 -dt 0.1
15 // ex16 -m ../data/amr-quad.mesh -o 4 -r 0
16 // ex16 -m ../data/amr-hex.mesh -o 2 -r 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  FiniteElementSpace &fespace;
52  Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
53 
54  BilinearForm *M;
55  BilinearForm *K;
56 
57  SparseMatrix Mmat, Kmat;
58  SparseMatrix *T; // T = M + dt K
59  double current_dt;
60 
61  CGSolver M_solver; // Krylov solver for inverting the mass matrix M
62  DSmoother M_prec; // Preconditioner for the mass matrix M
63 
64  CGSolver T_solver; // Implicit solver for T = M + dt K
65  DSmoother T_prec; // Preconditioner for the implicit solver
66 
67  double alpha, kappa;
68 
69  mutable Vector z; // auxiliary vector
70 
71 public:
72  ConductionOperator(FiniteElementSpace &f, double alpha, double kappa,
73  const Vector &u);
74 
75  virtual void Mult(const Vector &u, Vector &du_dt) const;
76  /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
77  This is the only requirement for high-order SDIRK implicit integration.*/
78  virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
79 
80  /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
81  void SetParameters(const Vector &u);
82 
83  virtual ~ConductionOperator();
84 };
85 
86 double InitialTemperature(const Vector &x);
87 
88 int main(int argc, char *argv[])
89 {
90  // 1. Parse command-line options.
91  const char *mesh_file = "../data/star.mesh";
92  int ref_levels = 2;
93  int order = 2;
94  int ode_solver_type = 3;
95  double t_final = 0.5;
96  double dt = 1.0e-2;
97  double alpha = 1.0e-2;
98  double kappa = 0.5;
99  bool visualization = true;
100  bool visit = false;
101  int vis_steps = 5;
102 
103  int precision = 8;
104  cout.precision(precision);
105 
106  OptionsParser args(argc, argv);
107  args.AddOption(&mesh_file, "-m", "--mesh",
108  "Mesh file to use.");
109  args.AddOption(&ref_levels, "-r", "--refine",
110  "Number of times to refine the mesh uniformly.");
111  args.AddOption(&order, "-o", "--order",
112  "Order (degree) of the finite elements.");
113  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
114  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
115  "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
116  args.AddOption(&t_final, "-tf", "--t-final",
117  "Final time; start time is 0.");
118  args.AddOption(&dt, "-dt", "--time-step",
119  "Time step.");
120  args.AddOption(&alpha, "-a", "--alpha",
121  "Alpha coefficient.");
122  args.AddOption(&kappa, "-k", "--kappa",
123  "Kappa coefficient offset.");
124  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
125  "--no-visualization",
126  "Enable or disable GLVis visualization.");
127  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
128  "--no-visit-datafiles",
129  "Save data files for VisIt (visit.llnl.gov) visualization.");
130  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
131  "Visualize every n-th timestep.");
132  args.Parse();
133  if (!args.Good())
134  {
135  args.PrintUsage(cout);
136  return 1;
137  }
138  args.PrintOptions(cout);
139 
140  // 2. Read the mesh from the given mesh file. We can handle triangular,
141  // quadrilateral, tetrahedral and hexahedral meshes with the same code.
142  Mesh *mesh = new Mesh(mesh_file, 1, 1);
143  int dim = mesh->Dimension();
144 
145  // 3. Define the ODE solver used for time integration. Several implicit
146  // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
147  // explicit Runge-Kutta methods are available.
148  ODESolver *ode_solver;
149  switch (ode_solver_type)
150  {
151  // Implicit L-stable methods
152  case 1: ode_solver = new BackwardEulerSolver; break;
153  case 2: ode_solver = new SDIRK23Solver(2); break;
154  case 3: ode_solver = new SDIRK33Solver; break;
155  // Explicit methods
156  case 11: ode_solver = new ForwardEulerSolver; break;
157  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
158  case 13: ode_solver = new RK3SSPSolver; break;
159  case 14: ode_solver = new RK4Solver; break;
160  // Implicit A-stable methods (not L-stable)
161  case 22: ode_solver = new ImplicitMidpointSolver; break;
162  case 23: ode_solver = new SDIRK23Solver; break;
163  case 24: ode_solver = new SDIRK34Solver; break;
164  default:
165  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
166  return 3;
167  }
168 
169  // 4. Refine the mesh to increase the resolution. In this example we do
170  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
171  // command-line parameter.
172  for (int lev = 0; lev < ref_levels; lev++)
173  {
174  mesh->UniformRefinement();
175  }
176 
177  // 5. Define the vector finite element space representing the current and the
178  // initial temperature, u_ref.
179  H1_FECollection fe_coll(order, dim);
180  FiniteElementSpace fespace(mesh, &fe_coll);
181 
182  int fe_size = fespace.GetTrueVSize();
183  cout << "Number of temperature unknowns: " << fe_size << endl;
184 
185  GridFunction u_gf(&fespace);
186 
187  // 6. Set the initial conditions for u. All boundaries are considered
188  // natural.
190  u_gf.ProjectCoefficient(u_0);
191  Vector u;
192  u_gf.GetTrueDofs(u);
193 
194  // 7. Initialize the conduction operator and the visualization.
195  ConductionOperator oper(fespace, alpha, kappa, u);
196 
197  u_gf.SetFromTrueDofs(u);
198  {
199  ofstream omesh("ex16.mesh");
200  omesh.precision(precision);
201  mesh->Print(omesh);
202  ofstream osol("ex16-init.gf");
203  osol.precision(precision);
204  u_gf.Save(osol);
205  }
206 
207  VisItDataCollection visit_dc("Example16", mesh);
208  visit_dc.RegisterField("temperature", &u_gf);
209  if (visit)
210  {
211  visit_dc.SetCycle(0);
212  visit_dc.SetTime(0.0);
213  visit_dc.Save();
214  }
215 
216  socketstream sout;
217  if (visualization)
218  {
219  char vishost[] = "localhost";
220  int visport = 19916;
221  sout.open(vishost, visport);
222  if (!sout)
223  {
224  cout << "Unable to connect to GLVis server at "
225  << vishost << ':' << visport << endl;
226  visualization = false;
227  cout << "GLVis visualization disabled.\n";
228  }
229  else
230  {
231  sout.precision(precision);
232  sout << "solution\n" << *mesh << u_gf;
233  sout << "pause\n";
234  sout << flush;
235  cout << "GLVis visualization paused."
236  << " Press space (in the GLVis window) to resume it.\n";
237  }
238  }
239 
240  // 8. Perform time-integration (looping over the time iterations, ti, with a
241  // time-step dt).
242  ode_solver->Init(oper);
243  double t = 0.0;
244 
245  bool last_step = false;
246  for (int ti = 1; !last_step; ti++)
247  {
248  if (t + dt >= t_final - dt/2)
249  {
250  last_step = true;
251  }
252 
253  ode_solver->Step(u, t, dt);
254 
255  if (last_step || (ti % vis_steps) == 0)
256  {
257  cout << "step " << ti << ", t = " << t << endl;
258 
259  u_gf.SetFromTrueDofs(u);
260  if (visualization)
261  {
262  sout << "solution\n" << *mesh << u_gf << flush;
263  }
264 
265  if (visit)
266  {
267  visit_dc.SetCycle(ti);
268  visit_dc.SetTime(t);
269  visit_dc.Save();
270  }
271  }
272  oper.SetParameters(u);
273  }
274 
275  // 9. Save the final solution. This output can be viewed later using GLVis:
276  // "glvis -m ex16.mesh -g ex16-final.gf".
277  {
278  ofstream osol("ex16-final.gf");
279  osol.precision(precision);
280  u_gf.Save(osol);
281  }
282 
283  // 10. Free the used memory.
284  delete ode_solver;
285  delete mesh;
286 
287  return 0;
288 }
289 
290 ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
291  double kap, const Vector &u)
292  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
293  T(NULL), current_dt(0.0), z(height)
294 {
295  const double rel_tol = 1e-8;
296 
297  M = new BilinearForm(&fespace);
298  M->AddDomainIntegrator(new MassIntegrator());
299  M->Assemble();
300  M->FormSystemMatrix(ess_tdof_list, Mmat);
301 
302  M_solver.iterative_mode = false;
303  M_solver.SetRelTol(rel_tol);
304  M_solver.SetAbsTol(0.0);
305  M_solver.SetMaxIter(30);
306  M_solver.SetPrintLevel(0);
307  M_solver.SetPreconditioner(M_prec);
308  M_solver.SetOperator(Mmat);
309 
310  alpha = al;
311  kappa = kap;
312 
313  T_solver.iterative_mode = false;
314  T_solver.SetRelTol(rel_tol);
315  T_solver.SetAbsTol(0.0);
316  T_solver.SetMaxIter(100);
317  T_solver.SetPrintLevel(0);
318  T_solver.SetPreconditioner(T_prec);
319 
320  SetParameters(u);
321 }
322 
323 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
324 {
325  // Compute:
326  // du_dt = M^{-1}*-K(u)
327  // for du_dt
328  Kmat.Mult(u, z);
329  z.Neg(); // z = -z
330  M_solver.Mult(z, du_dt);
331 }
332 
333 void ConductionOperator::ImplicitSolve(const double dt,
334  const Vector &u, Vector &du_dt)
335 {
336  // Solve the equation:
337  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
338  // for du_dt
339  if (!T)
340  {
341  T = Add(1.0, Mmat, dt, Kmat);
342  current_dt = dt;
343  T_solver.SetOperator(*T);
344  }
345  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
346  Kmat.Mult(u, z);
347  z.Neg();
348  T_solver.Mult(z, du_dt);
349 }
350 
351 void ConductionOperator::SetParameters(const Vector &u)
352 {
353  GridFunction u_alpha_gf(&fespace);
354  u_alpha_gf.SetFromTrueDofs(u);
355  for (int i = 0; i < u_alpha_gf.Size(); i++)
356  {
357  u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
358  }
359 
360  delete K;
361  K = new BilinearForm(&fespace);
362 
363  GridFunctionCoefficient u_coeff(&u_alpha_gf);
364 
365  K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
366  K->Assemble();
367  K->FormSystemMatrix(ess_tdof_list, Kmat);
368  delete T;
369  T = NULL; // re-compute T on the next ImplicitSolve
370 }
371 
372 ConductionOperator::~ConductionOperator()
373 {
374  delete T;
375  delete M;
376  delete K;
377 }
378 
379 double InitialTemperature(const Vector &x)
380 {
381  if (x.Norml2() < 0.5)
382  {
383  return 2.0;
384  }
385  else
386  {
387  return 1.0;
388  }
389 }
Conjugate gradient method.
Definition: solvers.hpp:111
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
double InitialTemperature(const Vector &x)
Definition: ex16.cpp:379
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
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:140
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:667
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.
STL namespace.
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:212
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2185
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2807
int dim
Definition: ex3.cpp:47
Data type sparse matrix.
Definition: sparsemat.hpp:38
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
Data collection with VisIt I/O routines.
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:280
int Dimension() const
Definition: mesh.hpp:611
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
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 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
virtual int GetTrueVSize()
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:166
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1232
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:36
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
The classical forward Euler method.
Definition: ode.hpp:101
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:295
int main(int argc, char *argv[])
Definition: ex16.cpp:88
virtual void Print(std::ostream &out=std::cout) const
Definition: mesh.hpp:988
bool Good() const
Definition: optparser.hpp:120