Split computation of Dedner Wave speeds into a separate file.

--HG--
branch : week-of-code

src/enzo/Make.config.objects

@@ -701,6 +701,7 @@ OBJS_HYDRO_RK = \
       	hydro_rk/AGNDiskInitialize.o \
         hydro_rk/Collapse1DInitialize.o \
         hydro_rk/Collapse3DInitialize.o \
+        hydro_rk/ComputeDednerWaveSpeeds.o \
         hydro_rk/EvolveLevel_RK2.o \
         hydro_rk/GalaxyDiskInitialize.o \
         hydro_rk/Grid_AddSelfGravity.o \

src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C

@@ -0,0 +1,77 @@
+/***********************************************************************
+/
+/  COMPUTE DEDNER WAVE SPEEDS
+/
+/  written by: Tom Abel
+/  date:       October 2009
+/
+/  ======================================================================= 
+/  This routine computes the wave speeds used for the Dedner formalism.
+/
+/       Reference: e.g.  Matsumoto, PASJ, 2007, 59, 905 
+/
+/  Output: C_h and C_p global variables defined in global_data.h
+/
+************************************************************************/ 
+
+#include <stdio.h>
+#include <math.h>
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "TopGridData.h"
+#include "Grid.h"
+#include "Hierarchy.h"
+#include "LevelHierarchy.h"
+#include "../hydro_rk/tools.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+
+int ComputeDednerWaveSpeeds(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], 
+			    int level, FLOAT dt0)
+{
+  /* Count the number of grids on this level. */
+
+  if (HydroMethod != MHD_RK) 
+    return SUCCESS;
+
+  float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, 
+    VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits;
+  double MassUnits;
+  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, 1.0);
+  MassUnits = DensityUnits*pow(LengthUnits,3);
+
+  int lmax;
+  LevelHierarchyEntry *Temp;
+  for (lmax = MAX_DEPTH_OF_HIERARCHY-1; lmax >= 0; lmax--) {
+    Temp = LevelArray[lmax];
+    if (Temp != NULL) 
+      break;
+  }
+
+  //      lmax = 0; // <- Pengs version had lmax = 6
+  FLOAT dx0, dy0, dz0, h_min, DivBDampingLength = 1.0;
+
+  dx0 = (DomainRightEdge[0] - DomainLeftEdge[0]) / MetaData->TopGridDims[0];
+  dy0 = (MetaData->TopGridRank > 1) ? 
+    (DomainRightEdge[1] - DomainLeftEdge[1]) / MetaData->TopGridDims[1] : 1e8;
+  dz0 = (MetaData->TopGridRank > 2) ? 
+    (DomainRightEdge[2] - DomainLeftEdge[2]) / MetaData->TopGridDims[2] : 1e8;
+  h_min = my_MIN(dx0, dy0, dz0);
+  h_min /= pow(RefineBy, lmax);
+  C_h = 0.5*MetaData->CourantSafetyNumber*(h_min/dt0);
+  C_h = min( C_h, 1e6/VelocityUnits); // never faster than __ cm/s (for very small dt0 a problems)
+  if (EOSType == 3)  // for isothermal runs just use the constant sound speed
+    C_h = EOSSoundSpeed;
+
+  C_p = sqrt(0.18*DivBDampingLength*C_h);
+
+  return SUCCESS;
+}
+
+

src/enzo/hydro_rk/EvolveLevel_RK2.C

@@ -42,7 +42,8 @@ int StarParticleFinalize(HierarchyEntry *Grids[], TopGridData *MetaData,
 			 int level, Star *&AllStars);
 int AdjustRefineRegion(LevelHierarchyEntry *LevelArray[], 
 		       TopGridData *MetaData, int EL_level);
-
+int ComputeDednerWaveSpeeds(TopGridData *MetaData,LevelHierarchyEntry *LevelArray[], 
+			    int level, FLOAT dt0);
 #ifdef TRANSFER
 int EvolvePhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 		  Star *AllStars, FLOAT GridTime, int level, int LoopTime = TRUE);
@@ -161,9 +162,6 @@ int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh);
 int FastSiblingLocatorInitializeStaticChainingMesh(ChainingMeshStructure *Mesh, int Rank,
 						   int TopGridDims[]); 
 
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
 double ReturnWallTime();
 int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
                int level);
@@ -223,12 +221,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
   FLOAT When;
 
-  float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, 
-    VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits;
-  double MassUnits;
-  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	   &TimeUnits, &VelocityUnits, 1.0);
-  MassUnits = DensityUnits*pow(LengthUnits,3);
 
   /* Create an array (Grids) of all the grids. */
 
@@ -366,40 +358,15 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
     } // end loop over grids (create Subgrid list)
 
-    /* compute wave speed
-       Reference: Matsumoto, PASJ, 2007, 59, 905 */
-
-    if (HydroMethod == MHD_RK) {
-      int lmax;
-      LevelHierarchyEntry *Temp;
-      for (lmax = MAX_DEPTH_OF_HIERARCHY-1; lmax >= 0; lmax--) {
-	Temp = LevelArray[lmax];
-	if (Temp != NULL) {
-	  break;
-	} 
-      }
-      //      lmax = 0; // <- Pengs version had lmax = 6
-      //      lmax = 1;
-      FLOAT dx0, dy0, dz0, h_min, DivBDampingLength = 1.0;
-
-      dx0 = (DomainRightEdge[0] - DomainLeftEdge[0]) / MetaData->TopGridDims[0];
-      dy0 = (MetaData->TopGridRank > 1) ? 
-	(DomainRightEdge[1] - DomainLeftEdge[1]) / MetaData->TopGridDims[1] : 1e8;
-      dz0 = (MetaData->TopGridRank > 2) ? 
-	(DomainRightEdge[2] - DomainLeftEdge[2]) / MetaData->TopGridDims[2] : 1e8;
-      h_min = my_MIN(dx0, dy0, dz0);
-      h_min /= pow(RefineBy, lmax);
-      C_h = 0.5*MetaData->CourantSafetyNumber*(h_min/dt0);
-      C_h = min( C_h, 1e6/VelocityUnits); // never faster than __ cm/s (for very small dt0 a problems)
-      if (EOSType == 3)  // for isothermal runs just use the constant sound speed
-	C_h = EOSSoundSpeed;
-
-      C_p = sqrt(0.18*DivBDampingLength*C_h);
-      //      C_p = sqrt(0.18*DivBDampingLength)*C_h;
-      if (debug) 
-	fprintf(stderr, "lengthscale %g timestep: %g  C_h: %g  C_p: %g\n ", 
-		h_min, dt0, C_h, C_p);
-    }
+    /* compute wave speed Reference: Matsumoto, PASJ, 2007, 59, 905 */
+
+    if (HydroMethod == MHD_RK) 
+      ComputeDednerWaveSpeeds(MetaData, LevelArray, level, dt0);
+
+    if (debug && HydroMethod == MHD_RK) 
+      fprintf(stderr, "wave speeds: timestep: %g  C_h: %g  C_p: %g\n ", 
+	       dt0, C_h, C_p);
+
 
     When = 0.5;
     RK2SecondStepBaryonDeposit = 0;