Use position hash as pseudo-random seed

geodynamics · Jun 24, 2017 · a3e0f29 · a3e0f29
1 parent b1b3b2b
commit a3e0f29
Show file tree

Hide file tree

Showing 11 changed files with 111 additions and 108 deletions.
diff --git a/doc/modules/changes/20170609_gassmoeller b/doc/modules/changes/20170609_gassmoeller
@@ -1,4 +1,4 @@
-New: A new initial temperature condition was added that adds a random
+New: A new initial temperature condition was added that adds random
 noise of user prescribed magnitude to the initial temperature field.
 The same type of plugin was added for initial compositions.
 <br>

diff --git a/include/aspect/initial_composition/random_perturbation.h b/include/aspect/initial_composition/random_perturbation.h
@@ -74,9 +74,13 @@ namespace aspect
         double magnitude;
 
         /**
-         * Random number generator.
+         * Whether to use a random seed for the random
+         * number generator. This parameter controls whether
+         * this plugin generates different or identical
+         * perturbations for subsequent model runs of
+         * the same setup.
          */
-        mutable boost::mt19937 random_number_generator;
+        bool use_random_seed;
     };
   }
 }

diff --git a/include/aspect/initial_temperature/random_perturbation.h b/include/aspect/initial_temperature/random_perturbation.h
@@ -73,9 +73,13 @@ namespace aspect
         double magnitude;
 
         /**
-         * Random number generator.
+         * Whether to use a random seed for the random
+         * number generator. This parameter controls whether
+         * this plugin generates different or identical
+         * perturbations for subsequent model runs of
+         * the same setup.
          */
-        mutable boost::mt19937 random_number_generator;
+        bool use_random_seed;
     };
   }
 }

diff --git a/source/initial_composition/random_perturbation.cc b/source/initial_composition/random_perturbation.cc
@@ -21,23 +21,44 @@
 
 #include <aspect/initial_composition/random_perturbation.h>
 
+#include <boost/functional/hash.hpp>
 #include <time.h>
 
 namespace aspect
 {
   namespace InitialComposition
   {
+    namespace
+    {
+      template<int dim>
+      std::size_t point_hash(const Point<dim> &position)
+      {
+        std::size_t hash;
+
+        for (unsigned int i = 0; i < dim; ++i)
+          boost::hash_combine(hash,position[i]);
+
+        return hash;
+      }
+    }
+
     template <int dim>
     double
     RandomPerturbation<dim>::
-    initial_composition (const Point<dim> &/*position*/,
-                         const unsigned int /*compositional_index*/) const
+    initial_composition (const Point<dim> &position,
+                         const unsigned int compositional_index) const
     {
-      // Uniform distribution on the interval [0,1]. This
-      // will be used to generate the random composition perturbation.
-      boost::uniform_01<double> uniform_distribution_01;
+      std::size_t seed = point_hash(position);
+      boost::hash_combine(seed,compositional_index);
+      if (use_random_seed)
+        boost::hash_combine(seed,time(NULL));
+
+      boost::mt19937 random_number_generator(seed);
 
-      return magnitude * 2.0 * (uniform_distribution_01(random_number_generator) - 0.5);
+      // Uniform distribution on the interval [-magnitude,magnitude). This
+      // will be used to generate the random composition perturbation.
+      boost::random::uniform_real_distribution<double> uniform_distribution(-magnitude,magnitude);
+      return uniform_distribution(random_number_generator);
     }
 
     template <int dim>
@@ -56,14 +77,8 @@ namespace aspect
                              "Whether to use a random seed for the random "
                              "number generator. This parameter controls whether "
                              "this plugin generates different or identical "
-                             "perturbations for different model runs.");
-          prm.declare_entry ("Random seed", "5432",
-                             Patterns::Integer(0),
-                             "The initial seed for the random perturbation. If "
-                             "'Use random seed' is set to 'false' model runs "
-                             "with identical seed will lead to identical "
-                             "perturbations, while model runs with different "
-                             "seeds will generate different perturbations.");
+                             "perturbations for subsequent model runs of"
+                             "the same setup.");
         }
         prm.leave_subsection ();
       }
@@ -80,11 +95,7 @@ namespace aspect
         prm.enter_subsection("Random perturbation");
         {
           magnitude = prm.get_double ("Magnitude");
-
-          if (prm.get_bool("Use random seed"))
-            random_number_generator.seed(time(NULL));
-          else
-            random_number_generator.seed(prm.get_integer("Random seed"));
+          use_random_seed = prm.get_bool("Use random seed");
         }
         prm.leave_subsection ();
       }

diff --git a/source/initial_temperature/random_perturbation.cc b/source/initial_temperature/random_perturbation.cc
@@ -21,22 +21,43 @@
 
 #include <aspect/initial_temperature/random_perturbation.h>
 
+#include <boost/functional/hash.hpp>
 #include <time.h>
 
 namespace aspect
 {
   namespace InitialTemperature
   {
+    namespace
+    {
+      template<int dim>
+      std::size_t point_hash(const Point<dim> &position)
+      {
+        std::size_t hash;
+
+        for (unsigned int i = 0; i < dim; ++i)
+          boost::hash_combine(hash,position[i]);
+
+        return hash;
+      }
+    }
+
     template <int dim>
     double
     RandomPerturbation<dim>::
-    initial_temperature (const Point<dim> &/*position*/) const
+    initial_temperature (const Point<dim> &position) const
     {
-      // Uniform distribution on the interval [0,1]. This
-      // will be used to generate the random temperature perturbation.
-      boost::uniform_01<double> uniform_distribution_01;
+      std::size_t seed = point_hash(position);
+
+      if (use_random_seed)
+        boost::hash_combine(seed,time(NULL));
 
-      return magnitude * 2.0 * (uniform_distribution_01(random_number_generator) - 0.5);
+      boost::mt19937 random_number_generator(seed);
+
+      // Uniform distribution on the interval [-magnitude,magnitude). This
+      // will be used to generate the random temperature perturbation.
+      boost::random::uniform_real_distribution<double> uniform_distribution(-magnitude,magnitude);
+      return uniform_distribution(random_number_generator);
     }
 
     template <int dim>
@@ -56,13 +77,6 @@ namespace aspect
                              "number generator. This parameter controls whether "
                              "this plugin generates different or identical "
                              "perturbations for different model runs.");
-          prm.declare_entry ("Random seed", "5432",
-                             Patterns::Integer(0),
-                             "The initial seed for the random perturbation. If "
-                             "'Use random seed' is set to 'false' model runs "
-                             "with identical seed will lead to identical "
-                             "perturbations, while model runs with different "
-                             "seeds will generate different perturbations.");
         }
         prm.leave_subsection ();
       }
@@ -79,11 +93,7 @@ namespace aspect
         prm.enter_subsection("Random perturbation");
         {
           magnitude = prm.get_double ("Magnitude");
-
-          if (prm.get_bool("Use random seed"))
-            random_number_generator.seed(time(NULL));
-          else
-            random_number_generator.seed(prm.get_integer("Random seed"));
+          use_random_seed = prm.get_bool("Use random seed");
         }
         prm.leave_subsection ();
       }

diff --git a/tests/random_composition_perturbation.prm b/tests/random_composition_perturbation.prm
@@ -1,69 +1,56 @@
-# Test to check the status of the random initial temperature perturbation
+# Test to check the status of the random initial composition perturbation
 # plugin.
 
-set End time                               = 0
+subsection Model settings
+  set Tangential velocity boundary indicators = left, right, bottom, top
+end
 
-subsection Boundary temperature model
-  set List of model names = box
+subsection Mesh refinement
+  set Initial global refinement = 3
 end
 
+subsection Postprocess
+  set List of postprocessors = composition statistics
+end
 
-subsection Gravity model
-  set Model name = vertical
+subsection Compositional fields
+  set Number of fields = 2
 end
 
+subsection Material model
+  set Model name = simpler
+end
 
+set End time = 0
 subsection Geometry model
   set Model name = box
-
-  subsection Box
-    set X extent = 1
-    set Y extent = 1
-  end
 end
 
+subsection Gravity model
+  set Model name = vertical
+end
 
 subsection Initial temperature model
   set List of model names = function
-
   subsection Function
     set Function expression = 0.5
   end
-end
 
-subsection Compositional fields
-  set Number of fields = 2
 end
 
 subsection Initial composition model
   set List of model names = function, random perturbation
-
   subsection Random perturbation
     set Magnitude = 0.2
   end
 
   subsection Function
     set Function expression = 0.5;0.0
   end
-end
-
 
-subsection Material model
-  set Model name = simpler
 end
 
-
-subsection Mesh refinement
-  set Initial adaptive refinement        = 0
-  set Initial global refinement          = 3
-end
-
-
-subsection Model settings
-  set Tangential velocity boundary indicators = left, right, bottom, top
-end
-
-subsection Postprocess
-  set List of postprocessors = composition statistics
+subsection Boundary temperature model
+  set List of model names = box
 end
 
diff --git a/tests/random_composition_perturbation/screen-output b/tests/random_composition_perturbation/screen-output
@@ -10,7 +10,7 @@ Number of degrees of freedom: 1,526 (578+81+289+289+289)
    Solving Stokes system... 2+0 iterations.
 
    Postprocessing:
-     Compositions min/max/mass: 0.3001/0.6984/0.4938 // -0.2/0.1996/-0.01184
+     Compositions min/max/mass: 0.3002/0.6997/0.5069 // -0.199/0.1999/-0.007404
 
 Termination requested by criterion: end time
 

diff --git a/tests/random_composition_perturbation/statistics b/tests/random_composition_perturbation/statistics
@@ -17,4 +17,4 @@
 # 17: Minimal value for composition C_2
 # 18: Maximal value for composition C_2
 # 19: Global mass for composition C_2
-0 0.000000000000e+00 0.000000000000e+00 64 659 289 578 0 0 0 2 3 2 3.00090365e-01 6.98414874e-01 4.93781411e-01 -1.99952361e-01 1.99564844e-01 -1.18366501e-02 
+0 0.000000000000e+00 0.000000000000e+00 64 659 289 578 0 0 0 2 3 2 3.00222459e-01 6.99744700e-01 5.06886136e-01 -1.99012748e-01 1.99867409e-01 -7.40357014e-03 
diff --git a/tests/random_temperature_perturbation.prm b/tests/random_temperature_perturbation.prm
@@ -1,57 +1,44 @@
 # Test to check the status of the random initial temperature perturbation
 # plugin.
 
-set End time                               = 0
-
-subsection Boundary temperature model
-  set List of model names = box
+subsection Model settings
+  set Tangential velocity boundary indicators = left, right, bottom, top
 end
 
+subsection Mesh refinement
+  set Initial global refinement = 3
+end
 
-subsection Gravity model
-  set Model name = vertical
+subsection Postprocess
+  set List of postprocessors = temperature statistics
 end
 
+subsection Material model
+  set Model name = simpler
+end
 
+set End time = 0
 subsection Geometry model
   set Model name = box
-
-  subsection Box
-    set X extent = 1
-    set Y extent = 1
-  end
 end
 
+subsection Gravity model
+  set Model name = vertical
+end
 
 subsection Initial temperature model
   set List of model names = function, random perturbation
-
   subsection Random perturbation
     set Magnitude = 0.2
   end
 
   subsection Function
     set Function expression = 0.5
   end
-end
-
-
-subsection Material model
-  set Model name = simpler
-end
-
-
-subsection Mesh refinement
-  set Initial adaptive refinement        = 0
-  set Initial global refinement          = 3
-end
-
 
-subsection Model settings
-  set Tangential velocity boundary indicators = left, right, bottom, top
 end
 
-subsection Postprocess
-  set List of postprocessors = temperature statistics
+subsection Boundary temperature model
+  set List of model names = box
 end
 
diff --git a/tests/random_temperature_perturbation/screen-output b/tests/random_temperature_perturbation/screen-output
@@ -5,10 +5,10 @@ Number of degrees of freedom: 948 (578+81+289)
 *** Timestep 0:  t=0 years
    Solving temperature system... 0 iterations.
    Rebuilding Stokes preconditioner...
-   Solving Stokes system... 17+0 iterations.
+   Solving Stokes system... 18+0 iterations.
 
    Postprocessing:
-     Temperature min/avg/max: 0.3001 K, 0.5056 K, 0.6973 K
+     Temperature min/avg/max: 0.3034 K, 0.4927 K, 0.6999 K
 
 Termination requested by criterion: end time
 

diff --git a/tests/random_temperature_perturbation/statistics b/tests/random_temperature_perturbation/statistics
@@ -11,4 +11,4 @@
 # 11: Minimal temperature (K)
 # 12: Average temperature (K)
 # 13: Maximal temperature (K)
-0 0.000000000000e+00 0.000000000000e+00 64 659 289 0 17 18 17 3.00129639e-01 5.05619087e-01 6.97305210e-01 
+0 0.000000000000e+00 0.000000000000e+00 64 659 289 0 18 19 18 3.03379361e-01 4.92662594e-01 6.99863972e-01