From 33c0c819b52472daaecaa451897cbd901dd02119 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Mon, 9 Dec 2013 11:35:24 +0000
Subject: [PATCH] Need to disable CellSimilarity when the number of threads is
 greater than one because otherwise we get different mapping data (on the
 order of roundoff) on cells that are translations of each other. When using
 WorkStream, this might result in different matrices that are eventually
 generated

git-svn-id: https://svn.dealii.org/trunk@31950 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/news/8.0.0-vs-8.1.0.h | 11 ++++++++++
 deal.II/source/fe/fe_values.cc    | 34 +++++++++++++++++++++++--------
 2 files changed, 37 insertions(+), 8 deletions(-)

diff --git a/deal.II/doc/news/8.0.0-vs-8.1.0.h b/deal.II/doc/news/8.0.0-vs-8.1.0.h
index acaaa648e3..4afef43f2f 100644
--- a/deal.II/doc/news/8.0.0-vs-8.1.0.h
+++ b/deal.II/doc/news/8.0.0-vs-8.1.0.h
@@ -259,6 +259,17 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: Since the introduction of ThreadLocalStorage in version 8.0, the
+  way in which FEValues objects visit cells in a parallel assembly loop is no
+  longer deterministic. Therefore, the detection of CellSimilarity that can
+  speed up computations of certain geometric quantities (shape gradients) on
+  cells that are translations is disabled when the number of threads is
+  greater than one. This produces somewhat slower code (usually not more than
+  a few percent) but ensures exact reproducibility of results.
+  <br>
+  (Martin Kronbichler, Wolfgang Bangerth, 2013/12/09)
+  </li>
+
   <li> Fixed: Several functions in namespace GridTools were not instantiated
   for parallel::distributed::Triangulation objects. This is now fixed.
   <br>
diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc
index a027938320..639d937fd9 100644
--- a/deal.II/source/fe/fe_values.cc
+++ b/deal.II/source/fe/fe_values.cc
@@ -15,6 +15,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/multithread_info.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
@@ -3256,15 +3257,33 @@ void
 FEValuesBase<dim,spacedim>::check_cell_similarity
 (const typename Triangulation<dim,spacedim>::cell_iterator &cell)
 {
-  // case that there has not been any cell
-  // before
+  // Unfortunately, the detection of simple geometries with CellSimilarity is
+  // sensitive to the first cell detected. When doing this with multiple
+  // threads, each thread will get its own scratch data object with an
+  // FEValues object in the implementation framework from late 2013, which is
+  // initialized to the first cell the thread sees. As this number might
+  // different between different runs (after all, the tasks are scheduled
+  // dynamically onto threads), this slight deviation leads to difference in
+  // roundoff errors that propagate through the program. Therefore, we need to
+  // disable CellSimilarity in case there is more than one thread in the
+  // problem. This will likely not affect many MPI test cases as there
+  // multithreading is disabled on default, but in many other situations
+  // because we rarely explicitly set the number of threads.
+  //
+  // TODO: Is it reasonable to introduce a flag "unsafe" in the constructor of
+  // FEValues to re-enable this feature?
+  if (multithread_info.n_threads() > 1)
+    {
+      cell_similarity = CellSimilarity::none;
+      return;
+    }
+
+  // case that there has not been any cell before
   if (this->present_cell.get() == 0)
     cell_similarity = CellSimilarity::none;
   else
-    // in MappingQ, data can have been
-    // modified during the last call. Then,
-    // we can't use that data on the new
-    // cell.
+    // in MappingQ, data can have been modified during the last call. Then, we
+    // can't use that data on the new cell.
     if (cell_similarity == CellSimilarity::invalid_next_cell)
       cell_similarity = CellSimilarity::none;
     else
@@ -3282,8 +3301,7 @@ FEValuesBase<dim,spacedim>::check_cell_similarity
           != cell->direction_flag() )
         cell_similarity =  CellSimilarity::inverted_translation;
     }
-  // TODO: here, one could implement other
-  // checks for similarity, e.g. for
+  // TODO: here, one could implement other checks for similarity, e.g. for
   // children of a parallelogram.
 }
 
-- 
2.39.5