1. Matthew Turk
  2. enzo-dev-answer-testing

Commits

dcollins4096  committed 821e70e Merge

merge with 293b81e1404c

  • Participants
  • Parent commits 8a5b935, 293b81e
  • Branches mhdctmerge

Comments (0)

Files changed (9)

File src/enzo/ProtostellarCollapseInitialize.C

View file
   
   /* set up the root grid */
 
-  if (MaximumRefinementLevel > 0 && startWithSubgrids)
+  if (MaximumRefinementLevel > 0 && startWithSubgrids){
     if (Subgrid[0]->GridData->ProjectSolutionToParentGrid(*(TopGrid.GridData))
 	== FAIL) {
             ENZO_FAIL("Error in ProjectSolutionToParentGrid.");
     }
 
-  else
+  }else{
     if (TopGrid.GridData->ProtostellarCollapseInitializeGrid(
                           ProtostellarCollapseCoreDensity,
 			  ProtostellarCollapseCoreEnergy,
 			  ProtostellarCollapseAngularVelocity) == FAIL) {
             ENZO_FAIL("Error in ProtostellarCollapseInitializeGrid.");
     }
+  }
 
   /* set up field names and units */
 

File src/enzo/kludges/EvolveHierarchy.C

-/***********************************************************************
-/
-/  EVOLVE HIERARCHY FUNCTION
-/
-/  written by: Greg Bryan
-/  date:       November, 1994
-/  modified1:  February, 1995 by GB
-/              Changed to reflect changes in EvolveGrid & EvolveTopGrid.
-/  modified2:  July, 1995 by GB
-/              Changed to reflect new routine EvolveLevel.
-/  modified3:  February, 2006 by Daniel Reynolds
-/              Updated call interface to ComputePotentialFieldLevelZero
-/  modified4:  February, 2007 by Robert Harkness
-/              Group and in-core i/o
-/  modified5:  December, 2007 by Robert Harkness
-/              Remove calls to ComputePotential to allow use of the
-/              Fast Sibling Locator to remedy Ncpu^2 scaling
-/  modified6:  February, 2008 by Robert Harkness
-/              Conditional calls to MPI_Barrier to force MPICH progress
-/  modified7:  October, 2009 by Ji-hoon Kim
-/              Added particle splitter routine
-/
-/  PURPOSE:
-/    This routine is responsible for the evolution of the grid hierarchy.
-/    It assumes the hierarchy is already constructed and the grids
-/    initialized.  The routine then loops over time until one of the
-/    stopping criteria is reached.  This routine also handles data dumps,
-/    history dumps and restarts dumps (although the later two have not
-/    yet been implemented).
-/
-************************************************************************/
-#include "preincludes.h"
- 
-#ifdef USE_MPI
-#include <mpi.h>
-#endif
- 
-#include <stdio.h>
- 
-#include "performance.h"
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-#include "LevelHierarchy.h"
-#include "CosmologyParameters.h"
-#include "communication.h"
-#include "CommunicationUtilities.h"
-#ifdef MHDCT
-#include "DaveTools.h"
-#endif //MHDCT 
-#ifdef TRANSFER
-#include "ImplicitProblemABC.h"
-#endif
- 
-// function prototypes
- 
-int RebuildHierarchy(TopGridData *MetaData,
-		     LevelHierarchyEntry *LevelArray[], int level);
-
-int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
-		int level, float dtLevelAbove, ExternalBoundary *Exterior
-#ifdef TRANSFER
-		, ImplicitProblemABC *ImplicitSolver
-#endif
-		);
-
-int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
-                    int level, float dtLevelAbove, ExternalBoundary *Exterior, 
-#ifdef TRANSFER
-		    ImplicitProblemABC *ImplicitSolver, 
-#endif
-		    FLOAT dt0);
-
-int WriteAllData(char *basename, int filenumber,
-		 HierarchyEntry *TopGrid, TopGridData &MetaData,
-		 ExternalBoundary *Exterior, 
-#ifdef TRANSFER
-		 ImplicitProblemABC *ImplicitSolver,
-#endif		 
-		 FLOAT WriteTime = -1);
-
-int Group_WriteAllData(char *basename, int filenumber,
-		 HierarchyEntry *TopGrid, TopGridData &MetaData,
-		 ExternalBoundary *Exterior, 
-#ifdef TRANSFER
-		 ImplicitProblemABC *ImplicitSolver,
-#endif		 
-		 FLOAT WriteTime = -1,
-		 int RestartDump = FALSE);
-
-int CopyOverlappingZones(grid* CurrentGrid, TopGridData *MetaData,
-			 LevelHierarchyEntry *LevelArray[], int level);
-int TestGravityCheckResults(LevelHierarchyEntry *LevelArray[]);
-int TestGravitySphereCheckResults(LevelHierarchyEntry *LevelArray[]);
-int CheckForOutput(HierarchyEntry *TopGrid, TopGridData &MetaData,
-		   ExternalBoundary *Exterior, 
-#ifdef TRANSFER
-		   ImplicitProblemABC *ImplicitSolver,
-#endif		 
-		   int Restart = FALSE);
-int CheckForTimeAction(LevelHierarchyEntry *LevelArray[],
-		       TopGridData &MetaData);
-int CheckForResubmit(TopGridData &MetaData, int &Stop);
-int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
-int OutputLevelInformation(FILE *fptr, TopGridData &MetaData,
-			   LevelHierarchyEntry *LevelArray[]);
-int PrepareGravitatingMassField(HierarchyEntry *Grid, TopGridData *MetaData,
-				LevelHierarchyEntry *LevelArray[], int level);
-int ReduceFragmentation(HierarchyEntry &TopGrid, TopGridData &MetaData,
-			ExternalBoundary *Exterior,
-			LevelHierarchyEntry *LevelArray[]);
-int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL,
-				int NumberOfSubgrids[] = NULL,
-				int FluxFlag = FALSE,
-				TopGridData* MetaData = NULL);
-double ReturnWallTime(void);
-int Enzo_Dims_create(int nnodes, int ndims, int *dims);
-int FOF(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], 
-	int WroteData, int FOFOnly=FALSE);
-int StarParticleCountOnly(LevelHierarchyEntry *LevelArray[]);
-int CommunicationLoadBalanceRootGrids(LevelHierarchyEntry *LevelArray[], 
-				      int TopGridRank, int CycleNumber);
-int CommunicationBroadcastValue(Eint32 *Value, int BroadcastProcessor);
-int CommunicationBroadcastValue(Eint64 *Value, int BroadcastProcessor);
-int ParticleSplitter(LevelHierarchyEntry *LevelArray[], int ThisLevel,
-		     TopGridData *MetaData); 
-int MagneticFieldResetter(LevelHierarchyEntry *LevelArray[], int ThisLevel,
-			  TopGridData *MetaData); 
-void PrintMemoryUsage(char *str);
-int SetEvolveRefineRegion(FLOAT time);
-
-#ifdef MEM_TRACE
-Eint64 mused(void);
-#endif
-#ifdef USE_PYTHON
-int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
-               int level, int from_topgrid);
-#endif
-
- 
- 
-#define NO_REDUCE_FRAGMENTATION
- 
-#ifdef MHDCT 
-int ExtraOutput(int output_flag, LevelHierarchyEntry *LevelArray[],TopGridData *MetaData, int level, ExternalBoundary *Exterior
-#ifdef TRANSFER
-			  , ImplicitProblemABC *ImplicitSolver
-#endif
-        );
-#endif //MHDCT
- 
- 
-int EvolveHierarchy(HierarchyEntry &TopGrid, TopGridData &MetaData,
-                    ExternalBoundary *Exterior,
-#ifdef TRANSFER
-		    ImplicitProblemABC *ImplicitSolver,
-#endif
-		    LevelHierarchyEntry *LevelArray[], float Initialdt)
-{
- 
-#ifdef MHDCT
-  fprintf(stderr,"PRGIO FORCED ON!!!\n");
-  ParallelRootGridIO = TRUE;
-#endif //MHDCT
-  float dt;
- 
-  int i, dim, Stop = FALSE;
-  int Restart = FALSE;
-  double tlev0, tlev1, treb0, treb1, tloop0, tloop1, tentry, texit;
-  LevelHierarchyEntry *Temp;
-  double LastCPUTime;
-
-  LCAPERF_BEGIN("EL");
-  LCAPERF_START("EvolveHierarchy");
-
-#ifdef USE_MPI
-  tentry = MPI_Wtime();
-#endif
- 
-  if (MetaData.Time        >= MetaData.StopTime ) Stop = TRUE;
-  if (MetaData.CycleNumber >= MetaData.StopCycle) Stop = TRUE;
-  MetaData.StartCPUTime = LastCPUTime = ReturnWallTime();
-  MetaData.LastCycleCPUTime = 0.0;
-
-  // Reset CPUTime, if it's very large (absolute UNIX time), which
-  // was the default from before.
-  if (MetaData.CPUTime > 1e2*MetaData.StopCPUTime)
-    MetaData.CPUTime = 0.0;
- 
-  /* Attach RandomForcingFields to BaryonFields temporarily to apply BCs */
- 
-  if (RandomForcing) { //AK
-    Temp = LevelArray[0];
-    while (Temp != NULL) {
-      Temp->GridData->AppendForcingToBaryonFields();
-      Temp = Temp->NextGridThisLevel;
-    }
-    Exterior->AppendForcingToBaryonFields();
-  }
- 
-  /* Set top grid boundary conditions. */
-
-  Temp = LevelArray[0];
-
-  PrintMemoryUsage("Enter EH");
-
-#ifdef FORCE_MSG_PROGRESS
-  CommunicationBarrier();
-#endif
-
-  CommunicationReceiveIndex = 0;
-  CommunicationDirection = COMMUNICATION_POST_RECEIVE;
-  CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE;
-
-  while (Temp != NULL) {
-    if (Temp->GridData->SetExternalBoundaryValues(Exterior) == FAIL) {
-      //      ENZO_FAIL("Error in grid->SetExternalBoundaryValues.\n");
-      Exterior->Prepare(Temp->GridData);
-
-    }
-    if (CopyOverlappingZones(Temp->GridData, &MetaData, LevelArray, 0)
-	== FAIL)
-      ENZO_FAIL("Error in CopyOverlappingZones.");
-    Temp = Temp->NextGridThisLevel;
-  }
-  
-  CommunicationDirection = COMMUNICATION_SEND;
-
-  Temp = LevelArray[0];
-  while (Temp != NULL) {
-    if (CopyOverlappingZones(Temp->GridData, &MetaData, LevelArray, 0)
-	== FAIL)
-      ENZO_FAIL("Error in CopyOverlappingZones.");
-    Temp = Temp->NextGridThisLevel;
-  }
-
-#ifdef FORCE_MSG_PROGRESS 
-  CommunicationBarrier();
-#endif
-
-  CommunicationReceiveHandler();
-
-#ifdef FORCE_MSG_PROGRESS
-  CommunicationBarrier();
-#endif
-
-  PrintMemoryUsage("Bdry set");
- 
-  /* Remove RandomForcingFields from BaryonFields when BCs are set. */
- 
-  if (RandomForcing) { //AK
-    LevelHierarchyEntry *Temp = LevelArray[0];
-    while (Temp != NULL) {
-      Temp->GridData->DetachForcingFromBaryonFields();
-      Temp = Temp->NextGridThisLevel;
-    }
-    Exterior->DetachForcingFromBaryonFields();
-  }
- 
-  /* Check for output. */
- 
-  CheckForOutput(&TopGrid, MetaData, Exterior, 
-#ifdef TRANSFER
-		 ImplicitSolver,
-#endif		 
-		 Restart);
-
-  PrintMemoryUsage("Output");
- 
-  /* Compute the acceleration field so ComputeTimeStep can find dtAccel.
-     (Actually, this is a huge pain-in-the-ass, so only do it if the
-      problem really requires it). */
- 
-/*
-  if (ProblemType == 21) {
-    PrepareGravitatingMassField(&TopGrid, &MetaData, LevelArray, 0);
-    ComputePotentialFieldLevelZero(&MetaData, Grids, NumberOfGrids);
-    TopGrid.GridData->ComputeAccelerationField(GRIDS);
-  }
-*/
-
-
-  /* Do the first grid regeneration. */
- 
-  if(CheckpointRestart == FALSE) {
-    RebuildHierarchy(&MetaData, LevelArray, 0);
-  }
-#ifdef MHDCT
-  ExtraOutput(0,LevelArray,&MetaData,0,Exterior
-#ifdef TRANSFER
-		      , ImplicitSolver
-#endif
-          );
-#endif //MHDCT
-
-  PrintMemoryUsage("1st rebuild");
- 
-  /* Particle Splitter. Split particles into 13 (=1+12) child particles */
-  
-  if (MetaData.FirstTimestepAfterRestart == TRUE &&
-      ParticleSplitterIterations > 0)
-    ParticleSplitter(LevelArray, 0, &MetaData);
-
-  /* Reset magnetic fields if requested. */
-  
-  if (MetaData.FirstTimestepAfterRestart == TRUE &&
-      ResetMagneticField == TRUE)
-    MagneticFieldResetter(LevelArray, 0, &MetaData);
-
-  /* Open the OutputLevelInformation file. */
- 
-  FILE *LevelInfofptr;
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    LevelInfofptr = fopen("OutputLevelInformation.out", "w");
-    fclose(LevelInfofptr);
-  }
-
-  /* For top-level timestepping with radiative star particles, we want
-     to restrict the timesteps.  Collect info here. */
-
-  StarParticleCountOnly(LevelArray);
- 
-#ifdef USE_LCAPERF
-  Eint32 lcaperf_iter;
-#endif
-
-  LCAPERF_STOP("EvolveHierarchy");
-  LCAPERF_END("EH");
-
-  /* ====== MAIN LOOP ===== */
-
-  bool FirstLoop = true;
-  while (!Stop) {
-
-#ifdef USE_LCAPERF
-    lcaperf_iter = MetaData.CycleNumber;
-    static bool isFirstCall = true;
-    if ((lcaperf_iter % LCAPERF_DUMP_FREQUENCY)==0 || isFirstCall) lcaperf.begin("EL");
-    isFirstCall = false;
-    lcaperf.attribute ("timestep",&lcaperf_iter, LCAPERF_INT);
-    lcaperf.start("EL");
-#endif
-
-#ifdef USE_MPI
-    tloop0 = MPI_Wtime();
-#endif
-
-#ifdef MEM_TRACE
-    fprintf(memtracePtr, "==== CYCLE %"ISYM" ====\n", MetaData.CycleNumber);
-#endif    
-    PrintMemoryUsage("Top");
-
-    /* Load balance the root grids if this isn't the initial call */
-
-    if ((CheckpointRestart == FALSE) && (!FirstLoop))
-      CommunicationLoadBalanceRootGrids(LevelArray, MetaData.TopGridRank, 
-					MetaData.CycleNumber);
-
-    /* Output level information to log file. */
- 
-    if (MyProcessorNumber == ROOT_PROCESSOR) {
-      LevelInfofptr = fopen("OutputLevelInformation.out", "a");
-      if (LevelInfofptr == NULL)
-        ENZO_FAIL("Can't open OutputLevelInformation.out!");
-    }
-
-    // OutputLevelInformation() only needs to be called by all processors
-    // when lcaperf is enabled.
-
-    OutputLevelInformation(LevelInfofptr, MetaData, LevelArray);
-
-    if (MyProcessorNumber == ROOT_PROCESSOR)
-      fclose(LevelInfofptr);
- 
-    /* Compute minimum timestep on the top level. */
- 
-    float dtProc   = huge_number;
-    Temp = LevelArray[0];
- 
-    // Start skipping
-    if(CheckpointRestart == FALSE) {
-      while (Temp != NULL) {
-        dtProc = min(dtProc, Temp->GridData->ComputeTimeStep());
-        Temp = Temp->NextGridThisLevel;
-      }
-
-      dt = RootGridCourantSafetyNumber*CommunicationMinValue(dtProc);
-      dt = min(MetaData.MaximumTopGridTimeStep, dt);
-
-      if (debug) fprintf(stderr, "dt, Initialdt: %g %g \n", dt, Initialdt);
-      if (Initialdt != 0) {
-      
-	dt = min(dt, Initialdt);
-	if (debug) fprintf(stderr, "dt, Initialdt: %g %g \n", dt, Initialdt);
-        Initialdt = 0;
-      }
-
-      /* Make sure timestep doesn't go past an output. */
-
-      if (ComovingCoordinates)
-        for (i = 0; i < MAX_NUMBER_OF_OUTPUT_REDSHIFTS; i++)
-          if (CosmologyOutputRedshift[i] != -1)
-            dt = min(1.0001*(CosmologyOutputRedshiftTime[i]-MetaData.Time), dt);
-      for (i = 0; i < MAX_TIME_ACTIONS; i++)
-        if (TimeActionTime[i] > 0 && TimeActionType[i] > 0)
-          dt = min(1.0001*(TimeActionTime[i] - MetaData.Time), dt);
-      if (MetaData.dtDataDump > 0.0) {
-        while (MetaData.TimeLastDataDump+MetaData.dtDataDump < MetaData.Time)
-          MetaData.TimeLastDataDump += MetaData.dtDataDump;
-        dt = min(1.0001*(MetaData.TimeLastDataDump + MetaData.dtDataDump -
-              MetaData.Time), dt);
-      }
-
-      /* Set the time step.  If it will cause Time += dt > StopTime, then
-         set dt = StopTime - Time */
-
-      dt = min(MetaData.StopTime - MetaData.Time, dt);
-    } else { 
-      dt = dtThisLevel[0]; 
-    }
-
-    /* Set the time step.  If it will cause Time += dt > StopTime, then
-       set dt = StopTime - Time */
- 
-    dt = min(MetaData.StopTime - MetaData.Time, dt);
-    Temp = LevelArray[0];
-    // Stop skipping
-
-    // Set dt from stored in CheckpointRestart
-    while (Temp != NULL) {
-      Temp->GridData->SetTimeStep(dt);
-      Temp = Temp->NextGridThisLevel;
-    }
- 
-#ifdef CONFIG_TESTING
-    if (MyProcessorNumber == ROOT_PROCESSOR) {
-      printf("enzo-test: MetaData.CycleNumber %"ISYM"\n", MetaData.CycleNumber);
-      printf("enzo-test: dt %.14g\n",dt);
-      printf("enzo-test: MetaData.Time %"GOUTSYM"\n", MetaData.Time);
-      fflush(stdout);
-    }
-#endif
-
-    if (MyProcessorNumber == ROOT_PROCESSOR) {
-      fprintf(stderr, "TopGrid dt = %"ESYM"     time = %"GOUTSYM"    cycle = %"ISYM,
-	     dt, MetaData.Time, MetaData.CycleNumber);
-
-      if (ComovingCoordinates) {
-	FLOAT a, dadt;
-	CosmologyComputeExpansionFactor(MetaData.Time, &a, &dadt);
-	fprintf(stderr, "    z = %"GOUTSYM, (1 + InitialRedshift)/a - 1);
-      }
-      fprintf(stderr, "\n");
-    }
-    //}
- 
-    /* Inline halo finder */
-
-    FOF(&MetaData, LevelArray, MetaData.WroteData);
-
-    /* If provided, set RefineRegion from evolving RefineRegion */
-    if ((RefineRegionTimeType == 1) || (RefineRegionTimeType == 0)) {
-        if (SetEvolveRefineRegion(MetaData.Time) == FAIL) {
-          fprintf(stderr, "Error in SetEvolveRefineRegion.\n");
-          return FAIL;
-        }
-    }
-    /* Evolve the top grid (and hence the entire hierarchy). */
-
-#ifdef USE_MPI 
-    CommunicationBarrier();
-    tlev0 = MPI_Wtime();
-#endif
-
-    /* Zeroing out the rootgrid Emissivity before EvolveLevel is called 
-       so when rootgrid emissivity values are calculated they are put in 
-       clean rootgrid array */
-#ifdef EMISSIVITY
-/*
-    if(StarMakerEmissivityField > 0){
-      LevelHierarchyEntry *RootTemp;
-      RootTemp = LevelArray[0];
-      while (RootTemp != NULL) {
-	RootTemp->GridData->ClearEmissivity();
-	RootTemp = RootTemp->NextGridThisLevel;
-      }
-    }
-*/
-#endif
-
-    if (HydroMethod == PPM_DirectEuler || HydroMethod == Zeus_Hydro || 
-	HydroMethod == PPM_LagrangeRemap || HydroMethod == HydroMethodUndefined ||
-#ifdef MHDCT
-	HydroMethod == MHD_Li ||
-#endif //MHDCT
-	HydroMethod < 0) {
-      if (EvolveLevel(&MetaData, LevelArray, 0, dt, Exterior
-#ifdef TRANSFER
-		      , ImplicitSolver
-#endif
-		      ) == FAIL) {
-        if (NumberOfProcessors == 1) {
-          fprintf(stderr, "Error in EvolveLevel.\n");
-          fprintf(stderr, "--> Dumping data (output number %d).\n",
-                  MetaData.DataDumpNumber);
-	Group_WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber,
-		     &TopGrid, MetaData, Exterior
-#ifdef TRANSFER
-		     , ImplicitSolver
-#endif		 
-		     );
-        }
-        return FAIL;
-      }
-    } else {
-      if (HydroMethod == HD_RK || HydroMethod == MHD_RK)
-	if (EvolveLevel_RK2(&MetaData, LevelArray, 0, dt, Exterior, 
-#ifdef TRANSFER
-			    ImplicitSolver, 
-#endif
-			    dt) == FAIL) {
-	  if (NumberOfProcessors == 1) {
-	    fprintf(stderr, "Error in EvolveLevel_RK2.\n");
-	    fprintf(stderr, "--> Dumping data (output number %d).\n",
-		    MetaData.DataDumpNumber);
-	    Group_WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber,
-			       &TopGrid, MetaData, Exterior
-#ifdef TRANSFER
-			       , ImplicitSolver
-#endif		 
-			       );
-	  }
-        return FAIL;
-      }
-    }
-
-
-
-#ifdef USE_MPI 
-    CommunicationBarrier();
-    tlev1 = MPI_Wtime();
-#endif
-  
-    /* Add time and check stopping criteria (steps #21 & #22)
-       (note the topgrid is also keeping its own time but this statement will
-       keep the two in synch). */
- 
-    MetaData.Time += dt;
-    MetaData.CycleNumber++;
-    MetaData.LastCycleCPUTime = ReturnWallTime() - LastCPUTime;
-    MetaData.CPUTime += MetaData.LastCycleCPUTime;
-    LastCPUTime = ReturnWallTime();
-
-    if (MyProcessorNumber == ROOT_PROCESSOR) {
-	
-    if (MetaData.Time >= MetaData.StopTime) {
-      if (MyProcessorNumber == ROOT_PROCESSOR)
-	printf("Stopping on top grid time limit.\n");
-      Stop = TRUE;
-    } else
-    if (MetaData.CycleNumber >= MetaData.StopCycle) {
-      if (MyProcessorNumber == ROOT_PROCESSOR)
-	printf("Stopping on top grid cycle limit.\n");
-      Stop = TRUE;
-    } else
-    if (MetaData.CPUTime >= MetaData.StopCPUTime) {
-      if (MyProcessorNumber == ROOT_PROCESSOR)
-	printf("Stopping on CPU time limit.\n");
-      Stop = TRUE;
-    } else
-    if ((ReturnWallTime() - MetaData.StartCPUTime >= MetaData.dtRestartDump &&
-	 MetaData.dtRestartDump > 0) ||
-	(MetaData.CycleNumber - MetaData.CycleLastRestartDump >= 
-	 MetaData.CycleSkipRestartDump &&
-	 MetaData.CycleSkipRestartDump > 0)) {
-      if (MyProcessorNumber == ROOT_PROCESSOR)
-	printf("Stopping to restart.\n");
-      Stop = TRUE;
-      Restart = TRUE;
-    }
-    } // ENDIF ROOT_PROCESSOR
-
-    CommunicationBroadcastValue(&Stop, ROOT_PROCESSOR);
-    CommunicationBroadcastValue(&Restart, ROOT_PROCESSOR);
-
-    /* If not restarting, rebuild the grids from level 0. */
-
-#ifdef USE_MPI
-    treb0 = MPI_Wtime();
-#endif
-
-    PrintMemoryUsage("Pre loop rebuild");
-#ifdef MHDCT
-  ExtraOutput(21,LevelArray,&MetaData,0,Exterior
-#ifdef TRANSFER
-		      , ImplicitSolver
-#endif
-          );
-#endif //MHDCT
- 
-    if (ProblemType != 25 && Restart == FALSE)
-      RebuildHierarchy(&MetaData, LevelArray, 0);
-#ifdef MHDCT
-  ExtraOutput(22,LevelArray,&MetaData,0,Exterior
-#ifdef TRANSFER
-		      , ImplicitSolver
-#endif
-          );
-#endif //MHDCT
-
-    PrintMemoryUsage("Post loop rebuild");
-
-#ifdef USE_MPI
-    treb1 = MPI_Wtime();
-#endif
- 
-    /* Check for time-actions. */
- 
-    CheckForTimeAction(LevelArray, MetaData);
- 
-    /* Check for output. */
- 
-    CheckForOutput(&TopGrid, MetaData, Exterior, 
-#ifdef TRANSFER
-		   ImplicitSolver,
-#endif		 
-		   Restart);
-
-    /* Call inline analysis. */
-#ifdef USE_PYTHON
-    CallPython(LevelArray, &MetaData, 0, 1);
-#endif
-
-    /* Check for resubmission */
-    
-    if (!Restart)
-      CheckForResubmit(MetaData, Stop);
-
-    /* If stopping, inline halo finder one more time */
-
-    if (Stop && !Restart)
-      FOF(&MetaData, LevelArray, TRUE);
-
-    /* Try to cut down on memory fragmentation. */
- 
-#ifdef REDUCE_FRAGMENTATION
- 
-    if (MetaData.WroteData && !Stop)
-      ReduceFragmentation(TopGrid, MetaData, Exterior, LevelArray);
- 
-#endif /* REDUCE_FRAGMENTATION */
-
-#ifdef USE_LCAPERF
-    lcaperf.stop("EL");
-    if (((lcaperf_iter+1) % LCAPERF_DUMP_FREQUENCY)==0) lcaperf.end("EL");
-#endif
-
-    PrintMemoryUsage("Bot");
-
-  for ( i = 0; i < MAX_NUMBER_OF_TASKS; i++ ) {
-    TaskMemory[i] = -1;
-  }
-
-#ifdef MEM_TRACE
-
-  /*
-  MPI_Datatype DataTypeInt = (sizeof(Eint64) == 4) ? MPI_INT : MPI_LONG_LONG_INT;
-  MPI_Arg ThisTask;
-  MPI_Arg TaskCount;
-  MPI_Arg Count = 1;
-  MPI_Arg stat;
-
-  stat = MPI_Comm_size(MPI_COMM_WORLD, &TaskCount);
-  stat = MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
-  stat = MPI_Allgather(&MemInUse, Count, DataTypeInt, TaskMemory, Count, DataTypeInt, MPI_COMM_WORLD);
-
-  if (ThisTask == 0 ) {
-    for ( i = 0; i < TaskCount; i++) {
-      fprintf(stderr, "TaskMemory : Task %"ISYM"  Memory %"ISYM"\n", i, TaskMemory[i]);
-    }
-  }
-*/
-
-#endif /* MEM_TRACE */
-#ifdef USE_MPI
-    tloop1 = MPI_Wtime();
-#endif
-
-    FILE *evlog;
-
-    if (MyProcessorNumber == ROOT_PROCESSOR) {
-      evlog = fopen("Evtime", "a");
-      fprintf(evlog, "%8"ISYM"  %16.9e  %16.9e  %16.9e\n", MetaData.CycleNumber, tlev1-tlev0, treb1-treb0, tloop1-tloop0);
-      fclose(evlog);
-    }
-
-#ifdef MEM_TRACE
-    Eint64 MemInUse;
-    if (MetaData.WroteData) {
-      MemInUse = mused();
-      MemInUse = CommunicationMaxValue(MemInUse);
-      if (MemInUse > MemoryLimit) {
-        if (MyProcessorNumber == ROOT_PROCESSOR)
-          printf("Stopping due to memory limit.\n");
-        Stop = TRUE;
-      }
-    }
-#endif
-
-    FirstLoop = false;
- 
-  } // ===== end of main loop ====
- 
-#ifdef USE_LCAPERF
-  if (((lcaperf_iter+1) % LCAPERF_DUMP_FREQUENCY)!=0) lcaperf.end("EL");
-  lcaperf.attribute ("timestep",0, LCAPERF_NULL);
-#endif
-
-  MetaData.CPUTime = ReturnWallTime() - MetaData.StartCPUTime;
- 
-  /* Done, so report on current time, etc. */
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    printf("Time     = %9"FSYM"   CycleNumber = %6"ISYM"    Wallclock   = %9"FSYM"\n",
-	   MetaData.Time, MetaData.CycleNumber, MetaData.CPUTime);
-    printf("StopTime = %9"FSYM"   StopCycle   = %6"ISYM"\n",
-	   MetaData.StopTime, MetaData.StopCycle);
-  }
- 
-  /* If we are running problem 23, TestGravity, then check the results. */
- 
-  if (ProblemType == 23)
-    TestGravityCheckResults(LevelArray);
-  if (ProblemType == 25 && NumberOfProcessors == 0)
-    TestGravitySphereCheckResults(LevelArray);
- 
-  /* if we are doing data dumps, then do one last one */
- 
-  if ((MetaData.dtDataDump != 0.0 || MetaData.CycleSkipDataDump != 0) &&
-      !MetaData.WroteData)
-    //#ifdef USE_HDF5_GROUPS
-    if (Group_WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber,
-			   &TopGrid, MetaData, Exterior, 
-#ifdef TRANSFER
-			   ImplicitSolver, 
-#endif		 
-			   -666) == FAIL)
-      ENZO_FAIL("Error in Group_WriteAllData.");
-// #else
-//     if (WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber,
-// 		     &TopGrid, MetaData, Exterior, 
-//#ifdef TRANSFER
-//		     ImplicitSolver, 
-//#endif		 
-//                   -666) == FAIL) {
-//       ENZO_FAIL("Error in WriteAllData.\n");
-//     }
-// #endif
- 
-  /* Write a file to indicate that we're finished. */
-
-  FILE *Exit_fptr;
-  if (!Restart && MyProcessorNumber == ROOT_PROCESSOR) {
-    if ((Exit_fptr = fopen("RunFinished", "w")) == NULL)
-      ENZO_FAIL("Error opening RunFinished.");
-    fprintf(Exit_fptr, "Finished on cycle %"ISYM"\n", MetaData.CycleNumber);
-    fclose(Exit_fptr);
-  }
-
-  if (NumberOfProcessors > 1)
-
-    printf("Communication: processor %"ISYM" CommunicationTime = %"FSYM"\n",
-	   MyProcessorNumber, CommunicationTime);
- 
-  /* done */
-
-#ifdef USE_MPI
-  texit = MPI_Wtime();
-#endif
- 
-  return SUCCESS;
-}

File src/enzo/kludges/FluxFix_Grid_CorrectForRefinedFluxes.C

-#ifdef FLUX_FIX
-/***********************************************************************
-/
-/  GRID CLASS (CORRECT SOLUTION GIVEN ORIGINAL AND REFINED FLUXES)
-/
-/  written by: Greg Bryan
-/  date:       November, 1994
-/  modified1:  Robert Harkness
-/  date:       January, 2003
-/	       Include extra fields beyond Metallicity!
-/  modified2:  David Collins & Rick Wagner
-/  date:       May, 2005
-/	       Include flux correction for outside grids.
-/              Re-instated CorrectBoundaryFlux code.
-/  modified2: David Collins, 2005
-/              Updated algebra so Cosmological Expansion is also
-/              conservative.  This fix also came with fixes to euler.src and
-/              Grid_GetProjectedBoundaryFluxes.C, so make sure you get those.
-/
-/  PURPOSE:    Ensures conservation of stuff.
-/
-/  RETURNS: SUCCESS or FAIL 
-/
-************************************************************************/
- 
-// Given both the original and refined fluxes for a subgrid in the current
-//   grid, correct the solution in the cells just outside the solution to
-//   reflect the new fluxes (which are determined from the solution of
-//   the subgrid).
-//   Note that if the subgrid is on the boundary of the current grid, we
-//     do not correct the values but instead replace the boundary fluxes
-//     for the current time step (BoundaryFluxesThisTimeStep).
-//   Also note that subgrids sharing a face with This Grid, but not proper subgrids,
-//   also need to be taken into account.
- 
-#include <stdio.h>
-#include <math.h>
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
- 
-/* function prototypes */
- 
-int FindField(int f, int farray[], int n);
-int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
- 
-int grid::CorrectForRefinedFluxes(fluxes *InitialFluxes,
-				  fluxes *RefinedFluxes,
-				  fluxes *BoundaryFluxesThisTimeStep,
-				  int SUBlingGrid,
-				  TopGridData *MetaData)
-{
- 
-
-  //fprintf(stderr,"CLOWN Kludge 1\n");
-  //return SUCCESS;
-  // Return if this doesn't concern us.
- 
-  if (ProcessorNumber != MyProcessorNumber || !UseHydro)
-    return SUCCESS;
-
-  if( SUBlingGrid == FALSE){
-    fprintf(stderr,"CFR p%d Only Sublings Today\n",MyProcessorNumber);
-    return SUCCESS;
-  }else{
-    fprintf(stderr,"CFR p%d Doing Sublings\n",MyProcessorNumber);
-  }
-  fprintf(stderr,"CFR: p%d SUBLING %d\n", MyProcessorNumber, SUBlingGrid);
- 
-  // declarations
-  int dtemp = 2;
-    fprintf(stderr,"CFR One p%d dtemp %d IF LeftStart %d RF RightStart %d\n",
-        MyProcessorNumber,
-        dtemp, InitialFluxes->LeftFluxStartGlobalIndex[dtemp][dtemp],
-        RefinedFluxes->RightFluxStartGlobalIndex[dtemp][dtemp]);
- 
-  int i1, i2, i, j, k, dim, field, ffield, index;
-  int FieldIndex, FluxIndex, GridFluxIndex, Offset, RefinedFluxIndex;
-  int End[MAX_DIMENSION], Start[MAX_DIMENSION];
- 
-  //Dimensions of Initial and Refined fluxes.  ("Dim" should be "InitialDim")
-  int Dim[MAX_DIMENSION],RefinedDim[MAX_DIMENSION] = {0,0,0};
- 
-  //Dimension and Start Index for the ParentGrid (this grid) boundary flux.
-  int GridFluxDim[MAX_DIMENSION], GridFluxStartIndex[MAX_DIMENSION];
- 
-  // Used for calculating position in the RefinedFlux and InitialFlux data structure.
-  int RefinedOffset[MAX_DIMENSION] ={0,0,0}, InitialOffset[MAX_DIMENSION] = {0,0,0};
- 
-  //Internal flags: Correct the BaryonField
-  int CorrectLeftBaryonField, CorrectRightBaryonField;
-  //For correction of the Parent Grid (this grid) boundary flux.
-  int CorrectLeftBoundaryFlux, CorrectRightBoundaryFlux;
-  float CorrectionAmountLeft, CorrectionAmountRight; 
- 
-  long_int GlobalDim;
- 
-  /* If there are no fields, don't do anything. */
- 
-  if (NumberOfBaryonFields > 0) {
- 
-    /* Find fields: density, total energy, velocity1-3. */
- 
-    int DensNum, GENum, Vel1Num, Vel2Num, Vel3Num, TENum, B1Num, B2Num, B3Num;
-    if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
-					 Vel3Num, TENum, B1Num, B2Num, B3Num) == FAIL) {
-            ENZO_FAIL("Error in grid->IdentifyPhysicalQuantities.");
-    }
-
-    //dcc kludge:  Just remove a(t)? 09/06/05 
-    /* If using comoving coordinates, compute a(t) because we'll need it
-       to multiply against the CellWidth. */
-
-    //  DC revision September 16th 2005 
-    //    FLOAT a = 1, dadt;
-    //    if (ComovingCoordinates)
-    //      if (CosmologyComputeExpansionFactor(Time, &a, &dadt) == FAIL) {
-    //        ENZO_FAIL("Error in CosmologyComputeExpansionFactors.\n");
-    //      }
- 
- 
-    /* Main loop over all faces. */
-   
-  int dtemp = 2;
-  
-    fprintf(stderr,"CFR Three p%d dtemp %d IF LeftStart %d RF RightStart %d\n",
-        MyProcessorNumber,
-        dtemp, InitialFluxes->LeftFluxStartGlobalIndex[dtemp][dtemp],
-        RefinedFluxes->RightFluxStartGlobalIndex[dtemp][dtemp]);
-    for (dim = 0; dim < GridRank; dim++) {
-      if (GridDimension[dim] > 1) {
-	/* Check that the dims of InitialFluxes & RefinedFluxes are the same */
- 
-	/* don't do this for SUBlings */
-	if( SUBlingGrid == FALSE ){
-	  for (j = 0; j < GridRank; j++)
-	    if ((InitialFluxes->LeftFluxStartGlobalIndex[dim][j] !=
-		 RefinedFluxes->LeftFluxStartGlobalIndex[dim][j])  ||
-		(InitialFluxes->LeftFluxEndGlobalIndex[dim][j] !=
-		 RefinedFluxes->LeftFluxEndGlobalIndex[dim][j])) {
-//	      printf("dim=%d / j=%d //// %d == %d :: %d == %d\n", 
-//		     dim, j, 
-//		     InitialFluxes->LeftFluxStartGlobalIndex[dim][j],
-//		     RefinedFluxes->LeftFluxStartGlobalIndex[dim][j],
-//		     InitialFluxes->LeftFluxEndGlobalIndex[dim][j],
-//		     RefinedFluxes->LeftFluxEndGlobalIndex[dim][j]);
-	      ENZO_FAIL("InitialFluxes & RefinedFluxes are different.\n");
-	    }
-	  /* Error check Fluxes to make sure they all exist. */
-	  for (field = 0; field < NumberOfBaryonFields; field++)
-	    if ((InitialFluxes->LeftFluxes[field][dim] == NULL) ||
-		(RefinedFluxes->LeftFluxes[field][dim] == NULL) ||
-		(InitialFluxes->RightFluxes[field][dim] == NULL) ||
-		(RefinedFluxes->RightFluxes[field][dim] == NULL)) {
-	      fprintf(stderr,"Some Flux data is not present.\n");
-	    return FAIL;
-	  }
-	}
- 
-	//by default, we want to correct the flux.
-	CorrectLeftBaryonField = CorrectRightBaryonField = TRUE;
-	if( SUBlingGrid == TRUE ){
-	
-	  /* calculate Global dimensions on this level */
-	  GlobalDim = 	nlongint(( DomainRightEdge[dim] - DomainLeftEdge[dim])
-				 / CellWidth[dim][0]);
-	
-	  /* get the dims of the refined fluxes to calculate
-	     array indices */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++){
-	    RefinedDim[i] = RefinedFluxes->LeftFluxEndGlobalIndex[dim][i] -
-	      RefinedFluxes->LeftFluxStartGlobalIndex[dim][i] + 1;
-	  }
-	
-	  /* check if SUBling left or right edge lies on
-	     Domain boundary, and if so, modulo the indices,
-	     otherwise, bump the indices to match the initial flux's.
-	  */
-	
-	  if( RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim] == 0 &&
-	      MetaData->LeftFaceBoundaryCondition[dim] == periodic ){
-	    RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim] = GlobalDim - 1;
-	    RefinedFluxes->LeftFluxEndGlobalIndex[dim][dim] = GlobalDim - 1;
-	  }else{
-	    RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim]--;
-	    RefinedFluxes->LeftFluxEndGlobalIndex[dim][dim]--;
-	  }
-	
-	
-	  if( RefinedFluxes->RightFluxStartGlobalIndex[dim][dim] == GlobalDim - 1 &&
-	      MetaData->RightFaceBoundaryCondition[dim] == periodic){
-	    RefinedFluxes->RightFluxStartGlobalIndex[dim][dim] = 0;
-	    RefinedFluxes->RightFluxEndGlobalIndex[dim][dim] = 0;
-	  }else{
-	    RefinedFluxes->RightFluxStartGlobalIndex[dim][dim]++;
-	    RefinedFluxes->RightFluxEndGlobalIndex[dim][dim]++;
-	  }
-	
-	  /* check to see if we're doing this dimension at all.
-	     only the dimension of contact needs to be checked,
-	     since SUBling grids can only have contact along a
-	     single axis. any corrections to this statement
-	     earns a beer */
-  }
-  int dtemp = 2;
-    fprintf(stderr,"CFR Four p%d dtemp %d IF LeftStart %d RF RightStart %d\n",
-        MyProcessorNumber,
-        dtemp, InitialFluxes->LeftFluxStartGlobalIndex[dtemp][dtemp],
-        RefinedFluxes->RightFluxStartGlobalIndex[dtemp][dtemp]);
-	
-  fprintf(stderr,"CFR: p%d IM OK1 %d %d\n",MyProcessorNumber,CorrectLeftBaryonField , CorrectRightBaryonField);
-	
-	  //note these are integers, so comparing them directly is ok.
-	  //Also note that here we don't do the complete check for SUBling-ness.
-	  //More logic is necessary for any two arbitrary grids, but it's been done
-	  //already in populating the SUBling list.  Here, all we need
-	  //to do is determine which face needs the correction, as we already know one exists.
-	
-	  if( InitialFluxes->RightFluxStartGlobalIndex[dim][dim] ==
-	      RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim] ){
-	    CorrectRightBaryonField = TRUE;
-	  }else{
-	    CorrectRightBaryonField = FALSE;
-	  }
-	
-	  if( InitialFluxes->LeftFluxStartGlobalIndex[dim][dim]==
-	      RefinedFluxes->RightFluxStartGlobalIndex[dim][dim] ){
-	    CorrectLeftBaryonField = TRUE;
-	  }else{
-	    CorrectLeftBaryonField = FALSE;
-	  }
-	
-    fprintf(stderr,"CFR IF dim %d LeftStart %d RF RightStart %d\n",
-        dim, InitialFluxes->LeftFluxStartGlobalIndex[dim][dim],
-        RefinedFluxes->RightFluxStartGlobalIndex[dim][dim]);
-  fprintf(stderr,"CFR: IM OK %d %d\n",CorrectLeftBaryonField , CorrectRightBaryonField);
-	  for (i = 0; i < MAX_DIMENSION; i++) {
-	    /* calculate the offset, so the index of the refined fluxes can
-	       be determined from the grid's index */
-	    RefinedOffset[i] = max(InitialFluxes->LeftFluxStartGlobalIndex[dim][i]-
-				   RefinedFluxes->LeftFluxStartGlobalIndex[dim][i],0);
-	  }
- 
-	  RefinedFluxes->LeftFluxStartGlobalIndex[dim][dim]=0;
-	  RefinedOffset[dim]=0;
- 
-	}//Subling == TRUE	
-  int dtemp = 2;
-    fprintf(stderr,"CFR Two p%d dtemp %d IF LeftStart %d RF RightStart %d\n",
-        MyProcessorNumber,
-        dtemp, InitialFluxes->LeftFluxStartGlobalIndex[dtemp][dtemp],
-        RefinedFluxes->RightFluxStartGlobalIndex[dtemp][dtemp]);
- 
-	if( CorrectLeftBaryonField || CorrectRightBaryonField){
-	
-	  /* Compute Start and end indicies of flux region (with respect to
-	     the current grid's flux region). */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++) {
-	    Start[i] = 0;
-	    End[i] = 0;
-	  }
-	
-	  /* start index = subgrid flux left edge global index -
-	     grid far left edge global index
-	     end index = subgrid flux right edge -
-	     grid far left edge global index. */
-	
-	  /* modified to account for different dimensions to the
-	     initial and refined fluxes */
-	
-	  for (i = 0; i < GridRank; i++) {
-	    Start[i] = max(InitialFluxes->LeftFluxStartGlobalIndex[dim][i],
-			   RefinedFluxes->LeftFluxStartGlobalIndex[dim][i]) -
-	      nlongint((CellLeftEdge[i][0] - DomainLeftEdge[i])/
-		       CellWidth[i][0]);
-	    End[i] = min(InitialFluxes->LeftFluxEndGlobalIndex[dim][i],
-			 RefinedFluxes->LeftFluxEndGlobalIndex[dim][i]) -
-	      nlongint((CellLeftEdge[i][0] - DomainLeftEdge[i])/
-		       CellWidth[i][0]);
-	
-	
-	    if (Start[i] < 0 || End[i] > GridDimension[i]) {
-	      fprintf(stderr, "Start/End[%"ISYM"] = %"ISYM"/%"ISYM"\n",
-		      dim, Start[i], End[i]);
-	      fprintf(stderr, "%"GOUTSYM" %"GOUTSYM" %lld\n",
-		      CellLeftEdge[i][0], CellWidth[i][0],
-		      InitialFluxes->LeftFluxStartGlobalIndex[dim][i]);
-	      ENZO_FAIL("Error in FluxFix_Grid_CorrectForRefinedFluxes!\n");
-	    }
-	  }
-	
-	  /* Correct vector to point at cells just outside the left face.
-	     Start[dim] and End[dim] should be the same because the
-	     layer to be corrected is but one cell thick. */
-	
-	  Start[dim] = max(Start[dim] - 1, 0);
-	  End[dim]   = Start[dim];
- 
-	  /* Compute Dimensions of InitialFluxes */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++)
-	    Dim[i] = End[i] - Start[i] + 1;
-	
-	  /* Compute Offset (in baryon field) for right side of face.
-	     The +2 is there because we want to correct the cells just the
-	     right face.*/
-	
-	  Offset = InitialFluxes->RightFluxStartGlobalIndex[dim][dim] -
-	    InitialFluxes->LeftFluxStartGlobalIndex[dim][dim] + 2;
-	  Offset = min(Offset, GridDimension[dim]-1);  // this isn't needed (?)
- 
- 
-	  //For SUBling grids, alter Offset, Start, and End to reflect that we're
-	  //adjusting the INNER edge of the grid if the SUBgrid is outside of it.
- 
-	  if( SUBlingGrid ){
-	    Offset -= 2;
-	    Start[dim]++;
-	    End[dim]++;
-	
-	    //Also correct Dim, the size of the Initial Flux: it comes from This Grid,
-	    //not the Subgrid.
-	
-	    for(i=0;i<GridRank;i++)
-	      if(i != dim){
-		Dim[i] = GridEndIndex[i]-GridStartIndex[i]+1;
-		InitialOffset[i] = max( RefinedFluxes->LeftFluxStartGlobalIndex[dim][i]-
-					InitialFluxes->LeftFluxStartGlobalIndex[dim][i],
-					0);
-	      }else{
-		Dim[i] = 1;
-		InitialOffset[i] = 0;
-	      }
-	  }
-	
-	
-	  /* Check to see if we should correct BoundaryFluxesThisTimeStep
-	     instead of the fields themselves. */
-	
-	  CorrectLeftBoundaryFlux = FALSE;
-	  CorrectRightBoundaryFlux = FALSE;
-	
-	  if (Start[dim] == GridStartIndex[dim]-1){
-	    CorrectLeftBoundaryFlux = TRUE;
-	    CorrectLeftBaryonField  = FALSE;
-	  }
-	  if (Start[dim] + Offset == GridEndIndex[dim]+1){
-	    CorrectRightBoundaryFlux = TRUE;
-	    CorrectRightBaryonField  = FALSE;
-	  }
- 
-	  /* Set GridFluxStartIndex to the starting position of the flux
-	     plane (i.e. exclude grid boundary zones), except for the direction
-	     of the flux is set such that GridStartIndex[dim] - Start[dim] = 0 */
-	
-	  for (i = 0; i < MAX_DIMENSION; i++) {
-	    GridFluxStartIndex[i] = GridStartIndex[i];
-	    GridFluxDim[i] = GridEndIndex[i] - GridStartIndex[i] + 1;
-	  }
-	
-	  GridFluxStartIndex[dim] = Start[dim];
-	  GridFluxDim[dim] = 1;
-	
-	  /* Turn Offset (in dim direction) to Offset (in field array) */
-	
-	  for (i = 0; i < dim; i++)
-	    Offset *= GridDimension[i];
-	
-	  /* Multiply faces by density to get conserved quantities
-	     (only multiply fields which we are going to correct) */
-	
-	  if (HydroMethod != Zeus_Hydro) {
-
-	    for (field = 0; field < NumberOfBaryonFields; field++) 
-	      if (FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
-		  FieldTypeIsDensity(FieldType[field]) == FALSE &&
-		  FieldTypeIsRadiation(FieldType[field]) == FALSE &&
-		  FieldType[field] != Bfield1 &&
-		  FieldType[field] != Bfield2 && FieldType[field] != Bfield3 &&
-		  FieldType[field] != PhiField && 
-		  FieldType[field] != DrivingField1 &&
-		  FieldType[field] != DrivingField2 &&
-		  FieldType[field] != DrivingField3 &&
-		  FieldType[field] != GravPotential &&
-		(RadiativeCooling == 0 || (FieldType[field] != TotalEnergy && 
-					   FieldType[field] != InternalEnergy))) {
-		for (k = Start[2]; k <= End[2]; k++) {
-		  for (j = Start[1]; j <= End[1]; j++) {
-		    index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		    for (i = Start[0]; i <= End[0]; i++, index++) {
-		      BaryonField[field][index] *= BaryonField[DensNum][index];
-		      BaryonField[field][index+Offset] *= 
-			BaryonField[DensNum][index+Offset];
-		    }
-		  }
-		}
-	      }
-	  }
-
-	  
-	
-	  /* Divide species by densities so that at the end we can multiply
-	     them by the new density (species are not otherwise modified --
-	     see the next comment).  This ensures that the species are changed
-	     to keep the same fractional density. */
-	
-	  for (field = 0; field < NumberOfBaryonFields; field++)
-	    if (FieldType[field] >= ElectronDensity &&
-		FieldType[field] < Metallicity &&
-		FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
-		FieldTypeIsRadiation(FieldType[field]) == FALSE)
-	      for (k = Start[2]; k <= End[2]; k++)
-		for (j = Start[1]; j <= End[1]; j++) {
-		  index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		  for (i = Start[0]; i <= End[0]; i++, index++) {
-		    BaryonField[field][index] /= BaryonField[DensNum][index];
-		    BaryonField[field][index+Offset] /=
-		      BaryonField[DensNum][index+Offset];
-		  }
-		}
-	
-	  /* Correct face for difference between refined and initial fluxes.
-	     (Don't do this for energy if radiative cooling is on because it's
-	     no longer conserved.  Similarly, don't do it for the species
-	     because they are not individually conserved either -- in fact,
-	     this could be done for the conserved quantities like charge,
-	     total number density summed over ionization, etc.) */
-	
-	  if (Coordinate == Cartesian) {
-	  for (field = 0; field < NumberOfBaryonFields; field++){
-	    if ((FieldTypeNoInterpolate(FieldType[field]) == FALSE) &&
-		(RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-					   FieldType[field] != InternalEnergy)) &&
-		(FieldType[field] < ElectronDensity) && 
-		FieldType[field] != DrivingField1 &&
-		FieldType[field] != DrivingField2 &&
-		FieldType[field] != DrivingField3 &&
-		FieldType[field] != GravPotential) {
-	      for (k = Start[2]; k <= End[2]; k++){
-		for (j = Start[1]; j <= End[1]; j++){
-		  for (i = Start[0]; i <= End[0]; i++) {
-		
-		    /* Compute indexes. */
-		
-		    FieldIndex = (k*GridDimension[1] + j)*GridDimension[0] + i;
-		    FluxIndex  = ((k - Start[2]+InitialOffset[2])*Dim[1] +
-				  (j - Start[1]+InitialOffset[1]))*Dim[0] +
-		      (i - Start[0]+InitialOffset[0]);
-		
-		    if( SUBlingGrid ){
-		      RefinedFluxIndex = ((k - Start[2] + RefinedOffset[2])*RefinedDim[1] +
-					  (j - Start[1] + RefinedOffset[1]))*RefinedDim[0] +
-			(i - Start[0] + RefinedOffset[0]);
-		
-		    }else{
-		      RefinedFluxIndex = FluxIndex;
-		    }
-		
-		    GridFluxIndex =
-		      (i - GridFluxStartIndex[0])
-		      + (j - GridFluxStartIndex[1])*GridFluxDim[0]
-		      + (k - GridFluxStartIndex[2])*GridFluxDim[1]*GridFluxDim[0];
-		
-		
-// 		    if (CorrectLeftBoundaryFlux)
-//		      BoundaryFluxesThisTimeStep->LeftFluxes[field][dim][GridFluxIndex] =
-//			RefinedFluxes->LeftFluxes[field][dim][FluxIndex];
-		
-		    if(CorrectLeftBaryonField){
-		
-		      if( SUBlingGrid == FALSE ){
-			CorrectionAmountLeft = 
-			  InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-			  RefinedFluxes->LeftFluxes[field][dim][FluxIndex];
-			BaryonField[field][FieldIndex] += CorrectionAmountLeft;
-			
-		      }else{ /* if( SUBlingGrid == False) */
-			
-            //fprintf(stderr,"CFR: Doin It Left\n");
-			CorrectionAmountLeft = 
-			  -(InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-			    RefinedFluxes->RightFluxes[field][dim][RefinedFluxIndex]);
-			BaryonField[field][FieldIndex] += CorrectionAmountLeft;
-			BaryonField[field][FieldIndex] = 500;
-			
-		      }
-		    }
-		
-// 		    if (CorrectRightBoundaryFlux)
-//		      BoundaryFluxesThisTimeStep->RightFluxes[field][dim] [GridFluxIndex] =
-//			RefinedFluxes->RightFluxes[field][dim][FluxIndex];
-		
-		    /* update only if necessary */
-		    if(CorrectRightBaryonField){
-		
-		      if( SUBlingGrid == FALSE ){
-
-			CorrectionAmountRight = 
-			  -(InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-			    RefinedFluxes->RightFluxes[field][dim][FluxIndex]);
-			  
-			BaryonField[field][FieldIndex + Offset] += CorrectionAmountRight;
-			
-		      }else{ /* if( SUBlingGrid == FALSE ){ */
-            fprintf(stderr,"CFR: Doin It Right\n");
-			  CorrectionAmountRight = 400;//
-			    InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-			    RefinedFluxes->LeftFluxes[field][dim][RefinedFluxIndex];
-			  BaryonField[field][FieldIndex + Offset] += CorrectionAmountRight;
-			
-		      } // else{ /* if( SUBlingGrid == FALSE ){ */
-		    } // if(CorrectRightBaryonField)
-		
-		    if ((FieldTypeIsDensity(FieldType[field]) == TRUE ||
-			 FieldType[field] == TotalEnergy ||
-			 FieldType[field] == InternalEnergy)) {
-
-		      /* If new density & energy is < 0 then undo the
-			 flux correction. */
-
-		      if (CorrectLeftBaryonField &&
-			  BaryonField[field][FieldIndex] <= 0) {
-			BaryonField[field][FieldIndex] -= CorrectionAmountLeft;
-
-			if (SUBlingGrid == FALSE) {
-			  if (debug)
-			    printf("P(%d) -- CFRFl warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->LeftFluxes[field][dim][FluxIndex],
-				   CorrectionAmountLeft,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->LeftFluxes[ffield][dim][FluxIndex] =
-			      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-			} else {
-			  if (debug)
-			    printf("P(%d) -- CFRFlS warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->RightFluxes[field][dim][RefinedFluxIndex],
-				   CorrectionAmountLeft,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->RightFluxes[ffield][dim][RefinedFluxIndex] =
-			      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-			} // ENDELSE (SUBlingGrid == FALSE)
-		      } // ENDIF CorrectLeftBaryonField
-
-		      if (CorrectRightBaryonField &&
-			  BaryonField[field][FieldIndex+Offset] <= 0) {
-			BaryonField[field][FieldIndex+Offset] -= CorrectionAmountRight;
-
-			if (SUBlingGrid == FALSE) {
-			  if (debug)
-			    printf("P(%d) -- CFRFr warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->RightFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->RightFluxes[field][dim][FluxIndex],
-				   CorrectionAmountRight,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->RightFluxes[ffield][dim][FluxIndex] =
-			      InitialFluxes->RightFluxes[ffield][dim][FluxIndex];
-			} else {
-			  if (debug)
-			    printf("P(%d) -- CFRFrS warn: %e %e %e %e %"ISYM
-				   " %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-				   MyProcessorNumber, BaryonField[field][FieldIndex],
-				   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-				   RefinedFluxes->RightFluxes[field][dim][RefinedFluxIndex],
-				   CorrectionAmountRight,
-				   i, j, k, dim, field);
-			  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-			    RefinedFluxes->LeftFluxes[ffield][dim][RefinedFluxIndex] =
-			      InitialFluxes->RightFluxes[ffield][dim][FluxIndex];
-			}
-		      } // ENDIF CorrectRightBaryonField
-		    }
-		  }// for (i = Start[0]; i <= End[0]; i++) {
-		} // for (j = Start[1]; j <= End[1]; j++){
-	      } // for (k = Start[2]; k <= End[2]; k++){
-	    }	// if ((RadiativeCooling == 0 || (FieldType[field] != TotalEnergy && etc
-	  } // for (field = 0; field < NumberOfBaryonFields; field++){
-	  }
-
-
-	if (Coordinate == Cylindrical) {
-	FLOAT xr, xl, xc, geofacr, geofacl;
-	for (field = 0; field < NumberOfBaryonFields; field++) {
-	  /*if ((RadiativeCooling == 0 || (FieldType[field] != TotalEnergy && 
-	    FieldType[field] != InternalEnergy))
-	    && (FieldType[field] < ElectronDensity) ) {*/
-	  if (FieldType[field] < ElectronDensity &&
-	      FieldType[field] != DrivingField1 &&
-	      FieldType[field] != DrivingField2 &&
-	      FieldType[field] != DrivingField3 &&
-	      FieldType[field] != GravPotential) {
-	  for (k = Start[2]; k <= End[2]; k++) {
-	    for (j = Start[1]; j <= End[1]; j++) {
-	      for (i = Start[0]; i <= End[0]; i++) {
-		/* Compute indexes. */
-		FieldIndex = (k*GridDimension[1] + j)*GridDimension[0] + i;
-		FluxIndex  = ((k - Start[2])*Dim[1] + (j - Start[1]))*Dim[0] +
-		              (i - Start[0]);
-		GridFluxIndex = 
-                     (i - GridFluxStartIndex[0]) 
-		   + (j - GridFluxStartIndex[1])*GridFluxDim[0]
-		   + (k - GridFluxStartIndex[2])*GridFluxDim[1]*GridFluxDim[0];
-		if (dim == 0) {
-		  /* Left side */
-		  xr = CellLeftEdge[0][i] + CellWidth[0][i];
-                  xc = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
-                  geofacr = xr/xc;
-		  BaryonField[field][FieldIndex] +=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] )*geofacr;
-
-		  /* Right side */
-		  xl = CellLeftEdge[0][i+Offset];
-                  xc = xl + 0.5*CellWidth[0][i+Offset];
-                  geofacl = xl/xc;
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] )*geofacl;
-		}		
-		if (dim == 1) {
-		  /* Left side */
-		  BaryonField[field][FieldIndex] +=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] );
-		  /* Right side */
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] );
-		}		
-		if (dim == 2) {
-
-		  /* Left side */
-		  xc = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; 
-		  BaryonField[field][FieldIndex] +=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] )*xc;
-		  /* Right side */
-		  xl = CellLeftEdge[0][i+Offset];
-                  xc = xl + 0.5*CellWidth[0][i+Offset];
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] )*xc;
-
-		}		
-		/* Check for posivtivity and undo flux correction if negative */
-		if ((FieldType[field] == Density || 
-		     FieldType[field] == TotalEnergy ||
-		     FieldType[field] == InternalEnergy) &&
-		    BaryonField[field][FieldIndex] <= 0) {
-		  /*if (debug) {
-		    printf("CFRFl warn: %e %e %e %d %d %d %d [%d]\n",
-			   BaryonField[field][FieldIndex],
-			   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-			   RefinedFluxes->LeftFluxes[field][dim][FluxIndex],
-			   i, j, k, dim, field);
-			   }*/
-		  /* If new density & energy is < 0 then undo the flux correction. */
-		  BaryonField[field][FieldIndex] -=
-		     (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		      RefinedFluxes->LeftFluxes[field][dim][FluxIndex] );
-		  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-		    RefinedFluxes->LeftFluxes[ffield][dim][FluxIndex] =
-		      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-		}
-		if ((FieldType[field] == Density || 
-		     FieldType[field] == TotalEnergy ||
-		     FieldType[field] == InternalEnergy) &&
-		    BaryonField[field][FieldIndex + Offset] <= 0.0) {
-		  /*if (debug) {
-		    printf("CFRFr warn: %e %e %e %d %d %d %d (%d) [%d]\n",
-			   BaryonField[field][FieldIndex + Offset],
-			   InitialFluxes->RightFluxes[field][dim][FluxIndex],
-			   RefinedFluxes->RightFluxes[field][dim][FluxIndex],
-			   i, j, k, dim, Offset, field);
-			   }*/
-		  /* If new density & energy is < 0 then undo the flux correction. */
-		  BaryonField[field][FieldIndex + Offset] +=
-		     (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		      RefinedFluxes->RightFluxes[field][dim][FluxIndex] );
-		  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-		    RefinedFluxes->RightFluxes[ffield][dim][FluxIndex] =
-		      InitialFluxes->RightFluxes[ffield][dim][FluxIndex];
-		}
-	      }
-	    }
-	  }
-	  
-	  }
-	}
-	}
-	
-
-	    /* Return faces to original quantity. */
-	
-	  if (HydroMethod != Zeus_Hydro)
-	    for (field = 0; field < NumberOfBaryonFields; field++)
-	      if (FieldTypeIsDensity(FieldType[field]) == FALSE &&
-		  FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
-		  FieldTypeIsRadiation(FieldType[field]) == FALSE &&
-		  (RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-					     FieldType[field] != InternalEnergy)) && 
-		  FieldType[field] != Bfield1 &&
-		  FieldType[field] != Bfield2 && FieldType[field] != Bfield3 &&
-		  FieldType[field] != PhiField &&
-		  FieldType[field] != DrivingField1 &&
-		  FieldType[field] != DrivingField2 &&
-		  FieldType[field] != DrivingField3 &&
-		  FieldType[field] != GravPotential)
-		for (k = Start[2]; k <= End[2]; k++)
-		  for (j = Start[1]; j <= End[1]; j++) {
-		    index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		    for (i = Start[0]; i <= End[0]; i++, index++) {
-		      BaryonField[field][index] /= BaryonField[DensNum][index];
-		      BaryonField[field][index+Offset] /=
-			BaryonField[DensNum][index+Offset];
-		    }
-		  }
-	
-	  /* If appropriate, restore consistency between total and internal
-	     energy in corrected faces. */
-	
-	  if (DualEnergyFormalism == TRUE && RadiativeCooling == FALSE){
-	    float B2 = 0.0;
-	    for (k = Start[2]; k <= End[2]; k++)
-	      for (j = Start[1]; j <= End[1]; j++) {
-		i1 = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		i2 = i1 + Offset;
-		for (i = Start[0]; i <= End[0]; i++, i1++, i2++) {
-		  BaryonField[GENum][i1] = max(BaryonField[GENum][i1],
-					       tiny_number);
-		  BaryonField[GENum][i2] = max(BaryonField[GENum][i2],
-					       tiny_number);
-		  BaryonField[TENum][i1] = BaryonField[GENum][i1] +
-		    0.5 * BaryonField[Vel1Num][i1] * BaryonField[Vel1Num][i1];
-		  BaryonField[TENum][i2] = BaryonField[GENum][i2] +
-		    0.5 * BaryonField[Vel1Num][i2] * BaryonField[Vel1Num][i2];
-		  if (GridRank > 1) {
-		    BaryonField[TENum][i1] +=
-		      0.5 * BaryonField[Vel2Num][i1] * BaryonField[Vel2Num][i1];
-		    BaryonField[TENum][i2] +=
-		      0.5 * BaryonField[Vel2Num][i2] * BaryonField[Vel2Num][i2];
-		  }
-		  if (GridRank > 2) {
-		    BaryonField[TENum][i1] +=
-		      0.5 * BaryonField[Vel3Num][i1] * BaryonField[Vel3Num][i1];
-		    BaryonField[TENum][i2] +=
-		      0.5 * BaryonField[Vel3Num][i2] * BaryonField[Vel3Num][i2];
-		  }
-
-		  if (HydroMethod == MHD_RK) {
-		    B2 = POW(BaryonField[B1Num][i1],2) + 
-		      POW(BaryonField[B2Num][i1],2) +
-		      POW(BaryonField[B3Num][i1],2);
-		    BaryonField[TENum][i1] += 0.5 * B2 / BaryonField[DensNum][i1];
-		    B2 = POW(BaryonField[B1Num][i2],2) + 
-		      POW(BaryonField[B2Num][i2],2) +
-		      POW(BaryonField[B3Num][i2],2);
-		    BaryonField[TENum][i2] += 0.5 * B2 / BaryonField[DensNum][i2];
-		}
-
-		
-		}		
-	      } // end: loop over faces
-	  } // end: if (DualEnergyFormalism)
-	
-	    /* Multiply species by density to return from fractional to real
-	       density. (see comments above regarding species). */
-	
-	  for (field = 0; field < NumberOfBaryonFields; field++)
-	    if (FieldType[field] >= ElectronDensity &&
-		FieldType[field] < Metallicity &&
-		FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
-		FieldTypeIsRadiation(FieldType[field]) == FALSE)
-	      for (k = Start[2]; k <= End[2]; k++)
-		for (j = Start[1]; j <= End[1]; j++) {
-		  index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		  for (i = Start[0]; i <= End[0]; i++, index++) {
-		    BaryonField[field][index] *= BaryonField[DensNum][index];
-		    BaryonField[field][index+Offset] *=
-		      BaryonField[DensNum][index+Offset];
-		  }
-		}
-	
-	} // if( CorrectLeftBaryonField || CorrectRightBaryonField){
-
-      } // end: if GridDimension[dim] > 1
-      /* delete Refined fluxes as they're not needed anymore. */
- 
-      for (field = 0; field < NumberOfBaryonFields; field++) {
-	delete [] RefinedFluxes->LeftFluxes[field][dim];
-	delete [] RefinedFluxes->RightFluxes[field][dim];
-	RefinedFluxes->LeftFluxes[field][dim] = NULL;
-	RefinedFluxes->RightFluxes[field][dim] = NULL;
-      }
- 
-    } // next dimension
-  } // Number of baryons fields > 0
- 
-  return SUCCESS;
- 
-}
-#endif 

File src/enzo/kludges/Grid_CorrectForRefinedFluxes.C

-#ifndef FLUX_FIX
-/***********************************************************************
-/
-/  GRID CLASS (CORRECT SOLUTION GIVEN ORIGINAL AND REFINED FLUXES)
-/
-/  written by: Greg Bryan
-/  date:       November, 1994
-/  modified1:  Robert Harkness
-/  date:       January, 2003
-/              Include extra fields beyond Metallicity!
-/  modified2: David Collins, 2005
-/              Updated algebra so Cosmological Expansion is also
-/              conservative.  This fix also came with fixes to euler.src and
-/              Grid_GetProjectedBoundaryFluxes.C, so make sure you get those.
-/
-/  PURPOSE:
-/
-/  RETURNS: SUCCESS or FAIL
-/
-************************************************************************/
- 
-// Given both the original and refined fluxes for a subgrid in the current
-//   grid, correct the solution in the cells just outside the solution to
-//   reflect the new fluxes (which are determined from the solution of
-//   the subgrid).
-//   Note that if the subgrid is on the boundary of the current grid, we
-//     do not correct the values but instead replace the boundary fluxes
-//     for the current time step (BoundaryFluxesThisTimeStep).
- 
-#include <stdio.h>
-#include <math.h>
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
- 
-/* function prototypes */
- 
-int FindField(int f, int farray[], int n);
-int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
- 
- 
-int grid::CorrectForRefinedFluxes(fluxes *InitialFluxes,
-				  fluxes *RefinedFluxes,
-				  fluxes *BoundaryFluxesThisTimeStep)
-{
-  /* Return if this doesn't concern us. */
- 
-  if (ProcessorNumber != MyProcessorNumber || !UseHydro)
-    return SUCCESS;
- 
-  /* declarations */
- 
-  int i1, i2, i, j, k, dim, field, ffield, index;
-  int FieldIndex, FluxIndex, GridFluxIndex;
-  int CorrectLeftBoundaryFlux, CorrectRightBoundaryFlux, Offset;
-  int Dim[MAX_DIMENSION], End[MAX_DIMENSION], Start[MAX_DIMENSION];
-  int GridFluxDim[MAX_DIMENSION], GridFluxStartIndex[MAX_DIMENSION];
-  float B2;
- 
-  /* If there are no fields, don't do anything. */
- 
-  if (NumberOfBaryonFields > 0) {
- 
-    /* Find fields: density, total energy, velocity1-3. */
- 
-    int DensNum, GENum, Vel1Num, Vel2Num, Vel3Num, TENum, B1Num, B2Num, B3Num;
-    if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num,
-					 Vel3Num, TENum, B1Num, B2Num, B3Num) == FAIL) {
-      ENZO_FAIL("Error in grid->IdentifyPhysicalQuantities.\n");
-    }
- 
-    /* If using comoving coordinates, compute a(t) because we'll need it
-       to multiply against the CellWidth. */
-
-//    DC revision 16th September 2005 
-//    FLOAT a = 1, dadt;
-//    if (ComovingCoordinates)
-//      if (CosmologyComputeExpansionFactor(Time, &a, &dadt) == FAIL) {
-//        ENZO_FAIL("Error in CosmologyComputeExpansionFactors.\n");
-//      }
- 
-    /* Main loop over all faces. */
- 
-    for (dim = 0; dim < GridRank; dim++) {
- 
-      /* If this dimension is flat, don't do any work. */
- 
-      if (GridDimension[dim] > 1) {
- 
-	/* Check that the dims of InitialFluxes & RefinedFluxes are the same */
- 
-	for (j = 0; j < GridRank; j++)
-	  if ((InitialFluxes->LeftFluxStartGlobalIndex[dim][j] !=
-	       RefinedFluxes->LeftFluxStartGlobalIndex[dim][j])  ||
-	      (InitialFluxes->LeftFluxEndGlobalIndex[dim][j] !=
-	       RefinedFluxes->LeftFluxEndGlobalIndex[dim][j])) {
-	    ENZO_FAIL("InitialFluxes & RefinedFluxes are different.\n");
-	  }
- 
-	/* Error check Fluxes to make sure they all exist. */
- 
-	for (field = 0; field < NumberOfBaryonFields; field++)
-	  if ((InitialFluxes->LeftFluxes[field][dim] == NULL) ||
-	      (RefinedFluxes->LeftFluxes[field][dim] == NULL) ||
-	      (InitialFluxes->RightFluxes[field][dim] == NULL) ||
-	      (RefinedFluxes->RightFluxes[field][dim] == NULL)) {
-	    ENZO_FAIL("Some Flux data is not present.\n");
-	  }
- 
-	/* Compute Start and end indicies of flux region (with respect to
-	   the current grid's flux region). */
- 
-	for (i = 0; i < MAX_DIMENSION; i++) {
-	  Start[i] = 0;
-	  End[i] = 0;
-	}
- 
-	/* start index = subgrid flux left edge global index -
-	                 grid far left edge global index
-	     end index = subgrid flux right edge -
-	                 grid far left edge global index. */
- 
-	for (i = 0; i < GridRank; i++) {
-	  Start[i] = InitialFluxes->LeftFluxStartGlobalIndex[dim][i] -
-	    nlongint((CellLeftEdge[i][0] - DomainLeftEdge[i])/CellWidth[i][0]);
-	  End[i] = InitialFluxes->LeftFluxEndGlobalIndex[dim][i] -
-	    nlongint((CellLeftEdge[i][0] - DomainLeftEdge[i])/CellWidth[i][0]);
-	  if (Start[i] < 0 || End[i] > GridDimension[i]) {
-	    fprintf(stderr, "Start/End[%"ISYM"] = %"ISYM"/%"ISYM"\n", dim, Start[i], End[i]);
-	    fprintf(stderr, "%"GOUTSYM" %"GOUTSYM" %ld\n",
-		    CellLeftEdge[i][0], CellWidth[i][0],
-		    InitialFluxes->LeftFluxStartGlobalIndex[dim][i]);
-	    ENZO_FAIL("Error in Grid_CorrectForRefinedFluxes!\n");
-	  }
-	}
- 
-	/* Correct vector to point at cells just outside the left face.
-	   Start[dim] and End[dim] should be the same because the
-	   layer to be corrected is but one cell thick. */
- 
-	Start[dim] = max(Start[dim] - 1, 0);
-	End[dim]   = Start[dim];
- 
-	/* Compute Dimensions of Fluxes */
- 
-	for (i = 0; i < MAX_DIMENSION; i++)
-	  Dim[i] = End[i] - Start[i] + 1;
- 
-	/* Compute Offset (in baryon field) for right side of face.
-	   The +2 is there because we want to correct the cells just the
-	   right face.*/
- 
-	Offset = InitialFluxes->RightFluxStartGlobalIndex[dim][dim] -
-	         InitialFluxes->LeftFluxStartGlobalIndex[dim][dim] + 2;
-	Offset = min(Offset, GridDimension[dim]-1);  // this isn't needed (?)
- 
-	/* Check to see if we should correct BoundaryFluxesThisTimeStep
-	   instead of the fields themselves. */
- 
-	CorrectLeftBoundaryFlux = FALSE;
-	CorrectRightBoundaryFlux = FALSE;
-#ifdef UNUSED
-	if (Start[dim] == GridStartIndex[dim]-1)
-	  CorrectLeftBoundaryFlux = TRUE;
-	if (Start[dim] + Offset == GridEndIndex[dim]+1)
-	  CorrectRightBoundaryFlux = TRUE;
-#endif /* UNUSED */
- 
-	/* Set GridFluxStartIndex to the starting position of the flux
-	   plane (i.e. exclude grid boundary zones), except for the direction
-	   of the flux is set such that GridStartIndex[dim] - Start[dim] = 0 */
-	
-	for (i = 0; i < MAX_DIMENSION; i++) {
-	  GridFluxStartIndex[i] = GridStartIndex[i];
-	  GridFluxDim[i] = GridEndIndex[i] - GridStartIndex[i] + 1;
-	}
-	GridFluxStartIndex[dim] = Start[dim];
-	GridFluxDim[dim] = 1;
-	
-	/* Turn Offset (in dim direction) to Offset (in field array) */
- 
-	for (i = 0; i < dim; i++)
-	  Offset *= GridDimension[i];
- 
-	/* Multiply faces by density to get conserved quantities
-	   (only multiply fields which we are going to correct) */
- 
-	if (HydroMethod != Zeus_Hydro)
-	  for (field = 0; field < NumberOfBaryonFields; field++)
-	    if (FieldTypeIsDensity(FieldType[field]) == FALSE &&
-		FieldTypeIsRadiation(FieldType[field]) == FALSE &&
-		FieldType[field] != Bfield1 &&
-		FieldType[field] != Bfield2 && FieldType[field] != Bfield3 &&
-		FieldType[field] != PhiField &&
-		FieldType[field] != DrivingField1 &&
-		FieldType[field] != DrivingField2 &&
-		FieldType[field] != DrivingField3 &&
-		FieldType[field] != GravPotential)  
-	      //		(RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-	      //	 			 FieldType[field] != InternalEnergy)))
-	      for (k = Start[2]; k <= End[2]; k++)
-		for (j = Start[1]; j <= End[1]; j++) {
-		  index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		  for (i = Start[0]; i <= End[0]; i++, index++) {
-		    BaryonField[field][index] *= BaryonField[DensNum][index];
-		    BaryonField[field][index+Offset] *=
-		      BaryonField[DensNum][index+Offset];
-		  }
-		}
- 
-	/* Divide species by densities so that at the end we can multiply
-           them by the new density (species are not otherwise modified --
-	   see the next comment).  This ensures that the species are changed
-	   to keep the same fractional density. */
- 
-	for (field = 0; field < NumberOfBaryonFields; field++)
-	  if (FieldType[field] >= ElectronDensity &&
-	      FieldType[field] <= Metallicity &&
-	      FieldTypeIsRadiation(FieldType[field]) == FALSE)
-	    for (k = Start[2]; k <= End[2]; k++)
-	      for (j = Start[1]; j <= End[1]; j++) {
-		index = (k*GridDimension[1] + j)*GridDimension[0] + Start[0];
-		for (i = Start[0]; i <= End[0]; i++, index++) {
-		  BaryonField[field][index] /= BaryonField[DensNum][index];
-		  BaryonField[field][index+Offset] /=
-		    BaryonField[DensNum][index+Offset];
-		}
-	      }
-	
-	/* Correct face for difference between refined and initial fluxes.
-	   (Don't do this for energy if radiative cooling is on because it's
-	    no longer conserved.  Similarly, don't do it for the species
-	    because they are not individually conserved either -- in fact,
-            this could be done for the conserved quantities like charge,
-	    total number density summed over ionization, etc.) */
- 
-	for (field = 0; field < NumberOfBaryonFields; field++)
- 	  if (FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
- 	      (RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-					 FieldType[field] != InternalEnergy))
-	      && (FieldType[field] < ElectronDensity)) {
-	  for (k = Start[2]; k <= End[2]; k++)
-	    for (j = Start[1]; j <= End[1]; j++)
-	      for (i = Start[0]; i <= End[0]; i++) {
- 
-		/* Compute indexes. */
- 
-		FieldIndex = (k*GridDimension[1] + j)*GridDimension[0] + i;
-		FluxIndex  = ((k - Start[2])*Dim[1] + (j - Start[1]))*Dim[0] +
-		              (i - Start[0]);
-		GridFluxIndex =
-                     (i - GridFluxStartIndex[0])
-		   + (j - GridFluxStartIndex[1])*GridFluxDim[0]
-		   + (k - GridFluxStartIndex[2])*GridFluxDim[1]*GridFluxDim[0];
- 
-		/* Left side */
- 
-		if (CorrectLeftBoundaryFlux)
-		  /*		  BoundaryFluxesThisTimeStep->LeftFluxes[field][dim][GridFluxIndex] =
-		    RefinedFluxes->LeftFluxes[field][dim][FluxIndex]; */
-		  ;
-		else
-		  BaryonField[field][FieldIndex] +=
-		     (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		      RefinedFluxes->LeftFluxes[field][dim][FluxIndex] );
- 
-		if ((FieldTypeIsDensity(FieldType[field]) == TRUE ||
-		     FieldType[field] == TotalEnergy ||
-		     FieldType[field] == InternalEnergy) &&
-		    BaryonField[field][FieldIndex] <= 0) {
-		  if (debug)
-		    printf("CFRFl warn: %e %e %e %"ISYM" %"ISYM" %"ISYM" %"ISYM" [%"ISYM"]\n",
-			   BaryonField[field][FieldIndex],
-			   InitialFluxes->LeftFluxes[field][dim][FluxIndex],
-			   RefinedFluxes->LeftFluxes[field][dim][FluxIndex],
-			   i, j, k, dim, field);
-		  
-		  /* If new density is < 0 then stop the flux correction. */
- 
-		  BaryonField[field][FieldIndex] -=
-		    (InitialFluxes->LeftFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->LeftFluxes[field][dim][FluxIndex] );
-
-		  for (ffield = 0; ffield < NumberOfBaryonFields; ffield++)
-		    RefinedFluxes->LeftFluxes[ffield][dim][FluxIndex] =
-		      InitialFluxes->LeftFluxes[ffield][dim][FluxIndex];
-		  
-		  //ENZO_FAIL("New density or energy is < 0!\n");
-		}
- 
-		/* Right side */
- 
-		if (CorrectRightBoundaryFlux)
-		  /*  BoundaryFluxesThisTimeStep->RightFluxes[field][dim] [GridFluxIndex] =
-		  RefinedFluxes->RightFluxes[field][dim][FluxIndex]; */
-		  ;
-		else
-		  BaryonField[field][FieldIndex + Offset] -=
-		    (InitialFluxes->RightFluxes[field][dim][FluxIndex] -
-		     RefinedFluxes->RightFluxes[field][dim][FluxIndex] );
-
- 
-		if ((FieldTypeIsDensity(FieldType[field]) == TRUE ||
-		     FieldType[field] == TotalEnergy ||
-		     FieldType[field] == InternalEnergy) &&
-		    BaryonField[field][FieldIndex + Offset] <= 0.0) {
-		  if (debug)
-		    printf("CFRFr warn: %e %e %e %"ISYM" %"ISYM" %"ISYM" %"ISYM" (%"ISYM") [%"ISYM"]\n",
-			   BaryonField[field][FieldIndex + Offset],
-			   InitialFluxes->RightFluxes[field][dim][FluxIndex],
-			   RefinedFluxes->RightFluxes[field][dim][FluxIndex],
-			   i, j, k, dim, Offset, field);