Skip to content

Commit

Permalink
Radiative transfer adaptive timestep is now an option. If not set to
Browse files Browse the repository at this point in the history 
1, a constant timestep will be used.  If undefined, it will calculate
dtPhoton on the first timestep.

--HG--
branch : week-of-code
  • Loading branch information
John Wise committed Oct 12, 2009
1 parent 9d1c097 commit db33dbc
Show file tree
Hide file tree
Showing 12 changed files with 126 additions and 88 deletions.
10 changes: 9 additions & 1 deletion src/enzo/Grid_TransportPhotonPackages.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {

Expand All @@ -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++;
Expand Down
10 changes: 8 additions & 2 deletions src/enzo/Grid_WalkPhotonPackage.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -789,6 +793,8 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP,

// are we done ?
if (((*PP)->CurrentTime) >= EndTime) {
(*PP)->Photons = -1;
DeleteMe = TRUE;
return SUCCESS;
}

Expand Down
91 changes: 46 additions & 45 deletions src/enzo/PhotonGrid_Methods.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

/************************************************************************/

Expand Down
17 changes: 12 additions & 5 deletions src/enzo/RadiativeTransferComputeTimestep.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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->
Expand Down Expand Up @@ -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;

Expand Down
5 changes: 3 additions & 2 deletions src/enzo/RadiativeTransferInitialize.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 6 additions & 0 deletions src/enzo/RadiativeTransferParameters.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
2 changes: 1 addition & 1 deletion src/enzo/RadiativeTransferPrepare.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */

Expand Down
7 changes: 6 additions & 1 deletion src/enzo/RadiativeTransferReadParameters.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -59,6 +59,7 @@ int RadiativeTransferReadParameters(FILE *fptr)
RadiativeTransferTimestepVelocityLimit = 100.0; // km/s
RadiativeTransferPeriodicBoundary = FALSE;
RadiativeTransferHIIRestrictedTimestep = FALSE;
RadiativeTransferAdaptiveTimestep = FALSE;

/* read input from file */

Expand Down Expand Up @@ -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 */

Expand Down
6 changes: 5 additions & 1 deletion src/enzo/RadiativeTransferWriteParameters.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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;
}
2 changes: 1 addition & 1 deletion src/enzo/ReadParameterFile.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
1 change: 1 addition & 0 deletions src/enzo/RestartPhotons.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
57 changes: 28 additions & 29 deletions src/enzo/WriteParameterFile.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit db33dbc

Please sign in to comment.