--- /dev/null
+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"
+ )
--- /dev/null
+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)
--- /dev/null
+#ifndef dealii__cdr_assemble_system_h
+#define dealii__cdr_assemble_system_h
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function_parser.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <array>
+#include <functional>
+
+#include <deal.II-cdr/parameters.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename MatrixType, typename VectorType>
+ void assemble_system
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const std::function<double(double, const Point<dim>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const VectorType ¤t_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ MatrixType &system_matrix,
+ VectorType &system_rhs);
+}
+#endif
--- /dev/null
+#ifndef dealii__cdr_assemble_system_templates_h
+#define dealii__cdr_assemble_system_templates_h
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <vector>
+
+#include <deal.II-cdr/assemble_system.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename MatrixType, typename VectorType>
+ void assemble_system
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const std::function<double(double, const Point<dim>)> 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<dim> fe_values(fe, quad, update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ Vector<double> cell_rhs(dofs_per_cell);
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+
+ Vector<double> current_fe_coefficients(dofs_per_cell);
+ std::vector<types::global_dof_index> 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
--- /dev/null
+#ifndef dealii__cdr_parameters_h
+#define dealii__cdr_parameters_h
+
+#include <fstream>
+#include <string>
+
+#include <deal.II/base/parameter_handler.h>
+
+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
--- /dev/null
+#ifndef dealii__cdr_system_matrix_h
+#define dealii__cdr_system_matrix_h
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II-cdr/parameters.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename Matrix>
+ void create_system_matrix
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ Matrix &system_matrix);
+
+ template<int dim, typename Matrix>
+ void create_system_matrix
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ const ConstraintMatrix &constraints,
+ Matrix &system_matrix);
+}
+#endif
--- /dev/null
+#ifndef dealii__cdr_system_matrix_templates_h
+#define dealii__cdr_system_matrix_templates_h
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II-cdr/system_matrix.h>
+
+#include <vector>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename UpdateFunction>
+ void internal_create_system_matrix
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> 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<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ FEValues<dim> fe_values(fe, quad, update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ std::vector<types::global_dof_index> 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<int dim, typename Matrix>
+ void create_system_matrix
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ const ConstraintMatrix &constraints,
+ Matrix &system_matrix)
+ {
+ internal_create_system_matrix<dim>
+ (dof_handler, quad, convection_function, parameters, time_step,
+ [&constraints, &system_matrix](const std::vector<types::global_dof_index> &local_indices,
+ const FullMatrix<double> &cell_matrix)
+ {
+ constraints.distribute_local_to_global
+ (cell_matrix, local_indices, system_matrix);
+ });
+ }
+
+ template<int dim, typename Matrix>
+ void create_system_matrix
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ Matrix &system_matrix)
+ {
+ internal_create_system_matrix<dim>
+ (dof_handler, quad, convection_function, parameters, time_step,
+ [&system_matrix](const std::vector<types::global_dof_index> &local_indices,
+ const FullMatrix<double> &cell_matrix)
+ {
+ system_matrix.add(local_indices, cell_matrix);
+ });
+ }
+}
+#endif
--- /dev/null
+#ifndef dealii__cdr_system_rhs_h
+#define dealii__cdr_system_rhs_h
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II-cdr/parameters.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename VectorType>
+ void create_system_rhs
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const std::function<double(double, const Point<dim>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const VectorType &previous_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ VectorType &system_rhs);
+}
+#endif
--- /dev/null
+#ifndef dealii__cdr_system_rhs_templates_h
+#define dealii__cdr_system_rhs_templates_h
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <vector>
+
+#include <deal.II-cdr/system_rhs.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename VectorType>
+ void create_system_rhs
+ (const DoFHandler<dim> &dof_handler,
+ const QGauss<dim> &quad,
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function,
+ const std::function<double(double, const Point<dim>)> 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<dim> fe_values(fe, quad, update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ Vector<double> cell_rhs(dofs_per_cell);
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+
+ Vector<double> current_fe_coefficients(dofs_per_cell);
+ std::vector<types::global_dof_index> 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
--- /dev/null
+#ifndef dealii__cdr_write_pvtu_output_h
+#define dealii__cdr_write_pvtu_output_h
+#include <deal.II/dofs/dof_handler.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ class WritePVTUOutput
+ {
+ public:
+ WritePVTUOutput(const unsigned int patch_level);
+
+ template<int dim, typename VectorType>
+ void write_output(const DoFHandler<dim> &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
--- /dev/null
+#ifndef dealii__cdr_write_pvtu_output_templates_h
+#define dealii__cdr_write_pvtu_output_templates_h
+#include <deal.II/base/data_out_base.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II-cdr/write_pvtu_output.h>
+
+#include <string>
+#include <fstream>
+#include <vector>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template<int dim, typename VectorType>
+ void WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
+ const VectorType &solution,
+ const unsigned int time_step_n,
+ const double current_time)
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "u");
+
+ Vector<float> 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<std::string> 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
--- /dev/null
+#ifndef dealii__cdr_write_xmdf_output_h
+#define dealii__cdr_write_xmdf_output_h
+#include <deal.II/base/data_out_base.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_component_interpretation.h>
+
+#include <string>
+#include <vector>
+
+
+namespace CDR
+{
+ using namespace dealii;
+
+ class WriteXDMFOutput
+ {
+ public:
+ WriteXDMFOutput(const unsigned int patch_level,
+ const bool update_mesh_at_each_step = true);
+
+ template<int dim, typename VectorType>
+ void write_output(const DoFHandler<dim> &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<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation;
+ std::vector<XDMFEntry> xdmf_entries;
+ bool write_mesh;
+ };
+}
+#endif
--- /dev/null
+#ifndef dealii__cdr_write_xmdf_output_templates_h
+#define dealii__cdr_write_xmdf_output_templates_h
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II-cdr/write_xdmf_output.h>
+
+namespace CDR
+{
+ template<int dim, typename VectorType>
+ void WriteXDMFOutput::write_output(const DoFHandler<dim> &dof_handler,
+ const VectorType &solution,
+ const unsigned int time_step_n,
+ const double current_time)
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "u", DataOut<dim>::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
--- /dev/null
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II-cdr/assemble_system.templates.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template
+ void assemble_system<2, SparseMatrix<double>, Vector<double>>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+ const std::function<double(double, const Point<2>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const Vector<double> ¤t_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ SparseMatrix<double> &system_matrix,
+ Vector<double> &system_rhs);
+
+ template
+ void assemble_system<3, SparseMatrix<double>, Vector<double>>
+ (const DoFHandler<3> &dof_handler,
+ const QGauss<3> &quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> convection_function,
+ const std::function<double(double, const Point<3>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const Vector<double> ¤t_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ SparseMatrix<double> &system_matrix,
+ Vector<double> &system_rhs);
+
+ template
+ void assemble_system<2, TrilinosWrappers::SparseMatrix,
+ TrilinosWrappers::MPI::Vector>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+ const std::function<double(double, const Point<2>)> 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<Tensor<1, 3>(const Point<3>)> convection_function,
+ const std::function<double(double, const Point<3>)> 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);
+}
--- /dev/null
+#include <deal.II-cdr/parameters.h>
+
+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();
+ }
+}
--- /dev/null
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+
+#include <deal.II-cdr/parameters.h>
+#include <deal.II-cdr/system_matrix.h>
+#include <deal.II-cdr/system_matrix.templates.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template
+ void create_system_matrix<2, SparseMatrix<double>>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ SparseMatrix<double> &system_matrix);
+
+ template
+ void create_system_matrix<3, SparseMatrix<double>>
+ (const DoFHandler<3> &dof_handler,
+ const QGauss<3> &quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ SparseMatrix<double> &system_matrix);
+
+ template
+ void create_system_matrix<2, SparseMatrix<double>>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ const ConstraintMatrix &constraints,
+ SparseMatrix<double> &system_matrix);
+
+ template
+ void create_system_matrix<3, SparseMatrix<double>>
+ (const DoFHandler<3> &dof_handler,
+ const QGauss<3> &quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ const ConstraintMatrix &constraints,
+ SparseMatrix<double> &system_matrix);
+
+ template
+ void create_system_matrix<2, TrilinosWrappers::SparseMatrix>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(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<Tensor<1, 3>(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<Tensor<1, 2>(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<Tensor<1, 3>(const Point<3>)> convection_function,
+ const CDR::Parameters ¶meters,
+ const double time_step,
+ const ConstraintMatrix &constraints,
+ TrilinosWrappers::SparseMatrix &system_matrix);
+}
--- /dev/null
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II-cdr/system_rhs.templates.h>
+
+namespace CDR
+{
+ using namespace dealii;
+
+ template
+ void create_system_rhs<2, Vector<double>>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+ const std::function<double(double, const Point<2>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const Vector<double> &previous_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ Vector<double> &system_rhs);
+
+ template
+ void create_system_rhs<3, Vector<double>>
+ (const DoFHandler<3> &dof_handler,
+ const QGauss<3> &quad,
+ const std::function<Tensor<1, 3>(const Point<3>)> convection_function,
+ const std::function<double(double, const Point<3>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const Vector<double> &previous_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ Vector<double> &system_rhs);
+
+ template
+ void create_system_rhs<2, TrilinosWrappers::MPI::Vector>
+ (const DoFHandler<2> &dof_handler,
+ const QGauss<2> &quad,
+ const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+ const std::function<double(double, const Point<2>)> 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<Tensor<1, 3>(const Point<3>)> convection_function,
+ const std::function<double(double, const Point<3>)> forcing_function,
+ const CDR::Parameters ¶meters,
+ const TrilinosWrappers::MPI::Vector &previous_solution,
+ const ConstraintMatrix &constraints,
+ const double current_time,
+ TrilinosWrappers::MPI::Vector &system_rhs);
+}
--- /dev/null
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II-cdr/write_pvtu_output.templates.h>
+
+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<double> &solution,
+ const unsigned int time_step_n,
+ const double current_time);
+
+ template
+ void WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
+ const Vector<double> &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);
+}
--- /dev/null
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II-cdr/write_xdmf_output.templates.h>
+
+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<double> &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<double> &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);
+}
--- /dev/null
+David R. Wells wellsd2@rpi.edu
--- /dev/null
+step-19
+step-40
--- /dev/null
+DEAL_II_WITH_CXX11
+DEAL_II_WITH_MPI
+DEAL_II_WITH_TRILINOS
+DEAL_II_WITH_ZLIB
--- /dev/null
+Parallel convection-diffusion-reaction (CDR) solver with theta time stepping utilizing Trilinos.
--- /dev/null
+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
--- /dev/null
+# 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.
--- /dev/null
+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)
--- /dev/null
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+// for distributed computations
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/solution_transfer.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <array>
+#include <functional>
+#include <chrono>
+#include <iostream>
+#include <vector>
+
+#include <deal.II-cdr/system_matrix.h>
+#include <deal.II-cdr/system_rhs.h>
+#include <deal.II-cdr/parameters.h>
+#include <deal.II-cdr/write_pvtu_output.h>
+
+using namespace dealii;
+
+constexpr int manifold_id {0};
+
+
+template<int dim>
+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<dim> fe;
+ QGauss<dim> quad;
+ const SphericalManifold<dim> boundary_description;
+ parallel::distributed::Triangulation<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+
+ const std::function<Tensor<1, dim>(const Point<dim>)> convection_function;
+ const std::function<double(double, const Point<dim>)> 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<int dim>
+CDRProblem<dim>::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<dim>()),
+ triangulation(mpi_communicator, typename Triangulation<dim>::MeshSmoothing
+ (Triangulation<dim>::smoothing_on_refinement |
+ Triangulation<dim>::smoothing_on_coarsening)),
+ dof_handler(triangulation),
+ convection_function
+{
+ [](const Point<dim> p) -> Tensor<1, dim>
+ {Tensor<1, dim> v; v[0] = -p[1]; v[1] = p[0]; return v;}},
+forcing_function
+{
+ [](double t, const Point<dim> 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<int dim>
+void CDRProblem<dim>::setup_geometry()
+{
+ const Point<dim> 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<int dim>
+void CDRProblem<dim>::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<int dim>
+void CDRProblem<dim>::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<dim>
+ (dof_handler, quad, convection_function, parameters, time_step, constraints,
+ system_matrix);
+ system_matrix.compress(VectorOperation::add);
+ preconditioner.initialize(system_matrix);
+}
+
+
+template<int dim>
+void CDRProblem<dim>::time_iterate()
+{
+ double current_time = parameters.start_time;
+ CDR::WritePVTUOutput pvtu_output(parameters.patch_level);
+ for (unsigned int time_step_n = 0; time_step_n < parameters.n_time_steps;
+ ++time_step_n)
+ {
+ current_time += time_step;
+
+ system_rhs = 0.0;
+ CDR::create_system_rhs<dim>
+ (dof_handler, quad, convection_function, forcing_function, parameters,
+ locally_relevant_solution, constraints, current_time, system_rhs);
+ system_rhs.compress(VectorOperation::add);
+
+ SolverControl solver_control(dof_handler.n_dofs(),
+ 1e-6*system_rhs.l2_norm(),
+ /*log_history = */ false,
+ /*log_result = */ false);
+ TrilinosWrappers::SolverGMRES solver(solver_control);
+ solver.solve(system_matrix, completely_distributed_solution, system_rhs,
+ preconditioner);
+ constraints.distribute(completely_distributed_solution);
+ locally_relevant_solution = completely_distributed_solution;
+
+ if (time_step_n % parameters.save_interval == 0)
+ {
+ pvtu_output.write_output(dof_handler, locally_relevant_solution,
+ time_step_n, current_time);
+ }
+
+ refine_mesh();
+ }
+}
+
+
+template<int dim>
+void CDRProblem<dim>::refine_mesh()
+{
+ parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
+ solution_transfer(dof_handler);
+
+ Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+ KellyErrorEstimator<dim>::estimate
+ (dof_handler, QGauss<dim - 1>(fe.degree + 1), typename FunctionMap<dim>::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<int dim>
+void CDRProblem<dim>::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<dim> cdr_problem(parameters);
+ cdr_problem.run();
+
+ auto t1 = std::chrono::high_resolution_clock::now();
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ {
+ std::cout << "time elapsed: "
+ << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count()
+ << " milliseconds."
+ << std::endl;
+ }
+
+ return 0;
+}
projects / code-gallery.git / commitdiff