step-35.h

Go to the documentation of this file.

 ((n % vel_update_prec == 0) || (n == 2)));
 *           verbose_cout << "  Projection Step" << std::endl;
 *           projection_step(/* reinit_prec= */ (n == 2));
 *           verbose_cout << "  Updating the Pressure" << std::endl;
 *           update_pressure(/* reinit_prec= */ (n == 2));
 *           vel_exact.advance_time(dt);
 *         }
 *       output_results(n_steps);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void NavierStokesProjection<dim>::interpolate_velocity()
 *     {
 *       for (unsigned int d = 0; d < dim; ++d)
 *         {
 *           u_star[d].equ(2., u_n[d]);
 *           u_star[d] -= u_n_minus_1[d];
 *         }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_35-codeNavierStokesProjectiondiffusion_stepcode"></a> 
 * <h4><code>NavierStokesProjection::diffusion_step</code></h4>
 * 
 
 * 
 * The implementation of a diffusion step. Note that the expensive operation
 * is the diffusion solve at the end of the function, which we have to do once
 * for each velocity component. To accelerate things a bit, we allow to do
 * this in %parallel, using the Threads::new_task function which makes sure
 * that the <code>dim</code> solves are all taken care of and are scheduled to
 * available processors: if your machine has more than one processor core and
 * no other parts of this program are using resources currently, then the
 * diffusion solves will run in %parallel. On the other hand, if your system
 * has only one processor core then running things in %parallel would be
 * inefficient (since it leads, for example, to cache congestion) and things
 * will be executed sequentially.
 * 
 * @code
 *     template <int dim>
 *     void NavierStokesProjection<dim>::diffusion_step(const bool reinit_prec)
 *     {
 *       pres_tmp.equ(-1., pres_n);
 *       pres_tmp.add(-4. / 3., phi_n, 1. / 3., phi_n_minus_1);
 *   
 *       assemble_advection_term();
 *   
 *       for (unsigned int d = 0; d < dim; ++d)
 *         {
 *           force[d] = 0.;
 *           v_tmp.equ(2. / dt, u_n[d]);
 *           v_tmp.add(-.5 / dt, u_n_minus_1[d]);
 *           vel_Mass.vmult_add(force[d], v_tmp);
 *   
 *           pres_Diff[d].vmult_add(force[d], pres_tmp);
 *           u_n_minus_1[d] = u_n[d];
 *   
 *           vel_it_matrix[d].copy_from(vel_Laplace_plus_Mass);
 *           vel_it_matrix[d].add(1., vel_Advection);
 *   
 *           vel_exact.set_component(d);
 *           boundary_values.clear();
 *           for (const auto &boundary_id : boundary_ids)
 *             {
 *               switch (boundary_id)
 *                 {
 *                   case 1:
 *                     VectorTools::interpolate_boundary_values(
 *                       dof_handler_velocity,
 *                       boundary_id,
 *                       Functions::ZeroFunction<dim>(),
 *                       boundary_values);
 *                     break;
 *                   case 2:
 *                     VectorTools::interpolate_boundary_values(dof_handler_velocity,
 *                                                              boundary_id,
 *                                                              vel_exact,
 *                                                              boundary_values);
 *                     break;
 *                   case 3:
 *                     if (d != 0)
 *                       VectorTools::interpolate_boundary_values(
 *                         dof_handler_velocity,
 *                         boundary_id,
 *                         Functions::ZeroFunction<dim>(),
 *                         boundary_values);
 *                     break;
 *                   case 4:
 *                     VectorTools::interpolate_boundary_values(
 *                       dof_handler_velocity,
 *                       boundary_id,
 *                       Functions::ZeroFunction<dim>(),
 *                       boundary_values);
 *                     break;
 *                   default:
 *                     DEAL_II_NOT_IMPLEMENTED();
 *                 }
 *             }
 *           MatrixTools::apply_boundary_values(boundary_values,
 *                                              vel_it_matrix[d],
 *                                              u_n[d],
 *                                              force[d]);
 *         }
 *   
 *   
 *       Threads::TaskGroup<void> tasks;
 *       for (unsigned int d = 0; d < dim; ++d)
 *         {
 *           if (reinit_prec)
 *             prec_velocity[d].initialize(vel_it_matrix[d],
 *                                         SparseILU<double>::AdditionalData(
 *                                           vel_diag_strength, vel_off_diagonals));
 *           tasks += Threads::new_task(
 *             &NavierStokesProjection<dim>::diffusion_component_solve, *this, d);
 *         }
 *       tasks.join_all();
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void
 *     NavierStokesProjection<dim>::diffusion_component_solve(const unsigned int d)
 *     {
 *       SolverControl solver_control(vel_max_its, vel_eps * force[d].l2_norm());
 *       SolverGMRES<Vector<double>> gmres(
 *         solver_control,
 *         SolverGMRES<Vector<double>>::AdditionalData(vel_Krylov_size));
 *       gmres.solve(vel_it_matrix[d], u_n[d], force[d], prec_velocity[d]);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_35-codeNavierStokesProjectionassemble_advection_termcode"></a> 
 * <h4> <code>NavierStokesProjection::assemble_advection_term</code> </h4>
 * 
 
 * 
 * The following few functions deal with assembling the advection terms, which
 * is the part of the system matrix for the diffusion step that changes at
 * every time step. As mentioned above, we will run the assembly loop over all
 * cells in %parallel, using the WorkStream class and other
 * facilities as described in the documentation topic on @ref threads.
 * 
 * @code
 *     template <int dim>
 *     void NavierStokesProjection<dim>::assemble_advection_term()
 *     {
 *       vel_Advection = 0.;
 *       AdvectionPerTaskData data(fe_velocity.n_dofs_per_cell());
 *       AdvectionScratchData scratch(fe_velocity,
 *                                    quadrature_velocity,
 *                                    update_values | update_JxW_values |
 *                                      update_gradients);
 *       WorkStream::run(
 *         dof_handler_velocity.begin_active(),
 *         dof_handler_velocity.end(),
 *         *this,
 *         &NavierStokesProjection<dim>::assemble_one_cell_of_advection,
 *         &NavierStokesProjection<dim>::copy_advection_local_to_global,
 *         scratch,
 *         data);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void NavierStokesProjection<dim>::assemble_one_cell_of_advection(
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       AdvectionScratchData                                 &scratch,
 *       AdvectionPerTaskData                                 &data)
 *     {
 *       scratch.fe_val.reinit(cell);
 *       cell->get_dof_indices(data.local_dof_indices);
 *       for (unsigned int d = 0; d < dim; ++d)
 *         {
 *           scratch.fe_val.get_function_values(u_star[d], scratch.u_star_tmp);
 *           for (unsigned int q = 0; q < scratch.nqp; ++q)
 *             scratch.u_star_local[q][d] = scratch.u_star_tmp[q];
 *         }
 *   
 *       for (unsigned int d = 0; d < dim; ++d)
 *         {
 *           scratch.fe_val.get_function_gradients(u_star[d], scratch.grad_u_star);
 *           for (unsigned int q = 0; q < scratch.nqp; ++q)
 *             {
 *               if (d == 0)
 *                 scratch.u_star_tmp[q] = 0.;
 *               scratch.u_star_tmp[q] += scratch.grad_u_star[q][d];
 *             }
 *         }
 *   
 *       data.local_advection = 0.;
 *       for (unsigned int q = 0; q < scratch.nqp; ++q)
 *         for (unsigned int i = 0; i < scratch.dpc; ++i)
 *           for (unsigned int j = 0; j < scratch.dpc; ++j)
 *             data.local_advection(i, j) += (scratch.u_star_local[q] *            
 *                                              scratch.fe_val.shape_grad(j, q) *  
 *                                              scratch.fe_val.shape_value(i, q)   
 *                                            +                                    
 *                                            0.5 *                                
 *                                              scratch.u_star_tmp[q] *            
 *                                              scratch.fe_val.shape_value(i, q) * 
 *                                              scratch.fe_val.shape_value(j, q))  
 *                                           * scratch.fe_val.JxW(q);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void NavierStokesProjection<dim>::copy_advection_local_to_global(
 *       const AdvectionPerTaskData &data)
 *     {
 *       for (unsigned int i = 0; i < fe_velocity.n_dofs_per_cell(); ++i)
 *         for (unsigned int j = 0; j < fe_velocity.n_dofs_per_cell(); ++j)
 *           vel_Advection.add(data.local_dof_indices[i],
 *                             data.local_dof_indices[j],
 *                             data.local_advection(i, j));
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_35-codeNavierStokesProjectionprojection_stepcode"></a> 
 * <h4><code>NavierStokesProjection::projection_step</code></h4>
 * 
 
 * 
 * This implements the projection step:
 * 
 * @code
 *     template <int dim>
 *     void NavierStokesProjection<dim>::projection_step(const bool reinit_prec)
 *     {
 *       pres_iterative.copy_from(pres_Laplace);
 *   
 *       pres_tmp = 0.;
 *       for (unsigned d = 0; d < dim; ++d)
 *         pres_Diff[d].Tvmult_add(pres_tmp, u_n[d]);
 *   
 *       phi_n_minus_1 = phi_n;
 *   
 *       static std::map<types::global_dof_index, double> bval;
 *       if (reinit_prec)
 *         VectorTools::interpolate_boundary_values(dof_handler_pressure,
 *                                                  3,
 *                                                  Functions::ZeroFunction<dim>(),
 *                                                  bval);
 *   
 *       MatrixTools::apply_boundary_values(bval, pres_iterative, phi_n, pres_tmp);
 *   
 *       if (reinit_prec)
 *         prec_pres_Laplace.initialize(pres_iterative,
 *                                      SparseILU<double>::AdditionalData(
 *                                        vel_diag_strength, vel_off_diagonals));
 *   
 *       SolverControl solvercontrol(vel_max_its, vel_eps * pres_tmp.l2_norm());
 *       SolverCG<Vector<double>> cg(solvercontrol);
 *       cg.solve(pres_iterative, phi_n, pres_tmp, prec_pres_Laplace);
 *   
 *       phi_n *= 1.5 / dt;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_35-codeNavierStokesProjectionupdate_pressurecode"></a> 
 * <h4> <code>NavierStokesProjection::update_pressure</code> </h4>
 * 
 
 * 
 * This is the pressure update step of the projection method. It implements
 * the standard formulation of the method, that is @f[ p^{n+1} = p^n +
 * \phi^{n+1}, @f] or the rotational form, which is @f[ p^{n+1} = p^n +
 * \phi^{n+1} - \frac{1}{Re} \nabla\cdot u^{n+1}. @f]
 * 
 * @code
 *     template <int dim>
 *     void NavierStokesProjection<dim>::update_pressure(const bool reinit_prec)
 *     {
 *       pres_n_minus_1 = pres_n;
 *       switch (type)
 *         {
 *           case RunTimeParameters::Method::standard:
 *             pres_n += phi_n;
 *             break;
 *           case RunTimeParameters::Method::rotational:
 *             if (reinit_prec)
 *               prec_mass.initialize(pres_Mass);
 *             pres_n = pres_tmp;
 *             prec_mass.solve(pres_n);
 *             pres_n.sadd(1. / Re, 1., pres_n_minus_1);
 *             pres_n += phi_n;
 *             break;
 *           default:
 *             DEAL_II_NOT_IMPLEMENTED();
 *         };
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_35-codeNavierStokesProjectionoutput_resultscode"></a> 
 * <h4> <code>NavierStokesProjection::output_results</code> </h4>
 * 
 
 * 
 * This method plots the current solution. The main difficulty is that we want
 * to create a single output file that contains the data for all velocity
 * components, the pressure, and also the vorticity of the flow. On the other
 * hand, velocities and the pressure live on separate DoFHandler objects, and
 * so can't be written to the same file using a single DataOut object. As a
 * consequence, we have to work a bit harder to get the various pieces of data
 * into a single DoFHandler object, and then use that to drive graphical
 * output.
 *   
 
 * 
 * We will not elaborate on this process here, but rather refer to @ref step_32 "step-32",
 * where a similar procedure is used (and is documented) to create a joint
 * DoFHandler object for all variables.
 *   
 
 * 
 * Let us also note that we here compute the vorticity as a scalar quantity in
 * a separate function, using the @f$L^2@f$ projection of the quantity
 * @f$\text{curl} u@f$ onto the finite element space used for the components of
 * the velocity. In principle, however, we could also have computed it as a
 * pointwise quantity from the velocity, and do so through the
 * DataPostprocessor mechanism discussed in @ref step_29 "step-29" and @ref step_33 "step-33".
 * 
 * @code
 *     template <int dim>
 *     void NavierStokesProjection<dim>::output_results(const unsigned int step)
 *     {
 *       assemble_vorticity((step == 1));
 *       const FESystem<dim> joint_fe(fe_velocity ^ dim, fe_pressure, fe_velocity);
 *       DoFHandler<dim>     joint_dof_handler(triangulation);
 *       joint_dof_handler.distribute_dofs(joint_fe);
 *       Assert(joint_dof_handler.n_dofs() ==
 *                ((dim + 1) * dof_handler_velocity.n_dofs() +
 *                 dof_handler_pressure.n_dofs()),
 *              ExcInternalError());
 *       Vector<double> joint_solution(joint_dof_handler.n_dofs());
 *       std::vector<types::global_dof_index> loc_joint_dof_indices(
 *         joint_fe.n_dofs_per_cell()),
 *         loc_vel_dof_indices(fe_velocity.n_dofs_per_cell()),
 *         loc_pres_dof_indices(fe_pressure.n_dofs_per_cell());
 *       typename DoFHandler<dim>::active_cell_iterator
 *         joint_cell = joint_dof_handler.begin_active(),
 *         joint_endc = joint_dof_handler.end(),
 *         vel_cell   = dof_handler_velocity.begin_active(),
 *         pres_cell  = dof_handler_pressure.begin_active();
 *       for (; joint_cell != joint_endc; ++joint_cell, ++vel_cell, ++pres_cell)
 *         {
 *           joint_cell->get_dof_indices(loc_joint_dof_indices);
 *           vel_cell->get_dof_indices(loc_vel_dof_indices);
 *           pres_cell->get_dof_indices(loc_pres_dof_indices);
 *           for (unsigned int i = 0; i < joint_fe.n_dofs_per_cell(); ++i)
 *             switch (joint_fe.system_to_base_index(i).first.first)
 *               {
 *                 case 0:
 *                   Assert(joint_fe.system_to_base_index(i).first.second < dim,
 *                          ExcInternalError());
 *                   joint_solution(loc_joint_dof_indices[i]) =
 *                     u_n[joint_fe.system_to_base_index(i).first.second](
 *                       loc_vel_dof_indices[joint_fe.system_to_base_index(i)
 *                                             .second]);
 *                   break;
 *                 case 1:
 *                   Assert(joint_fe.system_to_base_index(i).first.second == 0,
 *                          ExcInternalError());
 *                   joint_solution(loc_joint_dof_indices[i]) =
 *                     pres_n(loc_pres_dof_indices[joint_fe.system_to_base_index(i)
 *                                                   .second]);
 *                   break;
 *                 case 2:
 *                   Assert(joint_fe.system_to_base_index(i).first.second == 0,
 *                          ExcInternalError());
 *                   joint_solution(loc_joint_dof_indices[i]) = rot_u(
 *                     loc_vel_dof_indices[joint_fe.system_to_base_index(i).second]);
 *                   break;
 *                 default:
 *                   DEAL_II_ASSERT_UNREACHABLE();
 *               }
 *         }
 *       std::vector<std::string> joint_solution_names(dim, "v");
 *       joint_solution_names.emplace_back("p");
 *       joint_solution_names.emplace_back("rot_u");
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(joint_dof_handler);
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         component_interpretation(
 *           dim + 2, DataComponentInterpretation::component_is_part_of_vector);
 *       component_interpretation[dim] =
 *         DataComponentInterpretation::component_is_scalar;
 *       component_interpretation[dim + 1] =
 *         DataComponentInterpretation::component_is_scalar;
 *       data_out.add_data_vector(joint_solution,
 *                                joint_solution_names,
 *                                DataOut<dim>::type_dof_data,
 *                                component_interpretation);
 *       data_out.build_patches(deg + 1);
 *       std::ofstream output("solution-" + Utilities::int_to_string(step, 5) +
 *                            ".vtk");
 *       data_out.write_vtk(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * Following is the helper function that computes the vorticity by projecting
 * the term @f$\text{curl} u@f$ onto the finite element space used for the
 * components of the velocity. The function is only called whenever we
 * generate graphical output, so not very often, and as a consequence we
 * didn't bother parallelizing it using the WorkStream concept as we do for
 * the other assembly functions. That should not be overly complicated,
 * however, if needed. Moreover, the implementation that we have here only
 * works for 2d, so we bail if that is not the case.
 * 
 * @code
 *     template <int dim>
 *     void NavierStokesProjection<dim>::assemble_vorticity(const bool reinit_prec)
 *     {
 *       Assert(dim == 2, ExcNotImplemented());
 *       if (reinit_prec)
 *         prec_vel_mass.initialize(vel_Mass);
 *   
 *       FEValues<dim>      fe_val_vel(fe_velocity,
 *                                quadrature_velocity,
 *                                update_gradients | update_JxW_values |
 *                                  update_values);
 *       const unsigned int dpc = fe_velocity.n_dofs_per_cell(),
 *                          nqp = quadrature_velocity.size();
 *       std::vector<types::global_dof_index> ldi(dpc);
 *       Vector<double>                       loc_rot(dpc);
 *   
 *       std::vector<Tensor<1, dim>> grad_u1(nqp), grad_u2(nqp);
 *       rot_u = 0.;
 *   
 *       for (const auto &cell : dof_handler_velocity.active_cell_iterators())
 *         {
 *           fe_val_vel.reinit(cell);
 *           cell->get_dof_indices(ldi);
 *           fe_val_vel.get_function_gradients(u_n[0], grad_u1);
 *           fe_val_vel.get_function_gradients(u_n[1], grad_u2);
 *           loc_rot = 0.;
 *           for (unsigned int q = 0; q < nqp; ++q)
 *             for (unsigned int i = 0; i < dpc; ++i)
 *               loc_rot(i) += (grad_u2[q][0] - grad_u1[q][1]) * 
 *                             fe_val_vel.shape_value(i, q) *    
 *                             fe_val_vel.JxW(q);
 *   
 *           for (unsigned int i = 0; i < dpc; ++i)
 *             rot_u(ldi[i]) += loc_rot(i);
 *         }
 *   
 *       prec_vel_mass.solve(rot_u);
 *     }
 *   } // namespace Step35
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_35-Themainfunction"></a> 
 * <h3> The main function </h3>
 * 
 
 * 
 * The main function looks very much like in all the other tutorial programs, so
 * there is little to comment on here:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step35;
 *   
 *         RunTimeParameters::Data_Storage data;
 *         data.read_data("parameter-file.prm");
 *   
 *         deallog.depth_console(data.verbose ? 2 : 0);
 *   
 *         NavierStokesProjection<2> test(data);
 *         test.run(data.verbose, data.output_interval);
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     std::cout << "----------------------------------------------------"
 *               << std::endl
 *               << "Apparently everything went fine!" << std::endl
 *               << "Don't forget to brush your teeth :-)" << std::endl
 *               << std::endl;
 *     return 0;
 *   }
 * @endcode
<a name="step_35-Results"></a><h1>Results</h1>
 
 
<a name="step_35-Re100"></a><h3> Re = 100 </h3>
 
 
We run the code with the following <code>parameter-file.prm</code>, which can be found in the
same directory as the source:
@verbatim
  # First a global definition
  # the type of method we want to use
  set Method_Form = rotational
 
  subsection Physical data
    # In this subsection we declare the physical data
    # The initial and final time, and the Reynolds number
    set initial_time = 0.
    set final_time   = 25.
    set Reynolds     = 100
  end
 
  subsection Time step data
    # In this subsection we declare the data that is to be used for time discretization,
    # i.e. the time step dt
    set dt = 5e-3
  end
 
  subsection Space discretization
    # In this subsection we declare the data that is relevant to the space discretization
    # we set the number of global refines the triangulation must have
    # and the degree k of the pair Q_(k+1)--Q_k of velocity--pressure finite element spaces
    set n_of_refines = 3
    set pressure_fe_degree = 1
  end
 
  subsection Data solve velocity
    # In this section we declare the parameters that are going to control the solution process
    # for the velocity.
    set max_iterations = 1000  # maximal number of iterations that GMRES must make
    set eps            = 1e-6  # stopping criterion
    set Krylov_size    = 30    # size of the Krylov subspace to be used in GMRES
    set off_diagonals  = 60    # number of off diagonals that ILU must compute
    set diag_strength  = 0.01  # diagonal strengthening value
    set update_prec    = 10    # this number indicates how often the preconditioner must be updated
  end
 
  #The output frequency
  set output_interval = 50
 
  #Finally we set the verbosity level
  set verbose = false
@endverbatim
 
Since the <code>verbose</code> parameter is set to <code>false</code>,
we do not get any kind of output besides the number of the time step
the program is currently working on.
If we were to set it to <code>true</code> we would get information on what the program is doing and
how many steps each iterative process had to make to converge, etc.
 
Let us plot the obtained results for @f$t=1,5,12,20,25@f$ (i.e. time steps
200, 1000, 2400, 4000, and 5000), where in the left column we show the
vorticity and in the right the velocity field:
 
<table>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.vorticity.0.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.velocity.0.9.3.png" alt="" width="400"> </td>
  </tr>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.vorticity.1.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.velocity.1.9.3.png" alt="" width="400"> </td>
  </tr>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.vorticity.2.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.velocity.2.9.3.png" alt="" width="400"> </td>
  </tr>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.vorticity.3.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.velocity.3.9.3.png" alt="" width="400"> </td>
  </tr>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.vorticity.4.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_100.velocity.4.9.3.png" alt="" width="400"> </td>
  </tr>
</table>
 
The images show nicely the development and extension of a vortex chain
behind the obstacles, with the sign of the vorticity indicating
whether this is a left or right turning vortex.
 
 
<a name="step_35-Re500"></a><h3> Re = 500 </h3>
 
 
We can change the Reynolds number, @f$Re@f$, in the parameter file to a
value of @f$500@f$. Doing so, and reducing the time step somewhat as well,
yields the following images at times @f$t=20,40@f$:
 
<table>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_500.vorticity.0.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_500.velocity.0.9.3.png" alt="" width="400"> </td>
  </tr>
  <tr>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_500.vorticity.1.9.3.png" alt="" width="400"> </td>
    <td> <img src="https://dealii.org/images/steps/developer/step-35.Re_500.velocity.1.9.3.png" alt="" width="400"> </td>
  </tr>
</table>
 
For this larger Reynolds number, we observe unphysical oscillations, especially
for the vorticity. The discretization scheme has now difficulties in correctly
resolving the flow, which should still be laminar and well-organized.
These phenomena are typical of discretization schemes that lack robustness
in under-resolved scenarios, where under-resolved means that the Reynolds
number computed with the mesh size instead of the physical dimensions of
the geometry is large. We look at a zoom at the region behind the obstacle, and
the mesh size we have there:
 
 
<img src="https://dealii.org/images/steps/developer/step-35.Re_500.zoom.9.3.png" alt="" width="400">
 
We can easily test our hypothesis by re-running the simulation with one more
mesh refinement set in the parameter file :
 
<img src="https://dealii.org/images/steps/developer/step-35.Re_500.zoom_2.9.3.png" alt="" width="400">
 
Indeed, the vorticity field now looks much smoother. While we can expect that
further refining the mesh will suppress the remaining oscillations as well,
one should take measures to obtain a robust scheme in the limit of coarse
resolutions, as described below.
 
 
<a name="step-35-extensions"></a>
<a name="step_35-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
 
 
This program can be extended in the following directions:
<ul>
  <li> Adaptive mesh refinement: As we have seen, we computed everything on a single fixed mesh.
  Using adaptive mesh refinement can lead to increased accuracy while not significantly increasing the
  computational time.
 
  <li> Adaptive time-stepping: Although there apparently is currently no theory about
  projection methods with variable time step,
  practice shows that they perform very well.
 
  <li> High Reynolds %numbers: As we can see from the results, increasing the Reynolds number changes significantly
  the behavior of the discretization scheme. Using well-known stabilization techniques we could be able to
  compute the flow in this, or many other problems, when the Reynolds number is very large and where computational
  costs demand spatial resolutions for which the flow is only marginally resolved, especially for 3D turbulent
  flows.
 
  <li> Variable density incompressible flows: There are projection-like methods for the case of incompressible
  flows with variable density. Such flows play a role if fluids of different
  density mix, for example fresh water and salt water, or alcohol and water.
 
  <li> Compressible Navier-Stokes equations: These equations are relevant for
  cases where
  velocities are high enough so that the fluid becomes compressible, but not
  fast enough that we get into a regime where viscosity becomes negligible
  and the Navier-Stokes equations need to be replaced by the hyperbolic Euler
  equations of gas dynamics. Compressibility starts to become a factor if the
  velocity becomes greater than about one third of the speed of sound, so it
  is not a factor for almost all terrestrial vehicles. On the other hand,
  commercial jetliners fly at about 85 per cent of the speed of sound, and
  flow over the wings becomes significantly supersonic, a regime in which the
  compressible Navier-Stokes equations are not applicable any more
  either. There are significant applications for the range in between,
  however, such as for small aircraft or the fast trains in many European and
  East Asian countries.
</ul>
 *
 *
<a name="step_35-PlainProg"></a>
<h1> The plain program</h1>
@include "step-35.cc"
*/