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