From acca652fe069475e3ee9d8ceeecd165fdcaf8e93 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 17 Oct 2015 14:18:35 -0400 Subject: [PATCH] Added a parallel CDR solver. --- cdr/CMakeLists.txt | 39 +++ cdr/common/CMakeLists.txt | 7 + .../include/deal.II-cdr/assemble_system.h | 34 ++ .../deal.II-cdr/assemble_system.templates.h | 98 ++++++ cdr/common/include/deal.II-cdr/parameters.h | 39 +++ .../include/deal.II-cdr/system_matrix.h | 34 ++ .../deal.II-cdr/system_matrix.templates.h | 106 ++++++ cdr/common/include/deal.II-cdr/system_rhs.h | 27 ++ .../deal.II-cdr/system_rhs.templates.h | 94 ++++++ .../include/deal.II-cdr/write_pvtu_output.h | 24 ++ .../deal.II-cdr/write_pvtu_output.templates.h | 74 +++++ .../include/deal.II-cdr/write_xdmf_output.h | 40 +++ .../deal.II-cdr/write_xdmf_output.templates.h | 52 +++ cdr/common/source/assemble_system.cc | 65 ++++ cdr/common/source/parameters.cc | 110 +++++++ cdr/common/source/system_matrix.cc | 87 +++++ cdr/common/source/system_rhs.cc | 57 ++++ cdr/common/source/write_pvtu_output.cc | 37 +++ cdr/common/source/write_xdmf_output.cc | 42 +++ cdr/doc/author | 1 + cdr/doc/builds-on | 2 + cdr/doc/dependencies | 4 + cdr/doc/tooltip | 1 + cdr/parameters.prm | 27 ++ cdr/readme.md | 55 ++++ cdr/solver/CMakeLists.txt | 7 + cdr/solver/cdr.cc | 311 ++++++++++++++++++ 27 files changed, 1474 insertions(+) create mode 100644 cdr/CMakeLists.txt create mode 100644 cdr/common/CMakeLists.txt create mode 100644 cdr/common/include/deal.II-cdr/assemble_system.h create mode 100644 cdr/common/include/deal.II-cdr/assemble_system.templates.h create mode 100644 cdr/common/include/deal.II-cdr/parameters.h create mode 100644 cdr/common/include/deal.II-cdr/system_matrix.h create mode 100644 cdr/common/include/deal.II-cdr/system_matrix.templates.h create mode 100644 cdr/common/include/deal.II-cdr/system_rhs.h create mode 100644 cdr/common/include/deal.II-cdr/system_rhs.templates.h create mode 100644 cdr/common/include/deal.II-cdr/write_pvtu_output.h create mode 100644 cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h create mode 100644 cdr/common/include/deal.II-cdr/write_xdmf_output.h create mode 100644 cdr/common/include/deal.II-cdr/write_xdmf_output.templates.h create mode 100644 cdr/common/source/assemble_system.cc create mode 100644 cdr/common/source/parameters.cc create mode 100644 cdr/common/source/system_matrix.cc create mode 100644 cdr/common/source/system_rhs.cc create mode 100644 cdr/common/source/write_pvtu_output.cc create mode 100644 cdr/common/source/write_xdmf_output.cc create mode 100644 cdr/doc/author create mode 100644 cdr/doc/builds-on create mode 100644 cdr/doc/dependencies create mode 100644 cdr/doc/tooltip create mode 100644 cdr/parameters.prm create mode 100644 cdr/readme.md create mode 100644 cdr/solver/CMakeLists.txt create mode 100644 cdr/solver/cdr.cc diff --git a/cdr/CMakeLists.txt b/cdr/CMakeLists.txt new file mode 100644 index 0000000..8e54d00 --- /dev/null +++ b/cdr/CMakeLists.txt @@ -0,0 +1,39 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.4 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}) + +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +PROJECT(deal.II-cdr) +# Put the executable in the root directory (the same directory as this file) for +# easy access +SET(TARGET "cdr") +SET(EXECUTABLE_OUTPUT_PATH ${CMAKE_SOURCE_DIR}) + +ADD_SUBDIRECTORY("common") +ADD_SUBDIRECTORY("solver") + +ADD_CUSTOM_TARGET(debug + COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all + COMMENT "Switch CMAKE_BUILD_TYPE to Debug" + ) + +ADD_CUSTOM_TARGET(release + COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all + COMMENT "Switch CMAKE_BUILD_TYPE to Release" + ) + +ADD_CUSTOM_TARGET(run COMMAND ${CMAKE_SOURCE_DIR}/${TARGET} + COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration" + ) diff --git a/cdr/common/CMakeLists.txt b/cdr/common/CMakeLists.txt new file mode 100644 index 0000000..626efe9 --- /dev/null +++ b/cdr/common/CMakeLists.txt @@ -0,0 +1,7 @@ +INCLUDE_DIRECTORIES("include") + +FILE(GLOB_RECURSE DEAL_II_CDR_COMMON_SOURCES *.cc) + +ADD_LIBRARY(deal.II-cdr-common SHARED ${DEAL_II_CDR_COMMON_SOURCES}) + +DEAL_II_SETUP_TARGET(deal.II-cdr-common) diff --git a/cdr/common/include/deal.II-cdr/assemble_system.h b/cdr/common/include/deal.II-cdr/assemble_system.h new file mode 100644 index 0000000..849e673 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/assemble_system.h @@ -0,0 +1,34 @@ +#ifndef dealii__cdr_assemble_system_h +#define dealii__cdr_assemble_system_h +#include +#include + +#include + +#include + +#include + +#include +#include + +#include + +namespace CDR +{ + using namespace dealii; + + template + void assemble_system + (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 ¤t_solution, + const ConstraintMatrix &constraints, + const double current_time, + MatrixType &system_matrix, + VectorType &system_rhs); +} +#endif diff --git a/cdr/common/include/deal.II-cdr/assemble_system.templates.h b/cdr/common/include/deal.II-cdr/assemble_system.templates.h new file mode 100644 index 0000000..f76dc94 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/assemble_system.templates.h @@ -0,0 +1,98 @@ +#ifndef dealii__cdr_assemble_system_templates_h +#define dealii__cdr_assemble_system_templates_h +#include + +#include + +#include + +#include + +namespace CDR +{ + using namespace dealii; + + template + void assemble_system + (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 ¤t_solution, + const ConstraintMatrix &constraints, + const double current_time, + MatrixType &system_matrix, + 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); + + Vector cell_rhs(dofs_per_cell); + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + + Vector current_fe_coefficients(dofs_per_cell); + std::vector local_indices(dofs_per_cell); + + auto previous_time = current_time - time_step; + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_rhs = 0.0; + cell_matrix = 0.0; + cell->get_dof_indices(local_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + current_fe_coefficients[i] = current_solution[local_indices[i]]; + } + + for (unsigned int q = 0; q < quad.size(); ++q) + { + const auto current_convection = + convection_function(fe_values.quadrature_point(q)); + const auto current_forcing = + forcing_function(current_time, fe_values.quadrature_point(q)); + const auto 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 double mass_part + {fe_values.shape_value(i, q)*fe_values.shape_value(j, q)}; + const double system_part + {time_step/2.0* + // convection part + (fe_values.shape_value(i, q)*convection_contribution + // diffusion part + + parameters.diffusion_coefficient + *fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q) + // reaction part + + parameters.reaction_coefficient*mass_part)}; + cell_rhs(i) += fe_values.JxW(q) + *((mass_part - system_part)*current_fe_coefficients[j] + + time_step/2.0*fe_values.shape_value(i, q) + *(current_forcing + previous_forcing)); + cell_matrix(i, j) += + fe_values.JxW(q)*(mass_part + system_part); + } + } + } + constraints.distribute_local_to_global + (cell_matrix, cell_rhs, local_indices, system_matrix, system_rhs, true); + } + } + } +} +#endif diff --git a/cdr/common/include/deal.II-cdr/parameters.h b/cdr/common/include/deal.II-cdr/parameters.h new file mode 100644 index 0000000..4155338 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/parameters.h @@ -0,0 +1,39 @@ +#ifndef dealii__cdr_parameters_h +#define dealii__cdr_parameters_h + +#include +#include + +#include + +namespace CDR +{ + using namespace dealii; + + class Parameters + { + public: + double inner_radius; + double outer_radius; + + double diffusion_coefficient; + double reaction_coefficient; + bool time_dependent_forcing; + + unsigned int initial_refinement_level; + unsigned int max_refinement_level; + unsigned int fe_order; + + double start_time; + double stop_time; + unsigned int n_time_steps; + + unsigned int save_interval; + unsigned int patch_level; + + void read_parameter_file(std::string file_name); + private: + void configure_parameter_handler(ParameterHandler ¶meter_handler); + }; +} +#endif diff --git a/cdr/common/include/deal.II-cdr/system_matrix.h b/cdr/common/include/deal.II-cdr/system_matrix.h new file mode 100644 index 0000000..afd4f7f --- /dev/null +++ b/cdr/common/include/deal.II-cdr/system_matrix.h @@ -0,0 +1,34 @@ +#ifndef dealii__cdr_system_matrix_h +#define dealii__cdr_system_matrix_h +#include + +#include + +#include + +#include + +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, + Matrix &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 ConstraintMatrix &constraints, + Matrix &system_matrix); +} +#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 new file mode 100644 index 0000000..23629d0 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/system_matrix.templates.h @@ -0,0 +1,106 @@ +#ifndef dealii__cdr_system_matrix_templates_h +#define dealii__cdr_system_matrix_templates_h +#include +#include + +#include + +#include + +namespace CDR +{ + using namespace dealii; + + 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) + { + 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); + + std::vector local_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_matrix = 0.0; + cell->get_dof_indices(local_indices); + for (unsigned int q = 0; q < quad.size(); ++q) + { + const auto current_convection = + convection_function(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); + cell_matrix(i, j) += fe_values.JxW(q)* + // 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* + // convection part + (fe_values.shape_value(i, q)*convection_contribution + // Laplacian part + + parameters.diffusion_coefficient + *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q))) + ); + } + } + } + update_system_matrix(local_indices, cell_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 ConstraintMatrix &constraints, + Matrix &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, + Matrix &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); + }); + } +} +#endif diff --git a/cdr/common/include/deal.II-cdr/system_rhs.h b/cdr/common/include/deal.II-cdr/system_rhs.h new file mode 100644 index 0000000..ee921a8 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/system_rhs.h @@ -0,0 +1,27 @@ +#ifndef dealii__cdr_system_rhs_h +#define dealii__cdr_system_rhs_h +#include + +#include + +#include + +#include + +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 ConstraintMatrix &constraints, + const double current_time, + VectorType &system_rhs); +} +#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 new file mode 100644 index 0000000..7b5fd4b --- /dev/null +++ b/cdr/common/include/deal.II-cdr/system_rhs.templates.h @@ -0,0 +1,94 @@ +#ifndef dealii__cdr_system_rhs_templates_h +#define dealii__cdr_system_rhs_templates_h +#include +#include + +#include + +#include + +#include + +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 ConstraintMatrix &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); + + Vector cell_rhs(dofs_per_cell); + FullMatrix cell_matrix(dofs_per_cell, 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}; + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + cell_rhs = 0.0; + 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]]; + } + + for (unsigned int q = 0; q < quad.size(); ++q) + { + 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)); + 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); + + cell_rhs(i) += fe_values.JxW(q)* + // 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* + // convection part + (fe_values.shape_value(i, q)*convection_contribution + // Laplacian part + + parameters.diffusion_coefficient + *(fe_values.shape_grad(i, q)*fe_values.shape_grad(j, q)))) + *current_fe_coefficients[j] + // forcing parts + + 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); + } + } + } +} +#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 new file mode 100644 index 0000000..218d9ed --- /dev/null +++ b/cdr/common/include/deal.II-cdr/write_pvtu_output.h @@ -0,0 +1,24 @@ +#ifndef dealii__cdr_write_pvtu_output_h +#define dealii__cdr_write_pvtu_output_h +#include + +namespace CDR +{ + using namespace dealii; + + class WritePVTUOutput + { + 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); + private: + const unsigned int patch_level; + const unsigned int this_mpi_process; + }; +} +#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 new file mode 100644 index 0000000..b5b709c --- /dev/null +++ b/cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h @@ -0,0 +1,74 @@ +#ifndef dealii__cdr_write_pvtu_output_templates_h +#define dealii__cdr_write_pvtu_output_templates_h +#include +#include + +#include + +#include + +#include + +#include +#include +#include + +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) + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u"); + + Vector subdomain (dof_handler.get_tria().n_active_cells()); + for (auto &domain : subdomain) + { + domain = dof_handler.get_tria().locally_owned_subdomain(); + } + data_out.add_data_vector(subdomain, "subdomain"); + data_out.build_patches(patch_level); + + DataOutBase::VtkFlags flags; + flags.time = current_time; + flags.compression_level = DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; + data_out.set_flags(flags); + + unsigned int subdomain_n; + if (Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1) + { + subdomain_n = 0; + } + else + { + subdomain_n = dof_handler.get_tria().locally_owned_subdomain(); + } + + 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); + ++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"); + data_out.write_pvtu_record(master_output, filenames); + } + } +} +#endif diff --git a/cdr/common/include/deal.II-cdr/write_xdmf_output.h b/cdr/common/include/deal.II-cdr/write_xdmf_output.h new file mode 100644 index 0000000..6a1a696 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/write_xdmf_output.h @@ -0,0 +1,40 @@ +#ifndef dealii__cdr_write_xmdf_output_h +#define dealii__cdr_write_xmdf_output_h +#include + +#include + +#include + +#include + +#include +#include + + +namespace CDR +{ + using namespace dealii; + + class WriteXDMFOutput + { + public: + WriteXDMFOutput(const unsigned int patch_level, + const bool update_mesh_at_each_step = true); + + 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 bool update_mesh_at_each_step; + const std::string xdmf_file_name; + const std::vector + data_component_interpretation; + std::vector xdmf_entries; + bool write_mesh; + }; +} +#endif diff --git a/cdr/common/include/deal.II-cdr/write_xdmf_output.templates.h b/cdr/common/include/deal.II-cdr/write_xdmf_output.templates.h new file mode 100644 index 0000000..e944d92 --- /dev/null +++ b/cdr/common/include/deal.II-cdr/write_xdmf_output.templates.h @@ -0,0 +1,52 @@ +#ifndef dealii__cdr_write_xmdf_output_templates_h +#define dealii__cdr_write_xmdf_output_templates_h +#include + +#include + +#include + +namespace CDR +{ + template + void WriteXDMFOutput::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", DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches(patch_level); + + auto solution_file_name = "solution-" + + Utilities::int_to_string(time_step_n, 9) + ".h5"; + std::string mesh_file_name; + if (update_mesh_at_each_step) + { + mesh_file_name = "mesh-" + + Utilities::int_to_string(time_step_n, 9) + ".h5"; + } + else + { + mesh_file_name = "mesh.h5"; + } + + DataOutBase::DataOutFilter data_filter + (DataOutBase::DataOutFilterFlags(true, true)); + data_out.write_filtered_data(data_filter); + data_out.write_hdf5_parallel(data_filter, write_mesh, mesh_file_name, + solution_file_name, MPI_COMM_WORLD); + if (!update_mesh_at_each_step) + { + write_mesh = false; + } + + xdmf_entries.push_back(data_out.create_xdmf_entry + (data_filter, mesh_file_name, solution_file_name, + current_time, MPI_COMM_WORLD)); + data_out.write_xdmf_file(xdmf_entries, xdmf_file_name, MPI_COMM_WORLD); + } +} +#endif diff --git a/cdr/common/source/assemble_system.cc b/cdr/common/source/assemble_system.cc new file mode 100644 index 0000000..73c5bba --- /dev/null +++ b/cdr/common/source/assemble_system.cc @@ -0,0 +1,65 @@ +#include +#include + +#include + +#include + +namespace CDR +{ + using namespace dealii; + + template + void assemble_system<2, SparseMatrix, 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 ¤t_solution, + const ConstraintMatrix &constraints, + const double current_time, + SparseMatrix &system_matrix, + Vector &system_rhs); + + template + void assemble_system<3, SparseMatrix, 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 ¤t_solution, + const ConstraintMatrix &constraints, + const double current_time, + SparseMatrix &system_matrix, + Vector &system_rhs); + + template + void assemble_system<2, TrilinosWrappers::SparseMatrix, + 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 ¤t_solution, + const ConstraintMatrix &constraints, + const double current_time, + TrilinosWrappers::SparseMatrix &system_matrix, + TrilinosWrappers::MPI::Vector &system_rhs); + + template + void assemble_system<3, TrilinosWrappers::SparseMatrix, + 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 ¤t_solution, + const ConstraintMatrix &constraints, + const double current_time, + TrilinosWrappers::SparseMatrix &system_matrix, + TrilinosWrappers::MPI::Vector &system_rhs); +} diff --git a/cdr/common/source/parameters.cc b/cdr/common/source/parameters.cc new file mode 100644 index 0000000..f59e786 --- /dev/null +++ b/cdr/common/source/parameters.cc @@ -0,0 +1,110 @@ +#include + +namespace CDR +{ + 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.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.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.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.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.leave_subsection(); + } + + void Parameters::read_parameter_file(std::string file_name) + { + ParameterHandler parameter_handler; + { + std::ifstream file(file_name); + configure_parameter_handler(parameter_handler); + parameter_handler.read_input(file); + } + + parameter_handler.enter_subsection("Geometry"); + { + inner_radius = parameter_handler.get_double("inner_radius"); + outer_radius = parameter_handler.get_double("outer_radius"); + } + parameter_handler.leave_subsection(); + + 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"); + } + 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"); + 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"); + n_time_steps = parameter_handler.get_integer("n_time_steps"); + } + parameter_handler.leave_subsection(); + + parameter_handler.enter_subsection("Output"); + { + save_interval = parameter_handler.get_integer("save_interval"); + patch_level = parameter_handler.get_integer("patch_level"); + } + parameter_handler.leave_subsection(); + } +} diff --git a/cdr/common/source/system_matrix.cc b/cdr/common/source/system_matrix.cc new file mode 100644 index 0000000..a1934e0 --- /dev/null +++ b/cdr/common/source/system_matrix.cc @@ -0,0 +1,87 @@ +#include +#include + +#include +#include +#include + +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<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<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 ConstraintMatrix &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 ConstraintMatrix &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<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<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 ConstraintMatrix &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 ConstraintMatrix &constraints, + TrilinosWrappers::SparseMatrix &system_matrix); +} diff --git a/cdr/common/source/system_rhs.cc b/cdr/common/source/system_rhs.cc new file mode 100644 index 0000000..c3ee763 --- /dev/null +++ b/cdr/common/source/system_rhs.cc @@ -0,0 +1,57 @@ +#include +#include + +#include + +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 ConstraintMatrix &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 ConstraintMatrix &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 ConstraintMatrix &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 ConstraintMatrix &constraints, + const double current_time, + TrilinosWrappers::MPI::Vector &system_rhs); +} diff --git a/cdr/common/source/write_pvtu_output.cc b/cdr/common/source/write_pvtu_output.cc new file mode 100644 index 0000000..7bfe3de --- /dev/null +++ b/cdr/common/source/write_pvtu_output.cc @@ -0,0 +1,37 @@ +#include + +#include + +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)} + {} + + 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<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); +} diff --git a/cdr/common/source/write_xdmf_output.cc b/cdr/common/source/write_xdmf_output.cc new file mode 100644 index 0000000..4a6b3a7 --- /dev/null +++ b/cdr/common/source/write_xdmf_output.cc @@ -0,0 +1,42 @@ +#include + +#include + +namespace CDR +{ + using namespace dealii; + + WriteXDMFOutput::WriteXDMFOutput(const unsigned int patch_level, + const bool update_mesh_at_each_step) + : patch_level {patch_level}, + update_mesh_at_each_step {update_mesh_at_each_step}, + xdmf_file_name {"solution.xdmf"}, + data_component_interpretation + {DataComponentInterpretation::component_is_scalar}, + write_mesh {true} + {} + + template + void WriteXDMFOutput::write_output(const DoFHandler<2> &dof_handler, + const Vector &solution, + const unsigned int time_step_n, + const double current_time); + + template + void WriteXDMFOutput::write_output(const DoFHandler<2> &dof_handler, + const TrilinosWrappers::MPI::Vector &solution, + const unsigned int time_step_n, + const double current_time); + + template + void WriteXDMFOutput::write_output(const DoFHandler<3> &dof_handler, + const Vector &solution, + const unsigned int time_step_n, + const double current_time); + + template + void WriteXDMFOutput::write_output(const DoFHandler<3> &dof_handler, + const TrilinosWrappers::MPI::Vector &solution, + const unsigned int time_step_n, + const double current_time); +} diff --git a/cdr/doc/author b/cdr/doc/author new file mode 100644 index 0000000..751152f --- /dev/null +++ b/cdr/doc/author @@ -0,0 +1 @@ +David R. Wells wellsd2@rpi.edu diff --git a/cdr/doc/builds-on b/cdr/doc/builds-on new file mode 100644 index 0000000..abaea4c --- /dev/null +++ b/cdr/doc/builds-on @@ -0,0 +1,2 @@ +step-19 +step-40 diff --git a/cdr/doc/dependencies b/cdr/doc/dependencies new file mode 100644 index 0000000..b64b9f4 --- /dev/null +++ b/cdr/doc/dependencies @@ -0,0 +1,4 @@ +DEAL_II_WITH_CXX11 +DEAL_II_WITH_MPI +DEAL_II_WITH_TRILINOS +DEAL_II_WITH_ZLIB diff --git a/cdr/doc/tooltip b/cdr/doc/tooltip new file mode 100644 index 0000000..b8fe8a6 --- /dev/null +++ b/cdr/doc/tooltip @@ -0,0 +1 @@ +Parallel convection-diffusion-reaction (CDR) solver with theta time stepping utilizing Trilinos. diff --git a/cdr/parameters.prm b/cdr/parameters.prm new file mode 100644 index 0000000..a6f9fa2 --- /dev/null +++ b/cdr/parameters.prm @@ -0,0 +1,27 @@ +subsection Geometry + set inner_radius = 1.0 + set outer_radius = 2.0 +end + +subsection Physical Parameters + set diffusion_coefficient = 1.0e-3 + set reaction_coefficient = 1.0e-4 + set time_dependent_forcing = true +end + +subsection Finite Element + set initial_refinement_level = 3 + set max_refinement_level = 4 + set fe_order = 2 +end + +subsection Time Step + set start_time = 0.0 + set stop_time = 2.0 + set n_time_steps = 200 +end + +subsection Output + set save_interval = 1 + set patch_level = 1 +end diff --git a/cdr/readme.md b/cdr/readme.md new file mode 100644 index 0000000..6eabbed --- /dev/null +++ b/cdr/readme.md @@ -0,0 +1,55 @@ +# Gallery Entry: deal.II-cdr + + +## Overview +I started this project with the intent of better understanding adaptive mesh +refinement, parallel computing, and `CMake`. In particular, I started by writing +a uniform mesh, single process solver and then ultimately expanded it into +`solver/cdr.cc`. This example program might be useful to look at if you want to +see: +* A more complex `CMake` setup, which builds a shared object library and an + executable +* A simple parallel time stepping problem +* Use of `C++11` lambda functions + +The other solvers are available [here](http://www.github.com/drwells/dealii-cdr). + + +## Requirements +* A `C++-11` compliant compiler +* The development version of `deal.II` + + +## Compiling and running +Like the example programs, run +``` +cmake -DDEAL_II_DIR=/path/to/deal.II . +make +``` +in this directory. The solver may be run as +``` +make run +``` +or, for parallelism across `16` processes, +``` +mpirun -np 16 ./cdr +``` + + +## Why use convection-diffusion-reaction? +This equation exhibits very fine boundary layers (usually, from the literature, +these layers have width equal to the square root of the diffusion coefficient on +internal layers and the diffusion coefficient on boundary layers). A good way to +solve it is to use adaptive mesh refinement to refine the mesh only at interior +boundary layers. At the same time, this problem is linear (and does not have a +pressure term) so it is much simpler to solve than the Navier-Stokes equations +with comparable diffusion (Reynolds number). + +I use relatively large diffusion values so I can get away without any additional +scheme like SUPG. + + +## Recommended Literature +There are many good books and papers on numerical methods for this equation. A +good starting point is "Robust Numerical Methods for Singularly Perturbed +Problems" by Roos, Stynes, and Tobiska. diff --git a/cdr/solver/CMakeLists.txt b/cdr/solver/CMakeLists.txt new file mode 100644 index 0000000..27eecd5 --- /dev/null +++ b/cdr/solver/CMakeLists.txt @@ -0,0 +1,7 @@ +SET(TARGET_SRC ${TARGET}.cc) + +INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/common/include) + +ADD_EXECUTABLE(${TARGET} ${TARGET}.cc) +DEAL_II_SETUP_TARGET(${TARGET}) +TARGET_LINK_LIBRARIES(${TARGET} deal.II-cdr-common) diff --git a/cdr/solver/cdr.cc b/cdr/solver/cdr.cc new file mode 100644 index 0000000..6a24723 --- /dev/null +++ b/cdr/solver/cdr.cc @@ -0,0 +1,311 @@ +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +// for distributed computations +#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}; + + +template +class CDRProblem +{ +public: + CDRProblem(const CDR::Parameters ¶meters); + void run(); +private: + const CDR::Parameters parameters; + const double time_step; + double current_time; + + 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; + parallel::distributed::Triangulation triangulation; + DoFHandler dof_handler; + + const std::function(const Point)> convection_function; + const std::function)> forcing_function; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + ConstraintMatrix constraints; + bool first_run; + 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(); +}; + + +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(3*(2 + 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, + (Utilities::MPI::this_mpi_process(mpi_communicator) + == 0)) +{ + Assert(dim == 2, ExcNotImplemented()); +} + + +template +void CDRProblem::setup_geometry() +{ + const Point center; + 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()) + { + cell->set_all_manifold_ids(0); + } + triangulation.refine_global(parameters.initial_refinement_level); +} + + +template +void CDRProblem::setup_dofs() +{ + dof_handler.distribute_dofs(fe); + 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); + + 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); + constraints.close(); + + completely_distributed_solution.reinit + (locally_owned_dofs, mpi_communicator); + + locally_relevant_solution.reinit(locally_owned_dofs, locally_relevant_dofs, + mpi_communicator); +} + + +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.n_locally_owned_dofs_per_processor(), + 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 + (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() +{ + 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 + (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 +void CDRProblem::refine_mesh() +{ + parallel::distributed::SolutionTransfer + solution_transfer(dof_handler); + + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate + (dof_handler, QGauss(fe.degree + 1), typename FunctionMap::type(), + locally_relevant_solution, estimated_error_per_cell); + + // Poor man's version of refine and coarsen + 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(); + } + } + + triangulation.prepare_coarsening_and_refinement(); + solution_transfer.prepare_for_coarsening_and_refinement + (locally_relevant_solution); + triangulation.execute_coarsening_and_refinement(); + + setup_dofs(); + + TrilinosWrappers::MPI::Vector temporary + (locally_owned_dofs, mpi_communicator); + solution_transfer.interpolate(temporary); + locally_relevant_solution = temporary; + constraints.distribute(locally_relevant_solution); + setup_system(); +} + + +template +void CDRProblem::run() +{ + setup_geometry(); + setup_dofs(); + setup_system(); + time_iterate(); +} + + +constexpr int dim {2}; + + +int main(int argc, char *argv[]) +{ + 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 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(t1 - t0).count() + << " milliseconds." + << std::endl; + } + + return 0; +} -- 2.39.5