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
 
 
*/