cdr.h

Go to the documentation of this file.

 true);
 *   SparsityTools::distribute_sparsity_pattern(dynamic_sparsity_pattern,
 *                                              dof_handler.locally_owned_dofs(),
 *                                              mpi_communicator,
 *                                              locally_relevant_dofs);
 * 
 *   system_rhs.reinit(locally_owned_dofs, mpi_communicator);
 *   system_matrix.reinit(locally_owned_dofs,
 *                        dynamic_sparsity_pattern,
 *                        mpi_communicator);
 * 
 *   CDR::create_system_matrix<dim>(dof_handler,
 *                                  quad,
 *                                  convection_function,
 *                                  parameters,
 *                                  time_step,
 *                                  constraints,
 *                                  system_matrix);
 *   system_matrix.compress(VectorOperation::add);
 *   preconditioner.initialize(system_matrix);
 * }
 * 
 * 
 * template <int dim>
 * void
 * CDRProblem<dim>::time_iterate()
 * {
 *   double               current_time = parameters.start_time;
 *   CDR::WritePVTUOutput pvtu_output(parameters.patch_level);
 *   for (unsigned int time_step_n = 0; time_step_n < parameters.n_time_steps;
 *        ++time_step_n)
 *     {
 *       current_time += time_step;
 * 
 *       system_rhs = 0.0;
 *       CDR::create_system_rhs<dim>(dof_handler,
 *                                   quad,
 *                                   convection_function,
 *                                   forcing_function,
 *                                   parameters,
 *                                   locally_relevant_solution,
 *                                   constraints,
 *                                   current_time,
 *                                   system_rhs);
 *       system_rhs.compress(VectorOperation::add);
 * 
 *       SolverControl                 solver_control(dof_handler.n_dofs(),
 *                                    1e-6 * system_rhs.l2_norm(),
 *                                    /*log_history = */ false,
 *                                    /*log_result = */ false);
 *       TrilinosWrappers::SolverGMRES solver(solver_control);
 *       solver.solve(system_matrix,
 *                    completely_distributed_solution,
 *                    system_rhs,
 *                    preconditioner);
 *       constraints.distribute(completely_distributed_solution);
 *       locally_relevant_solution = completely_distributed_solution;
 * 
 *       if (time_step_n % parameters.save_interval == 0)
 *         {
 *           pvtu_output.write_output(dof_handler,
 *                                    locally_relevant_solution,
 *                                    time_step_n,
 *                                    current_time);
 *         }
 * 
 *       refine_mesh();
 *     }
 * }
 * 
 * 
 * template <int dim>
 * void
 * CDRProblem<dim>::refine_mesh()
 * {
 *   using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
 * 
 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   KellyErrorEstimator<dim>::estimate(dof_handler,
 *                                      QGauss<dim - 1>(fe.degree + 1),
 *                                      FunctionMap(),
 *                                      locally_relevant_solution,
 *                                      estimated_error_per_cell);
 * 
 * @endcode
 * 
 * This solver uses a crude refinement strategy where cells with relatively
 * high errors are refined and cells with relatively low errors are
 * coarsened. The maximum refinement level is capped to prevent run-away
 * refinement.
 * 
 * @code
 *   for (const auto &cell : triangulation.active_cell_iterators())
 *     {
 *       if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) >= 1e-3)
 *         {
 *           cell->set_refine_flag();
 *         }
 *       else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <=
 *                1e-5)
 *         {
 *           cell->set_coarsen_flag();
 *         }
 *     }
 * 
 *   if (triangulation.n_levels() > parameters.max_refinement_level)
 *     {
 *       for (const auto &cell : triangulation.cell_iterators_on_level(
 *              parameters.max_refinement_level))
 *         {
 *           cell->clear_refine_flag();
 *         }
 *     }
 * 
 * @endcode
 * 
 * Transferring the solution between different grids is ultimately just a
 * few function calls but they must be made in exactly the right order.
 * 
 * @code
 *   parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
 *     solution_transfer(dof_handler);
 * 
 *   triangulation.prepare_coarsening_and_refinement();
 *   solution_transfer.prepare_for_coarsening_and_refinement(
 *     locally_relevant_solution);
 *   triangulation.execute_coarsening_and_refinement();
 * 
 *   setup_dofs();
 * 
 * @endcode
 * 
 * The <code>solution_transfer</code> object stores a pointer to
 * <code>locally_relevant_solution</code>, so when
 * parallel::distributed::SolutionTransfer::interpolate is called it uses
 * those values to populate <code>temporary</code>.
 * 
 * @code
 *   TrilinosWrappers::MPI::Vector temporary(locally_owned_dofs, mpi_communicator);
 *   solution_transfer.interpolate(temporary);
 * @endcode
 * 
 * After <code>temporary</code> has the correct value, this call correctly
 * populates <code>completely_distributed_solution</code>, which had its
 * index set updated above with the call to <code>setup_dofs</code>.
 * 
 * @code
 *   completely_distributed_solution = temporary;
 * @endcode
 * 
 * Constraints cannot be applied to
 * @ref GlossGhostedVector "vectors with ghost entries" since the ghost
 * entries are write only, so this first goes through the completely
 * distributed vector.
 * 
 * @code
 *   constraints.distribute(completely_distributed_solution);
 *   locally_relevant_solution = completely_distributed_solution;
 *   setup_system();
 * }
 * 
 * 
 * template <int dim>
 * void
 * CDRProblem<dim>::run()
 * {
 *   setup_geometry();
 *   setup_dofs();
 *   setup_system();
 *   time_iterate();
 * }
 * 
 * 
 * constexpr int dim{2};
 * 
 * 
 * int
 * main(int argc, char *argv[])
 * {
 * @endcode
 * 
 * One of the new features in C++11 is the <code>chrono</code> component of
 * the standard library. This gives us an easy way to time the output.
 * 
 * @code
 *   auto t0 = std::chrono::high_resolution_clock::now();
 * 
 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
 *   CDR::Parameters                  parameters;
 *   parameters.read_parameter_file("parameters.prm");
 *   CDRProblem<dim> cdr_problem(parameters);
 *   cdr_problem.run();
 * 
 *   auto t1 = std::chrono::high_resolution_clock::now();
 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 *     {
 *       std::cout << "time elapsed: "
 *                 << std::chrono::duration_cast<std::chrono::milliseconds>(t1 -
 *                                                                          t0)
 *                      .count()
 *                 << " milliseconds." << std::endl;
 *     }
 * 
 *   return 0;
 * }
 * @endcode
 
 
*/