Fix 64bit GridRefinement overflow

GridRefinement::refine_and_coarsen_fixed_fraction causes overflow with
more than 2^32 global active cells. This is now fixed.

diff --git a/include/deal.II/distributed/grid_refinement.h b/include/deal.II/distributed/grid_refinement.h
index 2d6b6fd25da524c3f6bb4a773679bd11dd988ac3..4f96a88bc26626c2327aa369e723d5e4032267f1 100644 (file)
--- a/include/deal.II/distributed/grid_refinement.h
+++ b/include/deal.II/distributed/grid_refinement.h
@@ -57,7 +57,7 @@ namespace internal
           number
           compute_threshold(const dealii::Vector<number> &   criteria,
                             const std::pair<double, double> &global_min_and_max,
-                            const unsigned int               n_target_cells,
+                            const types::global_dof_index    n_target_cells,
                             MPI_Comm                         mpi_communicator);
         } // namespace RefineAndCoarsenFixedNumber
 

diff --git a/source/distributed/grid_refinement.cc b/source/distributed/grid_refinement.cc
index 661780ecafdfdc005de5c0287124de414f47b007..a1625db018f19dc755d2990c1f06b116f38be887 100644 (file)
--- a/source/distributed/grid_refinement.cc
+++ b/source/distributed/grid_refinement.cc
@@ -228,7 +228,7 @@ namespace internal
           number
           compute_threshold(const dealii::Vector<number> &   criteria,
                             const std::pair<double, double> &global_min_and_max,
-                            const unsigned int               n_target_cells,
+                            const types::global_dof_index    n_target_cells,
                             MPI_Comm                         mpi_communicator)
           {
             double interesting_range[2] = {global_min_and_max.first,
@@ -255,21 +255,22 @@ namespace internal
                      std::sqrt(interesting_range[0] * interesting_range[1]) :
                      (interesting_range[0] + interesting_range[1]) / 2);
 
-                // count how many of our own elements would be above this
-                // threshold and then add to it the number for all the others
-                unsigned int my_count =
+                // Count how many of our own elements would be above this
+                // threshold:
+                types::global_dof_index my_count =
                   std::count_if(criteria.begin(),
                                 criteria.end(),
                                 [test_threshold](const double c) {
                                   return c > test_threshold;
                                 });
 
-                unsigned int total_count = 0;
+                // Potentially accumulate in a 64bit int to avoid overflow:
+                types::global_dof_index total_count = 0;
 
                 ierr = MPI_Reduce(&my_count,
                                   &total_count,
                                   1,
-                                  MPI_UNSIGNED,
+                                  DEAL_II_DOF_INDEX_MPI_TYPE,
                                   MPI_SUM,
                                   master_mpi_rank,
                                   mpi_communicator);
@@ -477,8 +478,8 @@ namespace parallel
           GridRefinement::RefineAndCoarsenFixedNumber::compute_threshold(
             locally_owned_indicators,
             global_min_and_max,
-            static_cast<unsigned int>(adjusted_fractions.first *
-                                      tria.n_global_active_cells()),
+            static_cast<types::global_dof_index>(adjusted_fractions.first *
+                                                 tria.n_global_active_cells()),
             mpi_communicator);
 
         // compute bottom threshold only if necessary. otherwise use a threshold
@@ -488,7 +489,7 @@ namespace parallel
             GridRefinement::RefineAndCoarsenFixedNumber::compute_threshold(
               locally_owned_indicators,
               global_min_and_max,
-              static_cast<unsigned int>(
+              static_cast<types::global_dof_index>(
                 std::ceil((1 - adjusted_fractions.second) *
                           tria.n_global_active_cells())),
               mpi_communicator);
@@ -555,7 +556,7 @@ namespace parallel
             GridRefinement::RefineAndCoarsenFixedFraction::compute_threshold(
               locally_owned_indicators,
               global_min_and_max,
-              (1 - bottom_fraction_of_error) * total_error,
+              (1. - bottom_fraction_of_error) * total_error,
               mpi_communicator);
         else
           {

diff --git a/source/distributed/grid_refinement.inst.in b/source/distributed/grid_refinement.inst.in
index 3588a6f19f46db383ba5782753463354938a2f82..3e4477074f2e0d7f06b615a6cd8777bbb87a5d58 100644 (file)
--- a/source/distributed/grid_refinement.inst.in
+++ b/source/distributed/grid_refinement.inst.in
@@ -34,7 +34,7 @@ for (S : REAL_SCALARS)
               template S
               compute_threshold<S>(const dealii::Vector<S> &,
                                    const std::pair<double, double> &,
-                                   const unsigned int,
+                                   const types::global_dof_index,
                                    MPI_Comm);
             \}
             namespace RefineAndCoarsenFixedFraction