MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
39using namespace std;
40using 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 */
52class ConductionOperator : public TimeDependentOperator
53{
54protected:
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 real_t 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 real_t alpha, kappa;
72
73 mutable Vector z; // auxiliary vector
74
75public:
76 ConductionOperator(FiniteElementSpace &f, real_t alpha, real_t 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 real_t 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
91
92int 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 real_t t_final = 0.5;
100 real_t dt = 1.0e-2;
101 real_t alpha = 1.0e-2;
102 real_t 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 real_t 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
296ConductionOperator::ConductionOperator(FiniteElementSpace &f, real_t al,
297 real_t kap, const Vector &u)
298 : TimeDependentOperator(f.GetTrueVSize(), (real_t) 0.0), fespace(f),
299 M(NULL), K(NULL), T(NULL), current_dt(0.0), z(height)
300{
301 const real_t 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
329void 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
339void ConductionOperator::ImplicitSolve(const real_t 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
357void 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
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
378ConductionOperator::~ConductionOperator()
379{
380 delete T;
381 delete M;
382 delete K;
383}
384
386{
387 if (x.Norml2() < 0.5)
388 {
389 return 2.0;
390 }
391 else
392 {
393 return 1.0;
394 }
395}
Backward Euler ODE solver. L-stable.
Definition ode.hpp:412
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
Data type for scaled Jacobi-type smoother of sparse matrix.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
The classical forward Euler method.
Definition ode.hpp:118
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition gridfunc.cpp:381
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition gridfunc.cpp:366
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Implicit midpoint method. A-stable, not L-stable.
Definition ode.hpp:425
Mesh data type.
Definition mesh.hpp:56
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition ode.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &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 ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition ode.hpp:151
The classical explicit forth-order Runge-Kutta method, RK4.
Definition ode.hpp:164
Data type sparse matrix.
Definition sparsemat.hpp:51
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Base abstract class for first order time dependent operators.
Definition operator.hpp:317
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
const real_t alpha
Definition ex15.cpp:369
real_t InitialTemperature(const Vector &x)
Definition ex16.cpp:385
real_t kappa
Definition ex24.cpp:54
int dim
Definition ex24.cpp:53
int main()
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
const char vishost[]
RefCoord t[3]