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