]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
Added a parallel CDR solver. 2/head
authorDavid Wells <wellsd2@rpi.edu>
Sat, 17 Oct 2015 18:18:35 +0000 (14:18 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Mon, 19 Oct 2015 17:12:56 +0000 (13:12 -0400)
27 files changed:
cdr/CMakeLists.txt [new file with mode: 0644]
cdr/common/CMakeLists.txt [new file with mode: 0644]
cdr/common/include/deal.II-cdr/assemble_system.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/assemble_system.templates.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/parameters.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/system_matrix.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/system_matrix.templates.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/system_rhs.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/system_rhs.templates.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/write_pvtu_output.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/write_pvtu_output.templates.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/write_xdmf_output.h [new file with mode: 0644]
cdr/common/include/deal.II-cdr/write_xdmf_output.templates.h [new file with mode: 0644]
cdr/common/source/assemble_system.cc [new file with mode: 0644]
cdr/common/source/parameters.cc [new file with mode: 0644]
cdr/common/source/system_matrix.cc [new file with mode: 0644]
cdr/common/source/system_rhs.cc [new file with mode: 0644]
cdr/common/source/write_pvtu_output.cc [new file with mode: 0644]
cdr/common/source/write_xdmf_output.cc [new file with mode: 0644]
cdr/doc/author [new file with mode: 0644]
cdr/doc/builds-on [new file with mode: 0644]
cdr/doc/dependencies [new file with mode: 0644]
cdr/doc/tooltip [new file with mode: 0644]
cdr/parameters.prm [new file with mode: 0644]
cdr/readme.md [new file with mode: 0644]
cdr/solver/CMakeLists.txt [new file with mode: 0644]
cdr/solver/cdr.cc [new file with mode: 0644]

diff --git a/cdr/CMakeLists.txt b/cdr/CMakeLists.txt
new file mode 100644 (file)
index 0000000..8e54d00
--- /dev/null
@@ -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 (file)
index 0000000..626efe9
--- /dev/null
@@ -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 (file)
index 0000000..849e673
--- /dev/null
@@ -0,0 +1,34 @@
+#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                                 &parameters,
+   const VectorType                                      &current_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 (file)
index 0000000..f76dc94
--- /dev/null
@@ -0,0 +1,98 @@
+#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                                 &parameters,
+   const VectorType                                      &current_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
diff --git a/cdr/common/include/deal.II-cdr/parameters.h b/cdr/common/include/deal.II-cdr/parameters.h
new file mode 100644 (file)
index 0000000..4155338
--- /dev/null
@@ -0,0 +1,39 @@
+#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 &parameter_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 (file)
index 0000000..afd4f7f
--- /dev/null
@@ -0,0 +1,34 @@
+#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                                 &parameters,
+   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                                 &parameters,
+   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 (file)
index 0000000..23629d0
--- /dev/null
@@ -0,0 +1,106 @@
+#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                                 &parameters,
+   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                                 &parameters,
+   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                                 &parameters,
+   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
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 (file)
index 0000000..ee921a8
--- /dev/null
@@ -0,0 +1,27 @@
+#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                                 &parameters,
+   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 (file)
index 0000000..7b5fd4b
--- /dev/null
@@ -0,0 +1,94 @@
+#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                                 &parameters,
+   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
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 (file)
index 0000000..218d9ed
--- /dev/null
@@ -0,0 +1,24 @@
+#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
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 (file)
index 0000000..b5b709c
--- /dev/null
@@ -0,0 +1,74 @@
+#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
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 (file)
index 0000000..6a1a696
--- /dev/null
@@ -0,0 +1,40 @@
+#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
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 (file)
index 0000000..e944d92
--- /dev/null
@@ -0,0 +1,52 @@
+#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
diff --git a/cdr/common/source/assemble_system.cc b/cdr/common/source/assemble_system.cc
new file mode 100644 (file)
index 0000000..73c5bba
--- /dev/null
@@ -0,0 +1,65 @@
+#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                               &parameters,
+   const Vector<double>                                &current_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                               &parameters,
+   const Vector<double>                                &current_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                               &parameters,
+   const TrilinosWrappers::MPI::Vector                 &current_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                               &parameters,
+   const TrilinosWrappers::MPI::Vector                 &current_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 (file)
index 0000000..f59e786
--- /dev/null
@@ -0,0 +1,110 @@
+#include <deal.II-cdr/parameters.h>
+
+namespace CDR
+{
+  void Parameters::configure_parameter_handler(ParameterHandler &parameter_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 (file)
index 0000000..a1934e0
--- /dev/null
@@ -0,0 +1,87 @@
+#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                             &parameters,
+   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                             &parameters,
+   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                             &parameters,
+   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                             &parameters,
+   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                             &parameters,
+   const double                                      time_step,
+   TrilinosWrappers::SparseMatrix                    &system_matrix);
+
+  template
+  void create_system_matrix<3, TrilinosWrappers::SparseMatrix>
+  (const DoFHandler<3>                               &dof_handler,
+   const QGauss<3>                                   &quad,
+   const std::function<Tensor<1, 3>(const Point<3>)> convection_function,
+   const CDR::Parameters                             &parameters,
+   const double                                      time_step,
+   TrilinosWrappers::SparseMatrix                    &system_matrix);
+
+  template
+  void create_system_matrix<2, TrilinosWrappers::SparseMatrix>
+  (const DoFHandler<2>                               &dof_handler,
+   const QGauss<2>                                   &quad,
+   const std::function<Tensor<1, 2>(const Point<2>)> convection_function,
+   const CDR::Parameters                             &parameters,
+   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                             &parameters,
+   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 (file)
index 0000000..c3ee763
--- /dev/null
@@ -0,0 +1,57 @@
+#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                               &parameters,
+   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                               &parameters,
+   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                               &parameters,
+   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                               &parameters,
+   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 (file)
index 0000000..7bfe3de
--- /dev/null
@@ -0,0 +1,37 @@
+#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);
+}
diff --git a/cdr/common/source/write_xdmf_output.cc b/cdr/common/source/write_xdmf_output.cc
new file mode 100644 (file)
index 0000000..4a6b3a7
--- /dev/null
@@ -0,0 +1,42 @@
+#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);
+}
diff --git a/cdr/doc/author b/cdr/doc/author
new file mode 100644 (file)
index 0000000..751152f
--- /dev/null
@@ -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 (file)
index 0000000..abaea4c
--- /dev/null
@@ -0,0 +1,2 @@
+step-19
+step-40
diff --git a/cdr/doc/dependencies b/cdr/doc/dependencies
new file mode 100644 (file)
index 0000000..b64b9f4
--- /dev/null
@@ -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 (file)
index 0000000..b8fe8a6
--- /dev/null
@@ -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 (file)
index 0000000..a6f9fa2
--- /dev/null
@@ -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 (file)
index 0000000..6eabbed
--- /dev/null
@@ -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 (file)
index 0000000..27eecd5
--- /dev/null
@@ -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 (file)
index 0000000..6a24723
--- /dev/null
@@ -0,0 +1,311 @@
+#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 &parameters);
+  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 &parameters) :
+  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;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.