266 * <a name=
"time_dependent_navier_stokes.cc-Createthetriangulation"></a>
270 * code](https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc)
271 * with very few modifications.
276 * <a name="time_dependent_navier_stokes.cc-Helperfunction"></a>
277 * <h4>Helper function</h4>
280 * void create_triangulation_2d(Triangulation<2> &tria,
281 * bool compute_in_2d = true)
283 * SphericalManifold<2> boundary(Point<2>(0.5, 0.2));
284 * Triangulation<2> left, middle, right, tmp, tmp2;
285 * GridGenerator::subdivided_hyper_rectangle(
287 * std::vector<unsigned int>({3U, 4U}),
289 * Point<2>(0.3, 0.41),
291 * GridGenerator::subdivided_hyper_rectangle(
293 * std::vector<unsigned int>({18U, 4U}),
295 * Point<2>(2.5, 0.41),
300 * Create middle part first as a hyper shell.
303 * GridGenerator::hyper_shell(middle, Point<2>(0.5, 0.2), 0.05, 0.2, 4, true);
304 * middle.reset_all_manifolds();
305 * for (Triangulation<2>::cell_iterator cell = middle.begin();
306 * cell != middle.end();
308 * for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
310 * bool is_inner_rim = true;
311 * for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_face; ++v)
313 * Point<2> &vertex = cell->face(f)->vertex(v);
314 * if (std::abs(vertex.distance(Point<2>(0.5, 0.2)) - 0.05) > 1e-10)
316 * is_inner_rim = false;
321 * cell->face(f)->set_manifold_id(1);
323 * middle.set_manifold(1, boundary);
324 * middle.refine_global(1);
328 * Then move the vertices to the points where we want them to be to create a
329 * slightly asymmetric cube with a hole:
332 * for (Triangulation<2>::cell_iterator cell = middle.begin();
333 * cell != middle.end();
335 * for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
337 * Point<2> &vertex = cell->vertex(v);
338 * if (std::abs(vertex[0] - 0.7) < 1e-10 &&
339 * std::abs(vertex[1] - 0.2) < 1e-10)
340 * vertex = Point<2>(0.7, 0.205);
341 * else if (std::abs(vertex[0] - 0.6) < 1e-10 &&
342 * std::abs(vertex[1] - 0.3) < 1e-10)
343 * vertex = Point<2>(0.7, 0.41);
344 * else if (std::abs(vertex[0] - 0.6) < 1e-10 &&
345 * std::abs(vertex[1] - 0.1) < 1e-10)
346 * vertex = Point<2>(0.7, 0);
347 * else if (std::abs(vertex[0] - 0.5) < 1e-10 &&
348 * std::abs(vertex[1] - 0.4) < 1e-10)
349 * vertex = Point<2>(0.5, 0.41);
350 * else if (std::abs(vertex[0] - 0.5) < 1e-10 &&
351 * std::abs(vertex[1] - 0.0) < 1e-10)
352 * vertex = Point<2>(0.5, 0.0);
353 * else if (std::abs(vertex[0] - 0.4) < 1e-10 &&
354 * std::abs(vertex[1] - 0.3) < 1e-10)
355 * vertex = Point<2>(0.3, 0.41);
356 * else if (std::abs(vertex[0] - 0.4) < 1e-10 &&
357 * std::abs(vertex[1] - 0.1) < 1e-10)
358 * vertex = Point<2>(0.3, 0);
359 * else if (std::abs(vertex[0] - 0.3) < 1e-10 &&
360 * std::abs(vertex[1] - 0.2) < 1e-10)
361 * vertex = Point<2>(0.3, 0.205);
362 * else if (std::abs(vertex[0] - 0.56379) < 1e-4 &&
363 * std::abs(vertex[1] - 0.13621) < 1e-4)
364 * vertex = Point<2>(0.59, 0.11);
365 * else if (std::abs(vertex[0] - 0.56379) < 1e-4 &&
366 * std::abs(vertex[1] - 0.26379) < 1e-4)
367 * vertex = Point<2>(0.59, 0.29);
368 * else if (std::abs(vertex[0] - 0.43621) < 1e-4 &&
369 * std::abs(vertex[1] - 0.13621) < 1e-4)
370 * vertex = Point<2>(0.41, 0.11);
371 * else if (std::abs(vertex[0] - 0.43621) < 1e-4 &&
372 * std::abs(vertex[1] - 0.26379) < 1e-4)
373 * vertex = Point<2>(0.41, 0.29);
378 * Refine once to create the same level of refinement as in the
379 * neighboring domains:
382 * middle.refine_global(1);
386 * Must copy the triangulation because we cannot merge triangulations with
390 * GridGenerator::flatten_triangulation(middle, tmp2);
394 * Left domain is required in 3d only.
399 * GridGenerator::merge_triangulations(tmp2, right, tria);
403 * GridGenerator::merge_triangulations(left, tmp2, tmp);
404 * GridGenerator::merge_triangulations(tmp, right, tria);
411 * <a name="time_dependent_navier_stokes.cc-2Dflowaroundcylindertriangulation"></a>
412 * <h4>2D flow around cylinder triangulation</h4>
415 * void create_triangulation(Triangulation<2> &tria)
417 * create_triangulation_2d(tria);
420 * Set the left boundary (inflow) to 0, the right boundary (outflow) to 1,
421 * upper to 2, lower to 3 and the cylindrical surface to 4.
424 * for (Triangulation<2>::active_cell_iterator cell = tria.begin();
425 * cell != tria.end();
428 * for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
430 * if (cell->face(f)->at_boundary())
432 * if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
434 * cell->face(f)->set_all_boundary_ids(1);
436 * else if (std::abs(cell->face(f)->center()[0] - 0.3) < 1e-12)
438 * cell->face(f)->set_all_boundary_ids(0);
440 * else if (std::abs(cell->face(f)->center()[1] - 0.41) < 1e-12)
442 * cell->face(f)->set_all_boundary_ids(3);
444 * else if (std::abs(cell->face(f)->center()[1]) < 1e-12)
446 * cell->face(f)->set_all_boundary_ids(2);
450 * cell->face(f)->set_all_boundary_ids(4);
460 * <a name="time_dependent_navier_stokes.cc-3Dflowaroundcylindertriangulation"></a>
461 * <h4>3D flow around cylinder triangulation</h4>
464 * void create_triangulation(Triangulation<3> &tria)
466 * Triangulation<2> tria_2d;
467 * create_triangulation_2d(tria_2d, false);
468 * GridGenerator::extrude_triangulation(tria_2d, 5, 0.41, tria);
471 * Set the ids of the boundaries in x direction to 0 and 1; y direction to 2 and 3;
472 * z direction to 4 and 5; the cylindrical surface 6.
475 * for (Triangulation<3>::active_cell_iterator cell = tria.begin();
476 * cell != tria.end();
479 * for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
481 * if (cell->face(f)->at_boundary())
483 * if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
485 * cell->face(f)->set_all_boundary_ids(1);
487 * else if (std::abs(cell->face(f)->center()[0]) < 1e-12)
489 * cell->face(f)->set_all_boundary_ids(0);
491 * else if (std::abs(cell->face(f)->center()[1] - 0.41) < 1e-12)
493 * cell->face(f)->set_all_boundary_ids(3);
495 * else if (std::abs(cell->face(f)->center()[1]) < 1e-12)
497 * cell->face(f)->set_all_boundary_ids(2);
499 * else if (std::abs(cell->face(f)->center()[2] - 0.41) < 1e-12)
501 * cell->face(f)->set_all_boundary_ids(5);
503 * else if (std::abs(cell->face(f)->center()[2]) < 1e-12)
505 * cell->face(f)->set_all_boundary_ids(4);
509 * cell->face(f)->set_all_boundary_ids(6);
519 * <a name="time_dependent_navier_stokes.cc-Timestepping"></a>
520 * <h3>Time stepping</h3>
521 * This class is pretty much self-explanatory.
527 * Time(const double time_end,
528 * const double delta_t,
529 * const double output_interval,
530 * const double refinement_interval)
533 * time_end(time_end),
535 * output_interval(output_interval),
536 * refinement_interval(refinement_interval)
539 * double current() const { return time_current; }
540 * double end() const { return time_end; }
541 * double get_delta_t() const { return delta_t; }
542 * unsigned int get_timestep() const { return timestep; }
543 * bool time_to_output() const;
544 * bool time_to_refine() const;
548 * unsigned int timestep;
549 * double time_current;
550 * const double time_end;
551 * const double delta_t;
552 * const double output_interval;
553 * const double refinement_interval;
556 * bool Time::time_to_output() const
558 * unsigned int delta = static_cast<unsigned int>(output_interval / delta_t);
559 * return (timestep >= delta && timestep % delta == 0);
562 * bool Time::time_to_refine() const
564 * unsigned int delta = static_cast<unsigned int>(refinement_interval / delta_t);
565 * return (timestep >= delta && timestep % delta == 0);
568 * void Time::increment()
570 * time_current += delta_t;
577 * <a name="time_dependent_navier_stokes.cc-Boundaryvalues"></a>
578 * <h3>Boundary values</h3>
579 * Dirichlet boundary conditions for the velocity inlet and walls.
583 * class BoundaryValues : public Function<dim>
586 * BoundaryValues() : Function<dim>(dim + 1) {}
587 * virtual double value(const Point<dim> &p,
588 * const unsigned int component) const override;
590 * virtual void vector_value(const Point<dim> &p,
591 * Vector<double> &values) const override;
595 * double BoundaryValues<dim>::value(const Point<dim> &p,
596 * const unsigned int component) const
598 * Assert(component < this->n_components,
599 * ExcIndexRange(component, 0, this->n_components));
600 * double left_boundary = (dim == 2 ? 0.3 : 0.0);
601 * if (component == 0 && std::abs(p[0] - left_boundary) < 1e-10)
605 * For a parabolic velocity profile, @f$U_\mathrm{avg} = 2/3
607 * in 2D, and @f$U_\mathrm{avg} = 4/9 U_\mathrm{max}@f$ in 3D.
608 * If @f$\nu = 0.001@f$, @f$D = 0.1@f$, then @f$Re = 100 U_\mathrm{avg}@f$.
612 * double Umax = (dim == 2 ? 3 * Uavg / 2 : 9 * Uavg / 4);
613 * double value = 4 * Umax * p[1] * (0.41 - p[1]) / (0.41 * 0.41);
616 * value *= 4 * p[2] * (0.41 - p[2]) / (0.41 * 0.41);
624 * void BoundaryValues<dim>::vector_value(const Point<dim> &p,
625 * Vector<double> &values) const
627 * for (unsigned int c = 0; c < this->n_components; ++c)
628 * values(c) = BoundaryValues<dim>::value(p, c);
634 * <a name="time_dependent_navier_stokes.cc-Blockpreconditioner"></a>
635 * <h3>Block preconditioner</h3>
639 * The block Schur preconditioner can be written as the product of three
642 * P^{-1} = \begin{pmatrix} \tilde{A}^{-1} & 0\\ 0 & I\end{pmatrix}
643 * \begin{pmatrix} I & -B^T\\ 0 & I\end{pmatrix}
644 * \begin{pmatrix} I & 0\\ 0 & \tilde{S}^{-1}\end{pmatrix}
646 * @f$\tilde{A}@f$ is symmetric since the convection term is eliminated from the
648 * @f$\tilde{S}^{-1}@f$ is the inverse of the Schur complement of @f$\tilde{A}@f$,
649 * which consists of a reaction term, a diffusion term, a Grad-Div term
650 * and a convection term.
651 * In practice, the convection contribution is ignored, namely
652 * @f$\tilde{S}^{-1} = -(\nu + \gamma)M_p^{-1} -
653 * \frac{1}{\Delta{t}}{[B(diag(M_u))^{-1}B^T]}^{-1}@f$
654 * where @f$M_p@f$ is the pressure mass, and
655 * @f${[B(diag(M_u))^{-1}B^T]}@f$ is an approximation to the Schur complement of
656 * (velocity) mass matrix @f$BM_u^{-1}B^T@f$.
660 * Same as the tutorials, we define a vmult operation for the block
662 * instead of write it as a matrix. It can be seen from the above definition,
663 * the result of the vmult operation of the block preconditioner can be
665 * from the results of the vmult operations of @f$M_u^{-1}@f$, @f$M_p^{-1}@f$,
666 * @f$\tilde{A}^{-1}@f$, which can be transformed into solving three symmetric
671 * class BlockSchurPreconditioner : public Subscriptor
674 * BlockSchurPreconditioner(
675 * TimerOutput &timer,
679 * const std::vector<IndexSet> &owned_partitioning,
680 * const PETScWrappers::MPI::BlockSparseMatrix &system,
681 * const PETScWrappers::MPI::BlockSparseMatrix &mass,
682 * PETScWrappers::MPI::BlockSparseMatrix &schur);
684 * void vmult(PETScWrappers::MPI::BlockVector &dst,
685 * const PETScWrappers::MPI::BlockVector &src) const;
688 * TimerOutput &timer;
689 * const double gamma;
690 * const double viscosity;
693 * const SmartPointer<const PETScWrappers::MPI::BlockSparseMatrix>
695 * const SmartPointer<const PETScWrappers::MPI::BlockSparseMatrix> mass_matrix;
698 * As discussed, @f${[B(diag(M_u))^{-1}B^T]}@f$ and its inverse
699 * need to be computed.
700 * We can either explicitly compute it out as a matrix, or define
701 * it as a class with a vmult operation.
702 * The second approach saves some computation to construct the matrix,
703 * but leads to slow convergence in CG solver because it is impossible
704 * to apply a preconditioner. We go with the first route.
707 * const SmartPointer<PETScWrappers::MPI::BlockSparseMatrix> mass_schur;
713 * <a name="time_dependent_navier_stokes.cc-BlockSchurPreconditionerBlockSchurPreconditioner"></a>
714 * <h4>BlockSchurPreconditioner::BlockSchurPreconditioner</h4>
718 * Input parameters and system matrix, mass matrix as well as the mass schur
719 * matrix are needed in the preconditioner. In addition, we pass the
720 * partitioning information into this class because we need to create some
721 * temporary block vectors inside.
724 * BlockSchurPreconditioner::BlockSchurPreconditioner(
725 * TimerOutput &timer,
729 * const std::vector<IndexSet> &owned_partitioning,
730 * const PETScWrappers::MPI::BlockSparseMatrix &system,
731 * const PETScWrappers::MPI::BlockSparseMatrix &mass,
732 * PETScWrappers::MPI::BlockSparseMatrix &schur)
735 * viscosity(viscosity),
737 * system_matrix(&system),
738 * mass_matrix(&mass),
741 * TimerOutput::Scope timer_section(timer, "CG for Sm");
744 * The schur complemete of mass matrix is actually being computed here.
747 * PETScWrappers::MPI::BlockVector tmp1, tmp2;
748 * tmp1.reinit(owned_partitioning, mass_matrix->get_mpi_communicator());
749 * tmp2.reinit(owned_partitioning, mass_matrix->get_mpi_communicator());
754 * Jacobi preconditioner of matrix A is by definition @f${diag(A)}^{-1}@f$,
755 * this is exactly what we want to compute.
758 * PETScWrappers::PreconditionJacobi jacobi(mass_matrix->block(0, 0));
759 * jacobi.vmult(tmp2.block(0), tmp1.block(0));
760 * system_matrix->block(1, 0).mmult(
761 * mass_schur->block(1, 1), system_matrix->block(0, 1), tmp2.block(0));
767 * <a name="time_dependent_navier_stokes.cc-BlockSchurPreconditionervmult"></a>
768 * <h4>BlockSchurPreconditioner::vmult</h4>
772 * The vmult operation strictly follows the definition of
773 * BlockSchurPreconditioner
774 * introduced above. Conceptually it computes @f$u = P^{-1}v@f$.
777 * void BlockSchurPreconditioner::vmult(
778 * PETScWrappers::MPI::BlockVector &dst,
779 * const PETScWrappers::MPI::BlockVector &src) const
786 * PETScWrappers::MPI::Vector utmp(src.block(0));
787 * PETScWrappers::MPI::Vector tmp(src.block(1));
791 * This block computes @f$u_1 = \tilde{S}^{-1} v_1@f$,
792 * where CG solvers are used for @f$M_p^{-1}@f$ and @f$S_m^{-1}@f$.
796 * TimerOutput::Scope timer_section(timer, "CG for Mp");
797 * SolverControl mp_control(src.block(1).size(),
798 * 1e-6 * src.block(1).l2_norm());
799 * PETScWrappers::SolverCG cg_mp(mp_control,
800 * mass_schur->get_mpi_communicator());
803 * @f$-(\nu + \gamma)M_p^{-1}v_1@f$
806 * PETScWrappers::PreconditionBlockJacobi Mp_preconditioner;
807 * Mp_preconditioner.initialize(mass_matrix->block(1, 1));
809 * mass_matrix->block(1, 1), tmp, src.block(1), Mp_preconditioner);
810 * tmp *= -(viscosity + gamma);
814 * @f$-\frac{1}{dt}S_m^{-1}v_1@f$
818 * TimerOutput::Scope timer_section(timer, "CG for Sm");
819 * SolverControl sm_control(src.block(1).size(),
820 * 1e-6 * src.block(1).l2_norm());
821 * PETScWrappers::SolverCG cg_sm(sm_control,
822 * mass_schur->get_mpi_communicator());
825 * PreconditionBlockJacobi works find on Sm if we do not refine the mesh.
826 * Because after refine_mesh is called, zero entries will be created on
827 * the diagonal (not sure why), which prevents PreconditionBlockJacobi
831 * PETScWrappers::PreconditionNone Sm_preconditioner;
832 * Sm_preconditioner.initialize(mass_schur->block(1, 1));
834 * mass_schur->block(1, 1), dst.block(1), src.block(1), Sm_preconditioner);
835 * dst.block(1) *= -1 / dt;
839 * Adding up these two, we get @f$\tilde{S}^{-1}v_1@f$.
842 * dst.block(1) += tmp;
845 * Compute @f$v_0 - B^T\tilde{S}^{-1}v_1@f$ based on @f$u_1@f$.
848 * system_matrix->block(0, 1).vmult(utmp, dst.block(1));
850 * utmp += src.block(0);
853 * Finally, compute the product of @f$\tilde{A}^{-1}@f$ and utmp
854 * using another CG solver.
858 * TimerOutput::Scope timer_section(timer, "CG for A");
859 * SolverControl a_control(src.block(0).size(),
860 * 1e-6 * src.block(0).l2_norm());
861 * PETScWrappers::SolverCG cg_a(a_control,
862 * mass_schur->get_mpi_communicator());
865 * We do not use any preconditioner for this block, which is of course
867 * only because the performance of the only two preconditioners available
868 * PreconditionBlockJacobi and PreconditionBoomerAMG are even worse than
872 * PETScWrappers::PreconditionNone A_preconditioner;
873 * A_preconditioner.initialize(system_matrix->block(0, 0));
875 * system_matrix->block(0, 0), dst.block(0), utmp, A_preconditioner);
882 * <a name="time_dependent_navier_stokes.cc-TheincompressibleNavierStokessolver"></a>
883 * <h3>The incompressible Navier-Stokes solver</h3>
887 * Parallel incompressible Navier Stokes equation solver using
888 * implicit-explicit time scheme.
889 * This program is built upon dealii tutorials @ref step_57 "step-57", @ref step_40 "step-40", @ref step_22 "step-22",
890 * and @ref step_20 "step-20".
891 * The system equation is written in the incremental form, and we treat
892 * the convection term explicitly. Therefore the system equation is linear
893 * and symmetric, which does not need to be solved with Newton's
iteration.
904 *
~InsIMEX() { timer.print_summary(); }
913 *
void refine_mesh(
const unsigned int,
const unsigned int);
917 *
const unsigned int degree;
1001 *
std::shared_ptr<BlockSchurPreconditioner> preconditioner;
1010 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXInsIMEX"></a>
1014 *
template <
int dim>
1016 *
: viscosity(0.001),
1026 *
time(1
e0, 1e-3, 1e-2, 1e-2),
1035 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXsetup_dofs"></a>
1039 *
template <
int dim>
1047 *
dof_handler.distribute_dofs(fe);
1076 *
pcout <<
" Number of active fluid cells: "
1078 *
<<
" Number of degrees of freedom: " << dof_handler.n_dofs() <<
" ("
1079 *
<<
dof_u <<
'+' <<
dof_p <<
')' << std::endl;
1085 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXmake_constraints"></a>
1089 *
template <
int dim>
1139 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXinitialize_system"></a>
1143 *
template <
int dim>
1146 *
preconditioner.reset();
1147 *
system_matrix.
clear();
1153 *
sparsity_pattern.copy_from(
dsp);
1156 *
dof_handler.locally_owned_dofs(),
1171 *
schur_dsp.block(1, 1).compute_mmult_pattern(sparsity_pattern.block(1, 0),
1172 *
sparsity_pattern.block(0, 1));
1202 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXassemble"></a>
1203 * <
h4>InsIMEX::assemble</
h4>
1217 *
template <
int dim>
1225 *
system_matrix = 0;
1240 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1250 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1257 *
std::vector<double>
div_phi_u(dofs_per_cell);
1258 *
std::vector<Tensor<1, dim>>
phi_u(dofs_per_cell);
1259 *
std::vector<Tensor<2, dim>>
grad_phi_u(dofs_per_cell);
1260 *
std::vector<double>
phi_p(dofs_per_cell);
1262 *
for (
auto cell = dof_handler.begin_active(); cell != dof_handler.end();
1265 *
if (cell->is_locally_owned())
1267 *
fe_values.reinit(cell);
1279 *
fe_values[
velocities].get_function_gradients(
1282 *
fe_values[
velocities].get_function_divergences(
1294 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1296 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
1304 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1308 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1335 *
cell->get_dof_indices(local_dof_indices);
1343 *
local_dof_indices,
1368 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXsolve"></a>
1369 * <
h4>InsIMEX::solve</
h4>
1377 *
template <
int dim>
1378 *
std::pair<unsigned int, double>
1386 *
time.get_delta_t(),
1394 *
system_matrix.m(), 1e-8 *
system_rhs.l2_norm(),
true);
1417 *
return {solver_control.last_step(), solver_control.last_value()};
1423 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXrun"></a>
1424 * <
h4>InsIMEX::run</
h4>
1427 *
template <
int dim>
1430 *
pcout <<
"Running with PETSc on "
1432 *
<<
" MPI rank(s)..." << std::endl;
1445 *
while (time.end() - time.current() > 1e-12)
1447 *
if (time.get_timestep() == 0)
1452 *
std::cout.precision(6);
1453 *
std::cout.width(12);
1454 *
pcout << std::string(96,
'*') << std::endl
1455 *
<<
"Time step = " << time.get_timestep()
1456 *
<<
", at t = " << std::scientific << time.current() << std::endl;
1490 *
pcout << std::scientific << std::left <<
" GMRES_ITR = " << std::setw(3)
1491 *
<< state.first <<
" GMRES_RES = " << state.second << std::endl;
1497 *
if (time.time_to_output())
1501 *
if (time.time_to_refine())
1512 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXoutput_result"></a>
1519 *
template <
int dim>
1523 *
pcout <<
"Writing results..." << std::endl;
1527 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
1528 *
data_component_interpretation(
1530 *
data_component_interpretation.push_back(
1533 *
data_out.attach_dof_handler(dof_handler);
1542 *
data_component_interpretation);
1550 *
for (
unsigned int i = 0; i <
subdomain.size(); ++i)
1554 *
data_out.add_data_vector(
subdomain,
"subdomain");
1556 *
data_out.build_patches(degree + 1);
1558 *
std::string basename =
1567 *
data_out.write_vtu(output);
1572 *
for (
unsigned int i = 0;
1580 *
std::ofstream
pvd_output(
"navierstokes.pvd");
1588 * <a name=
"time_dependent_navier_stokes.cc-InsIMEXrefine_mesh"></a>
1595 *
template <
int dim>
1600 *
pcout <<
"Refining mesh..." << std::endl;
1618 *
cell->clear_refine_flag();
1625 *
cell->clear_coarsen_flag();
1635 *
trans(dof_handler);
1665 *
trans.interpolate(tmp);
1673 * <a name=
"time_dependent_navier_stokes.cc-mainfunction"></a>
1684 *
using namespace dealii;
1685 *
using namespace fluid;
1693 *
catch (std::exception &exc)
1695 *
std::cerr << std::endl
1697 *
<<
"----------------------------------------------------"
1699 *
std::cerr <<
"Exception on processing: " << std::endl
1700 *
<< exc.what() << std::endl
1701 *
<<
"Aborting!" << std::endl
1702 *
<<
"----------------------------------------------------"
1708 *
std::cerr << std::endl
1710 *
<<
"----------------------------------------------------"
1712 *
std::cerr <<
"Unknown exception!" << std::endl
1713 *
<<
"Aborting!" << std::endl
1714 *
<<
"----------------------------------------------------"
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
@ component_is_part_of_vector
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > ×_and_names)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void refine_and_coarsen_fixed_fraction(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::L1_norm)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
void advance(std::tuple< I1, I2 > &t, const unsigned int n)