1 /* ---------------------------------------------------------------------
2 *
3 * Copyright (C) 2009 - 2020 by the deal.II authors
4 *
5 * This file is part of the deal.II library.
6 *
7 * The deal.II library is free software; you can use it, redistribute
8 * it, and/or modify it under the terms of the GNU Lesser General
9 * Public License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 * The full text of the license can be found in the file LICENSE.md at
12 * the top level directory of deal.II.
13 *
14 * ---------------------------------------------------------------------
15
16 *
17 * Author: Abner Salgado, Texas A&M University 2009
18 */
19
20
21 // @sect3{Include files}
22
23 // We start by including all the necessary deal.II header files and some C++
24 // related ones. Each one of them has been discussed in previous tutorial
25 // programs, so we will not get into details here.
26 #include <deal.II/base/parameter_handler.h>
27 #include <deal.II/base/point.h>
28 #include <deal.II/base/function.h>
29 #include <deal.II/base/quadrature_lib.h>
30 #include <deal.II/base/multithread_info.h>
31 #include <deal.II/base/thread_management.h>
32 #include <deal.II/base/work_stream.h>
33 #include <deal.II/base/parallel.h>
34 #include <deal.II/base/utilities.h>
35 #include <deal.II/base/conditional_ostream.h>
36
37 #include <deal.II/lac/vector.h>
38 #include <deal.II/lac/sparse_matrix.h>
39 #include <deal.II/lac/dynamic_sparsity_pattern.h>
40 #include <deal.II/lac/solver_cg.h>
41 #include <deal.II/lac/precondition.h>
42 #include <deal.II/lac/solver_gmres.h>
43 #include <deal.II/lac/sparse_ilu.h>
44 #include <deal.II/lac/sparse_direct.h>
45 #include <deal.II/lac/affine_constraints.h>
46
47 #include <deal.II/grid/tria.h>
48 #include <deal.II/grid/grid_generator.h>
49 #include <deal.II/grid/grid_refinement.h>
50 #include <deal.II/grid/tria_accessor.h>
51 #include <deal.II/grid/tria_iterator.h>
52 #include <deal.II/grid/grid_in.h>
53
54 #include <deal.II/dofs/dof_handler.h>
55 #include <deal.II/dofs/dof_accessor.h>
56 #include <deal.II/dofs/dof_tools.h>
57 #include <deal.II/dofs/dof_renumbering.h>
58
59 #include <deal.II/fe/fe_q.h>
60 #include <deal.II/fe/fe_values.h>
61 #include <deal.II/fe/fe_tools.h>
62 #include <deal.II/fe/fe_system.h>
63
64 #include <deal.II/numerics/matrix_tools.h>
65 #include <deal.II/numerics/vector_tools.h>
66 #include <deal.II/numerics/data_out.h>
67
68 #include <fstream>
69 #include <cmath>
70 #include <iostream>
71
72 // Finally this is as in all previous programs:
73 namespace Step35
74 {
75 using namespace dealii;
76
77 // @sect3{Run time parameters}
78 //
79 // Since our method has several parameters that can be fine-tuned we put them
80 // into an external file, so that they can be determined at run-time.
81 //
82 // This includes, in particular, the formulation of the equation for the
83 // auxiliary variable $\phi$, for which we declare an <code>enum</code>. Next,
84 // we declare a class that is going to read and store all the parameters that
85 // our program needs to run.
86 namespace RunTimeParameters
87 {
88 enum class Method
89 {
90 standard,
91 rotational
92 };
93
94 class Data_Storage
95 {
96 public:
97 Data_Storage();
98
99 void read_data(const std::string &filename);
100
101 Method form;
102
103 double dt;
104 double initial_time;
105 double final_time;
106
107 double Reynolds;
108
109 unsigned int n_global_refines;
110
111 unsigned int pressure_degree;
112
113 unsigned int vel_max_iterations;
114 unsigned int vel_Krylov_size;
115 unsigned int vel_off_diagonals;
116 unsigned int vel_update_prec;
117 double vel_eps;
118 double vel_diag_strength;
119
120 bool verbose;
121 unsigned int output_interval;
122
123 protected:
124 ParameterHandler prm;
125 };
126
127 // In the constructor of this class we declare all the parameters. The
128 // details of how this works have been discussed elsewhere, for example in
129 // step-29.
Data_Storage()130 Data_Storage::Data_Storage()
131 : form(Method::rotational)
132 , dt(5e-4)
133 , initial_time(0.)
134 , final_time(1.)
135 , Reynolds(1.)
136 , n_global_refines(0)
137 , pressure_degree(1)
138 , vel_max_iterations(1000)
139 , vel_Krylov_size(30)
140 , vel_off_diagonals(60)
141 , vel_update_prec(15)
142 , vel_eps(1e-12)
143 , vel_diag_strength(0.01)
144 , verbose(true)
145 , output_interval(15)
146 {
147 prm.declare_entry("Method_Form",
148 "rotational",
149 Patterns::Selection("rotational|standard"),
150 " Used to select the type of method that we are going "
151 "to use. ");
152 prm.enter_subsection("Physical data");
153 {
154 prm.declare_entry("initial_time",
155 "0.",
156 Patterns::Double(0.),
157 " The initial time of the simulation. ");
158 prm.declare_entry("final_time",
159 "1.",
160 Patterns::Double(0.),
161 " The final time of the simulation. ");
162 prm.declare_entry("Reynolds",
163 "1.",
164 Patterns::Double(0.),
165 " The Reynolds number. ");
166 }
167 prm.leave_subsection();
168
169 prm.enter_subsection("Time step data");
170 {
171 prm.declare_entry("dt",
172 "5e-4",
173 Patterns::Double(0.),
174 " The time step size. ");
175 }
176 prm.leave_subsection();
177
178 prm.enter_subsection("Space discretization");
179 {
180 prm.declare_entry("n_of_refines",
181 "0",
182 Patterns::Integer(0, 15),
183 " The number of global refines we do on the mesh. ");
184 prm.declare_entry("pressure_fe_degree",
185 "1",
186 Patterns::Integer(1, 5),
187 " The polynomial degree for the pressure space. ");
188 }
189 prm.leave_subsection();
190
191 prm.enter_subsection("Data solve velocity");
192 {
193 prm.declare_entry(
194 "max_iterations",
195 "1000",
196 Patterns::Integer(1, 1000),
197 " The maximal number of iterations GMRES must make. ");
198 prm.declare_entry("eps",
199 "1e-12",
200 Patterns::Double(0.),
201 " The stopping criterion. ");
202 prm.declare_entry("Krylov_size",
203 "30",
204 Patterns::Integer(1),
205 " The size of the Krylov subspace to be used. ");
206 prm.declare_entry("off_diagonals",
207 "60",
208 Patterns::Integer(0),
209 " The number of off-diagonal elements ILU must "
210 "compute. ");
211 prm.declare_entry("diag_strength",
212 "0.01",
213 Patterns::Double(0.),
214 " Diagonal strengthening coefficient. ");
215 prm.declare_entry("update_prec",
216 "15",
217 Patterns::Integer(1),
218 " This number indicates how often we need to "
219 "update the preconditioner");
220 }
221 prm.leave_subsection();
222
223 prm.declare_entry("verbose",
224 "true",
225 Patterns::Bool(),
226 " This indicates whether the output of the solution "
227 "process should be verbose. ");
228
229 prm.declare_entry("output_interval",
230 "1",
231 Patterns::Integer(1),
232 " This indicates between how many time steps we print "
233 "the solution. ");
234 }
235
236
237
read_data(const std::string & filename)238 void Data_Storage::read_data(const std::string &filename)
239 {
240 std::ifstream file(filename);
241 AssertThrow(file, ExcFileNotOpen(filename));
242
243 prm.parse_input(file);
244
245 if (prm.get("Method_Form") == std::string("rotational"))
246 form = Method::rotational;
247 else
248 form = Method::standard;
249
250 prm.enter_subsection("Physical data");
251 {
252 initial_time = prm.get_double("initial_time");
253 final_time = prm.get_double("final_time");
254 Reynolds = prm.get_double("Reynolds");
255 }
256 prm.leave_subsection();
257
258 prm.enter_subsection("Time step data");
259 {
260 dt = prm.get_double("dt");
261 }
262 prm.leave_subsection();
263
264 prm.enter_subsection("Space discretization");
265 {
266 n_global_refines = prm.get_integer("n_of_refines");
267 pressure_degree = prm.get_integer("pressure_fe_degree");
268 }
269 prm.leave_subsection();
270
271 prm.enter_subsection("Data solve velocity");
272 {
273 vel_max_iterations = prm.get_integer("max_iterations");
274 vel_eps = prm.get_double("eps");
275 vel_Krylov_size = prm.get_integer("Krylov_size");
276 vel_off_diagonals = prm.get_integer("off_diagonals");
277 vel_diag_strength = prm.get_double("diag_strength");
278 vel_update_prec = prm.get_integer("update_prec");
279 }
280 prm.leave_subsection();
281
282 verbose = prm.get_bool("verbose");
283
284 output_interval = prm.get_integer("output_interval");
285 }
286 } // namespace RunTimeParameters
287
288
289
290 // @sect3{Equation data}
291
292 // In the next namespace, we declare the initial and boundary conditions:
293 namespace EquationData
294 {
295 // As we have chosen a completely decoupled formulation, we will not take
296 // advantage of deal.II's capabilities to handle vector valued problems. We
297 // do, however, want to use an interface for the equation data that is
298 // somehow dimension independent. To be able to do that, our functions
299 // should be able to know on which spatial component we are currently
300 // working, and we should be able to have a common interface to do that. The
301 // following class is an attempt in that direction.
302 template <int dim>
303 class MultiComponentFunction : public Function<dim>
304 {
305 public:
306 MultiComponentFunction(const double initial_time = 0.);
307 void set_component(const unsigned int d);
308
309 protected:
310 unsigned int comp;
311 };
312
313 template <int dim>
MultiComponentFunction(const double initial_time)314 MultiComponentFunction<dim>::MultiComponentFunction(
315 const double initial_time)
316 : Function<dim>(1, initial_time)
317 , comp(0)
318 {}
319
320
321 template <int dim>
set_component(const unsigned int d)322 void MultiComponentFunction<dim>::set_component(const unsigned int d)
323 {
324 Assert(d < dim, ExcIndexRange(d, 0, dim));
325 comp = d;
326 }
327
328
329 // With this class defined, we declare classes that describe the boundary
330 // conditions for velocity and pressure:
331 template <int dim>
332 class Velocity : public MultiComponentFunction<dim>
333 {
334 public:
335 Velocity(const double initial_time = 0.0);
336
337 virtual double value(const Point<dim> & p,
338 const unsigned int component = 0) const override;
339
340 virtual void value_list(const std::vector<Point<dim>> &points,
341 std::vector<double> & values,
342 const unsigned int component = 0) const override;
343 };
344
345
346 template <int dim>
Velocity(const double initial_time)347 Velocity<dim>::Velocity(const double initial_time)
348 : MultiComponentFunction<dim>(initial_time)
349 {}
350
351
352 template <int dim>
value_list(const std::vector<Point<dim>> & points,std::vector<double> & values,const unsigned int) const353 void Velocity<dim>::value_list(const std::vector<Point<dim>> &points,
354 std::vector<double> & values,
355 const unsigned int) const
356 {
357 const unsigned int n_points = points.size();
358 Assert(values.size() == n_points,
359 ExcDimensionMismatch(values.size(), n_points));
360 for (unsigned int i = 0; i < n_points; ++i)
361 values[i] = Velocity<dim>::value(points[i]);
362 }
363
364
365 template <int dim>
value(const Point<dim> & p,const unsigned int) const366 double Velocity<dim>::value(const Point<dim> &p, const unsigned int) const
367 {
368 if (this->comp == 0)
369 {
370 const double Um = 1.5;
371 const double H = 4.1;
372 return 4. * Um * p(1) * (H - p(1)) / (H * H);
373 }
374 else
375 return 0.;
376 }
377
378
379
380 template <int dim>
381 class Pressure : public Function<dim>
382 {
383 public:
384 Pressure(const double initial_time = 0.0);
385
386 virtual double value(const Point<dim> & p,
387 const unsigned int component = 0) const override;
388
389 virtual void value_list(const std::vector<Point<dim>> &points,
390 std::vector<double> & values,
391 const unsigned int component = 0) const override;
392 };
393
394 template <int dim>
Pressure(const double initial_time)395 Pressure<dim>::Pressure(const double initial_time)
396 : Function<dim>(1, initial_time)
397 {}
398
399
400 template <int dim>
value(const Point<dim> & p,const unsigned int component) const401 double Pressure<dim>::value(const Point<dim> & p,
402 const unsigned int component) const
403 {
404 (void)component;
405 AssertIndexRange(component, 1);
406 return 25. - p(0);
407 }
408
409 template <int dim>
value_list(const std::vector<Point<dim>> & points,std::vector<double> & values,const unsigned int component) const410 void Pressure<dim>::value_list(const std::vector<Point<dim>> &points,
411 std::vector<double> & values,
412 const unsigned int component) const
413 {
414 (void)component;
415 AssertIndexRange(component, 1);
416 const unsigned int n_points = points.size();
417 Assert(values.size() == n_points,
418 ExcDimensionMismatch(values.size(), n_points));
419 for (unsigned int i = 0; i < n_points; ++i)
420 values[i] = Pressure<dim>::value(points[i]);
421 }
422 } // namespace EquationData
423
424
425
426 // @sect3{The <code>NavierStokesProjection</code> class}
427
428 // Now for the main class of the program. It implements the various versions
429 // of the projection method for Navier-Stokes equations. The names for all the
430 // methods and member variables should be self-explanatory, taking into
431 // account the implementation details given in the introduction.
432 template <int dim>
433 class NavierStokesProjection
434 {
435 public:
436 NavierStokesProjection(const RunTimeParameters::Data_Storage &data);
437
438 void run(const bool verbose = false, const unsigned int n_plots = 10);
439
440 protected:
441 RunTimeParameters::Method type;
442
443 const unsigned int deg;
444 const double dt;
445 const double t_0;
446 const double T;
447 const double Re;
448
449 EquationData::Velocity<dim> vel_exact;
450 std::map<types::global_dof_index, double> boundary_values;
451 std::vector<types::boundary_id> boundary_ids;
452
453 Triangulation<dim> triangulation;
454
455 FE_Q<dim> fe_velocity;
456 FE_Q<dim> fe_pressure;
457
458 DoFHandler<dim> dof_handler_velocity;
459 DoFHandler<dim> dof_handler_pressure;
460
461 QGauss<dim> quadrature_pressure;
462 QGauss<dim> quadrature_velocity;
463
464 SparsityPattern sparsity_pattern_velocity;
465 SparsityPattern sparsity_pattern_pressure;
466 SparsityPattern sparsity_pattern_pres_vel;
467
468 SparseMatrix<double> vel_Laplace_plus_Mass;
469 SparseMatrix<double> vel_it_matrix[dim];
470 SparseMatrix<double> vel_Mass;
471 SparseMatrix<double> vel_Laplace;
472 SparseMatrix<double> vel_Advection;
473 SparseMatrix<double> pres_Laplace;
474 SparseMatrix<double> pres_Mass;
475 SparseMatrix<double> pres_Diff[dim];
476 SparseMatrix<double> pres_iterative;
477
478 Vector<double> pres_n;
479 Vector<double> pres_n_minus_1;
480 Vector<double> phi_n;
481 Vector<double> phi_n_minus_1;
482 Vector<double> u_n[dim];
483 Vector<double> u_n_minus_1[dim];
484 Vector<double> u_star[dim];
485 Vector<double> force[dim];
486 Vector<double> v_tmp;
487 Vector<double> pres_tmp;
488 Vector<double> rot_u;
489
490 SparseILU<double> prec_velocity[dim];
491 SparseILU<double> prec_pres_Laplace;
492 SparseDirectUMFPACK prec_mass;
493 SparseDirectUMFPACK prec_vel_mass;
494
495 DeclException2(ExcInvalidTimeStep,
496 double,
497 double,
498 << " The time step " << arg1 << " is out of range."
499 << std::endl
500 << " The permitted range is (0," << arg2 << "]");
501
502 void create_triangulation_and_dofs(const unsigned int n_refines);
503
504 void initialize();
505
506 void interpolate_velocity();
507
508 void diffusion_step(const bool reinit_prec);
509
510 void projection_step(const bool reinit_prec);
511
512 void update_pressure(const bool reinit_prec);
513
514 private:
515 unsigned int vel_max_its;
516 unsigned int vel_Krylov_size;
517 unsigned int vel_off_diagonals;
518 unsigned int vel_update_prec;
519 double vel_eps;
520 double vel_diag_strength;
521
522 void initialize_velocity_matrices();
523
524 void initialize_pressure_matrices();
525
526 // The next few structures and functions are for doing various things in
527 // parallel. They follow the scheme laid out in @ref threads, using the
528 // WorkStream class. As explained there, this requires us to declare two
529 // structures for each of the assemblers, a per-task data and a scratch data
530 // structure. These are then handed over to functions that assemble local
531 // contributions and that copy these local contributions to the global
532 // objects.
533 //
534 // One of the things that are specific to this program is that we don't just
535 // have a single DoFHandler object that represents both the velocities and
536 // the pressure, but we use individual DoFHandler objects for these two
537 // kinds of variables. We pay for this optimization when we want to assemble
538 // terms that involve both variables, such as the divergence of the velocity
539 // and the gradient of the pressure, times the respective test functions.
540 // When doing so, we can't just anymore use a single FEValues object, but
541 // rather we need two, and they need to be initialized with cell iterators
542 // that point to the same cell in the triangulation but different
543 // DoFHandlers.
544 //
545 // To do this in practice, we declare a "synchronous" iterator -- an object
546 // that internally consists of several (in our case two) iterators, and each
547 // time the synchronous iteration is moved forward one step, each of the
548 // iterators stored internally is moved forward one step as well, thereby
549 // always staying in sync. As it so happens, there is a deal.II class that
550 // facilitates this sort of thing. (What is important here is to know that
551 // two DoFHandler objects built on the same triangulation will walk over the
552 // cells of the triangulation in the same order.)
553 using IteratorTuple =
554 std::tuple<typename DoFHandler<dim>::active_cell_iterator,
555 typename DoFHandler<dim>::active_cell_iterator>;
556
557 using IteratorPair = SynchronousIterators<IteratorTuple>;
558
559 void initialize_gradient_operator();
560
561 struct InitGradPerTaskData
562 {
563 unsigned int d;
564 unsigned int vel_dpc;
565 unsigned int pres_dpc;
566 FullMatrix<double> local_grad;
567 std::vector<types::global_dof_index> vel_local_dof_indices;
568 std::vector<types::global_dof_index> pres_local_dof_indices;
569
InitGradPerTaskDataStep35::NavierStokesProjection::InitGradPerTaskData570 InitGradPerTaskData(const unsigned int dd,
571 const unsigned int vdpc,
572 const unsigned int pdpc)
573 : d(dd)
574 , vel_dpc(vdpc)
575 , pres_dpc(pdpc)
576 , local_grad(vdpc, pdpc)
577 , vel_local_dof_indices(vdpc)
578 , pres_local_dof_indices(pdpc)
579 {}
580 };
581
582 struct InitGradScratchData
583 {
584 unsigned int nqp;
585 FEValues<dim> fe_val_vel;
586 FEValues<dim> fe_val_pres;
InitGradScratchDataStep35::NavierStokesProjection::InitGradScratchData587 InitGradScratchData(const FE_Q<dim> & fe_v,
588 const FE_Q<dim> & fe_p,
589 const QGauss<dim> &quad,
590 const UpdateFlags flags_v,
591 const UpdateFlags flags_p)
592 : nqp(quad.size())
593 , fe_val_vel(fe_v, quad, flags_v)
594 , fe_val_pres(fe_p, quad, flags_p)
595 {}
InitGradScratchDataStep35::NavierStokesProjection::InitGradScratchData596 InitGradScratchData(const InitGradScratchData &data)
597 : nqp(data.nqp)
598 , fe_val_vel(data.fe_val_vel.get_fe(),
599 data.fe_val_vel.get_quadrature(),
600 data.fe_val_vel.get_update_flags())
601 , fe_val_pres(data.fe_val_pres.get_fe(),
602 data.fe_val_pres.get_quadrature(),
603 data.fe_val_pres.get_update_flags())
604 {}
605 };
606
607 void assemble_one_cell_of_gradient(const IteratorPair & SI,
608 InitGradScratchData &scratch,
609 InitGradPerTaskData &data);
610
611 void copy_gradient_local_to_global(const InitGradPerTaskData &data);
612
613 // The same general layout also applies to the following classes and
614 // functions implementing the assembly of the advection term:
615 void assemble_advection_term();
616
617 struct AdvectionPerTaskData
618 {
619 FullMatrix<double> local_advection;
620 std::vector<types::global_dof_index> local_dof_indices;
AdvectionPerTaskDataStep35::NavierStokesProjection::AdvectionPerTaskData621 AdvectionPerTaskData(const unsigned int dpc)
622 : local_advection(dpc, dpc)
623 , local_dof_indices(dpc)
624 {}
625 };
626
627 struct AdvectionScratchData
628 {
629 unsigned int nqp;
630 unsigned int dpc;
631 std::vector<Point<dim>> u_star_local;
632 std::vector<Tensor<1, dim>> grad_u_star;
633 std::vector<double> u_star_tmp;
634 FEValues<dim> fe_val;
AdvectionScratchDataStep35::NavierStokesProjection::AdvectionScratchData635 AdvectionScratchData(const FE_Q<dim> & fe,
636 const QGauss<dim> &quad,
637 const UpdateFlags flags)
638 : nqp(quad.size())
639 , dpc(fe.n_dofs_per_cell())
640 , u_star_local(nqp)
641 , grad_u_star(nqp)
642 , u_star_tmp(nqp)
643 , fe_val(fe, quad, flags)
644 {}
645
AdvectionScratchDataStep35::NavierStokesProjection::AdvectionScratchData646 AdvectionScratchData(const AdvectionScratchData &data)
647 : nqp(data.nqp)
648 , dpc(data.dpc)
649 , u_star_local(nqp)
650 , grad_u_star(nqp)
651 , u_star_tmp(nqp)
652 , fe_val(data.fe_val.get_fe(),
653 data.fe_val.get_quadrature(),
654 data.fe_val.get_update_flags())
655 {}
656 };
657
658 void assemble_one_cell_of_advection(
659 const typename DoFHandler<dim>::active_cell_iterator &cell,
660 AdvectionScratchData & scratch,
661 AdvectionPerTaskData & data);
662
663 void copy_advection_local_to_global(const AdvectionPerTaskData &data);
664
665 // The final few functions implement the diffusion solve as well as
666 // postprocessing the output, including computing the curl of the velocity:
667 void diffusion_component_solve(const unsigned int d);
668
669 void output_results(const unsigned int step);
670
671 void assemble_vorticity(const bool reinit_prec);
672 };
673
674
675
676 // @sect4{ <code>NavierStokesProjection::NavierStokesProjection</code> }
677
678 // In the constructor, we just read all the data from the
679 // <code>Data_Storage</code> object that is passed as an argument, verify that
680 // the data we read is reasonable and, finally, create the triangulation and
681 // load the initial data.
682 template <int dim>
NavierStokesProjection(const RunTimeParameters::Data_Storage & data)683 NavierStokesProjection<dim>::NavierStokesProjection(
684 const RunTimeParameters::Data_Storage &data)
685 : type(data.form)
686 , deg(data.pressure_degree)
687 , dt(data.dt)
688 , t_0(data.initial_time)
689 , T(data.final_time)
690 , Re(data.Reynolds)
691 , vel_exact(data.initial_time)
692 , fe_velocity(deg + 1)
693 , fe_pressure(deg)
694 , dof_handler_velocity(triangulation)
695 , dof_handler_pressure(triangulation)
696 , quadrature_pressure(deg + 1)
697 , quadrature_velocity(deg + 2)
698 , vel_max_its(data.vel_max_iterations)
699 , vel_Krylov_size(data.vel_Krylov_size)
700 , vel_off_diagonals(data.vel_off_diagonals)
701 , vel_update_prec(data.vel_update_prec)
702 , vel_eps(data.vel_eps)
703 , vel_diag_strength(data.vel_diag_strength)
704 {
705 if (deg < 1)
706 std::cout
707 << " WARNING: The chosen pair of finite element spaces is not stable."
708 << std::endl
709 << " The obtained results will be nonsense" << std::endl;
710
711 AssertThrow(!((dt <= 0.) || (dt > .5 * T)), ExcInvalidTimeStep(dt, .5 * T));
712
713 create_triangulation_and_dofs(data.n_global_refines);
714 initialize();
715 }
716
717
718 // @sect4{<code>NavierStokesProjection::create_triangulation_and_dofs</code>}
719
720 // The method that creates the triangulation and refines it the needed number
721 // of times. After creating the triangulation, it creates the mesh dependent
722 // data, i.e. it distributes degrees of freedom and renumbers them, and
723 // initializes the matrices and vectors that we will use.
724 template <int dim>
create_triangulation_and_dofs(const unsigned int n_refines)725 void NavierStokesProjection<dim>::create_triangulation_and_dofs(
726 const unsigned int n_refines)
727 {
728 GridIn<dim> grid_in;
729 grid_in.attach_triangulation(triangulation);
730
731 {
732 std::string filename = "nsbench2.inp";
733 std::ifstream file(filename);
734 Assert(file, ExcFileNotOpen(filename.c_str()));
735 grid_in.read_ucd(file);
736 }
737
738 std::cout << "Number of refines = " << n_refines << std::endl;
739 triangulation.refine_global(n_refines);
740 std::cout << "Number of active cells: " << triangulation.n_active_cells()
741 << std::endl;
742
743 boundary_ids = triangulation.get_boundary_ids();
744
745 dof_handler_velocity.distribute_dofs(fe_velocity);
746 DoFRenumbering::boost::Cuthill_McKee(dof_handler_velocity);
747 dof_handler_pressure.distribute_dofs(fe_pressure);
748 DoFRenumbering::boost::Cuthill_McKee(dof_handler_pressure);
749
750 initialize_velocity_matrices();
751 initialize_pressure_matrices();
752 initialize_gradient_operator();
753
754 pres_n.reinit(dof_handler_pressure.n_dofs());
755 pres_n_minus_1.reinit(dof_handler_pressure.n_dofs());
756 phi_n.reinit(dof_handler_pressure.n_dofs());
757 phi_n_minus_1.reinit(dof_handler_pressure.n_dofs());
758 pres_tmp.reinit(dof_handler_pressure.n_dofs());
759 for (unsigned int d = 0; d < dim; ++d)
760 {
761 u_n[d].reinit(dof_handler_velocity.n_dofs());
762 u_n_minus_1[d].reinit(dof_handler_velocity.n_dofs());
763 u_star[d].reinit(dof_handler_velocity.n_dofs());
764 force[d].reinit(dof_handler_velocity.n_dofs());
765 }
766 v_tmp.reinit(dof_handler_velocity.n_dofs());
767 rot_u.reinit(dof_handler_velocity.n_dofs());
768
769 std::cout << "dim (X_h) = " << (dof_handler_velocity.n_dofs() * dim) //
770 << std::endl //
771 << "dim (M_h) = " << dof_handler_pressure.n_dofs() //
772 << std::endl //
773 << "Re = " << Re << std::endl //
774 << std::endl;
775 }
776
777
778 // @sect4{ <code>NavierStokesProjection::initialize</code> }
779
780 // This method creates the constant matrices and loads the initial data
781 template <int dim>
initialize()782 void NavierStokesProjection<dim>::initialize()
783 {
784 vel_Laplace_plus_Mass = 0.;
785 vel_Laplace_plus_Mass.add(1. / Re, vel_Laplace);
786 vel_Laplace_plus_Mass.add(1.5 / dt, vel_Mass);
787
788 EquationData::Pressure<dim> pres(t_0);
789 VectorTools::interpolate(dof_handler_pressure, pres, pres_n_minus_1);
790 pres.advance_time(dt);
791 VectorTools::interpolate(dof_handler_pressure, pres, pres_n);
792 phi_n = 0.;
793 phi_n_minus_1 = 0.;
794 for (unsigned int d = 0; d < dim; ++d)
795 {
796 vel_exact.set_time(t_0);
797 vel_exact.set_component(d);
798 VectorTools::interpolate(dof_handler_velocity,
799 vel_exact,
800 u_n_minus_1[d]);
801 vel_exact.advance_time(dt);
802 VectorTools::interpolate(dof_handler_velocity, vel_exact, u_n[d]);
803 }
804 }
805
806
807 // @sect4{ <code>NavierStokesProjection::initialize_*_matrices</code> }
808
809 // In this set of methods we initialize the sparsity patterns, the constraints
810 // (if any) and assemble the matrices that do not depend on the timestep
811 // <code>dt</code>. Note that for the Laplace and mass matrices, we can use
812 // functions in the library that do this. Because the expensive operations of
813 // this function -- creating the two matrices -- are entirely independent, we
814 // could in principle mark them as tasks that can be worked on in %parallel
815 // using the Threads::new_task functions. We won't do that here since these
816 // functions internally already are parallelized, and in particular because
817 // the current function is only called once per program run and so does not
818 // incur a cost in each time step. The necessary modifications would be quite
819 // straightforward, however.
820 template <int dim>
initialize_velocity_matrices()821 void NavierStokesProjection<dim>::initialize_velocity_matrices()
822 {
823 {
824 DynamicSparsityPattern dsp(dof_handler_velocity.n_dofs(),
825 dof_handler_velocity.n_dofs());
826 DoFTools::make_sparsity_pattern(dof_handler_velocity, dsp);
827 sparsity_pattern_velocity.copy_from(dsp);
828 }
829 vel_Laplace_plus_Mass.reinit(sparsity_pattern_velocity);
830 for (unsigned int d = 0; d < dim; ++d)
831 vel_it_matrix[d].reinit(sparsity_pattern_velocity);
832 vel_Mass.reinit(sparsity_pattern_velocity);
833 vel_Laplace.reinit(sparsity_pattern_velocity);
834 vel_Advection.reinit(sparsity_pattern_velocity);
835
836 MatrixCreator::create_mass_matrix(dof_handler_velocity,
837 quadrature_velocity,
838 vel_Mass);
839 MatrixCreator::create_laplace_matrix(dof_handler_velocity,
840 quadrature_velocity,
841 vel_Laplace);
842 }
843
844 // The initialization of the matrices that act on the pressure space is
845 // similar to the ones that act on the velocity space.
846 template <int dim>
initialize_pressure_matrices()847 void NavierStokesProjection<dim>::initialize_pressure_matrices()
848 {
849 {
850 DynamicSparsityPattern dsp(dof_handler_pressure.n_dofs(),
851 dof_handler_pressure.n_dofs());
852 DoFTools::make_sparsity_pattern(dof_handler_pressure, dsp);
853 sparsity_pattern_pressure.copy_from(dsp);
854 }
855
856 pres_Laplace.reinit(sparsity_pattern_pressure);
857 pres_iterative.reinit(sparsity_pattern_pressure);
858 pres_Mass.reinit(sparsity_pattern_pressure);
859
860 MatrixCreator::create_laplace_matrix(dof_handler_pressure,
861 quadrature_pressure,
862 pres_Laplace);
863 MatrixCreator::create_mass_matrix(dof_handler_pressure,
864 quadrature_pressure,
865 pres_Mass);
866 }
867
868
869 // For the gradient operator, we start by initializing the sparsity pattern
870 // and compressing it. It is important to notice here that the gradient
871 // operator acts from the pressure space into the velocity space, so we have
872 // to deal with two different finite element spaces. To keep the loops
873 // synchronized, we use the alias that we have defined before, namely
874 // <code>PairedIterators</code> and <code>IteratorPair</code>.
875 template <int dim>
initialize_gradient_operator()876 void NavierStokesProjection<dim>::initialize_gradient_operator()
877 {
878 {
879 DynamicSparsityPattern dsp(dof_handler_velocity.n_dofs(),
880 dof_handler_pressure.n_dofs());
881 DoFTools::make_sparsity_pattern(dof_handler_velocity,
882 dof_handler_pressure,
883 dsp);
884 sparsity_pattern_pres_vel.copy_from(dsp);
885 }
886
887 InitGradPerTaskData per_task_data(0,
888 fe_velocity.n_dofs_per_cell(),
889 fe_pressure.n_dofs_per_cell());
890 InitGradScratchData scratch_data(fe_velocity,
891 fe_pressure,
892 quadrature_velocity,
893 update_gradients | update_JxW_values,
894 update_values);
895
896 for (unsigned int d = 0; d < dim; ++d)
897 {
898 pres_Diff[d].reinit(sparsity_pattern_pres_vel);
899 per_task_data.d = d;
900 WorkStream::run(
901 IteratorPair(IteratorTuple(dof_handler_velocity.begin_active(),
902 dof_handler_pressure.begin_active())),
903 IteratorPair(IteratorTuple(dof_handler_velocity.end(),
904 dof_handler_pressure.end())),
905 *this,
906 &NavierStokesProjection<dim>::assemble_one_cell_of_gradient,
907 &NavierStokesProjection<dim>::copy_gradient_local_to_global,
908 scratch_data,
909 per_task_data);
910 }
911 }
912
913 template <int dim>
assemble_one_cell_of_gradient(const IteratorPair & SI,InitGradScratchData & scratch,InitGradPerTaskData & data)914 void NavierStokesProjection<dim>::assemble_one_cell_of_gradient(
915 const IteratorPair & SI,
916 InitGradScratchData &scratch,
917 InitGradPerTaskData &data)
918 {
919 scratch.fe_val_vel.reinit(std::get<0>(*SI));
920 scratch.fe_val_pres.reinit(std::get<1>(*SI));
921
922 std::get<0>(*SI)->get_dof_indices(data.vel_local_dof_indices);
923 std::get<1>(*SI)->get_dof_indices(data.pres_local_dof_indices);
924
925 data.local_grad = 0.;
926 for (unsigned int q = 0; q < scratch.nqp; ++q)
927 {
928 for (unsigned int i = 0; i < data.vel_dpc; ++i)
929 for (unsigned int j = 0; j < data.pres_dpc; ++j)
930 data.local_grad(i, j) +=
931 -scratch.fe_val_vel.JxW(q) *
932 scratch.fe_val_vel.shape_grad(i, q)[data.d] *
933 scratch.fe_val_pres.shape_value(j, q);
934 }
935 }
936
937
938 template <int dim>
copy_gradient_local_to_global(const InitGradPerTaskData & data)939 void NavierStokesProjection<dim>::copy_gradient_local_to_global(
940 const InitGradPerTaskData &data)
941 {
942 for (unsigned int i = 0; i < data.vel_dpc; ++i)
943 for (unsigned int j = 0; j < data.pres_dpc; ++j)
944 pres_Diff[data.d].add(data.vel_local_dof_indices[i],
945 data.pres_local_dof_indices[j],
946 data.local_grad(i, j));
947 }
948
949
950 // @sect4{ <code>NavierStokesProjection::run</code> }
951
952 // This is the time marching function, which starting at <code>t_0</code>
953 // advances in time using the projection method with time step <code>dt</code>
954 // until <code>T</code>.
955 //
956 // Its second parameter, <code>verbose</code> indicates whether the function
957 // should output information what it is doing at any given moment: for
958 // example, it will say whether we are working on the diffusion, projection
959 // substep; updating preconditioners etc. Rather than implementing this
960 // output using code like
961 // @code
962 // if (verbose) std::cout << "something";
963 // @endcode
964 // we use the ConditionalOStream class to do that for us. That
965 // class takes an output stream and a condition that indicates whether the
966 // things you pass to it should be passed through to the given output
967 // stream, or should just be ignored. This way, above code simply becomes
968 // @code
969 // verbose_cout << "something";
970 // @endcode
971 // and does the right thing in either case.
972 template <int dim>
run(const bool verbose,const unsigned int output_interval)973 void NavierStokesProjection<dim>::run(const bool verbose,
974 const unsigned int output_interval)
975 {
976 ConditionalOStream verbose_cout(std::cout, verbose);
977
978 const auto n_steps = static_cast<unsigned int>((T - t_0) / dt);
979 vel_exact.set_time(2. * dt);
980 output_results(1);
981 for (unsigned int n = 2; n <= n_steps; ++n)
982 {
983 if (n % output_interval == 0)
984 {
985 verbose_cout << "Plotting Solution" << std::endl;
986 output_results(n);
987 }
988 std::cout << "Step = " << n << " Time = " << (n * dt) << std::endl;
989 verbose_cout << " Interpolating the velocity " << std::endl;
990
991 interpolate_velocity();
992 verbose_cout << " Diffusion Step" << std::endl;
993 if (n % vel_update_prec == 0)
994 verbose_cout << " With reinitialization of the preconditioner"
995 << std::endl;
996 diffusion_step((n % vel_update_prec == 0) || (n == 2));
997 verbose_cout << " Projection Step" << std::endl;
998 projection_step((n == 2));
999 verbose_cout << " Updating the Pressure" << std::endl;
1000 update_pressure((n == 2));
1001 vel_exact.advance_time(dt);
1002 }
1003 output_results(n_steps);
1004 }
1005
1006
1007
1008 template <int dim>
interpolate_velocity()1009 void NavierStokesProjection<dim>::interpolate_velocity()
1010 {
1011 for (unsigned int d = 0; d < dim; ++d)
1012 {
1013 u_star[d].equ(2., u_n[d]);
1014 u_star[d] -= u_n_minus_1[d];
1015 }
1016 }
1017
1018
1019 // @sect4{<code>NavierStokesProjection::diffusion_step</code>}
1020
1021 // The implementation of a diffusion step. Note that the expensive operation
1022 // is the diffusion solve at the end of the function, which we have to do once
1023 // for each velocity component. To accelerate things a bit, we allow to do
1024 // this in %parallel, using the Threads::new_task function which makes sure
1025 // that the <code>dim</code> solves are all taken care of and are scheduled to
1026 // available processors: if your machine has more than one processor core and
1027 // no other parts of this program are using resources currently, then the
1028 // diffusion solves will run in %parallel. On the other hand, if your system
1029 // has only one processor core then running things in %parallel would be
1030 // inefficient (since it leads, for example, to cache congestion) and things
1031 // will be executed sequentially.
1032 template <int dim>
diffusion_step(const bool reinit_prec)1033 void NavierStokesProjection<dim>::diffusion_step(const bool reinit_prec)
1034 {
1035 pres_tmp.equ(-1., pres_n);
1036 pres_tmp.add(-4. / 3., phi_n, 1. / 3., phi_n_minus_1);
1037
1038 assemble_advection_term();
1039
1040 for (unsigned int d = 0; d < dim; ++d)
1041 {
1042 force[d] = 0.;
1043 v_tmp.equ(2. / dt, u_n[d]);
1044 v_tmp.add(-.5 / dt, u_n_minus_1[d]);
1045 vel_Mass.vmult_add(force[d], v_tmp);
1046
1047 pres_Diff[d].vmult_add(force[d], pres_tmp);
1048 u_n_minus_1[d] = u_n[d];
1049
1050 vel_it_matrix[d].copy_from(vel_Laplace_plus_Mass);
1051 vel_it_matrix[d].add(1., vel_Advection);
1052
1053 vel_exact.set_component(d);
1054 boundary_values.clear();
1055 for (const auto &boundary_id : boundary_ids)
1056 {
1057 switch (boundary_id)
1058 {
1059 case 1:
1060 VectorTools::interpolate_boundary_values(
1061 dof_handler_velocity,
1062 boundary_id,
1063 Functions::ZeroFunction<dim>(),
1064 boundary_values);
1065 break;
1066 case 2:
1067 VectorTools::interpolate_boundary_values(dof_handler_velocity,
1068 boundary_id,
1069 vel_exact,
1070 boundary_values);
1071 break;
1072 case 3:
1073 if (d != 0)
1074 VectorTools::interpolate_boundary_values(
1075 dof_handler_velocity,
1076 boundary_id,
1077 Functions::ZeroFunction<dim>(),
1078 boundary_values);
1079 break;
1080 case 4:
1081 VectorTools::interpolate_boundary_values(
1082 dof_handler_velocity,
1083 boundary_id,
1084 Functions::ZeroFunction<dim>(),
1085 boundary_values);
1086 break;
1087 default:
1088 Assert(false, ExcNotImplemented());
1089 }
1090 }
1091 MatrixTools::apply_boundary_values(boundary_values,
1092 vel_it_matrix[d],
1093 u_n[d],
1094 force[d]);
1095 }
1096
1097
1098 Threads::TaskGroup<void> tasks;
1099 for (unsigned int d = 0; d < dim; ++d)
1100 {
1101 if (reinit_prec)
1102 prec_velocity[d].initialize(vel_it_matrix[d],
1103 SparseILU<double>::AdditionalData(
1104 vel_diag_strength, vel_off_diagonals));
1105 tasks += Threads::new_task(
1106 &NavierStokesProjection<dim>::diffusion_component_solve, *this, d);
1107 }
1108 tasks.join_all();
1109 }
1110
1111
1112
1113 template <int dim>
1114 void
diffusion_component_solve(const unsigned int d)1115 NavierStokesProjection<dim>::diffusion_component_solve(const unsigned int d)
1116 {
1117 SolverControl solver_control(vel_max_its, vel_eps * force[d].l2_norm());
1118 SolverGMRES<Vector<double>> gmres(
1119 solver_control,
1120 SolverGMRES<Vector<double>>::AdditionalData(vel_Krylov_size));
1121 gmres.solve(vel_it_matrix[d], u_n[d], force[d], prec_velocity[d]);
1122 }
1123
1124
1125 // @sect4{ <code>NavierStokesProjection::assemble_advection_term</code> }
1126
1127 // The following few functions deal with assembling the advection terms, which
1128 // is the part of the system matrix for the diffusion step that changes at
1129 // every time step. As mentioned above, we will run the assembly loop over all
1130 // cells in %parallel, using the WorkStream class and other
1131 // facilities as described in the documentation module on @ref threads.
1132 template <int dim>
assemble_advection_term()1133 void NavierStokesProjection<dim>::assemble_advection_term()
1134 {
1135 vel_Advection = 0.;
1136 AdvectionPerTaskData data(fe_velocity.n_dofs_per_cell());
1137 AdvectionScratchData scratch(fe_velocity,
1138 quadrature_velocity,
1139 update_values | update_JxW_values |
1140 update_gradients);
1141 WorkStream::run(
1142 dof_handler_velocity.begin_active(),
1143 dof_handler_velocity.end(),
1144 *this,
1145 &NavierStokesProjection<dim>::assemble_one_cell_of_advection,
1146 &NavierStokesProjection<dim>::copy_advection_local_to_global,
1147 scratch,
1148 data);
1149 }
1150
1151
1152
1153 template <int dim>
assemble_one_cell_of_advection(const typename DoFHandler<dim>::active_cell_iterator & cell,AdvectionScratchData & scratch,AdvectionPerTaskData & data)1154 void NavierStokesProjection<dim>::assemble_one_cell_of_advection(
1155 const typename DoFHandler<dim>::active_cell_iterator &cell,
1156 AdvectionScratchData & scratch,
1157 AdvectionPerTaskData & data)
1158 {
1159 scratch.fe_val.reinit(cell);
1160 cell->get_dof_indices(data.local_dof_indices);
1161 for (unsigned int d = 0; d < dim; ++d)
1162 {
1163 scratch.fe_val.get_function_values(u_star[d], scratch.u_star_tmp);
1164 for (unsigned int q = 0; q < scratch.nqp; ++q)
1165 scratch.u_star_local[q](d) = scratch.u_star_tmp[q];
1166 }
1167
1168 for (unsigned int d = 0; d < dim; ++d)
1169 {
1170 scratch.fe_val.get_function_gradients(u_star[d], scratch.grad_u_star);
1171 for (unsigned int q = 0; q < scratch.nqp; ++q)
1172 {
1173 if (d == 0)
1174 scratch.u_star_tmp[q] = 0.;
1175 scratch.u_star_tmp[q] += scratch.grad_u_star[q][d];
1176 }
1177 }
1178
1179 data.local_advection = 0.;
1180 for (unsigned int q = 0; q < scratch.nqp; ++q)
1181 for (unsigned int i = 0; i < scratch.dpc; ++i)
1182 for (unsigned int j = 0; j < scratch.dpc; ++j)
1183 data.local_advection(i, j) += (scratch.u_star_local[q] * //
1184 scratch.fe_val.shape_grad(j, q) * //
1185 scratch.fe_val.shape_value(i, q) //
1186 + //
1187 0.5 * //
1188 scratch.u_star_tmp[q] * //
1189 scratch.fe_val.shape_value(i, q) * //
1190 scratch.fe_val.shape_value(j, q)) //
1191 * scratch.fe_val.JxW(q);
1192 }
1193
1194
1195
1196 template <int dim>
copy_advection_local_to_global(const AdvectionPerTaskData & data)1197 void NavierStokesProjection<dim>::copy_advection_local_to_global(
1198 const AdvectionPerTaskData &data)
1199 {
1200 for (unsigned int i = 0; i < fe_velocity.n_dofs_per_cell(); ++i)
1201 for (unsigned int j = 0; j < fe_velocity.n_dofs_per_cell(); ++j)
1202 vel_Advection.add(data.local_dof_indices[i],
1203 data.local_dof_indices[j],
1204 data.local_advection(i, j));
1205 }
1206
1207
1208
1209 // @sect4{<code>NavierStokesProjection::projection_step</code>}
1210
1211 // This implements the projection step:
1212 template <int dim>
projection_step(const bool reinit_prec)1213 void NavierStokesProjection<dim>::projection_step(const bool reinit_prec)
1214 {
1215 pres_iterative.copy_from(pres_Laplace);
1216
1217 pres_tmp = 0.;
1218 for (unsigned d = 0; d < dim; ++d)
1219 pres_Diff[d].Tvmult_add(pres_tmp, u_n[d]);
1220
1221 phi_n_minus_1 = phi_n;
1222
1223 static std::map<types::global_dof_index, double> bval;
1224 if (reinit_prec)
1225 VectorTools::interpolate_boundary_values(dof_handler_pressure,
1226 3,
1227 Functions::ZeroFunction<dim>(),
1228 bval);
1229
1230 MatrixTools::apply_boundary_values(bval, pres_iterative, phi_n, pres_tmp);
1231
1232 if (reinit_prec)
1233 prec_pres_Laplace.initialize(pres_iterative,
1234 SparseILU<double>::AdditionalData(
1235 vel_diag_strength, vel_off_diagonals));
1236
1237 SolverControl solvercontrol(vel_max_its, vel_eps * pres_tmp.l2_norm());
1238 SolverCG<Vector<double>> cg(solvercontrol);
1239 cg.solve(pres_iterative, phi_n, pres_tmp, prec_pres_Laplace);
1240
1241 phi_n *= 1.5 / dt;
1242 }
1243
1244
1245 // @sect4{ <code>NavierStokesProjection::update_pressure</code> }
1246
1247 // This is the pressure update step of the projection method. It implements
1248 // the standard formulation of the method, that is @f[ p^{n+1} = p^n +
1249 // \phi^{n+1}, @f] or the rotational form, which is @f[ p^{n+1} = p^n +
1250 // \phi^{n+1} - \frac{1}{Re} \nabla\cdot u^{n+1}. @f]
1251 template <int dim>
update_pressure(const bool reinit_prec)1252 void NavierStokesProjection<dim>::update_pressure(const bool reinit_prec)
1253 {
1254 pres_n_minus_1 = pres_n;
1255 switch (type)
1256 {
1257 case RunTimeParameters::Method::standard:
1258 pres_n += phi_n;
1259 break;
1260 case RunTimeParameters::Method::rotational:
1261 if (reinit_prec)
1262 prec_mass.initialize(pres_Mass);
1263 pres_n = pres_tmp;
1264 prec_mass.solve(pres_n);
1265 pres_n.sadd(1. / Re, 1., pres_n_minus_1);
1266 pres_n += phi_n;
1267 break;
1268 default:
1269 Assert(false, ExcNotImplemented());
1270 };
1271 }
1272
1273
1274 // @sect4{ <code>NavierStokesProjection::output_results</code> }
1275
1276 // This method plots the current solution. The main difficulty is that we want
1277 // to create a single output file that contains the data for all velocity
1278 // components, the pressure, and also the vorticity of the flow. On the other
1279 // hand, velocities and the pressure live on separate DoFHandler objects, and
1280 // so can't be written to the same file using a single DataOut object. As a
1281 // consequence, we have to work a bit harder to get the various pieces of data
1282 // into a single DoFHandler object, and then use that to drive graphical
1283 // output.
1284 //
1285 // We will not elaborate on this process here, but rather refer to step-32,
1286 // where a similar procedure is used (and is documented) to create a joint
1287 // DoFHandler object for all variables.
1288 //
1289 // Let us also note that we here compute the vorticity as a scalar quantity in
1290 // a separate function, using the $L^2$ projection of the quantity
1291 // $\text{curl} u$ onto the finite element space used for the components of
1292 // the velocity. In principle, however, we could also have computed as a
1293 // pointwise quantity from the velocity, and do so through the
1294 // DataPostprocessor mechanism discussed in step-29 and step-33.
1295 template <int dim>
output_results(const unsigned int step)1296 void NavierStokesProjection<dim>::output_results(const unsigned int step)
1297 {
1298 assemble_vorticity((step == 1));
1299 const FESystem<dim> joint_fe(
1300 fe_velocity, dim, fe_pressure, 1, fe_velocity, 1);
1301 DoFHandler<dim> joint_dof_handler(triangulation);
1302 joint_dof_handler.distribute_dofs(joint_fe);
1303 Assert(joint_dof_handler.n_dofs() ==
1304 ((dim + 1) * dof_handler_velocity.n_dofs() +
1305 dof_handler_pressure.n_dofs()),
1306 ExcInternalError());
1307 Vector<double> joint_solution(joint_dof_handler.n_dofs());
1308 std::vector<types::global_dof_index> loc_joint_dof_indices(
1309 joint_fe.n_dofs_per_cell()),
1310 loc_vel_dof_indices(fe_velocity.n_dofs_per_cell()),
1311 loc_pres_dof_indices(fe_pressure.n_dofs_per_cell());
1312 typename DoFHandler<dim>::active_cell_iterator
1313 joint_cell = joint_dof_handler.begin_active(),
1314 joint_endc = joint_dof_handler.end(),
1315 vel_cell = dof_handler_velocity.begin_active(),
1316 pres_cell = dof_handler_pressure.begin_active();
1317 for (; joint_cell != joint_endc; ++joint_cell, ++vel_cell, ++pres_cell)
1318 {
1319 joint_cell->get_dof_indices(loc_joint_dof_indices);
1320 vel_cell->get_dof_indices(loc_vel_dof_indices);
1321 pres_cell->get_dof_indices(loc_pres_dof_indices);
1322 for (unsigned int i = 0; i < joint_fe.n_dofs_per_cell(); ++i)
1323 switch (joint_fe.system_to_base_index(i).first.first)
1324 {
1325 case 0:
1326 Assert(joint_fe.system_to_base_index(i).first.second < dim,
1327 ExcInternalError());
1328 joint_solution(loc_joint_dof_indices[i]) =
1329 u_n[joint_fe.system_to_base_index(i).first.second](
1330 loc_vel_dof_indices[joint_fe.system_to_base_index(i)
1331 .second]);
1332 break;
1333 case 1:
1334 Assert(joint_fe.system_to_base_index(i).first.second == 0,
1335 ExcInternalError());
1336 joint_solution(loc_joint_dof_indices[i]) =
1337 pres_n(loc_pres_dof_indices[joint_fe.system_to_base_index(i)
1338 .second]);
1339 break;
1340 case 2:
1341 Assert(joint_fe.system_to_base_index(i).first.second == 0,
1342 ExcInternalError());
1343 joint_solution(loc_joint_dof_indices[i]) = rot_u(
1344 loc_vel_dof_indices[joint_fe.system_to_base_index(i).second]);
1345 break;
1346 default:
1347 Assert(false, ExcInternalError());
1348 }
1349 }
1350 std::vector<std::string> joint_solution_names(dim, "v");
1351 joint_solution_names.emplace_back("p");
1352 joint_solution_names.emplace_back("rot_u");
1353 DataOut<dim> data_out;
1354 data_out.attach_dof_handler(joint_dof_handler);
1355 std::vector<DataComponentInterpretation::DataComponentInterpretation>
1356 component_interpretation(
1357 dim + 2, DataComponentInterpretation::component_is_part_of_vector);
1358 component_interpretation[dim] =
1359 DataComponentInterpretation::component_is_scalar;
1360 component_interpretation[dim + 1] =
1361 DataComponentInterpretation::component_is_scalar;
1362 data_out.add_data_vector(joint_solution,
1363 joint_solution_names,
1364 DataOut<dim>::type_dof_data,
1365 component_interpretation);
1366 data_out.build_patches(deg + 1);
1367 std::ofstream output("solution-" + Utilities::int_to_string(step, 5) +
1368 ".vtk");
1369 data_out.write_vtk(output);
1370 }
1371
1372
1373
1374 // Following is the helper function that computes the vorticity by projecting
1375 // the term $\text{curl} u$ onto the finite element space used for the
1376 // components of the velocity. The function is only called whenever we
1377 // generate graphical output, so not very often, and as a consequence we
1378 // didn't bother parallelizing it using the WorkStream concept as we do for
1379 // the other assembly functions. That should not be overly complicated,
1380 // however, if needed. Moreover, the implementation that we have here only
1381 // works for 2d, so we bail if that is not the case.
1382 template <int dim>
assemble_vorticity(const bool reinit_prec)1383 void NavierStokesProjection<dim>::assemble_vorticity(const bool reinit_prec)
1384 {
1385 Assert(dim == 2, ExcNotImplemented());
1386 if (reinit_prec)
1387 prec_vel_mass.initialize(vel_Mass);
1388
1389 FEValues<dim> fe_val_vel(fe_velocity,
1390 quadrature_velocity,
1391 update_gradients | update_JxW_values |
1392 update_values);
1393 const unsigned int dpc = fe_velocity.n_dofs_per_cell(),
1394 nqp = quadrature_velocity.size();
1395 std::vector<types::global_dof_index> ldi(dpc);
1396 Vector<double> loc_rot(dpc);
1397
1398 std::vector<Tensor<1, dim>> grad_u1(nqp), grad_u2(nqp);
1399 rot_u = 0.;
1400
1401 for (const auto &cell : dof_handler_velocity.active_cell_iterators())
1402 {
1403 fe_val_vel.reinit(cell);
1404 cell->get_dof_indices(ldi);
1405 fe_val_vel.get_function_gradients(u_n[0], grad_u1);
1406 fe_val_vel.get_function_gradients(u_n[1], grad_u2);
1407 loc_rot = 0.;
1408 for (unsigned int q = 0; q < nqp; ++q)
1409 for (unsigned int i = 0; i < dpc; ++i)
1410 loc_rot(i) += (grad_u2[q][0] - grad_u1[q][1]) * //
1411 fe_val_vel.shape_value(i, q) * //
1412 fe_val_vel.JxW(q);
1413
1414 for (unsigned int i = 0; i < dpc; ++i)
1415 rot_u(ldi[i]) += loc_rot(i);
1416 }
1417
1418 prec_vel_mass.solve(rot_u);
1419 }
1420 } // namespace Step35
1421
1422
1423 // @sect3{ The main function }
1424
1425 // The main function looks very much like in all the other tutorial programs, so
1426 // there is little to comment on here:
main()1427 int main()
1428 {
1429 try
1430 {
1431 using namespace Step35;
1432
1433 RunTimeParameters::Data_Storage data;
1434 data.read_data("parameter-file.prm");
1435
1436 deallog.depth_console(data.verbose ? 2 : 0);
1437
1438 NavierStokesProjection<2> test(data);
1439 test.run(data.verbose, data.output_interval);
1440 }
1441 catch (std::exception &exc)
1442 {
1443 std::cerr << std::endl
1444 << std::endl
1445 << "----------------------------------------------------"
1446 << std::endl;
1447 std::cerr << "Exception on processing: " << std::endl
1448 << exc.what() << std::endl
1449 << "Aborting!" << std::endl
1450 << "----------------------------------------------------"
1451 << std::endl;
1452 return 1;
1453 }
1454 catch (...)
1455 {
1456 std::cerr << std::endl
1457 << std::endl
1458 << "----------------------------------------------------"
1459 << std::endl;
1460 std::cerr << "Unknown exception!" << std::endl
1461 << "Aborting!" << std::endl
1462 << "----------------------------------------------------"
1463 << std::endl;
1464 return 1;
1465 }
1466 std::cout << "----------------------------------------------------"
1467 << std::endl
1468 << "Apparently everything went fine!" << std::endl
1469 << "Don't forget to brush your teeth :-)" << std::endl
1470 << std::endl;
1471 return 0;
1472 }
1473