Commits

Cameron Hummels  committed 75e0662

Modifying KH ICs to use convention of BaryonField[TENum]. Removing reference to level in Grid_KHInitializeGridRamp.

  • Participants
  • Parent commits 17852b0

Comments (0)

Files changed (4)

File src/enzo/Grid.h

                            float KHInnerVx, float KHOuterVx,
                            float KHInnerPressure,
                            float KHOuterPressure,
-                           float KHRampWidth,
-                           int   level);
+                           float KHRampWidth);
 
   /* Initialize a grid and set boundary for the 2D/3D Noh problem. */
 

File src/enzo/Grid_KHInitializeGrid.C

   for (dim = 0; dim < GridRank; dim++)
     size *= GridDimension[dim];
 
-    /* set fields in "inner" region: .25 < y < .75 */
+  /* Find fields: density, total energy, velocity1-3. */
+
+  int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num;
+  if (this-IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
+                                      Vel3Num, TENum) == FAIL) {
+    ENZO_FAIL("Error in IdentifyPhysicalQuantities.\n");
+  }
+
+  /* set fields in "inner" region: .25 < y < .75 */
 
   int index, jndex, i;
   for (i = 0; i < size; i++) {
     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;
+    BaryonField[Vel2Num][i] = KHPerturbationAmplitude * rand_y;
+    BaryonField[TENum][i] = POW(BaryonField[Vel2Num][i],2) / 2.0;
       
     if (jndex >= GridDimension[1]/4 && jndex < 3*GridDimension[1]/4) { //AK
-      BaryonField[0][i]  = KHInnerDensity;
-      BaryonField[2][i]  = KHInnerVx + KHPerturbationAmplitude * rand_x;
-      BaryonField[1][i] += KHInnerInternalEnergy + 
-                           POW(BaryonField[2][i],2) / 2.0;
+      BaryonField[DensNum][i]  = KHInnerDensity;
+      BaryonField[Vel1Num][i]  = KHInnerVx + KHPerturbationAmplitude * rand_x;
+      BaryonField[TENum][i] += KHInnerInternalEnergy + 
+                               POW(BaryonField[Vel1Num][i],2) / 2.0;
     }
     else {
-      BaryonField[2][i]  = KHOuterVx + KHPerturbationAmplitude * rand_x;
-      BaryonField[1][i] += KHOuterInternalEnergy + 
-                           POW(BaryonField[2][i],2) / 2.0;
+      BaryonField[Vel1Num][i]  = KHOuterVx + KHPerturbationAmplitude * rand_x;
+      BaryonField[TENum][i] += KHOuterInternalEnergy + 
+                               POW(BaryonField[Vel1Num][i],2) / 2.0;
     }
   }
       
   if (debug) {
     fprintf(stderr, "GKHIG: BF[2][0]:%"FSYM" BF[2][size/2]: %"FSYM"\n", 
-            BaryonField[2][0], BaryonField[2][size/2]);
+            BaryonField[Vel1Num][0], BaryonField[Vel1Num][size/2]);
 
     fprintf(stderr, "GKHIG: BF[3][0]:%"FSYM" BF[3][size/2]: %"FSYM"\n", 
-            BaryonField[3][0], BaryonField[3][size/2]);
+            BaryonField[Vel2Num][0], BaryonField[Vel2Num][size/2]);
   }
   return SUCCESS;
 }

File src/enzo/Grid_KHInitializeGridRamp.C

                                float KHInnerVx, float KHOuterVx,
                                float KHInnerPressure,
                                float KHOuterPressure,
-                               float KHRampWidth,
-                               int   level)
+                               float KHRampWidth)
 {
 
   if (ProcessorNumber != MyProcessorNumber)
   for (dim = 0; dim < GridRank; dim++)
     size *= GridDimension[dim];
 
+  /* Find fields: density, total energy, velocity1-3. */
+
+  int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num;
+  if (this-IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num,
+                                      Vel3Num, TENum) == FAIL) {
+    ENZO_FAIL("Error in IdentifyPhysicalQuantities.\n");
+  }
+
   int index, jndex, i, j, k, n = 0;
   FLOAT x, y, KHRampWidth2;
   KHRampWidth2 = KHRampWidth*KHRampWidth;
     y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j];
 
     /* everywhere set the sinusoidal perturbation in v_y */
-    BaryonField[3][n] = KHPerturbationAmplitude * sin(4.0 * M_PI * x);
+    BaryonField[Vel2Num][n] = KHPerturbationAmplitude * sin(4.0 * M_PI * x);
 
     /* set density and velocity fields in "inner region" */
     if (y > 0.25 && y < 0.75) { 
-      BaryonField[0][n]  = KHInnerDensity;
-      BaryonField[2][n]  = KHInnerVx;
+      BaryonField[DensNum][n]  = KHInnerDensity;
+      BaryonField[Vel1Num][n]  = KHInnerVx;
 
     /* modify the top half of the inner fluid to account for ramp */
       if (y >= 0.5) {
-        BaryonField[0][n] -= exp( (-0.5/KHRampWidth2)*
+        BaryonField[DensNum][n] -= exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.75 - 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
-        BaryonField[2][n] -= exp( (-0.5/KHRampWidth2)*
+        BaryonField[Vel1Num][n] -= exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.75 - 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
       } else {
 
     /* modify the bottom half of the inner fluid to account for ramp */
-        BaryonField[0][n] -= exp( (-0.5/KHRampWidth2)*
+        BaryonField[DensNum][n] -= exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.25 + 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
-        BaryonField[2][n] -= exp( (-0.5/KHRampWidth2)*
+        BaryonField[Vel1Num][n] -= exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.25 + 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
       }
 
     /* set the energy in the "inner region" to be thermal + kinetic */
-      BaryonField[1][n] = KHInnerPressure/((Gamma - 1.0)*BaryonField[0][n]) + 
-                       (POW(BaryonField[2][n],2) +
-                        POW(BaryonField[3][n],2)) / 2.0;
+      BaryonField[TENum][n] = KHInnerPressure/((Gamma - 1.0)*BaryonField[DensNum][n]) + 
+                       (POW(BaryonField[Vel1Num][n],2) +
+                        POW(BaryonField[Vel2Num][n],2)) / 2.0;
 
     /* set density and velocity fields in "outer region" */
     } else {
-      BaryonField[0][n]  = KHOuterDensity;
-      BaryonField[2][n]  = KHOuterVx;
+      BaryonField[DensNum][n]  = KHOuterDensity;
+      BaryonField[Vel1Num][n]  = KHOuterVx;
 
     /* modify the top half of the outer fluid to account for ramp */
       if (y >= 0.5) {
-        BaryonField[0][n] += exp( (-0.5/KHRampWidth2)*
+        BaryonField[DensNum][n] += exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.75 + 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
-        BaryonField[2][n] += exp( (-0.5/KHRampWidth2)*
+        BaryonField[Vel1Num][n] += exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.75 + 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
       } else {
 
     /* modify the bottom half of the outer fluid to account for ramp */
-        BaryonField[0][n] += exp( (-0.5/KHRampWidth2)*
+        BaryonField[DensNum][n] += exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.25 - 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
-        BaryonField[2][n] += exp( (-0.5/KHRampWidth2)*
+        BaryonField[Vel1Num][n] += exp( (-0.5/KHRampWidth2)*
                                   pow(y-0.25 - 
                                       sqrt(-2.0*KHRampWidth2*log(0.5)),2));
       }
 
     /* set the energy in the "outer region" to be thermal + kinetic */
-      BaryonField[1][n] = KHOuterPressure/((Gamma - 1.0)*BaryonField[0][n]) + 
-                       (POW(BaryonField[2][n],2) +
-                        POW(BaryonField[3][n],2)) / 2.0;
+      BaryonField[TENum][n] = KHOuterPressure/((Gamma - 1.0)*BaryonField[DensNum][n]) + 
+                       (POW(BaryonField[Vel1Num][n],2) +
+                        POW(BaryonField[Vel2Num][n],2)) / 2.0;
     }
   }
   return SUCCESS;

File src/enzo/KHInitialize.C

                                                KHOuterVelocity[0],
                                                KHInnerPressure,
                                                KHOuterPressure,
-                                               KHRampWidth,
-                                               0) 
+                                               KHRampWidth)
         == FAIL) {
       ENZO_FAIL("Error in KHInitializeGridRamp.");
     }
                                                    KHOuterVelocity[0],
                                                    KHInnerPressure,
                                                    KHOuterPressure,
-                                                   KHRampWidth, 
-                                                   level+1) 
+                                                   KHRampWidth)
               == FAIL) {
             ENZO_FAIL("Error in KHInitializeGridRamp.");
           }