diff --git a/input/hm12_photorates.dat b/input/hm12_photorates.dat index 482762782..ace3dc594 100644 --- a/input/hm12_photorates.dat +++ b/input/hm12_photorates.dat @@ -1,27 +1,26 @@ - # Optically thin H and He photoionization and heating rates + Compton heating rate - # - # The following 11 lines are model identification flags - # CUBA V6.3 Nov 2011 - # 0.1000E+01 - # 0.4000E+00 0.1300E+01 0.2000E+00 0.4400E+00 0.1570E+01 0.9000E+00 - # 0.1200E+06 0.1000E+05 0.5000E+04 0.1300E+04 - # 0.3000E+00 0.7000E+00 0.8900E-01 0.7000E+00 - # 0.1000E+01 - # Inputs/salp_SFRD1parfit.hiz.dat - # Inputs/column_dist.dat - # 0.1000E+01 - # 0.1800E-03 0.3400E+01 - # 0.1000E+02 - # - # The following 7 lines are column identification labels - # z=redshift, HI_i=hydrogen phoionization rate, HI_h=hydrogen photoheating rate - # HeI_i,h=photoionization and photoheating rates for neutral helium - # HeII_i,h=photoionization and photoheating rates for singly ionized helium - # Compton= Compton heating rate. - # All photoionization and photoheating rates are *per ion* and have units of 1/s and eV/s, respectively. - # The Compton heating rate is *per electron* and has units of eV/s. - # z HI_i HI_h HeI_i HeI_h HeII_i HeII_h Compton - +# Optically thin H and He photoionization and heating rates + Compton heating rate +# +# The following 11 lines are model identification flags +# CUBA V6.3 Nov 2011 +# 0.1000E+01 +# 0.4000E+00 0.1300E+01 0.2000E+00 0.4400E+00 0.1570E+01 0.9000E+00 +# 0.1200E+06 0.1000E+05 0.5000E+04 0.1300E+04 +# 0.3000E+00 0.7000E+00 0.8900E-01 0.7000E+00 +# 0.1000E+01 +# Inputs/salp_SFRD1parfit.hiz.dat +# Inputs/column_dist.dat +# 0.1000E+01 +# 0.1800E-03 0.3400E+01 +# 0.1000E+02 +# +# The following 7 lines are column identification labels +# z=redshift, HI_i=hydrogen phoionization rate, HI_h=hydrogen photoheating rate +# HeI_i,h=photoionization and photoheating rates for neutral helium +# HeII_i,h=photoionization and photoheating rates for singly ionized helium +# Compton= Compton heating rate. +# All photoionization and photoheating rates are *per ion* and have units of 1/s and eV/s, respectively. +# The Compton heating rate is *per electron* and has units of eV/s. +# z HI_i HI_h HeI_i HeI_h HeII_i HeII_h Compton .000E+00 .228E-13 .889E-13 .124E-13 .112E-12 .555E-15 .114E-13 .690E-19 .491E-01 .284E-13 .111E-12 .157E-13 .140E-12 .676E-15 .138E-13 .827E-19 .101E+00 .354E-13 .139E-12 .196E-13 .174E-12 .823E-15 .168E-13 .985E-19 diff --git a/src/enzo/InitializeHM12Photorates.C b/src/enzo/InitializeHM12Photorates.C new file mode 100644 index 000000000..823e6897b --- /dev/null +++ b/src/enzo/InitializeHM12Photorates.C @@ -0,0 +1,104 @@ +/*********************************************************************** +/ +/ INITIALIZE THE TABULATED HM12 Photo-ionization/heating rates +/ +/ written by: Gabriel Altay +/ date: April, 2013 +/ modified1: +/ +/ PURPOSE: +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + +#include +#include +#include +#include "ErrorExceptions.h" +#include "macros_and_parameters.h" +#include "typedefs.h" +#include "global_data.h" +#include "CosmologyParameters.h" + +int InitializeHM12Photorates() +{ + + + /* Open input file for data. */ + + FILE *fptr = fopen("hm12_photorates.dat", "r"); + if (fptr == NULL) { + ENZO_FAIL("Error opening hm12_photorates.dat\n"); + } + + /* Read rate data, skipping over comments (count first). + Note that photoheating rates are per ion in eV/s */ + + int index = 0; + char line[MAX_LINE_LENGTH]; + double ev2erg; + + ev2erg = 1.60217657e-12; + + RateData.HM12NumberOfRedshiftBins = 0; + + while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) { + if (line[0] != '#') RateData.HM12NumberOfRedshiftBins++; + } + + RateData.HM12Redshifts = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12GH1 = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12GHe1 = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12GHe2 = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12GhH1 = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12GhHe1 = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12GhHe2 = new float[RateData.HM12NumberOfRedshiftBins]; + RateData.HM12Compton = new float[RateData.HM12NumberOfRedshiftBins]; + + rewind(fptr); + while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) { + if (line[0] != '#') + if (sscanf(line, "%"ESYM" %"ESYM" %"ESYM" %"ESYM" %"ESYM" %"ESYM" %"ESYM" %"ESYM, + &RateData.HM12Redshifts[index], + &RateData.HM12GH1[index], + &RateData.HM12GhH1[index], + &RateData.HM12GHe1[index], + &RateData.HM12GhHe1[index], + &RateData.HM12GHe2[index], + &RateData.HM12GhHe2[index], + &RateData.HM12Compton[index]) == 8) { + + RateData.HM12GH1[index] = log10( RateData.HM12GH1[index] ); + RateData.HM12GhH1[index] = log10( RateData.HM12GhH1[index]*ev2erg ); + + RateData.HM12GHe1[index] = log10( RateData.HM12GHe1[index] ); + RateData.HM12GhHe1[index] = log10( RateData.HM12GhHe1[index]*ev2erg ); + + RateData.HM12GHe2[index] = log10( RateData.HM12GHe2[index] ); + RateData.HM12GhHe2[index] = log10( RateData.HM12GhHe2[index]*ev2erg ); + + RateData.HM12Compton[index] = log10( RateData.HM12Compton[index]*ev2erg ); + + index++; + + } + } + + fclose(fptr); + + RateData.HM12RedshiftLo = RateData.HM12Redshifts[0]; + RateData.HM12RedshiftHi = RateData.HM12Redshifts[RateData.HM12NumberOfRedshiftBins-1]; + + if (debug) { + printf("InitializeHM12Photorates: HM12NumberOfRedshiftBins = %"ISYM"\n", + RateData.HM12NumberOfRedshiftBins); + printf("InitializeHM12Photorates: HM12RedshiftLo = %"ESYM"\n", + RateData.HM12RedshiftLo); + printf("InitializeHM12Photorates: HM12RedshiftHi = %"ESYM"\n", + RateData.HM12RedshiftHi); + } + + + return SUCCESS; +} diff --git a/src/enzo/InitializeRateData.C b/src/enzo/InitializeRateData.C index 93da75ef2..d7b880948 100644 --- a/src/enzo/InitializeRateData.C +++ b/src/enzo/InitializeRateData.C @@ -31,6 +31,7 @@ int GetUnits(float *DensityUnits, float *LengthUnits, int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt); int InitializeCloudyCooling(FLOAT Time); int InitializeLymanWernerTable(); +int InitializeHM12Photorates(); int ReadMetalCoolingRates(float TemperatureUnits, float LengthUnits, float aUnits, float DensityUnits, float TimeUnits, float aye); @@ -76,6 +77,10 @@ int InitializeRateData(FLOAT Time) if (debug) printf("InitializeRateData: NumberOfTemperatureBins = %"ISYM"\n", CoolData.NumberOfTemperatureBins); + if (debug) printf("InitializeRateData: RadiationFieldType = %"ISYM"\n", + RadiationFieldType); + + /* Allocate CoolData space for rates. */ CoolData.ceHI = new float[CoolData.NumberOfTemperatureBins]; @@ -220,6 +225,13 @@ int InitializeRateData(FLOAT Time) } } + /* If using tabulated HM12 photo-ionization/heating rates, initialize. */ + if (RadiationFieldType == 15) { + if (InitializeHM12Photorates() == FAIL) { + ENZO_FAIL("Error in InitializeHM12Photorates.\n"); + } + } + /* Initialize Cloudy cooling, even if not being used. */ /* If not used, this will just initialize some data structues. */ if (InitializeCloudyCooling(Time) == FAIL) { diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects index 538ea8b11..781f4edd0 100644 --- a/src/enzo/Make.config.objects +++ b/src/enzo/Make.config.objects @@ -608,6 +608,7 @@ OBJS_CONFIG_LIB = \ InitializeGadgetEquilibriumCoolData.o \ InitializeLocal.o \ InitializeLymanWernerTable.o \ + InitializeHM12Photorates.o \ InitializeMovieFile.o \ InitializeNew.o \ InitializePythonInterface.o \ diff --git a/src/enzo/Make.mach.honest-mistake b/src/enzo/Make.mach.honest-mistake new file mode 100644 index 000000000..2b14c599c --- /dev/null +++ b/src/enzo/Make.mach.honest-mistake @@ -0,0 +1,105 @@ +#======================================================================= +# +# FILE: Make.mach.honest-mistake +# +# DESCRIPTION: Makefile settings for my laptop +# +# AUTHOR: Gabriel Altay (gabriel.altay@gmail.com) +# +# DATE: 2013-04-11 +# +# This configuration assumes that build-essentials, gfortran, +# OpenMPI and HDF5 have been installed using apt-get. +# +#======================================================================= + +MACH_TEXT = Generic Ubuntu +MACH_VALID = 1 +MACH_FILE = Make.mach.honest-mistake + +#----------------------------------------------------------------------- +# Install paths (local variables) +#----------------------------------------------------------------------- + +# LOCAL_HDF5_INSTALL = /usr/local/hdf5/1.8.2s +LOCAL_HDF5_INSTALL = + +#----------------------------------------------------------------------- +# Compiler settings +#----------------------------------------------------------------------- + +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_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 + +# Without MPI + +MACH_CC_NOMPI = gcc # C compiler when not using MPI +MACH_CXX_NOMPI = g++ # C++ compiler when not using MPI +MACH_FC_NOMPI = gfortran # Fortran 77 compiler when not using MPI +MACH_F90_NOMPI = gfortran # Fortran 90 compiler when not using MPI +MACH_LD_NOMPI = g++ # Linker when not using MPI + +#----------------------------------------------------------------------- +# Machine-dependent defines +#----------------------------------------------------------------------- + +MACH_DEFINES = -DLINUX -DH5_USE_16_API + +#----------------------------------------------------------------------- +# Compiler flag settings +#----------------------------------------------------------------------- + + +MACH_CPPFLAGS = -P -traditional +MACH_CFLAGS = +MACH_CXXFLAGS = +MACH_FFLAGS = -fno-second-underscore -ffixed-line-length-132 +MACH_F90FLAGS = -fno-second-underscore +MACH_LDFLAGS = + +#----------------------------------------------------------------------- +# Optimization flags +#----------------------------------------------------------------------- + +MACH_OPT_WARN = -Wall -g +MACH_OPT_DEBUG = -g +MACH_OPT_HIGH = -O2 +MACH_OPT_AGGRESSIVE = -O3 -g + +#----------------------------------------------------------------------- +# Includes +#----------------------------------------------------------------------- + +LOCAL_INCLUDES_MPI = # MPI includes +LOCAL_INCLUDES_HDF5 = -I$(LOCAL_HDF5_INSTALL)/include # HDF5 includes +LOCAL_INCLUDES_HYPRE = # hypre includes +LOCAL_INCLUDES_PAPI = # PAPI includes + +MACH_INCLUDES = $(LOCAL_INCLUDES_HDF5) + +MACH_INCLUDES_MPI = $(LOCAL_INCLUDES_MPI) +MACH_INCLUDES_HYPRE = $(LOCAL_INCLUDES_HYPRE) +MACH_INCLUDES_PAPI = $(LOCAL_INCLUDES_PAPI) + +#----------------------------------------------------------------------- +# Libraries +#----------------------------------------------------------------------- + +LOCAL_LIBS_MPI = # MPI libraries +LOCAL_LIBS_HDF5 = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz # HDF5 libraries +LOCAL_LIBS_HYPRE = # hypre libraries +LOCAL_LIBS_PAPI = # PAPI libraries + +LOCAL_LIBS_MACH = -lgfortran # Machine-dependent libraries + +MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) +MACH_LIBS_MPI = $(LOCAL_LIBS_MPI) +MACH_LIBS_HYPRE = $(LOCAL_LIBS_HYPRE) +MACH_LIBS_PAPI = $(LOCAL_LIBS_PAPI) diff --git a/src/enzo/RadiationFieldCalculateRates.C b/src/enzo/RadiationFieldCalculateRates.C index 95b384848..3877efc94 100644 --- a/src/enzo/RadiationFieldCalculateRates.C +++ b/src/enzo/RadiationFieldCalculateRates.C @@ -87,7 +87,6 @@ int RadiationFieldCalculateRates(FLOAT Time) if (ComovingCoordinates) { CosmologyComputeExpansionFactor(Time, &a, &dadt); - aUnits = 1.0/(1.0 + InitialRedshift); Redshift = 1.0/(a*aUnits) - 1; } else { @@ -97,6 +96,7 @@ int RadiationFieldCalculateRates(FLOAT Time) CoolData.RadiationRedshiftFullOn = RadiationFieldRedshift+0.1; CoolData.RadiationRedshiftDropOff = 0.0; } + double tbase1 = TimeUnits; double xbase1 = LengthUnits/(a*aUnits); @@ -509,8 +509,89 @@ int RadiationFieldCalculateRates(FLOAT Time) RateData.k31 = 1.13e8 * CoolData.f3 * TimeUnits; } + + /* ------------------------------------------------------------------ */ + /* 15) tabulated photoionization and photoheating rates from + Haardt and Madau 2012. rates are read into memory in + InitializeHM12Photorates.C and interpolated here. */ + + if (RadiationFieldType == 15) { + + double zhi, delta, frac; + int ilo, ihi; + + + /* first find indices that bracket the input redshift */ + if ( Redshift <= RateData.HM12RedshiftLo ) { + ilo = 0; + ihi = ilo + 1; + frac = 0.0; + } + else if ( Redshift >= RateData.HM12RedshiftHi ) { + ihi = RateData.HM12NumberOfRedshiftBins - 1; + ilo = ihi - 1; + frac = 1.0; + } + else { + ihi = 0; + zhi = RateData.HM12Redshifts[ihi]; + while ( zhi < Redshift ) { + ihi = ihi + 1; + zhi = RateData.HM12Redshifts[ihi]; + } + ilo = ihi-1; + delta = RateData.HM12Redshifts[ihi] - RateData.HM12Redshifts[ilo]; + frac = ( Redshift - RateData.HM12Redshifts[ilo] ) / delta; + } + + /* now interpolate the rates */ + delta = RateData.HM12GH1[ihi] - RateData.HM12GH1[ilo]; + RateData.k24 = POW( 10.0, RateData.HM12GH1[ilo] + frac * delta ); + RateData.k24 *= TimeUnits * Ramp; + + delta = RateData.HM12GHe1[ihi] - RateData.HM12GHe1[ilo]; + RateData.k25 = POW( 10.0, RateData.HM12GHe1[ilo] + frac * delta ); + RateData.k25 *= TimeUnits * Ramp; + + delta = RateData.HM12GHe2[ihi] - RateData.HM12GHe2[ilo]; + RateData.k26 = POW( 10.0, RateData.HM12GHe2[ilo] + frac * delta ); + RateData.k26 *= TimeUnits * Ramp; + + delta = RateData.HM12GhH1[ihi] - RateData.HM12GhH1[ilo]; + CoolData.piHI = POW( 10.0, RateData.HM12GhH1[ilo] + frac * delta ); + CoolData.piHI = CoolData.piHI / CoolingUnits * Ramp; + + delta = RateData.HM12GhHe1[ihi] - RateData.HM12GhHe1[ilo]; + CoolData.piHeI = POW( 10.0, RateData.HM12GhHe1[ilo] + frac * delta ); + CoolData.piHeI = CoolData.piHeI / CoolingUnits * Ramp; + + delta = RateData.HM12GhHe2[ihi] - RateData.HM12GhHe2[ilo]; + CoolData.piHeII = POW( 10.0, RateData.HM12GhHe2[ilo] + frac * delta ); + CoolData.piHeII = CoolData.piHeII / CoolingUnits * Ramp; + + /* + printf( "Altay: Redshift = %"ESYM"\n", Redshift); + printf( "Altay: RadiationRedshiftOn = %"ESYM"\n", CoolData.RadiationRedshiftOn); + printf( "Altay: RadiationRedshiftFullOn = %"ESYM"\n", CoolData.RadiationRedshiftFullOn); + printf( "Altay: Ramp = %"ESYM"\n", Ramp); + + printf( "Altay: TimeUnits = %"ESYM"\n", TimeUnits); + printf( "Altay: CoolingUnits = %"ESYM"\n", CoolingUnits); + printf( "Altay: k24 = %"ESYM"\n",RateData.k24); + printf( "Altay: k25 = %"ESYM"\n",RateData.k25); + printf( "Altay: k26 = %"ESYM"\n",RateData.k26); + printf( "Altay: piHI = %"ESYM"\n",CoolData.piHI); + printf( "Altay: piHeI = %"ESYM"\n",CoolData.piHeI); + printf( "Altay: piHeII = %"ESYM"\n",CoolData.piHeII); + + ENZO_FAIL("Altay: In Radiation Field Calculate \n"); + */ + + } + + /* ------------------------------------------------------------------ */ - if (RadiationFieldType < 0 || RadiationFieldType > 14) { + if (RadiationFieldType < 0 || RadiationFieldType > 15) { ENZO_VFAIL("RadiationFieldType %"ISYM" not recognized.\n", RadiationFieldType) } @@ -598,5 +679,18 @@ int RadiationFieldCalculateRates(FLOAT Time) fclose(fp); */ + if (MyProcessorNumber == ROOT_PROCESSOR) { + + printf( "Altay: TimeUnits = %"ESYM"\n", TimeUnits); + printf( "Altay: CoolingUnits = %"ESYM"\n", CoolingUnits); + printf( "Altay: k24 = %"ESYM"\n",RateData.k24); + printf( "Altay: k25 = %"ESYM"\n",RateData.k25); + printf( "Altay: k26 = %"ESYM"\n",RateData.k26); + printf( "Altay: piHI = %"ESYM"\n",CoolData.piHI); + printf( "Altay: piHeI = %"ESYM"\n",CoolData.piHeI); + printf( "Altay: piHeII = %"ESYM"\n",CoolData.piHeII); + + } + return SUCCESS; } diff --git a/src/enzo/RateData.h b/src/enzo/RateData.h index 7d8491e55..28a11b392 100644 --- a/src/enzo/RateData.h +++ b/src/enzo/RateData.h @@ -79,4 +79,28 @@ struct RateDataType float *n_cr_n; float *n_cr_d1; float *n_cr_d2; + + + + /* Haardt and Madau 2012 photorates */ + int HM12NumberOfRedshiftBins; + float HM12RedshiftLo; + float HM12RedshiftHi; + + float *HM12Redshifts; + + /* hydrogen / helium photo-ionization rates (functions of z) */ + float *HM12GH1; + float *HM12GHe1; + float *HM12GHe2; + + /* hydrogen / helium photo-heating rates (functions of z) */ + float *HM12GhH1; + float *HM12GhHe1; + float *HM12GhHe2; + + /* Compton heating rate */ + float *HM12Compton; + + };