From 99a77d5d8bc29b678ccab31b9c33a07c1706bbb1 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Fri, 17 Sep 2021 08:36:46 -0400 Subject: [PATCH 01/10] Adding PopIII_Rotating parameter --- src/enzo/ReadParameterFile.C | 3 +++ src/enzo/SetDefaultGlobalValues.C | 3 +++ src/enzo/WriteParameterFile.C | 3 +++ src/enzo/global_data.h | 3 +++ 4 files changed, 12 insertions(+) diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C index 4996fa6db..ed29869ed 100644 --- a/src/enzo/ReadParameterFile.C +++ b/src/enzo/ReadParameterFile.C @@ -1368,6 +1368,9 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt) ret += sscanf(line,"MagneticSupernovaEnergy = %"FSYM, &MagneticSupernovaEnergy); ret += sscanf(line,"MagneticSupernovaDuration = %"FSYM, &MagneticSupernovaDuration); + // Rotating Pop III Models + ret += sscanf(line, "PopIII_Rotating = %"ISYM, &PopIII_Rotating); + /* If the dummy char space was used, then make another. */ if (*dummy != 0) { diff --git a/src/enzo/SetDefaultGlobalValues.C b/src/enzo/SetDefaultGlobalValues.C index 6cae8caa1..8cb4879ec 100644 --- a/src/enzo/SetDefaultGlobalValues.C +++ b/src/enzo/SetDefaultGlobalValues.C @@ -1050,5 +1050,8 @@ int SetDefaultGlobalValues(TopGridData &MetaData) MagneticSupernovaDuration = 5e4 ; // Total duration of magnetic feedback in years MagneticSupernovaEnergy = 1.0e51; // Total energy (ergs) injected per star particle (supernova) + /* Rotating Pop III Stars Model */ + PopIII_Rotating = 0; // 0 = off + return SUCCESS; } diff --git a/src/enzo/WriteParameterFile.C b/src/enzo/WriteParameterFile.C index e4ce7134b..4e459eb87 100644 --- a/src/enzo/WriteParameterFile.C +++ b/src/enzo/WriteParameterFile.C @@ -1274,6 +1274,9 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData, char *name = NULL) fprintf(fptr,"MagneticSupernovaEnergy = %"GSYM"\n",MagneticSupernovaEnergy); fprintf(fptr,"MagneticSupernovaDuration = %"GSYM"\n",MagneticSupernovaDuration); + /* Rotating PopIII Stars Model */ + fprintf(fptr, "PopIII_Rotating = %"ISYM"\n", PopIII_Rotating); + /* Output current time */ time_t ID; ID = time(NULL); diff --git a/src/enzo/global_data.h b/src/enzo/global_data.h index 92fa5a51c..1be6c0872 100644 --- a/src/enzo/global_data.h +++ b/src/enzo/global_data.h @@ -1205,4 +1205,7 @@ EXTERN float MagneticSupernovaRadius; EXTERN float MagneticSupernovaDuration; EXTERN float MagneticSupernovaEnergy; +/* Rotating Pop III Stars Model */ +EXTERN int PopIII_Rotating; + #endif From 24c706396fad942558a82f85eb51b93285b49ce8 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Fri, 1 Oct 2021 15:07:44 -0400 Subject: [PATCH 02/10] Making Test Star Particle a radiating P3 star. --- src/enzo/Grid.h | 13 +- .../Grid_TestStarParticleInitializeGrid.C | 306 +++++++++++++++++- src/enzo/Make.config.objects | 3 +- src/enzo/Make.mach.linux-gnu | 21 +- src/enzo/PopIIIRotationModel.C | 138 ++++++++ src/enzo/Star.h | 3 + src/enzo/Star_ComputePhotonRates.C | 47 ++- src/enzo/TestStarParticleInitialize.C | 133 +++++--- 8 files changed, 597 insertions(+), 67 deletions(-) create mode 100644 src/enzo/PopIIIRotationModel.C diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 4757742a1..14831002c 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -2070,7 +2070,18 @@ int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], int TestStarParticleInitializeGrid(float TestStarParticleStarMass, float *Initialdt, FLOAT TestStarParticleStarVelocity[], - FLOAT TestStarParticleStarPosition[]); + FLOAT TestStarParticleStarPosition[], + float UniformDensity, + float UniformTotalEnergy, + float UniformVelocity[], + float UniformBField[], + float InitialTemperature, + float PhotonTestInitialFractionHII, + float PhotonTestInitialFractionHeII, + float PhotonTestInitialFractionHeIII, + float PhotonTestInitialFractionHM, + float PhotonTestInitialFractionH2I, + float PhotonTestInitialFractionH2II); /* Gravity Test: initialize grid. */ diff --git a/src/enzo/Grid_TestStarParticleInitializeGrid.C b/src/enzo/Grid_TestStarParticleInitializeGrid.C index 8c3e2cbb9..a6900148c 100644 --- a/src/enzo/Grid_TestStarParticleInitializeGrid.C +++ b/src/enzo/Grid_TestStarParticleInitializeGrid.C @@ -14,6 +14,7 @@ #include #include +#include #include #include "ErrorExceptions.h" #include "macros_and_parameters.h" @@ -29,27 +30,144 @@ int GetUnits(float *DensityUnits, float *LengthUnits, float *TemperatureUnits, float *TimeUnits, float *VelocityUnits, double *MassUnits, FLOAT Time); +// Used to compute Bonner-Ebert density profile +double ph_BE(double r); +double ph_q(double r); +double ph_Ang(double a1, double a2, double R, double r); + +// Returns random velocity from Maxwellian distribution +double ph_Maxwellian(double c_tilda, double vel_unit, double mu, double gamma); + +static int PhotonTestParticleCount = 0; + int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, float *Initialdt, FLOAT TestStarParticleStarVelocity[], - FLOAT TestStarParticleStarPosition[]) + FLOAT TestStarParticleStarPosition[], + float TestStarParticleDensity, + float TestStarParticleEnergy, + float TestStarParticleVelocity[], + float TestStarParticleBField[], + float InitialTemperature, + float PhotonTestInitialFractionHII, + float PhotonTestInitialFractionHeII, + float PhotonTestInitialFractionHeIII, + float PhotonTestInitialFractionHM, + float PhotonTestInitialFractionH2I, + float PhotonTestInitialFractionH2II) { /* declarations */ float CentralMass = 1.0; - int i, dim; - float TestInitialdt = *Initialdt; + int dim, i, j, k, m, field, size, active_size, index, cindex; + + int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, + DINum, DIINum, HDINum, MetalNum, MetalIaNum, B1Num, B2Num, B3Num, PhiNum, CRNum; + + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum, kdissH2IINum, kphHMNum, RPresNum1, RPresNum2, RPresNum3; + + int ExtraField[2]; + + float TestInitialdt = *Initialdt, *density_field = NULL, *HII_field = NULL, *HeII_field = NULL, + *HeIII_field = NULL, *Temperature_field = NULL; /* Return if this doesn't concern us. */ if (ProcessorNumber != MyProcessorNumber) return SUCCESS; + /* create fields */ + + NumberOfBaryonFields = 0; + FieldType[NumberOfBaryonFields++] = Density; + FieldType[NumberOfBaryonFields++] = TotalEnergy; + if (DualEnergyFormalism) + FieldType[NumberOfBaryonFields++] = InternalEnergy; + int vel = NumberOfBaryonFields; + FieldType[NumberOfBaryonFields++] = Velocity1; + if (GridRank > 1 || HydroMethod > 2) + FieldType[NumberOfBaryonFields++] = Velocity2; + if (GridRank > 2 || HydroMethod > 2) + FieldType[NumberOfBaryonFields++] = Velocity3; + + if ( CRModel ) { + CRNum = NumberOfBaryonFields; + FieldType[NumberOfBaryonFields++] = CRDensity; + } + + if (WritePotential) + FieldType[NumberOfBaryonFields++] = GravPotential; + + + int colorfields = NumberOfBaryonFields; + + // Enzo's standard multispecies (primordial chemistry - H, D, He) + if (TestProblemData.MultiSpecies) { + FieldType[DeNum = NumberOfBaryonFields++] = ElectronDensity; + FieldType[HINum = NumberOfBaryonFields++] = HIDensity; + FieldType[HIINum = NumberOfBaryonFields++] = HIIDensity; + FieldType[HeINum = NumberOfBaryonFields++] = HeIDensity; + FieldType[HeIINum = NumberOfBaryonFields++] = HeIIDensity; + FieldType[HeIIINum = NumberOfBaryonFields++] = HeIIIDensity; + if (TestProblemData.MultiSpecies > 1) { + FieldType[HMNum = NumberOfBaryonFields++] = HMDensity; + FieldType[H2INum = NumberOfBaryonFields++] = H2IDensity; + FieldType[H2IINum = NumberOfBaryonFields++] = H2IIDensity; + } + if (TestProblemData.MultiSpecies > 2) { + FieldType[DINum = NumberOfBaryonFields++] = DIDensity; + FieldType[DIINum = NumberOfBaryonFields++] = DIIDensity; + FieldType[HDINum = NumberOfBaryonFields++] = HDIDensity; + } + } + + // Metal fields, including the standard 'metallicity' as well + // as two extra fields + if (TestProblemData.UseMetallicityField) { + FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity; + + if (StarMakerTypeIaSNe) + FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity; + + if(TestProblemData.MultiMetals){ + FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0; + FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1; + } + } + + if (RadiativeTransfer && (MultiSpecies < 1)) { + ENZO_FAIL("Grid_TestStarParticleInitialize: Radiative Transfer but not MultiSpecies set"); + } + + // Allocate fields for photo ionization and heating rates + if (RadiativeTransfer) + if (MultiSpecies) { + FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; + FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; + if (RadiativeTransferHydrogenOnly == FALSE) { + FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; + FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; + } + if (MultiSpecies > 1) { + FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; + FieldType[kdissH2IINum = NumberOfBaryonFields++] = kdissH2II; + FieldType[kphHMNum = NumberOfBaryonFields++] = kphHM; + } + } + + if (RadiationPressure && RadiativeTransfer) { + FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0; + FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1; + FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2; + } + + NumberOfPhotonPackages = 0; + PhotonPackages-> NextPackage= NULL; /* Get Units. */ float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, - VelocityUnits = 1, TimeUnits = 1; + VelocityUnits = 1, TimeUnits = 1, CriticalDensity = 1, BoxLength = 1, mu = 0.6, mu_data; double MassUnits = 1; if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, @@ -70,12 +188,22 @@ int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, this->AllocateNewParticles(NumberOfParticles); printf("Allocated %d particles\n", NumberOfParticles); + this->AllocateGrids(); + + /* Initialize radiation fields */ // need this + + if (this->InitializeRadiativeTransferFields() == FAIL) { + ENZO_FAIL("\nError in InitializeRadiativeTransferFields.\n"); + } + /* Set particle IDs and types */ for (i = 0; i < NumberOfParticles; i++) { ParticleNumber[i] = i; - ParticleType[i] = PARTICLE_TYPE_STAR; + ParticleType[i] = -PopIII; + ParticleAttribute[0][i] = Time + 1e-7; } + ParticleMass[0] = CentralMass; /* Set central particle. */ for (dim = 0; dim < GridRank; dim++) { @@ -83,8 +211,7 @@ int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, (DomainLeftEdge[dim]+DomainRightEdge[dim]) + 0.5*CellWidth[0][0]; ParticleVelocity[dim][0] = TestStarParticleStarVelocity[dim]*1e5*TimeUnits/LengthUnits; } - ParticleMass[0] = CentralMass; - ParticleAttribute[0][0] = Time+1e-7; //creation time:make sure it is non-zero + if (STARFEED_METHOD(UNIGRID_STAR)) ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; if (STARFEED_METHOD(MOM_STAR)) if(StarMakerExplosionDelayTime >= 0.0) @@ -95,6 +222,171 @@ int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, ParticleAttribute[2][0] = 0.0; // Metal fraction ParticleAttribute[3][0] = 0.0; // metalfSNIa + if (this->FindNewStarParticles(GridLevel) == FAIL) { + ENZO_FAIL("Error in grid::FindNewStarParticles."); + } + + /* Reset particle type to be positive*/ + for (i = 0; i < NumberOfParticles; i++) { + ParticleNumber[i] = i; + ParticleType[i] = PopIII; + } + Star *cstar; + for (cstar = Stars; cstar; cstar = cstar->NextStar) + cstar->type = PopIII; + + + /* Set up the baryon field. */ // need thid + /* compute size of fields */ + active_size = 1; + int ActiveDims[MAX_DIMENSION]; + for (dim = 0; dim < GridRank; dim++) { + ActiveDims[dim] = GridEndIndex[dim] - GridStartIndex[dim] + 1; + active_size *= ActiveDims[dim]; + } + + /* Loop over the mesh. */ + float density, dens1, Velocity[MAX_DIMENSION], + temperature, temp1, sigma, sigma1, colour, outer_radius; + float HII_Fraction, HeII_Fraction, HeIII_Fraction, H2I_Fraction; + FLOAT r, x, y = 0, z = 0; + int n = 0; + + for (k = 0; k < GridDimension[2]; k++) + for (j = 0; j < GridDimension[1]; j++) + for (i = 0; i < GridDimension[0]; i++, n++) { + + /* Compute position */ // need this + + x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; + if (GridRank > 1) + y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j]; + if (GridRank > 2) + z = CellLeftEdge[2][k] + 0.5*CellWidth[2][k]; + + if (i >= GridStartIndex[0] && i <= GridEndIndex[0] && + j >= GridStartIndex[1] && j <= GridEndIndex[1] && + k >= GridStartIndex[2] && k <= GridEndIndex[2]) { + cindex = (i-GridStartIndex[0]) + ActiveDims[0] * + ((j-GridStartIndex[1]) + (k-GridStartIndex[2])*ActiveDims[1]); + if (density_field != NULL) + density = max(density_field[cindex], 1e-6); + else + density = 1.0; + if (HII_field != NULL) + HII_Fraction = HII_field[cindex]; + else + HII_Fraction = PhotonTestInitialFractionHII; + if (HeII_field != NULL) + HeII_Fraction = HeII_field[cindex]; + else + HeII_Fraction = PhotonTestInitialFractionHeII; + if (HeIII_field != NULL) + HeIII_Fraction = HeIII_field[cindex]; + else + HeIII_Fraction = PhotonTestInitialFractionHeIII; + if (Temperature_field != NULL) + temperature = temp1 = Temperature_field[cindex]; + else + temperature = temp1 = InitialTemperature; + } + + /* set density, total energy */ // check + BaryonField[0][n] = TestStarParticleDensity; + BaryonField[1][n] = TestStarParticleEnergy; + + /* set velocities */ + + for (dim = 0; dim < GridRank; dim++) + BaryonField[vel+dim][n] = Velocity[dim] + TestStarParticleVelocity[dim]; + + /* Set internal energy if necessary. */ //check + + if (DualEnergyFormalism) + BaryonField[2][n] = BaryonField[1][n]; + + /* If doing multi-species (HI, etc.), set these. */ + + if (MultiSpecies > 0) { + BaryonField[HIINum][n] = TestProblemData.HII_Fraction * + TestProblemData.HydrogenFractionByMass * TestStarParticleDensity; + + BaryonField[HeIINum][n] = TestProblemData.HeII_Fraction * + TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); + + BaryonField[HeIIINum][n] = TestProblemData.HeIII_Fraction * + TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); + + BaryonField[HeINum][n] = + (1.0 - TestProblemData.HydrogenFractionByMass)*TestStarParticleDensity - + BaryonField[HeIINum][n] - BaryonField[HeIIINum][n]; + } + + if (MultiSpecies > 1) { + BaryonField[HMNum][n] = TestProblemData.HM_Fraction * + BaryonField[HIINum][n]; + + BaryonField[H2INum][n] = TestProblemData.H2I_Fraction * + BaryonField[0][n] * TestProblemData.HydrogenFractionByMass; + + BaryonField[H2IINum][n] = TestProblemData.H2II_Fraction * 2.0 * + BaryonField[HIINum][n]; + + BaryonField[kdissH2INum][n] = 0.0; + BaryonField[kphHMNum][n] = 0.0; + BaryonField[kdissH2IINum][n] = 0.0; + } + + // HI density is calculated by subtracting off the various ionized fractions + // from the total + BaryonField[HINum][n] = TestProblemData.HydrogenFractionByMass*BaryonField[0][n] + - BaryonField[HIINum][n]; + + if (MultiSpecies > 1) + BaryonField[HINum][n] -= (BaryonField[HMNum][n] + BaryonField[H2IINum][n] + + BaryonField[H2INum][n]); + + BaryonField[DeNum][n] = BaryonField[HIINum][n] + + 0.25*BaryonField[HeIINum][n] + 0.5*BaryonField[HeIIINum][n]; + + if (MultiSpecies > 1) + BaryonField[DeNum][n] += 0.5*BaryonField[H2IINum][n] - + BaryonField[HMNum][n]; + + /* Set Deuterium species (assumed to be negligible). */ + + if (MultiSpecies > 2) { + BaryonField[DINum ][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HINum][n]; + BaryonField[DIINum][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HIINum][n]; + BaryonField[HDINum][n] = 0.75 * TestProblemData.DeuteriumToHydrogenRatio * BaryonField[H2INum][n]; + } + + /* Set energy (thermal and then total if necessary). */ + + if (MultiSpecies) { + mu_data = + 0.25*(BaryonField[HeINum][n] + BaryonField[HeIINum][n] + + BaryonField[HeIIINum][n] ) + + BaryonField[HINum][n] + BaryonField[HIINum][n] + + BaryonField[DeNum][n]; + if (MultiSpecies > 1) + mu_data += BaryonField[HMNum][n] + + 0.5*(BaryonField[H2INum][n] + BaryonField[H2IINum][n]); + mu_data = BaryonField[0][n] / mu_data; + } else + mu_data = mu; + + BaryonField[1][n] = temperature/TemperatureUnits/ + ((Gamma-1.0)*mu_data); + + } + + delete [] density_field; + delete [] HII_field; + delete [] HeII_field; + delete [] HeIII_field; + delete [] Temperature_field; + return SUCCESS; } diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects index 0f05a6669..ddcfd9244 100644 --- a/src/enzo/Make.config.objects +++ b/src/enzo/Make.config.objects @@ -1035,7 +1035,8 @@ OBJS_CONFIG_LIB = \ rk4.o \ SchrodingerAddPotential.o \ FDMCollapse.o \ - Grid_FDMCollapse.o + Grid_FDMCollapse.o \ + PopIIIRotationModel.o #----------------------------------------------------------------------- diff --git a/src/enzo/Make.mach.linux-gnu b/src/enzo/Make.mach.linux-gnu index 3d583d4e0..f9ec4422e 100644 --- a/src/enzo/Make.mach.linux-gnu +++ b/src/enzo/Make.mach.linux-gnu @@ -36,8 +36,8 @@ MACH_FILE = Make.mach.linux-gnu # Install paths (local variables) #----------------------------------------------------------------------- -LOCAL_HDF5_INSTALL = /PATH/TO/SERIAL-HDF5/INSTALL # mandatory -LOCAL_GRACKLE_INSTALL = /PATH/TO/GRACKLE/INSTALL # optional +LOCAL_HDF5_INSTALL = /home/dskinner6/software/yt-conda/pkgs/hdf5-1.10.2-hc401514_3# mandatory +LOCAL_GRACKLE_INSTALL = /home/dskinner6/software/grackle # optional LOCAL_HYPRE_INSTALL = /PATH/TO/HYPRE/INSTALL # optional #----------------------------------------------------------------------- @@ -49,10 +49,10 @@ MACH_CPP = cpp # C preprocessor command # With MPI MACH_CC_MPI = mpicc # C compiler when using MPI -MACH_CXX_MPI = mpic++ # C++ compiler when using MPI +MACH_CXX_MPI = mpicxx # C++ compiler when using MPI MACH_FC_MPI = gfortran # Fortran 77 compiler when using MPI MACH_F90_MPI = gfortran # Fortran 90 compiler when using MPI -MACH_LD_MPI = mpic++ # Linker when using MPI +MACH_LD_MPI = mpicxx # Linker when using MPI # Without MPI @@ -66,14 +66,15 @@ MACH_LD_NOMPI = g++ # Linker when not using MPI # Machine-dependent defines #----------------------------------------------------------------------- -MACH_DEFINES = -DLINUX -DH5_USE_16_API +MACH_DEFINES = -DLINUX -DH5_USE_16_API -fPIC #----------------------------------------------------------------------- # Compiler flag settings #----------------------------------------------------------------------- -MACH_CPPFLAGS = -P -traditional +MACH_OMPFLAGS = -fopenmp +MACH_CPPFLAGS = -P -traditional MACH_CFLAGS = MACH_CXXFLAGS = MACH_FFLAGS = -fno-second-underscore -ffixed-line-length-132 @@ -86,7 +87,7 @@ MACH_LDFLAGS = MACH_OPT_WARN = -Wall -g MACH_OPT_DEBUG = -g -MACH_OPT_HIGH = -O2 +MACH_OPT_HIGH = -O2 -g MACH_OPT_AGGRESSIVE = -O3 -g #----------------------------------------------------------------------- @@ -94,10 +95,10 @@ MACH_OPT_AGGRESSIVE = -O3 -g #----------------------------------------------------------------------- LOCAL_INCLUDES_MPI = # MPI includes -LOCAL_INCLUDES_HDF5 = -I$(LOCAL_HDF5_INSTALL) # HDF5 includes +LOCAL_INCLUDES_HDF5 = -I$(LOCAL_HDF5_INSTALL)/include # HDF5 includes LOCAL_INCLUDES_HYPRE = -I$(LOCAL_HYPRE_INSTALL)/include LOCAL_INCLUDES_PAPI = # PAPI includes -LOCAL_INCLUDES_GRACKLE = -I$(LOCAL_GRACKLE_INSTALL)/include +LOCAL_INCLUDES_GRACKLE = -I/home/dskinner6/software/grackle/include MACH_INCLUDES = $(LOCAL_INCLUDES_HDF5) MACH_INCLUDES_MPI = $(LOCAL_INCLUDES_MPI) @@ -114,7 +115,7 @@ LOCAL_LIBS_HDF5 = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz LOCAL_LIBS_HYPRE = -L$(LOCAL_HYPRE_INSTALL)/lib -lHYPRE LOCAL_LIBS_PAPI = # PAPI libraries LOCAL_LIBS_MACH = -lgfortran # Machine-dependent libraries -LOCAL_LIBS_GRACKLE = -L$(LOCAL_GRACKLE_INSTALL)/lib -lgrackle +LOCAL_LIBS_GRACKLE = -L/home/dskinner6/software/grackle/lib -lgrackle MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) MACH_LIBS_MPI = $(LOCAL_LIBS_MPI) diff --git a/src/enzo/PopIIIRotationModel.C b/src/enzo/PopIIIRotationModel.C new file mode 100644 index 000000000..bf5bfe843 --- /dev/null +++ b/src/enzo/PopIIIRotationModel.C @@ -0,0 +1,138 @@ +/*********************************************************************** +/ +/ COMPUTE PHOTON PRODUCTION RATES FOR (NON)ROTATING POPIII STARS +/ +/ written by: Danielle Skinner +/ date: September, 2021 +/ modified1: +/ +/ ---------- SPECIES -------- +/ 0 : HI +/ 1 : HeI +/ 2 : HeII +/ +/ Model comes from Murphy et al. 2021. Only valid for stars with masses +/ between 9 - 120 Msun. H2 not included. The fits here are up until the +/ star is at 80% of its lifetime, after which the photon production rate +/ is assumed to be constant. +/ +************************************************************************/ + +#include +#include +#include +#include "ErrorExceptions.h" +#include "phys_constants.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" + +int j, bin, nvalue; +float frac, a, b, c, Mass, age; +double Q, Q_lower, Q_upper; + +const int nentries = 9; +const float mass_list[] = {9.0, 12.0, 15.0, 20.0, 30.0, 40.0, 60.0, 85.0, 120.0}; + +int search_lower_bound(float *arr, float value, int low, int high, + int total); + +int find_element(float element, const float array[], int n){ + nvalue = -1; + for (j = 0; j < n; j++){ + if(array[j] == element){ + nvalue = j; + break; + } + } + return nvalue; +} + +double InterpolateQ(float Mass, float age, float *parameters_a, float *parameters_b, float *parameters_c){ + float age2 = age*age; + + if (find_element(Mass, (float*)mass_list, nentries) != -1) { + Q = (parameters_a[nvalue] * age2 + parameters_b[nvalue] * age + parameters_c[nvalue]); + return Q; + } + else { + bin = search_lower_bound((float*)mass_list, Mass, 0, nentries, nentries); + frac = (Mass - mass_list[bin]) / (mass_list[bin+1] - mass_list[bin]); + Q_lower = ( (parameters_a[bin] * age2) + (parameters_b[bin] * age) + parameters_c[bin] ); + Q_upper = ( (parameters_a[bin+1] * age2) + (parameters_b[bin+1] * age) + parameters_c[bin+1] ); + + Q = Q_lower + frac * (Q_upper - Q_lower); + return Q; + } +} + +double CalculateRotationalPhotonRates(float Mass, float age, int species){ + + if (species == 0) { + static float parameters_a[] = {-0.3014, 0.03601, 0.1357, 0.2404, 0.2116, 0.1807, 0.1118, 0.0901, 0.0683}; + static float parameters_b[] = {0.66228, 0.40061, 0.26264, 0.15633, 0.17335, 0.177001, 0.21019, 0.19481, 0.19028}; + static float parameters_c[] = {47.42321, 47.82289, 48.15684, 48.52598, 48.96249, 49.23924, 49.58557, 49.85428, 50.09629}; + + Q = InterpolateQ(Mass, age, parameters_a, parameters_b, parameters_c); + + } + + else if (species == 1) { + static float parameters_a[] = {-0.9481, -0.6126, -0.3906, -0.1218, -0.1451, -0.1122, -0.4014, -0.4259, -0.4637}; + static float parameters_b[] = {1.19688, 0.99590, 0.65148, 0.32491, 0.28988, 0.23664, 0.37460, 0.33801, 0.31139}; + static float parameters_c[] = {46.48224, 46.98303, 47.47341, 47.98433, 48.52012, 48.84545, 49.22893, 49.52746, 49.79099}; + + Q = InterpolateQ(Mass, age, parameters_a, parameters_b, parameters_c); + + } + + else if (species == 2) { + static float parameters_a[] = {-2.8125, -2.5140, -1.9560, -1.2121, -1.2286, -1.0128, -1.9777, -2.0174, -2.0865}; + static float parameters_b[] = {2.73813, 2.74062, 1.80823, 0.83063, 0.63568, 0.41093, 0.86035, 0.75522, 0.63998}; + static float parameters_c[] = {43.43041, 44.22116, 45.17038, 46.10727, 46.95021, 47.42831, 47.93338, 48.32898, 48.66499}; + + Q = InterpolateQ(Mass, age, parameters_a, parameters_b, parameters_c); + + } + return Q; +} + +double CalculateNonRotationalPhotonRates(float Mass, float age, int species){ + + if (species == 0) { + static float parameters_a[] = {-0.35938, -0.03421, 0.05179, 0.17966, 0.18772, 0.17035, 0.15019, 0.12873, 0.12268}; + static float parameters_b[] = {0.64365, 0.38223, 0.26300, 0.13555, 0.13217, 0.14523, 0.14671, 0.14156, 0.12346}; + static float parameters_c[] = {47.444944, 47.830253, 48.174224, 48.544643, 48.979915, 49.254436, 49.601645, 49.868288, 50.109797}; + + Q = InterpolateQ(Mass, age, parameters_a, parameters_b, parameters_c); + + } + + else if (species == 1) { + static float parameters_a[] = {-1.05193, -0.53651, -0.46799, -0.15949, -0.08082, -0.11586, -0.17175, -0.22875, -0.27835}; + static float parameters_b[] = {1.19633, 0.95788, 0.65528, 0.30155, 0.21199, 0.22530, 0.23914, 0.24920, 0.25542}; + static float parameters_c[] = {46.522791, 46.981776, 47.505882, 48.018253, 48.553852, 48.872077, 49.263237, 49.555736, 49.814517}; + + Q = InterpolateQ(Mass, age, parameters_a, parameters_b, parameters_c); + + } + + else if (species == 2) { + static float parameters_a[] = {-3.05209, -2.01521, -2.01975, -1.18593, -0.90444, -1.00130, -1.17825, -1.35364, -1.54739}; + static float parameters_b[] = {2.79275, 2.65228, 1.82661, 0.80239, 0.45104, 0.46494, 0.51727, 0.57479, 0.65844}; + static float parameters_c[] = {43.524763, 44.193070, 45.247383, 46.187693, 47.034429, 47.491437, 48.024623, 48.402114, 48.718482}; + + Q = InterpolateQ(Mass, age, parameters_a, parameters_b, parameters_c); + + } + + return Q; + +} \ No newline at end of file diff --git a/src/enzo/Star.h b/src/enzo/Star.h index 7ee4127b7..9e02d1bf0 100644 --- a/src/enzo/Star.h +++ b/src/enzo/Star.h @@ -68,6 +68,8 @@ class Star Star operator+=(Star a); Star* copy(void); + FLOAT Time; + // Routines star_type ReturnType(void) { return type; }; int ReturnID(void) { return Identifier; }; @@ -75,6 +77,7 @@ class Star float ReturnBirthTime(void) { return BirthTime; }; double ReturnFinalMass(void) { return FinalMass; }; void AssignFinalMass(double value) { FinalMass = value; }; + float ReturnAge(void) { return ( Time - BirthTime ); }; float ReturnLifetime(void) { return LifeTime; }; float ReturnBirthtime(void) { return BirthTime; }; int ReturnLevel(void) { return level; }; diff --git a/src/enzo/Star_ComputePhotonRates.C b/src/enzo/Star_ComputePhotonRates.C index 542678794..fb1fe898c 100644 --- a/src/enzo/Star_ComputePhotonRates.C +++ b/src/enzo/Star_ComputePhotonRates.C @@ -30,19 +30,22 @@ #include "LevelHierarchy.h" float ReturnValuesFromSpectrumTable(float ColumnDensity, float dColumnDensity, int mode); +float CalculateRotationalPhotonRates(float Mass, float age, int species); +float CalculateNonRotationalPhotonRates(float Mass, float age, int species); int Star::ComputePhotonRates(const float TimeUnits, int &nbins, float E[], double Q[]) { int i; - double L_UV, cgs_convert, _mass; + double L_UV, cgs_convert, _mass, _age; float x, x2, EnergyFractionLW, MeanEnergy, XrayLuminosityFraction; float Mform, EnergyFractionHeI, EnergyFractionHeII; if (this->Mass < 0.1) // Not "born" yet _mass = this->FinalMass; else _mass = this->Mass; + _age = this->ReturnAge()/this->LifeTime; //Percentage of lifetime x = log10((float)(_mass)); x2 = x*x; @@ -61,17 +64,45 @@ int Star::ComputePhotonRates(const float TimeUnits, int &nbins, float E[], doubl E[2] = 58.0; E[3] = 12.8; _mass = max(min((float)(_mass), 500), 5); - if (_mass > 9 && _mass <= 500) { - Q[0] = pow(10.0, 43.61 + 4.9*x - 0.83*x2); - Q[1] = pow(10.0, 42.51 + 5.69*x - 1.01*x2); - Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); - Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); - } else if (_mass > 5 && _mass <= 9) { + if (_mass > 5 && _mass < 9) { Q[0] = pow(10.0, 39.29 + 8.55*x); Q[1] = pow(10.0, 29.24 + 18.49*x); Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); - } // ENDELSE + } + else if (_mass >= 9 && _mass <= 120) { + if (PopIII_Rotating == 1) { + if (_age <= 0.8){ + Q[0] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 0)); + Q[1] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 1)); + Q[2] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 2)); + } + else { + Q[0] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 0)); + Q[1] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 1)); + Q[2] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 2)); + } + } + else { + if (_age <= 0.8){ + Q[0] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 0)); + Q[1] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 1)); + Q[2] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 2)); + } + else { + Q[0] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 0)); + Q[1] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 1)); + Q[2] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 2)); + } + } + Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); + } + else if (_mass > 120 && _mass < 500) { + Q[0] = pow(10.0, 43.61 + 4.9*x - 0.83*x2); + Q[1] = pow(10.0, 42.51 + 5.69*x - 1.01*x2); + Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); + Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); + } else { for (i = 0; i < nbins; i++) Q[i] = 0.0; } diff --git a/src/enzo/TestStarParticleInitialize.C b/src/enzo/TestStarParticleInitialize.C index 1a8147dfc..43c1b8d1e 100644 --- a/src/enzo/TestStarParticleInitialize.C +++ b/src/enzo/TestStarParticleInitialize.C @@ -30,28 +30,53 @@ #include "Hierarchy.h" #include "TopGridData.h" +static float PhotonTestInitialFractionHII = 1.2e-5; +static float PhotonTestInitialFractionHeII = 1.0e-14; +static float PhotonTestInitialFractionHeIII = 1.0e-17; +static float PhotonTestInitialFractionHM = 2.0e-9; +static float PhotonTestInitialFractionH2I = 2.0e-20; +static float PhotonTestInitialFractionH2II = 3.0e-14; + int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData,float *Initialdt) { - char *DensName = "Density"; - char *TEName = "TotalEnergy"; - char *GEName = "GasEnergy"; - char *Vel1Name = "x-velocity"; - char *Vel2Name = "y-velocity"; - char *Vel3Name = "z-velocity"; - char *MetalName = "Metal_Density"; - char *ElectronName = "Electron_Density"; - char *HIName = "HI_Density"; - char *HIIName = "HII_Density"; - char *HeIName = "HeI_Density"; - char *HeIIName = "HeII_Density"; - char *HeIIIName = "HeIII_Density"; + const char *DensName = "Density"; /* make const */ + const char *TEName = "TotalEnergy"; + const char *GEName = "GasEnergy"; + const char *Vel1Name = "x-velocity"; + const char *Vel2Name = "y-velocity"; + const char *Vel3Name = "z-velocity"; + const char *MetalName = "Metal_Density"; + const char *ElectronName = "Electron_Density"; + const char *HIName = "HI_Density"; + const char *HIIName = "HII_Density"; + const char *HeIName = "HeI_Density"; + const char *HeIIName = "HeII_Density"; + const char *HeIIIName = "HeIII_Density"; + const char *HMName = "HM_Density"; + const char *H2IName = "H2I_Density"; + const char *H2IIName = "H2II_Density"; + const char *DIName = "DI_Density"; + const char *DIIName = "DII_Density"; + const char *HDIName = "HDI_Density"; + const char *kphHIName = "HI_kph"; + const char *gammaName = "PhotoGamma"; + const char *kphHeIName = "HeI_kph"; + const char *kphHeIIName = "HeII_kph"; + const char *kphHMName = "HM_kph"; + const char *kdissH2IIName = "H2II_kdiss"; + const char *kdissH2IName = "H2I_kdiss"; + const char *RadAccel1Name = "x-RadPressure"; + const char *RadAccel2Name = "y-RadPressure"; + const char *RadAccel3Name = "z-RadPressure"; /* declarations */ char line[MAX_LINE_LENGTH]; - int dim, ret; + char *dummy = new char[MAX_LINE_LENGTH]; + int dim, ret, i, source; + dummy[0] = 0; /* Error check. */ @@ -67,7 +92,7 @@ int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGri float TestStarParticleStarMass = 100.0; int TestProblemUseMetallicityField = 1; float TestProblemInitialMetallicityFraction = 2e-3; // 0.1 Zsun - + float PhotonTestInitialTemperature = 1000; @@ -122,42 +147,70 @@ int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGri /* set up uniform grid as of before explosion */ - - - if (TopGrid.GridData->InitializeUniformGrid(TestStarParticleDensity, - TestStarParticleEnergy, - TestStarParticleEnergy, - TestStarParticleVelocity, - TestStarParticleBField) == FAIL) - ENZO_FAIL("Error in InitializeUniformGrid."); - if (TopGrid.GridData-> TestStarParticleInitializeGrid(TestStarParticleStarMass, Initialdt, TestStarParticleStarVelocity, - TestStarParticleStarPosition) == FAIL) - ENZO_FAIL("Error in TestStarParticleInitializeGrid.\n"); + TestStarParticleStarPosition, + TestStarParticleDensity, + TestStarParticleEnergy, + TestStarParticleVelocity, + TestStarParticleBField, + PhotonTestInitialTemperature, + PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII, + PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM, + PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II) == FAIL) + ENZO_FAIL("Error in TestStarParticleInitializeGrid.\n"); /* set up field names and units */ int count = 0; - DataLabel[count++] = DensName; - DataLabel[count++] = TEName; + DataLabel[count++] = (char*) DensName; + DataLabel[count++] = (char*) TEName; if (DualEnergyFormalism) - DataLabel[count++] = GEName; - DataLabel[count++] = Vel1Name; - DataLabel[count++] = Vel2Name; - DataLabel[count++] = Vel3Name; + DataLabel[count++] = (char*) GEName; + DataLabel[count++] = (char*) Vel1Name; + DataLabel[count++] = (char*) Vel2Name; + DataLabel[count++] = (char*) Vel3Name; if (TestProblemData.MultiSpecies){ - DataLabel[count++] = ElectronName; - DataLabel[count++] = HIName; - DataLabel[count++] = HIIName; - DataLabel[count++] = HeIName; - DataLabel[count++] = HeIIName; - DataLabel[count++] = HeIIIName; - } + DataLabel[count++] = (char*) ElectronName; + DataLabel[count++] = (char*) HIName; + DataLabel[count++] = (char*) HIIName; + DataLabel[count++] = (char*) HeIName; + DataLabel[count++] = (char*) HeIIName; + DataLabel[count++] = (char*) HeIIIName; + if (MultiSpecies > 1) { + DataLabel[count++] = (char*) HMName; + DataLabel[count++] = (char*) H2IName; + DataLabel[count++] = (char*) H2IIName; + } + if (MultiSpecies > 2) { + DataLabel[count++] = (char*) DIName; + DataLabel[count++] = (char*) DIIName; + DataLabel[count++] = (char*) HDIName; + } + } // if Multispecies if (TestProblemData.UseMetallicityField) - DataLabel[count++] = MetalName; + DataLabel[count++] = (char*) MetalName; + + if (RadiativeTransfer) + if (MultiSpecies) { + DataLabel[count++] = (char*) kphHIName; + DataLabel[count++] = (char*) gammaName; + DataLabel[count++] = (char*) kphHeIName; + DataLabel[count++] = (char*) kphHeIIName; + if (MultiSpecies > 1) { + DataLabel[count++]= (char*) kdissH2IName; + DataLabel[count++]= (char*) kdissH2IIName; + DataLabel[count++]= (char*) kphHMName; + } + } // if RadiativeTransfer + + if (RadiationPressure) { + DataLabel[count++] = (char*) RadAccel1Name; + DataLabel[count++] = (char*) RadAccel2Name; + DataLabel[count++] = (char*) RadAccel3Name; + } int j; From 232bff73e9c0395e0846082f243c64e5b27dbbaa Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Thu, 21 Jul 2022 09:08:55 -0400 Subject: [PATCH 03/10] Separating TestStarParticle and TestRadiatingStarParticle. New test added. --- src/enzo/Grid.h | 6 + ..._TestRadiatingStarParticleInitializeGrid.C | 393 ++++++++++++++++++ .../Grid_TestStarParticleInitializeGrid.C | 306 +------------- src/enzo/InitializeNew.C | 8 +- src/enzo/Make.config.objects | 2 + src/enzo/Make.mach.linux-gnu | 2 +- .../TestRadiatingStarParticleInitialize.C | 247 +++++++++++ src/enzo/TestStarParticleInitialize.C | 128 ++---- 8 files changed, 699 insertions(+), 393 deletions(-) create mode 100644 src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C create mode 100644 src/enzo/TestRadiatingStarParticleInitialize.C diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 14831002c..7072bf612 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -2068,6 +2068,12 @@ int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], /* Star Particle test: initialize particle */ int TestStarParticleInitializeGrid(float TestStarParticleStarMass, + float *Initialdt, + FLOAT TestStarParticleStarVelocity[], + FLOAT TestStarParticleStarPosition[]); + +/* Radiating Star Particle test: initialize particle */ + int TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass, float *Initialdt, FLOAT TestStarParticleStarVelocity[], FLOAT TestStarParticleStarPosition[], diff --git a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C new file mode 100644 index 000000000..c480da812 --- /dev/null +++ b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C @@ -0,0 +1,393 @@ +/*********************************************************************** +/ +/ GRID CLASS (INITIALIZE THE GRID FOR A RADIATING STAR PARTICLE TEST) +/ +/ written by: Greg Bryan +/ date: June, 2012 +/ modified1: Danielle Skinner +/ date: July, 2022 +/ +/ PURPOSE: +/ +/ RETURNS: FAIL or SUCCESS +/ +************************************************************************/ + +#include +#include +#include +#include +#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 "phys_constants.h" + +int GetUnits(float *DensityUnits, float *LengthUnits, + float *TemperatureUnits, float *TimeUnits, + float *VelocityUnits, double *MassUnits, FLOAT Time); + +// Used to compute Bonner-Ebert density profile +double ph_BE(double r); +double ph_q(double r); +double ph_Ang(double a1, double a2, double R, double r); + +// Returns random velocity from Maxwellian distribution +double ph_Maxwellian(double c_tilda, double vel_unit, double mu, double gamma); + +static int PhotonTestParticleCount = 0; + +int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass, + float *Initialdt, + FLOAT TestStarParticleStarVelocity[], + FLOAT TestStarParticleStarPosition[], + float TestStarParticleDensity, + float TestStarParticleEnergy, + float TestStarParticleVelocity[], + float TestStarParticleBField[], + float InitialTemperature, + float PhotonTestInitialFractionHII, + float PhotonTestInitialFractionHeII, + float PhotonTestInitialFractionHeIII, + float PhotonTestInitialFractionHM, + float PhotonTestInitialFractionH2I, + float PhotonTestInitialFractionH2II) +{ + /* declarations */ + + float CentralMass = 1.0; + int dim, i, j, k, m, field, size, active_size, index, cindex; + + int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, + DINum, DIINum, HDINum, MetalNum, MetalIaNum, B1Num, B2Num, B3Num, PhiNum, CRNum; + + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum, kdissH2IINum, kphHMNum, RPresNum1, RPresNum2, RPresNum3; + + int ExtraField[2]; + + float TestInitialdt = *Initialdt, *density_field = NULL, *HII_field = NULL, *HeII_field = NULL, + *HeIII_field = NULL, *Temperature_field = NULL; + + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) + return SUCCESS; + + /* create fields */ + + NumberOfBaryonFields = 0; + FieldType[NumberOfBaryonFields++] = Density; + FieldType[NumberOfBaryonFields++] = TotalEnergy; + if (DualEnergyFormalism) + FieldType[NumberOfBaryonFields++] = InternalEnergy; + int vel = NumberOfBaryonFields; + FieldType[NumberOfBaryonFields++] = Velocity1; + if (GridRank > 1 || HydroMethod > 2) + FieldType[NumberOfBaryonFields++] = Velocity2; + if (GridRank > 2 || HydroMethod > 2) + FieldType[NumberOfBaryonFields++] = Velocity3; + + if ( CRModel ) { + CRNum = NumberOfBaryonFields; + FieldType[NumberOfBaryonFields++] = CRDensity; + } + + if (WritePotential) + FieldType[NumberOfBaryonFields++] = GravPotential; + + + int colorfields = NumberOfBaryonFields; + + // Enzo's standard multispecies (primordial chemistry - H, D, He) + if (TestProblemData.MultiSpecies) { + FieldType[DeNum = NumberOfBaryonFields++] = ElectronDensity; + FieldType[HINum = NumberOfBaryonFields++] = HIDensity; + FieldType[HIINum = NumberOfBaryonFields++] = HIIDensity; + FieldType[HeINum = NumberOfBaryonFields++] = HeIDensity; + FieldType[HeIINum = NumberOfBaryonFields++] = HeIIDensity; + FieldType[HeIIINum = NumberOfBaryonFields++] = HeIIIDensity; + if (TestProblemData.MultiSpecies > 1) { + FieldType[HMNum = NumberOfBaryonFields++] = HMDensity; + FieldType[H2INum = NumberOfBaryonFields++] = H2IDensity; + FieldType[H2IINum = NumberOfBaryonFields++] = H2IIDensity; + } + if (TestProblemData.MultiSpecies > 2) { + FieldType[DINum = NumberOfBaryonFields++] = DIDensity; + FieldType[DIINum = NumberOfBaryonFields++] = DIIDensity; + FieldType[HDINum = NumberOfBaryonFields++] = HDIDensity; + } + } + + // Metal fields, including the standard 'metallicity' as well + // as two extra fields + if (TestProblemData.UseMetallicityField) { + FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity; + + if (StarMakerTypeIaSNe) + FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity; + + if(TestProblemData.MultiMetals){ + FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0; + FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1; + } + } + + if (RadiativeTransfer && (MultiSpecies < 1)) { + ENZO_FAIL("Grid_TestRadiatingStarParticleInitialize: Radiative Transfer but not MultiSpecies set"); + } + + // Allocate fields for photo ionization and heating rates + if (RadiativeTransfer) + if (MultiSpecies) { + FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; + FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; + if (RadiativeTransferHydrogenOnly == FALSE) { + FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; + FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; + } + if (MultiSpecies > 1) { + FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; + FieldType[kdissH2IINum = NumberOfBaryonFields++] = kdissH2II; + FieldType[kphHMNum = NumberOfBaryonFields++] = kphHM; + } + } + + if (RadiationPressure && RadiativeTransfer) { + FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0; + FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1; + FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2; + } + + NumberOfPhotonPackages = 0; + PhotonPackages-> NextPackage= NULL; + + /* Get Units. */ + + float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, + VelocityUnits = 1, TimeUnits = 1, CriticalDensity = 1, BoxLength = 1, mu = 0.6, mu_data; + double MassUnits = 1; + + if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, + &TimeUnits, &VelocityUnits, &MassUnits, Time) == FAIL) { + ENZO_FAIL("Error in GetUnits.\n"); + } + + /* Set Central Mass in simulation units */ + + CentralMass = TestStarParticleStarMass*1.99e33* pow(LengthUnits*CellWidth[0][0],-3.0)/DensityUnits; + + printf("Central Mass: %f \n",CentralMass); + + /* Set number of particles for this grid and allocate space. */ + + NumberOfParticles = 1; + NumberOfParticleAttributes = 4; + this->AllocateNewParticles(NumberOfParticles); + printf("Allocated %d particles\n", NumberOfParticles); + + this->AllocateGrids(); + + /* Initialize radiation fields */ // need this + + if (this->InitializeRadiativeTransferFields() == FAIL) { + ENZO_FAIL("\nError in InitializeRadiativeTransferFields.\n"); + } + + /* Set particle IDs and types */ + + for (i = 0; i < NumberOfParticles; i++) { + ParticleNumber[i] = i; + ParticleType[i] = -PopIII; + ParticleAttribute[0][i] = Time + 1e-7; + } + ParticleMass[0] = CentralMass; + + /* Set central particle. */ + for (dim = 0; dim < GridRank; dim++) { + ParticlePosition[dim][0] = TestStarParticleStarPosition[dim]* + (DomainLeftEdge[dim]+DomainRightEdge[dim]) + 0.5*CellWidth[0][0]; + ParticleVelocity[dim][0] = TestStarParticleStarVelocity[dim]*1e5*TimeUnits/LengthUnits; + } + + if (STARFEED_METHOD(UNIGRID_STAR)) ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; + if (STARFEED_METHOD(MOM_STAR)) + if(StarMakerExplosionDelayTime >= 0.0) + ParticleAttribute[1][0] = 1.0; + else + ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; + + ParticleAttribute[2][0] = 0.0; // Metal fraction + ParticleAttribute[3][0] = 0.0; // metalfSNIa + + if (this->FindNewStarParticles(GridLevel) == FAIL) { + ENZO_FAIL("Error in grid::FindNewStarParticles."); + } + + /* Reset particle type to be positive*/ + for (i = 0; i < NumberOfParticles; i++) { + ParticleNumber[i] = i; + ParticleType[i] = PopIII; + } + Star *cstar; + for (cstar = Stars; cstar; cstar = cstar->NextStar) + cstar->type = PopIII; + + + /* Set up the baryon field. */ // need thid + /* compute size of fields */ + active_size = 1; + int ActiveDims[MAX_DIMENSION]; + for (dim = 0; dim < GridRank; dim++) { + ActiveDims[dim] = GridEndIndex[dim] - GridStartIndex[dim] + 1; + active_size *= ActiveDims[dim]; + } + + /* Loop over the mesh. */ + float density, dens1, Velocity[MAX_DIMENSION], + temperature, temp1, sigma, sigma1, colour, outer_radius; + float HII_Fraction, HeII_Fraction, HeIII_Fraction, H2I_Fraction; + FLOAT r, x, y = 0, z = 0; + int n = 0; + + for (k = 0; k < GridDimension[2]; k++) + for (j = 0; j < GridDimension[1]; j++) + for (i = 0; i < GridDimension[0]; i++, n++) { + + /* Compute position */ // need this + + x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; + if (GridRank > 1) + y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j]; + if (GridRank > 2) + z = CellLeftEdge[2][k] + 0.5*CellWidth[2][k]; + + if (i >= GridStartIndex[0] && i <= GridEndIndex[0] && + j >= GridStartIndex[1] && j <= GridEndIndex[1] && + k >= GridStartIndex[2] && k <= GridEndIndex[2]) { + cindex = (i-GridStartIndex[0]) + ActiveDims[0] * + ((j-GridStartIndex[1]) + (k-GridStartIndex[2])*ActiveDims[1]); + if (density_field != NULL) + density = max(density_field[cindex], 1e-6); + else + density = 1.0; + if (HII_field != NULL) + HII_Fraction = HII_field[cindex]; + else + HII_Fraction = PhotonTestInitialFractionHII; + if (HeII_field != NULL) + HeII_Fraction = HeII_field[cindex]; + else + HeII_Fraction = PhotonTestInitialFractionHeII; + if (HeIII_field != NULL) + HeIII_Fraction = HeIII_field[cindex]; + else + HeIII_Fraction = PhotonTestInitialFractionHeIII; + if (Temperature_field != NULL) + temperature = temp1 = Temperature_field[cindex]; + else + temperature = temp1 = InitialTemperature; + } + + /* set density, total energy */ // check + BaryonField[0][n] = TestStarParticleDensity; + BaryonField[1][n] = TestStarParticleEnergy; + + /* set velocities */ + + for (dim = 0; dim < GridRank; dim++) + BaryonField[vel+dim][n] = Velocity[dim] + TestStarParticleVelocity[dim]; + + /* Set internal energy if necessary. */ //check + + if (DualEnergyFormalism) + BaryonField[2][n] = BaryonField[1][n]; + + /* If doing multi-species (HI, etc.), set these. */ + + if (MultiSpecies > 0) { + BaryonField[HIINum][n] = TestProblemData.HII_Fraction * + TestProblemData.HydrogenFractionByMass * TestStarParticleDensity; + + BaryonField[HeIINum][n] = TestProblemData.HeII_Fraction * + TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); + + BaryonField[HeIIINum][n] = TestProblemData.HeIII_Fraction * + TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); + + BaryonField[HeINum][n] = + (1.0 - TestProblemData.HydrogenFractionByMass)*TestStarParticleDensity - + BaryonField[HeIINum][n] - BaryonField[HeIIINum][n]; + } + + if (MultiSpecies > 1) { + BaryonField[HMNum][n] = TestProblemData.HM_Fraction * + BaryonField[HIINum][n]; + + BaryonField[H2INum][n] = TestProblemData.H2I_Fraction * + BaryonField[0][n] * TestProblemData.HydrogenFractionByMass; + + BaryonField[H2IINum][n] = TestProblemData.H2II_Fraction * 2.0 * + BaryonField[HIINum][n]; + + BaryonField[kdissH2INum][n] = 0.0; + BaryonField[kphHMNum][n] = 0.0; + BaryonField[kdissH2IINum][n] = 0.0; + } + + // HI density is calculated by subtracting off the various ionized fractions + // from the total + BaryonField[HINum][n] = TestProblemData.HydrogenFractionByMass*BaryonField[0][n] + - BaryonField[HIINum][n]; + + if (MultiSpecies > 1) + BaryonField[HINum][n] -= (BaryonField[HMNum][n] + BaryonField[H2IINum][n] + + BaryonField[H2INum][n]); + + BaryonField[DeNum][n] = BaryonField[HIINum][n] + + 0.25*BaryonField[HeIINum][n] + 0.5*BaryonField[HeIIINum][n]; + + if (MultiSpecies > 1) + BaryonField[DeNum][n] += 0.5*BaryonField[H2IINum][n] - + BaryonField[HMNum][n]; + + /* Set Deuterium species (assumed to be negligible). */ + + if (MultiSpecies > 2) { + BaryonField[DINum ][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HINum][n]; + BaryonField[DIINum][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HIINum][n]; + BaryonField[HDINum][n] = 0.75 * TestProblemData.DeuteriumToHydrogenRatio * BaryonField[H2INum][n]; + } + + /* Set energy (thermal and then total if necessary). */ + + if (MultiSpecies) { + mu_data = + 0.25*(BaryonField[HeINum][n] + BaryonField[HeIINum][n] + + BaryonField[HeIIINum][n] ) + + BaryonField[HINum][n] + BaryonField[HIINum][n] + + BaryonField[DeNum][n]; + if (MultiSpecies > 1) + mu_data += BaryonField[HMNum][n] + + 0.5*(BaryonField[H2INum][n] + BaryonField[H2IINum][n]); + mu_data = BaryonField[0][n] / mu_data; + } else + mu_data = mu; + + BaryonField[1][n] = temperature/TemperatureUnits/ + ((Gamma-1.0)*mu_data); + + } + + delete [] density_field; + delete [] HII_field; + delete [] HeII_field; + delete [] HeIII_field; + delete [] Temperature_field; + + return SUCCESS; +} + diff --git a/src/enzo/Grid_TestStarParticleInitializeGrid.C b/src/enzo/Grid_TestStarParticleInitializeGrid.C index a6900148c..d470d1654 100644 --- a/src/enzo/Grid_TestStarParticleInitializeGrid.C +++ b/src/enzo/Grid_TestStarParticleInitializeGrid.C @@ -14,7 +14,6 @@ #include #include -#include #include #include "ErrorExceptions.h" #include "macros_and_parameters.h" @@ -30,144 +29,26 @@ int GetUnits(float *DensityUnits, float *LengthUnits, float *TemperatureUnits, float *TimeUnits, float *VelocityUnits, double *MassUnits, FLOAT Time); -// Used to compute Bonner-Ebert density profile -double ph_BE(double r); -double ph_q(double r); -double ph_Ang(double a1, double a2, double R, double r); - -// Returns random velocity from Maxwellian distribution -double ph_Maxwellian(double c_tilda, double vel_unit, double mu, double gamma); - -static int PhotonTestParticleCount = 0; - int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, float *Initialdt, FLOAT TestStarParticleStarVelocity[], - FLOAT TestStarParticleStarPosition[], - float TestStarParticleDensity, - float TestStarParticleEnergy, - float TestStarParticleVelocity[], - float TestStarParticleBField[], - float InitialTemperature, - float PhotonTestInitialFractionHII, - float PhotonTestInitialFractionHeII, - float PhotonTestInitialFractionHeIII, - float PhotonTestInitialFractionHM, - float PhotonTestInitialFractionH2I, - float PhotonTestInitialFractionH2II) + FLOAT TestStarParticleStarPosition[]) { /* declarations */ float CentralMass = 1.0; - int dim, i, j, k, m, field, size, active_size, index, cindex; - - int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, - DINum, DIINum, HDINum, MetalNum, MetalIaNum, B1Num, B2Num, B3Num, PhiNum, CRNum; - - int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum, kdissH2IINum, kphHMNum, RPresNum1, RPresNum2, RPresNum3; - - int ExtraField[2]; - - float TestInitialdt = *Initialdt, *density_field = NULL, *HII_field = NULL, *HeII_field = NULL, - *HeIII_field = NULL, *Temperature_field = NULL; + int i, dim; + float TestInitialdt = *Initialdt; /* Return if this doesn't concern us. */ if (ProcessorNumber != MyProcessorNumber) return SUCCESS; - - /* create fields */ - - NumberOfBaryonFields = 0; - FieldType[NumberOfBaryonFields++] = Density; - FieldType[NumberOfBaryonFields++] = TotalEnergy; - if (DualEnergyFormalism) - FieldType[NumberOfBaryonFields++] = InternalEnergy; - int vel = NumberOfBaryonFields; - FieldType[NumberOfBaryonFields++] = Velocity1; - if (GridRank > 1 || HydroMethod > 2) - FieldType[NumberOfBaryonFields++] = Velocity2; - if (GridRank > 2 || HydroMethod > 2) - FieldType[NumberOfBaryonFields++] = Velocity3; - - if ( CRModel ) { - CRNum = NumberOfBaryonFields; - FieldType[NumberOfBaryonFields++] = CRDensity; - } - - if (WritePotential) - FieldType[NumberOfBaryonFields++] = GravPotential; - - - int colorfields = NumberOfBaryonFields; - - // Enzo's standard multispecies (primordial chemistry - H, D, He) - if (TestProblemData.MultiSpecies) { - FieldType[DeNum = NumberOfBaryonFields++] = ElectronDensity; - FieldType[HINum = NumberOfBaryonFields++] = HIDensity; - FieldType[HIINum = NumberOfBaryonFields++] = HIIDensity; - FieldType[HeINum = NumberOfBaryonFields++] = HeIDensity; - FieldType[HeIINum = NumberOfBaryonFields++] = HeIIDensity; - FieldType[HeIIINum = NumberOfBaryonFields++] = HeIIIDensity; - if (TestProblemData.MultiSpecies > 1) { - FieldType[HMNum = NumberOfBaryonFields++] = HMDensity; - FieldType[H2INum = NumberOfBaryonFields++] = H2IDensity; - FieldType[H2IINum = NumberOfBaryonFields++] = H2IIDensity; - } - if (TestProblemData.MultiSpecies > 2) { - FieldType[DINum = NumberOfBaryonFields++] = DIDensity; - FieldType[DIINum = NumberOfBaryonFields++] = DIIDensity; - FieldType[HDINum = NumberOfBaryonFields++] = HDIDensity; - } - } - - // Metal fields, including the standard 'metallicity' as well - // as two extra fields - if (TestProblemData.UseMetallicityField) { - FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity; - - if (StarMakerTypeIaSNe) - FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity; - - if(TestProblemData.MultiMetals){ - FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0; - FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1; - } - } - - if (RadiativeTransfer && (MultiSpecies < 1)) { - ENZO_FAIL("Grid_TestStarParticleInitialize: Radiative Transfer but not MultiSpecies set"); - } - - // Allocate fields for photo ionization and heating rates - if (RadiativeTransfer) - if (MultiSpecies) { - FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; - FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; - if (RadiativeTransferHydrogenOnly == FALSE) { - FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; - FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; - } - if (MultiSpecies > 1) { - FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; - FieldType[kdissH2IINum = NumberOfBaryonFields++] = kdissH2II; - FieldType[kphHMNum = NumberOfBaryonFields++] = kphHM; - } - } - - if (RadiationPressure && RadiativeTransfer) { - FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0; - FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1; - FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2; - } - - NumberOfPhotonPackages = 0; - PhotonPackages-> NextPackage= NULL; /* Get Units. */ float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, - VelocityUnits = 1, TimeUnits = 1, CriticalDensity = 1, BoxLength = 1, mu = 0.6, mu_data; + VelocityUnits = 1, TimeUnits = 1; double MassUnits = 1; if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, @@ -188,22 +69,12 @@ int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, this->AllocateNewParticles(NumberOfParticles); printf("Allocated %d particles\n", NumberOfParticles); - this->AllocateGrids(); - - /* Initialize radiation fields */ // need this - - if (this->InitializeRadiativeTransferFields() == FAIL) { - ENZO_FAIL("\nError in InitializeRadiativeTransferFields.\n"); - } - /* Set particle IDs and types */ for (i = 0; i < NumberOfParticles; i++) { ParticleNumber[i] = i; - ParticleType[i] = -PopIII; - ParticleAttribute[0][i] = Time + 1e-7; + ParticleType[i] = PARTICLE_TYPE_STAR; } - ParticleMass[0] = CentralMass; /* Set central particle. */ for (dim = 0; dim < GridRank; dim++) { @@ -211,6 +82,8 @@ int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, (DomainLeftEdge[dim]+DomainRightEdge[dim]) + 0.5*CellWidth[0][0]; ParticleVelocity[dim][0] = TestStarParticleStarVelocity[dim]*1e5*TimeUnits/LengthUnits; } + ParticleMass[0] = CentralMass; + ParticleAttribute[0][0] = Time+1e-7; if (STARFEED_METHOD(UNIGRID_STAR)) ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; if (STARFEED_METHOD(MOM_STAR)) @@ -222,171 +95,6 @@ int grid::TestStarParticleInitializeGrid(float TestStarParticleStarMass, ParticleAttribute[2][0] = 0.0; // Metal fraction ParticleAttribute[3][0] = 0.0; // metalfSNIa - if (this->FindNewStarParticles(GridLevel) == FAIL) { - ENZO_FAIL("Error in grid::FindNewStarParticles."); - } - - /* Reset particle type to be positive*/ - for (i = 0; i < NumberOfParticles; i++) { - ParticleNumber[i] = i; - ParticleType[i] = PopIII; - } - Star *cstar; - for (cstar = Stars; cstar; cstar = cstar->NextStar) - cstar->type = PopIII; - - - /* Set up the baryon field. */ // need thid - /* compute size of fields */ - active_size = 1; - int ActiveDims[MAX_DIMENSION]; - for (dim = 0; dim < GridRank; dim++) { - ActiveDims[dim] = GridEndIndex[dim] - GridStartIndex[dim] + 1; - active_size *= ActiveDims[dim]; - } - - /* Loop over the mesh. */ - float density, dens1, Velocity[MAX_DIMENSION], - temperature, temp1, sigma, sigma1, colour, outer_radius; - float HII_Fraction, HeII_Fraction, HeIII_Fraction, H2I_Fraction; - FLOAT r, x, y = 0, z = 0; - int n = 0; - - for (k = 0; k < GridDimension[2]; k++) - for (j = 0; j < GridDimension[1]; j++) - for (i = 0; i < GridDimension[0]; i++, n++) { - - /* Compute position */ // need this - - x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; - if (GridRank > 1) - y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j]; - if (GridRank > 2) - z = CellLeftEdge[2][k] + 0.5*CellWidth[2][k]; - - if (i >= GridStartIndex[0] && i <= GridEndIndex[0] && - j >= GridStartIndex[1] && j <= GridEndIndex[1] && - k >= GridStartIndex[2] && k <= GridEndIndex[2]) { - cindex = (i-GridStartIndex[0]) + ActiveDims[0] * - ((j-GridStartIndex[1]) + (k-GridStartIndex[2])*ActiveDims[1]); - if (density_field != NULL) - density = max(density_field[cindex], 1e-6); - else - density = 1.0; - if (HII_field != NULL) - HII_Fraction = HII_field[cindex]; - else - HII_Fraction = PhotonTestInitialFractionHII; - if (HeII_field != NULL) - HeII_Fraction = HeII_field[cindex]; - else - HeII_Fraction = PhotonTestInitialFractionHeII; - if (HeIII_field != NULL) - HeIII_Fraction = HeIII_field[cindex]; - else - HeIII_Fraction = PhotonTestInitialFractionHeIII; - if (Temperature_field != NULL) - temperature = temp1 = Temperature_field[cindex]; - else - temperature = temp1 = InitialTemperature; - } - - /* set density, total energy */ // check - BaryonField[0][n] = TestStarParticleDensity; - BaryonField[1][n] = TestStarParticleEnergy; - - /* set velocities */ - - for (dim = 0; dim < GridRank; dim++) - BaryonField[vel+dim][n] = Velocity[dim] + TestStarParticleVelocity[dim]; - - /* Set internal energy if necessary. */ //check - - if (DualEnergyFormalism) - BaryonField[2][n] = BaryonField[1][n]; - - /* If doing multi-species (HI, etc.), set these. */ - - if (MultiSpecies > 0) { - BaryonField[HIINum][n] = TestProblemData.HII_Fraction * - TestProblemData.HydrogenFractionByMass * TestStarParticleDensity; - - BaryonField[HeIINum][n] = TestProblemData.HeII_Fraction * - TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); - - BaryonField[HeIIINum][n] = TestProblemData.HeIII_Fraction * - TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); - - BaryonField[HeINum][n] = - (1.0 - TestProblemData.HydrogenFractionByMass)*TestStarParticleDensity - - BaryonField[HeIINum][n] - BaryonField[HeIIINum][n]; - } - - if (MultiSpecies > 1) { - BaryonField[HMNum][n] = TestProblemData.HM_Fraction * - BaryonField[HIINum][n]; - - BaryonField[H2INum][n] = TestProblemData.H2I_Fraction * - BaryonField[0][n] * TestProblemData.HydrogenFractionByMass; - - BaryonField[H2IINum][n] = TestProblemData.H2II_Fraction * 2.0 * - BaryonField[HIINum][n]; - - BaryonField[kdissH2INum][n] = 0.0; - BaryonField[kphHMNum][n] = 0.0; - BaryonField[kdissH2IINum][n] = 0.0; - } - - // HI density is calculated by subtracting off the various ionized fractions - // from the total - BaryonField[HINum][n] = TestProblemData.HydrogenFractionByMass*BaryonField[0][n] - - BaryonField[HIINum][n]; - - if (MultiSpecies > 1) - BaryonField[HINum][n] -= (BaryonField[HMNum][n] + BaryonField[H2IINum][n] - + BaryonField[H2INum][n]); - - BaryonField[DeNum][n] = BaryonField[HIINum][n] + - 0.25*BaryonField[HeIINum][n] + 0.5*BaryonField[HeIIINum][n]; - - if (MultiSpecies > 1) - BaryonField[DeNum][n] += 0.5*BaryonField[H2IINum][n] - - BaryonField[HMNum][n]; - - /* Set Deuterium species (assumed to be negligible). */ - - if (MultiSpecies > 2) { - BaryonField[DINum ][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HINum][n]; - BaryonField[DIINum][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HIINum][n]; - BaryonField[HDINum][n] = 0.75 * TestProblemData.DeuteriumToHydrogenRatio * BaryonField[H2INum][n]; - } - - /* Set energy (thermal and then total if necessary). */ - - if (MultiSpecies) { - mu_data = - 0.25*(BaryonField[HeINum][n] + BaryonField[HeIINum][n] + - BaryonField[HeIIINum][n] ) + - BaryonField[HINum][n] + BaryonField[HIINum][n] + - BaryonField[DeNum][n]; - if (MultiSpecies > 1) - mu_data += BaryonField[HMNum][n] + - 0.5*(BaryonField[H2INum][n] + BaryonField[H2IINum][n]); - mu_data = BaryonField[0][n] / mu_data; - } else - mu_data = mu; - - BaryonField[1][n] = temperature/TemperatureUnits/ - ((Gamma-1.0)*mu_data); - - } - - delete [] density_field; - delete [] HII_field; - delete [] HeII_field; - delete [] HeIII_field; - delete [] Temperature_field; - return SUCCESS; } diff --git a/src/enzo/InitializeNew.C b/src/enzo/InitializeNew.C index e47f8491d..c6fb9b267 100644 --- a/src/enzo/InitializeNew.C +++ b/src/enzo/InitializeNew.C @@ -83,6 +83,8 @@ int StratifiedMediumExplosionInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr TopGridData &MetaData); int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData, float *Initialdt); +int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, + TopGridData &MetaData, float *Initialdt); int KHInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData); int NohInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, @@ -748,7 +750,11 @@ int InitializeNew(char *filename, HierarchyEntry &TopGrid, #endif // Insert new problem intializer here... - + + // 251) Test a star particle explosion + if (ProblemType == 251) + ret = TestRadiatingStarParticleInitialize(fptr, Outfptr, TopGrid, MetaData, + Initialdt); if (ret == INT_UNDEFINED) { ENZO_VFAIL("Problem Type %"ISYM" undefined.\n", ProblemType) diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects index ddcfd9244..271c12986 100644 --- a/src/enzo/Make.config.objects +++ b/src/enzo/Make.config.objects @@ -650,6 +650,7 @@ OBJS_CONFIG_LIB = \ Grid_TestGravitySphereInitializeGrid.o \ Grid_TestOrbitInitializeGrid.o \ Grid_TestStarParticleInitializeGrid.o \ + Grid_TestRadiatingStarParticleInitializeGrid.o \ Grid_TracerParticleCreateParticles.o \ Grid_TracerParticleOutputData.o \ Grid_TracerParticleSetVelocity.o \ @@ -938,6 +939,7 @@ OBJS_CONFIG_LIB = \ TestGravitySphereInitialize.o \ TestOrbitInitialize.o \ TestStarParticleInitialize.o \ + TestRadiatingStarParticleInitialize.o \ TracerParticleCreation.o \ TransposeRegionOverlap.o \ tscint1d.o \ diff --git a/src/enzo/Make.mach.linux-gnu b/src/enzo/Make.mach.linux-gnu index f9ec4422e..d3a8bc446 100644 --- a/src/enzo/Make.mach.linux-gnu +++ b/src/enzo/Make.mach.linux-gnu @@ -114,7 +114,7 @@ LOCAL_LIBS_MPI = # MPI libraries LOCAL_LIBS_HDF5 = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz LOCAL_LIBS_HYPRE = -L$(LOCAL_HYPRE_INSTALL)/lib -lHYPRE LOCAL_LIBS_PAPI = # PAPI libraries -LOCAL_LIBS_MACH = -lgfortran # Machine-dependent libraries +LOCAL_LIBS_MACH = -lgfortran -lquadmath # Machine-dependent libraries LOCAL_LIBS_GRACKLE = -L/home/dskinner6/software/grackle/lib -lgrackle MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) diff --git a/src/enzo/TestRadiatingStarParticleInitialize.C b/src/enzo/TestRadiatingStarParticleInitialize.C new file mode 100644 index 000000000..edcc6bc7b --- /dev/null +++ b/src/enzo/TestRadiatingStarParticleInitialize.C @@ -0,0 +1,247 @@ +/*********************************************************************** +/ +/ INITIALIZE A RADIATING STAR PARTICLE TEST +/ +/ written by: Greg Bryan +/ date: June, 2012 +/ modified1: Danielle Skinner +/ date: July, 2022 +/ +/ PURPOSE: +/ Initialize a radiating star particle test in a uniform medium. +/ Based on the Single Star Particle Test and Photon Test. +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + +// This routine intializes a new simulation based on the parameter file. +// + +#include +#include +#include +#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" + +static float PhotonTestInitialFractionHII = 1.2e-5; +static float PhotonTestInitialFractionHeII = 1.0e-14; +static float PhotonTestInitialFractionHeIII = 1.0e-17; +static float PhotonTestInitialFractionHM = 2.0e-9; +static float PhotonTestInitialFractionH2I = 2.0e-20; +static float PhotonTestInitialFractionH2II = 3.0e-14; + +int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, + TopGridData &MetaData,float *Initialdt) +{ + const char *DensName = "Density"; /* make const */ + const char *TEName = "TotalEnergy"; + const char *GEName = "GasEnergy"; + const char *Vel1Name = "x-velocity"; + const char *Vel2Name = "y-velocity"; + const char *Vel3Name = "z-velocity"; + const char *MetalName = "Metal_Density"; + const char *ElectronName = "Electron_Density"; + const char *HIName = "HI_Density"; + const char *HIIName = "HII_Density"; + const char *HeIName = "HeI_Density"; + const char *HeIIName = "HeII_Density"; + const char *HeIIIName = "HeIII_Density"; + const char *HMName = "HM_Density"; + const char *H2IName = "H2I_Density"; + const char *H2IIName = "H2II_Density"; + const char *DIName = "DI_Density"; + const char *DIIName = "DII_Density"; + const char *HDIName = "HDI_Density"; + const char *kphHIName = "HI_kph"; + const char *gammaName = "PhotoGamma"; + const char *kphHeIName = "HeI_kph"; + const char *kphHeIIName = "HeII_kph"; + const char *kphHMName = "HM_kph"; + const char *kdissH2IIName = "H2II_kdiss"; + const char *kdissH2IName = "H2I_kdiss"; + const char *RadAccel1Name = "x-RadPressure"; + const char *RadAccel2Name = "y-RadPressure"; + const char *RadAccel3Name = "z-RadPressure"; + + + /* declarations */ + + char line[MAX_LINE_LENGTH]; + char *dummy = new char[MAX_LINE_LENGTH]; + int dim, ret, i, source; + dummy[0] = 0; + + /* Error check. */ + + + /* set default parameters */ + + float TestStarParticleDensity = 1.0; + float TestStarParticleEnergy = 1.0; + float TestStarParticleVelocity[3] = {0.0, 0.0, 0.0}; + FLOAT TestStarParticleStarVelocity[3] = {0.0, 0.0, 0.0}; + FLOAT TestStarParticleStarPosition[3] = {0.5, 0.5, 0.5}; + float TestStarParticleBField[3] = {0.0, 0.0, 0.0}; + float TestStarParticleStarMass = 100.0; + int TestProblemUseMetallicityField = 1; + float TestProblemInitialMetallicityFraction = 2e-3; // 0.1 Zsun + float PhotonTestInitialTemperature = 1000; + + + + + TestProblemData.MultiSpecies = MultiSpecies; + TestProblemData.UseMetallicityField = TestProblemUseMetallicityField; + TestProblemData.MetallicityField_Fraction = TestProblemInitialMetallicityFraction; + + /* read input from file */ + + + while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) { + + ret = 0; + + /* read parameters */ + + ret += sscanf(line, "TestStarParticleDensity = %"FSYM, + &TestStarParticleDensity); + ret += sscanf(line, "TestStarParticleEnergy = %"FSYM, + &TestStarParticleEnergy); + ret += sscanf(line, "TestStarParticleStarMass = %"FSYM, + &TestStarParticleStarMass); + ret += sscanf(line,"TestStarParticleStarVelocity = %"PSYM" %"PSYM" %"PSYM, + &TestStarParticleStarVelocity[0], + &TestStarParticleStarVelocity[1], + &TestStarParticleStarVelocity[2]); + ret += sscanf(line,"TestStarParticleStarPosition = %"PSYM" %"PSYM" %"PSYM, + &TestStarParticleStarPosition[0], + &TestStarParticleStarPosition[1], + &TestStarParticleStarPosition[2]); + + + ret += sscanf(line, "TestProblemUseMetallicityField = %"ISYM, &TestProblemData.UseMetallicityField); + ret += sscanf(line, "TestProblemInitialMetallicityFraction = %"FSYM, &TestProblemData.MetallicityField_Fraction); + + ret += sscanf(line, "TestProblemInitialHIFraction = %"FSYM, &TestProblemData.HI_Fraction); + ret += sscanf(line, "TestProblemInitialHIIFraction = %"FSYM, &TestProblemData.HII_Fraction); + ret += sscanf(line, "TestProblemInitialHeIFraction = %"FSYM, &TestProblemData.HeI_Fraction); + ret += sscanf(line, "TestProblemInitialHeIIFraction = %"FSYM, &TestProblemData.HeII_Fraction); + ret += sscanf(line, "TestProblemInitialHeIIIIFraction = %"FSYM, &TestProblemData.HeIII_Fraction); + ret += sscanf(line, "TestProblemHydrogenFractionByMass = %"FSYM, &TestProblemData.HydrogenFractionByMass); + + + /* if the line is suspicious, issue a warning */ + + if (ret == 0 && strstr(line, "=") && strstr(line, "TestStarParticle") + && line[0] != '#') + fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line); + + } // end input from parameter file + + /* set up uniform grid as of before explosion */ + + if (TopGrid.GridData-> + TestRadiatingStarParticleInitializeGrid(TestStarParticleStarMass, + Initialdt, + TestStarParticleStarVelocity, + TestStarParticleStarPosition, + TestStarParticleDensity, + TestStarParticleEnergy, + TestStarParticleVelocity, + TestStarParticleBField, + PhotonTestInitialTemperature, + PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII, + PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM, + PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II) == FAIL) + ENZO_FAIL("Error in TestRadiatingStarParticleInitializeGrid.\n"); + + /* set up field names and units */ + + int count = 0; + DataLabel[count++] = (char*) DensName; + DataLabel[count++] = (char*) TEName; + if (DualEnergyFormalism) + DataLabel[count++] = (char*) GEName; + DataLabel[count++] = (char*) Vel1Name; + DataLabel[count++] = (char*) Vel2Name; + DataLabel[count++] = (char*) Vel3Name; + if (TestProblemData.MultiSpecies){ + DataLabel[count++] = (char*) ElectronName; + DataLabel[count++] = (char*) HIName; + DataLabel[count++] = (char*) HIIName; + DataLabel[count++] = (char*) HeIName; + DataLabel[count++] = (char*) HeIIName; + DataLabel[count++] = (char*) HeIIIName; + if (MultiSpecies > 1) { + DataLabel[count++] = (char*) HMName; + DataLabel[count++] = (char*) H2IName; + DataLabel[count++] = (char*) H2IIName; + } + if (MultiSpecies > 2) { + DataLabel[count++] = (char*) DIName; + DataLabel[count++] = (char*) DIIName; + DataLabel[count++] = (char*) HDIName; + } + } // if Multispecies + if (TestProblemData.UseMetallicityField) + DataLabel[count++] = (char*) MetalName; + + if (RadiativeTransfer) + if (MultiSpecies) { + DataLabel[count++] = (char*) kphHIName; + DataLabel[count++] = (char*) gammaName; + DataLabel[count++] = (char*) kphHeIName; + DataLabel[count++] = (char*) kphHeIIName; + if (MultiSpecies > 1) { + DataLabel[count++]= (char*) kdissH2IName; + DataLabel[count++]= (char*) kdissH2IIName; + DataLabel[count++]= (char*) kphHMName; + } + } // if RadiativeTransfer + + if (RadiationPressure) { + DataLabel[count++] = (char*) RadAccel1Name; + DataLabel[count++] = (char*) RadAccel2Name; + DataLabel[count++] = (char*) RadAccel3Name; + } + + + int j; + for(j=0; j < count; j++) + DataUnits[j] = NULL; + + + /* Write parameters to parameter output file */ + + if (MyProcessorNumber == ROOT_PROCESSOR) { + + fprintf(Outfptr, "TestStarParticleDensity = %"FSYM"\n", + TestStarParticleDensity); + fprintf(Outfptr, "TestStarParticleEnergy = %"FSYM"\n", + TestStarParticleEnergy); + fprintf(Outfptr, "MetallicityField_Fraction = %"FSYM"\n", + TestProblemData.MetallicityField_Fraction); + } + + fprintf(stderr, "TestStarParticleDensity = %"FSYM"\n", + TestStarParticleDensity); + fprintf(stderr, "TestStarParticleEnergy = %"FSYM"\n", + TestStarParticleEnergy); + fprintf(stderr, "MetallicityField_Fraction = %"FSYM"\n", + TestProblemData.MetallicityField_Fraction); + + + + return SUCCESS; + +} + diff --git a/src/enzo/TestStarParticleInitialize.C b/src/enzo/TestStarParticleInitialize.C index 43c1b8d1e..d84e0f7c6 100644 --- a/src/enzo/TestStarParticleInitialize.C +++ b/src/enzo/TestStarParticleInitialize.C @@ -30,53 +30,28 @@ #include "Hierarchy.h" #include "TopGridData.h" -static float PhotonTestInitialFractionHII = 1.2e-5; -static float PhotonTestInitialFractionHeII = 1.0e-14; -static float PhotonTestInitialFractionHeIII = 1.0e-17; -static float PhotonTestInitialFractionHM = 2.0e-9; -static float PhotonTestInitialFractionH2I = 2.0e-20; -static float PhotonTestInitialFractionH2II = 3.0e-14; - int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData,float *Initialdt) { - const char *DensName = "Density"; /* make const */ - const char *TEName = "TotalEnergy"; - const char *GEName = "GasEnergy"; - const char *Vel1Name = "x-velocity"; - const char *Vel2Name = "y-velocity"; - const char *Vel3Name = "z-velocity"; - const char *MetalName = "Metal_Density"; - const char *ElectronName = "Electron_Density"; - const char *HIName = "HI_Density"; - const char *HIIName = "HII_Density"; - const char *HeIName = "HeI_Density"; - const char *HeIIName = "HeII_Density"; - const char *HeIIIName = "HeIII_Density"; - const char *HMName = "HM_Density"; - const char *H2IName = "H2I_Density"; - const char *H2IIName = "H2II_Density"; - const char *DIName = "DI_Density"; - const char *DIIName = "DII_Density"; - const char *HDIName = "HDI_Density"; - const char *kphHIName = "HI_kph"; - const char *gammaName = "PhotoGamma"; - const char *kphHeIName = "HeI_kph"; - const char *kphHeIIName = "HeII_kph"; - const char *kphHMName = "HM_kph"; - const char *kdissH2IIName = "H2II_kdiss"; - const char *kdissH2IName = "H2I_kdiss"; - const char *RadAccel1Name = "x-RadPressure"; - const char *RadAccel2Name = "y-RadPressure"; - const char *RadAccel3Name = "z-RadPressure"; + char *DensName = "Density"; + char *TEName = "TotalEnergy"; + char *GEName = "GasEnergy"; + char *Vel1Name = "x-velocity"; + char *Vel2Name = "y-velocity"; + char *Vel3Name = "z-velocity"; + char *MetalName = "Metal_Density"; + char *ElectronName = "Electron_Density"; + char *HIName = "HI_Density"; + char *HIIName = "HII_Density"; + char *HeIName = "HeI_Density"; + char *HeIIName = "HeII_Density"; + char *HeIIIName = "HeIII_Density"; /* declarations */ char line[MAX_LINE_LENGTH]; - char *dummy = new char[MAX_LINE_LENGTH]; - int dim, ret, i, source; - dummy[0] = 0; + int dim, ret; /* Error check. */ @@ -92,7 +67,6 @@ int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGri float TestStarParticleStarMass = 100.0; int TestProblemUseMetallicityField = 1; float TestProblemInitialMetallicityFraction = 2e-3; // 0.1 Zsun - float PhotonTestInitialTemperature = 1000; @@ -147,70 +121,40 @@ int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGri /* set up uniform grid as of before explosion */ + if (TopGrid.GridData->InitializeUniformGrid(TestStarParticleDensity, + TestStarParticleEnergy, + TestStarParticleEnergy, + TestStarParticleVelocity, + TestStarParticleBField) == FAIL) + ENZO_FAIL("Error in InitializeUniformGrid."); + if (TopGrid.GridData-> TestStarParticleInitializeGrid(TestStarParticleStarMass, Initialdt, TestStarParticleStarVelocity, - TestStarParticleStarPosition, - TestStarParticleDensity, - TestStarParticleEnergy, - TestStarParticleVelocity, - TestStarParticleBField, - PhotonTestInitialTemperature, - PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII, - PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM, - PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II) == FAIL) + TestStarParticleStarPosition) == FAIL) ENZO_FAIL("Error in TestStarParticleInitializeGrid.\n"); /* set up field names and units */ int count = 0; - DataLabel[count++] = (char*) DensName; - DataLabel[count++] = (char*) TEName; + DataLabel[count++] = DensName; + DataLabel[count++] = TEName; if (DualEnergyFormalism) - DataLabel[count++] = (char*) GEName; - DataLabel[count++] = (char*) Vel1Name; - DataLabel[count++] = (char*) Vel2Name; - DataLabel[count++] = (char*) Vel3Name; + DataLabel[count++] = GEName; + DataLabel[count++] = Vel1Name; + DataLabel[count++] = Vel2Name; + DataLabel[count++] = Vel3Name; if (TestProblemData.MultiSpecies){ - DataLabel[count++] = (char*) ElectronName; - DataLabel[count++] = (char*) HIName; - DataLabel[count++] = (char*) HIIName; - DataLabel[count++] = (char*) HeIName; - DataLabel[count++] = (char*) HeIIName; - DataLabel[count++] = (char*) HeIIIName; - if (MultiSpecies > 1) { - DataLabel[count++] = (char*) HMName; - DataLabel[count++] = (char*) H2IName; - DataLabel[count++] = (char*) H2IIName; - } - if (MultiSpecies > 2) { - DataLabel[count++] = (char*) DIName; - DataLabel[count++] = (char*) DIIName; - DataLabel[count++] = (char*) HDIName; - } - } // if Multispecies - if (TestProblemData.UseMetallicityField) - DataLabel[count++] = (char*) MetalName; - - if (RadiativeTransfer) - if (MultiSpecies) { - DataLabel[count++] = (char*) kphHIName; - DataLabel[count++] = (char*) gammaName; - DataLabel[count++] = (char*) kphHeIName; - DataLabel[count++] = (char*) kphHeIIName; - if (MultiSpecies > 1) { - DataLabel[count++]= (char*) kdissH2IName; - DataLabel[count++]= (char*) kdissH2IIName; - DataLabel[count++]= (char*) kphHMName; - } - } // if RadiativeTransfer - - if (RadiationPressure) { - DataLabel[count++] = (char*) RadAccel1Name; - DataLabel[count++] = (char*) RadAccel2Name; - DataLabel[count++] = (char*) RadAccel3Name; + DataLabel[count++] = ElectronName; + DataLabel[count++] = HIName; + DataLabel[count++] = HIIName; + DataLabel[count++] = HeIName; + DataLabel[count++] = HeIIName; + DataLabel[count++] = HeIIIName; } + if (TestProblemData.UseMetallicityField) + DataLabel[count++] = MetalName; int j; From 28dec0bae86528288f3126aabf0ca280eb8adf48 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Thu, 21 Jul 2022 11:35:57 -0400 Subject: [PATCH 04/10] Adding param file for new test. --- .../TestRadiatingStarParticleSingle.enzo | 92 +++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo diff --git a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo new file mode 100644 index 000000000..79543e182 --- /dev/null +++ b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo @@ -0,0 +1,92 @@ +# +# AMR PROBLEM DEFINITION FILE: TestRadiatingStarParticle +# +# define problem +# +ProblemType = 251 // TestRadiatingStarParaticle +TopGridRank = 3 +TopGridDimensions = 50 50 50 +StarParticleFeedback = 16384 // Population III Stars +TestStarParticleDensity = 1.0 +TestStarParticleEnergy = 6.55410093179e-06 +TestProblemInitialMetallicityFraction = 2e-3 +TestStarParticleStarMass = 20.0 +StarEnergyToThermalFeedback = 5.59e-6 +StarMassEjectionFraction = 0.25 +StarFeedbackKineticFraction = 0.3 +StarMakerExplosionDelayTime = 0.0 +OutputTemperature = 1 +OutputCoolingTime = 1 +Initialdt = 0.0010017410642 +# +# Units +# +LengthUnits = 1.2344e+21 +TimeUnits = 3.15e13 // 10 Myr +DensityUnits = 1.0e-24 // 1 part/cc +# +# set I/O and stop/start parameters +# +StopTime = 1 +dtDataDump = 0.1 +# +# set Hydro parameters +# +HydroMethod = 0 +DualEnergyFormalism = 1 +CourantSafetyNumber = 0.4 +PPMDiffusionParameter = 0 // diffusion off +PPMFlatteningParameter = 0 // flattening on +PPMSteepeningParameter = 0 // steepening on +Gamma = 1.6666667 +# +# set grid refinement parameters +# +StaticHierarchy = 1 // static hierarchy +RefineBy = 2 // refinement factor +# +# Radiation +RadiativeTransferRaysPerCell = 5.100000 +RadiativeTransferInitialHEALPixLevel = 1 +RadiativeTransferHydrogenOnly = 0 +RadiativeTransferOpticallyThinH2 = 1 +RadiativeTransferPeriodicBoundary = 0 +RadiativeTransferAdaptiveTimestep = 1 +RadiativeTransferRadiationPressure = 1 +RadiativeTransferPhotonMergeRadius = 3.0 +RadiativeTransferSourceClustering = 1 + +# +RadiativeTransfer = 1 + +MaximumRefinementLevel = 0 +CellFlaggingMethod = 0 +# +# Cooling +MultiSpecies = 1 +RadiativeCooling = 0 +#MultiSpecies = 1 +#MetalCooling = 3 +#CloudyCoolingGridFile = solar_2008_3D_metals.h5 +#CMBTemperatureFloor = 0 +#TestProblemMultiSpecies = 1 +#TestProblemUseMetallicityField = 1 +#TestProblemHydrogenFractionByMass = 0.75 +#TestProblemInitialMetallicityFraction = 2e-3 +#TestProblemInitialHIFraction = 0.998 +#TestProblemInitialHeIFraction = 1.0 +#TestProblemInitialHeIIFraction = 1.0e-20 +#TestProblemInitialHeIIIIFraction = 1.0e-20 +# +# set some misc global parameters +# + +MixSpeciesAndColors = 1 +PopIIISupernovaUseColour = 1 + + +# +# Rotating Pop III Stars +# + +PopIII_Rotating = 1 \ No newline at end of file From e7b547c6e0cd8f0ae6180c1955ab626970e90c98 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Mon, 25 Jul 2022 15:05:10 -0400 Subject: [PATCH 05/10] Making requested changes. --- doc/manual/source/parameters/index.rst | 6 ++ .../TestRadiatingStarParticleSingle.enzo | 2 +- .../TestRadiatingStarParticleSingle.enzotest | 2 + src/enzo/Grid.h | 12 +-- ..._TestRadiatingStarParticleInitializeGrid.C | 22 +++--- src/enzo/Make.mach.linux-gnu | 21 +++-- src/enzo/ReadParameterFile.C | 2 +- src/enzo/SetDefaultGlobalValues.C | 2 +- src/enzo/Star_ComputePhotonRates.C | 77 +++++++++++-------- .../TestRadiatingStarParticleInitialize.C | 24 +++--- src/enzo/WriteParameterFile.C | 2 +- 11 files changed, 95 insertions(+), 77 deletions(-) create mode 100644 run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest diff --git a/doc/manual/source/parameters/index.rst b/doc/manual/source/parameters/index.rst index a0aedf421..cab6975d9 100644 --- a/doc/manual/source/parameters/index.rst +++ b/doc/manual/source/parameters/index.rst @@ -2338,6 +2338,12 @@ The parameters below are considered in ``StarParticleCreation`` method 3. Above this density, a Pop III "color" particle forms, and it will populate the surrounding region with a color field. Units: mean density. Default: 1e6 ``PopIIIColorMass`` (external) A Pop III "color" particle will populate the surrounding region with a mass of PopIIIColorMass. Units: solar masses. Default: 1e6 +``PopIIIRotating`` (external) + Updates ionizing photon rates to the rotational and non-rotational rates given by Murphy et al. (2021). Default: 0 + + 0 - Model is off, Schaerer (2002) photon rates are used. + 1 - Rotating model is on + 2 - Non-rotating model is on .. _radiative_star_cluster_formation_parameters: diff --git a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo index 79543e182..35aef4ecb 100644 --- a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo +++ b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo @@ -6,7 +6,7 @@ ProblemType = 251 // TestRadiatingStarParaticle TopGridRank = 3 TopGridDimensions = 50 50 50 -StarParticleFeedback = 16384 // Population III Stars +StarParticleFeedback = 8 // Population III Stars TestStarParticleDensity = 1.0 TestStarParticleEnergy = 6.55410093179e-06 TestProblemInitialMetallicityFraction = 2e-3 diff --git a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest new file mode 100644 index 000000000..7867f0e89 --- /dev/null +++ b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest @@ -0,0 +1,2 @@ +name = 'TestRadiatingStarParticleSingle' +radiation = 'ray' \ No newline at end of file diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 4c9592172..185ba4c05 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -2101,12 +2101,12 @@ int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], float UniformVelocity[], float UniformBField[], float InitialTemperature, - float PhotonTestInitialFractionHII, - float PhotonTestInitialFractionHeII, - float PhotonTestInitialFractionHeIII, - float PhotonTestInitialFractionHM, - float PhotonTestInitialFractionH2I, - float PhotonTestInitialFractionH2II); + float TestStarParticleInitialFractionHII, + float TestStarParticleInitialFractionHeII, + float TestStarParticleInitialFractionHeIII, + float TestStarParticleInitialFractionHM, + float TestStarParticleInitialFractionH2I, + float TestStarParticleInitialFractionH2II); /* Gravity Test: initialize grid. */ diff --git a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C index c480da812..c879983f8 100644 --- a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C +++ b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C @@ -1,6 +1,6 @@ /*********************************************************************** / -/ GRID CLASS (INITIALIZE THE GRID FOR A RADIATING STAR PARTICLE TEST) +/ GRID CLASS (INITIALIZE THE GRID FOR A RADIATING POP III STAR PARTICLE TEST) / / written by: Greg Bryan / date: June, 2012 @@ -39,7 +39,7 @@ double ph_Ang(double a1, double a2, double R, double r); // Returns random velocity from Maxwellian distribution double ph_Maxwellian(double c_tilda, double vel_unit, double mu, double gamma); -static int PhotonTestParticleCount = 0; +static int TestStarParticleParticleCount = 0; int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass, float *Initialdt, @@ -50,12 +50,12 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass float TestStarParticleVelocity[], float TestStarParticleBField[], float InitialTemperature, - float PhotonTestInitialFractionHII, - float PhotonTestInitialFractionHeII, - float PhotonTestInitialFractionHeIII, - float PhotonTestInitialFractionHM, - float PhotonTestInitialFractionH2I, - float PhotonTestInitialFractionH2II) + float TestStarParticleInitialFractionHII, + float TestStarParticleInitialFractionHeII, + float TestStarParticleInitialFractionHeIII, + float TestStarParticleInitialFractionHM, + float TestStarParticleInitialFractionH2I, + float TestStarParticleInitialFractionH2II) { /* declarations */ @@ -277,15 +277,15 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass if (HII_field != NULL) HII_Fraction = HII_field[cindex]; else - HII_Fraction = PhotonTestInitialFractionHII; + HII_Fraction = TestStarParticleInitialFractionHII; if (HeII_field != NULL) HeII_Fraction = HeII_field[cindex]; else - HeII_Fraction = PhotonTestInitialFractionHeII; + HeII_Fraction = TestStarParticleInitialFractionHeII; if (HeIII_field != NULL) HeIII_Fraction = HeIII_field[cindex]; else - HeIII_Fraction = PhotonTestInitialFractionHeIII; + HeIII_Fraction = TestStarParticleInitialFractionHeIII; if (Temperature_field != NULL) temperature = temp1 = Temperature_field[cindex]; else diff --git a/src/enzo/Make.mach.linux-gnu b/src/enzo/Make.mach.linux-gnu index d3a8bc446..1998aac0b 100644 --- a/src/enzo/Make.mach.linux-gnu +++ b/src/enzo/Make.mach.linux-gnu @@ -36,8 +36,8 @@ MACH_FILE = Make.mach.linux-gnu # Install paths (local variables) #----------------------------------------------------------------------- -LOCAL_HDF5_INSTALL = /home/dskinner6/software/yt-conda/pkgs/hdf5-1.10.2-hc401514_3# mandatory -LOCAL_GRACKLE_INSTALL = /home/dskinner6/software/grackle # optional +LOCAL_HDF5_INSTALL = /PATH/TO/SERIAL-HDF5/INSTALL # mandatory +LOCAL_GRACKLE_INSTALL = /PATH/TO/GRACKLE/INSTALL # optional LOCAL_HYPRE_INSTALL = /PATH/TO/HYPRE/INSTALL # optional #----------------------------------------------------------------------- @@ -49,10 +49,10 @@ MACH_CPP = cpp # C preprocessor command # With MPI MACH_CC_MPI = mpicc # C compiler when using MPI -MACH_CXX_MPI = mpicxx # C++ compiler when using MPI +MACH_CXX_MPI = mpic++ # C++ compiler when using MPI MACH_FC_MPI = gfortran # Fortran 77 compiler when using MPI MACH_F90_MPI = gfortran # Fortran 90 compiler when using MPI -MACH_LD_MPI = mpicxx # Linker when using MPI +MACH_LD_MPI = mpic++ # Linker when using MPI # Without MPI @@ -66,15 +66,14 @@ MACH_LD_NOMPI = g++ # Linker when not using MPI # Machine-dependent defines #----------------------------------------------------------------------- -MACH_DEFINES = -DLINUX -DH5_USE_16_API -fPIC +MACH_DEFINES = -DLINUX -DH5_USE_16_API #----------------------------------------------------------------------- # Compiler flag settings #----------------------------------------------------------------------- -MACH_OMPFLAGS = -fopenmp -MACH_CPPFLAGS = -P -traditional +MACH_CPPFLAGS = -P -traditional MACH_CFLAGS = MACH_CXXFLAGS = MACH_FFLAGS = -fno-second-underscore -ffixed-line-length-132 @@ -87,7 +86,7 @@ MACH_LDFLAGS = MACH_OPT_WARN = -Wall -g MACH_OPT_DEBUG = -g -MACH_OPT_HIGH = -O2 -g +MACH_OPT_HIGH = -O2 MACH_OPT_AGGRESSIVE = -O3 -g #----------------------------------------------------------------------- @@ -98,7 +97,7 @@ LOCAL_INCLUDES_MPI = # MPI includes LOCAL_INCLUDES_HDF5 = -I$(LOCAL_HDF5_INSTALL)/include # HDF5 includes LOCAL_INCLUDES_HYPRE = -I$(LOCAL_HYPRE_INSTALL)/include LOCAL_INCLUDES_PAPI = # PAPI includes -LOCAL_INCLUDES_GRACKLE = -I/home/dskinner6/software/grackle/include +LOCAL_INCLUDES_GRACKLE = -I$(LOCAL_GRACKLE_INSTALL)/include MACH_INCLUDES = $(LOCAL_INCLUDES_HDF5) MACH_INCLUDES_MPI = $(LOCAL_INCLUDES_MPI) @@ -114,8 +113,8 @@ LOCAL_LIBS_MPI = # MPI libraries LOCAL_LIBS_HDF5 = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz LOCAL_LIBS_HYPRE = -L$(LOCAL_HYPRE_INSTALL)/lib -lHYPRE LOCAL_LIBS_PAPI = # PAPI libraries -LOCAL_LIBS_MACH = -lgfortran -lquadmath # Machine-dependent libraries -LOCAL_LIBS_GRACKLE = -L/home/dskinner6/software/grackle/lib -lgrackle +LOCAL_LIBS_MACH = -lgfortran # Machine-dependent libraries +LOCAL_LIBS_GRACKLE = -L$(LOCAL_GRACKLE_INSTALL)/lib -lgrackle MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) MACH_LIBS_MPI = $(LOCAL_LIBS_MPI) diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C index 3842d5842..5c5557203 100644 --- a/src/enzo/ReadParameterFile.C +++ b/src/enzo/ReadParameterFile.C @@ -1374,7 +1374,7 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt) ret += sscanf(line,"MagneticSupernovaDuration = %"FSYM, &MagneticSupernovaDuration); // Rotating Pop III Models - ret += sscanf(line, "PopIII_Rotating = %"ISYM, &PopIII_Rotating); + ret += sscanf(line, "PopIIIRotating = %"ISYM, &PopIIIRotating); /* If the dummy char space was used, then make another. */ diff --git a/src/enzo/SetDefaultGlobalValues.C b/src/enzo/SetDefaultGlobalValues.C index 4bb32817e..8e0838936 100644 --- a/src/enzo/SetDefaultGlobalValues.C +++ b/src/enzo/SetDefaultGlobalValues.C @@ -1055,7 +1055,7 @@ int SetDefaultGlobalValues(TopGridData &MetaData) MagneticSupernovaEnergy = 1.0e51; // Total energy (ergs) injected per star particle (supernova) /* Rotating Pop III Stars Model */ - PopIII_Rotating = 0; // 0 = off + PopIIIRotating = 0; // 0 = off; 1 = rotating; 2 = non-rotating return SUCCESS; } diff --git a/src/enzo/Star_ComputePhotonRates.C b/src/enzo/Star_ComputePhotonRates.C index fb1fe898c..fed5fc805 100644 --- a/src/enzo/Star_ComputePhotonRates.C +++ b/src/enzo/Star_ComputePhotonRates.C @@ -64,44 +64,55 @@ int Star::ComputePhotonRates(const float TimeUnits, int &nbins, float E[], doubl E[2] = 58.0; E[3] = 12.8; _mass = max(min((float)(_mass), 500), 5); - if (_mass > 5 && _mass < 9) { - Q[0] = pow(10.0, 39.29 + 8.55*x); - Q[1] = pow(10.0, 29.24 + 18.49*x); - Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); - Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); + + if (PopIIIRotating = 0) { + if (_mass > 9 && _mass <= 500) { + Q[0] = pow(10.0, 43.61 + 4.9*x - 0.83*x2); + Q[1] = pow(10.0, 42.51 + 5.69*x - 1.01*x2); + Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); + Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); + } else if (_mass > 5 && _mass <= 9) { + Q[0] = pow(10.0, 39.29 + 8.55*x); + Q[1] = pow(10.0, 29.24 + 18.49*x); + Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); + Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); + } // ENDELSE } - else if (_mass >= 9 && _mass <= 120) { - if (PopIII_Rotating == 1) { - if (_age <= 0.8){ - Q[0] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 0)); - Q[1] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 1)); - Q[2] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 2)); - } - else { - Q[0] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 0)); - Q[1] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 1)); - Q[2] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 2)); - } + + else { + if (_mass > 5 && _mass < 9) { + Q[0] = pow(10.0, 39.29 + 8.55*x); + Q[1] = pow(10.0, 29.24 + 18.49*x); + Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); + Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); } - else { - if (_age <= 0.8){ - Q[0] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 0)); - Q[1] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 1)); - Q[2] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 2)); + else if (_mass >= 9 && _mass <= 120) { + if (PopIIIRotating == 1) { + if (_age <= 0.8){ + Q[0] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 0)); + Q[1] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 1)); + Q[2] = pow(10.0, CalculateRotationalPhotonRates(_mass, _age, 2)); + } + else { + Q[0] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 0)); + Q[1] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 1)); + Q[2] = pow(10.0, CalculateRotationalPhotonRates(_mass, 0.8, 2)); + } } - else { - Q[0] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 0)); - Q[1] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 1)); - Q[2] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 2)); + else if (PopIIIRotating == 2) { + if (_age <= 0.8){ + Q[0] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 0)); + Q[1] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 1)); + Q[2] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, _age, 2)); + } + else { + Q[0] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 0)); + Q[1] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 1)); + Q[2] = pow(10.0, CalculateNonRotationalPhotonRates(_mass, 0.8, 2)); + } } + Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); } - Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); - } - else if (_mass > 120 && _mass < 500) { - Q[0] = pow(10.0, 43.61 + 4.9*x - 0.83*x2); - Q[1] = pow(10.0, 42.51 + 5.69*x - 1.01*x2); - Q[2] = pow(10.0, 26.71 + 18.14*x - 3.58*x2); - Q[3] = pow(10.0, 44.03 + 4.59*x - 0.77*x2); } else { for (i = 0; i < nbins; i++) Q[i] = 0.0; diff --git a/src/enzo/TestRadiatingStarParticleInitialize.C b/src/enzo/TestRadiatingStarParticleInitialize.C index edcc6bc7b..6fc0ed417 100644 --- a/src/enzo/TestRadiatingStarParticleInitialize.C +++ b/src/enzo/TestRadiatingStarParticleInitialize.C @@ -8,7 +8,7 @@ / date: July, 2022 / / PURPOSE: -/ Initialize a radiating star particle test in a uniform medium. +/ Initialize a radiating Pop III star particle test in a uniform medium. / Based on the Single Star Particle Test and Photon Test. / / RETURNS: SUCCESS or FAIL @@ -32,12 +32,12 @@ #include "Hierarchy.h" #include "TopGridData.h" -static float PhotonTestInitialFractionHII = 1.2e-5; -static float PhotonTestInitialFractionHeII = 1.0e-14; -static float PhotonTestInitialFractionHeIII = 1.0e-17; -static float PhotonTestInitialFractionHM = 2.0e-9; -static float PhotonTestInitialFractionH2I = 2.0e-20; -static float PhotonTestInitialFractionH2II = 3.0e-14; +static float TestStarParticleInitialFractionHII = 1.2e-5; +static float TestStarParticleInitialFractionHeII = 1.0e-14; +static float TestStarParticleInitialFractionHeIII = 1.0e-17; +static float TestStarParticleInitialFractionHM = 2.0e-9; +static float TestStarParticleInitialFractionH2I = 2.0e-20; +static float TestStarParticleInitialFractionH2II = 3.0e-14; int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData,float *Initialdt) @@ -94,7 +94,7 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr float TestStarParticleStarMass = 100.0; int TestProblemUseMetallicityField = 1; float TestProblemInitialMetallicityFraction = 2e-3; // 0.1 Zsun - float PhotonTestInitialTemperature = 1000; + float TestStarParticleInitialTemperature = 1000; @@ -158,10 +158,10 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr TestStarParticleEnergy, TestStarParticleVelocity, TestStarParticleBField, - PhotonTestInitialTemperature, - PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII, - PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM, - PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II) == FAIL) + TestStarParticleInitialTemperature, + TestStarParticleInitialFractionHII, TestStarParticleInitialFractionHeII, + TestStarParticleInitialFractionHeIII, TestStarParticleInitialFractionHM, + TestStarParticleInitialFractionH2I, TestStarParticleInitialFractionH2II) == FAIL) ENZO_FAIL("Error in TestRadiatingStarParticleInitializeGrid.\n"); /* set up field names and units */ diff --git a/src/enzo/WriteParameterFile.C b/src/enzo/WriteParameterFile.C index 26c15b13a..cb76e0207 100644 --- a/src/enzo/WriteParameterFile.C +++ b/src/enzo/WriteParameterFile.C @@ -1280,7 +1280,7 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData, char *name = NULL) fprintf(fptr,"MagneticSupernovaDuration = %"GSYM"\n",MagneticSupernovaDuration); /* Rotating PopIII Stars Model */ - fprintf(fptr, "PopIII_Rotating = %"ISYM"\n", PopIII_Rotating); + fprintf(fptr, "PopIIIRotating = %"ISYM"\n", PopIIIRotating); /* Output current time */ time_t ID; From 56acb9b11d44143ff8fe107317a2405d6aa5c305 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Mon, 25 Jul 2022 15:13:03 -0400 Subject: [PATCH 06/10] Fixing parameter name --- src/enzo/global_data.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/enzo/global_data.h b/src/enzo/global_data.h index c402c4cab..17ad9a9b5 100644 --- a/src/enzo/global_data.h +++ b/src/enzo/global_data.h @@ -1219,6 +1219,6 @@ EXTERN float MagneticSupernovaDuration; EXTERN float MagneticSupernovaEnergy; /* Rotating Pop III Stars Model */ -EXTERN int PopIII_Rotating; +EXTERN int PopIIIRotating; #endif From 444dbe2d0d46b7ba8e742be6b26a7c7ac16c78f3 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Mon, 25 Jul 2022 15:26:23 -0400 Subject: [PATCH 07/10] Fixing an else if statement. --- src/enzo/Star_ComputePhotonRates.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/enzo/Star_ComputePhotonRates.C b/src/enzo/Star_ComputePhotonRates.C index fed5fc805..7b04a5659 100644 --- a/src/enzo/Star_ComputePhotonRates.C +++ b/src/enzo/Star_ComputePhotonRates.C @@ -79,7 +79,7 @@ int Star::ComputePhotonRates(const float TimeUnits, int &nbins, float E[], doubl } // ENDELSE } - else { + else if (PopIIIRotating == 1 || PopIIIRotating == 2) { if (_mass > 5 && _mass < 9) { Q[0] = pow(10.0, 39.29 + 8.55*x); Q[1] = pow(10.0, 29.24 + 18.49*x); From 221ab20f558f87e88ed41c7250ce6ab9c8ef4583 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Fri, 29 Jul 2022 16:30:39 -0400 Subject: [PATCH 08/10] Fixing radiating star particle to run on multiple cores and compile with photon-no. --- .../TestRadiatingStarParticleSingle.enzo | 30 +- .../TestRadiatingStarParticleSingle.enzotest | 2 +- src/enzo/Grid.h | 17 - src/enzo/Grid_FindNewStarParticles.C | 4 +- ..._TestRadiatingStarParticleInitializeGrid.C | 373 +++++++----------- src/enzo/InitializeNew.C | 10 +- src/enzo/Make.config.objects | 4 +- src/enzo/PhotonGrid_Methods.h | 8 + src/enzo/Star_ComputePhotonRates.C | 2 +- .../TestRadiatingStarParticleInitialize.C | 28 +- 10 files changed, 188 insertions(+), 290 deletions(-) diff --git a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo index 35aef4ecb..f8535a89b 100644 --- a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo +++ b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzo @@ -3,18 +3,20 @@ # # define problem # -ProblemType = 251 // TestRadiatingStarParaticle +ProblemType = 252 // TestRadiatingStarParaticle TopGridRank = 3 TopGridDimensions = 50 50 50 StarParticleFeedback = 8 // Population III Stars TestStarParticleDensity = 1.0 -TestStarParticleEnergy = 6.55410093179e-06 -TestProblemInitialMetallicityFraction = 2e-3 -TestStarParticleStarMass = 20.0 +TestStarParticleEnergy = 6e-10 +TestStarParticleStarPosition = .5 .5 .5 +TestStarParticleStarVelocity = 0 0 0 +#TestProblemInitialMetallicityFraction = 2e-3 +TestStarParticleStarMass = 20.0 // Doesn't immediately work since this is a PopIII star StarEnergyToThermalFeedback = 5.59e-6 StarMassEjectionFraction = 0.25 StarFeedbackKineticFraction = 0.3 -StarMakerExplosionDelayTime = 0.0 +StarMakerExplosionDelayTime = 0.5 // Lifetime of star OutputTemperature = 1 OutputCoolingTime = 1 Initialdt = 0.0010017410642 @@ -44,8 +46,12 @@ Gamma = 1.6666667 # StaticHierarchy = 1 // static hierarchy RefineBy = 2 // refinement factor +MaximumRefinementLevel = 0 +CellFlaggingMethod = 2 + # # Radiation +RadiativeTransfer = 1 RadiativeTransferRaysPerCell = 5.100000 RadiativeTransferInitialHEALPixLevel = 1 RadiativeTransferHydrogenOnly = 0 @@ -55,16 +61,14 @@ RadiativeTransferAdaptiveTimestep = 1 RadiativeTransferRadiationPressure = 1 RadiativeTransferPhotonMergeRadius = 3.0 RadiativeTransferSourceClustering = 1 - +RadiativeCooling = 1 # -RadiativeTransfer = 1 - MaximumRefinementLevel = 0 CellFlaggingMethod = 0 # # Cooling MultiSpecies = 1 -RadiativeCooling = 0 +#RadiativeCooling = 0 #MultiSpecies = 1 #MetalCooling = 3 #CloudyCoolingGridFile = solar_2008_3D_metals.h5 @@ -74,6 +78,7 @@ RadiativeCooling = 0 #TestProblemHydrogenFractionByMass = 0.75 #TestProblemInitialMetallicityFraction = 2e-3 #TestProblemInitialHIFraction = 0.998 +#TestProblemInitialHIIFraction #TestProblemInitialHeIFraction = 1.0 #TestProblemInitialHeIIFraction = 1.0e-20 #TestProblemInitialHeIIIIFraction = 1.0e-20 @@ -81,12 +86,13 @@ RadiativeCooling = 0 # set some misc global parameters # -MixSpeciesAndColors = 1 -PopIIISupernovaUseColour = 1 +#MixSpeciesAndColors = 1 +#PopIIISupernovaUseColour = 1 +#PopIIIHeliumIonization = 1 # # Rotating Pop III Stars # -PopIII_Rotating = 1 \ No newline at end of file +PopIIIRotating = 0 \ No newline at end of file diff --git a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest index 7867f0e89..aa9bd8f43 100644 --- a/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest +++ b/run/StarParticle/RadiatingStarParticleSingleTest/TestRadiatingStarParticleSingle.enzotest @@ -1,2 +1,2 @@ name = 'TestRadiatingStarParticleSingle' -radiation = 'ray' \ No newline at end of file +radiation = 'ray' diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 185ba4c05..d54d9ef1c 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -2091,23 +2091,6 @@ int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], FLOAT TestStarParticleStarVelocity[], FLOAT TestStarParticleStarPosition[]); -/* Radiating Star Particle test: initialize particle */ - int TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass, - float *Initialdt, - FLOAT TestStarParticleStarVelocity[], - FLOAT TestStarParticleStarPosition[], - float UniformDensity, - float UniformTotalEnergy, - float UniformVelocity[], - float UniformBField[], - float InitialTemperature, - float TestStarParticleInitialFractionHII, - float TestStarParticleInitialFractionHeII, - float TestStarParticleInitialFractionHeIII, - float TestStarParticleInitialFractionHM, - float TestStarParticleInitialFractionH2I, - float TestStarParticleInitialFractionH2II); - /* Gravity Test: initialize grid. */ int TestGravityInitializeGrid(float CentralDensity, diff --git a/src/enzo/Grid_FindNewStarParticles.C b/src/enzo/Grid_FindNewStarParticles.C index a8f6ce3e2..1e23d1e42 100644 --- a/src/enzo/Grid_FindNewStarParticles.C +++ b/src/enzo/Grid_FindNewStarParticles.C @@ -65,8 +65,10 @@ int grid::FindNewStarParticles(int level) calls to the IMF. */ if (ParticleType[i] == -PARTICLE_TYPE_SINGLE_STAR) - if (PopIIIInitialMassFunction == FALSE) + if (PopIIIInitialMassFunction == FALSE){ + if (ProblemType != 252) NewStar->AssignFinalMass(PopIIIStarMass); + } if (ParticleType[i] == -PARTICLE_TYPE_SIMPLE_SOURCE) NewStar->AssignFinalMass(PopIIIStarMass); InsertStarAfter(Stars, NewStar); diff --git a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C index c879983f8..434a01900 100644 --- a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C +++ b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C @@ -31,14 +31,6 @@ int GetUnits(float *DensityUnits, float *LengthUnits, float *TemperatureUnits, float *TimeUnits, float *VelocityUnits, double *MassUnits, FLOAT Time); -// Used to compute Bonner-Ebert density profile -double ph_BE(double r); -double ph_q(double r); -double ph_Ang(double a1, double a2, double R, double r); - -// Returns random velocity from Maxwellian distribution -double ph_Maxwellian(double c_tilda, double vel_unit, double mu, double gamma); - static int TestStarParticleParticleCount = 0; int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass, @@ -47,23 +39,15 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass FLOAT TestStarParticleStarPosition[], float TestStarParticleDensity, float TestStarParticleEnergy, - float TestStarParticleVelocity[], - float TestStarParticleBField[], - float InitialTemperature, - float TestStarParticleInitialFractionHII, - float TestStarParticleInitialFractionHeII, - float TestStarParticleInitialFractionHeIII, - float TestStarParticleInitialFractionHM, - float TestStarParticleInitialFractionH2I, - float TestStarParticleInitialFractionH2II) + float TestStarParticleVelocity[]) { /* declarations */ - float CentralMass = 1.0; + float CentralMass = 1.0, temperature; int dim, i, j, k, m, field, size, active_size, index, cindex; int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, - DINum, DIINum, HDINum, MetalNum, MetalIaNum, B1Num, B2Num, B3Num, PhiNum, CRNum; + DINum, DIINum, HDINum, MetalNum, MetalIaNum; int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum, kdissH2IINum, kphHMNum, RPresNum1, RPresNum2, RPresNum3; @@ -71,11 +55,6 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass float TestInitialdt = *Initialdt, *density_field = NULL, *HII_field = NULL, *HeII_field = NULL, *HeIII_field = NULL, *Temperature_field = NULL; - - /* Return if this doesn't concern us. */ - - if (ProcessorNumber != MyProcessorNumber) - return SUCCESS; /* create fields */ @@ -86,56 +65,33 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass FieldType[NumberOfBaryonFields++] = InternalEnergy; int vel = NumberOfBaryonFields; FieldType[NumberOfBaryonFields++] = Velocity1; - if (GridRank > 1 || HydroMethod > 2) + if (GridRank > 1) FieldType[NumberOfBaryonFields++] = Velocity2; - if (GridRank > 2 || HydroMethod > 2) + if (GridRank > 2) FieldType[NumberOfBaryonFields++] = Velocity3; - - if ( CRModel ) { - CRNum = NumberOfBaryonFields; - FieldType[NumberOfBaryonFields++] = CRDensity; - } - - if (WritePotential) - FieldType[NumberOfBaryonFields++] = GravPotential; - - - int colorfields = NumberOfBaryonFields; - - // Enzo's standard multispecies (primordial chemistry - H, D, He) - if (TestProblemData.MultiSpecies) { - FieldType[DeNum = NumberOfBaryonFields++] = ElectronDensity; - FieldType[HINum = NumberOfBaryonFields++] = HIDensity; - FieldType[HIINum = NumberOfBaryonFields++] = HIIDensity; - FieldType[HeINum = NumberOfBaryonFields++] = HeIDensity; - FieldType[HeIINum = NumberOfBaryonFields++] = HeIIDensity; - FieldType[HeIIINum = NumberOfBaryonFields++] = HeIIIDensity; - if (TestProblemData.MultiSpecies > 1) { - FieldType[HMNum = NumberOfBaryonFields++] = HMDensity; - FieldType[H2INum = NumberOfBaryonFields++] = H2IDensity; - FieldType[H2IINum = NumberOfBaryonFields++] = H2IIDensity; + if (MultiSpecies) { + FieldType[DeNum = NumberOfBaryonFields++] = ElectronDensity; + FieldType[HINum = NumberOfBaryonFields++] = HIDensity; + FieldType[HIINum = NumberOfBaryonFields++] = HIIDensity; + FieldType[HeINum = NumberOfBaryonFields++] = HeIDensity; + FieldType[HeIINum = NumberOfBaryonFields++] = HeIIDensity; + FieldType[HeIIINum = NumberOfBaryonFields++] = HeIIIDensity; + if (MultiSpecies > 1) { + FieldType[HMNum = NumberOfBaryonFields++] = HMDensity; + FieldType[H2INum = NumberOfBaryonFields++] = H2IDensity; + FieldType[H2IINum = NumberOfBaryonFields++] = H2IIDensity; } - if (TestProblemData.MultiSpecies > 2) { + if (MultiSpecies > 2) { FieldType[DINum = NumberOfBaryonFields++] = DIDensity; FieldType[DIINum = NumberOfBaryonFields++] = DIIDensity; FieldType[HDINum = NumberOfBaryonFields++] = HDIDensity; } } - // Metal fields, including the standard 'metallicity' as well - // as two extra fields - if (TestProblemData.UseMetallicityField) { - FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity; - - if (StarMakerTypeIaSNe) - FieldType[MetalIaNum = NumberOfBaryonFields++] = MetalSNIaDensity; - - if(TestProblemData.MultiMetals){ - FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0; - FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1; - } + if (TestProblemData.UseMetallicityField){ + MetalNum = NumberOfBaryonFields; + FieldType[NumberOfBaryonFields++] = Metallicity; } - if (RadiativeTransfer && (MultiSpecies < 1)) { ENZO_FAIL("Grid_TestRadiatingStarParticleInitialize: Radiative Transfer but not MultiSpecies set"); } @@ -146,13 +102,13 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; if (RadiativeTransferHydrogenOnly == FALSE) { - FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; - FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; + FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; + FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; } if (MultiSpecies > 1) { - FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; - FieldType[kdissH2IINum = NumberOfBaryonFields++] = kdissH2II; - FieldType[kphHMNum = NumberOfBaryonFields++] = kphHM; + FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; + FieldType[kdissH2IINum = NumberOfBaryonFields++] = kdissH2II; + FieldType[kphHMNum = NumberOfBaryonFields++] = kphHM; } } @@ -165,6 +121,11 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass NumberOfPhotonPackages = 0; PhotonPackages-> NextPackage= NULL; + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) + return SUCCESS; + /* Get Units. */ float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, @@ -176,20 +137,23 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass ENZO_FAIL("Error in GetUnits.\n"); } - /* Set Central Mass in simulation units */ - - CentralMass = TestStarParticleStarMass*1.99e33* pow(LengthUnits*CellWidth[0][0],-3.0)/DensityUnits; - - printf("Central Mass: %f \n",CentralMass); - /* Set number of particles for this grid and allocate space. */ NumberOfParticles = 1; NumberOfParticleAttributes = 4; this->AllocateNewParticles(NumberOfParticles); printf("Allocated %d particles\n", NumberOfParticles); + + /* Set up the baryon field. */ + /* compute size of fields */ + + size = 1; + for (dim = 0; dim < GridRank; dim++) + size *= GridDimension[dim]; + + /* allocate fields */ - this->AllocateGrids(); + this->AllocateGrids(); /* Initialize radiation fields */ // need this @@ -197,190 +161,129 @@ int grid::TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass ENZO_FAIL("\nError in InitializeRadiativeTransferFields.\n"); } - /* Set particle IDs and types */ - - for (i = 0; i < NumberOfParticles; i++) { - ParticleNumber[i] = i; - ParticleType[i] = -PopIII; - ParticleAttribute[0][i] = Time + 1e-7; - } - ParticleMass[0] = CentralMass; - - /* Set central particle. */ - for (dim = 0; dim < GridRank; dim++) { - ParticlePosition[dim][0] = TestStarParticleStarPosition[dim]* - (DomainLeftEdge[dim]+DomainRightEdge[dim]) + 0.5*CellWidth[0][0]; - ParticleVelocity[dim][0] = TestStarParticleStarVelocity[dim]*1e5*TimeUnits/LengthUnits; + for (i = 0; i < size; i++) { + BaryonField[0][i] = TestStarParticleDensity; + BaryonField[1][i] = TestStarParticleEnergy; } - if (STARFEED_METHOD(UNIGRID_STAR)) ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; - if (STARFEED_METHOD(MOM_STAR)) - if(StarMakerExplosionDelayTime >= 0.0) - ParticleAttribute[1][0] = 1.0; - else - ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; - - ParticleAttribute[2][0] = 0.0; // Metal fraction - ParticleAttribute[3][0] = 0.0; // metalfSNIa - - if (this->FindNewStarParticles(GridLevel) == FAIL) { - ENZO_FAIL("Error in grid::FindNewStarParticles."); - } - - /* Reset particle type to be positive*/ - for (i = 0; i < NumberOfParticles; i++) { - ParticleNumber[i] = i; - ParticleType[i] = PopIII; - } - Star *cstar; - for (cstar = Stars; cstar; cstar = cstar->NextStar) - cstar->type = PopIII; - - - /* Set up the baryon field. */ // need thid - /* compute size of fields */ - active_size = 1; - int ActiveDims[MAX_DIMENSION]; - for (dim = 0; dim < GridRank; dim++) { - ActiveDims[dim] = GridEndIndex[dim] - GridStartIndex[dim] + 1; - active_size *= ActiveDims[dim]; - } - - /* Loop over the mesh. */ - float density, dens1, Velocity[MAX_DIMENSION], - temperature, temp1, sigma, sigma1, colour, outer_radius; - float HII_Fraction, HeII_Fraction, HeIII_Fraction, H2I_Fraction; - FLOAT r, x, y = 0, z = 0; - int n = 0; - - for (k = 0; k < GridDimension[2]; k++) - for (j = 0; j < GridDimension[1]; j++) - for (i = 0; i < GridDimension[0]; i++, n++) { - - /* Compute position */ // need this - - x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; - if (GridRank > 1) - y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j]; - if (GridRank > 2) - z = CellLeftEdge[2][k] + 0.5*CellWidth[2][k]; - - if (i >= GridStartIndex[0] && i <= GridEndIndex[0] && - j >= GridStartIndex[1] && j <= GridEndIndex[1] && - k >= GridStartIndex[2] && k <= GridEndIndex[2]) { - cindex = (i-GridStartIndex[0]) + ActiveDims[0] * - ((j-GridStartIndex[1]) + (k-GridStartIndex[2])*ActiveDims[1]); - if (density_field != NULL) - density = max(density_field[cindex], 1e-6); - else - density = 1.0; - if (HII_field != NULL) - HII_Fraction = HII_field[cindex]; - else - HII_Fraction = TestStarParticleInitialFractionHII; - if (HeII_field != NULL) - HeII_Fraction = HeII_field[cindex]; - else - HeII_Fraction = TestStarParticleInitialFractionHeII; - if (HeIII_field != NULL) - HeIII_Fraction = HeIII_field[cindex]; - else - HeIII_Fraction = TestStarParticleInitialFractionHeIII; - if (Temperature_field != NULL) - temperature = temp1 = Temperature_field[cindex]; - else - temperature = temp1 = InitialTemperature; - } - - /* set density, total energy */ // check - BaryonField[0][n] = TestStarParticleDensity; - BaryonField[1][n] = TestStarParticleEnergy; - - /* set velocities */ - + /* set velocities */ + for (dim = 0; dim < GridRank; dim++) - BaryonField[vel+dim][n] = Velocity[dim] + TestStarParticleVelocity[dim]; - - /* Set internal energy if necessary. */ //check + for (i = 0; i < size; i++) + BaryonField[vel+dim][i] = TestStarParticleVelocity[dim]; + /* Set internal energy if necessary. */ + if (DualEnergyFormalism) - BaryonField[2][n] = BaryonField[1][n]; - - /* If doing multi-species (HI, etc.), set these. */ + for (i = 0; i < size; i++) + BaryonField[2][i] = BaryonField[1][i]; - if (MultiSpecies > 0) { - BaryonField[HIINum][n] = TestProblemData.HII_Fraction * - TestProblemData.HydrogenFractionByMass * TestStarParticleDensity; + for (i = 0; i < size; i++){ - BaryonField[HeIINum][n] = TestProblemData.HeII_Fraction * - TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); + if(MultiSpecies > 0){ - BaryonField[HeIIINum][n] = TestProblemData.HeIII_Fraction * - TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); + BaryonField[HIINum][i] = TestProblemData.HII_Fraction * + TestProblemData.HydrogenFractionByMass * TestStarParticleDensity; + + BaryonField[HeIINum][i] = TestProblemData.HeII_Fraction * + TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); - BaryonField[HeINum][n] = - (1.0 - TestProblemData.HydrogenFractionByMass)*TestStarParticleDensity - - BaryonField[HeIINum][n] - BaryonField[HeIIINum][n]; - } + BaryonField[HeIIINum][i] = TestProblemData.HeIII_Fraction * + TestStarParticleDensity * (1.0-TestProblemData.HydrogenFractionByMass); - if (MultiSpecies > 1) { - BaryonField[HMNum][n] = TestProblemData.HM_Fraction * - BaryonField[HIINum][n]; + BaryonField[HeINum][i] = + (1.0 - TestProblemData.HydrogenFractionByMass)*TestStarParticleDensity - + BaryonField[HeIINum][i] - BaryonField[HeIIINum][i]; - BaryonField[H2INum][n] = TestProblemData.H2I_Fraction * - BaryonField[0][n] * TestProblemData.HydrogenFractionByMass; + if(MultiSpecies > 1){ + BaryonField[HMNum][i] = TestProblemData.HM_Fraction * + BaryonField[HIINum][i]; - BaryonField[H2IINum][n] = TestProblemData.H2II_Fraction * 2.0 * - BaryonField[HIINum][n]; + BaryonField[H2INum][i] = TestProblemData.H2I_Fraction * + BaryonField[0][i] * TestProblemData.HydrogenFractionByMass; - BaryonField[kdissH2INum][n] = 0.0; - BaryonField[kphHMNum][n] = 0.0; - BaryonField[kdissH2IINum][n] = 0.0; - } + BaryonField[H2IINum][i] = TestProblemData.H2II_Fraction * 2.0 * + BaryonField[HIINum][i]; - // HI density is calculated by subtracting off the various ionized fractions - // from the total - BaryonField[HINum][n] = TestProblemData.HydrogenFractionByMass*BaryonField[0][n] - - BaryonField[HIINum][n]; - - if (MultiSpecies > 1) - BaryonField[HINum][n] -= (BaryonField[HMNum][n] + BaryonField[H2IINum][n] - + BaryonField[H2INum][n]); - - BaryonField[DeNum][n] = BaryonField[HIINum][n] + - 0.25*BaryonField[HeIINum][n] + 0.5*BaryonField[HeIIINum][n]; + BaryonField[kdissH2INum][i] = 0.0; + BaryonField[kphHMNum][i] = 0.0; + BaryonField[kdissH2IINum][i] = 0.0; + } - if (MultiSpecies > 1) - BaryonField[DeNum][n] += 0.5*BaryonField[H2IINum][n] - - BaryonField[HMNum][n]; + BaryonField[HINum][i] = TestProblemData.HydrogenFractionByMass*BaryonField[0][i] + - BaryonField[HIINum][i]; + if (MultiSpecies > 1) + BaryonField[HINum][i] -= (BaryonField[HMNum][i] + BaryonField[H2IINum][i] + + BaryonField[H2INum][i]); + + BaryonField[DeNum][i] = BaryonField[HIINum][i] + + 0.25*BaryonField[HeIINum][i] + 0.5*BaryonField[HeIIINum][i]; + if (MultiSpecies > 1) + BaryonField[DeNum][i] += 0.5*BaryonField[H2IINum][i] - + BaryonField[HMNum][i]; + + if(MultiSpecies > 2){ + BaryonField[DINum ][i] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HINum][i]; + BaryonField[DIINum][i] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HIINum][i]; + BaryonField[HDINum][i] = 0.75 * TestProblemData.DeuteriumToHydrogenRatio * BaryonField[H2INum][i]; + } - /* Set Deuterium species (assumed to be negligible). */ + } // if(TestProblemData.MultiSpecies) - if (MultiSpecies > 2) { - BaryonField[DINum ][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HINum][n]; - BaryonField[DIINum][n] = TestProblemData.DeuteriumToHydrogenRatio * BaryonField[HIINum][n]; - BaryonField[HDINum][n] = 0.75 * TestProblemData.DeuteriumToHydrogenRatio * BaryonField[H2INum][n]; - } + // metallicity fields + if(TestProblemData.UseMetallicityField){ + BaryonField[MetalNum][i] = TestProblemData.MetallicityField_Fraction* TestStarParticleDensity; - /* Set energy (thermal and then total if necessary). */ + } - if (MultiSpecies) { +if (MultiSpecies) { mu_data = - 0.25*(BaryonField[HeINum][n] + BaryonField[HeIINum][n] + - BaryonField[HeIIINum][n] ) + - BaryonField[HINum][n] + BaryonField[HIINum][n] + - BaryonField[DeNum][n]; + 0.25*(BaryonField[HeINum][i] + BaryonField[HeIINum][i] + + BaryonField[HeIIINum][i] ) + + BaryonField[HINum][i] + BaryonField[HIINum][i] + + BaryonField[DeNum][i]; if (MultiSpecies > 1) - mu_data += BaryonField[HMNum][n] + - 0.5*(BaryonField[H2INum][n] + BaryonField[H2IINum][n]); - mu_data = BaryonField[0][n] / mu_data; + mu_data += BaryonField[HMNum][i] + + 0.5*(BaryonField[H2INum][i] + BaryonField[H2IINum][i]); + mu_data = BaryonField[0][i] / mu_data; } else mu_data = mu; - BaryonField[1][n] = temperature/TemperatureUnits/ - ((Gamma-1.0)*mu_data); - } + } // for (i = 0; i < size; i++) + + /* Set Central Mass in simulation units */ + + CentralMass = TestStarParticleStarMass*1.99e33* pow(LengthUnits*CellWidth[0][0],-3.0)/DensityUnits; + + printf("Central Mass: %f \n",CentralMass); + + /* Set particle IDs and types */ + + for (i = 0; i < NumberOfParticles; i++) { + ParticleNumber[i] = i; + ParticleType[i] = -PopIII; + ParticleAttribute[0][i] = Time + 1e-7; + ParticleAttribute[1][i] = StarMakerExplosionDelayTime; + } + ParticleMass[0] = CentralMass; + + /* Set central particle. */ + for (dim = 0; dim < GridRank; dim++) { + ParticlePosition[dim][0] = TestStarParticleStarPosition[dim]* + (DomainLeftEdge[dim]+DomainRightEdge[dim]) + 0.5*CellWidth[0][0]; + ParticleVelocity[dim][0] = TestStarParticleStarVelocity[dim]*1e5*TimeUnits/LengthUnits; + } + + if (STARFEED_METHOD(UNIGRID_STAR)) ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; + if (STARFEED_METHOD(MOM_STAR)) + if(StarMakerExplosionDelayTime >= 0.0) + ParticleAttribute[1][0] = 1.0; + else + ParticleAttribute[1][0] = 10.0 * Myr_s/TimeUnits; + + ParticleAttribute[2][0] = 0.0; // Metal fraction + ParticleAttribute[3][0] = 0.0; // metalfSNIa delete [] density_field; delete [] HII_field; diff --git a/src/enzo/InitializeNew.C b/src/enzo/InitializeNew.C index 2a58f3c46..9937815c9 100644 --- a/src/enzo/InitializeNew.C +++ b/src/enzo/InitializeNew.C @@ -89,8 +89,6 @@ int StratifiedMediumExplosionInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr TopGridData &MetaData); int TestStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData, float *Initialdt); -int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, - TopGridData &MetaData, float *Initialdt); int KHInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData); int NohInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, @@ -201,6 +199,8 @@ int RHIonizationSteepInitialize(FILE *fptr, FILE *Outfptr, int CosmoIonizationInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData, int local); +int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, + TopGridData &MetaData, float *Initialdt); #endif /* TRANSFER */ @@ -765,10 +765,12 @@ int InitializeNew(char *filename, HierarchyEntry &TopGrid, // Insert new problem intializer here... - // 251) Test a star particle explosion - if (ProblemType == 251) + // 252) Test a star particle explosion + #ifdef TRANSFER + if (ProblemType == 252) ret = TestRadiatingStarParticleInitialize(fptr, Outfptr, TopGrid, MetaData, Initialdt); + #endif /* TRANSFER */ if (ret == INT_UNDEFINED) { ENZO_VFAIL("Problem Type %"ISYM" undefined.\n", ProblemType) diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects index d33672629..48472677a 100644 --- a/src/enzo/Make.config.objects +++ b/src/enzo/Make.config.objects @@ -655,7 +655,6 @@ OBJS_CONFIG_LIB = \ Grid_TestGravitySphereInitializeGrid.o \ Grid_TestOrbitInitializeGrid.o \ Grid_TestStarParticleInitializeGrid.o \ - Grid_TestRadiatingStarParticleInitializeGrid.o \ Grid_TracerParticleCreateParticles.o \ Grid_TracerParticleOutputData.o \ Grid_TracerParticleSetVelocity.o \ @@ -944,7 +943,6 @@ OBJS_CONFIG_LIB = \ TestGravitySphereInitialize.o \ TestOrbitInitialize.o \ TestStarParticleInitialize.o \ - TestRadiatingStarParticleInitialize.o \ TracerParticleCreation.o \ TransposeRegionOverlap.o \ tscint1d.o \ @@ -1107,6 +1105,7 @@ POBJS_CONFIG_LIB = \ Grid_SetSubgridMarkerFromSubgrid.o \ Grid_SubgridMarkerPostParallel.o \ Grid_Shine.o \ + Grid_TestRadiatingStarParticleInitializeGrid.o \ Grid_TransportPhotonPackages.o \ Grid_WalkPhotonPackage.o \ LinkedListRoutines.o \ @@ -1138,6 +1137,7 @@ POBJS_CONFIG_LIB = \ StarParticleRadTransfer.o \ TemperatureFieldToolsForComptonHeating.o \ TemperatureFieldToolsH2Shielding.o \ + TestRadiatingStarParticleInitialize.o \ WritePhotonSources.o OBJS_HYDRO_RK = \ diff --git a/src/enzo/PhotonGrid_Methods.h b/src/enzo/PhotonGrid_Methods.h index f4d53b8e4..b3df0d13b 100644 --- a/src/enzo/PhotonGrid_Methods.h +++ b/src/enzo/PhotonGrid_Methods.h @@ -509,6 +509,14 @@ int PhotonTestInitializeGrid(int NumberOfSpheres, char *HeIIIFractionFilename, char *TemperatureFilename); +int TestRadiatingStarParticleInitializeGrid(float TestStarParticleStarMass, + float *Initialdt, + FLOAT TestStarParticleStarVelocity[], + FLOAT TestStarParticleStarPosition[], + float TestStarParticleDensity, + float TestStarParticleEnergy, + float TestStarParticleVelocity[]); + /************************************************************************/ int StarParticlesToRadSources(FLOAT Time, double ConversionFactor); diff --git a/src/enzo/Star_ComputePhotonRates.C b/src/enzo/Star_ComputePhotonRates.C index 7b04a5659..920db2fa5 100644 --- a/src/enzo/Star_ComputePhotonRates.C +++ b/src/enzo/Star_ComputePhotonRates.C @@ -65,7 +65,7 @@ int Star::ComputePhotonRates(const float TimeUnits, int &nbins, float E[], doubl E[3] = 12.8; _mass = max(min((float)(_mass), 500), 5); - if (PopIIIRotating = 0) { + if (PopIIIRotating == 0) { if (_mass > 9 && _mass <= 500) { Q[0] = pow(10.0, 43.61 + 4.9*x - 0.83*x2); Q[1] = pow(10.0, 42.51 + 5.69*x - 1.01*x2); diff --git a/src/enzo/TestRadiatingStarParticleInitialize.C b/src/enzo/TestRadiatingStarParticleInitialize.C index 6fc0ed417..0285b58dd 100644 --- a/src/enzo/TestRadiatingStarParticleInitialize.C +++ b/src/enzo/TestRadiatingStarParticleInitialize.C @@ -42,7 +42,7 @@ static float TestStarParticleInitialFractionH2II = 3.0e-14; int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData,float *Initialdt) { - const char *DensName = "Density"; /* make const */ + const char *DensName = "Density"; const char *TEName = "TotalEnergy"; const char *GEName = "GasEnergy"; const char *Vel1Name = "x-velocity"; @@ -90,11 +90,10 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr float TestStarParticleVelocity[3] = {0.0, 0.0, 0.0}; FLOAT TestStarParticleStarVelocity[3] = {0.0, 0.0, 0.0}; FLOAT TestStarParticleStarPosition[3] = {0.5, 0.5, 0.5}; - float TestStarParticleBField[3] = {0.0, 0.0, 0.0}; + float TestStarParticleBField[3] = {0.0, 0.0, 0.0}; float TestStarParticleStarMass = 100.0; int TestProblemUseMetallicityField = 1; float TestProblemInitialMetallicityFraction = 2e-3; // 0.1 Zsun - float TestStarParticleInitialTemperature = 1000; @@ -105,7 +104,7 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr /* read input from file */ - + rewind(fptr); while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) { ret = 0; @@ -131,12 +130,12 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr ret += sscanf(line, "TestProblemUseMetallicityField = %"ISYM, &TestProblemData.UseMetallicityField); ret += sscanf(line, "TestProblemInitialMetallicityFraction = %"FSYM, &TestProblemData.MetallicityField_Fraction); - ret += sscanf(line, "TestProblemInitialHIFraction = %"FSYM, &TestProblemData.HI_Fraction); - ret += sscanf(line, "TestProblemInitialHIIFraction = %"FSYM, &TestProblemData.HII_Fraction); - ret += sscanf(line, "TestProblemInitialHeIFraction = %"FSYM, &TestProblemData.HeI_Fraction); - ret += sscanf(line, "TestProblemInitialHeIIFraction = %"FSYM, &TestProblemData.HeII_Fraction); - ret += sscanf(line, "TestProblemInitialHeIIIIFraction = %"FSYM, &TestProblemData.HeIII_Fraction); - ret += sscanf(line, "TestProblemHydrogenFractionByMass = %"FSYM, &TestProblemData.HydrogenFractionByMass); + //ret += sscanf(line, "TestProblemInitialHIFraction = %"FSYM, &TestProblemData.HI_Fraction); + //ret += sscanf(line, "TestProblemInitialHIIFraction = %"FSYM, &TestProblemData.HII_Fraction); + //ret += sscanf(line, "TestProblemInitialHeIFraction = %"FSYM, &TestProblemData.HeI_Fraction); + ////ret += sscanf(line, "TestProblemInitialHeIIFraction = %"FSYM, &TestProblemData.HeII_Fraction); + //ret += sscanf(line, "TestProblemInitialHeIIIIFraction = %"FSYM, &TestProblemData.HeIII_Fraction); + //ret += sscanf(line, "TestProblemHydrogenFractionByMass = %"FSYM, &TestProblemData.HydrogenFractionByMass); /* if the line is suspicious, issue a warning */ @@ -156,12 +155,7 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr TestStarParticleStarPosition, TestStarParticleDensity, TestStarParticleEnergy, - TestStarParticleVelocity, - TestStarParticleBField, - TestStarParticleInitialTemperature, - TestStarParticleInitialFractionHII, TestStarParticleInitialFractionHeII, - TestStarParticleInitialFractionHeIII, TestStarParticleInitialFractionHM, - TestStarParticleInitialFractionH2I, TestStarParticleInitialFractionH2II) == FAIL) + TestStarParticleVelocity) == FAIL) ENZO_FAIL("Error in TestRadiatingStarParticleInitializeGrid.\n"); /* set up field names and units */ @@ -174,7 +168,7 @@ int TestRadiatingStarParticleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntr DataLabel[count++] = (char*) Vel1Name; DataLabel[count++] = (char*) Vel2Name; DataLabel[count++] = (char*) Vel3Name; - if (TestProblemData.MultiSpecies){ + if (MultiSpecies){ DataLabel[count++] = (char*) ElectronName; DataLabel[count++] = (char*) HIName; DataLabel[count++] = (char*) HIIName; From 8e7cac45b885e1fd7b380412e2c2d1dcb6899134 Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Mon, 15 May 2023 12:31:27 -0400 Subject: [PATCH 09/10] Update src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C Co-authored-by: Greg Bryan --- src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C index 434a01900..aeaec466c 100644 --- a/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C +++ b/src/enzo/Grid_TestRadiatingStarParticleInitializeGrid.C @@ -254,7 +254,7 @@ if (MultiSpecies) { /* Set Central Mass in simulation units */ - CentralMass = TestStarParticleStarMass*1.99e33* pow(LengthUnits*CellWidth[0][0],-3.0)/DensityUnits; + CentralMass = TestStarParticleStarMass*SolarMass* pow(LengthUnits*CellWidth[0][0],-3.0)/DensityUnits; printf("Central Mass: %f \n",CentralMass); From c60970f36f7365a38469bf828f8e9e034600a3cb Mon Sep 17 00:00:00 2001 From: Danielle Skinner Date: Mon, 15 May 2023 12:48:07 -0400 Subject: [PATCH 10/10] Addressing comments made by Greg Bryan. --- src/enzo/Grid_FindNewStarParticles.C | 8 +++++--- src/enzo/Star_ComputePhotonRates.C | 4 +++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/enzo/Grid_FindNewStarParticles.C b/src/enzo/Grid_FindNewStarParticles.C index 1e23d1e42..6caa80e5e 100644 --- a/src/enzo/Grid_FindNewStarParticles.C +++ b/src/enzo/Grid_FindNewStarParticles.C @@ -64,11 +64,13 @@ int grid::FindNewStarParticles(int level) merging new (massless) star particles to avoid unnecessary calls to the IMF. */ - if (ParticleType[i] == -PARTICLE_TYPE_SINGLE_STAR) + if (ParticleType[i] == -PARTICLE_TYPE_SINGLE_STAR){ if (PopIIIInitialMassFunction == FALSE){ - if (ProblemType != 252) - NewStar->AssignFinalMass(PopIIIStarMass); + if (ProblemType != 252){ + NewStar->AssignFinalMass(PopIIIStarMass); + } } + } if (ParticleType[i] == -PARTICLE_TYPE_SIMPLE_SOURCE) NewStar->AssignFinalMass(PopIIIStarMass); InsertStarAfter(Stars, NewStar); diff --git a/src/enzo/Star_ComputePhotonRates.C b/src/enzo/Star_ComputePhotonRates.C index 920db2fa5..0f9b7a895 100644 --- a/src/enzo/Star_ComputePhotonRates.C +++ b/src/enzo/Star_ComputePhotonRates.C @@ -45,7 +45,9 @@ int Star::ComputePhotonRates(const float TimeUnits, int &nbins, float E[], doubl _mass = this->FinalMass; else _mass = this->Mass; - _age = this->ReturnAge()/this->LifeTime; //Percentage of lifetime + + _age = this->ReturnAge()/this->LifeTime; //Percentage of lifetime + x = log10((float)(_mass)); x2 = x*x;