Commits

Cameron Hummels  committed 17852b0

Adding mersennes twister random number generator to KH test for reproducible results in old method (with random y-velocity perturbations).

  • Participants
  • Parent commits dcfbebb

Comments (0)

Files changed (3)

File src/enzo/Grid.h

                        float KHPerturbationAmplitude,
                        float KHInnerVx, float KHOuterVx,
                        float KHInnerPressure,
-                       float KHOuterPressure);
+                       float KHOuterPressure,
+                       int   KHRandomSeed);
 
   /* Initialize a grid for the KH instability problem including a ramp. */
 

File src/enzo/Grid_KHInitializeGrid.C

 /  modified1:  Alexei Kritsuk, December 2004.  
 /  modified2:  Gregg Dobrowalski, Feb 2005.  
 /  modified3:  Alexei Kritsuk, April 2005. v_y perturbations + more parameters.
+/  modified4:  Cameron Hummels, December 2013. mersennes twister RNG
 /
 /  PURPOSE: Sets the field variables in the domain.
 /          
 #include "ExternalBoundary.h"
 #include "Grid.h"
 
+/* including mersennes twister random number generator; allows
+   reproducible behavior with same seed */
+
+void mt_init(unsigned_int seed);
+
+unsigned_long_int mt_random();
+
 int grid::KHInitializeGrid(float KHInnerDensity, 
                            float KHInnerInternalEnergy,
                            float KHOuterInternalEnergy,
                            float KHPerturbationAmplitude,
                            float KHInnerVx, float KHOuterVx,
                            float KHInnerPressure,
-                           float KHOuterPressure)
+                           float KHOuterPressure,
+                           int   KHRandomSeed)
 {
 
   if (ProcessorNumber != MyProcessorNumber)
   /* declarations */
 
   int size = 1, dim;
+  float rand_x, rand_y;
+  unsigned_long_int random_int;
+
+  /* initialize random number generator from seed */
+  mt_init(((unsigned_int) KHRandomSeed)); 
+
   for (dim = 0; dim < GridRank; dim++)
     size *= GridDimension[dim];
 
     index = i % GridDimension[0];
     jndex = (i-index)/GridDimension[0];
 
-    BaryonField[3][i] = KHPerturbationAmplitude * 
-                        ((float)rand()/(float)(RAND_MAX) - 0.5); //AK
+    rand_x = (float)(mt_random()%32768)/(32768.0) - 0.5;
+    rand_y = (float)(mt_random()%32768)/(32768.0) - 0.5;
+
+    BaryonField[3][i] = KHPerturbationAmplitude * rand_y;
     BaryonField[1][i] = POW(BaryonField[3][i],2) / 2.0;
       
     if (jndex >= GridDimension[1]/4 && jndex < 3*GridDimension[1]/4) { //AK
       BaryonField[0][i]  = KHInnerDensity;
-      BaryonField[2][i]  = KHInnerVx + 
-                           KHPerturbationAmplitude * 
-                           ((float)rand()/(float)(RAND_MAX) - 0.5);
+      BaryonField[2][i]  = KHInnerVx + KHPerturbationAmplitude * rand_x;
       BaryonField[1][i] += KHInnerInternalEnergy + 
                            POW(BaryonField[2][i],2) / 2.0;
     }
     else {
-      BaryonField[2][i]  = KHOuterVx + 
-                           KHPerturbationAmplitude * 
-                           ((float)rand()/(float)(RAND_MAX) - 0.5);
+      BaryonField[2][i]  = KHOuterVx + KHPerturbationAmplitude * rand_x;
       BaryonField[1][i] += KHOuterInternalEnergy + 
                            POW(BaryonField[2][i],2) / 2.0;
     }

File src/enzo/KHInitialize.C

   float KHBulkVelocity          = 0.0;
   int   KHRamp                  = 1;    // Convergent ICs with Ramp
   float KHRampWidth             = 0.05;
+  int   KHRandomSeed            = 123456789;
 
   float KHInnerInternalEnergy, KHOuterInternalEnergy;
 
     ret += sscanf(line, "KHBulkVelocity  = %"FSYM, &KHBulkVelocity);
     ret += sscanf(line, "KHRamp = %"ISYM, &KHRamp);
     ret += sscanf(line, "KHRampWidth     = %"FSYM, &KHRampWidth);
+    ret += sscanf(line, "KHRandomSeed    = %"ISYM, &KHRandomSeed);
 
     /* if the line is suspicious, issue a warning */
 
                                            KHInnerVelocity[0], 
                                            KHOuterVelocity[0],
                                            KHInnerPressure,
-                                           KHOuterPressure)
+                                           KHOuterPressure,
+                                           KHRandomSeed)
         == FAIL) {
       ENZO_FAIL("Error in KHInitializeGrid.");
     }