1. Britton Smith
  2. enzo-dev-coolingtime-rh

Commits

Britton Smith  committed ddaf458

Changed dynamic hierarchy rebuilding to use a time instead of a number
of cycles to go without rebuilding.
Changed Grid_ComputeConductionTimeStep to return the timestep instead
of setting a pointer value.

  • Participants
  • Parent commits 7e4b29c
  • Branches week-of-code

Comments (0)

Files changed (9)

File src/enzo/EvolveLevel.C

View file
     SetLevelTimeStep(Grids, NumberOfGrids, level, 
         &dtThisLevelSoFar[level], &dtThisLevel[level], dtLevelAbove);
 
+    TimeSinceRebuildHierarchy[level] += dtThisLevel[level];
+
     /* If StarFormationOncePerRootGridTimeStep, stars are only created
     once per root grid time step and only on MaximumRefinementLevel
     grids. The following sets the MakeStars flag for all

File src/enzo/Grid.h

View file
 /* Member functions for dealing with thermal conduction */
    int ComputeHeat(float dedt[]);	     /* Compute Heat */
    int ConductHeat();			     /* Conduct Heat */
-   int ComputeConductionTimeStep(float &dt); /* Estimate conduction time-step */
+   float ComputeConductionTimeStep(); /* Estimate conduction time-step */
 
 /* Baryons: Copy current solution to Old solution (returns success/fail)
     (for step #16) */

File src/enzo/Grid_ComputeConductionTimeStep.C

View file
 /  on a given grid patch.
 /
 /  RETURNS:
-/    SUCCESS or FAIL
+/    The conduction timestep.
 /
 ************************************************************************/
 
 	      float *VelocityUnits, double *MassUnits, FLOAT Time);
 
 // Member functions
-int grid::ComputeConductionTimeStep (float &dt) {
-  if (ProcessorNumber != MyProcessorNumber) {return SUCCESS;}
-  if (NumberOfBaryonFields == 0) {return SUCCESS;}
+float grid::ComputeConductionTimeStep () {
+  if (ProcessorNumber != MyProcessorNumber) {return huge_number;}
+  if (NumberOfBaryonFields == 0) {return huge_number;}
   this->DebugCheck("ComputeConductionTimeStep");
 
   // Some locals
   FLOAT a = 1.0, dadt;
   double MassUnits = 1.0;
   float *rho;
-  float dt_est;
+  float dt, dt_est;
   double all_units;
 
   float SpitzerFraction;
     SpitzerFraction = AnisotropicConductionSpitzerFraction;
   }
   else {
-    return SUCCESS;
+    return huge_number;
   }
 
   int size = 1, grid_index, right_side_index; 
     ( 6.0e-7 * SpitzerFraction * mh * TimeUnits );
   
   dt *= float(all_units);
+  dt *= ConductionCourantSafetyNumber;  // for stability, this has to be < 0.5
 
   delete [] Temp;
 
-  return SUCCESS;
+  return dt;
  
 }

File src/enzo/Grid_ComputeTimeStep.C

View file
   /* 5) Calculate minimum dt due to thermal conduction. */
 
   if(IsotropicConduction || AnisotropicConduction){
-    if (this->ComputeConductionTimeStep(dtConduction) == FAIL) {
-      fprintf(stderr, "Error in ComputeConductionTimeStep.\n");
-      return FAIL;
-    }
-    dtConduction *= ConductionCourantSafetyNumber;  // for stability
-    dtConduction *= float(DEFAULT_GHOST_ZONES);     // for subcycling 
+    dtConduction = this->ComputeConductionTimeStep();
+    dtConduction *= float(DEFAULT_GHOST_ZONES);  // for subcycling
   }
 
   /* 6) GasDrag time step */

File src/enzo/Grid_ConductHeat.C

View file
     }
 
     // compute this subcycle timestep
-    if (this->ComputeConductionTimeStep(dtSubcycle) == FAIL) {
-      ENZO_FAIL("Error in ComputeConductionTimeStep.");
-    }
-    dtSubcycle *= ConductionCourantSafetyNumber;  // for stability, this has to be < 0.5
+    dtSubcycle = this->ComputeConductionTimeStep();
 
     // make sure we don't extend past dtFixed
     dtSubcycle = min(dtSubcycle, dtFixed-dtSoFar);

File src/enzo/RebuildHierarchy.C

View file
     return SUCCESS;
   }
 
+  if (ConductionDynamicRebuildHierarchy) {
+    if (TimeSinceRebuildHierarchy[level] < dtRebuildHierarchy[level]) {
+      return SUCCESS;
+    }
+    else {
+      for (int i = level;i <= MaximumRefinementLevel;i++) {
+        dtRebuildHierarchy[i] = -1.0;
+        TimeSinceRebuildHierarchy[i] = 0.0;
+      }
+    }
+  }
+
   double tt0, tt1, tt2, tt3;
  
   /* declarations */

File src/enzo/SetDefaultGlobalValues.C

View file
 
   IsotropicConduction = FALSE;
   AnisotropicConduction = FALSE;
-  IsotropicConductionSpitzerFraction = 1.0;
-  AnisotropicConductionSpitzerFraction = 1.0;
+  IsotropicConductionSpitzerFraction = 0.0;
+  AnisotropicConductionSpitzerFraction = 0.0;
   ConductionCourantSafetyNumber = 0.5;
 
   PythonTopGridSkip                = 0;
 
   /* Some stateful variables for EvolveLevel */
   for(i = 0; i < MAX_DEPTH_OF_HIERARCHY; i++) {
-    LevelCycleCount[i] = 0;
-    LevelSubCycleCount[i] = 0;
+    LevelCycleCount[i] = LevelSubCycleCount[i] = 0;
+    dtRebuildHierarchy[i] = -1.0;
+    TimeSinceRebuildHierarchy[i] = 0.0;
     dtThisLevelSoFar[i] = dtThisLevel[i] = 0.0;
   }
 

File src/enzo/SetLevelTimeStep.C

View file
  
   } else {
  
+    /* Calculate timestep without conduction and get conduction separately later. */
+
+    int my_isotropic_conduction = IsotropicConduction;
+    int my_anisotropic_conduction = AnisotropicConduction;
+    int dynamic_hierarchy_rebuild = ConductionDynamicRebuildHierarchy &&
+      dtRebuildHierarchy[level] <= 0.0;
+
+    if (dynamic_hierarchy_rebuild) {
+      IsotropicConduction = AnisotropicConduction = FALSE;      
+    }
+
     /* Compute the mininum timestep for all grids. */
  
     *dtThisLevel = huge_number;
     }
     *dtThisLevel = CommunicationMinValue(*dtThisLevel);
 
-    dtActual = *dtThisLevel;
-
-    /* Compute timestep without conduction and use to set the number 
+    /* Compute conduction timestep and use to set the number 
        of iterations without rebuiding the hierarchy. */
 
-    if (ConductionDynamicRebuildHierarchy && LevelSubCycleCount[level] == 0) {
-      float dt_no_conduction, my_dt_no_conduction;
-      int my_isotropic_conduction = IsotropicConduction;
-      int my_anisotropic_conduction = AnisotropicConduction;
-      IsotropicConduction = AnisotropicConduction = FALSE;
+    if (dynamic_hierarchy_rebuild) {
 
-      dt_no_conduction = huge_number;
-      for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
-        my_dt_no_conduction = Grids[grid1]->GridData->ComputeTimeStep();
-        dt_no_conduction = min(dt_no_conduction, my_dt_no_conduction);
-      }
-      dt_no_conduction = CommunicationMinValue(dt_no_conduction);
-      int my_cycle_skip = (int) (dt_no_conduction / dtActual);
-      RebuildHierarchyCycleSkip[level] = max(1, my_cycle_skip);
-      if (debug) fprintf(stderr, "AdaptiveRebuildHierarchyCycleSkip[%"ISYM"] = %"ISYM"\n",
-                         level, RebuildHierarchyCycleSkip[level]);
-
+      /* Return conduction parameters to original values. */
       IsotropicConduction = my_isotropic_conduction;
       AnisotropicConduction = my_anisotropic_conduction;
+
+      float dt_conduction;
+      dt_conduction = huge_number;
+      for (grid1 = 0; grid1 < NumberOfGrids; grid1++) {
+        dt_conduction = min(dt_conduction, 
+                            Grids[grid1]->GridData->ComputeConductionTimeStep());
+      }
+      dt_conduction = CommunicationMinValue(dt_conduction);
+      dt_conduction *= float(DEFAULT_GHOST_ZONES);  // for subcycling
+
+      int my_cycle_skip = max(1, (int) (*dtThisLevel / dt_conduction));
+      dtRebuildHierarchy[level] = *dtThisLevel;
+      if (debug)
+        fprintf(stderr, "Conduction dt[%"ISYM"] = %"GSYM", will rebuild hierarchy in about %"ISYM" cycles.\n",
+                         level, dt_conduction, my_cycle_skip);
+
+      /* Set actual timestep correctly. */
+      *dtThisLevel = min(*dtThisLevel, dt_conduction);
+
     }
 
+    dtActual = *dtThisLevel;
+
 #ifdef USE_DT_LIMIT
 
     //    dtLimit = LevelZeroDeltaT/(4.0)/POW(RefineBy,level);

File src/enzo/global_data.h

View file
 
 EXTERN int LevelCycleCount[MAX_DEPTH_OF_HIERARCHY];
 EXTERN int LevelSubCycleCount[MAX_DEPTH_OF_HIERARCHY];
+EXTERN float dtRebuildHierarchy[MAX_DEPTH_OF_HIERARCHY];
+EXTERN float TimeSinceRebuildHierarchy[MAX_DEPTH_OF_HIERARCHY];
 EXTERN float dtThisLevelSoFar[MAX_DEPTH_OF_HIERARCHY];
 EXTERN float dtThisLevel[MAX_DEPTH_OF_HIERARCHY];