From 18b844182ac3b24d58ff196d36848f270c599a24 Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 25 Jul 2022 11:28:18 -0400 Subject: [PATCH] Run clang-format on everything --- cdr/common/include/deal.II-cdr/parameters.h | 15 +- .../include/deal.II-cdr/system_matrix.h | 38 +-- .../deal.II-cdr/system_matrix.templates.h | 138 +++++---- cdr/common/include/deal.II-cdr/system_rhs.h | 25 +- .../deal.II-cdr/system_rhs.templates.h | 100 ++++--- .../include/deal.II-cdr/write_pvtu_output.h | 14 +- .../deal.II-cdr/write_pvtu_output.templates.h | 42 +-- cdr/common/source/parameters.cc | 111 ++++--- cdr/common/source/system_matrix.cc | 138 ++++----- cdr/common/source/system_rhs.cc | 90 +++--- cdr/common/source/write_pvtu_output.cc | 46 +-- cdr/solver/cdr.cc | 272 ++++++++++-------- 12 files changed, 568 insertions(+), 461 deletions(-) diff --git a/cdr/common/include/deal.II-cdr/parameters.h b/cdr/common/include/deal.II-cdr/parameters.h index 83e21bf..4a88d3b 100644 --- a/cdr/common/include/deal.II-cdr/parameters.h +++ b/cdr/common/include/deal.II-cdr/parameters.h @@ -23,22 +23,25 @@ namespace CDR double diffusion_coefficient; double reaction_coefficient; - bool time_dependent_forcing; + bool time_dependent_forcing; unsigned int initial_refinement_level; unsigned int max_refinement_level; unsigned int fe_order; - double start_time; - double stop_time; + double start_time; + double stop_time; unsigned int n_time_steps; unsigned int save_interval; unsigned int patch_level; - void read_parameter_file(const std::string &file_name); + void + read_parameter_file(const std::string &file_name); + private: - void configure_parameter_handler(ParameterHandler ¶meter_handler); + void + configure_parameter_handler(ParameterHandler ¶meter_handler); }; -} +} // namespace CDR #endif diff --git a/cdr/common/include/deal.II-cdr/system_matrix.h b/cdr/common/include/deal.II-cdr/system_matrix.h index 41fd50b..489497f 100644 --- a/cdr/common/include/deal.II-cdr/system_matrix.h +++ b/cdr/common/include/deal.II-cdr/system_matrix.h @@ -18,23 +18,25 @@ namespace CDR { using namespace dealii; - template - void create_system_matrix - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - MatrixType &system_matrix); + template + void + create_system_matrix( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + MatrixType & system_matrix); - template - void create_system_matrix - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - const AffineConstraints &constraints, - MatrixType &system_matrix); -} + template + void + create_system_matrix( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + const AffineConstraints & constraints, + MatrixType & system_matrix); +} // namespace CDR #endif diff --git a/cdr/common/include/deal.II-cdr/system_matrix.templates.h b/cdr/common/include/deal.II-cdr/system_matrix.templates.h index 12340b5..8e4cabe 100644 --- a/cdr/common/include/deal.II-cdr/system_matrix.templates.h +++ b/cdr/common/include/deal.II-cdr/system_matrix.templates.h @@ -4,11 +4,11 @@ #include -#include - #include #include +#include + #include #include @@ -22,20 +22,23 @@ namespace CDR // This is the actual implementation of the create_system_matrix // function described in the header file. It is similar to the system matrix // assembly routine in step-40. - template - void internal_create_system_matrix - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - UpdateFunction update_system_matrix) + template + void + internal_create_system_matrix( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + UpdateFunction update_system_matrix) { - auto &fe = dof_handler.get_fe(); - const auto dofs_per_cell = fe.dofs_per_cell; + auto & fe = dof_handler.get_fe(); + const auto dofs_per_cell = fe.dofs_per_cell; FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - FEValues fe_values(fe, quad, update_values | update_gradients | - update_quadrature_points | update_JxW_values); + FEValues fe_values(fe, + quad, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); std::vector local_indices(dofs_per_cell); @@ -55,19 +58,23 @@ namespace CDR { for (unsigned int j = 0; j < dofs_per_cell; ++j) { - const auto convection_contribution = current_convection - *fe_values.shape_grad(j, q); - cell_matrix(i, j) += fe_values.JxW(q)* - // Here are the time step, mass, and reaction parts: - ((1.0 + time_step/2.0*parameters.reaction_coefficient) - *fe_values.shape_value(i, q)*fe_values.shape_value(j, q) - + time_step/2.0* - // and the convection part: - (fe_values.shape_value(i, q)*convection_contribution - // and, finally, the diffusion part: - + parameters.diffusion_coefficient - *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q))) - ); + const auto convection_contribution = + current_convection * fe_values.shape_grad(j, q); + cell_matrix(i, j) += + fe_values.JxW(q) * + // Here are the time step, mass, and reaction parts: + ((1.0 + + time_step / 2.0 * parameters.reaction_coefficient) * + fe_values.shape_value(i, q) * + fe_values.shape_value(j, q) + + time_step / 2.0 * + // and the convection part: + (fe_values.shape_value(i, q) * + convection_contribution + // and, finally, the diffusion part: + + parameters.diffusion_coefficient * + (fe_values.shape_grad(i, q) * + fe_values.shape_grad(j, q)))); } } } @@ -76,42 +83,53 @@ namespace CDR } } - template - void create_system_matrix - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - const AffineConstraints &constraints, - MatrixType &system_matrix) + template + void + create_system_matrix( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + const AffineConstraints & constraints, + MatrixType & system_matrix) { - internal_create_system_matrix - (dof_handler, quad, convection_function, parameters, time_step, - [&constraints, &system_matrix](const std::vector &local_indices, - const FullMatrix &cell_matrix) - { - constraints.distribute_local_to_global - (cell_matrix, local_indices, system_matrix); - }); + internal_create_system_matrix( + dof_handler, + quad, + convection_function, + parameters, + time_step, + [&constraints, &system_matrix]( + const std::vector &local_indices, + const FullMatrix & cell_matrix) { + constraints.distribute_local_to_global(cell_matrix, + local_indices, + system_matrix); + }); } - template - void create_system_matrix - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - MatrixType &system_matrix) + template + void + create_system_matrix( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + MatrixType & system_matrix) { - internal_create_system_matrix - (dof_handler, quad, convection_function, parameters, time_step, - [&system_matrix](const std::vector &local_indices, - const FullMatrix &cell_matrix) - { - system_matrix.add(local_indices, cell_matrix); - }); + internal_create_system_matrix( + dof_handler, + quad, + convection_function, + parameters, + time_step, + [&system_matrix]( + const std::vector &local_indices, + const FullMatrix & cell_matrix) { + system_matrix.add(local_indices, cell_matrix); + }); } -} +} // namespace CDR #endif diff --git a/cdr/common/include/deal.II-cdr/system_rhs.h b/cdr/common/include/deal.II-cdr/system_rhs.h index 986b5a1..a21d043 100644 --- a/cdr/common/include/deal.II-cdr/system_rhs.h +++ b/cdr/common/include/deal.II-cdr/system_rhs.h @@ -18,16 +18,17 @@ namespace CDR { using namespace dealii; - template - void create_system_rhs - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const std::function)> &forcing_function, - const CDR::Parameters ¶meters, - const VectorType &previous_solution, - const AffineConstraints &constraints, - const double current_time, - VectorType &system_rhs); -} + template + void + create_system_rhs( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const std::function)> &forcing_function, + const CDR::Parameters & parameters, + const VectorType & previous_solution, + const AffineConstraints & constraints, + const double current_time, + VectorType & system_rhs); +} // namespace CDR #endif diff --git a/cdr/common/include/deal.II-cdr/system_rhs.templates.h b/cdr/common/include/deal.II-cdr/system_rhs.templates.h index e7e2ae6..34f1378 100644 --- a/cdr/common/include/deal.II-cdr/system_rhs.templates.h +++ b/cdr/common/include/deal.II-cdr/system_rhs.templates.h @@ -22,32 +22,35 @@ namespace CDR { using namespace dealii; - template - void create_system_rhs - (const DoFHandler &dof_handler, - const QGauss &quad, - const std::function(const Point)> &convection_function, - const std::function)> &forcing_function, - const CDR::Parameters ¶meters, - const VectorType &previous_solution, - const AffineConstraints &constraints, - const double current_time, - VectorType &system_rhs) + template + void + create_system_rhs( + const DoFHandler & dof_handler, + const QGauss & quad, + const std::function(const Point)> &convection_function, + const std::function)> &forcing_function, + const CDR::Parameters & parameters, + const VectorType & previous_solution, + const AffineConstraints & constraints, + const double current_time, + VectorType & system_rhs) { - auto &fe = dof_handler.get_fe(); - const auto dofs_per_cell = fe.dofs_per_cell; - const double time_step = (parameters.stop_time - parameters.start_time) - /parameters.n_time_steps; - FEValues fe_values(fe, quad, update_values | update_gradients | - update_quadrature_points | update_JxW_values); + auto & fe = dof_handler.get_fe(); + const auto dofs_per_cell = fe.dofs_per_cell; + const double time_step = + (parameters.stop_time - parameters.start_time) / parameters.n_time_steps; + FEValues fe_values(fe, + quad, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); - Vector cell_rhs(dofs_per_cell); + Vector cell_rhs(dofs_per_cell); FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Vector current_fe_coefficients(dofs_per_cell); + Vector current_fe_coefficients(dofs_per_cell); std::vector local_indices(dofs_per_cell); - const double previous_time {current_time - time_step}; + const double previous_time{current_time - time_step}; for (const auto &cell : dof_handler.active_cell_iterators()) { @@ -58,7 +61,8 @@ namespace CDR cell->get_dof_indices(local_indices); for (unsigned int i = 0; i < dofs_per_cell; ++i) { - current_fe_coefficients[i] = previous_solution[local_indices[i]]; + current_fe_coefficients[i] = + previous_solution[local_indices[i]]; } for (unsigned int q = 0; q < quad.size(); ++q) @@ -66,38 +70,46 @@ namespace CDR const auto current_convection = convection_function(fe_values.quadrature_point(q)); - const double current_forcing = forcing_function - (current_time, fe_values.quadrature_point(q)); - const double previous_forcing = forcing_function - (previous_time, fe_values.quadrature_point(q)); + const double current_forcing = + forcing_function(current_time, fe_values.quadrature_point(q)); + const double previous_forcing = + forcing_function(previous_time, + fe_values.quadrature_point(q)); for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { - const auto convection_contribution = current_convection - *fe_values.shape_grad(j, q); + const auto convection_contribution = + current_convection * fe_values.shape_grad(j, q); - cell_rhs(i) += fe_values.JxW(q)* - // Here are the mass and reaction part: - (((1.0 - time_step/2.0*parameters.reaction_coefficient) - *fe_values.shape_value(i, q)*fe_values.shape_value(j, q) - - time_step/2.0* - // the convection part: - (fe_values.shape_value(i, q)*convection_contribution - // the diffusion part: - + parameters.diffusion_coefficient - *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q)))) - *current_fe_coefficients[j] - // and, finally, the forcing function part: - + time_step/2.0* - (current_forcing + previous_forcing) - *fe_values.shape_value(i, q)); + cell_rhs(i) += + fe_values.JxW(q) * + // Here are the mass and reaction part: + (((1.0 - time_step / 2.0 * + parameters.reaction_coefficient) * + fe_values.shape_value(i, q) * + fe_values.shape_value(j, q) - + time_step / 2.0 * + // the convection part: + (fe_values.shape_value(i, q) * + convection_contribution + // the diffusion part: + + parameters.diffusion_coefficient * + (fe_values.shape_grad(i, q) * + fe_values.shape_grad(j, q)))) * + current_fe_coefficients[j] + // and, finally, the forcing function part: + + time_step / 2.0 * + (current_forcing + previous_forcing) * + fe_values.shape_value(i, q)); } } } - constraints.distribute_local_to_global(cell_rhs, local_indices, system_rhs); + constraints.distribute_local_to_global(cell_rhs, + local_indices, + system_rhs); } } } -} +} // namespace CDR #endif diff --git a/cdr/common/include/deal.II-cdr/write_pvtu_output.h b/cdr/common/include/deal.II-cdr/write_pvtu_output.h index 5b0b1c2..9530d88 100644 --- a/cdr/common/include/deal.II-cdr/write_pvtu_output.h +++ b/cdr/common/include/deal.II-cdr/write_pvtu_output.h @@ -12,14 +12,16 @@ namespace CDR public: WritePVTUOutput(const unsigned int patch_level); - template - void write_output(const DoFHandler &dof_handler, - const VectorType &solution, - const unsigned int time_step_n, - const double current_time); + template + void + write_output(const DoFHandler &dof_handler, + const VectorType & solution, + const unsigned int time_step_n, + const double current_time); + private: const unsigned int patch_level; const unsigned int this_mpi_process; }; -} +} // namespace CDR #endif diff --git a/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h b/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h index 15e5154..a0af9ba 100644 --- a/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h +++ b/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h @@ -11,8 +11,8 @@ #include -#include #include +#include #include // Here is the implementation of the important function. This is similar to @@ -21,18 +21,19 @@ namespace CDR { using namespace dealii; - template - void WritePVTUOutput::write_output(const DoFHandler &dof_handler, - const VectorType &solution, - const unsigned int time_step_n, - const double current_time) + template + void + WritePVTUOutput::write_output(const DoFHandler &dof_handler, + const VectorType & solution, + const unsigned int time_step_n, + const double current_time) { DataOut data_out; data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(solution, "u"); - const auto &triangulation = dof_handler.get_triangulation(); - Vector subdomain (triangulation.n_active_cells()); + const auto & triangulation = dof_handler.get_triangulation(); + Vector subdomain(triangulation.n_active_cells()); for (auto &domain : subdomain) { domain = triangulation.locally_owned_subdomain(); @@ -44,7 +45,8 @@ namespace CDR flags.time = current_time; // While the default flag is for the best compression level, using // best_speed makes this function much faster. - flags.compression_level = DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; + flags.compression_level = + DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; data_out.set_flags(flags); unsigned int subdomain_n; @@ -57,25 +59,25 @@ namespace CDR subdomain_n = triangulation.locally_owned_subdomain(); } - std::ofstream output - ("solution-" + Utilities::int_to_string(time_step_n) + "." - + Utilities::int_to_string(subdomain_n, 4) - + ".vtu"); + std::ofstream output("solution-" + Utilities::int_to_string(time_step_n) + + "." + Utilities::int_to_string(subdomain_n, 4) + + ".vtu"); data_out.write_vtu(output); if (this_mpi_process == 0) { std::vector filenames; - for (unsigned int i = 0; i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i) - filenames.push_back - ("solution-" + Utilities::int_to_string (time_step_n) + "." - + Utilities::int_to_string (i, 4) + ".vtu"); - std::ofstream master_output - ("solution-" + Utilities::int_to_string(time_step_n) + ".pvtu"); + filenames.push_back("solution-" + + Utilities::int_to_string(time_step_n) + "." + + Utilities::int_to_string(i, 4) + ".vtu"); + std::ofstream master_output( + "solution-" + Utilities::int_to_string(time_step_n) + ".pvtu"); data_out.write_pvtu_record(master_output, filenames); } } -} +} // namespace CDR #endif diff --git a/cdr/common/source/parameters.cc b/cdr/common/source/parameters.cc index daefef0..16be678 100644 --- a/cdr/common/source/parameters.cc +++ b/cdr/common/source/parameters.cc @@ -5,64 +5,90 @@ namespace CDR { - void Parameters::configure_parameter_handler(ParameterHandler ¶meter_handler) + void + Parameters::configure_parameter_handler(ParameterHandler ¶meter_handler) { parameter_handler.enter_subsection("Geometry"); { - parameter_handler.declare_entry - ("inner_radius", "1.0", Patterns::Double(0.0), "Inner radius."); - parameter_handler.declare_entry - ("outer_radius", "2.0", Patterns::Double(0.0), "Outer radius."); + parameter_handler.declare_entry("inner_radius", + "1.0", + Patterns::Double(0.0), + "Inner radius."); + parameter_handler.declare_entry("outer_radius", + "2.0", + Patterns::Double(0.0), + "Outer radius."); } parameter_handler.leave_subsection(); parameter_handler.enter_subsection("Physical Parameters"); { - parameter_handler.declare_entry - ("diffusion_coefficient", "1.0", Patterns::Double(0.0), "Diffusion coefficient."); - parameter_handler.declare_entry - ("reaction_coefficient", "1.0", Patterns::Double(0.0), "Reaction coefficient."); - parameter_handler.declare_entry - ("time_dependent_forcing", "true", Patterns::Bool(), "Whether or not " - "the forcing function depends on time."); + parameter_handler.declare_entry("diffusion_coefficient", + "1.0", + Patterns::Double(0.0), + "Diffusion coefficient."); + parameter_handler.declare_entry("reaction_coefficient", + "1.0", + Patterns::Double(0.0), + "Reaction coefficient."); + parameter_handler.declare_entry("time_dependent_forcing", + "true", + Patterns::Bool(), + "Whether or not " + "the forcing function depends on time."); } parameter_handler.leave_subsection(); parameter_handler.enter_subsection("Finite Element"); { - parameter_handler.declare_entry - ("initial_refinement_level", "1", Patterns::Integer(1), - "Initial number of levels in the mesh."); - parameter_handler.declare_entry - ("max_refinement_level", "1", Patterns::Integer(1), - "Maximum number of levels in the mesh."); - parameter_handler.declare_entry - ("fe_order", "1", Patterns::Integer(1), "Finite element order."); + parameter_handler.declare_entry("initial_refinement_level", + "1", + Patterns::Integer(1), + "Initial number of levels in the mesh."); + parameter_handler.declare_entry("max_refinement_level", + "1", + Patterns::Integer(1), + "Maximum number of levels in the mesh."); + parameter_handler.declare_entry("fe_order", + "1", + Patterns::Integer(1), + "Finite element order."); } parameter_handler.leave_subsection(); parameter_handler.enter_subsection("Time Step"); { - parameter_handler.declare_entry - ("start_time", "0.0", Patterns::Double(0.0), "Start time."); - parameter_handler.declare_entry - ("stop_time", "1.0", Patterns::Double(1.0), "Stop time."); - parameter_handler.declare_entry - ("n_time_steps", "1", Patterns::Integer(1), "Number of time steps."); + parameter_handler.declare_entry("start_time", + "0.0", + Patterns::Double(0.0), + "Start time."); + parameter_handler.declare_entry("stop_time", + "1.0", + Patterns::Double(1.0), + "Stop time."); + parameter_handler.declare_entry("n_time_steps", + "1", + Patterns::Integer(1), + "Number of time steps."); } parameter_handler.leave_subsection(); parameter_handler.enter_subsection("Output"); { - parameter_handler.declare_entry - ("save_interval", "10", Patterns::Integer(1), "Save interval."); - parameter_handler.declare_entry - ("patch_level", "2", Patterns::Integer(0), "Patch level."); + parameter_handler.declare_entry("save_interval", + "10", + Patterns::Integer(1), + "Save interval."); + parameter_handler.declare_entry("patch_level", + "2", + Patterns::Integer(0), + "Patch level."); } parameter_handler.leave_subsection(); } - void Parameters::read_parameter_file(const std::string &file_name) + void + Parameters::read_parameter_file(const std::string &file_name) { ParameterHandler parameter_handler; { @@ -80,25 +106,30 @@ namespace CDR parameter_handler.enter_subsection("Physical Parameters"); { - diffusion_coefficient = parameter_handler.get_double("diffusion_coefficient"); - reaction_coefficient = parameter_handler.get_double("reaction_coefficient"); - time_dependent_forcing = parameter_handler.get_bool("time_dependent_forcing"); + diffusion_coefficient = + parameter_handler.get_double("diffusion_coefficient"); + reaction_coefficient = + parameter_handler.get_double("reaction_coefficient"); + time_dependent_forcing = + parameter_handler.get_bool("time_dependent_forcing"); } parameter_handler.leave_subsection(); parameter_handler.enter_subsection("Finite Element"); { - initial_refinement_level = parameter_handler.get_integer("initial_refinement_level"); - max_refinement_level = parameter_handler.get_integer("max_refinement_level"); + initial_refinement_level = + parameter_handler.get_integer("initial_refinement_level"); + max_refinement_level = + parameter_handler.get_integer("max_refinement_level"); fe_order = parameter_handler.get_integer("fe_order"); } parameter_handler.leave_subsection(); parameter_handler.enter_subsection("Time Step"); { - start_time = parameter_handler.get_double("start_time"); - stop_time = parameter_handler.get_double("stop_time"); + start_time = parameter_handler.get_double("start_time"); + stop_time = parameter_handler.get_double("stop_time"); n_time_steps = parameter_handler.get_integer("n_time_steps"); } parameter_handler.leave_subsection(); @@ -106,8 +137,8 @@ namespace CDR parameter_handler.enter_subsection("Output"); { save_interval = parameter_handler.get_integer("save_interval"); - patch_level = parameter_handler.get_integer("patch_level"); + patch_level = parameter_handler.get_integer("patch_level"); } parameter_handler.leave_subsection(); } -} +} // namespace CDR diff --git a/cdr/common/source/system_matrix.cc b/cdr/common/source/system_matrix.cc index 4e18902..bc03457 100644 --- a/cdr/common/source/system_matrix.cc +++ b/cdr/common/source/system_matrix.cc @@ -13,79 +13,79 @@ namespace CDR { using namespace dealii; - template - void create_system_matrix<2, SparseMatrix> - (const DoFHandler<2> &dof_handler, - const QGauss<2> &quad, - const std::function(const Point<2>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - SparseMatrix &system_matrix); + template void + create_system_matrix<2, SparseMatrix>( + const DoFHandler<2> & dof_handler, + const QGauss<2> & quad, + const std::function(const Point<2>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + SparseMatrix & system_matrix); - template - void create_system_matrix<3, SparseMatrix> - (const DoFHandler<3> &dof_handler, - const QGauss<3> &quad, - const std::function(const Point<3>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - SparseMatrix &system_matrix); + template void + create_system_matrix<3, SparseMatrix>( + const DoFHandler<3> & dof_handler, + const QGauss<3> & quad, + const std::function(const Point<3>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + SparseMatrix & system_matrix); - template - void create_system_matrix<2, SparseMatrix> - (const DoFHandler<2> &dof_handler, - const QGauss<2> &quad, - const std::function(const Point<2>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - const AffineConstraints &constraints, - SparseMatrix &system_matrix); + template void + create_system_matrix<2, SparseMatrix>( + const DoFHandler<2> & dof_handler, + const QGauss<2> & quad, + const std::function(const Point<2>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + const AffineConstraints & constraints, + SparseMatrix & system_matrix); - template - void create_system_matrix<3, SparseMatrix> - (const DoFHandler<3> &dof_handler, - const QGauss<3> &quad, - const std::function(const Point<3>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - const AffineConstraints &constraints, - SparseMatrix &system_matrix); + template void + create_system_matrix<3, SparseMatrix>( + const DoFHandler<3> & dof_handler, + const QGauss<3> & quad, + const std::function(const Point<3>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + const AffineConstraints & constraints, + SparseMatrix & system_matrix); - template - void create_system_matrix<2, TrilinosWrappers::SparseMatrix> - (const DoFHandler<2> &dof_handler, - const QGauss<2> &quad, - const std::function(const Point<2>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - TrilinosWrappers::SparseMatrix &system_matrix); + template void + create_system_matrix<2, TrilinosWrappers::SparseMatrix>( + const DoFHandler<2> & dof_handler, + const QGauss<2> & quad, + const std::function(const Point<2>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + TrilinosWrappers::SparseMatrix & system_matrix); - template - void create_system_matrix<3, TrilinosWrappers::SparseMatrix> - (const DoFHandler<3> &dof_handler, - const QGauss<3> &quad, - const std::function(const Point<3>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - TrilinosWrappers::SparseMatrix &system_matrix); + template void + create_system_matrix<3, TrilinosWrappers::SparseMatrix>( + const DoFHandler<3> & dof_handler, + const QGauss<3> & quad, + const std::function(const Point<3>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + TrilinosWrappers::SparseMatrix & system_matrix); - template - void create_system_matrix<2, TrilinosWrappers::SparseMatrix> - (const DoFHandler<2> &dof_handler, - const QGauss<2> &quad, - const std::function(const Point<2>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - const AffineConstraints &constraints, - TrilinosWrappers::SparseMatrix &system_matrix); + template void + create_system_matrix<2, TrilinosWrappers::SparseMatrix>( + const DoFHandler<2> & dof_handler, + const QGauss<2> & quad, + const std::function(const Point<2>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + const AffineConstraints & constraints, + TrilinosWrappers::SparseMatrix & system_matrix); - template - void create_system_matrix<3, TrilinosWrappers::SparseMatrix> - (const DoFHandler<3> &dof_handler, - const QGauss<3> &quad, - const std::function(const Point<3>)> &convection_function, - const CDR::Parameters ¶meters, - const double time_step, - const AffineConstraints &constraints, - TrilinosWrappers::SparseMatrix &system_matrix); -} + template void + create_system_matrix<3, TrilinosWrappers::SparseMatrix>( + const DoFHandler<3> & dof_handler, + const QGauss<3> & quad, + const std::function(const Point<3>)> &convection_function, + const CDR::Parameters & parameters, + const double time_step, + const AffineConstraints & constraints, + TrilinosWrappers::SparseMatrix & system_matrix); +} // namespace CDR diff --git a/cdr/common/source/system_rhs.cc b/cdr/common/source/system_rhs.cc index df70485..d70e86e 100644 --- a/cdr/common/source/system_rhs.cc +++ b/cdr/common/source/system_rhs.cc @@ -11,51 +11,51 @@ namespace CDR { using namespace dealii; - template - void create_system_rhs<2, Vector> - (const DoFHandler<2> &dof_handler, - const QGauss<2> &quad, - const std::function(const Point<2>)> &convection_function, - const std::function)> &forcing_function, - const CDR::Parameters ¶meters, - const Vector &previous_solution, - const AffineConstraints &constraints, - const double current_time, - Vector &system_rhs); + template void + create_system_rhs<2, Vector>( + const DoFHandler<2> & dof_handler, + const QGauss<2> & quad, + const std::function(const Point<2>)> & convection_function, + const std::function)> &forcing_function, + const CDR::Parameters & parameters, + const Vector & previous_solution, + const AffineConstraints & constraints, + const double current_time, + Vector & system_rhs); - template - void create_system_rhs<3, Vector> - (const DoFHandler<3> &dof_handler, - const QGauss<3> &quad, - const std::function(const Point<3>)> &convection_function, - const std::function)> &forcing_function, - const CDR::Parameters ¶meters, - const Vector &previous_solution, - const AffineConstraints &constraints, - const double current_time, - Vector &system_rhs); + template void + create_system_rhs<3, Vector>( + const DoFHandler<3> & dof_handler, + const QGauss<3> & quad, + const std::function(const Point<3>)> & convection_function, + const std::function)> &forcing_function, + const CDR::Parameters & parameters, + const Vector & previous_solution, + const AffineConstraints & constraints, + const double current_time, + Vector & system_rhs); - template - void create_system_rhs<2, TrilinosWrappers::MPI::Vector> - (const DoFHandler<2> &dof_handler, - const QGauss<2> &quad, - const std::function(const Point<2>)> &convection_function, - const std::function)> &forcing_function, - const CDR::Parameters ¶meters, - const TrilinosWrappers::MPI::Vector &previous_solution, - const AffineConstraints &constraints, - const double current_time, - TrilinosWrappers::MPI::Vector &system_rhs); + template void + create_system_rhs<2, TrilinosWrappers::MPI::Vector>( + const DoFHandler<2> & dof_handler, + const QGauss<2> & quad, + const std::function(const Point<2>)> & convection_function, + const std::function)> &forcing_function, + const CDR::Parameters & parameters, + const TrilinosWrappers::MPI::Vector & previous_solution, + const AffineConstraints & constraints, + const double current_time, + TrilinosWrappers::MPI::Vector & system_rhs); - template - void create_system_rhs<3, TrilinosWrappers::MPI::Vector> - (const DoFHandler<3> &dof_handler, - const QGauss<3> &quad, - const std::function(const Point<3>)> &convection_function, - const std::function)> &forcing_function, - const CDR::Parameters ¶meters, - const TrilinosWrappers::MPI::Vector &previous_solution, - const AffineConstraints &constraints, - const double current_time, - TrilinosWrappers::MPI::Vector &system_rhs); -} + template void + create_system_rhs<3, TrilinosWrappers::MPI::Vector>( + const DoFHandler<3> & dof_handler, + const QGauss<3> & quad, + const std::function(const Point<3>)> & convection_function, + const std::function)> &forcing_function, + const CDR::Parameters & parameters, + const TrilinosWrappers::MPI::Vector & previous_solution, + const AffineConstraints & constraints, + const double current_time, + TrilinosWrappers::MPI::Vector & system_rhs); +} // namespace CDR diff --git a/cdr/common/source/write_pvtu_output.cc b/cdr/common/source/write_pvtu_output.cc index f0efe08..b918ebb 100644 --- a/cdr/common/source/write_pvtu_output.cc +++ b/cdr/common/source/write_pvtu_output.cc @@ -10,31 +10,31 @@ namespace CDR using namespace dealii; WritePVTUOutput::WritePVTUOutput(const unsigned int patch_level) - : patch_level {patch_level}, - this_mpi_process {Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)} + : patch_level{patch_level} + , this_mpi_process{Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)} {} - template - void WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler, - const Vector &solution, - const unsigned int time_step_n, - const double current_time); + template void + WritePVTUOutput::write_output(const DoFHandler<2> & dof_handler, + const Vector &solution, + const unsigned int time_step_n, + const double current_time); - template - void WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler, - const Vector &solution, - const unsigned int time_step_n, - const double current_time); + template void + WritePVTUOutput::write_output(const DoFHandler<3> & dof_handler, + const Vector &solution, + const unsigned int time_step_n, + const double current_time); - template - void WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler, - const TrilinosWrappers::MPI::Vector &solution, - const unsigned int time_step_n, - const double current_time); + template void + WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler, + const TrilinosWrappers::MPI::Vector &solution, + const unsigned int time_step_n, + const double current_time); - template - void WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler, - const TrilinosWrappers::MPI::Vector &solution, - const unsigned int time_step_n, - const double current_time); -} + template void + WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler, + const TrilinosWrappers::MPI::Vector &solution, + const unsigned int time_step_n, + const double current_time); +} // namespace CDR diff --git a/cdr/solver/cdr.cc b/cdr/solver/cdr.cc index f5b35d3..4b62f1e 100644 --- a/cdr/solver/cdr.cc +++ b/cdr/solver/cdr.cc @@ -9,59 +9,62 @@ #include #include -#include #include +#include #include // These headers are for distributed computations: -#include #include -#include +#include + #include #include -#include +#include +#include +#include #include #include -#include #include -#include -#include -#include - +#include #include #include -#include #include +#include +#include +#include + using namespace dealii; -constexpr int manifold_id {0}; +constexpr int manifold_id{0}; // This is the actual solver class which performs time iteration and calls the // appropriate library functions to do it. -template +template class CDRProblem { public: CDRProblem(const CDR::Parameters ¶meters); - void run(); + void + run(); + private: const CDR::Parameters parameters; - const double time_step; - double current_time; + const double time_step; + double current_time; - MPI_Comm mpi_communicator; + MPI_Comm mpi_communicator; const unsigned int n_mpi_processes; const unsigned int this_mpi_process; - FE_Q fe; - QGauss quad; - const SphericalManifold boundary_description; + FE_Q fe; + QGauss quad; + const SphericalManifold boundary_description; parallel::distributed::Triangulation triangulation; - DoFHandler dof_handler; + DoFHandler dof_handler; const std::function(const Point)> convection_function; const std::function)> forcing_function; @@ -70,7 +73,7 @@ private: IndexSet locally_relevant_dofs; AffineConstraints constraints; - bool first_run; + bool first_run; // As is usual in parallel programs, I keep two copies of parts of the // complete solution: locally_relevant_solution contains both @@ -78,64 +81,70 @@ private: // boundary (the @ref GlossGhostCell "ghost cells") while // completely_distributed_solution only contains the parts of // the solution computed on the current @ref GlossMPIProcess "MPI process". - TrilinosWrappers::MPI::Vector locally_relevant_solution; - TrilinosWrappers::MPI::Vector completely_distributed_solution; - TrilinosWrappers::MPI::Vector system_rhs; - TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector completely_distributed_solution; + TrilinosWrappers::MPI::Vector system_rhs; + TrilinosWrappers::SparseMatrix system_matrix; TrilinosWrappers::PreconditionAMG preconditioner; ConditionalOStream pcout; - void setup_geometry(); - void setup_system(); - void setup_dofs(); - void refine_mesh(); - void time_iterate(); + void + setup_geometry(); + void + setup_system(); + void + setup_dofs(); + void + refine_mesh(); + void + time_iterate(); }; -template -CDRProblem::CDRProblem(const CDR::Parameters ¶meters) : - parameters(parameters), - time_step {(parameters.stop_time - parameters.start_time) - /parameters.n_time_steps -}, -current_time {parameters.start_time}, -mpi_communicator (MPI_COMM_WORLD), -n_mpi_processes {Utilities::MPI::n_mpi_processes(mpi_communicator)}, -this_mpi_process {Utilities::MPI::this_mpi_process(mpi_communicator)}, -fe(parameters.fe_order), -quad(parameters.fe_order + 2), -boundary_description(Point()), -triangulation(mpi_communicator, typename Triangulation::MeshSmoothing - (Triangulation::smoothing_on_refinement | - Triangulation::smoothing_on_coarsening)), -dof_handler(triangulation), -convection_function -{ - [](const Point p) -> Tensor<1, dim> - {Tensor<1, dim> v; v[0] = -p[1]; v[1] = p[0]; return v;} -}, -forcing_function -{ - [](double t, const Point p) -> double - { - return std::exp(-8*t)*std::exp(-40*Utilities::fixed_power<6>(p[0] - 1.5)) - *std::exp(-40*Utilities::fixed_power<6>(p[1])); - } -}, -first_run {true}, -pcout (std::cout, this_mpi_process == 0) +template +CDRProblem::CDRProblem(const CDR::Parameters ¶meters) + : parameters(parameters) + , time_step{(parameters.stop_time - parameters.start_time) / + parameters.n_time_steps} + , current_time{parameters.start_time} + , mpi_communicator(MPI_COMM_WORLD) + , n_mpi_processes{Utilities::MPI::n_mpi_processes(mpi_communicator)} + , this_mpi_process{Utilities::MPI::this_mpi_process(mpi_communicator)} + , fe(parameters.fe_order) + , quad(parameters.fe_order + 2) + , boundary_description(Point()) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , dof_handler(triangulation) + , convection_function{[](const Point p) -> Tensor<1, dim> { + Tensor<1, dim> v; + v[0] = -p[1]; + v[1] = p[0]; + return v; + }} + , forcing_function{[](double t, const Point p) -> double { + return std::exp(-8 * t) * + std::exp(-40 * Utilities::fixed_power<6>(p[0] - 1.5)) * + std::exp(-40 * Utilities::fixed_power<6>(p[1])); + }} + , first_run{true} + , pcout(std::cout, this_mpi_process == 0) { Assert(dim == 2, ExcNotImplemented()); } -template -void CDRProblem::setup_geometry() +template +void +CDRProblem::setup_geometry() { const Point center; - GridGenerator::hyper_shell(triangulation, center, parameters.inner_radius, + GridGenerator::hyper_shell(triangulation, + center, + parameters.inner_radius, parameters.outer_radius); triangulation.set_manifold(manifold_id, boundary_description); for (const auto &cell : triangulation.active_cell_iterators()) @@ -146,12 +155,12 @@ void CDRProblem::setup_geometry() } -template -void CDRProblem::setup_dofs() +template +void +CDRProblem::setup_dofs() { dof_handler.distribute_dofs(fe); - pcout << "Number of degrees of freedom: " - << dof_handler.n_dofs() + pcout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; locally_owned_dofs = dof_handler.locally_owned_dofs(); DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); @@ -159,43 +168,55 @@ void CDRProblem::setup_dofs() constraints.clear(); constraints.reinit(locally_relevant_dofs); DoFTools::make_hanging_node_constraints(dof_handler, constraints); - DoFTools::make_zero_boundary_constraints(dof_handler, manifold_id, constraints); + DoFTools::make_zero_boundary_constraints(dof_handler, + manifold_id, + constraints); constraints.close(); - completely_distributed_solution.reinit - (locally_owned_dofs, mpi_communicator); + completely_distributed_solution.reinit(locally_owned_dofs, mpi_communicator); - locally_relevant_solution.reinit(locally_owned_dofs, locally_relevant_dofs, + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, mpi_communicator); } -template -void CDRProblem::setup_system() +template +void +CDRProblem::setup_system() { DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern, - constraints, /*keep_constrained_dofs*/true); - SparsityTools::distribute_sparsity_pattern - (dynamic_sparsity_pattern, dof_handler.locally_owned_dofs(), - mpi_communicator, locally_relevant_dofs); + DoFTools::make_sparsity_pattern(dof_handler, + dynamic_sparsity_pattern, + constraints, + /*keep_constrained_dofs*/ 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, + system_matrix.reinit(locally_owned_dofs, + dynamic_sparsity_pattern, mpi_communicator); - CDR::create_system_matrix - (dof_handler, quad, convection_function, parameters, time_step, constraints, - system_matrix); + CDR::create_system_matrix(dof_handler, + quad, + convection_function, + parameters, + time_step, + constraints, + system_matrix); system_matrix.compress(VectorOperation::add); preconditioner.initialize(system_matrix); } -template -void CDRProblem::time_iterate() +template +void +CDRProblem::time_iterate() { - double current_time = parameters.start_time; + 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) @@ -203,25 +224,35 @@ void CDRProblem::time_iterate() current_time += time_step; system_rhs = 0.0; - CDR::create_system_rhs - (dof_handler, quad, convection_function, forcing_function, parameters, - locally_relevant_solution, constraints, current_time, system_rhs); + CDR::create_system_rhs(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(), + 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, + 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); + pvtu_output.write_output(dof_handler, + locally_relevant_solution, + time_step_n, + current_time); } refine_mesh(); @@ -229,16 +260,18 @@ void CDRProblem::time_iterate() } -template -void CDRProblem::refine_mesh() +template +void +CDRProblem::refine_mesh() { - using FunctionMap = - std::map *>; + using FunctionMap = std::map *>; Vector estimated_error_per_cell(triangulation.n_active_cells()); - KellyErrorEstimator::estimate - (dof_handler, QGauss(fe.degree + 1), FunctionMap(), - locally_relevant_solution, estimated_error_per_cell); + KellyErrorEstimator::estimate(dof_handler, + QGauss(fe.degree + 1), + FunctionMap(), + locally_relevant_solution, + estimated_error_per_cell); // This solver uses a crude refinement strategy where cells with relatively // high errors are refined and cells with relatively low errors are @@ -250,7 +283,8 @@ void CDRProblem::refine_mesh() { cell->set_refine_flag(); } - else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <= 1e-5) + else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <= + 1e-5) { cell->set_coarsen_flag(); } @@ -258,8 +292,8 @@ void CDRProblem::refine_mesh() if (triangulation.n_levels() > parameters.max_refinement_level) { - for (const auto &cell : - triangulation.cell_iterators_on_level(parameters.max_refinement_level)) + for (const auto &cell : triangulation.cell_iterators_on_level( + parameters.max_refinement_level)) { cell->clear_refine_flag(); } @@ -268,11 +302,11 @@ void CDRProblem::refine_mesh() // Transferring the solution between different grids is ultimately just a // few function calls but they must be made in exactly the right order. parallel::distributed::SolutionTransfer - solution_transfer(dof_handler); + solution_transfer(dof_handler); triangulation.prepare_coarsening_and_refinement(); - solution_transfer.prepare_for_coarsening_and_refinement - (locally_relevant_solution); + solution_transfer.prepare_for_coarsening_and_refinement( + locally_relevant_solution); triangulation.execute_coarsening_and_refinement(); setup_dofs(); @@ -281,8 +315,7 @@ void CDRProblem::refine_mesh() // locally_relevant_solution, so when // parallel::distributed::SolutionTransfer::interpolate is called it uses // those values to populate temporary. - TrilinosWrappers::MPI::Vector temporary - (locally_owned_dofs, mpi_communicator); + TrilinosWrappers::MPI::Vector temporary(locally_owned_dofs, mpi_communicator); solution_transfer.interpolate(temporary); // After temporary has the correct value, this call correctly // populates completely_distributed_solution, which had its @@ -298,8 +331,9 @@ void CDRProblem::refine_mesh() } -template -void CDRProblem::run() +template +void +CDRProblem::run() { setup_geometry(); setup_dofs(); @@ -308,17 +342,18 @@ void CDRProblem::run() } -constexpr int dim {2}; +constexpr int dim{2}; -int main(int argc, char *argv[]) +int +main(int argc, char *argv[]) { // One of the new features in C++11 is the chrono component of // the standard library. This gives us an easy way to time the output. auto t0 = std::chrono::high_resolution_clock::now(); Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); - CDR::Parameters parameters; + CDR::Parameters parameters; parameters.read_parameter_file("parameters.prm"); CDRProblem cdr_problem(parameters); cdr_problem.run(); @@ -327,9 +362,10 @@ int main(int argc, char *argv[]) if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) { std::cout << "time elapsed: " - << std::chrono::duration_cast(t1 - t0).count() - << " milliseconds." - << std::endl; + << std::chrono::duration_cast(t1 - + t0) + .count() + << " milliseconds." << std::endl; } return 0; -- 2.39.5