From 0bd559bfb701a6f83aded1b158085c8f0fe60c8c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 23 Jul 2015 11:00:51 -0500 Subject: [PATCH] Make sure that compute refinement thresholds are no larger than the max indicator. This could happen because we search the threshold in an interval larger than the actual indicators. If we end up with a threshold that is larger than the largest indicator, we need to cut things back. In the process of debugging, I also found a few oddities that have to do with the fact that we do a lot of arithmetic in double precision, but often pass in criteria as vectors of float. Fix some of those as well. --- source/distributed/grid_refinement.cc | 39 +++- .../refine_and_coarsen_fixed_fraction_07.cc | 221 ++++++++++++++++++ ..._coarsen_fixed_fraction_07.mpirun=1.output | 3 + 3 files changed, 255 insertions(+), 8 deletions(-) create mode 100644 tests/mpi/refine_and_coarsen_fixed_fraction_07.cc create mode 100644 tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output diff --git a/source/distributed/grid_refinement.cc b/source/distributed/grid_refinement.cc index 714f255962..f4ab173345 100644 --- a/source/distributed/grid_refinement.cc +++ b/source/distributed/grid_refinement.cc @@ -34,6 +34,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -76,7 +77,7 @@ namespace * others get a pair of zeros. */ template - std::pair + std::pair compute_global_min_and_max_at_root (const Vector &criteria, MPI_Comm mpi_communicator) { @@ -153,7 +154,7 @@ namespace void get_locally_owned_indicators (const parallel::distributed::Triangulation &tria, const Vector &criteria, - dealii::Vector &locally_owned_indicators) + dealii::Vector &locally_owned_indicators) { Assert (locally_owned_indicators.size() == tria.n_locally_owned_active_cells(), ExcInternalError()); @@ -412,7 +413,22 @@ namespace 0, mpi_communicator); if (interesting_range[0] == interesting_range[1]) - return interesting_range[0]; + { + // so we have found our threshold. since we adjust + // the range at the top of the function to be slightly + // larger than the actual extremes of the refinement + // criteria values, we can end up in a situation where + // the threshold is in fact larger than the maximal + // refinement indicator. in such cases, we get no + // refinement at all. thus, cap the threshold by the + // actual largest value + double final_threshold = std::min (interesting_range[0], + global_min_and_max.second); + MPI_Bcast (&final_threshold, 1, MPI_DOUBLE, + 0, mpi_communicator); + + return final_threshold; + } const double test_threshold = (interesting_range[0] > 0 @@ -500,7 +516,15 @@ namespace 0, mpi_communicator); if (interesting_range[0] == interesting_range[1]) - return interesting_range[0]; + { + // we got the final value. the master adjusts the + // value one more time, so receive it from there + double final_threshold = -1; + MPI_Bcast (&final_threshold, 1, MPI_DOUBLE, + 0, mpi_communicator); + + return final_threshold; + } // count how many elements // there are that are bigger @@ -569,7 +593,7 @@ namespace parallel // vector of indicators the // ones that correspond to // cells that we locally own - dealii::Vector + dealii::Vector locally_owned_indicators (tria.n_locally_owned_active_cells()); get_locally_owned_indicators (tria, criteria, @@ -583,7 +607,7 @@ namespace parallel // need it here, but it's a // collective communication // call - const std::pair global_min_and_max + const std::pair global_min_and_max = compute_global_min_and_max_at_root (locally_owned_indicators, mpi_communicator); @@ -681,7 +705,7 @@ namespace parallel // vector of indicators the // ones that correspond to // cells that we locally own - dealii::Vector + dealii::Vector locally_owned_indicators (tria.n_locally_owned_active_cells()); get_locally_owned_indicators (tria, criteria, @@ -718,7 +742,6 @@ namespace parallel top_fraction_of_error * total_error, mpi_communicator); - // compute bottom // threshold only if // necessary. otherwise diff --git a/tests/mpi/refine_and_coarsen_fixed_fraction_07.cc b/tests/mpi/refine_and_coarsen_fixed_fraction_07.cc new file mode 100644 index 0000000000..5612c11f7e --- /dev/null +++ b/tests/mpi/refine_and_coarsen_fixed_fraction_07.cc @@ -0,0 +1,221 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// A testcase by Andrea Bonito (in modified form): Given a particular +// set of refinement indicators, we get no refinement at all because +// the threshold is computed to be *larger* than the largest +// refinement indicator + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include + + +void test() +{ + const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + // create a mesh with 120 cells (because that's what Andrea's + // original testcase had) + parallel::distributed::Triangulation<2> triangulation(MPI_COMM_WORLD); + std::vector subdivisions(2); + subdivisions[0] = 120; + subdivisions[1] = 1; + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + Point<2>(), + Point<2>(1,1)); + // initialize the refinement indicators with a set of particular values from the original + // testcase + const double values[] = + { + 1.48589e-08, + 3.31859e-06, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 0.0025918, + 0.0025918, + 3.31859e-06, + 1.48589e-08, + 0.0025918, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 0.0025918, + 1.48589e-08, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 0.0025918, + 0.0025918, + 3.31859e-06, + 3.31859e-06, + 1.48589e-08, + 0.0025918, + 0.0025918, + 3.31859e-06, + 3.31859e-06, + 0.0025918, + 0.0025918, + 1.48589e-08, + 3.31859e-06, + 3.31859e-06, + 0.0025918, + 0.0025918, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 1.48589e-08, + 0.0025918, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 0.0025918, + 1.48589e-08, + 3.31859e-06, + 0.0025918, + 0.0025918, + 3.31859e-06, + 0.0025918, + 3.31859e-06, + 3.31859e-06, + 1.48589e-08, + 0.0120662, + 0.0446999, + 0.0446999, + 0.0644361, + 0.0446999, + 0.0644361, + 0.0644361, + 3.14051, + 0.0446999, + 0.0120662, + 0.0644361, + 0.0446999, + 0.0644361, + 0.0446999, + 3.14051, + 0.0644361, + 0.0446999, + 0.0644361, + 0.0120662, + 0.0446999, + 0.0644361, + 3.14051, + 0.0446999, + 0.0644361, + 0.0644361, + 0.0446999, + 0.0446999, + 0.0120662, + 3.14051, + 0.0644361, + 0.0644361, + 0.0446999, + 0.0446999, + 0.0644361, + 0.0644361, + 3.14051, + 0.0120662, + 0.0446999, + 0.0446999, + 0.0644361, + 0.0644361, + 0.0446999, + 3.14051, + 0.0644361, + 0.0446999, + 0.0120662, + 0.0644361, + 0.0446999, + 0.0644361, + 3.14051, + 0.0446999, + 0.0644361, + 0.0446999, + 0.0644361, + 0.0120662, + 0.0446999, + 3.14051, + 0.0644361, + 0.0644361, + 0.0446999, + 0.0644361, + 0.0446999, + 0.0446999, + 0.0120662 + }; + + Vector estimated_error_per_cell(&values[0], + &values[0] + 120); + + deallog << "Number of cells before refinement: " << triangulation.n_active_cells() + << std::endl; + + parallel::distributed::GridRefinement:: + refine_and_coarsen_fixed_fraction (triangulation, + estimated_error_per_cell, + 0.3, 0.0); + triangulation.execute_coarsening_and_refinement (); + + deallog << "Number of cells after refinement: " << triangulation.n_active_cells() + << std::endl; + +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test(); + } + else + test(); + +} diff --git a/tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output b/tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output new file mode 100644 index 0000000000..5bf1b8a10e --- /dev/null +++ b/tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output @@ -0,0 +1,3 @@ + +DEAL:0::Number of cells before refinement: 120 +DEAL:0::Number of cells after refinement: 144 -- 2.39.5