diff --git a/src/enzo/Grid_TransportPhotonPackages.C b/src/enzo/Grid_TransportPhotonPackages.C index 109bf63bf..405a5b80f 100644 --- a/src/enzo/Grid_TransportPhotonPackages.C +++ b/src/enzo/Grid_TransportPhotonPackages.C @@ -163,7 +163,14 @@ int grid::TransportPhotonPackages(int level, ListOfPhotonsToMove **PhotonsToMove int AdvancePhotonPointer; int DeleteMe, DeltaLevel, PauseMe; - FLOAT EndTime = PhotonTime+dtPhoton-ROUNDOFF; + const float clight = 2.9979e10; + float LightCrossingTime = VelocityUnits / + (clight * RadiativeTransferPropagationSpeedFraction); + FLOAT EndTime; + if (RadiativeTransferAdaptiveTimestep) + EndTime = PhotonTime+LightCrossingTime; + else + EndTime = PhotonTime+dtPhoton-ROUNDOFF; while (PP != NULL) { @@ -179,6 +186,7 @@ int grid::TransportPhotonPackages(int level, ListOfPhotonsToMove **PhotonsToMove kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3, DeleteMe, PauseMe, DeltaLevel, + LightCrossingTime, DensityUnits, TemperatureUnits, VelocityUnits, LengthUnits, TimeUnits); tcount++; diff --git a/src/enzo/Grid_WalkPhotonPackage.C b/src/enzo/Grid_WalkPhotonPackage.C index 0beac67be..3034bbeab 100644 --- a/src/enzo/Grid_WalkPhotonPackage.C +++ b/src/enzo/Grid_WalkPhotonPackage.C @@ -57,7 +57,7 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, int gammaHeINum, int kphHeIINum, int gammaHeIINum, int kdissH2INum, int RPresNum1, int RPresNum2, int RPresNum3, int &DeleteMe, int &PauseMe, - int &DeltaLevel, + int &DeltaLevel, float LightCrossingTime, float DensityUnits, float TemperatureUnits, float VelocityUnits, float LengthUnits, float TimeUnits) { @@ -248,7 +248,11 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, FLOAT nH, nHI, nHeI, nHeII, fH, nH2I, nHI_inv, nHeI_inv, nHeII_inv, xe; FLOAT thisDensity, inverse_rho; float shield1, shield2, solid_angle, filling_factor, midpoint, nearest_edge; - EndTime = PhotonTime+dtPhoton; + + if (RadiativeTransferAdaptiveTimestep) + EndTime = PhotonTime + LightCrossingTime; + else + EndTime = PhotonTime + dtPhoton; int dummy; FLOAT dr; @@ -789,6 +793,8 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, // are we done ? if (((*PP)->CurrentTime) >= EndTime) { + (*PP)->Photons = -1; + DeleteMe = TRUE; return SUCCESS; } diff --git a/src/enzo/PhotonGrid_Methods.h b/src/enzo/PhotonGrid_Methods.h index 8b650de96..5ef83d0af 100644 --- a/src/enzo/PhotonGrid_Methods.h +++ b/src/enzo/PhotonGrid_Methods.h @@ -260,61 +260,62 @@ int MergePausedPhotonPackages(void); /* Trace a line thorugh the grid */ - int TraceRay(int NumberOfSegments, - FLOAT r, - FLOAT x, FLOAT y, FLOAT z, - FLOAT ux, FLOAT uy, FLOAT uz, - FLOAT dr[], - long cindex[], int ci[], int cj[], int ck[]); +int TraceRay(int NumberOfSegments, + FLOAT r, + FLOAT x, FLOAT y, FLOAT z, + FLOAT ux, FLOAT uy, FLOAT uz, + FLOAT dr[], + long cindex[], int ci[], int cj[], int ck[]); /* Walk Photon Package one by one */ - int WalkPhotonPackage(PhotonPackageEntry **PP, - grid **MoveToGrid, grid *ParentGrid, grid *CurrentGrid, - grid **Grids0, int nGrids0, - int DensNum, int HINum, int HeINum, - int HeIINum, int H2INum, - int kphHINum, int gammaHINum, - int kphHeINum, int gammaHeINum, - int kphHeIINum, int gammaHeIINum, - int kdissH2INum, int RPresNum1, int RPresNum2, - int RPresNum3, int &DeleteMe, int &PauseMe, - int &DeltaLevel, float DensityUnits, - float TemperatureUnits, float VelocityUnits, - float LengthUnits, float TimeUnits); +int WalkPhotonPackage(PhotonPackageEntry **PP, + grid **MoveToGrid, grid *ParentGrid, grid *CurrentGrid, + grid **Grids0, int nGrids0, + int DensNum, int HINum, int HeINum, + int HeIINum, int H2INum, + int kphHINum, int gammaHINum, + int kphHeINum, int gammaHeINum, + int kphHeIINum, int gammaHeIINum, + int kdissH2INum, int RPresNum1, int RPresNum2, + int RPresNum3, int &DeleteMe, int &PauseMe, + int &DeltaLevel, float LightCrossingTime, + float DensityUnits, + float TemperatureUnits, float VelocityUnits, + float LengthUnits, float TimeUnits); /* Create PhotonPackages for a given radiation sources */ - int Shine(RadiationSourceEntry *RadiationSource); +int Shine(RadiationSourceEntry *RadiationSource); /* PhotonTest: Initialize grid allowing for up to ten sources */ #define MAX_SPHERES 10 - int PhotonTestInitializeGrid(int NumberOfSpheres, - float SphereRadius[MAX_SPHERES], - float SphereCoreRadius[MAX_SPHERES], - float SphereDensity[MAX_SPHERES], - float SphereTemperature[MAX_SPHERES], - FLOAT SpherePosition[MAX_SPHERES][MAX_DIMENSION], - float SphereVelocity[MAX_SPHERES][MAX_DIMENSION], - float SphereFracKeplarianRot[MAX_SPHERES], - float SphereTurbulence[MAX_SPHERES], - float SphereCutOff[MAX_SPHERES], - float SphereAng1[MAX_SPHERES], - float SphereAng2[MAX_SPHERES], - int SphereNumShells[MAX_SPHERES], - int SphereType[MAX_SPHERES], - int SphereUseParticles, - float UniformVelocity[MAX_DIMENSION], - int SphereUseColour, - float InitialTemperature, int level, - float PhotonTestInitialFractionHII, - float PhotonTestInitialFractionHeII, - float PhotonTestInitialFractionHeIII, - float PhotonTestInitialFractionHM, - float PhotonTestInitialFractionH2I, - float PhotonTestInitialFractionH2II, - int RefineByOpticalDepth); +int PhotonTestInitializeGrid(int NumberOfSpheres, + float SphereRadius[MAX_SPHERES], + float SphereCoreRadius[MAX_SPHERES], + float SphereDensity[MAX_SPHERES], + float SphereTemperature[MAX_SPHERES], + FLOAT SpherePosition[MAX_SPHERES][MAX_DIMENSION], + float SphereVelocity[MAX_SPHERES][MAX_DIMENSION], + float SphereFracKeplarianRot[MAX_SPHERES], + float SphereTurbulence[MAX_SPHERES], + float SphereCutOff[MAX_SPHERES], + float SphereAng1[MAX_SPHERES], + float SphereAng2[MAX_SPHERES], + int SphereNumShells[MAX_SPHERES], + int SphereType[MAX_SPHERES], + int SphereUseParticles, + float UniformVelocity[MAX_DIMENSION], + int SphereUseColour, + float InitialTemperature, int level, + float PhotonTestInitialFractionHII, + float PhotonTestInitialFractionHeII, + float PhotonTestInitialFractionHeIII, + float PhotonTestInitialFractionHM, + float PhotonTestInitialFractionH2I, + float PhotonTestInitialFractionH2II, + int RefineByOpticalDepth); /************************************************************************/ diff --git a/src/enzo/RadiativeTransferComputeTimestep.C b/src/enzo/RadiativeTransferComputeTimestep.C index 8186baaeb..00e3b898c 100644 --- a/src/enzo/RadiativeTransferComputeTimestep.C +++ b/src/enzo/RadiativeTransferComputeTimestep.C @@ -42,6 +42,9 @@ int RadiativeTransferComputeTimestep(LevelHierarchyEntry *LevelArray[], int level) { + if (!RadiativeTransferAdaptiveTimestep && dtPhoton != FLOAT_UNDEFINED) + return SUCCESS; + const int MaxStepsPerHydroStep = 8; const float PhotonCourantFactor = 1.0; @@ -69,14 +72,14 @@ int RadiativeTransferComputeTimestep(LevelHierarchyEntry *LevelArray[], dtPhoton = 10*huge_number; + float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, + VelocityUnits = 1, TimeUnits = 1; + GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, + &VelocityUnits, LevelArray[level]->GridData->ReturnTime()); + // Calculate timestep by limiting to a 50% max change in HII if (RadiativeTransferHIIRestrictedTimestep) { - float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, - VelocityUnits = 1, TimeUnits = 1; - GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, - &VelocityUnits, LevelArray[level]->GridData->ReturnTime()); - for (l = 0; l < MAX_DEPTH_OF_HIERARCHY-1; l++) for (Temp = LevelArray[l]; Temp; Temp = Temp->NextGridThisLevel) { ThisPhotonDT = Temp->GridData-> @@ -119,6 +122,10 @@ int RadiativeTransferComputeTimestep(LevelHierarchyEntry *LevelArray[], if (InitialTimestep && !MetaData->FirstTimestepAfterRestart) dtPhoton = min(dtPhoton, dtLevelAbove); + + if (!RadiativeTransferAdaptiveTimestep && debug) + printf("RadiativeTransfer: Setting dtPhoton = %g = %g years\n", + dtPhoton, dtPhoton*TimeUnits/3.1557e7); return SUCCESS; diff --git a/src/enzo/RadiativeTransferInitialize.C b/src/enzo/RadiativeTransferInitialize.C index 2f5207cc1..69186bae3 100644 --- a/src/enzo/RadiativeTransferInitialize.C +++ b/src/enzo/RadiativeTransferInitialize.C @@ -140,8 +140,9 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, "from %"ISYM" to %"ISYM"\n", OldNumberOfBaryonFields, OldNumberOfBaryonFields+FieldsToAdd); - if (OldNumberOfBaryonFields+FieldsToAdd > MAX_DEPTH_OF_HIERARCHY) - ENZO_FAIL("Exceeds MAX_DEPTH_OF_HIERARCHY. Please increase and re-compile."); + if (OldNumberOfBaryonFields+FieldsToAdd > MAX_NUMBER_OF_BARYON_FIELDS) + ENZO_FAIL("Exceeds MAX_NUMBER_OF_BARYON_FIELDS. " + "Please increase and re-compile."); for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel) diff --git a/src/enzo/RadiativeTransferParameters.h b/src/enzo/RadiativeTransferParameters.h index f1ca71b32..2540979da 100644 --- a/src/enzo/RadiativeTransferParameters.h +++ b/src/enzo/RadiativeTransferParameters.h @@ -62,3 +62,9 @@ EXTERN int RadiativeTransferPeriodicBoundary; /* Flag to use timestepping to restrict HII fraction change to 50% */ EXTERN int RadiativeTransferHIIRestrictedTimestep; + +/* Flag to use adaptive timestepping for ray tracing. However this + removes the time-derivative from the radiative transfer + equation. */ + +EXTERN int RadiativeTransferAdaptiveTimestep; diff --git a/src/enzo/RadiativeTransferPrepare.C b/src/enzo/RadiativeTransferPrepare.C index 12b2ec8bb..a1c9d68ed 100644 --- a/src/enzo/RadiativeTransferPrepare.C +++ b/src/enzo/RadiativeTransferPrepare.C @@ -47,7 +47,7 @@ int RadiativeTransferPrepare(LevelHierarchyEntry *LevelArray[], int level, dt = LevelArray[level]->GridData->ReturnTimeStep(); //if (dtPhoton >= 0.0 && GridTime+dt >= PhotonTime) { - if (GridTime+dt >= PhotonTime) { + if (GridTime+dt >= PhotonTime || MetaData->FirstTimestepAfterRestart) { /* Determine the photon timestep */ diff --git a/src/enzo/RadiativeTransferReadParameters.C b/src/enzo/RadiativeTransferReadParameters.C index 7872d1fb8..91513b88e 100644 --- a/src/enzo/RadiativeTransferReadParameters.C +++ b/src/enzo/RadiativeTransferReadParameters.C @@ -32,7 +32,7 @@ int RadiativeTransferReadParameters(FILE *fptr) RadiationPressure = FALSE; // off PhotonTime = 0; - dtPhoton = 0.1; + dtPhoton = FLOAT_UNDEFINED; for (i = 0; i < 4; i++) { EscapedPhotonCount[i] = 0.0; TotalEscapedPhotonCount[i] = 0.0; @@ -59,6 +59,7 @@ int RadiativeTransferReadParameters(FILE *fptr) RadiativeTransferTimestepVelocityLimit = 100.0; // km/s RadiativeTransferPeriodicBoundary = FALSE; RadiativeTransferHIIRestrictedTimestep = FALSE; + RadiativeTransferAdaptiveTimestep = FALSE; /* read input from file */ @@ -100,6 +101,10 @@ int RadiativeTransferReadParameters(FILE *fptr) &RadiativeTransferPhotonMergeRadius); ret += sscanf(line, "RadiativeTransferHIIRestrictedTimestep = %"ISYM, &RadiativeTransferHIIRestrictedTimestep); + ret += sscanf(line, "RadiativeTransferAdaptiveTimestep = %"ISYM, + &RadiativeTransferAdaptiveTimestep); + + ret += sscanf(line, "dtPhoton = %"FSYM, dtPhoton); /* if the line is suspicious, issue a warning */ diff --git a/src/enzo/RadiativeTransferWriteParameters.C b/src/enzo/RadiativeTransferWriteParameters.C index 07a8be870..32ab6eff5 100644 --- a/src/enzo/RadiativeTransferWriteParameters.C +++ b/src/enzo/RadiativeTransferWriteParameters.C @@ -25,6 +25,8 @@ int RadiativeTransferWriteParameters(FILE *fptr) { + fprintf(fptr, "dtPhoton = %"GOUTSYM"\n", + dtPhoton); fprintf(fptr, "RadiativeTransferRadiationPressure = %"ISYM"\n", RadiationPressure); fprintf(fptr, "RadiativeTransferSourceRadius = %"GSYM"\n", @@ -55,8 +57,10 @@ int RadiativeTransferWriteParameters(FILE *fptr) RadiativeTransferSourceClustering); fprintf(fptr, "RadiativeTransferPhotonMergeRadius = %"FSYM"\n", RadiativeTransferPhotonMergeRadius); - fprintf(fptr, "RadiativeTransferHIIRestrictedTimestep = %"ISYM"\n\n", + fprintf(fptr, "RadiativeTransferHIIRestrictedTimestep = %"ISYM"\n", RadiativeTransferHIIRestrictedTimestep); + fprintf(fptr, "RadiativeTransferAdaptiveTimestep = %"ISYM"\n\n", + RadiativeTransferAdaptiveTimestep); return SUCCESS; } diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C index df938494d..4e49c73d3 100644 --- a/src/enzo/ReadParameterFile.C +++ b/src/enzo/ReadParameterFile.C @@ -1146,7 +1146,7 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt) my_exit(EXIT_SUCCESS); #endif } - printf("Initialdt in ReadParameterFiled %e\n", *Initialdt); + if (debug) printf("Initialdt in ReadParameterFile = %e\n", *Initialdt); CheckShearingBoundaryConsistency(MetaData); return SUCCESS; } diff --git a/src/enzo/RestartPhotons.C b/src/enzo/RestartPhotons.C index 56c8ba91b..41098896b 100644 --- a/src/enzo/RestartPhotons.C +++ b/src/enzo/RestartPhotons.C @@ -73,6 +73,7 @@ int RestartPhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], float LightCrossingTime = VelocityUnits / (clight * RadiativeTransferPropagationSpeedFraction); FLOAT SavedPhotonTime = PhotonTime; + FLOAT SavedPhotonTimestep = dtPhoton; PhotonTime -= LightCrossingTime; if (debug) diff --git a/src/enzo/WriteParameterFile.C b/src/enzo/WriteParameterFile.C index 24387e44b..ba11a5e54 100644 --- a/src/enzo/WriteParameterFile.C +++ b/src/enzo/WriteParameterFile.C @@ -583,39 +583,38 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData) /* Most Stanford additions: */ - fprintf(fptr, "Theta_Limiter = %f\n", Theta_Limiter); - fprintf(fptr, "RiemannSolver = %d\n", RiemannSolver); - fprintf(fptr, "ReconstructionMethod = %d\n", ReconstructionMethod); - fprintf(fptr, "RKOrder = %d\n", RKOrder); - fprintf(fptr, "UsePhysicalUnit = %d\n", UsePhysicalUnit); - fprintf(fptr, "UseFloor = %d\n", UseFloor); - fprintf(fptr, "UseViscosity = %d\n", UseViscosity); - fprintf(fptr, "UseAmbipolarDiffusion = %d\n", UseAmbipolarDiffusion); - fprintf(fptr, "UseResistivity = %d\n", UseResistivity); - fprintf(fptr, "SmallRho = %g\n", SmallRho*rhou); - fprintf(fptr, "SmallP = %g\n", SmallP*presu); - fprintf(fptr, "SmallT = %g\n", SmallT*tempu); - fprintf(fptr, "MaximumAlvenSpeed = %g\n", MaximumAlvenSpeed*velu); - fprintf(fptr, "Coordinate = %d\n", Coordinate); - fprintf(fptr, "EOSType = %d\n", EOSType); - fprintf(fptr, "EOSSoundSpeed = %g\n", EOSSoundSpeed); - fprintf(fptr, "EOSCriticalDensity = %g\n", EOSCriticalDensity); - fprintf(fptr, "EOSGamma = %g\n", EOSGamma); - fprintf(fptr, "Mu = %g\n", Mu); - fprintf(fptr, "CoolingCutOffDensity1 = %g\n", CoolingCutOffDensity1); - fprintf(fptr, "CoolingCutOffDensity2 = %g\n", CoolingCutOffDensity2); - fprintf(fptr, "CoolingCutOffTemperature = %g\n", CoolingCutOffTemperature); + fprintf(fptr, "Theta_Limiter = %f\n", Theta_Limiter); + fprintf(fptr, "RiemannSolver = %d\n", RiemannSolver); + fprintf(fptr, "ReconstructionMethod = %d\n", ReconstructionMethod); + fprintf(fptr, "RKOrder = %d\n", RKOrder); + fprintf(fptr, "UsePhysicalUnit = %d\n", UsePhysicalUnit); + fprintf(fptr, "UseFloor = %d\n", UseFloor); + fprintf(fptr, "UseViscosity = %d\n", UseViscosity); + fprintf(fptr, "UseAmbipolarDiffusion = %d\n", UseAmbipolarDiffusion); + fprintf(fptr, "UseResistivity = %d\n", UseResistivity); + fprintf(fptr, "SmallRho = %g\n", SmallRho*rhou); + fprintf(fptr, "SmallP = %g\n", SmallP*presu); + fprintf(fptr, "SmallT = %g\n", SmallT*tempu); + fprintf(fptr, "MaximumAlvenSpeed = %g\n", MaximumAlvenSpeed*velu); + fprintf(fptr, "Coordinate = %d\n", Coordinate); + fprintf(fptr, "EOSType = %d\n", EOSType); + fprintf(fptr, "EOSSoundSpeed = %g\n", EOSSoundSpeed); + fprintf(fptr, "EOSCriticalDensity = %g\n", EOSCriticalDensity); + fprintf(fptr, "EOSGamma = %g\n", EOSGamma); + fprintf(fptr, "Mu = %g\n", Mu); + fprintf(fptr, "CoolingCutOffDensity1 = %g\n", CoolingCutOffDensity1); + fprintf(fptr, "CoolingCutOffDensity2 = %g\n", CoolingCutOffDensity2); + fprintf(fptr, "CoolingCutOffTemperature = %g\n", CoolingCutOffTemperature); fprintf(fptr, "CoolingPowerCutOffDensity1 = %g\n", CoolingPowerCutOffDensity1); fprintf(fptr, "CoolingPowerCutOffDensity2 = %g\n", CoolingPowerCutOffDensity2); - fprintf(fptr, "UseConstantAcceleration = %d\n", UseConstantAcceleration); - fprintf(fptr, "ConstantAcceleration = %g %g %g\n", ConstantAcceleration[0], + fprintf(fptr, "UseConstantAcceleration = %d\n", UseConstantAcceleration); + fprintf(fptr, "ConstantAcceleration = %g %g %g\n", ConstantAcceleration[0], ConstantAcceleration[1], ConstantAcceleration[2]); - - fprintf(fptr, "AngularVelocity = %g\n", AngularVelocity); - fprintf(fptr, "VelocityGradient = %g\n", VelocityGradient); - fprintf(fptr, "UseDrivingField = %d\n", UseDrivingField); - fprintf(fptr, "DrivingEfficiency = %f\n", DrivingEfficiency); + fprintf(fptr, "AngularVelocity = %g\n", AngularVelocity); + fprintf(fptr, "VelocityGradient = %g\n", VelocityGradient); + fprintf(fptr, "UseDrivingField = %d\n", UseDrivingField); + fprintf(fptr, "DrivingEfficiency = %f\n", DrivingEfficiency); #ifdef ECUDA fprintf(fptr, "UseCUDA = %f\n", UseCUDA); #endif