From 9a0cc75502be3963a897fff637a78ae5f3c24462 Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Fri, 11 Sep 2009 11:06:00 -0400 Subject: [PATCH 01/19] Merged with tip --HG-- branch : week-of-code --- src/enzo/Make.mach.scinet | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/src/enzo/Make.mach.scinet b/src/enzo/Make.mach.scinet index a8b181461..604e32697 100644 --- a/src/enzo/Make.mach.scinet +++ b/src/enzo/Make.mach.scinet @@ -23,31 +23,38 @@ MACH_FILE = Make.mach.scinet # Install paths (local variables) #----------------------------------------------------------------------- +#LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-intel-v11.0-ofed LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-gcc-v4.4.0-ofed LOCAL_HDF5_INSTALL = /scinet/gpc/Libraries/HDF5-1.8.3-v18-openmpi -LOCAL_PYTHON_INSTALL = /usr/bin/python +LOCAL_PYTHON_INSTALL = /scinet/gpc/tools/Python/Python262/ #----------------------------------------------------------------------- # Compiler settings #----------------------------------------------------------------------- -MACH_CPP = /usr/bin/cpp +MACH_CPP = mpiCC # With MPI -MACH_CC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpicc -MACH_CXX_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC -MACH_FC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif77 -MACH_F90_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif90 -MACH_LD_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC +MACH_CC_MPI = mpicc +MACH_CXX_MPI = mpiCC +MACH_FC_MPI = mpif77 +MACH_F90_MPI = mpif90 +MACH_LD_MPI = mpiCC +#MACH_CC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpicc +#MACH_CXX_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC +#MACH_FC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif77 +#MACH_F90_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif90 +#MACH_LD_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC + # Without MPI -MACH_CC_NOMPI = gcc -MACH_CXX_NOMPI = gcc -MACH_FC_NOMPI = gfortran -MACH_F90_NOMPI = gfortran -MACH_LD_NOMPI = gcc +MACH_CC_NOMPI = icc +MACH_CXX_NOMPI = icc +MACH_FC_NOMPI = ifort +MACH_F90_NOMPI = ifort +MACH_LD_NOMPI = icc #----------------------------------------------------------------------- # Machine-dependent defines From 009f3eacae1d296f68e2a46f8a310acb25a3981e Mon Sep 17 00:00:00 2001 From: Fen Zhao Date: Wed, 28 Oct 2009 20:14:11 -0700 Subject: [PATCH 02/19] Added 2D Shearing test problems. Added vorticity output for enzo writing. --HG-- branch : week-of-code --- src/enzo/Grid.h | 8 + src/enzo/Grid_Group_WriteGrid.C | 191 +++++++++++- src/enzo/Grid_ShearingBox2DInitializeGrid.C | 294 ++++++++++++++++++ src/enzo/Grid_ShearingBoxInitializeGrid.C | 2 + ...Grid_ShearingBoxStratifiedInitializeGrid.C | 213 +++++++++++++ src/enzo/Grid_WriteGrid.C | 220 ++++++++++++- src/enzo/InitializeNew.C | 11 +- src/enzo/InterpretCommandLine.C | 6 +- src/enzo/Make.config.objects | 4 + src/enzo/Make.mach.orange | 3 +- src/enzo/Makefile | 2 +- src/enzo/SetDefaultGlobalValues.C | 2 +- src/enzo/ShearingBox2DInitialize.C | 176 +++++++++++ src/enzo/enzo.C | 89 +++++- src/enzo/global_data.h | 5 + src/enzo/hydro_rk/Grid_MHDSourceTerms.C | 30 +- src/enzo/hydro_rk/Grid_SourceTerms.C | 23 +- 17 files changed, 1238 insertions(+), 41 deletions(-) create mode 100644 src/enzo/Grid_ShearingBox2DInitializeGrid.C create mode 100644 src/enzo/Grid_ShearingBoxStratifiedInitializeGrid.C create mode 100644 src/enzo/ShearingBox2DInitialize.C diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 1cc69184d..da8238edc 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -1736,6 +1736,14 @@ int CollapseTestInitializeGrid(int NumberOfSpheres, int ShearingBoxInitializeGrid(float ThermalMagneticRatio, float fraction, float ShearingGeometry, int InitialMagneticFieldConfiguration); + + int ShearingBox2DInitializeGrid(float ThermalMagneticRatio, float fraction, + float ShearingGeometry, + int InitialMagneticFieldConfiguration); + + int ShearingBoxStratifiedInitializeGrid(float ThermalMagneticRatio, float fraction, + float ShearingGeometry, + int InitialMagneticFieldConfiguration); // ------------------------------------------------------------------------- // Analysis functions for AnalysisBaseClass and it's derivatives. // diff --git a/src/enzo/Grid_Group_WriteGrid.C b/src/enzo/Grid_Group_WriteGrid.C index 147aa9132..bb249331b 100644 --- a/src/enzo/Grid_Group_WriteGrid.C +++ b/src/enzo/Grid_Group_WriteGrid.C @@ -300,7 +300,7 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f /* 2c) Loop over fields, writing each one. */ for (field = 0; field < NumberOfBaryonFields; field++) { - + /* copy active part of field into grid */ for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) @@ -332,6 +332,7 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f DataUnits[field] = "none"; } + printf("OutPut Field2: %d \n", field); WriteStringAttr(dset_id, "Label", DataLabel[field], log_fptr); WriteStringAttr(dset_id, "Units", DataUnits[field], log_fptr); WriteStringAttr(dset_id, "Format", "e10.4", log_fptr); @@ -415,6 +416,194 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f } // end: if (OutputTemperature) + + + + + if (VelAnyl) { + int Vel1Num, Vel2Num, Vel3Num; + + for (int i=0; i< NumberOfBaryonFields; i++){ + switch (FieldType[i]){ + case Velocity1: + Vel1Num=i; break; + case Velocity2: + Vel2Num=i; break; + case Velocity3: + Vel3Num=i; break; + } + } + + + io_type *curlx; io_type *curly; + + + if(GridRank==3){ + curlx = new io_type [size]; + curly = new io_type [size];} + + + io_type *curlz = new io_type [size]; + io_type *div = new io_type [size]; + + FLOAT dx = CellWidth[0][0], + dy = CellWidth[1][0], dz; + if (GridRank>2) + dz = CellWidth[2][0]; + + /* Copy active part of field into grid */ + int igrid, igridyp1, igridym1, igridzp1, igridzm1; + for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) { + for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { + for (int i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) { + + igrid = (k*GridDimension[1] + j)*GridDimension[0] + i; + igridyp1 = (k*GridDimension[1] + j + 1)*GridDimension[0] + i; + igridym1 = (k*GridDimension[1] + j - 1)*GridDimension[0] + i; + igridzp1 = ((k+1)*GridDimension[1]+j)*GridDimension[0] + i; + igridzm1 = ((k-1)*GridDimension[1]+j)*GridDimension[0] + i; + + + + if (GridRank==3){ + + div[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx + + 0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy + + 0.5*(BaryonField[Vel3Num][igridzp1]-BaryonField[Vel3Num][igridzm1])/dz) + ); + curlx[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel3Num][igridyp1]-BaryonField[Vel3Num][igridym1])/dy - + 0.5*(BaryonField[Vel2Num][igridzp1]-BaryonField[Vel2Num][igridzm1])/dz) + ); + curly[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel1Num][igridzp1]-BaryonField[Vel1Num][igridzm1])/dz - + 0.5*(BaryonField[Vel3Num][igrid+1]-BaryonField[Vel3Num][igrid-1])/dx) + ); + curlz[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx - + 0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy) + ); + } + + if (GridRank==2){ + + div[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx + + 0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy + )); + curlz[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx - + 0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy) + ); + } + + + + } + } + } + + /* output */ + + char *DataLabelN[4]; + if (GridRank==2) { + DataLabelN[0]="Velocity_Div"; + DataLabelN[1]="Velocity_Vorticity"; + } + if (GridRank==3) { + DataLabelN[0]="Velocity_Div"; + DataLabelN[1]="Velocity_Vorticity1"; + DataLabelN[2]="Velocity_Vorticity2"; + DataLabelN[3]="Velocity_Vorticity3"; + } + + + + int tFields=4; + if (GridRank==2) tFields=2; + + for (int field=0; field<=tFields; field++){ + + file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL); + if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id); + if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);} + + if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n",DataLabelN[field]); + + dset_id = H5Dcreate(file_id, DataLabelN[field], file_type_id, file_dsp_id, H5P_DEFAULT); + if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id); + if( dset_id == h5_error ){my_exit(EXIT_FAILURE);} + + /* set datafield name and units, etc. */ + + if ( DataUnits[field] == NULL ) + { + DataUnits[field] = "none"; + } + + + WriteStringAttr(dset_id, "Label", DataLabelN[field], log_fptr); + WriteStringAttr(dset_id, "Units", NULL, log_fptr); + WriteStringAttr(dset_id, "Format", "e10.4", log_fptr); + WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr); + + switch (field) { + case 0: + temp=div; + case 1: + temp=curlz; + case 2: + temp=curly; + case 3: + temp=curlx; + + } + + h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp); + if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + + h5_status = H5Sclose(file_dsp_id); + if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + + h5_status = H5Dclose(dset_id); + if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + + + } + + + + + if(GridRank==3){ + delete curlx; + delete curly;} + delete curlz; + delete div; + } + + if (OutputCoolingTime) { /* Allocate field and compute cooling time. */ diff --git a/src/enzo/Grid_ShearingBox2DInitializeGrid.C b/src/enzo/Grid_ShearingBox2DInitializeGrid.C new file mode 100644 index 000000000..d7af4c16b --- /dev/null +++ b/src/enzo/Grid_ShearingBox2DInitializeGrid.C @@ -0,0 +1,294 @@ +/*********************************************************************** +/ +/ INITIALIZE A SHEARING BOX TEST ON A GRID +/ +/ written by: Fen Zhao +/ date: June, 2009 +/ modified1: +/ +/ PURPOSE: +/ Set up exither an advecting sphere or the standard shearing box simluation +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + +#include +#include +#include +#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 "CosmologyParameters.h" +#include "hydro_rk/EOS.h" + +int GetUnits(float *DensityUnits, float *LengthUnits, + float *TemperatureUnits, float *TimeUnits, + float *VelocityUnits, FLOAT Time); +double Gaussian(double cs); + +int grid::ShearingBox2DInitializeGrid(float ThermalMagneticRatio, float fraction, float ShearingGeometry, int InitialMagneticFieldConfiguration) +{ + + // 1 to 5 :2D shearing box + // 1 = hydro shearing box + // 2 = KH instablity (v_y= Omega*sin(2.*3.1415927*ShearingGeometry*y)*fraction) + // 3 = sheared sphere (radius= ShearingGeometry) + // 4 = rotating sphere, uniform density (radius= ShearingGeometry) + // 5 = rotating sphere, ramped density (radius= ShearingGeometry) + + + float Pi= 3.14159; + + + + /* declarations */ + + int phip_num; + NumberOfBaryonFields = 0; + FieldType[iden=NumberOfBaryonFields++] = Density; + FieldType[ivx=NumberOfBaryonFields++] = Velocity1; + FieldType[ivy=NumberOfBaryonFields++] = Velocity2; + + FieldType[ietot=NumberOfBaryonFields++] = TotalEnergy; + if (DualEnergyFormalism) { + FieldType[ieint=NumberOfBaryonFields++] = InternalEnergy; + } + if (HydroMethod == MHD_RK) { + FieldType[iBx=NumberOfBaryonFields++] = Bfield1; + FieldType[iBy=NumberOfBaryonFields++] = Bfield2; + + FieldType[NumberOfBaryonFields++] = PhiField; + } + + if(UseDivergenceCleaning) { + FieldType[NumberOfBaryonFields++] = Phi_pField; + } + + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) { + return SUCCESS; + } + + int size = 1; + for (int dim = 0; dim < GridRank; dim++) { + size *= GridDimension[dim]; + } + + for (int field = 0; field < NumberOfBaryonFields; field++) { + if (BaryonField[field] == NULL) { + BaryonField[field] = new float[size]; + for (int i = 0; i < size; i++) + BaryonField[field][i] = 0.0; + } + } + + float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1.0, + TimeUnits = 1.0, VelocityUnits = 1.0; + if (UsePhysicalUnit) + GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, Time); + float MagneticUnits = sqrt(4.0*Pi*DensityUnits)*VelocityUnits; + + /* Problem parameters */ + float rho = 1.0; + float cs = 1e-3; + float pres = rho*cs*cs; + float B0 = 0.0; + float c_s=1e-3; + + float rhou, lenu, tempu, tu, velu; + GetUnits(&rhou, &lenu, &tempu, &tu, &velu, Time); + float PressureUnits = rhou*pow(velu,2); + float bunit=sqrt(4.0*3.14159*rhou*velu*velu); + + const float q = VelocityGradient; + const float Omega = AngularVelocity; + + + + + const FLOAT Lx = DomainRightEdge[0] - DomainLeftEdge[0]; + const FLOAT Ly = DomainRightEdge[1] - DomainLeftEdge[1]; + /* Set up the background */ + + int n = 0; + for (int k = 0; k < GridDimension[2]; k++) { + for (int j = 0; j < GridDimension[1]; j++) { + for (int i = 0; i < GridDimension[0]; i++, n++) { + + FLOAT x; + if (ShearingBoxProblemType==1) + x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; + else + x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]-Lx/2; + + FLOAT y; + if (ShearingBoxProblemType==3 || ShearingBoxProblemType==4 || + ShearingBoxProblemType==5) + y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j] -Ly/2; + else + y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j]; + + + + BaryonField[iden][n] = rho; + + + float vx = 0.0; + float vy = -x*q*Omega; + float vz = 0.0; + + BaryonField[ivx][n] = vx; + BaryonField[ivy][n] = vy; + BaryonField[ivz][n] = vz; + + + + if (ShearingBoxProblemType==2){ //KH instability + float deltaRho=1.0; + FLOAT x1 = 0.3*Lx/2.0; + FLOAT delx = x1 * 0.1; + FLOAT ramp = 1./((1.+exp(-2/delx*(x+x1)))*(1.+exp(-2/delx*(x1-x)))); + float drho=deltaRho*ramp; + + BaryonField[ivy][n]+=(x*q*Omega)*ramp; + BaryonField[ivx][n]+=Omega*sin(2.*3.1415927*ShearingGeometry*y)*fraction; + BaryonField[iden][n]+= drho; + + } + + + if (ShearingBoxProblemType==3){ //sheared sphere + float radius= ShearingGeometry; + float deltaRho=1.0; + FLOAT delx = radius * 0.1; + float r=sqrt(x*x+y*y); + FLOAT ramp = 1./(1.+exp(-2/delx*(radius-r))); + float drho=deltaRho*ramp; + + BaryonField[iden][n] += drho; + + } + + + if (ShearingBoxProblemType==4){ //rotating sphere, no density difference + float radius= ShearingGeometry; + //float deltaRho=1.0; + FLOAT delx = radius * 0.1; + float r=sqrt(x*x+y*y); + float theta=atan2(y,x); + //if (theta<0) theta=theta+3.1415927*2; + + float vx=r*sin(theta); + float vy=r*cos(theta); + + FLOAT ramp = 1./(1.+exp(-2/delx*(radius-r))); + float dvx=vx*ramp; + float dvy=vy*ramp; + + BaryonField[ivx][n] += dvx; + BaryonField[ivy][n] += dvy; + + + } + + if (ShearingBoxProblemType==5){ //rotating sphere, density difference + FLOAT radius=(FLOAT) ShearingGeometry; + float deltaRho=1.0; + FLOAT delx = radius * 0.1; + + float r=sqrt(x*x+y*y); + float theta=atan(y/x); + + float vx=r*sin(theta); + float vy=r*cos(theta); + + FLOAT ramp = 1./(1.+exp(-2/delx*(radius-r))); + float dvx=vx*ramp; + float dvy=vy*ramp; + float drho=deltaRho*ramp; + + BaryonField[ivx][n] += dvx; + BaryonField[ivy][n] += dvy; + BaryonField[iden][n] += drho; + + + } + + + vx= BaryonField[ivx][n]; + vy =BaryonField[ivy][n]; + vz= BaryonField[ivz][n]; + + float eint, h, cs_temp, dpdrho, dpde; + EOS(pres, BaryonField[iden][n], eint, h, cs_temp, dpdrho, dpde, EOSType, 1); + //EOS(pres, rho, eint, h, cs_temp, dpdrho, dpde, EOSType, 1); + + BaryonField[ietot][n] = eint + 0.5*(vx*vx + vy*vy + vz*vz); + + if (DualEnergyFormalism) { + BaryonField[ieint][n] = eint; + } + + if (HydroMethod == MHD_RK) { + + float rhoActual=BaryonField[iden ][n]; + float pressure=c_s*c_s*rhoActual/Gamma; + float realpressure=pressure*PressureUnits; + float InitialBField=sqrt((8*3.14159*realpressure/(ThermalMagneticRatio)))/bunit; + if (InitialMagneticFieldConfiguration == 0) B0 = InitialBField; + else if (InitialMagneticFieldConfiguration == 1) B0 = InitialBField*sin(2*3.14159*x); + + + BaryonField[iBz][n] = B0; + BaryonField[ietot][n] += 0.5 * pow(B0,2) / rho; + } + + } // end loop over grid + } + } + + + + + if (ShearingBoxProblemType == 1){ + /* ProblemType 1: Vortex wave. + Reference: B. M. Johnson & C. F. Gammie, ApJ, 2005, 626, 978. */ + + const FLOAT kx0 = (-8.0*2.0*Pi/Lx)/(8.0); + const FLOAT ky = 2.0*2.0*Pi/Ly; + const float vx0 = 1e-4; // in unit of cs + n = 0; + for (int k = 0; k < GridDimension[2]; k++) { + for (int j = 0; j < GridDimension[1]; j++) { + for (int i = 0; i < GridDimension[0]; i++, n++) { + int igrid = i + (j + k * GridDimension[1]) * GridDimension[0]; + FLOAT x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; + FLOAT y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j]; + + float dvx = vx0 * cs * cos(kx0 * x + ky * y); + float dvy = -kx0 / ky * dvx; + float drho = rho/(cs*ky) * (-2.0*kx0*kx0*q*Omega + 2.0*(q-1.0)*Omega) * vx0 * sin(kx0*x + ky*y); + + //BaryonField[iden][igrid] += drho; + + BaryonField[ivx ][igrid] += dvx; + BaryonField[ivy ][igrid] += dvy; + + } + }}} + + + + //PrintToScreenBoundaries(BaryonField[ivy], "Vy (Initial)\n" ,2, 0); + // PrintToScreenBoundaries(BaryonField[iden], "Dens (Initial)\n", 2, 0); + + + return SUCCESS; + +} diff --git a/src/enzo/Grid_ShearingBoxInitializeGrid.C b/src/enzo/Grid_ShearingBoxInitializeGrid.C index 4c7d4c166..1da55e403 100644 --- a/src/enzo/Grid_ShearingBoxInitializeGrid.C +++ b/src/enzo/Grid_ShearingBoxInitializeGrid.C @@ -34,6 +34,8 @@ double Gaussian(double cs); int grid::ShearingBoxInitializeGrid(float ThermalMagneticRatio, float fraction, float ShearingGeometry, int InitialMagneticFieldConfiguration) { + // + /* declarations */ diff --git a/src/enzo/Grid_ShearingBoxStratifiedInitializeGrid.C b/src/enzo/Grid_ShearingBoxStratifiedInitializeGrid.C new file mode 100644 index 000000000..2ba46996c --- /dev/null +++ b/src/enzo/Grid_ShearingBoxStratifiedInitializeGrid.C @@ -0,0 +1,213 @@ +/*********************************************************************** +/ +/ INITIALIZE A SHEARING BOX TEST ON A GRID +/ +/ written by: Fen Zhao +/ date: June, 2009 +/ modified1: +/ +/ PURPOSE: +/ Set up exither an advecting sphere or the standard shearing box simluation +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + +#include +#include +#include +#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 "CosmologyParameters.h" +#include "hydro_rk/EOS.h" + +int GetUnits(float *DensityUnits, float *LengthUnits, + float *TemperatureUnits, float *TimeUnits, + float *VelocityUnits, FLOAT Time); +double Gaussian(double cs); + +int grid::ShearingBoxStratifiedInitializeGrid(float ThermalMagneticRatio, float fraction, float ShearingGeometry, int InitialMagneticFieldConfiguration) +{ + + + /* declarations */ + + + + + int phip_num; + NumberOfBaryonFields = 0; + FieldType[iden=NumberOfBaryonFields++] = Density; + FieldType[ivx=NumberOfBaryonFields++] = Velocity1; + FieldType[ivy=NumberOfBaryonFields++] = Velocity2; + FieldType[ivz=NumberOfBaryonFields++] = Velocity3; + + FieldType[ietot=NumberOfBaryonFields++] = TotalEnergy; + if (DualEnergyFormalism) { + FieldType[ieint=NumberOfBaryonFields++] = InternalEnergy; + } + if (useMHD){ + FieldType[iBx=NumberOfBaryonFields++] = Bfield1; + FieldType[iBy=NumberOfBaryonFields++] = Bfield2; + FieldType[iBz=NumberOfBaryonFields++] = Bfield3; + FieldType[NumberOfBaryonFields++] = PhiField; + } + + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) { + return SUCCESS; + } + + + + + + +// int iBx, iBy, iBz; + int iB[3]={-1,-1,-1}; + if (useMHD){ + iB[0]=iBx; + iB[1]=iBy; + if (GridRank==3){ + iB[2]=iBz; + } + } + + int size = 1; + for (int dim = 0; dim < GridRank; dim++) { + size *= GridDimension[dim]; + } + + for (int field = 0; field < NumberOfBaryonFields; field++) { + if (BaryonField[field] == NULL) { + BaryonField[field] = new float[size]; + } + } + + srand(110182*ProcessorNumber); + + + int i,j,k; + int n=0; + float eint, v2, vx,vy,vz; + + float rhou, lenu, tempu, tu, velu; + GetUnits(&rhou, &lenu, &tempu, &tu, &velu, Time); + float PressureUnits = rhou*pow(velu,2); + + + float magnitude=AngularVelocity*fraction; + float rho=1.0; + + float c_s=1e-3; + + float lengthx=DomainRightEdge[0]-DomainLeftEdge[0]; + float lengthy=DomainRightEdge[1]-DomainLeftEdge[1]; + float lengthz=DomainRightEdge[2]-DomainLeftEdge[2]; + + float h, cs, dpdrho, dpde, H, pressure; + float bunit=sqrt(4.0*3.14159*rhou*velu*velu); + + FLOAT x,y,z; + + + for (k = 0; k < GridDimension[2]; k++) { + for (j = 0; j < GridDimension[1]; j++) { + for (i = 0; i < GridDimension[0]; i++, n++) { + + FLOAT xPos[3] = {CellLeftEdge[0][i] + 0.5*CellWidth[0][i]-lengthx/2.0, + CellLeftEdge[1][j] + 0.5*CellWidth[1][j]-lengthy/2.0, + 0.0}; + + + x=xPos[0]; y=xPos[1]; z=0.0; + + if (GridRank==3) { + xPos[2]=CellLeftEdge[2][k] + 0.5*CellWidth[2][k]-lengthz/2.0; + z=xPos[2]; + } + + + float xVel[3] ={0,0,0}; + + if (ShearingBoxProblemType == 0){ + if (x*x+y*y+z*z<0.25*ShearingGeometry*ShearingGeometry){ + BaryonField[iden ][n]=500.0*rho; + + } + else if (x*x+y*y+z*z2) + dz = CellWidth[2][0]; + + /* Copy active part of field into grid */ + int igrid, igridyp1, igridym1, igridzp1, igridzm1; + for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) { + for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { + for (int i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) { + + igrid = (k*GridDimension[1] + j)*GridDimension[0] + i; + igridyp1 = (k*GridDimension[1] + j + 1)*GridDimension[0] + i; + igridym1 = (k*GridDimension[1] + j - 1)*GridDimension[0] + i; + igridzp1 = ((k+1)*GridDimension[1]+j)*GridDimension[0] + i; + igridzm1 = ((k-1)*GridDimension[1]+j)*GridDimension[0] + i; + + + + if (GridRank==3){ + + div[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx + + 0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy + + 0.5*(BaryonField[Vel3Num][igridzp1]-BaryonField[Vel3Num][igridzm1])/dz) + ); + curlx[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel3Num][igridyp1]-BaryonField[Vel3Num][igridym1])/dy - + 0.5*(BaryonField[Vel2Num][igridzp1]-BaryonField[Vel2Num][igridzm1])/dz) + ); + curly[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel1Num][igridzp1]-BaryonField[Vel1Num][igridzm1])/dz - + 0.5*(BaryonField[Vel3Num][igrid+1]-BaryonField[Vel3Num][igrid-1])/dx) + ); + curlz[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx - + 0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy) + ); + } + + if (GridRank==2){ + + div[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx + + 0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy + )); + curlz[(i-GridStartIndex[0])+ + (j-GridStartIndex[1])*ActiveDim[0]+ + (k-GridStartIndex[2])*ActiveDim[0]*ActiveDim[1] ] = + io_type( + (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx - + 0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy) + ); + } + + + + } + } + } + + /* output */ + + char *DataLabelN[4]; + if (GridRank==2) { + DataLabelN[0]="Velocity_Div"; + DataLabelN[1]="Velocity_Vorticity"; + } + if (GridRank==3) { + DataLabelN[0]="Velocity_Div"; + DataLabelN[1]="Velocity_Vorticity1"; + DataLabelN[2]="Velocity_Vorticity2"; + DataLabelN[3]="Velocity_Vorticity3"; + } + + + + int tFields=3; + if (GridRank==2) tFields=1; + fprintf(stderr,"Inside 0 \n"); + + for (int field=0; field<=tFields; field++){ + fprintf(stderr,"Inside 1 %d \n", field); + file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL); + if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id); + if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);} + + if (io_log) fprintf(log_fptr,"H5Dcreate with Name = %s\n",DataLabelN[field]); + + dset_id = H5Dcreate(file_id, DataLabelN[field], file_type_id, file_dsp_id, H5P_DEFAULT); + if (io_log) fprintf(log_fptr, "H5Dcreate id: %"ISYM"\n", dset_id); + if( dset_id == h5_error ){my_exit(EXIT_FAILURE);} + + /* set datafield name and units, etc. */ + + if ( DataUnits[field] == NULL ) + { + DataUnits[field] = "none"; + } + + fprintf(stderr,"Inside 2 %d \n", field); + fprintf(stderr, DataLabelN[field]); + + WriteStringAttr(dset_id, "Label", DataLabelN[field], log_fptr); + fprintf(stderr,"Inside 3 %d \n", field); + WriteStringAttr(dset_id, "Units", "VortUnits", log_fptr); + fprintf(stderr,"Inside 4 %d \n", field); + WriteStringAttr(dset_id, "Format", "e10.4", log_fptr); + fprintf(stderr,"Inside 5 %d \n", field); + WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr); + fprintf(stderr,"Inside 6 %d \n", field); + + switch (field) { + case 0: + temp=div; + case 1: + temp=curlz; + case 2: + temp=curly; + case 3: + temp=curlx; + + } + + h5_status = H5Dwrite(dset_id, float_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) temp); + if (io_log) fprintf(log_fptr, "H5Dwrite: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + + h5_status = H5Sclose(file_dsp_id); + if (io_log) fprintf(log_fptr, "H5Sclose: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + + h5_status = H5Dclose(dset_id); + if (io_log) fprintf(log_fptr, "H5Dclose: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + + + } + + fprintf(stderr,"Inside sorta \n"); + + + if(GridRank==3){ + delete curlx; + delete curly;} + delete curlz; + delete div; + } + + fprintf(stderr,"Outside 10 \n"); + if (OutputCoolingTime) { /* Allocate field and compute cooling time. */ @@ -454,6 +654,7 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) delete [] cooling_time; } // if (OutputCoolingTime) + fprintf(stderr,"Outside 11 \n"); /* Make sure that there is a copy of dark matter field to save (and at the right resolution). */ @@ -853,6 +1054,7 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) } + /* clean up */ delete [] temp; @@ -860,22 +1062,31 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) } // end: if (MyProcessorNumber...) } // end: if (NumberOfParticles > 0) + + fprintf(stderr,"Outside 11.5 \n"); /* Close HDF file. */ if (MyProcessorNumber == ProcessorNumber) { + h5_status = H5Fclose(file_id); - if (io_log) fprintf(log_fptr, "H5Fclose: %"ISYM"\n", h5_status); - if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + fprintf(stderr,"Outside 11.51 %d %d\n", h5_status, h5_error); + if (io_log) fprintf(log_fptr, "H5Fclose: %"ISYM"\n", h5_status); + if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} + fprintf(stderr,"Outside 11.52 \n"); } if (MyProcessorNumber == ProcessorNumber) - { + { fprintf(stderr,"Outside 11.53 \n"); if (io_log) fclose(log_fptr); + fprintf(stderr,"Outside 11.54 \n"); } /* 4) Save Gravity info. */ + + fprintf(stderr,"Outside 11.6 \n"); + if (MyProcessorNumber == ROOT_PROCESSOR) if (SelfGravity) @@ -883,7 +1094,10 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) /* Clean up. */ + fprintf(stderr,"Outside 12 \n"); delete [] name; + + fprintf(stderr,"Outside 13 \n"); return SUCCESS; diff --git a/src/enzo/InitializeNew.C b/src/enzo/InitializeNew.C index 8d7832425..4750edf81 100644 --- a/src/enzo/InitializeNew.C +++ b/src/enzo/InitializeNew.C @@ -126,7 +126,10 @@ int TracerParticleCreation(FILE *fptr, HierarchyEntry &TopGrid, int ShearingBoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData); - +int ShearingBox2DInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, + TopGridData &MetaData); +int ShearingBoxStratifiedInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, + TopGridData &MetaData); #ifdef TRANSFER int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData); @@ -435,7 +438,11 @@ int InitializeNew(char *filename, HierarchyEntry &TopGrid, // 35) Shearing Box Simulation if (ProblemType == 35) ret = ShearingBoxInitialize(fptr, Outfptr, TopGrid, MetaData); - + if (ProblemType == 36) + ret = ShearingBox2DInitialize(fptr, Outfptr, TopGrid, MetaData); + if (ProblemType == 37) + ret = ShearingBoxStratifiedInitialize(fptr, Outfptr, TopGrid, MetaData); + // 40) Supernova Explosion from restart diff --git a/src/enzo/InterpretCommandLine.C b/src/enzo/InterpretCommandLine.C index 6cf3467cd..9d05e7aa4 100644 --- a/src/enzo/InterpretCommandLine.C +++ b/src/enzo/InterpretCommandLine.C @@ -28,7 +28,7 @@ int InterpretCommandLine(int argc, char *argv[], char *myname, int &InformationOutput, int &OutputAsParticleData, int &project, int &ProjectionDimension, - int &ProjectionSmooth, + int &ProjectionSmooth, int &velanyl, char *ParameterFile[], int RegionStart[], int RegionEnd[], FLOAT RegionStartCoordinate[], @@ -224,6 +224,10 @@ int InterpretCommandLine(int argc, char *argv[], char *myname, extract = TRUE; break; + case 'v': + velanyl = TRUE; + break; + /* Unknown */ default: diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects index 476384124..eb8dcf3a0 100644 --- a/src/enzo/Make.config.objects +++ b/src/enzo/Make.config.objects @@ -385,6 +385,8 @@ OBJS_CONFIG_LIB = \ Grid_SetMinimumSupport.o \ Grid_SetParticleMassFlaggingField.o \ Grid_ShearingBoxInitializeGrid.o \ + Grid_ShearingBox2DInitializeGrid.o \ + Grid_ShearingBoxStratifiedInitializeGrid.o \ Grid_ShockTubeInitializeGrid.o \ Grid_SolveForPotential.o \ Grid_SolveHydroEquations.o \ @@ -553,6 +555,8 @@ OBJS_CONFIG_LIB = \ SetLevelTimeStep.o \ sgi_st1_fft64.o \ ShearingBoxInitialize.o \ + ShearingBox2DInitialize.o \ + ShearingBoxStratifiedInitialize.o \ ShockInABoxInitialize.o \ ShockPoolInitialize.o \ ShockTubeInitialize.o \ diff --git a/src/enzo/Make.mach.orange b/src/enzo/Make.mach.orange index 41ff4f4db..bc42c3d5a 100644 --- a/src/enzo/Make.mach.orange +++ b/src/enzo/Make.mach.orange @@ -58,7 +58,8 @@ MACH_LD_NOMPI = icc #MACH_DEFINES = -DXT3 -DNO_IO_LOG -DSYSCALL -DENZO_ANALYSIS #MACH_DEFINES = -DXT3 -DNO_IO_LOG -DSYSCALL -DHAVE_SPRNG #MACH_DEFINES = -DXT3 -DNO_IO_LOG -DSYSCALL -DSFGEN_PERF -DHAVE_SPRNG -DUSE_STOCHASTIC_FORCING -MACH_DEFINES = -DLINUX -DH5_USE_16_API +MACH_DEFINES = -DLINUX -DH5_USE_16_API -DIO_LOG + #----------------------------------------------------------------------- # Compiler flag settings diff --git a/src/enzo/Makefile b/src/enzo/Makefile index 90a9a2f3d..a37e9344a 100644 --- a/src/enzo/Makefile +++ b/src/enzo/Makefile @@ -46,7 +46,7 @@ ENZO_DIR = . MODULES = # please remove this -SVN = /u/ki/mturk/Research/local/yt-x86_64/bin/hg +SVN = hg #----------------------------------------------------------------------- # Make.config.settings is used for setting default values to all compile-time diff --git a/src/enzo/SetDefaultGlobalValues.C b/src/enzo/SetDefaultGlobalValues.C index e233013cd..8efdb62b5 100644 --- a/src/enzo/SetDefaultGlobalValues.C +++ b/src/enzo/SetDefaultGlobalValues.C @@ -574,6 +574,6 @@ int SetDefaultGlobalValues(TopGridData &MetaData) useMHD=0; MoveParticlesBetweenSiblings = FALSE; - + VelAnyl = FALSE; return SUCCESS; } diff --git a/src/enzo/ShearingBox2DInitialize.C b/src/enzo/ShearingBox2DInitialize.C new file mode 100644 index 000000000..21f25fa09 --- /dev/null +++ b/src/enzo/ShearingBox2DInitialize.C @@ -0,0 +1,176 @@ +/*********************************************************************** +/ +/ INITIALIZE A SHEARING BOX TEST +/ +/ written by: Fen Zhao +/ date: June, 2009 +/ modified1: +/ +/ PURPOSE: +/ Set up exither an advecting sphere or the standard shearing box simluation +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + + + +#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 "Hierarchy.h" +#include "LevelHierarchy.h" +#include "TopGridData.h" + +// void WriteListOfFloats(FILE *fptr, int N, float floats[]); +// void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]); +// void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); +// int RebuildHierarchy(TopGridData *MetaData, +// LevelHierarchyEntry *LevelArray[], int level); +// void WriteListOfFloats(FILE *fptr, int N, float floats[]); +// void WriteListOfFloats(FILE *fptr, int N, EFLOAT floats[]); +// void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); +// int RebuildHierarchy(TopGridData *MetaData, +// LevelHierarchyEntry *LevelArray[], int level); +// int GetUnits(float *DensityUnits, float *LengthUnits, +// float *TemperatureUnits, float *TimeUnits, +// float *VelocityUnits, EFLOAT Time); + +void WriteListOfFloats(FILE *fptr, int N, float floats[]); +void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]); +void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); +int RebuildHierarchy(TopGridData *MetaData, + LevelHierarchyEntry *LevelArray[], int level); + + +int ShearingBox2DInitialize (FILE *fptr, FILE *Outfptr, + HierarchyEntry &TopGrid, TopGridData &MetaData) + +{ + + //Lets Assume the shear is in the y direction and the x direction is radial + + const char *DensName = "Density"; + 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 *BxName = "Bx"; + const char *ByName = "By"; + // const char *BzName = "Bz"; + const char *PhiName = "Phi"; + const char *DebugName = "Debug"; + const char *Phi_pName = "Phip"; + + /* declarations */ + + char line[MAX_LINE_LENGTH]; + + + /* set default parameters */ + + + /* read input from file */ + + + float ThermalMagneticRatio=400; + float FluctuationAmplitudeFraction=0.1; + int ShearingBoxRefineAtStart = FALSE; + float ShearingGeometry=2.0; + int InitialMagneticFieldConfiguration=0; + int RefineAtStart=1; + + int ret; + + + + while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) { + + ret = 0; + +/* read parameters */ + ret += sscanf(line, "ShearingBoxRefineAtStart = %"ISYM, &RefineAtStart); + ret += sscanf(line, "ShearingBoxThermalMagneticRatio= %"FSYM, &ThermalMagneticRatio); + ret += sscanf(line, "ShearingBoxFluctuationAmplitudeFraction = %"FSYM, &FluctuationAmplitudeFraction); + ret += sscanf(line, "ShearingBoxGeometry = %"FSYM, &ShearingGeometry); + ret += sscanf(line, "ShearingBoxInitialMagneticFieldConfiguration = %"ISYM, &InitialMagneticFieldConfiguration); + + } + + + + if (TopGrid.GridData->ShearingBox2DInitializeGrid(ThermalMagneticRatio, FluctuationAmplitudeFraction, ShearingGeometry, InitialMagneticFieldConfiguration) + == FAIL) + ENZO_FAIL("Error in ShearingBoxInitializeGrid.\n"); + + + + if (RefineAtStart) { + /* Declare, initialize and fill out the LevelArray. */ + + LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY]; + for (int level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) + LevelArray[level] = NULL; + AddLevel(LevelArray, &TopGrid, 0); + + /* Add levels to the maximum depth or until no new levels are created, + and re-initialize the level after it is created. */ + + for (int level = 0; level < MaximumRefinementLevel; level++) { + if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) { + fprintf(stderr, "Error in RebuildHierarchy.\n"); + ENZO_FAIL(""); + } + if (LevelArray[level+1] == NULL) + break; + LevelHierarchyEntry *Temp = LevelArray[level+1]; + while (Temp != NULL) { + if (Temp->GridData->ShearingBox2DInitializeGrid(ThermalMagneticRatio, FluctuationAmplitudeFraction, ShearingGeometry, InitialMagneticFieldConfiguration) + == FAIL) + ENZO_FAIL("Error in ShearingBoxInitializeGrid.\n"); + Temp = Temp->NextGridThisLevel; + } + } // end: loop over levels + } // end: if (CollapseTestRefineAtStart) + + + + + /* set up field names and units */ + + int count = 0; + DataLabel[count++] = (char*) DensName; + DataLabel[count++] = (char*) Vel1Name; + DataLabel[count++] = (char*) Vel2Name; + //DataLabel[count++] = (char*) Vel3Name; + DataLabel[count++] = (char*) TEName; + if (DualEnergyFormalism) { + DataLabel[count++] = (char*) GEName; + } + if(useMHD){ + DataLabel[count++] = (char*) BxName; + DataLabel[count++] = (char*) ByName; + // DataLabel[count++] = (char*) BzName; + DataLabel[count++] = (char*) PhiName; + } + + for (int i = 0; i < count; i++) { + DataUnits[i] = NULL; + } + + + + + return SUCCESS; + +} diff --git a/src/enzo/enzo.C b/src/enzo/enzo.C index a74448df3..a6cbfe777 100644 --- a/src/enzo/enzo.C +++ b/src/enzo/enzo.C @@ -101,7 +101,7 @@ int InterpretCommandLine(int argc, char *argv[], char *myname, int &InformationOutput, int &OutputAsParticleDataFlag, int &project, int &ProjectionDimension, - int &ProjectionSmooth, + int &ProjectionSmooth, int &velanyl, char *ParameterFile[], int RegionStart[], int RegionEnd[], FLOAT RegionStartCoordinates[], @@ -109,6 +109,21 @@ int InterpretCommandLine(int argc, char *argv[], char *myname, int &Level, int MyProcessorNumber); void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); int SetDefaultGlobalValues(TopGridData &MetaData); + +int WriteAllData(char *basename, int filenumber, HierarchyEntry *TopGrid, + TopGridData &MetaData, ExternalBoundary *Exterior, + FLOAT WriteTime = -1); +int FastSiblingLocatorInitialize(ChainingMeshStructure *Mesh, int Rank, + int TopGridDims[]); +int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh); +int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, + SiblingGridList SiblingList[], + int level, TopGridData *MetaData, + ExternalBoundary *Exterior, LevelHierarchyEntry *Level); + +int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level, + HierarchyEntry **Grids[]); + int CommunicationInitialize(Eint32 *argc, char **argv[]); int CommunicationFinalize(); @@ -243,11 +258,12 @@ Eint32 main(Eint32 argc, char *argv[]) // Initialize int restart = FALSE, - OutputAsParticleDataFlag = FALSE, - InformationOutput = FALSE, - project = FALSE, - ProjectionDimension = INT_UNDEFINED, - ProjectionSmooth = FALSE; + OutputAsParticleDataFlag = FALSE, + InformationOutput = FALSE, + project = FALSE, + ProjectionDimension = INT_UNDEFINED, + ProjectionSmooth = FALSE, + velanyl = FALSE; debug = FALSE; extract = FALSE; flow_trace_on = FALSE; @@ -329,7 +345,7 @@ Eint32 main(Eint32 argc, char *argv[]) if (InterpretCommandLine(int_argc, argv, myname, restart, debug, extract, InformationOutput, OutputAsParticleDataFlag, - project, ProjectionDimension, ProjectionSmooth, + project, ProjectionDimension, ProjectionSmooth, velanyl, &ParameterFile, RegionStart, RegionEnd, RegionStartCoordinates, RegionEndCoordinates, @@ -344,7 +360,7 @@ Eint32 main(Eint32 argc, char *argv[]) // If we need to read the parameter file as a restart file, do it now - if (restart || OutputAsParticleDataFlag || extract || InformationOutput || project) { + if (restart || OutputAsParticleDataFlag || extract || InformationOutput || project || velanyl) { SetDefaultGlobalValues(MetaData); @@ -454,6 +470,63 @@ Eint32 main(Eint32 argc, char *argv[]) } // ENDFOR dim } + + + /* Do vector analysis */ + + if (velanyl) { + VelAnyl = TRUE; + + for (int level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) { + HierarchyEntry **Grids; + int NumberOfGrids = GenerateGridArray(LevelArray, level, &Grids); + if (LevelArray[level] != NULL) { + /* Initialize the chaining mesh used in the FastSiblingLocator. */ + + ChainingMeshStructure ChainingMesh; + FastSiblingLocatorInitialize(&ChainingMesh, MetaData.TopGridRank, + MetaData.TopGridDims); + SiblingGridList *SiblingList = new SiblingGridList[NumberOfGrids]; + if (SiblingList == NULL) { + fprintf(stderr, "Error allocating SiblingList\n"); + return FAIL; + } + + /* Add all the grids to the chaining mesh. */ + + for (int grid1 = 0; grid1 < NumberOfGrids; grid1++) + Grids[grid1]->GridData->FastSiblingLocatorAddGrid(&ChainingMesh); + + /* For each grid, get a list of possible siblings from the chaining mesh. */ + + for (int grid1 = 0; grid1 < NumberOfGrids; grid1++) + if (Grids[grid1]->GridData->FastSiblingLocatorFindSiblings( + &ChainingMesh, &SiblingList[grid1], + MetaData.LeftFaceBoundaryCondition, + MetaData.RightFaceBoundaryCondition) == FAIL) { + fprintf(stderr, "Error in grid->FastSiblingLocatorFindSiblings.\n"); + return FAIL; + } + + /* Clean up the chaining mesh. */ + + FastSiblingLocatorFinalize(&ChainingMesh); + + if (SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, &MetaData, + &Exterior, LevelArray[level]) == FAIL) { + printf("error setboundary"); + } + } + } + if (WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber-1, + &TopGrid, MetaData, &Exterior) == FAIL) { + fprintf(stderr, "Error in WriteAllData.\n"); + return FAIL; + } + my_exit(EXIT_SUCCESS); + } + + // Normal start: Open and read parameter file #ifdef MEM_TRACE diff --git a/src/enzo/global_data.h b/src/enzo/global_data.h index c39f7d302..e1ce8b625 100644 --- a/src/enzo/global_data.h +++ b/src/enzo/global_data.h @@ -704,4 +704,9 @@ EXTERN int RefineByJeansLengthUnits; EXTERN int MoveParticlesBetweenSiblings; + +/* Vorticity Calculations */ + +EXTERN int VelAnyl; + #endif diff --git a/src/enzo/hydro_rk/Grid_MHDSourceTerms.C b/src/enzo/hydro_rk/Grid_MHDSourceTerms.C index 5127cbb6a..d0d14afe6 100644 --- a/src/enzo/hydro_rk/Grid_MHDSourceTerms.C +++ b/src/enzo/hydro_rk/Grid_MHDSourceTerms.C @@ -321,9 +321,10 @@ int grid::MHDSourceTerms(float **dU) /* Add centrifugal force for the shearing box */ - if (ProblemType == 35 && ShearingBoxProblemType !=0) { - - int igrid; + if ((ProblemType == 35 || ProblemType == 36 ||ProblemType == 37) && ShearingBoxProblemType !=0) { + + + int igrid; float rho, gx, gy, gz; FLOAT xPos[3]; float vels[3]; @@ -332,7 +333,8 @@ int grid::MHDSourceTerms(float **dU) int iden=FindField(Density, FieldType, NumberOfBaryonFields); int ivx=FindField(Velocity1, FieldType, NumberOfBaryonFields); int ivy=FindField(Velocity2, FieldType, NumberOfBaryonFields); - int ivz=FindField(Velocity3, FieldType, NumberOfBaryonFields); + int ivz; + if (GridRank==3) ivz=FindField(Velocity3, FieldType, NumberOfBaryonFields); int indexNumbers[3]={iS1,iS2,iS3}; @@ -341,9 +343,10 @@ int grid::MHDSourceTerms(float **dU) float lengthx=DomainRightEdge[0]-DomainLeftEdge[0]; float lengthy=DomainRightEdge[1]-DomainLeftEdge[1]; - float lengthz=DomainRightEdge[2]-DomainLeftEdge[2]; - - + float lengthz; + if (GridRank==3) lengthz=DomainRightEdge[2]-DomainLeftEdge[2]; + else lengthz-0.0; + for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) { for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { @@ -353,19 +356,19 @@ int grid::MHDSourceTerms(float **dU) rho = BaryonField[iden][igrid]; xPos[0] = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]-lengthx/2.0; xPos[1] = CellLeftEdge[1][i] + 0.5*CellWidth[1][i]-lengthy/2.0; - xPos[2] = CellLeftEdge[2][i] + 0.5*CellWidth[2][i]-lengthz/2.0; - + if (GridRank==3) xPos[2] = CellLeftEdge[2][i] + 0.5*CellWidth[2][i]-lengthz/2.0; + else xPos[2]=0; + vels[0] = BaryonField[ivx][igrid]; vels[1] = BaryonField[ivy][igrid]; - vels[2] = BaryonField[ivz][igrid]; - + if (GridRank==3) vels[2] = BaryonField[ivz][igrid]; + else vels[2]=0; //Omega cross V - FLOAT temp= dU[indexNumbers[1]][n]; dU[indexNumbers[0]][n] -= dtFixed*2.0*rho*(A[1]*vels[2]-A[2]*vels[1]); dU[indexNumbers[1]][n] -= dtFixed*2.0*rho*(A[2]*vels[0]-A[0]*vels[2]); - dU[indexNumbers[2]][n] -= dtFixed*2.0*rho*(A[0]*vels[1]-A[1]*vels[0]); + if (GridRank==3) dU[indexNumbers[2]][n] -= dtFixed*2.0*rho*(A[0]*vels[1]-A[1]*vels[0]); dU[indexNumbers[ShearingBoundaryDirection]][n] += dtFixed*2.0*rho*VelocityGradient*AngularVelocity*AngularVelocity*xPos[ShearingBoundaryDirection]; @@ -375,6 +378,7 @@ int grid::MHDSourceTerms(float **dU) + } } diff --git a/src/enzo/hydro_rk/Grid_SourceTerms.C b/src/enzo/hydro_rk/Grid_SourceTerms.C index dd054bb67..9fa8ade94 100644 --- a/src/enzo/hydro_rk/Grid_SourceTerms.C +++ b/src/enzo/hydro_rk/Grid_SourceTerms.C @@ -244,7 +244,7 @@ int grid::SourceTerms(float **dU) /* Add centrifugal force for the shearing box */ - if (ProblemType == 35 && ShearingBoxProblemType !=0) { + if ((ProblemType == 35 || ProblemType == 36 ||ProblemType == 37) && ShearingBoxProblemType !=0) { int igrid; float rho, gx, gy, gz; @@ -255,7 +255,8 @@ int grid::SourceTerms(float **dU) int iden=FindField(Density, FieldType, NumberOfBaryonFields); int ivx=FindField(Velocity1, FieldType, NumberOfBaryonFields); int ivy=FindField(Velocity2, FieldType, NumberOfBaryonFields); - int ivz=FindField(Velocity3, FieldType, NumberOfBaryonFields); + int ivz; + if (GridRank==3) ivz=FindField(Velocity3, FieldType, NumberOfBaryonFields); int indexNumbers[3]={iS1,iS2,iS3}; @@ -264,9 +265,10 @@ int grid::SourceTerms(float **dU) float lengthx=DomainRightEdge[0]-DomainLeftEdge[0]; float lengthy=DomainRightEdge[1]-DomainLeftEdge[1]; - float lengthz=DomainRightEdge[2]-DomainLeftEdge[2]; - - + float lengthz; + if (GridRank==3) lengthz=DomainRightEdge[2]-DomainLeftEdge[2]; + else lengthz-0.0; + for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) { for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { @@ -276,18 +278,19 @@ int grid::SourceTerms(float **dU) rho = BaryonField[iden][igrid]; xPos[0] = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]-lengthx/2.0; xPos[1] = CellLeftEdge[1][i] + 0.5*CellWidth[1][i]-lengthy/2.0; - xPos[2] = CellLeftEdge[2][i] + 0.5*CellWidth[2][i]-lengthz/2.0; - + if (GridRank==3) xPos[2] = CellLeftEdge[2][i] + 0.5*CellWidth[2][i]-lengthz/2.0; + else xPos[2]=0; + vels[0] = BaryonField[ivx][igrid]; vels[1] = BaryonField[ivy][igrid]; - vels[2] = BaryonField[ivz][igrid]; - + if (GridRank==3) vels[2] = BaryonField[ivz][igrid]; + else vels[2]=0; //Omega cross V dU[indexNumbers[0]][n] -= dtFixed*2.0*rho*(A[1]*vels[2]-A[2]*vels[1]); dU[indexNumbers[1]][n] -= dtFixed*2.0*rho*(A[2]*vels[0]-A[0]*vels[2]); - dU[indexNumbers[2]][n] -= dtFixed*2.0*rho*(A[0]*vels[1]-A[1]*vels[0]); + if (GridRank==3) dU[indexNumbers[2]][n] -= dtFixed*2.0*rho*(A[0]*vels[1]-A[1]*vels[0]); dU[indexNumbers[ShearingBoundaryDirection]][n] += dtFixed*2.0*rho*VelocityGradient*AngularVelocity*AngularVelocity*xPos[ShearingBoundaryDirection]; From 63e8004d9a2c53550efbff6e5c073c69754beeba Mon Sep 17 00:00:00 2001 From: Fen Zhao Date: Wed, 28 Oct 2009 23:43:02 -0700 Subject: [PATCH 03/19] stratified template --HG-- branch : week-of-code --- src/enzo/ShearingBoxStratifiedInitialize.C | 173 +++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 src/enzo/ShearingBoxStratifiedInitialize.C diff --git a/src/enzo/ShearingBoxStratifiedInitialize.C b/src/enzo/ShearingBoxStratifiedInitialize.C new file mode 100644 index 000000000..3adffd20b --- /dev/null +++ b/src/enzo/ShearingBoxStratifiedInitialize.C @@ -0,0 +1,173 @@ +/*********************************************************************** +/ +/ INITIALIZE A SHEARING BOX TEST +/ +/ written by: Fen Zhao +/ date: June, 2009 +/ modified1: +/ +/ PURPOSE: +/ Set up exither an advecting sphere or the standard shearing box simluation +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + + + +#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 "Hierarchy.h" +#include "LevelHierarchy.h" +#include "TopGridData.h" + +// void WriteListOfFloats(FILE *fptr, int N, float floats[]); +// void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]); +// void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); +// int RebuildHierarchy(TopGridData *MetaData, +// LevelHierarchyEntry *LevelArray[], int level); +// void WriteListOfFloats(FILE *fptr, int N, float floats[]); +// void WriteListOfFloats(FILE *fptr, int N, EFLOAT floats[]); +// void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); +// int RebuildHierarchy(TopGridData *MetaData, +// LevelHierarchyEntry *LevelArray[], int level); +// int GetUnits(float *DensityUnits, float *LengthUnits, +// float *TemperatureUnits, float *TimeUnits, +// float *VelocityUnits, EFLOAT Time); + +void WriteListOfFloats(FILE *fptr, int N, float floats[]); +void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]); +void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level); +int RebuildHierarchy(TopGridData *MetaData, + LevelHierarchyEntry *LevelArray[], int level); + + +int ShearingBoxStratifiedInitialize (FILE *fptr, FILE *Outfptr, + HierarchyEntry &TopGrid, TopGridData &MetaData) + +{ + const char *DensName = "Density"; + 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 *BxName = "Bx"; + const char *ByName = "By"; + const char *BzName = "Bz"; + const char *PhiName = "Phi"; + const char *DebugName = "Debug"; + const char *Phi_pName = "Phip"; + + /* declarations */ + + char line[MAX_LINE_LENGTH]; + + + /* set default parameters */ + + + /* read input from file */ + + + float ThermalMagneticRatio=400; + float FluctuationAmplitudeFraction=0.1; + int ShearingBoxRefineAtStart = FALSE; + float ShearingGeometry=2.0; + int InitialMagneticFieldConfiguration=0; + int RefineAtStart=1; + + int ret; + + + + while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) { + + ret = 0; + +/* read parameters */ + ret += sscanf(line, "ShearingBoxRefineAtStart = %"ISYM, &RefineAtStart); + ret += sscanf(line, "ShearingBoxThermalMagneticRatio= %"FSYM, &ThermalMagneticRatio); + ret += sscanf(line, "ShearingBoxFluctuationAmplitudeFraction = %"FSYM, &FluctuationAmplitudeFraction); + ret += sscanf(line, "ShearingBoxGeometry = %"FSYM, &ShearingGeometry); + ret += sscanf(line, "ShearingBoxInitialMagneticFieldConfiguration = %"ISYM, &InitialMagneticFieldConfiguration); + + } + + + + if (TopGrid.GridData->ShearingBoxInitializeGrid(ThermalMagneticRatio, FluctuationAmplitudeFraction, ShearingGeometry, InitialMagneticFieldConfiguration) + == FAIL) + ENZO_FAIL("Error in ShearingBoxInitializeGrid.\n"); + + + + if (RefineAtStart) { + /* Declare, initialize and fill out the LevelArray. */ + + LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY]; + for (int level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) + LevelArray[level] = NULL; + AddLevel(LevelArray, &TopGrid, 0); + + /* Add levels to the maximum depth or until no new levels are created, + and re-initialize the level after it is created. */ + + for (int level = 0; level < MaximumRefinementLevel; level++) { + if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) { + fprintf(stderr, "Error in RebuildHierarchy.\n"); + ENZO_FAIL(""); + } + if (LevelArray[level+1] == NULL) + break; + LevelHierarchyEntry *Temp = LevelArray[level+1]; + while (Temp != NULL) { + if (Temp->GridData->ShearingBoxInitializeGrid(ThermalMagneticRatio, FluctuationAmplitudeFraction, ShearingGeometry, InitialMagneticFieldConfiguration) + == FAIL) + ENZO_FAIL("Error in ShearingBoxInitializeGrid.\n"); + Temp = Temp->NextGridThisLevel; + } + } // end: loop over levels + } // end: if (CollapseTestRefineAtStart) + + + + + /* set up field names and units */ + + int count = 0; + DataLabel[count++] = (char*) DensName; + DataLabel[count++] = (char*) Vel1Name; + DataLabel[count++] = (char*) Vel2Name; + DataLabel[count++] = (char*) Vel3Name; + DataLabel[count++] = (char*) TEName; + if (DualEnergyFormalism) { + DataLabel[count++] = (char*) GEName; + } + if(useMHD){ + DataLabel[count++] = (char*) BxName; + DataLabel[count++] = (char*) ByName; + DataLabel[count++] = (char*) BzName; + DataLabel[count++] = (char*) PhiName; + } + + for (int i = 0; i < count; i++) { + DataUnits[i] = NULL; + } + + + + + return SUCCESS; + +} From a30bed46038563790df20e0035e953ea43bcb09e Mon Sep 17 00:00:00 2001 From: Fen Zhao Date: Fri, 30 Oct 2009 22:58:51 -0700 Subject: [PATCH 04/19] added vorticity write to new io routines --HG-- branch : week-of-code --- src/enzo/Grid_WriteGrid.C | 37 ++++++----- src/enzo/New_Grid_WriteGrid.C | 114 ++++++++++++++++++++++++++++++++++ src/enzo/enzo.C | 24 ++++--- 3 files changed, 149 insertions(+), 26 deletions(-) diff --git a/src/enzo/Grid_WriteGrid.C b/src/enzo/Grid_WriteGrid.C index bcfa99717..fffd82aee 100644 --- a/src/enzo/Grid_WriteGrid.C +++ b/src/enzo/Grid_WriteGrid.C @@ -49,6 +49,7 @@ int GetUnits(float *DensityUnits, float *LengthUnits, float *TemperatureUnits, float *TimeUnits, float *VelocityUnits, FLOAT Time); + int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) { @@ -512,10 +513,10 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) int tFields=3; if (GridRank==2) tFields=1; - fprintf(stderr,"Inside 0 \n"); + for (int field=0; field<=tFields; field++){ - fprintf(stderr,"Inside 1 %d \n", field); + file_dsp_id = H5Screate_simple((Eint32) GridRank, OutDims, NULL); if (io_log) fprintf(log_fptr, "H5Screate file_dsp_id: %"ISYM"\n", file_dsp_id); if( file_dsp_id == h5_error ){my_exit(EXIT_FAILURE);} @@ -533,17 +534,17 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) DataUnits[field] = "none"; } - fprintf(stderr,"Inside 2 %d \n", field); + fprintf(stderr, DataLabelN[field]); WriteStringAttr(dset_id, "Label", DataLabelN[field], log_fptr); - fprintf(stderr,"Inside 3 %d \n", field); + WriteStringAttr(dset_id, "Units", "VortUnits", log_fptr); - fprintf(stderr,"Inside 4 %d \n", field); + WriteStringAttr(dset_id, "Format", "e10.4", log_fptr); - fprintf(stderr,"Inside 5 %d \n", field); + WriteStringAttr(dset_id, "Geometry", "Cartesian", log_fptr); - fprintf(stderr,"Inside 6 %d \n", field); + switch (field) { case 0: @@ -582,7 +583,7 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) delete div; } - fprintf(stderr,"Outside 10 \n"); + if (OutputCoolingTime) { @@ -654,7 +655,7 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) delete [] cooling_time; } // if (OutputCoolingTime) - fprintf(stderr,"Outside 11 \n"); + /* Make sure that there is a copy of dark matter field to save (and at the right resolution). */ @@ -1063,7 +1064,7 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) } // end: if (MyProcessorNumber...) } // end: if (NumberOfParticles > 0) - fprintf(stderr,"Outside 11.5 \n"); + /* Close HDF file. */ @@ -1071,21 +1072,21 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) { h5_status = H5Fclose(file_id); - fprintf(stderr,"Outside 11.51 %d %d\n", h5_status, h5_error); + if (io_log) fprintf(log_fptr, "H5Fclose: %"ISYM"\n", h5_status); if( h5_status == h5_error ){my_exit(EXIT_FAILURE);} - fprintf(stderr,"Outside 11.52 \n"); + } if (MyProcessorNumber == ProcessorNumber) - { fprintf(stderr,"Outside 11.53 \n"); + { if (io_log) fclose(log_fptr); - fprintf(stderr,"Outside 11.54 \n"); + } /* 4) Save Gravity info. */ - fprintf(stderr,"Outside 11.6 \n"); + if (MyProcessorNumber == ROOT_PROCESSOR) @@ -1094,11 +1095,13 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) /* Clean up. */ - fprintf(stderr,"Outside 12 \n"); + delete [] name; - fprintf(stderr,"Outside 13 \n"); + return SUCCESS; } + + diff --git a/src/enzo/New_Grid_WriteGrid.C b/src/enzo/New_Grid_WriteGrid.C index f0ebd5b79..61ea7f1d3 100644 --- a/src/enzo/New_Grid_WriteGrid.C +++ b/src/enzo/New_Grid_WriteGrid.C @@ -268,6 +268,120 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f } } + if (VelAnyl){ + + int Vel1Num, Vel2Num, Vel3Num; + Vel1Num=FindField(Velocity1, FieldType, NumberOfBaryonFields); + Vel2Num=FindField(Velocity2, FieldType, NumberOfBaryonFields); + Vel3Num=FindField(Velocity3, FieldType, NumberOfBaryonFields); + + float *curlx; float *curly; + + if(GridRank==3){ + curlx = new float [size]; + curly = new float [size];} + + float *curlz = new float [size]; + float *div = new float [size]; + + FLOAT dx = CellWidth[0][0], + dy = CellWidth[1][0], dz; + if (GridRank>2) + dz = CellWidth[2][0]; + + /* Copy active part of field into grid */ + int igrid, igridyp1, igridym1, igridzp1, igridzm1; + for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) { + for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { + for (int i = GridStartIndex[0]; i <= GridEndIndex[0]; i++) { + + igrid = (k*GridDimension[1] + j)*GridDimension[0] + i; + igridyp1 = (k*GridDimension[1] + j + 1)*GridDimension[0] + i; + igridym1 = (k*GridDimension[1] + j - 1)*GridDimension[0] + i; + igridzp1 = ((k+1)*GridDimension[1]+j)*GridDimension[0] + i; + igridzm1 = ((k-1)*GridDimension[1]+j)*GridDimension[0] + i; + + + + if (GridRank==3){ + + div[igrid ] = + float( + (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx + + 0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy + + 0.5*(BaryonField[Vel3Num][igridzp1]-BaryonField[Vel3Num][igridzm1])/dz) + ); + curlx[igrid ] = + float( + (0.5*(BaryonField[Vel3Num][igridyp1]-BaryonField[Vel3Num][igridym1])/dy - + 0.5*(BaryonField[Vel2Num][igridzp1]-BaryonField[Vel2Num][igridzm1])/dz) + ); + curly[igrid ] = + float( + (0.5*(BaryonField[Vel1Num][igridzp1]-BaryonField[Vel1Num][igridzm1])/dz - + 0.5*(BaryonField[Vel3Num][igrid+1]-BaryonField[Vel3Num][igrid-1])/dx) + ); + curlz[igrid ] = + float( + (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx - + 0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy) + ); + } + + if (GridRank==2){ + + div[igrid ] = + float( + (0.5*(BaryonField[Vel1Num][igrid+1]-BaryonField[Vel1Num][igrid-1])/dx + + 0.5*(BaryonField[Vel2Num][igridyp1]-BaryonField[Vel2Num][igridym1])/dy + )); + + + curlz[igrid] = + float( + (0.5*(BaryonField[Vel2Num][igrid+1]-BaryonField[Vel2Num][igrid-1])/dx - + 0.5*(BaryonField[Vel1Num][igridyp1]-BaryonField[Vel1Num][igridym1])/dy) + ); + } + + + + } + } + } + + + if (GridRank==3){ + + this->write_dataset(GridRank, OutDims, "Velocity_Div", + group_id, file_type_id, (VOIDP) div, TRUE, temp); + this->write_dataset(GridRank, OutDims, "Velocity_Vorticity1", + group_id, file_type_id, (VOIDP) curlx, TRUE, temp); + this->write_dataset(GridRank, OutDims, "Velocity_Vorticity2", + group_id, file_type_id, (VOIDP) curly, TRUE, temp); + this->write_dataset(GridRank, OutDims, "Velocity_Vorticity3", + group_id, file_type_id, (VOIDP) curlz, TRUE, temp); + + } + + if (GridRank==2){ + + this->write_dataset(GridRank, OutDims, "Velocity_Div", + group_id, file_type_id, (VOIDP) div, TRUE, temp); + this->write_dataset(GridRank, OutDims, "Velocity_Vorticity", + group_id, file_type_id, (VOIDP) curlz, TRUE, temp); + + } + + if(GridRank==3){ + delete [] curlx; + delete [] curly;} + delete [] curlz; + delete [] div; + + + } + /* If this is cosmology, compute the temperature field as well since its such a pain to compute after the fact. */ diff --git a/src/enzo/enzo.C b/src/enzo/enzo.C index 5b4ce6a57..1472b2153 100644 --- a/src/enzo/enzo.C +++ b/src/enzo/enzo.C @@ -113,6 +113,12 @@ int SetDefaultGlobalValues(TopGridData &MetaData); int WriteAllData(char *basename, int filenumber, HierarchyEntry *TopGrid, TopGridData &MetaData, ExternalBoundary *Exterior, FLOAT WriteTime = -1); + +int Group_WriteAllData(char *basename, int filenumber, + HierarchyEntry *TopGrid, TopGridData &MetaData, + ExternalBoundary *Exterior, FLOAT WriteTime = -1, + int RestartDump = FALSE); + int FastSiblingLocatorInitialize(ChainingMeshStructure *Mesh, int Rank, int TopGridDims[]); int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh); @@ -192,14 +198,14 @@ Eint32 main(Eint32 argc, char *argv[]) int int_argc; int_argc = argc; - if (MyProcessorNumber == ROOT_PROCESSOR && - ENZO_SVN_REVISION != 0) { - printf("=========================\n"); - printf("Enzo SVN Branch %s\n",ENZO_SVN_BRANCH); - printf("Enzo SVN Revision %s\n",ENZO_SVN_REVISION); - printf("=========================\n"); - fflush(stdout); - } + // if (MyProcessorNumber == ROOT_PROCESSOR && +// ENZO_SVN_REVISION != 0) { +// printf("=========================\n"); +// printf("Enzo SVN Branch %s\n",ENZO_SVN_BRANCH); +// printf("Enzo SVN Revision %s\n",ENZO_SVN_REVISION); +// printf("=========================\n"); +// fflush(stdout); +// } // Performance Monitoring #ifdef USE_MPI @@ -518,7 +524,7 @@ Eint32 main(Eint32 argc, char *argv[]) } } } - if (WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber-1, + if (Group_WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber-1, &TopGrid, MetaData, &Exterior) == FAIL) { fprintf(stderr, "Error in WriteAllData.\n"); return FAIL; From 1e3ea6fbf860c872b304d5d4680525c4ab6775b6 Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Mon, 2 Nov 2009 11:46:43 -0500 Subject: [PATCH 05/19] Merged with Tip --HG-- branch : week-of-code --- src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C | 2 +- src/enzo/star_maker8.C | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C index 2053a1c09..03f2f4d01 100644 --- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C @@ -584,7 +584,7 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL if (PutSink == 1 && level == 0) { // set it up on level zero and make it mustrefine // double mass_p = 20.0*1.989e33; - double mass_p = 0.01*1.989e33; + double mass_p = 10.0*1.989e33; mass_p /= MassUnits; double dx = CellWidth[0][0]; double den_p = mass_p / pow(dx,3); diff --git a/src/enzo/star_maker8.C b/src/enzo/star_maker8.C index dd60fc0b7..cafd69105 100644 --- a/src/enzo/star_maker8.C +++ b/src/enzo/star_maker8.C @@ -84,6 +84,8 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float printf("Star Maker 8 running - SinkMergeDistance = %g\n", SinkMergeDistance); printf("Star Maker 8: massthresh=%g, jlrefine=%g\n", *massthresh,*jlrefine); + printf("Star Maker 8: time = %g\n", *t); + /* Compute Units. */ From 70572146d475a8b6287ee6cde962e000aeaa5890 Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Mon, 2 Nov 2009 16:42:54 -0500 Subject: [PATCH 06/19] changes to SciNet makefile --HG-- branch : week-of-code --- src/enzo/Make.mach.scinet | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/enzo/Make.mach.scinet b/src/enzo/Make.mach.scinet index 604e32697..b7c8ce0c1 100644 --- a/src/enzo/Make.mach.scinet +++ b/src/enzo/Make.mach.scinet @@ -23,8 +23,8 @@ MACH_FILE = Make.mach.scinet # Install paths (local variables) #----------------------------------------------------------------------- -#LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-intel-v11.0-ofed -LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-gcc-v4.4.0-ofed +LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-intel-v11.0-ofed +#LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-gcc-v4.4.0-ofed LOCAL_HDF5_INSTALL = /scinet/gpc/Libraries/HDF5-1.8.3-v18-openmpi LOCAL_PYTHON_INSTALL = /scinet/gpc/tools/Python/Python262/ From a5bcd85fa2ff0680e0d2139d27c55b33bafb7315 Mon Sep 17 00:00:00 2001 From: Fen Zhao Date: Mon, 2 Nov 2009 22:10:38 -0800 Subject: [PATCH 07/19] changed the loop for SBC --HG-- branch : week-of-code --- src/enzo/SetBoundaryConditions.C | 15 +++++++++----- src/enzo/enzo.C | 35 +++++++++++++++++++++++++++----- 2 files changed, 40 insertions(+), 10 deletions(-) diff --git a/src/enzo/SetBoundaryConditions.C b/src/enzo/SetBoundaryConditions.C index 7e8400ee7..d811eb90d 100644 --- a/src/enzo/SetBoundaryConditions.C +++ b/src/enzo/SetBoundaryConditions.C @@ -66,6 +66,7 @@ int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, int loopEnd = (ShearingBoundaryDirection != -1) ? 2 : 1; + int grid1, grid2, StartGrid, EndGrid, loop; LCAPERF_START("SetBoundaryConditions"); @@ -76,8 +77,9 @@ int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, CommunicationBarrier(); #endif - TIME_MSG("Interpolating boundaries from parent"); - + if (loop == 0) { + TIME_MSG("Interpolating boundaries from parent"); + for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { if (traceMPI) fprintf(tracePtr, "SBC loop\n"); @@ -88,6 +90,8 @@ int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, /* Here, we just generate the calls to generate the receive buffers, without actually doing anything. */ + + CommunicationDirection = COMMUNICATION_POST_RECEIVE; CommunicationReceiveIndex = 0; @@ -98,16 +102,17 @@ int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; - if (loop == 0) { + if (level == 0) { Grids[grid1]->GridData->SetExternalBoundaryValues(Exterior); } else { Grids[grid1]->GridData->InterpolateBoundaryFromParent (Grids[grid1]->ParentGrid->GridData); } - } // ENDIF loop == 0 + } // ENDFOR grids + /* -------------- SECOND PASS ----------------- */ /* Now we generate all the sends, and do all the computation @@ -136,7 +141,7 @@ int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, if (CommunicationReceiveHandler() == FAIL) ENZO_FAIL(""); - } // ENDFOR grid batches + }} // ENDFOR grid batches TIME_MSG("Copying zones in SetBoundaryConditions"); for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { diff --git a/src/enzo/enzo.C b/src/enzo/enzo.C index 1472b2153..56f025223 100644 --- a/src/enzo/enzo.C +++ b/src/enzo/enzo.C @@ -119,13 +119,35 @@ int Group_WriteAllData(char *basename, int filenumber, ExternalBoundary *Exterior, FLOAT WriteTime = -1, int RestartDump = FALSE); -int FastSiblingLocatorInitialize(ChainingMeshStructure *Mesh, int Rank, - int TopGridDims[]); -int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh); + + +#ifdef FAST_SIB int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, SiblingGridList SiblingList[], - int level, TopGridData *MetaData, + int level, TopGridData *MetaData, + ExternalBoundary *Exterior, LevelHierarchyEntry *Level); +#else +int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, + int level, TopGridData *MetaData, ExternalBoundary *Exterior, LevelHierarchyEntry *Level); +#endif + +int FastSiblingLocatorInitialize(ChainingMeshStructure *Mesh, int Rank, + int TopGridDims[]); +int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh); + +int RebuildHierarchy(TopGridData *MetaData, + LevelHierarchyEntry *LevelArray[], int level); +int CopyOverlappingZones(grid* CurrentGrid, TopGridData *MetaData, + LevelHierarchyEntry *LevelArray[], int level); +int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL, + int NumberOfSubgrids[] = NULL, + int FluxFlag = FALSE, + TopGridData* MetaData = NULL); +int CreateSiblingList(HierarchyEntry ** Grids, int NumberOfGrids, SiblingGridList *SiblingList, + int StaticLevelZero,TopGridData * MetaData,int level); + + int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level, HierarchyEntry **Grids[]); @@ -482,6 +504,8 @@ Eint32 main(Eint32 argc, char *argv[]) if (velanyl) { VelAnyl = TRUE; + RebuildHierarchy(&MetaData, LevelArray, 0); + for (int level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) { HierarchyEntry **Grids; @@ -518,12 +542,13 @@ Eint32 main(Eint32 argc, char *argv[]) FastSiblingLocatorFinalize(&ChainingMesh); - if (SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, &MetaData, + if (SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, &MetaData, &Exterior, LevelArray[level]) == FAIL) { printf("error setboundary"); } } } + if (Group_WriteAllData(MetaData.DataDumpName, MetaData.DataDumpNumber-1, &TopGrid, MetaData, &Exterior) == FAIL) { fprintf(stderr, "Error in WriteAllData.\n"); From 6557bb0ae56571800c49ad0fc34c741fab241fdd Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Wed, 11 Nov 2009 13:09:59 -0500 Subject: [PATCH 08/19] debugging isotropic stellar wind (in progress) --HG-- branch : week-of-code --- src/enzo/Grid_StarParticleHandler.C | 6 ++-- src/enzo/ReadParameterFile.C | 1 + src/enzo/SetDefaultGlobalValues.C | 1 + src/enzo/WriteParameterFile.C | 2 ++ src/enzo/global_data.h | 1 + .../hydro_rk/Grid_TurbulenceInitializeGrid.C | 10 +++--- src/enzo/star_maker8.C | 33 +++++++++++++++---- 7 files changed, 41 insertions(+), 13 deletions(-) diff --git a/src/enzo/Grid_StarParticleHandler.C b/src/enzo/Grid_StarParticleHandler.C index 017d795e2..07fe0e618 100644 --- a/src/enzo/Grid_StarParticleHandler.C +++ b/src/enzo/Grid_StarParticleHandler.C @@ -293,6 +293,8 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) if (!StarParticleCreation && !StarParticleFeedback) return SUCCESS; + //printf("\n XXXX StarParticleHandler Called XXXX \n \n"); + if (MyProcessorNumber != ProcessorNumber) return SUCCESS; @@ -446,7 +448,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) h2field[index] = BaryonField[H2INum][index] + BaryonField[H2IINum][index]; } } - + printf("Star type \n"); /* Set the units. */ float DensityUnits = 1, LengthUnits = 1, TemperatureUnits = 1, @@ -847,7 +849,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) ENZO_FAIL("Error in star_maker8.\n"); } } else { - printf("Grid_StarParticleHandler 784 - sink maker called\n"); + printf("Grid_StarParticleHandler 784 - sink maker called (NOT STAR_MAKER8)\n"); if (sink_maker(GridDimension, GridDimension+1, GridDimension+2, &size, BaryonField[DensNum], BaryonField[Vel1Num], BaryonField[Vel2Num], BaryonField[Vel3Num], diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C index 6546259dd..cc5deaae4 100644 --- a/src/enzo/ReadParameterFile.C +++ b/src/enzo/ReadParameterFile.C @@ -714,6 +714,7 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt) ret += sscanf(line, "SinkMergeMass = %"FSYM, &SinkMergeMass); ret += sscanf(line, "StellarWindFeedback = %"ISYM, &StellarWindFeedback); ret += sscanf(line, "StellarWindTurnOnMass = %"FSYM, &StellarWindTurnOnMass); + ret += sscanf(line, "MSStellarWindTurnOnMass = %"FSYM, &MSStellarWindTurnOnMass); // ret += sscanf(line, "VelAnyl = %"ISYM, &VelAnyl); diff --git a/src/enzo/SetDefaultGlobalValues.C b/src/enzo/SetDefaultGlobalValues.C index 3d2abc2d6..2236d7cc6 100644 --- a/src/enzo/SetDefaultGlobalValues.C +++ b/src/enzo/SetDefaultGlobalValues.C @@ -403,6 +403,7 @@ int SetDefaultGlobalValues(TopGridData &MetaData) TotalSinkMass = 0.0; StellarWindFeedback = 0; StellarWindTurnOnMass = 0.1; + MSStellarWindTurnOnMass = 10.0; UseHydro = 1; Coordinate = Cartesian; diff --git a/src/enzo/WriteParameterFile.C b/src/enzo/WriteParameterFile.C index 69478c39b..b486e1ff2 100644 --- a/src/enzo/WriteParameterFile.C +++ b/src/enzo/WriteParameterFile.C @@ -503,6 +503,8 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData) StellarWindFeedback); fprintf(fptr, "StellarWindTurnOnMass = %"FSYM"\n", StellarWindTurnOnMass); + fprintf(fptr, "MSStellarWindTurnOnMass = %"FSYM"\n", + MSStellarWindTurnOnMass); fprintf(fptr, "StarMakerOverDensityThreshold = %"GSYM"\n", diff --git a/src/enzo/global_data.h b/src/enzo/global_data.h index 0d63417bd..88a135ee7 100644 --- a/src/enzo/global_data.h +++ b/src/enzo/global_data.h @@ -623,6 +623,7 @@ EXTERN float SinkMergeMass; EXTERN float TotalSinkMass; EXTERN int StellarWindFeedback; EXTERN float StellarWindTurnOnMass; +EXTERN float MSStellarWindTurnOnMass; EXTERN int NBodyDirectSummation; /* Turbulence simulation parameters */ diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C index 03f2f4d01..d5b5b1dda 100644 --- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C @@ -584,7 +584,7 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL if (PutSink == 1 && level == 0) { // set it up on level zero and make it mustrefine // double mass_p = 20.0*1.989e33; - double mass_p = 10.0*1.989e33; + double mass_p = 30.0*1.989e33; mass_p /= MassUnits; double dx = CellWidth[0][0]; double den_p = mass_p / pow(dx,3); @@ -604,15 +604,15 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL ParticleMass[0] = den_p; ParticleNumber[0] = 0; ParticleType[0] = PARTICLE_TYPE_MUST_REFINE; - ParticlePosition[0][0] = 0.5;//+0.5*dx; - ParticlePosition[1][0] = 0.5;//+0.5*dx; - ParticlePosition[2][0] = 0.5;//+0.5*dx; + ParticlePosition[0][0] = 0.5+0.5*dx; + ParticlePosition[1][0] = 0.5+0.5*dx; + ParticlePosition[2][0] = 0.5+0.5*dx; ParticleVelocity[0][0] = 0.0; ParticleVelocity[1][0] = 0.0; ParticleVelocity[2][0] = 0.0; ParticleAttribute[0][0] = 0.0; // creation time - ParticleAttribute[1][0] = 0; + ParticleAttribute[1][0] = 0.0; ParticleAttribute[2][0] = mass_p; if (StellarWindFeedback) { diff --git a/src/enzo/star_maker8.C b/src/enzo/star_maker8.C index cafd69105..68b0b7c38 100644 --- a/src/enzo/star_maker8.C +++ b/src/enzo/star_maker8.C @@ -11,6 +11,7 @@ SW = 1 - magnetic field protostellar jets SW = 2 - random direction protostellar jets SW = 3 - isotropic main sequence stellar wind + SW = 4 - protostellar and main sequence winds with an accretion disc (in progress) INPUTS: d - density field u,v,w - velocity fields @@ -83,7 +84,7 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float double Pi = 3.1415926; printf("Star Maker 8 running - SinkMergeDistance = %g\n", SinkMergeDistance); - printf("Star Maker 8: massthresh=%g, jlrefine=%g\n", *massthresh,*jlrefine); + // printf("Star Maker 8: massthresh=%g, jlrefine=%g\n", *massthresh,*jlrefine); printf("Star Maker 8: time = %g\n", *t); @@ -547,18 +548,30 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float /* StellarWindFeedback 3: Isotropic wind */ float mdot_wind = 1e-5*(*dt)*(*t1)/(3.1557e7*umass); /* 10^-5 solar mases per year - this is in code units: density x length^3*/ - v_wind = 2.0e6/(*v1); - mdot_wind = mdot_wind/(4.0*Pi); /* mass Per solid angle */ //printf("Adding Stellar wind 3: dt =%e, mdot =%e, Vwind =%e, rho_wind =%e \n",dt,mdot_wind*umass/(*t1),v_wind*(*v1),rho_wind*(*d1)); FLOAT radius_cell[MAX_SUPERCELL_NUMBER]; FLOAT radius2_cell[MAX_SUPERCELL_NUMBER]; float SolidAngle; + FLOAT mdot_wind1, mdot_wind2; if (StellarWindFeedback == 3) { - printf("mdotwind = %e\n",mdot_wind); + //printf("mdotwind = %e\n",mdot_wind); // printf("STELLAR WIND FEEDBACK = 3\n"); for (n = 0; n < nsinks; n++) { bb = sink_index[n]; + + if (mpold[bb]*pow(*dx,3)*umass < MSStellarWindTurnOnMass) { + continue; + } + + v_wind = 1.e5*(-355.554+892.32*log10(mpold[bb]*pow(*dx,3)*umass - 5.24765))/(*v1); + mdot_wind1 = (pow(10,-9.47)*pow(mpold[bb]*pow(*dx,3)*umass,2.2427))*((*t1)*(*dx)/v_wind)/(3.1557e7*umass); + mdot_wind2 = (pow(10,-9.47)*pow(mpold[bb]*pow(*dx,3)*umass,2.2427))*(*dt)*(*t1)/(3.1557e7*umass); + mdot_wind = max(mdot_wind1,mdot_wind2); + //(pow(10,-9.47)*pow(mpold[bb]*pow(*dx,3)*umass,2.2427)) - mass loss in Msun/yr (from N. Smith 2006 table 1) + mdot_wind = mdot_wind/(4.0*Pi);/* mass Per solid angle */ + + printf("\n Mass star = %e Msun, Mdot_wind = %e Msun/yr, and v_wind = %e cm/s \n \n",mpold[bb]*pow(*dx,3)*umass,(4.0*Pi)*mdot_wind/((*dt)*(*t1)/(3.1557e7*umass)), v_wind*(*v1)); i = (xpold[bb] - *xstart)/(*dx); j = (ypold[bb] - *ystart)/(*dx); @@ -576,7 +589,6 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float m_cell = 0; /* Caclualte the total mass of the supercell */ - for (int kk = -3; kk <= 3; kk++) { for (int jj = -3; jj <= 3; jj++) { for (int ii = -3; ii <= 3; ii++) { @@ -610,6 +622,8 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float /* Do the feedback */ + + //printf("n_cell = %i \n", n_cell); float m_wind = 0.0; float cells_volume = 0.0; for (int ic = 0; ic < n_cell; ic++) { @@ -631,7 +645,7 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float else if (radius2_cell[ic] == 27.0) SolidAngle = 0.0302870901; else { SolidAngle = 4.*3.1415926/n_cell; - //printf("star_maker8.C line 373: Radius squared is wrong?!? radius =%f, n_cell = %i\n",radius2_cell[ic],n_cell); + printf("star_maker8.C line 373: Radius squared is wrong?!? radius =%f, n_cell = %i\n",radius2_cell[ic],n_cell); } rho_wind = mdot_wind*SolidAngle/(pow((*dx),3)); cells_volume += pow((*dx),3); @@ -646,6 +660,7 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float te[ind_cell[ic]] += 0.5*v_wind*v_wind; } + if (m_wind > 0.0) printf(" density in cells changed, dt = %e, t = %e\n",(*dt),(*t)); /* Substract the ejected mass and set dm to be zero */ /*printf("Iso-wind injected: dt = %e s, vwind=%g, n_cell=%"ISYM", xp=(%g, %g, %g),m_star = %e, m_cell=%e, m_wind=%e, rho=%e, umass = %e\n", (*dt)*(*t1),v_wind*(*v1), n_cell, xpold[bb], ypold[bb], zpold[bb],mpold[bb]*umass , m_cell*umass,m_wind*umass, rho_wind*(*d1),umass); @@ -661,6 +676,12 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float } + if (StellarWindFeedback == 4 && bx == NULL) { /*protostellar jets and MS stellar wind*/ + + } + + + /* Loop over grid looking for a cell with mass larger than massthres */ From 20556e1de517441a9462cff86fa32452ada8cd43 Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Wed, 11 Nov 2009 14:07:56 -0500 Subject: [PATCH 09/19] broken?? --HG-- branch : week-of-code --- src/enzo/Grid_StarParticleHandler.C | 8 +++++--- src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C | 8 ++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/enzo/Grid_StarParticleHandler.C b/src/enzo/Grid_StarParticleHandler.C index 07fe0e618..6f9ad7267 100644 --- a/src/enzo/Grid_StarParticleHandler.C +++ b/src/enzo/Grid_StarParticleHandler.C @@ -448,7 +448,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) h2field[index] = BaryonField[H2INum][index] + BaryonField[H2IINum][index]; } } - printf("Star type \n"); + //printf("Star type \n"); /* Set the units. */ float DensityUnits = 1, LengthUnits = 1, TemperatureUnits = 1, @@ -465,7 +465,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) // if (StarParticleCreation > 0 && level == MaximumRefinementLevel) { if (StarParticleCreation > 0) { - + /* Generate a fake grid to keep the particles in. */ grid *tg = new grid; @@ -797,6 +797,8 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) /* This creates sink particles which suck up mass off the grid. */ + if (STARMAKE_METHOD(SINK_PARTICLE)) printf(" Sink Particle\n"); + if (level == MaximumRefinementLevel) printf(" Max Refinement\n"); if (STARMAKE_METHOD(SINK_PARTICLE) && level == MaximumRefinementLevel) { /* Set the density threshold by using the mass in a cell which would have caused another refinement. */ @@ -822,7 +824,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) ny_jet = ParticleAttribute[4]; nz_jet = ParticleAttribute[5]; } - + printf(" test line\n"); if (star_maker8(GridDimension, GridDimension+1, GridDimension+2, &size, BaryonField[DensNum], BaryonField[TENum], BaryonField[GENum], BaryonField[Vel1Num], BaryonField[Vel2Num], BaryonField[Vel3Num], diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C index d5b5b1dda..32e02a405 100644 --- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C @@ -584,7 +584,7 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL if (PutSink == 1 && level == 0) { // set it up on level zero and make it mustrefine // double mass_p = 20.0*1.989e33; - double mass_p = 30.0*1.989e33; + double mass_p = 20.0*1.989e33; mass_p /= MassUnits; double dx = CellWidth[0][0]; double den_p = mass_p / pow(dx,3); @@ -604,9 +604,9 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL ParticleMass[0] = den_p; ParticleNumber[0] = 0; ParticleType[0] = PARTICLE_TYPE_MUST_REFINE; - ParticlePosition[0][0] = 0.5+0.5*dx; - ParticlePosition[1][0] = 0.5+0.5*dx; - ParticlePosition[2][0] = 0.5+0.5*dx; + ParticlePosition[0][0] = 0.5; //+0.5*dx; + ParticlePosition[1][0] = 0.5; //+0.5*dx; + ParticlePosition[2][0] = 0.5; //+0.5*dx; ParticleVelocity[0][0] = 0.0; ParticleVelocity[1][0] = 0.0; From 540357cf83fda0f44f4c03e852fbbbcd8c6eac19 Mon Sep 17 00:00:00 2001 From: Matthew Turk Date: Wed, 11 Nov 2009 14:24:46 -0800 Subject: [PATCH 10/19] This was removed some revisions ago and the removal breaks compilation on any system where ImageMagick is in the path. --HG-- branch : week-of-code --- src/enzo/Makefile | 1 + 1 file changed, 1 insertion(+) diff --git a/src/enzo/Makefile b/src/enzo/Makefile index 3dff4f709..5a4038634 100644 --- a/src/enzo/Makefile +++ b/src/enzo/Makefile @@ -47,6 +47,7 @@ MODULES = VERBOSE = 0 # please remove this +SVN = hg #----------------------------------------------------------------------- # Make.config.settings is used for setting default values to all compile-time From 831b26bff14213b6a20b2cc2482e936e48090986 Mon Sep 17 00:00:00 2001 From: Matthew Turk Date: Wed, 11 Nov 2009 14:30:19 -0800 Subject: [PATCH 11/19] ErrorsExceptions has to be before macros_and_parameters for a 64-bit compilation to succeed... (hooray #defines of built-in types!) --HG-- branch : week-of-code --- src/enzo/hydro_rk/MHDLine.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/enzo/hydro_rk/MHDLine.C b/src/enzo/hydro_rk/MHDLine.C index fef1e47df..ff69cfd39 100644 --- a/src/enzo/hydro_rk/MHDLine.C +++ b/src/enzo/hydro_rk/MHDLine.C @@ -14,10 +14,10 @@ #include #include +#include "ErrorExceptions.h" #include "macros_and_parameters.h" #include "typedefs.h" #include "global_data.h" -#include "ErrorExceptions.h" int HLL_PLM_MHD(float **prim, float **priml, float **primr, float **species, float **colors, float **FluxLine, int ActiveSize, From 8ca1643884f67128485eb2c49dce0042bed68d20 Mon Sep 17 00:00:00 2001 From: John Wise Date: Wed, 11 Nov 2009 18:51:14 -0500 Subject: [PATCH 12/19] Fixed make system to report compiler errors to stdout when breaking. --HG-- branch : week-of-code --- src/enzo/Makefile | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/enzo/Makefile b/src/enzo/Makefile index 3dff4f709..a6d5d395a 100644 --- a/src/enzo/Makefile +++ b/src/enzo/Makefile @@ -157,10 +157,10 @@ verbose: $(EXE).exe $(FC) -c -o $@ $(FFLAGS) $*.f >& $(OUTPUT) ; \ else \ $(FC) -c -o $@ $(FFLAGS) $*.f >> $(OUTPUT) 2>&1 ; \ - fi) - @(if [ -e $@ ]; then \ + fi ; \ + if [ -e $@ ]; then \ rm -f $*.f; \ - else \ + else \ echo; \ echo "$(CPP) $(DEFINES) $(CPPFLAGS) $*.src > $*.f"; \ echo "$(FC) -c -o $@ $(FFLAGS) $*.f"; \ @@ -169,7 +169,7 @@ verbose: $(EXE).exe $(FC) -c -o $@ $(FFLAGS) $*.f; \ echo; \ exit 1; \ - fi) + fi) .src90.f90: @echo "Compiling $<" @@ -180,10 +180,10 @@ verbose: $(EXE).exe $(F90) -c -o $@ $(F90FLAGS) $*.f90 >& $(OUTPUT) ; \ else \ $(F90) -c -o $@ $(F90FLAGS) $*.f90 >> $(OUTPUT) 2>&1 ; \ - fi) - @(if [ -e $@ ]; then \ + fi ; \ + if [ -e $@ ]; then \ rm -f $*.f90; \ - else \ + else \ echo; \ echo "$(CPP) $(DEFINES) $(CPPFLAGS) $*.src90 > $*.f90"; \ echo "$(F90) -c -o $@ $(F90FLAGS) $*.f90"; \ @@ -192,7 +192,7 @@ verbose: $(EXE).exe $(F90) -c -o $@ $(F90FLAGS) $*.f90; \ echo; \ exit 1; \ - fi) + fi) .c.o: @rm -f $@ @@ -203,15 +203,15 @@ verbose: $(EXE).exe else \ $(CC) -c -o $@ $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c \ >> $(OUTPUT) 2>&1 ; \ - fi) - @(if [ ! -e $@ ]; then \ + fi ; \ + if [ ! -e $@ ]; then \ echo; \ echo "$(CC) -c -o $@ $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c"; \ echo; \ $(CC) -c -o $@ $(DEFINES) $(CFLAGS) $(INCLUDES) $*.c;\ echo; \ exit 1; \ - fi) + fi) .C.o: @rm -f $@ @@ -222,15 +222,15 @@ verbose: $(EXE).exe else \ $(CXX) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C \ >> $(OUTPUT) 2>&1 ; \ - fi) - @(if [ ! -e $@ ]; then \ + fi ; \ + if [ ! -e $@ ]; then \ echo; \ echo "$(CXX) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C"; \ echo; \ $(CXX) -c -o $@ $(DEFINES) $(CXXFLAGS) $(INCLUDES) $*.C;\ echo; \ exit 1; \ - fi) + fi) .cu.o: @rm -f $@ @@ -241,15 +241,15 @@ verbose: $(EXE).exe else \ $(CUDACOMPILER) -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) \ $(INCLUDES) $*.cu >> $(OUTPUT) 2>&1 ; \ - fi) - @(if [ ! -e $@ ]; then \ + fi ; \ + if [ ! -e $@ ]; then \ echo; \ echo "$(CUDACOMPILER) -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) $(INCLUDES) $*.cu"; \ echo; \ $(CUDACOMPILER) -c -o $@ $(DEFINES) $(CUDACOMPFLAGS) $(INCLUDES) $*.cu;\ echo; \ exit 1; \ - fi) + fi) #----------------------------------------------------------------------- # Generate all make-generated source files From fd7e81e78499eaac54da67bb7bafb5d6c95f269a Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Wed, 11 Nov 2009 18:53:23 -0500 Subject: [PATCH 13/19] Fixed particle attribute label definitions so all the same (and with 7). This should be improved on later so not redefined so muchgmake install --HG-- branch : week-of-code --- src/enzo/AMRH5writer.C | 14 ++++++++++---- src/enzo/Grid_ConvertToNumpy.C | 4 +++- src/enzo/Grid_Group_ReadGrid.C | 4 +++- src/enzo/Grid_ReadGrid.C | 4 +++- src/enzo/Grid_StarParticleHandler.C | 6 +++--- src/enzo/Grid_WriteCube.C | 7 ++++--- src/enzo/Grid_WriteGrid.C | 4 +++- src/enzo/Grid_WriteGridX.C | 4 +++- src/enzo/New_Grid_ReadGrid.C | 4 +++- src/enzo/New_Grid_WriteGrid.C | 4 +++- 10 files changed, 38 insertions(+), 17 deletions(-) diff --git a/src/enzo/AMRH5writer.C b/src/enzo/AMRH5writer.C index 6bc47787d..3ffa5afcf 100644 --- a/src/enzo/AMRH5writer.C +++ b/src/enzo/AMRH5writer.C @@ -57,8 +57,11 @@ void AMRHDF5Writer::AMRHDF5Create( const char* fileName, {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; const char *ParticleOtherLabel[] = {"particle_type", "particle_number", "particle_mass"}; - const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + /* const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */ + char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + int i; for (i=0; i<3; i++) { @@ -422,8 +425,11 @@ herr_t AMRHDF5Writer::writeParticles ( const int nPart, {"particle_position_x", "particle_position_y", "particle_position_z"}; const char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; - const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + + /* const char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction", "p5", "p6"}; */ sprintf(gridDataName, "/grid-%d", gridId); if (nBaryonFields > 0) diff --git a/src/enzo/Grid_ConvertToNumpy.C b/src/enzo/Grid_ConvertToNumpy.C index f1a8f4263..4cdfca5b9 100644 --- a/src/enzo/Grid_ConvertToNumpy.C +++ b/src/enzo/Grid_ConvertToNumpy.C @@ -32,7 +32,9 @@ void grid::ConvertToNumpy(int GridID, PyArrayObject *container[], int ParentID, char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction", "p5", "p6"};*/ this->DebugCheck("Converting to NumPy arrays"); diff --git a/src/enzo/Grid_Group_ReadGrid.C b/src/enzo/Grid_Group_ReadGrid.C index 181bb068c..dba8c91da 100644 --- a/src/enzo/Grid_Group_ReadGrid.C +++ b/src/enzo/Grid_Group_ReadGrid.C @@ -85,7 +85,9 @@ int grid::Group_ReadGrid(FILE *fptr, int GridID, HDF5_hid_t file_id, char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /*char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction", "p5", "p6"};*/ #ifdef IO_LOG int io_log = 1; diff --git a/src/enzo/Grid_ReadGrid.C b/src/enzo/Grid_ReadGrid.C index 7a9cd9f37..9c12fe095 100644 --- a/src/enzo/Grid_ReadGrid.C +++ b/src/enzo/Grid_ReadGrid.C @@ -78,7 +78,9 @@ int grid::ReadGrid(FILE *fptr, int GridID, char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction"};*/ #ifdef USE_HDF4 Eint32 TempIntArray2[MAX_DIMENSION]; diff --git a/src/enzo/Grid_StarParticleHandler.C b/src/enzo/Grid_StarParticleHandler.C index 6f9ad7267..e07c0fd8d 100644 --- a/src/enzo/Grid_StarParticleHandler.C +++ b/src/enzo/Grid_StarParticleHandler.C @@ -797,8 +797,8 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) /* This creates sink particles which suck up mass off the grid. */ - if (STARMAKE_METHOD(SINK_PARTICLE)) printf(" Sink Particle\n"); - if (level == MaximumRefinementLevel) printf(" Max Refinement\n"); + //if (STARMAKE_METHOD(SINK_PARTICLE)) printf(" Sink Particle\n"); + //if (level == MaximumRefinementLevel) printf(" Max Refinement\n"); if (STARMAKE_METHOD(SINK_PARTICLE) && level == MaximumRefinementLevel) { /* Set the density threshold by using the mass in a cell which would have caused another refinement. */ @@ -824,7 +824,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) ny_jet = ParticleAttribute[4]; nz_jet = ParticleAttribute[5]; } - printf(" test line\n"); + //printf(" test line\n"); if (star_maker8(GridDimension, GridDimension+1, GridDimension+2, &size, BaryonField[DensNum], BaryonField[TENum], BaryonField[GENum], BaryonField[Vel1Num], BaryonField[Vel2Num], BaryonField[Vel3Num], diff --git a/src/enzo/Grid_WriteCube.C b/src/enzo/Grid_WriteCube.C index 06a953852..935b016c1 100644 --- a/src/enzo/Grid_WriteCube.C +++ b/src/enzo/Grid_WriteCube.C @@ -118,9 +118,10 @@ int grid::WriteCube(char *base_name, int grid_id, int TGdims[]) char *ParticleMassLabel = "particle_mass"; char *ParticleTypeLabel = "particle_type"; char *ParticleIndexLabel = "particle_index"; - - char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction"};*/ #ifdef IO_LOG int io_log = 1; #else diff --git a/src/enzo/Grid_WriteGrid.C b/src/enzo/Grid_WriteGrid.C index eca030624..2c9b0ad66 100644 --- a/src/enzo/Grid_WriteGrid.C +++ b/src/enzo/Grid_WriteGrid.C @@ -83,7 +83,9 @@ int grid::WriteGrid(FILE *fptr, char *base_name, int grid_id) char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction"};*/ char *SmoothedDMLabel[] = {"Dark_Matter_Density", "Velocity_Dispersion", "Particle_x-velocity", "Particle_y-velocity", "Particle_z-velocity"}; diff --git a/src/enzo/Grid_WriteGridX.C b/src/enzo/Grid_WriteGridX.C index fd09b63d2..44147e668 100644 --- a/src/enzo/Grid_WriteGridX.C +++ b/src/enzo/Grid_WriteGridX.C @@ -73,7 +73,9 @@ int grid::WriteGridX(FILE *fptr, char *base_name, int grid_id) char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction"};*/ #ifdef IO_LOG int io_log = 1; #else diff --git a/src/enzo/New_Grid_ReadGrid.C b/src/enzo/New_Grid_ReadGrid.C index 860146b7b..1cd562162 100644 --- a/src/enzo/New_Grid_ReadGrid.C +++ b/src/enzo/New_Grid_ReadGrid.C @@ -86,7 +86,9 @@ int grid::Group_ReadGrid(FILE *fptr, int GridID, HDF5_hid_t file_id, char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction"};*/ if(ReadText){ diff --git a/src/enzo/New_Grid_WriteGrid.C b/src/enzo/New_Grid_WriteGrid.C index 1fd222955..424a43b2f 100644 --- a/src/enzo/New_Grid_WriteGrid.C +++ b/src/enzo/New_Grid_WriteGrid.C @@ -88,7 +88,9 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f char *ParticleVelocityLabel[] = {"particle_velocity_x", "particle_velocity_y", "particle_velocity_z"}; char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", - "metallicity_fraction", "alpha_fraction"}; + "metallicity_fraction", "particle_jet_x", "particle_jet_y", "particle_jet_z", "alpha_fraction"}; + /* char *ParticleAttributeLabel[] = {"creation_time", "dynamical_time", + "metallicity_fraction", "alpha_fraction"};*/ char *SmoothedDMLabel[] = {"Dark_Matter_Density", "Velocity_Dispersion", "Particle_x-velocity", "Particle_y-velocity", "Particle_z-velocity"}; From 02570ee6cadbef8f2e3a7c8ce1e33cbab5cba5dc Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Thu, 12 Nov 2009 11:45:09 -0500 Subject: [PATCH 14/19] SciNet makefile tests --HG-- branch : week-of-code --- src/enzo/Make.mach.scinet | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/src/enzo/Make.mach.scinet b/src/enzo/Make.mach.scinet index d93cc33f6..09b86f110 100644 --- a/src/enzo/Make.mach.scinet +++ b/src/enzo/Make.mach.scinet @@ -32,7 +32,7 @@ LOCAL_PYTHON_INSTALL = /scinet/gpc/tools/Python/Python262/ # Compiler settings #----------------------------------------------------------------------- -MACH_CPP = mpiCC +MACH_CPP = cpp # With MPI @@ -70,11 +70,11 @@ MACH_DEFINES = -DLINUX -DH5_USE_16_API -DTRANSFER -DWINDS -DIO_LOG #----------------------------------------------------------------------- MACH_CPPFLAGS = -P -traditional -MACH_CFLAGS = -c -MACH_CXXFLAGS = -c -MACH_FFLAGS = -c -132 -MACH_F90FLAGS = -MACH_LDFLAGS = -Wl +MACH_CFLAGS = +MACH_CXXFLAGS = +MACH_FFLAGS = -fno-second-underscore +MACH_F90FLAGS = -fno-second-underscore +MACH_LDFLAGS = #,-static #----------------------------------------------------------------------- @@ -101,8 +101,8 @@ MACH_FFLAGS_REAL_64 = -r8 # # *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING *** -MACH_OPT_WARN = -MACH_OPT_DEBUG = -g +MACH_OPT_WARN = -Wall # Flags for verbose compiler warnings +MACH_OPT_DEBUG = -O0 -g MACH_OPT_HIGH = -O2 MACH_OPT_AGGRESSIVE = -O3 @@ -129,14 +129,11 @@ MACH_INCLUDES_HYPRE = $(LOCAL_INCLUDES_HYPRE) # variables will be properly set # -LOCAL_LIBS_MPI = -L$(LOCAL_MPI_INSTALL)/lib -LOCAL_LIBS_HDF5 = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz +LOCAL_LIBS_MPI = $(LOCAL_MPI_INSTALL)/lib +LOCAL_LIBS_HDF5 = $(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz LOCAL_LIBS_PYTHON = $(LOCAL_PYTHON_INSTALL)/lib/python2.6/config/libpython2.6.a -LOCAL_LIBS_MACH = -L$(INT_DIR)/fc/10.1.018/@sys/lib/ \ - -limf -lifcore -lifport \ - -L$(INT_DIR)/cc/10.1.018/@sys/lib/ \ - -lstdc++ -lg2c +LOCAL_LIBS_MACH = -limf -lifcore -lifport -lstdc++ -lg2c MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) MACH_LIBS_MPI = $(LOCAL_LIBS_MPI) From 2c88991e4f0fc55ef5a64eb7c813565fa0763c7a Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Thu, 12 Nov 2009 12:11:52 -0500 Subject: [PATCH 15/19] Linked libraries found for SciNet --HG-- branch : week-of-code --- src/enzo/Make.mach.scinet | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/enzo/Make.mach.scinet b/src/enzo/Make.mach.scinet index 09b86f110..6a0f5b57a 100644 --- a/src/enzo/Make.mach.scinet +++ b/src/enzo/Make.mach.scinet @@ -22,6 +22,7 @@ MACH_FILE = Make.mach.scinet #----------------------------------------------------------------------- # Install paths (local variables) #----------------------------------------------------------------------- +INT_DIR=/scinet/gpc/intel/Compiler/11.1/056 LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-intel-v11.0-ofed #LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-gcc-v4.4.0-ofed @@ -129,11 +130,16 @@ MACH_INCLUDES_HYPRE = $(LOCAL_INCLUDES_HYPRE) # variables will be properly set # -LOCAL_LIBS_MPI = $(LOCAL_MPI_INSTALL)/lib +LOCAL_LIBS_MPI = $(LOCAL_MPI_INSTALL)/lib -lmyriexpress LOCAL_LIBS_HDF5 = $(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz LOCAL_LIBS_PYTHON = $(LOCAL_PYTHON_INSTALL)/lib/python2.6/config/libpython2.6.a -LOCAL_LIBS_MACH = -limf -lifcore -lifport -lstdc++ -lg2c +#LOCAL_LIBS_MACH = -limf -lifcore -lifport -lstdc++ -lg2c +LOCAL_LIBS_MACH = -L$(INT_DIR)/lib/intel64 \ + -limf -lifcore -lifport \ + -lstdc++ -lg2c + + MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) MACH_LIBS_MPI = $(LOCAL_LIBS_MPI) From 2a2fc268c31412da9d8e1d8810a09c79ea21e763 Mon Sep 17 00:00:00 2001 From: Elizabeth Harper-Clark Date: Thu, 12 Nov 2009 12:25:32 -0500 Subject: [PATCH 16/19] Alternative hdf libraries (no mpi) added to Make.mach.scinet --HG-- branch : week-of-code --- src/enzo/Make.mach.scinet | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/enzo/Make.mach.scinet b/src/enzo/Make.mach.scinet index 6a0f5b57a..9de6c2b15 100644 --- a/src/enzo/Make.mach.scinet +++ b/src/enzo/Make.mach.scinet @@ -27,6 +27,7 @@ INT_DIR=/scinet/gpc/intel/Compiler/11.1/056 LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-intel-v11.0-ofed #LOCAL_MPI_INSTALL = /scinet/gpc/mpi/openmpi/1.3.2-gcc-v4.4.0-ofed LOCAL_HDF5_INSTALL = /scinet/gpc/Libraries/HDF5-1.8.3-v18-openmpi +#LOCAL_HDF5_INSTALL = /scinet/gpc/Libraries/HDF5-1.8.3-v18-mvapich2 LOCAL_PYTHON_INSTALL = /scinet/gpc/tools/Python/Python262/ #----------------------------------------------------------------------- @@ -37,16 +38,16 @@ MACH_CPP = cpp # With MPI -MACH_CC_MPI = mpicc -MACH_CXX_MPI = mpiCC -MACH_FC_MPI = mpif77 -MACH_F90_MPI = mpif90 -MACH_LD_MPI = mpiCC -#MACH_CC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpicc -#MACH_CXX_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC -#MACH_FC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif77 -#MACH_F90_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif90 -#MACH_LD_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC +#MACH_CC_MPI = mpicc +#MACH_CXX_MPI = mpiCC +#MACH_FC_MPI = mpif77 +#MACH_F90_MPI = mpif90 +#MACH_LD_MPI = mpiCC +MACH_CC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpicc +MACH_CXX_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC +MACH_FC_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif77 +MACH_F90_MPI = $(LOCAL_MPI_INSTALL)/bin/mpif90 +MACH_LD_MPI = $(LOCAL_MPI_INSTALL)/bin/mpiCC # Without MPI @@ -73,8 +74,10 @@ MACH_DEFINES = -DLINUX -DH5_USE_16_API -DTRANSFER -DWINDS -DIO_LOG MACH_CPPFLAGS = -P -traditional MACH_CFLAGS = MACH_CXXFLAGS = -MACH_FFLAGS = -fno-second-underscore -MACH_F90FLAGS = -fno-second-underscore +MACH_FFLAGS = +#-fno-second-underscore +MACH_F90FLAGS = +#-fno-second-underscore MACH_LDFLAGS = #,-static From 61f579108afd048885c9f63a06446f9f6533525d Mon Sep 17 00:00:00 2001 From: John Wise Date: Fri, 13 Nov 2009 08:51:22 -0500 Subject: [PATCH 17/19] Consolidated all photo-heating fields into one. Added option for hydrogen only ionization (RadiativeTransferHydrogenOnly) so we don't have to allocate the kphHeI and kphHeII fields in those sims. --HG-- branch : week-of-code --- src/enzo/ExternalBoundary.h | 2 + src/enzo/ExternalBoundary_AddField.C | 6 +- src/enzo/Grid.h | 7 +- src/enzo/Grid_AddH2Dissociation.C | 10 +-- src/enzo/Grid_AllocateInterpolatedRadiation.C | 35 +++++----- src/enzo/Grid_ComputeCoolingTime.C | 15 ++-- src/enzo/Grid_ComputeLuminosity.C | 16 ++--- src/enzo/Grid_ComputePhotonTimestepHII.C | 7 +- .../Grid_ConvertToCellCenteredRadiation.C | 35 ++++------ .../Grid_CorrectRadiationIncompleteness.C | 13 ++-- src/enzo/Grid_DetectIonizationFrontApprox.C | 12 ++-- src/enzo/Grid_ElectronFractionEstimate.C | 29 +++----- src/enzo/Grid_FinalizeRadiationFields.C | 33 ++++----- .../Grid_FlagCellsToBeRefinedByOpticalDepth.C | 11 +-- src/enzo/Grid_Group_WriteGrid.C | 2 +- .../Grid_IdentifyRadiativeTransferFields.C | 30 ++++---- .../Grid_InitializeRadiativeTransferFields.C | 40 ++++++----- src/enzo/Grid_PhotonTestInitializeGrid.C | 14 ++-- src/enzo/Grid_Shine.C | 1 + src/enzo/Grid_SolveCoupledRateEquations.C | 17 ++--- ...Grid_SolveHighDensityPrimordialChemistry.C | 9 +-- src/enzo/Grid_SolveRadiativeCooling.C | 15 ++-- src/enzo/Grid_SolveRateAndCoolEquations.C | 16 ++--- src/enzo/Grid_SolveRateEquations.C | 17 ++--- src/enzo/Grid_StarParticleHandler.C | 2 +- src/enzo/Grid_TransportPhotonPackages.C | 24 ++----- src/enzo/Grid_WalkPhotonPackage.C | 15 ++-- src/enzo/Make.config.objects | 2 + src/enzo/PhotonGrid_Methods.h | 30 ++++---- src/enzo/PhotonTestInitialize.C | 8 +-- src/enzo/RadiativeTransferInitialize.C | 50 ++++++++----- src/enzo/RadiativeTransferParameters.h | 2 + src/enzo/RadiativeTransferReadParameters.C | 3 + src/enzo/RadiativeTransferWriteParameters.C | 4 +- src/enzo/StarParticleAddFeedback.C | 2 +- src/enzo/cool1d_multi.src | 13 ++-- src/enzo/cool1d_sep.src | 4 +- src/enzo/cool_multi_lum.src | 6 +- src/enzo/cool_multi_time.src | 6 +- .../hydro_rk/Grid_Collapse3DInitializeGrid.C | 8 +-- .../hydro_rk/Grid_TurbulenceInitializeGrid.C | 8 +-- src/enzo/multi_cool.src | 9 +-- src/enzo/solve_rate.src | 43 ++++++++---- src/enzo/solve_rate_cool.src | 70 ++++++++++++------- src/enzo/typedefs.h | 2 +- 45 files changed, 334 insertions(+), 369 deletions(-) diff --git a/src/enzo/ExternalBoundary.h b/src/enzo/ExternalBoundary.h index eef41c7c5..c574ed655 100644 --- a/src/enzo/ExternalBoundary.h +++ b/src/enzo/ExternalBoundary.h @@ -126,6 +126,8 @@ class ExternalBoundary int DetachForcingFromBaryonFields(); int AddField(int FieldType); + int DeleteObsoleteFields(int *ObsoleteFields, + int NumberOfObsoleteFields); }; diff --git a/src/enzo/ExternalBoundary_AddField.C b/src/enzo/ExternalBoundary_AddField.C index d3b9e446c..7521759ba 100644 --- a/src/enzo/ExternalBoundary_AddField.C +++ b/src/enzo/ExternalBoundary_AddField.C @@ -1,9 +1,9 @@ /*********************************************************************** / -/ EXTERNAL BOUNDARY CLASS (READ THE EXTERNAL BOUNDARY VALUES) +/ EXTERNAL BOUNDARY CLASS (ADD EXTERNAL BOUNDARY VALUES FOR A NEW FIELD) / -/ written by: Greg Bryan -/ date: November, 1994 +/ written by: John Wise +/ date: March, 2009 / modified1: / / PURPOSE: diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 0c6e5149a..e0b1437e7 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -1732,6 +1732,8 @@ int CollapseTestInitializeGrid(int NumberOfSpheres, int ReadRandomForcingFields(FILE *main_file_pointer); int AddFields(int TypesToAdd[], int NumberOfFields); + int DeleteObsoleteFields(int *ObsoleteFields, + int NumberOfObsoleteFields); inline bool isLocal () {return MyProcessorNumber == ProcessorNumber; }; @@ -1968,9 +1970,8 @@ int CollapseTestInitializeGrid(int NumberOfSpheres, // Radiative transfer methods that don't fit in the TRANSFER define //------------------------------------------------------------------------ - int IdentifyRadiativeTransferFields(int &kphHINum, int &gammaHINum, - int &kphHeINum, int &gammaHeINum, - int &kphHeIINum, int &gammaHeIINum, + int IdentifyRadiativeTransferFields(int &kphHINum, int &gammaNum, + int &kphHeINum, int &kphHeIINum, int &kdissH2INum); #ifdef TRANSFER diff --git a/src/enzo/Grid_AddH2Dissociation.C b/src/enzo/Grid_AddH2Dissociation.C index 7fc07cb2d..a41632c8b 100644 --- a/src/enzo/Grid_AddH2Dissociation.C +++ b/src/enzo/Grid_AddH2Dissociation.C @@ -60,13 +60,9 @@ int grid::AddH2Dissociation(Star *AllStars) /* Get photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, kphHeIINum, + kdissH2INum); /* For now, initialize H2 photo-dissociation field. */ diff --git a/src/enzo/Grid_AllocateInterpolatedRadiation.C b/src/enzo/Grid_AllocateInterpolatedRadiation.C index aac873cf9..48aea1e77 100644 --- a/src/enzo/Grid_AllocateInterpolatedRadiation.C +++ b/src/enzo/Grid_AllocateInterpolatedRadiation.C @@ -30,14 +30,9 @@ int grid::AllocateInterpolatedRadiation() if (MyProcessorNumber != ProcessorNumber) return SUCCESS; - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); int i,j,k, index, field, dim, vcsize = 1; @@ -50,28 +45,32 @@ int grid::AllocateInterpolatedRadiation() for (field = 0; field < NumberOfBaryonFields; field++) if (FieldsToInterpolate[field] == TRUE) { + rkph = FieldUndefined; + rgamma = FieldUndefined; switch (FieldType[field]) { case HIDensity: rkph = kphHINum; - rgamma = gammaHINum; + rgamma = gammaNum; break; case HeIDensity: - rkph = kphHeINum; - rgamma = gammaHeINum; + if (RadiativeTransferHydrogenOnly == FALSE) + rkph = kphHeINum; break; case HeIIDensity: - rkph = kphHeIINum; - rgamma = gammaHeIINum; + if (RadiativeTransferHydrogenOnly == FALSE) + rkph = kphHeIINum; break; } // ENDSWITCH field - if (InterpolatedField[rkph] == NULL) + if (InterpolatedField[rkph] == NULL && rkph != FieldUndefined) InterpolatedField[rkph] = new float[vcsize]; - if (InterpolatedField[rgamma] == NULL) + if (InterpolatedField[rgamma] == NULL && rgamma != FieldUndefined) InterpolatedField[rgamma] = new float[vcsize]; - for (i = 0; i < vcsize; i++) { - InterpolatedField[rkph][i] = + for (i = 0; i < vcsize; i++) + InterpolatedField[rkph][i] = 0.0f; + if (rgamma != FieldUndefined) + for (i = 0; i < vcsize; i++) InterpolatedField[rgamma][i] = 0.0f; - } + } // ENDIF fields to interpolate return SUCCESS; diff --git a/src/enzo/Grid_ComputeCoolingTime.C b/src/enzo/Grid_ComputeCoolingTime.C index 286337767..8d963c2b0 100644 --- a/src/enzo/Grid_ComputeCoolingTime.C +++ b/src/enzo/Grid_ComputeCoolingTime.C @@ -80,7 +80,7 @@ extern "C" void FORTRAN_NAME(cool_multi_time)( float *metala, int *n_xe, float *xe_start, float *xe_end, float *inutot, int *iradfield, int *nfreq, int *imetalregen, int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p, - int *iradtrans, float *gammaHI, float *gammaHeI, float *gammaHeII, + int *iradtrans, float *photogamma, int *icmbTfloor, int *iClHeat, int *iClMMW, float *clMetNorm, float *clEleFra, int *clGridRank, int *clGridDim, float *clPar1, float *clPar2, float *clPar3, float *clPar4, float *clPar5, @@ -141,13 +141,9 @@ int grid::ComputeCoolingTime(float *cooling_time) /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Get easy to handle pointers for each variable. */ @@ -254,8 +250,7 @@ int grid::ComputeCoolingTime(float *cooling_time) &RadiationData.NumberOfFrequencyBins, &RadiationFieldRecomputeMetalRates, &RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor, - &RadiativeTransfer, BaryonField[gammaHINum], - BaryonField[gammaHeINum], BaryonField[gammaHeIINum], + &RadiativeTransfer, BaryonField[gammaNum], &CloudyCoolingData.CMBTemperatureFloor, &CloudyCoolingData.IncludeCloudyHeating, &CloudyCoolingData.IncludeCloudyMMW, &CloudyCoolingData.CloudyMetallicityNormalization, diff --git a/src/enzo/Grid_ComputeLuminosity.C b/src/enzo/Grid_ComputeLuminosity.C index 9f4ae68f6..11dae5bc6 100644 --- a/src/enzo/Grid_ComputeLuminosity.C +++ b/src/enzo/Grid_ComputeLuminosity.C @@ -67,7 +67,7 @@ extern "C" void FORTRAN_NAME(cool_multi_lum)( float *metala, int *n_xe, float *xe_start, float *xe_end, float *inutot, int *iradfield, int *nfreq, int *imetalregen, int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p, - int *iradtrans, float *gammaHI, float *gammaHeI, float *gammaHeII); + int *iradtrans, float *photogamma); extern "C" void FORTRAN_NAME(cool_time)( float *d, float *e, float *ge, float *u, float *v, float *w, float *cooltime, @@ -115,14 +115,9 @@ int grid::ComputeLuminosity(float *luminosity, int NumberOfLuminosityFields) /* Find photo-ionization fields */ - int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, kphHeINum, kphHeIINum, kdissH2INum, gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Get easy to handle pointers for each variable. */ @@ -239,8 +234,7 @@ int grid::ComputeLuminosity(float *luminosity, int NumberOfLuminosityFields) &RadiationData.NumberOfFrequencyBins, &RadiationFieldRecomputeMetalRates, &RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor, - &RadiativeTransfer, BaryonField[gammaHINum], - BaryonField[gammaHeINum], BaryonField[gammaHeIINum]); + &RadiativeTransfer, BaryonField[gammaNum]); else { #ifdef UNUSED FORTRAN_NAME(cool_time)( diff --git a/src/enzo/Grid_ComputePhotonTimestepHII.C b/src/enzo/Grid_ComputePhotonTimestepHII.C index b24c11b24..8e0ca5f05 100644 --- a/src/enzo/Grid_ComputePhotonTimestepHII.C +++ b/src/enzo/Grid_ComputePhotonTimestepHII.C @@ -76,10 +76,9 @@ float grid::ComputePhotonTimestepHII(float DensityUnits, float LengthUnits, /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum); + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Compute temperature field */ diff --git a/src/enzo/Grid_ConvertToCellCenteredRadiation.C b/src/enzo/Grid_ConvertToCellCenteredRadiation.C index 276041ba6..def4e4b9f 100644 --- a/src/enzo/Grid_ConvertToCellCenteredRadiation.C +++ b/src/enzo/Grid_ConvertToCellCenteredRadiation.C @@ -38,31 +38,28 @@ int grid::ConvertToCellCenteredRadiation() /* Find radiative transfer fields. */ - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); int rkph, rgamma, field; for (field = 0; field < NumberOfBaryonFields; field++) if (FieldsToInterpolate[field] == TRUE) { + rkph = FieldUndefined; + rgamma = FieldUndefined; switch (FieldType[field]) { case HIDensity: rkph = kphHINum; - rgamma = gammaHINum; + rgamma = gammaNum; break; case HeIDensity: - rkph = kphHeINum; - rgamma = gammaHeINum; + if (RadiativeTransferHydrogenOnly == FALSE) + rkph = kphHeINum; break; case HeIIDensity: - rkph = kphHeIINum; - rgamma = gammaHeIINum; + if (RadiativeTransferHydrogenOnly == FALSE) + rkph = kphHeIINum; break; } // ENDSWITCH field @@ -75,15 +72,11 @@ int grid::ConvertToCellCenteredRadiation() ENZO_FAIL(""); } - if (this->ComputeCellCenteredField(rkph) == FAIL) { - fprintf(stderr, "Error in grid->ComputeCellCenteredField.\n"); - ENZO_FAIL(""); - } + if (rkph != FieldUndefined) + this->ComputeCellCenteredField(rkph); - if (this->ComputeCellCenteredField(rgamma) == FAIL) { - fprintf(stderr, "Error in grid->ComputeCellCenteredField.\n"); - ENZO_FAIL(""); - } + if (rgamma != FieldUndefined) + this->ComputeCellCenteredField(rgamma); } // ENDIF fields to interpolate diff --git a/src/enzo/Grid_CorrectRadiationIncompleteness.C b/src/enzo/Grid_CorrectRadiationIncompleteness.C index 89a6c4b68..44beeadbb 100644 --- a/src/enzo/Grid_CorrectRadiationIncompleteness.C +++ b/src/enzo/Grid_CorrectRadiationIncompleteness.C @@ -40,14 +40,9 @@ int grid::CorrectRadiationIncompleteness(void) /* Find radiative transfer fields. */ - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); int nsrc; @@ -55,7 +50,7 @@ int grid::CorrectRadiationIncompleteness(void) if (BaryonField[kphHeIINum][i] > 0) { nsrc = max(int(BaryonField[kphHeIINum][i]-0.1), 1); BaryonField[kphHINum][i] /= BaryonField[kphHeIINum][i] / nsrc; - BaryonField[gammaHINum][i] /= BaryonField[kphHeIINum][i] / nsrc; + BaryonField[gammaNum][i] /= BaryonField[kphHeIINum][i] / nsrc; } #endif /* TRANSFER */ diff --git a/src/enzo/Grid_DetectIonizationFrontApprox.C b/src/enzo/Grid_DetectIonizationFrontApprox.C index 4bf21983f..27508bff9 100644 --- a/src/enzo/Grid_DetectIonizationFrontApprox.C +++ b/src/enzo/Grid_DetectIonizationFrontApprox.C @@ -37,20 +37,16 @@ int grid::DetectIonizationFrontApprox(float TemperatureUnits) float maxDT, dT; float TOLERANCE = 5e3; // if abs[dT(i+1) - dT(i)] > TOLERANCE, then we return TRUE - int eNum, kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, - gammaHeIINum, kdissH2INum; + int eNum, kphHINum, gammaNum, kphHeINum, kphHeIINum, + kdissH2INum; if (DualEnergyFormalism) eNum = FindField(InternalEnergy, FieldType, NumberOfBaryonFields); else eNum = FindField(TotalEnergy, FieldType, NumberOfBaryonFields); - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); maxDT = -1e20; TOLERANCE /= TemperatureUnits; diff --git a/src/enzo/Grid_ElectronFractionEstimate.C b/src/enzo/Grid_ElectronFractionEstimate.C index 0c8239b36..6ebe7bf57 100644 --- a/src/enzo/Grid_ElectronFractionEstimate.C +++ b/src/enzo/Grid_ElectronFractionEstimate.C @@ -47,27 +47,16 @@ int grid::ElectronFractionEstimate(float dt) int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num; int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, DINum, DIINum, HDINum; - int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; + int kphHINum, kphHeINum, kphHeIINum, kdissH2INum, gammaNum; - if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, - Vel3Num, TENum) == FAIL) { - fprintf(stderr, "Error in IdentifyPhysicalQuantities.\n"); - ENZO_FAIL(""); - } + IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, + Vel3Num, TENum); - if (IdentifySpeciesFields(DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, - HMNum, H2INum, H2IINum, DINum, DIINum, HDINum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifySpeciesFields.\n"); - ENZO_FAIL(""); - } + IdentifySpeciesFields(DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, + HMNum, H2INum, H2IINum, DINum, DIINum, HDINum); - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* If using cosmology, get units. */ @@ -223,7 +212,7 @@ int grid::ElectronFractionEstimate(float dt) // printf("edot[0] = %"GSYM"\n", edot); - edotplus = CoolData.ipiht * BaryonField[gammaHINum][index] * rtunits + edotplus = CoolData.ipiht * BaryonField[gammaNum][index] * rtunits * proper_hi / dom; edotplus = min(edotplus, max_edotplus); @@ -238,7 +227,7 @@ int grid::ElectronFractionEstimate(float dt) // printf("(%"ISYM" %"ISYM" %"ISYM") tem=%.2g, gg=%.2g, edot=%.2g, " // "ge=%.2g\n", -// i, j, k, temperature, BaryonField[gammaHINum][index], +// i, j, k, temperature, BaryonField[gammaNum][index], // edot, BaryonField[GENum][index]); /* Estimate ionization fraction */ diff --git a/src/enzo/Grid_FinalizeRadiationFields.C b/src/enzo/Grid_FinalizeRadiationFields.C index 856392f09..9bd199512 100644 --- a/src/enzo/Grid_FinalizeRadiationFields.C +++ b/src/enzo/Grid_FinalizeRadiationFields.C @@ -55,14 +55,9 @@ int grid::FinalizeRadiationFields(void) /* Find radiative transfer fields. */ - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Get units. */ @@ -76,21 +71,23 @@ int grid::FinalizeRadiationFields(void) float DensityConversion = DensityUnits / 1.673e-24; float factor = DensityConversion * CellVolume; - float invrho; for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { index = GRIDINDEX_NOGHOST(GridStartIndex[0],j,k); - for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) { - invrho = 1.0 / (factor * BaryonField[HINum][index]); - BaryonField[kphHINum][index] *= invrho; - BaryonField[gammaHINum][index] *= invrho; - BaryonField[kphHeINum][index] /= (factor * BaryonField[HeINum][index]); - BaryonField[gammaHeINum][index] /= (factor * BaryonField[HeINum][index]); - //BaryonField[kphHeIINum][index] /= (factor * BaryonField[HeIINum][index]); - BaryonField[gammaHeIINum][index] /= (factor * BaryonField[HeIINum][index]); - } // ENDFOR i + for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) + BaryonField[kphHINum][index] /= (factor * BaryonField[HINum][index]); } // ENDFOR j + + if (RadiativeTransferHydrogenOnly == FALSE) + for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) + for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { + index = GRIDINDEX_NOGHOST(GridStartIndex[0],j,k); + for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) { + BaryonField[kphHeINum][index] /= (factor * BaryonField[HeINum][index]); + BaryonField[kphHeIINum][index] /= (factor * BaryonField[HeIINum][index]); + } // ENDFOR i + } // ENDFOR j if (RadiativeTransferHIIRestrictedTimestep && this->IndexOfMaximumkph >= 0) diff --git a/src/enzo/Grid_FlagCellsToBeRefinedByOpticalDepth.C b/src/enzo/Grid_FlagCellsToBeRefinedByOpticalDepth.C index 062237c01..b23fb17c2 100644 --- a/src/enzo/Grid_FlagCellsToBeRefinedByOpticalDepth.C +++ b/src/enzo/Grid_FlagCellsToBeRefinedByOpticalDepth.C @@ -64,14 +64,9 @@ int grid::FlagCellsToBeRefinedByOpticalDepth() /* Find radiative transfer fields. */ - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Get density units. */ diff --git a/src/enzo/Grid_Group_WriteGrid.C b/src/enzo/Grid_Group_WriteGrid.C index 8cd888487..bb3501b4f 100644 --- a/src/enzo/Grid_Group_WriteGrid.C +++ b/src/enzo/Grid_Group_WriteGrid.C @@ -333,7 +333,7 @@ int grid::Group_WriteGrid(FILE *fptr, char *base_name, int grid_id, HDF5_hid_t f DataUnits[field] = "none"; } - printf("OutPut Field2: %d \n", field); + //printf("OutPut Field2: %d \n", field); WriteStringAttr(dset_id, "Label", DataLabel[field], log_fptr); WriteStringAttr(dset_id, "Units", DataUnits[field], log_fptr); WriteStringAttr(dset_id, "Format", "e10.4", log_fptr); diff --git a/src/enzo/Grid_IdentifyRadiativeTransferFields.C b/src/enzo/Grid_IdentifyRadiativeTransferFields.C index a376146ba..bb625ac73 100644 --- a/src/enzo/Grid_IdentifyRadiativeTransferFields.C +++ b/src/enzo/Grid_IdentifyRadiativeTransferFields.C @@ -27,33 +27,31 @@ int FindField(int f, int farray[], int n); -int grid::IdentifyRadiativeTransferFields(int &kphHINum, int &gammaHINum, - int &kphHeINum, int &gammaHeINum, - int &kphHeIINum, int &gammaHeIINum, +int grid::IdentifyRadiativeTransferFields(int &kphHINum, int &gammaNum, + int &kphHeINum, + int &kphHeIINum, int &kdissH2INum) { - kphHINum = gammaHINum = kphHeINum = gammaHeINum = kphHeIINum = gammaHeIINum = kdissH2INum = 0; + kphHINum = gammaNum = kphHeINum = kphHeIINum = kdissH2INum = 0; if (RadiativeTransfer) { kphHINum = FindField(kphHI, FieldType, NumberOfBaryonFields); - gammaHINum = FindField(gammaHI, FieldType, NumberOfBaryonFields); - kphHeINum = FindField(kphHeI, FieldType, NumberOfBaryonFields); - gammaHeINum = FindField(gammaHeI, FieldType, NumberOfBaryonFields); - kphHeIINum = FindField(kphHeII, FieldType, NumberOfBaryonFields); - gammaHeIINum= FindField(gammaHeII, FieldType, NumberOfBaryonFields); + gammaNum = FindField(PhotoGamma, FieldType, NumberOfBaryonFields); + if (RadiativeTransferHydrogenOnly == FALSE) { + kphHeINum = FindField(kphHeI, FieldType, NumberOfBaryonFields); + kphHeIINum = FindField(kphHeII, FieldType, NumberOfBaryonFields); + } if (MultiSpecies > 1) kdissH2INum = FindField(kdissH2I, FieldType, NumberOfBaryonFields); - - if (kphHINum<0 || gammaHINum<0 || kphHeINum<0 || gammaHeINum<0 || - kphHeIINum<0 || gammaHeIINum<0 || kdissH2INum<0) { + if (kphHINum<0 || gammaNum<0 || kphHeINum<0 || kphHeIINum<0 || + kdissH2INum<0) { fprintf(stderr, "Grid_IdentifyRadiativeTransferFields: failed\n"); - fprintf(stderr, "kphHI=%"ISYM", gammaHI=%"ISYM", kphHeI=%"ISYM", gammaHeI=%"ISYM", " - "kphHeII=%"ISYM", gammaHeII=%"ISYM", kdissH2I=%"ISYM"\n", - kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, - gammaHeIINum, kdissH2INum); + fprintf(stderr, "kphHI=%"ISYM", gamma=%"ISYM", kphHeI=%"ISYM", " + "kphHeII=%"ISYM", kdissH2I=%"ISYM"\n", + kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum); ENZO_FAIL(""); } } diff --git a/src/enzo/Grid_InitializeRadiativeTransferFields.C b/src/enzo/Grid_InitializeRadiativeTransferFields.C index 365582815..07895bb55 100644 --- a/src/enzo/Grid_InitializeRadiativeTransferFields.C +++ b/src/enzo/Grid_InitializeRadiativeTransferFields.C @@ -33,33 +33,39 @@ int grid::InitializeRadiativeTransferFields() // if (HasRadiation == FALSE) // return SUCCESS; - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); int i,j,k, index; /* Initialize photo and heating rates and compute number densities */ + for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { index = (k*GridDimension[1] + j)*GridDimension[0] + GridStartIndex[0]; - for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) { + for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) BaryonField[kphHINum][index] = - BaryonField[gammaHINum][index] = - BaryonField[kphHeINum][index] = - BaryonField[gammaHeINum][index] = - BaryonField[kphHeIINum][index] = - BaryonField[gammaHeIINum][index] = 0.; - if (MultiSpecies > 1 && !RadiativeTransferOpticallyThinH2) - BaryonField[kdissH2INum][index] = 0.; - } + BaryonField[gammaNum][index] = 0.0; } // loop over grid + if (RadiativeTransferHydrogenOnly == FALSE) + for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) + for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { + index = (k*GridDimension[1] + j)*GridDimension[0] + GridStartIndex[0]; + for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) + BaryonField[kphHeINum][index] = + BaryonField[kphHeIINum][index] = 0.0; + } // loop over grid + + if (MultiSpecies > 1 && !RadiativeTransferOpticallyThinH2) + for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) + for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) { + index = (k*GridDimension[1] + j)*GridDimension[0] + GridStartIndex[0]; + for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) + BaryonField[kdissH2INum][index] = 0.; + } // loop over grid + if (RadiationPressure) { int RPresNum1, RPresNum2, RPresNum3; diff --git a/src/enzo/Grid_PhotonTestInitializeGrid.C b/src/enzo/Grid_PhotonTestInitializeGrid.C index 12fb0572d..214f4d280 100644 --- a/src/enzo/Grid_PhotonTestInitializeGrid.C +++ b/src/enzo/Grid_PhotonTestInitializeGrid.C @@ -79,8 +79,8 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, int dim, i, j, k, m, field, sphere, size, index; int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, - DINum, DIINum, HDINum, kphHINum, gammaHINum, kphHeINum, gammaHeINum, - kphHeIINum, gammaHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; + DINum, DIINum, HDINum, kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; /* create fields */ @@ -126,11 +126,11 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, if (RadiativeTransfer) if (MultiSpecies) { FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; - FieldType[gammaHINum = NumberOfBaryonFields++] = gammaHI; - FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; - FieldType[gammaHeINum = NumberOfBaryonFields++] = gammaHeI; - FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; - FieldType[gammaHeIINum= NumberOfBaryonFields++] = gammaHeII; + FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; + if (RadiativeTransferHydrogenOnly == FALSE) { + FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; + FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; + } if (MultiSpecies > 1) FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; } diff --git a/src/enzo/Grid_Shine.C b/src/enzo/Grid_Shine.C index 393cb5e30..614aaf371 100644 --- a/src/enzo/Grid_Shine.C +++ b/src/enzo/Grid_Shine.C @@ -48,6 +48,7 @@ int grid::Shine(RadiationSourceEntry *RadiationSource) stype = 1; #endif if (MultiSpecies>1 && !RadiativeTransferOpticallyThinH2) stype++; + if (RadiativeTransferHydrogenOnly == TRUE) stype = 1; /* At most how many new Photon Packages should be allocated and created? */ diff --git a/src/enzo/Grid_SolveCoupledRateEquations.C b/src/enzo/Grid_SolveCoupledRateEquations.C index 67abd468f..6050b2f05 100644 --- a/src/enzo/Grid_SolveCoupledRateEquations.C +++ b/src/enzo/Grid_SolveCoupledRateEquations.C @@ -76,8 +76,9 @@ extern "C" void FORTRAN_NAME(solve_rate_cool)( float *inutot, int *iradtype, int *nfreq, int *imetalregen, int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p, int *iradtrans, int *iradcoupled, int *iradstep, int *ierr, + int *irt_honly, float *kphHI, float *kphHeI, float *kphHeII, - float *kdissH2I, float *gammaHI, float *gammaHeI, float *gammaHeII, + float *kdissH2I, float *photogamma, int *icmbTfloor, int *iClHeat, int *iClMMW, float *clMetNorm, float *clEleFra, int *clGridRank, int *clGridDim, float *clPar1, float *clPar2, float *clPar3, float *clPar4, float *clPar5, @@ -123,13 +124,9 @@ int grid::SolveCoupledRateEquations() /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Get easy to handle pointers for each variable. */ @@ -267,9 +264,9 @@ int grid::SolveCoupledRateEquations() &RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor, &RadiativeTransfer, &RadiativeTransferCoupledRateSolver, &RTCoupledSolverIntermediateStep, &ierr, + &RadiativeTransferHydrogenOnly, BaryonField[kphHINum], BaryonField[kphHeINum], BaryonField[kphHeIINum], - BaryonField[kdissH2INum], BaryonField[gammaHINum], BaryonField[gammaHeINum], - BaryonField[gammaHeIINum], + BaryonField[kdissH2INum], BaryonField[gammaNum], &CloudyCoolingData.CMBTemperatureFloor, &CloudyCoolingData.IncludeCloudyHeating, &CloudyCoolingData.IncludeCloudyMMW, &CloudyCoolingData.CloudyMetallicityNormalization, diff --git a/src/enzo/Grid_SolveHighDensityPrimordialChemistry.C b/src/enzo/Grid_SolveHighDensityPrimordialChemistry.C index 43cf45e10..e0e603aac 100644 --- a/src/enzo/Grid_SolveHighDensityPrimordialChemistry.C +++ b/src/enzo/Grid_SolveHighDensityPrimordialChemistry.C @@ -121,12 +121,9 @@ int grid::SolveHighDensityPrimordialChemistry() /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - ENZO_FAIL("Error in grid->IdentifyRadiativeTransferFields."); -} + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Compute size of the current grid. */ diff --git a/src/enzo/Grid_SolveRadiativeCooling.C b/src/enzo/Grid_SolveRadiativeCooling.C index 1fe861646..fce1b1dd2 100644 --- a/src/enzo/Grid_SolveRadiativeCooling.C +++ b/src/enzo/Grid_SolveRadiativeCooling.C @@ -76,7 +76,7 @@ extern "C" void FORTRAN_NAME(multi_cool)( float *metala, int *n_xe, float *xe_start, float *xe_end, float *inutot, int *iradtype, int *nfreq, int *imetalregen, int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p, - int *iradtrans, float *gammaHI, float *gammaHeI, float *gammaHeII, + int *iradtrans, float *photogamma, int *icmbTfloor, int *iClHeat, int *iClMMW, float *clMetNorm, float *clEleFra, int *clGridRank, int *clGridDim, float *clPar1, float *clPar2, float *clPar3, float *clPar4, float *clPar5, @@ -145,13 +145,9 @@ int grid::SolveRadiativeCooling() /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Get easy to handle pointers for each variable. */ @@ -287,8 +283,7 @@ int grid::SolveRadiativeCooling() &RadiationData.NumberOfFrequencyBins, &RadiationFieldRecomputeMetalRates, &RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor, - &RadiativeTransfer, BaryonField[gammaHINum], BaryonField[gammaHeINum], - BaryonField[gammaHeIINum], + &RadiativeTransfer, BaryonField[gammaNum], &CloudyCoolingData.CMBTemperatureFloor, &CloudyCoolingData.IncludeCloudyHeating, &CloudyCoolingData.IncludeCloudyMMW, &CloudyCoolingData.CloudyMetallicityNormalization, diff --git a/src/enzo/Grid_SolveRateAndCoolEquations.C b/src/enzo/Grid_SolveRateAndCoolEquations.C index 3484783ae..c9261e23b 100644 --- a/src/enzo/Grid_SolveRateAndCoolEquations.C +++ b/src/enzo/Grid_SolveRateAndCoolEquations.C @@ -76,8 +76,9 @@ extern "C" void FORTRAN_NAME(solve_rate_cool)( float *inutot, int *iradtype, int *nfreq, int *imetalregen, int *iradshield, float *avgsighp, float *avgsighep, float *avgsighe2p, int *iradtrans, int *iradcoupled, int *iradstep, int *ierr, + int *irt_honly, float *kphHI, float *kphHeI, float *kphHeII, - float *kdissH2I, float *gammaHI, float *gammaHeI, float *gammaHeII, + float *kdissH2I, float *photogamma, int *icmbTfloor, int *iClHeat, int *iClMMW, float *clMetNorm, float *clEleFra, int *clGridRank, int *clGridDim, float *clPar1, float *clPar2, float *clPar3, float *clPar4, float *clPar5, @@ -124,12 +125,9 @@ int grid::SolveRateAndCoolEquations() /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - ENZO_FAIL("Error in grid->IdentifyRadiativeTransferFields."); -} + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Compute size of the current grid. */ @@ -284,9 +282,9 @@ int grid::SolveRateAndCoolEquations() &RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor, &RadiativeTransfer, &RadiativeTransferCoupledRateSolver, &RTCoupledSolverIntermediateStep, &ierr, + &RadiativeTransferHydrogenOnly, BaryonField[kphHINum], BaryonField[kphHeINum], BaryonField[kphHeIINum], - BaryonField[kdissH2INum], BaryonField[gammaHINum], BaryonField[gammaHeINum], - BaryonField[gammaHeIINum], + BaryonField[kdissH2INum], BaryonField[gammaNum], &CloudyCoolingData.CMBTemperatureFloor, &CloudyCoolingData.IncludeCloudyHeating, &CloudyCoolingData.IncludeCloudyMMW, &CloudyCoolingData.CloudyMetallicityNormalization, diff --git a/src/enzo/Grid_SolveRateEquations.C b/src/enzo/Grid_SolveRateEquations.C index 14e10581f..a2c127fd4 100644 --- a/src/enzo/Grid_SolveRateEquations.C +++ b/src/enzo/Grid_SolveRateEquations.C @@ -53,8 +53,8 @@ extern "C" void FORTRAN_NAME(solve_rate)( float *HM, float *H2I, float *H2II, float *DI, float *DII, float *HDI, int *iradshield, float *avgsigh, float *avgsighe, float *avgsighe2, int *iradfield, float *piHI, float *piHeI, - int *iradtrans, float *kphHI, float *kphHeI, float *kphHeII, - float *kdissH2I); + int *iradtrans, int *irt_honly, float *kphHI, float *kphHeI, + float *kphHeII, float *kdissH2I); int grid::SolveRateEquations() @@ -91,13 +91,9 @@ int grid::SolveRateEquations() /* Find photo-ionization fields */ int kphHINum, kphHeINum, kphHeIINum, kdissH2INum; - int gammaHINum, gammaHeINum, gammaHeIINum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stderr, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int gammaNum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); /* Find the density field. */ @@ -197,7 +193,8 @@ int grid::SolveRateEquations() BaryonField[DINum], BaryonField[DIINum], BaryonField[HDINum], &RadiationShield, &HIShieldFactor, &HeIShieldFactor, &HeIIShieldFactor, &RadiationFieldType, &CoolData.piHI, &CoolData.piHeI, - &RadiativeTransfer, BaryonField[kphHINum], BaryonField[kphHeINum], + &RadiativeTransfer, &RadiativeTransferHydrogenOnly, + BaryonField[kphHINum], BaryonField[kphHeINum], BaryonField[kphHeIINum], BaryonField[kdissH2INum]); /* deallocate temporary space for solver */ diff --git a/src/enzo/Grid_StarParticleHandler.C b/src/enzo/Grid_StarParticleHandler.C index 07fe0e618..8b3b76709 100644 --- a/src/enzo/Grid_StarParticleHandler.C +++ b/src/enzo/Grid_StarParticleHandler.C @@ -448,7 +448,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level) h2field[index] = BaryonField[H2INum][index] + BaryonField[H2IINum][index]; } } - printf("Star type \n"); + //printf("Star type \n"); /* Set the units. */ float DensityUnits = 1, LengthUnits = 1, TemperatureUnits = 1, diff --git a/src/enzo/Grid_TransportPhotonPackages.C b/src/enzo/Grid_TransportPhotonPackages.C index c61333f8f..0e1c45341 100644 --- a/src/enzo/Grid_TransportPhotonPackages.C +++ b/src/enzo/Grid_TransportPhotonPackages.C @@ -93,23 +93,13 @@ int grid::TransportPhotonPackages(int level, ListOfPhotonsToMove **PhotonsToMove /* Find radiative transfer fields. */ - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { - fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); - ENZO_FAIL(""); - } + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum); int RPresNum1, RPresNum2, RPresNum3; - if (RadiationPressure) { - if (IdentifyRadiationPressureFields(RPresNum1, RPresNum2, RPresNum3) - == FAIL) { - fprintf(stdout, "Error in IdentifyRadiationPressureFields.\n"); - ENZO_FAIL(""); - } - } + if (RadiationPressure) + IdentifyRadiationPressureFields(RPresNum1, RPresNum2, RPresNum3); /* Get units. */ @@ -183,8 +173,8 @@ int grid::TransportPhotonPackages(int level, ListOfPhotonsToMove **PhotonsToMove WalkPhotonPackage(&PP, &MoveToGrid, ParentGrid, CurrentGrid, Grids0, nGrids0, DensNum, HINum, HeINum, HeIINum, H2INum, - kphHINum, gammaHINum, kphHeINum, gammaHeINum, - kphHeIINum, gammaHeIINum, kdissH2INum, RPresNum1, + kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3, DeleteMe, PauseMe, DeltaLevel, LightCrossingTime, DensityUnits, TemperatureUnits, VelocityUnits, diff --git a/src/enzo/Grid_WalkPhotonPackage.C b/src/enzo/Grid_WalkPhotonPackage.C index 13f03500c..b77fd56ce 100644 --- a/src/enzo/Grid_WalkPhotonPackage.C +++ b/src/enzo/Grid_WalkPhotonPackage.C @@ -51,8 +51,8 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, grid **MoveToGrid, grid *ParentGrid, grid *CurrentGrid, grid **Grids0, int nGrids0, int DensNum, int HINum, int HeINum, int HeIINum, int H2INum, - int kphHINum, int gammaHINum, int kphHeINum, - int gammaHeINum, int kphHeIINum, int gammaHeIINum, + int kphHINum, int gammaNum, int kphHeINum, + int kphHeIINum, int kdissH2INum, int RPresNum1, int RPresNum2, int RPresNum3, int &DeleteMe, int &PauseMe, int &DeltaLevel, float LightCrossingTime, @@ -66,7 +66,6 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, const float PopulationFractions[] = {1.0, 0.25, 0.25, 1.0}; const float EscapeRadiusFractions[] = {0.5, 1.0, 2.0}; const int kphNum[] = {kphHINum, kphHeINum, kphHeIINum}; - const int gammaNum[] = {gammaHINum, gammaHeINum, gammaHeIINum}; float ConvertToProperNumberDensity = DensityUnits/1.673e-24f; @@ -492,9 +491,9 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, // contributions to the photoionization rate is over whole timestep BaryonField[kphNum[type]][index] += dP1*factor1; - // the heating rate is just the number of photo ionizations times - // the excess energy units here are eV/s/cm^3 *TimeUnits. - BaryonField[gammaNum[type]][index] += dP1*factor2[0]; + // the heating rate is just the number of photo ionizations + // times the excess energy units here are eV/s *TimeUnits. + BaryonField[gammaNum][index] += dP1*factor2[0]; break; @@ -571,8 +570,8 @@ int grid::WalkPhotonPackage(PhotonPackageEntry **PP, BaryonField[kphNum[i]][index] += dP1 * factor1 * ion2_factor[i]; // the heating rate is just the number of photo ionizations times - // the excess energy units here are eV/s/cm^3 *TimeUnits. - BaryonField[gammaNum[i]][index] += dP1 * factor2[i] * heat_factor; + // the excess energy units here are eV/s/cm^3 *TimeUnits. + BaryonField[gammaNum][index] += dP1 * factor2[i] * heat_factor; } // ENDFOR absorber diff --git a/src/enzo/Make.config.objects b/src/enzo/Make.config.objects index 250f29d52..5a4b39bad 100644 --- a/src/enzo/Make.config.objects +++ b/src/enzo/Make.config.objects @@ -128,6 +128,7 @@ OBJS_CONFIG_LIB = \ ExposeGridHierarchy.o \ ExternalBoundary_AppendForcingToBaryonFields.o \ ExternalBoundary_constructor.o \ + ExternalBoundary_DeleteObsoleteFields.o \ ExternalBoundary_DetachForcingFromBaryonFields.o \ ExternalBoundary_IdentifyPhysicalQuantities.o \ ExternalBoundary_InitializeExternalBoundaryFaceIO.o \ @@ -280,6 +281,7 @@ OBJS_CONFIG_LIB = \ Grid_FlagCellsToBeRefinedByMetallicity.o \ Grid_DepositMustRefineParticles.o \ Grid_FlagCellsToBeRefinedByMustRefineRegion.o \ + Grid_DeleteObsoleteFields.o \ Grid_DepositBaryons.o \ Grid_DepositParticlePositions.o \ Grid_DepositParticlePositionsLocal.o \ diff --git a/src/enzo/PhotonGrid_Methods.h b/src/enzo/PhotonGrid_Methods.h index 69b353ae4..ca72c3ac5 100644 --- a/src/enzo/PhotonGrid_Methods.h +++ b/src/enzo/PhotonGrid_Methods.h @@ -202,11 +202,9 @@ int CountRadiationCells(void) { ncells = 0; - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + if (IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum) == FAIL) { fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); return 0; } @@ -229,11 +227,9 @@ float Max_kph(int &ncells) { ncells = 0; - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + if (IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum) == FAIL) { fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); return 0; } @@ -259,11 +255,9 @@ float Min_kph(int &ncells) { ncells = 0; - int kphHINum, gammaHINum, kphHeINum, gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum; - if (IdentifyRadiativeTransferFields(kphHINum, gammaHINum, kphHeINum, - gammaHeINum, kphHeIINum, gammaHeIINum, - kdissH2INum) == FAIL) { + int kphHINum, gammaNum, kphHeINum, kphHeIINum, kdissH2INum; + if (IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum) == FAIL) { fprintf(stdout, "Error in grid->IdentifyRadiativeTransferFields.\n"); return 0; } @@ -314,9 +308,9 @@ int WalkPhotonPackage(PhotonPackageEntry **PP, grid **Grids0, int nGrids0, int DensNum, int HINum, int HeINum, int HeIINum, int H2INum, - int kphHINum, int gammaHINum, - int kphHeINum, int gammaHeINum, - int kphHeIINum, int gammaHeIINum, + int kphHINum, int gammaNum, + int kphHeINum, + int kphHeIINum, int kdissH2INum, int RPresNum1, int RPresNum2, int RPresNum3, int &DeleteMe, int &PauseMe, int &DeltaLevel, float LightCrossingTime, diff --git a/src/enzo/PhotonTestInitialize.C b/src/enzo/PhotonTestInitialize.C index a4e1ba6f5..4725887ea 100644 --- a/src/enzo/PhotonTestInitialize.C +++ b/src/enzo/PhotonTestInitialize.C @@ -71,11 +71,9 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, const char *DIIName = "DII_Density"; const char *HDIName = "HDI_Density"; const char *kphHIName = "HI_kph"; - const char *gammaHIName = "HI_gamma"; + const char *gammaName = "PhotoGamma"; const char *kphHeIName = "HeI_kph"; - const char *gammaHeIName = "HeI_gamma"; const char *kphHeIIName = "HeII_kph"; - const char *gammaHeIIName= "HeII_gamma"; const char *kdissH2IName = "H2I_kdiss"; const char *RadAccel1Name = "x-RadPressure"; const char *RadAccel2Name = "y-RadPressure"; @@ -453,11 +451,9 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, if (RadiativeTransfer) if (MultiSpecies) { DataLabel[count++] = (char*) kphHIName; - DataLabel[count++] = (char*) gammaHIName; + DataLabel[count++] = (char*) gammaName; DataLabel[count++] = (char*) kphHeIName; - DataLabel[count++] = (char*) gammaHeIName; DataLabel[count++] = (char*) kphHeIIName; - DataLabel[count++] = (char*) gammaHeIIName; if (MultiSpecies > 1) DataLabel[count++]= (char*) kdissH2IName; } // if RadiativeTransfer diff --git a/src/enzo/RadiativeTransferInitialize.C b/src/enzo/RadiativeTransferInitialize.C index 69186bae3..d9504c7a9 100644 --- a/src/enzo/RadiativeTransferInitialize.C +++ b/src/enzo/RadiativeTransferInitialize.C @@ -38,11 +38,9 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, { const char *kphHIName = "HI_kph"; - const char *gammaHIName = "HI_gamma"; + const char *gammaName = "PhotoGamma"; const char *kphHeIName = "HeI_kph"; - const char *gammaHeIName = "HeI_gamma"; const char *kphHeIIName = "HeII_kph"; - const char *gammaHeIIName = "HeII_gamma"; const char *kdissH2IName = "H2I_kdiss"; const char *RadAccel1Name = "RadAccel1"; const char *RadAccel2Name = "RadAccel2"; @@ -101,8 +99,12 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, ExistingTypes[i] = FieldUndefined; if (RadiativeTransfer) { - for (i = kphHI; i <= gammaHeII; i++) - TypesToAdd[FieldsToAdd++] = i; + TypesToAdd[FieldsToAdd++] = kphHI; + TypesToAdd[FieldsToAdd++] = PhotoGamma; + if (RadiativeTransferHydrogenOnly == FALSE) { + TypesToAdd[FieldsToAdd++] = kphHeI; + TypesToAdd[FieldsToAdd++] = kphHeII; + } if (MultiSpecies > 1) TypesToAdd[FieldsToAdd++] = kdissH2I; if (RadiationPressure) @@ -146,10 +148,7 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel) - if (Temp->GridData->AddFields(TypesToAdd, FieldsToAdd) == FAIL) { - fprintf(stderr, "Error in grid::AddFields\n"); - ENZO_FAIL(""); - } + Temp->GridData->AddFields(TypesToAdd, FieldsToAdd); /* Add external boundaries */ @@ -157,6 +156,27 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, Exterior.AddField(TypesToAdd[i]); } // ENDFOR fields + /* Check for old gammaHeI and gammaHeII fields. Delete if they + exist. */ + + int NumberOfObsoleteFields = 2; + int ObsoleteFields[MAX_NUMBER_OF_BARYON_FIELDS]; + + ObsoleteFields[0] = gammaHeI; + ObsoleteFields[1] = gammaHeII; + if (RadiativeTransferHydrogenOnly == TRUE) { + NumberOfObsoleteFields += 2; + ObsoleteFields[2] = kphHeI; + ObsoleteFields[3] = kphHeII; + } + + for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) + for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel) + Temp->GridData->DeleteObsoleteFields(ObsoleteFields, + NumberOfObsoleteFields); + + Exterior.DeleteObsoleteFields(ObsoleteFields, NumberOfObsoleteFields); + /* Assign the radiation field DataLabels */ for (i = 0; i < FieldsToAdd; i++) { @@ -164,21 +184,15 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, case kphHI: DataLabel[OldNumberOfBaryonFields+i] = (char*) kphHIName; break; - case gammaHI: - DataLabel[OldNumberOfBaryonFields+i] = (char*) gammaHIName; + case PhotoGamma: + DataLabel[OldNumberOfBaryonFields+i] = (char*) gammaName; break; case kphHeI: DataLabel[OldNumberOfBaryonFields+i] = (char*) kphHeIName; break; - case gammaHeI: - DataLabel[OldNumberOfBaryonFields+i] = (char*) gammaHeIName; - break; case kphHeII: DataLabel[OldNumberOfBaryonFields+i] = (char*) kphHeIIName; break; - case gammaHeII: - DataLabel[OldNumberOfBaryonFields+i] = (char*) gammaHeIIName; - break; case kdissH2I: DataLabel[OldNumberOfBaryonFields+i] = (char*) kdissH2IName; break; @@ -217,4 +231,4 @@ int RadiativeTransferInitialize(char *ParameterFile, TopGridData &MetaData, return SUCCESS; - } +} diff --git a/src/enzo/RadiativeTransferParameters.h b/src/enzo/RadiativeTransferParameters.h index 3f3e0165b..c168f0d4b 100644 --- a/src/enzo/RadiativeTransferParameters.h +++ b/src/enzo/RadiativeTransferParameters.h @@ -70,3 +70,5 @@ EXTERN int RadiativeTransferHIIRestrictedTimestep; EXTERN int RadiativeTransferAdaptiveTimestep; EXTERN float GlobalMaximumkphIfront; + +EXTERN int RadiativeTransferHydrogenOnly; diff --git a/src/enzo/RadiativeTransferReadParameters.C b/src/enzo/RadiativeTransferReadParameters.C index 9c223855f..744c1297c 100644 --- a/src/enzo/RadiativeTransferReadParameters.C +++ b/src/enzo/RadiativeTransferReadParameters.C @@ -59,6 +59,7 @@ int RadiativeTransferReadParameters(FILE *fptr) RadiativeTransferPeriodicBoundary = FALSE; RadiativeTransferHIIRestrictedTimestep = FALSE; RadiativeTransferAdaptiveTimestep = FALSE; + RadiativeTransferHydrogenOnly = FALSE; /* read input from file */ @@ -102,6 +103,8 @@ int RadiativeTransferReadParameters(FILE *fptr) &RadiativeTransferHIIRestrictedTimestep); ret += sscanf(line, "RadiativeTransferAdaptiveTimestep = %"ISYM, &RadiativeTransferAdaptiveTimestep); + ret += sscanf(line, "RadiativeTransferHydrogenOnly = %"ISYM, + &RadiativeTransferHydrogenOnly); ret += sscanf(line, "dtPhoton = %"FSYM, &dtPhoton); /* if the line is suspicious, issue a warning */ diff --git a/src/enzo/RadiativeTransferWriteParameters.C b/src/enzo/RadiativeTransferWriteParameters.C index 32ab6eff5..0ab0f53fa 100644 --- a/src/enzo/RadiativeTransferWriteParameters.C +++ b/src/enzo/RadiativeTransferWriteParameters.C @@ -59,8 +59,10 @@ int RadiativeTransferWriteParameters(FILE *fptr) RadiativeTransferPhotonMergeRadius); fprintf(fptr, "RadiativeTransferHIIRestrictedTimestep = %"ISYM"\n", RadiativeTransferHIIRestrictedTimestep); - fprintf(fptr, "RadiativeTransferAdaptiveTimestep = %"ISYM"\n\n", + fprintf(fptr, "RadiativeTransferAdaptiveTimestep = %"ISYM"\n", RadiativeTransferAdaptiveTimestep); + fprintf(fptr, "RadiativeTransferHydrogenOnly = %"ISYM"\n\n", + RadiativeTransferHydrogenOnly); return SUCCESS; } diff --git a/src/enzo/StarParticleAddFeedback.C b/src/enzo/StarParticleAddFeedback.C index c39c4a2e9..6cfdafe7e 100644 --- a/src/enzo/StarParticleAddFeedback.C +++ b/src/enzo/StarParticleAddFeedback.C @@ -96,7 +96,7 @@ int StarParticleAddFeedback(TopGridData *MetaData, (!cstar->ApplyFeedbackTrue(SNe_dt))) continue; - float dtForThisStar = Temp->GridData->ReturnTimeStep(); + float dtForThisStar = LevelArray[level]->GridData->ReturnTimeStep(); /* Compute some parameters */ cstar->CalculateFeedbackParameters(influenceRadius, RootCellWidth, diff --git a/src/enzo/cool1d_multi.src b/src/enzo/cool1d_multi.src index 7bb2fc95f..c5ac6b570 100644 --- a/src/enzo/cool1d_multi.src +++ b/src/enzo/cool1d_multi.src @@ -30,7 +30,7 @@ & tgas, tgasold, p2d, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, + & iradtrans, photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, @@ -75,7 +75,7 @@ real HM(in,jn,kn), H2I(in,jn,kn), H2II(in,jn,kn), & DI(in,jn,kn), DII(in,jn,kn), HDI(in,jn,kn), & metal(in,jn,kn) - real gammaHI(in,jn,kn), gammaHeI(in,jn,kn), gammaHeII(in,jn,kn) + real photogamma(in,jn,kn) real hyd01ka(nratec), h2k01a(nratec), vibha(nratec), & rotha(nratec), rotla(nratec), gpldla(nratec), & gphdla(nratec), hdltea(nratec), hdlowa(nratec), @@ -630,14 +630,11 @@ c for radiative transfer, convert from eV/s*TimeUnits to the coolunits use rtunits = 1.60217646e-12/utim/coolunit do i = is+1, ie+1 if (itmask(i)) then - edot(i) = edot(i) + float(ipiht)*( - & + gammaHI (i,j,k) * rtunits * HI (i,j,k) - & + gammaHeI (i,j,k) * rtunits * HeI (i,j,k)*0.25 - & + gammaHeII(i,j,k) * rtunits * HeII(i,j,k)*0.25 - & )/dom + edot(i) = edot(i) + float(ipiht) * photogamma(i,j,k) * + & rtunits if (edot(i) .ne. edot(i)) then write(6,*) 'NaN in edot[2]: ', i, j, k, edot(i), - & gammaHI(i,j,k), HI(i,j,k), de(i,j,k), d(i,j,k), + & photogamma(i,j,k),HI(i,j,k), de(i,j,k), d(i,j,k), & ge(i,j,k), p2d(i), tgas(i), dom, urho, aye, mh ERROR_MESSAGE endif diff --git a/src/enzo/cool1d_sep.src b/src/enzo/cool1d_sep.src index 23c478c7b..a35369916 100644 --- a/src/enzo/cool1d_sep.src +++ b/src/enzo/cool1d_sep.src @@ -27,7 +27,7 @@ c & tgas, tgasold, p2d, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, itmask) + & iradtrans, photogamma, itmask) c c CALCULATE COOLING LUMINOSTIIES c @@ -65,7 +65,7 @@ c real HM(in,jn,kn), H2I(in,jn,kn), H2II(in,jn,kn), & DI(in,jn,kn), DII(in,jn,kn), HDI(in,jn,kn), & metal(in,jn,kn) - real gammaHI(in,jn,kn), gammaHeI(in,jn,kn), gammaHeII(in,jn,kn) + real photogamma(in,jn,kn) real hyd01ka(nratec), h2k01a(nratec), vibha(nratec), & rotha(nratec), rotla(nratec), gpldla(nratec), & gphdla(nratec), hdltea(nratec), hdlowa(nratec), diff --git a/src/enzo/cool_multi_lum.src b/src/enzo/cool_multi_lum.src index 3e2a935f6..381a28006 100644 --- a/src/enzo/cool_multi_lum.src +++ b/src/enzo/cool_multi_lum.src @@ -21,7 +21,7 @@ c & metala, n_xe, xe_start, xe_end, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII) + & iradtrans, photogamma) c c SOLVE RADIATIVE COOLING/HEATING EQUATIONS c @@ -69,7 +69,7 @@ c & reHeII2a(nratec), reHeIIIa(nratec), brema(nratec) real compa, piHI, piHeI, piHeII, comp_xraya, comp_temp, & inutot(nfreq), avgsighp, avgsighep, avgsighe2p - real gammaHI(in,jn,kn), gammaHeI(in,jn,kn), gammaHeII(in,jn,kn) + real photogamma(in,jn,kn) c c Parameters c @@ -189,7 +189,7 @@ c & tgas, tgasold, p2d, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, itmask) + & iradtrans, photogamma, itmask) c c Compute the cooling luminosity on the slice. Multiply by c density because edot is specific energy (erg/s), not per unit diff --git a/src/enzo/cool_multi_time.src b/src/enzo/cool_multi_time.src index 0a60e6d95..2e72e7aea 100644 --- a/src/enzo/cool_multi_time.src +++ b/src/enzo/cool_multi_time.src @@ -22,7 +22,7 @@ & metala, n_xe, xe_start, xe_end, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, + & iradtrans, photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, @@ -65,7 +65,7 @@ real HM(in,jn,kn), H2I(in,jn,kn), H2II(in,jn,kn), & DI(in,jn,kn), DII(in,jn,kn), HDI(in,jn,kn) real metal(in,jn,kn) - real gammaHI(in,jn,kn), gammaHeI(in,jn,kn), gammaHeII(in,jn,kn) + real photogamma(in,jn,kn) real hyd01ka(nratec), h2k01a(nratec), vibha(nratec), & rotha(nratec), rotla(nratec), gpldla(nratec), & gphdla(nratec), hdltea(nratec), hdlowa(nratec) @@ -195,7 +195,7 @@ & tgas, tgasold, p2d, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, + & iradtrans, photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, diff --git a/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C b/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C index 9adb96b29..763224fbe 100644 --- a/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C @@ -37,8 +37,8 @@ int grid::Collapse3DInitializeGrid(int n_sphere, int dim, i, j, k, m, field, sphere, size; int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, - DINum, DIINum, HDINum, kphHINum, gammaHINum, kphHeINum, gammaHeINum, - kphHeIINum, gammaHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; + DINum, DIINum, HDINum, kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; NumberOfBaryonFields = 0; @@ -84,11 +84,9 @@ int grid::Collapse3DInitializeGrid(int n_sphere, if (RadiativeTransfer) if (MultiSpecies) { FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; - FieldType[gammaHINum = NumberOfBaryonFields++] = gammaHI; + FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; - FieldType[gammaHeINum = NumberOfBaryonFields++] = gammaHeI; FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; - FieldType[gammaHeIINum= NumberOfBaryonFields++] = gammaHeII; if (MultiSpecies > 1) { FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; } diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C index d5b5b1dda..17614dba1 100644 --- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C @@ -40,8 +40,8 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL int dim, i, j, k, m, n, field, sphere, size, igrid, activesize; int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum, - DINum, DIINum, HDINum, kphHINum, gammaHINum, kphHeINum, gammaHeINum, - kphHeIINum, gammaHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; + DINum, DIINum, HDINum, kphHINum, gammaNum, kphHeINum, + kphHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; NumberOfBaryonFields = 0; @@ -91,11 +91,9 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL if (RadiativeTransfer) { if (MultiSpecies) { FieldType[kphHINum = NumberOfBaryonFields++] = kphHI; - FieldType[gammaHINum = NumberOfBaryonFields++] = gammaHI; + FieldType[gammaNum = NumberOfBaryonFields++] = PhotoGamma; FieldType[kphHeINum = NumberOfBaryonFields++] = kphHeI; - FieldType[gammaHeINum = NumberOfBaryonFields++] = gammaHeI; FieldType[kphHeIINum = NumberOfBaryonFields++] = kphHeII; - FieldType[gammaHeIINum= NumberOfBaryonFields++] = gammaHeII; if (MultiSpecies > 1) { FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; } diff --git a/src/enzo/multi_cool.src b/src/enzo/multi_cool.src index 2aa147ed5..351a6488b 100644 --- a/src/enzo/multi_cool.src +++ b/src/enzo/multi_cool.src @@ -24,8 +24,7 @@ & metala, n_xe, xe_start, xe_end, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, kphHI, kphHeI, kphHeII, kdissH2I, - & gammaHI, gammaHeI, gammaHeII, + & iradtrans, photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, @@ -69,9 +68,7 @@ real HM(in,jn,kn), H2I(in,jn,kn), H2II(in,jn,kn), & DI(in,jn,kn), DII(in,jn,kn), HDI(in,jn,kn), & metal(in,jn,kn) - real kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn), - & kdissH2I(in,jn,kn) - real gammaHI(in,jn,kn), gammaHeI(in,jn,kn), gammaHeII(in,jn,kn) + real photogamma(in,jn,kn) real hyd01ka(nratec), h2k01a(nratec), vibha(nratec), & rotha(nratec), rotla(nratec), gpldla(nratec), & gphdla(nratec), hdltea(nratec), hdlowa(nratec) @@ -233,7 +230,7 @@ & tgas, tgasold, p2d, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, + & iradtrans, photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, diff --git a/src/enzo/solve_rate.src b/src/enzo/solve_rate.src index 5896d5d91..9c6c08fb6 100644 --- a/src/enzo/solve_rate.src +++ b/src/enzo/solve_rate.src @@ -18,7 +18,8 @@ c & HM, H2I, H2II, DI, DII, HDI, & iradshield, avgsigh, avgsighe, avgsighe2, & iradtype, piHI, piHeI, - & iradtrans, kphHI, kphHeI, kphHeII, kdissH2I) + & iradtrans, irt_honly, kphHI, kphHeI, kphHeII, + & kdissH2I) c c SOLVE MULTI-SPECIES RATE EQUATIONS c @@ -44,7 +45,7 @@ c c Arguments c integer in, jn, kn, is, js, ks, ie, je, ke, nratec, ispecies, - & iradshield, iradtype, iradtrans + & iradshield, iradtype, iradtrans, irt_honly real dt, aye, temstart, temend, utim, uaye, urho, uxyz, & fh, dtoh real de(in,jn,kn), HI(in,jn,kn), HII(in,jn,kn), @@ -414,12 +415,19 @@ c c Add photo-ionization rates if necessary c if (iradtrans.eq.1) then - do i = is+1, ie+1 - HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) - dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) - $ + kphHeI(i,j,k) * HeI(i,j,k) / 4.0 - $ + kphHeII(i,j,k) *HeII(i,j,k) / 4.0 - enddo + if (irt_honly.eq.0) then + do i = is+1, ie+1 + HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) + dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) + $ + kphHeI(i,j,k) * HeI(i,j,k) / 4.0 + $ + kphHeII(i,j,k) *HeII(i,j,k) / 4.0 + enddo + else + do i = is+1, ie+1 + HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) + dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) + enddo + endif endif c c Find timestep that keeps relative changes below 10% @@ -519,10 +527,12 @@ c & + k25shield(i)*HeII(i,j,k)/4. & + k26shield(i)*HeI(i,j,k)/4. #endif /* RADIATION */ - if (iradtrans .eq. 1) + if (iradtrans .eq. 1 .and. irt_honly.eq.0) $ scoef = scoef + kphHI(i,j,k) * HI(i,j,k) $ + kphHeI(i,j,k) * HeI(i,j,k)/4. $ + kphHeII(i,j,k) * HeII(i,j,k)/4. + if (iradtrans .eq. 1 .and. irt_honly.eq.1) + $ scoef = scoef + kphHI(i,j,k) * HI(i,j,k) acoef = -(k1(i)*HI(i,j,k) - k2(i)*HII(i,j,k) & + k3(i)*HeI(i,j,k)/4. - k6(i)*HeIII(i,j,k)/4. & + k5(i)*HeII(i,j,k)/4. - k4(i)*HeII(i,j,k)/4. @@ -545,7 +555,8 @@ c #ifdef RADIATION & + k26shield(i) #endif /* RADIATION */ - if (iradtrans.eq.1) acoef = acoef + kphHeI(i,j,k) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ acoef = acoef + kphHeI(i,j,k) HeIp(i) = ( scoef*dtit(i) + HeI(i,j,k) ) & / ( 1. + acoef*dtit(i) ) c @@ -556,12 +567,14 @@ c #ifdef RADIATION & + k26shield(i)*HeIp(i) #endif /* RADIATION */ - if (iradtrans.eq.1) scoef = scoef + kphHeI(i,j,k)*HeIp(i) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ scoef = scoef + kphHeI(i,j,k)*HeIp(i) acoef = k4(i)*de(i,j,k) + k5(i)*de(i,j,k) #ifdef RADIATION & + k25shield(i) #endif /* RADIATION */ - if (iradtrans.eq.1) acoef = acoef + kphHeII(i,j,k) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ acoef = acoef + kphHeII(i,j,k) HeIIp(i) = ( scoef*dtit(i) + HeII(i,j,k) ) & / ( 1. + acoef*dtit(i) ) c @@ -571,7 +584,7 @@ c #ifdef RADIATION & + k25shield(i)*HeIIp(i) #endif /* RADIATION */ - if (iradtrans.eq.1) + if (iradtrans.eq.1 .and. irt_honly.eq.0) $ scoef = scoef + kphHeII(i,j,k) * HeIIp(i) acoef = k6(i)*de(i,j,k) HeIIIp(i) = ( scoef*dtit(i) + HeIII(i,j,k) ) @@ -643,10 +656,12 @@ c & + k25shield(i)*HeIIp(i)/4. & + k26shield(i)*HeIp(i)/4. #endif /* RADIATION */ - if (iradtrans .eq. 1) + if (iradtrans .eq. 1 .and. irt_honly.eq.0) $ scoef = scoef + kphHI(i,j,k) * HIp(i) $ + kphHeI(i,j,k) * HeIp(i)/4. $ + kphHeII(i,j,k) * HeIIp(i)/4. + if (iradtrans .eq. 1 .and. irt_honly.eq.1) + $ scoef = scoef + kphHI(i,j,k) * HIp(i) acoef = - (k1(i) *HI(i,j,k) - k2(i)*HII(i,j,k) & + k3(i) *HeI(i,j,k)/4. - k6(i)*HeIII(i,j,k)/4. & + k5(i) *HeII(i,j,k)/4.- k4(i)*HeII(i,j,k)/4. diff --git a/src/enzo/solve_rate_cool.src b/src/enzo/solve_rate_cool.src index 76cf12274..0902a4df4 100644 --- a/src/enzo/solve_rate_cool.src +++ b/src/enzo/solve_rate_cool.src @@ -29,8 +29,9 @@ & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, & iradtrans, iradcoupled, iradstep, ierr, + & irt_honly, & kphHI, kphHeI, kphHeII, kdissH2I, - & gammaHI, gammaHeI, gammaHeII, + & photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, @@ -125,7 +126,7 @@ integer in, jn, kn, is, js, ks, ie, je, ke, nratec, imethod, & idual, iexpand, ih2co, ipiht, ispecies, imetal, idim, & iradtype, nfreq, imetalregen, iradshield, iradtrans, - & iradcoupled, iradstep, n_xe, ierr, imcool + & iradcoupled, iradstep, n_xe, ierr, imcool, irt_honly real dt, aye, temstart, temend, eta1, eta2, gamma, & utim, uxyz, uaye, urho, utem, fh, dtoh, xe_start, xe_end @@ -143,7 +144,7 @@ real kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn), & kdissH2I(in,jn,kn) - real gammaHI(in,jn,kn), gammaHeI(in,jn,kn), gammaHeII(in,jn,kn) + real photogamma(in,jn,kn) ! Cooling tables (coolings rates as a function of temperature) @@ -341,7 +342,7 @@ & tgas, tgasold, p2d, & inutot, iradtype, nfreq, imetalregen, & iradshield, avgsighp, avgsighep, avgsighe2p, - & iradtrans, gammaHI, gammaHeI, gammaHeII, + & iradtrans, photogamma, & icmbTfloor, iClHeat, iClMMW, & clMetNorm, clEleFra, clGridRank, clGridDim, & clPar1, clPar2, clPar3, clPar4, clPar5, @@ -379,7 +380,7 @@ & k24, k25, k26, k27, k28, k29, k30, k31, & k50, k51, k52, k53, k54, k55, & k56, k24shield, k25shield, k26shield, - & iradtrans, kphHI, kphHeI, kphHeII, + & iradtrans, irt_honly, kphHI, kphHeI, kphHeII, & kdissH2I, itmask) ! Find timestep that keeps relative chemical changes below 10% @@ -549,8 +550,8 @@ & HIp, HIIp, HeIp, HeIIp, HeIIIp, dep, & HMp, H2Ip, H2IIp, DIp, DIIp, HDIp, & dedot_prev, HIdot_prev, - & iradtrans, kphHI, kphHeI, kphHeII, kdissH2I, - & itmask) + & iradtrans, irt_honly, kphHI, kphHeI, kphHeII, + & kdissH2I, itmask) ! Add the timestep to the elapsed time for each cell and find ! minimum elapsed time step in this row @@ -944,7 +945,7 @@ c ------------------------------------------------------------------- & k24, k25, k26, k27, k28, k29, k30, k31, & k50, k51, k52, k53, k54, k55, & k56, k24shield, k25shield, k26shield, - & iradtrans, kphHI, kphHeI, kphHeII, + & iradtrans, irt_honly, kphHI, kphHeI, kphHeII, & kdissH2I, itmask) ! ------------------------------------------------------------------- @@ -953,7 +954,7 @@ c ------------------------------------------------------------------- ! arguments - integer ispecies, is, ie, j, k, in, jn, kn, iradtrans + integer ispecies, is, ie, j, k, in, jn, kn, iradtrans, irt_honly real dedot(in), HIdot(in) logical itmask(in) @@ -1068,14 +1069,23 @@ c ------------------------------------------------------------------- ! Add photo-ionization rates if needed if (iradtrans .eq. 1) then - do i = is+1, ie+1 - if (itmask(i)) then - HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) - dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) - & + kphHeI(i,j,k) * HeI(i,j,k) / 4.0 - & + kphHeII(i,j,k) *HeII(i,j,k) / 4.0 - endif - enddo + if (irt_honly .eq. 0) then + do i = is+1, ie+1 + if (itmask(i)) then + HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) + dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) + & + kphHeI(i,j,k) * HeI(i,j,k) / 4.0 + & + kphHeII(i,j,k) *HeII(i,j,k) / 4.0 + endif + enddo + else + do i = is+1, ie+1 + if (itmask(i)) then + HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) + dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) + endif + enddo + endif endif return @@ -1097,15 +1107,15 @@ c ------------------------------------------------------------------- & HIp, HIIp, HeIp, HeIIp, HeIIIp, dep, & HMp, H2Ip, H2IIp, DIp, DIIp, HDIp, & dedot_prev, HIdot_prev, - & iradtrans, kphHI, kphHeI, kphHeII, kdissH2I, - & itmask) + & iradtrans, irt_honly, kphHI, kphHeI, kphHeII, + & kdissH2I, itmask) c ------------------------------------------------------------------- implicit NONE ! arguments - integer ispecies, in, jn, kn, is, ie, j, k, iradtrans + integer ispecies, in, jn, kn, is, ie, j, k, iradtrans, irt_honly real dtit(in), dedot_prev(in), HIdot_prev(in) logical itmask(in) @@ -1184,10 +1194,12 @@ c & + k25shield(i)*HeII(i,j,k)/4. & + k26shield(i)*HeI(i,j,k)/4. #endif /* RADIATION */ - if (iradtrans .eq. 1) + if (iradtrans .eq. 1 .and. irt_honly .eq. 0) & scoef = scoef + kphHI(i,j,k) * HI(i,j,k) & + kphHeI(i,j,k) * HeI(i,j,k)/4. & + kphHeII(i,j,k) * HeII(i,j,k)/4. + if (iradtrans .eq. 1 .and. irt_honly .eq. 1) + & scoef = scoef + kphHI(i,j,k) * HI(i,j,k) acoef = -(k1(i)*HI(i,j,k) - k2(i)*HII(i,j,k) & + k3(i)*HeI(i,j,k)/4. - k6(i)*HeIII(i,j,k)/4.0 & + k5(i)*HeII(i,j,k)/4. - k4(i)*HeII(i,j,k)/4.0) @@ -1211,7 +1223,8 @@ c #ifdef RADIATION & + k26shield(i) #endif /* RADIATION */ - if (iradtrans.eq.1) acoef = acoef + kphHeI(i,j,k) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ acoef = acoef + kphHeI(i,j,k) HeIp(i) = ( scoef*dtit(i) + HeI(i,j,k) ) & / ( 1. + acoef*dtit(i) ) @@ -1222,12 +1235,14 @@ c #ifdef RADIATION & + k26shield(i)*HeIp(i) #endif /* RADIATION */ - if (iradtrans.eq.1) scoef = scoef + kphHeI(i,j,k)*HeIp(i) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ scoef = scoef + kphHeI(i,j,k)*HeIp(i) acoef = k4(i)*de(i,j,k) + k5(i)*de(i,j,k) #ifdef RADIATION & + k25shield(i) #endif /* RADIATION */ - if (iradtrans.eq.1) acoef = acoef + kphHeII(i,j,k) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ acoef = acoef + kphHeII(i,j,k) HeIIp(i) = ( scoef*dtit(i) + HeII(i,j,k) ) & / ( 1. + acoef*dtit(i) ) @@ -1237,7 +1252,8 @@ c #ifdef RADIATION & + k25shield(i)*HeIIp(i) #endif /* RADIATION */ - if (iradtrans.eq.1) scoef = scoef + kphHeII(i,j,k) * HeIIp(i) + if (iradtrans.eq.1 .and. irt_honly.eq.0) + $ scoef = scoef + kphHeII(i,j,k) * HeIIp(i) acoef = k6(i)*de(i,j,k) HeIIIp(i) = ( scoef*dtit(i) + HeIII(i,j,k) ) & / ( 1. + acoef*dtit(i) ) @@ -1315,10 +1331,12 @@ c --- (C) Now do extra 3-species for molecular hydrogen --- & + k25shield(i)*HeIIp(i)/4. & + k26shield(i)*HeIp(i)/4. #endif /* RADIATION */ - if (iradtrans .eq. 1) + if (iradtrans .eq. 1 .and. irt_honly.eq.0) & scoef = scoef + kphHI(i,j,k) * HIp(i) & + kphHeI(i,j,k) * HeIp(i)/4. & + kphHeII(i,j,k) * HeIIp(i)/4. + if (iradtrans .eq. 1 .and. irt_honly.eq.1) + & scoef = scoef + kphHI(i,j,k) * HIp(i) acoef = - (k1(i) *HI(i,j,k) - k2(i)*HII(i,j,k) & + k3(i) *HeI(i,j,k)/4. - k6(i)*HeIII(i,j,k)/4. & + k5(i) *HeII(i,j,k)/4.- k4(i)*HeII(i,j,k)/4. diff --git a/src/enzo/typedefs.h b/src/enzo/typedefs.h index 81fd3b6c3..aa01df098 100644 --- a/src/enzo/typedefs.h +++ b/src/enzo/typedefs.h @@ -63,7 +63,7 @@ const field_type ExtraType0 = 21, ExtraType1 = 22, kphHI = 23, - gammaHI = 24, + PhotoGamma = 24, kphHeI = 25, gammaHeI = 26, kphHeII = 27, From e33479ed7d0fff83f84de6bcaf5e28b48d1a5cf9 Mon Sep 17 00:00:00 2001 From: John Wise Date: Fri, 13 Nov 2009 12:59:55 -0500 Subject: [PATCH 18/19] We need to call CommunicationCollectParticles before *and* after CTP in RebuildHierarchy. Plus, CTP assumes that particles are distributed across processors, so I added an argument to CCP to not sync NumberOfParticles if needed. --HG-- branch : week-of-code --- src/enzo/CommunicationCollectParticles.C | 34 ++++++++++----- src/enzo/CommunicationTransferParticles.C | 40 +++++++++-------- src/enzo/CommunicationTransferParticlesOpt.C | 46 ++++++++++---------- src/enzo/CommunicationTransferStars.C | 23 ++++++++++ src/enzo/CommunicationTransferStarsOpt.C | 28 +++++++++--- src/enzo/Grid_MoveAllParticles.C | 4 +- src/enzo/RebuildHierarchy.C | 26 ++++++++--- 7 files changed, 133 insertions(+), 68 deletions(-) diff --git a/src/enzo/CommunicationCollectParticles.C b/src/enzo/CommunicationCollectParticles.C index 47b99739f..7d65dc120 100644 --- a/src/enzo/CommunicationCollectParticles.C +++ b/src/enzo/CommunicationCollectParticles.C @@ -53,7 +53,7 @@ int CommunicationShareStars(int *NumberToMove, star_data* &SendList, int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], int level, bool ParticlesAreLocal, - int CollectMode) + bool SyncNumberOfParticles, int CollectMode) { /* Create pointer arrays and count grids */ @@ -117,15 +117,19 @@ int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], #ifdef DEBUG_CCP for (i = 0; i < NumberOfGrids; i++) - printf("CCP[P%"ISYM"B]: grid %"ISYM", %"ISYM" proc, %"ISYM" particles\n", + printf("CCP[P%"ISYM"B]: grid %"ISYM", %"ISYM" proc, %"ISYM" particles, " + "%"ISYM" stars\n", MyProcessorNumber, i, GridHierarchyPointer[i]->GridData->ReturnProcessorNumber(), - GridHierarchyPointer[i]->GridData->ReturnNumberOfParticles()); + GridHierarchyPointer[i]->GridData->ReturnNumberOfParticles(), + GridHierarchyPointer[i]->GridData->ReturnNumberOfStars()); for (i = 0; i < NumberOfSubgrids; i++) - printf("CCP[P%"ISYM"B]: subgrid %"ISYM", %"ISYM" proc, %"ISYM" particles\n", + printf("CCP[P%"ISYM"B]: subgrid %"ISYM", %"ISYM" proc, %"ISYM" particles, " + "%"ISYM" stars\n", MyProcessorNumber, i, SubgridHierarchyPointer[i]->GridData->ReturnProcessorNumber(), - SubgridHierarchyPointer[i]->GridData->ReturnNumberOfParticles()); + SubgridHierarchyPointer[i]->GridData->ReturnNumberOfParticles(), + SubgridHierarchyPointer[i]->GridData->ReturnNumberOfStars()); #endif /* DEBUG_CCP */ for (i = 0; i < NumberOfProcessors; i++) { @@ -246,7 +250,8 @@ int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], processor, set number of particles so everybody agrees. ************************************************************************/ - if ((!KeepLocal && NumberOfProcessors > 1) || ParticlesAreLocal) { + if ((!KeepLocal && NumberOfProcessors > 1) || + (ParticlesAreLocal && SyncNumberOfParticles)) { CommunicationSyncNumberOfParticles(GridHierarchyPointer, NumberOfGrids); CommunicationSyncNumberOfParticles(SubgridHierarchyPointer, NumberOfSubgrids); @@ -255,15 +260,19 @@ int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], #ifdef DEBUG_CCP for (i = 0; i < NumberOfGrids; i++) - printf("CCP[P%"ISYM"A]: grid %"ISYM", %"ISYM" proc, %"ISYM" particles\n", + printf("CCP[P%"ISYM"A]: grid %"ISYM", %"ISYM" proc, %"ISYM" particles, " + "%"ISYM" stars\n", MyProcessorNumber, i, GridHierarchyPointer[i]->GridData->ReturnProcessorNumber(), - GridHierarchyPointer[i]->GridData->ReturnNumberOfParticles()); + GridHierarchyPointer[i]->GridData->ReturnNumberOfParticles(), + GridHierarchyPointer[i]->GridData->ReturnNumberOfStars()); for (i = 0; i < NumberOfSubgrids; i++) - printf("CCP[P%"ISYM"A]: subgrid %"ISYM", %"ISYM" proc, %"ISYM" particles\n", + printf("CCP[P%"ISYM"A]: subgrid %"ISYM", %"ISYM" proc, %"ISYM" particles, " + "%"ISYM" stars\n", MyProcessorNumber, i, SubgridHierarchyPointer[i]->GridData->ReturnProcessorNumber(), - SubgridHierarchyPointer[i]->GridData->ReturnNumberOfParticles()); + SubgridHierarchyPointer[i]->GridData->ReturnNumberOfParticles(), + SubgridHierarchyPointer[i]->GridData->ReturnNumberOfStars()); #endif /* DEBUG_CCP */ /* Cleanup. */ @@ -444,8 +453,9 @@ int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], If the particles and stars are only on the grid's host processor, set number of particles so everybody agrees. ************************************************************************/ - - CommunicationSyncNumberOfParticles(GridHierarchyPointer, NumberOfGrids); + + if (SyncNumberOfParticles) + CommunicationSyncNumberOfParticles(GridHierarchyPointer, NumberOfGrids); #ifdef DEBUG_CCP for (i = StartGrid; i < EndGrid; i++) diff --git a/src/enzo/CommunicationTransferParticles.C b/src/enzo/CommunicationTransferParticles.C index 294392f7b..42b7f6df7 100644 --- a/src/enzo/CommunicationTransferParticles.C +++ b/src/enzo/CommunicationTransferParticles.C @@ -344,25 +344,27 @@ int CommunicationTransferParticles(grid *GridPointer[], int NumberOfGrids) } // end: if grid is on my processor /* Set number of particles so everybody agrees. */ - -// if (NumberOfProcessors > 1) { -// int *Changes = new int[NumberOfGrids]; -// for (j = 0; j < NumberOfGrids; j++) -// Changes[j] = 0; -// for (j = 0; j < NumberOfGrids; j++) -// for (i = 0; i < 6; i++) -// if (SharedList[j].ToGrid[i] != -1) { -// Changes[SharedList[j].FromGrid] -= SharedList[j].NumberToMove[i]; -// Changes[SharedList[j].ToGrid[i]] += SharedList[j].NumberToMove[i]; -// } -// for (j = 0; j < NumberOfGrids; j++) { -// if (GridPointer[j]->ReturnProcessorNumber() != MyProcessorNumber) -// GridPointer[j]->SetNumberOfParticles( -// GridPointer[j]->ReturnNumberOfParticles()+Changes[j]); -// // printf("Pb(%"ISYM") CTP grid[%"ISYM"] = %"ISYM"\n", MyProcessorNumber, j, GridPointer[j]->ReturnNumberOfParticles()); -// } -// delete [] Changes; -// } + +#ifdef UNUSED + if (NumberOfProcessors > 1) { + int *Changes = new int[NumberOfGrids]; + for (j = 0; j < NumberOfGrids; j++) + Changes[j] = 0; + for (j = 0; j < NumberOfGrids; j++) + for (i = 0; i < 6; i++) + if (SharedList[j].ToGrid[i] != -1) { + Changes[SharedList[j].FromGrid] -= SharedList[j].NumberToMove[i]; + Changes[SharedList[j].ToGrid[i]] += SharedList[j].NumberToMove[i]; + } + for (j = 0; j < NumberOfGrids; j++) { + if (GridPointer[j]->ReturnProcessorNumber() != MyProcessorNumber) + GridPointer[j]->SetNumberOfParticles( + GridPointer[j]->ReturnNumberOfParticles()+Changes[j]); + // printf("Pb(%"ISYM") CTP grid[%"ISYM"] = %"ISYM"\n", MyProcessorNumber, j, GridPointer[j]->ReturnNumberOfParticles()); + } + delete [] Changes; + } +#endif /* CleanUp. */ diff --git a/src/enzo/CommunicationTransferParticlesOpt.C b/src/enzo/CommunicationTransferParticlesOpt.C index 1f9837ac4..0aa072a7a 100644 --- a/src/enzo/CommunicationTransferParticlesOpt.C +++ b/src/enzo/CommunicationTransferParticlesOpt.C @@ -31,17 +31,18 @@ #include "TopGridData.h" #include "Hierarchy.h" #include "LevelHierarchy.h" +#include "CommunicationUtilities.h" void my_exit(int status); // function prototypes Eint32 compare_grid(const void *a, const void *b); int Enzo_Dims_create(int nnodes, int ndims, int *dims); -int CommunicationAllSumIntegerValues(int *Values, int Number); int CommunicationShareParticles(int *NumberToMove, particle_data* &SendList, int &NumberOfReceives, particle_data* &SharedList); +#define NO_DEBUG_CTP #define KEEP_PARTICLES_LOCAL int CommunicationTransferParticles(grid *GridPointer[], int NumberOfGrids) @@ -190,23 +191,26 @@ int CommunicationTransferParticles(grid *GridPointer[], int NumberOfGrids) /* Set number of particles so everybody agrees. */ - // (JHW, May 2009) No longer needed because we synchronize in - // CommunicationCollectParticles now. - -// if (NumberOfProcessors > 1) { -// int *AllNumberOfParticles = new int[NumberOfGrids]; -// for (j = 0; j < NumberOfGrids; j++) -// if (GridPointer[j]->ReturnProcessorNumber() == MyProcessorNumber) -// AllNumberOfParticles[j] = GridPointer[j]->ReturnNumberOfParticles(); -// else -// AllNumberOfParticles[j] = 0; -// -// CommunicationAllSumIntegerValues(AllNumberOfParticles, NumberOfGrids); -// for (j = 0; j < NumberOfGrids; j++) -// GridPointer[j]->SetNumberOfParticles(AllNumberOfParticles[j]); -// -// delete [] AllNumberOfParticles; -// } +#ifdef UNUSED + if (NumberOfProcessors > 1) { + int *AllNumberOfParticles = new int[NumberOfGrids]; + for (j = 0; j < NumberOfGrids; j++) + if (GridPointer[j]->ReturnProcessorNumber() == MyProcessorNumber) + AllNumberOfParticles[j] = GridPointer[j]->ReturnNumberOfParticles(); + else + AllNumberOfParticles[j] = 0; + + CommunicationAllSumValues(AllNumberOfParticles, NumberOfGrids); + printf("P%d: npart = %d %d (%d %d)\n", MyProcessorNumber, + GridPointer[0]->ReturnNumberOfParticles(), + GridPointer[1]->ReturnNumberOfParticles(), + AllNumberOfParticles[0], AllNumberOfParticles[1]); + for (j = 0; j < NumberOfGrids; j++) + GridPointer[j]->SetNumberOfParticles(AllNumberOfParticles[j]); + + delete [] AllNumberOfParticles; + } +#endif /* UNUSED */ /* Cleanup. */ @@ -216,11 +220,7 @@ int CommunicationTransferParticles(grid *GridPointer[], int NumberOfGrids) delete [] NumberToMove; delete [] GridMap; -#ifdef USE_MPI - int temp_int = TotalNumberToMove; - MPI_Reduce(&temp_int, &TotalNumberToMove, 1, IntDataType, MPI_SUM, - ROOT_PROCESSOR, MPI_COMM_WORLD); -#endif + CommunicationSumValues(&TotalNumberToMove, 1); if (debug) printf("CommunicationTransferParticles: moved = %"ISYM"\n", TotalNumberToMove); diff --git a/src/enzo/CommunicationTransferStars.C b/src/enzo/CommunicationTransferStars.C index a13489ba3..a8c54d205 100644 --- a/src/enzo/CommunicationTransferStars.C +++ b/src/enzo/CommunicationTransferStars.C @@ -335,6 +335,29 @@ int CommunicationTransferStars(grid *GridPointer[], int NumberOfGrids) LocalNumberToMove, LocalPointer, COPY_IN); } // end: if grid is on my processor + + /* Set number of particles so everybody agrees. */ + +#ifdef UNUSED + if (NumberOfProcessors > 1) { + int *Changes = new int[NumberOfGrids]; + for (j = 0; j < NumberOfGrids; j++) + Changes[j] = 0; + for (j = 0; j < NumberOfGrids; j++) + for (i = 0; i < 6; i++) + if (SharedList[j].ToGrid[i] != -1) { + Changes[SharedList[j].FromGrid] -= SharedList[j].NumberToMove[i]; + Changes[SharedList[j].ToGrid[i]] += SharedList[j].NumberToMove[i]; + } + for (j = 0; j < NumberOfGrids; j++) { + if (GridPointer[j]->ReturnProcessorNumber() != MyProcessorNumber) + GridPointer[j]->SetNumberOfStars( + GridPointer[j]->ReturnNumberOfStars()+Changes[j]); + // printf("Pb(%"ISYM") CTP grid[%"ISYM"] = %"ISYM"\n", MyProcessorNumber, j, GridPointer[j]->ReturnNumberOfStars()); + } + delete [] Changes; + } +#endif /* UNUSED */ /* CleanUp. */ diff --git a/src/enzo/CommunicationTransferStarsOpt.C b/src/enzo/CommunicationTransferStarsOpt.C index 18dec4d1a..4ad3b569b 100644 --- a/src/enzo/CommunicationTransferStarsOpt.C +++ b/src/enzo/CommunicationTransferStarsOpt.C @@ -32,13 +32,13 @@ #include "TopGridData.h" #include "Hierarchy.h" #include "LevelHierarchy.h" +#include "CommunicationUtilities.h" void my_exit(int status); // function prototypes Eint32 compare_star_grid(const void *a, const void *b); int Enzo_Dims_create(int nnodes, int ndims, int *dims); -int CommunicationAllSumIntegerValues(int *Values, int Number); // Remove define. This method will always be used. //#define KEEP_PARTICLES_LOCAL @@ -137,6 +137,26 @@ int CommunicationTransferStars(grid *GridPointer[], int NumberOfGrids) } // ENDFOR grids } // ENDIF NumberOfRecieves > 0 + + /* Set number of stars so everybody agrees. */ + +#ifdef UNUSED + if (NumberOfProcessors > 1) { + int *AllNumberOfStars = new int[NumberOfGrids]; + for (j = 0; j < NumberOfGrids; j++) + if (GridPointer[j]->ReturnProcessorNumber() == MyProcessorNumber) + AllNumberOfStars[j] = GridPointer[j]->ReturnNumberOfStars(); + else + AllNumberOfStars[j] = 0; + + CommunicationAllSumValues(AllNumberOfStars, NumberOfGrids); + for (j = 0; j < NumberOfGrids; j++) + GridPointer[j]->SetNumberOfStars(AllNumberOfStars[j]); + + delete [] AllNumberOfStars; + } +#endif + /* Cleanup. */ if (SendList != SharedList) @@ -145,11 +165,7 @@ int CommunicationTransferStars(grid *GridPointer[], int NumberOfGrids) delete [] NumberToMove; delete [] GridMap; -#ifdef USE_MPI - int temp_int = TotalNumberToMove; - MPI_Reduce(&temp_int, &TotalNumberToMove, 1, IntDataType, MPI_SUM, - ROOT_PROCESSOR, MPI_COMM_WORLD); -#endif + CommunicationSumValues(&TotalNumberToMove, 1); if (debug) printf("CommunicationTransferStars: moved = %"ISYM"\n", TotalNumberToMove); diff --git a/src/enzo/Grid_MoveAllParticles.C b/src/enzo/Grid_MoveAllParticles.C index e80af4c05..94feed243 100644 --- a/src/enzo/Grid_MoveAllParticles.C +++ b/src/enzo/Grid_MoveAllParticles.C @@ -52,8 +52,8 @@ int grid::MoveAllParticles(int NumberOfGrids, grid* FromGrid[]) /* Debugging info. */ - if (debug) printf("MoveAllParticles: %"ISYM" (before: ThisGrid = %"ISYM").\n", - TotalNumberOfParticles, NumberOfParticles); + if (debug1) printf("MoveAllParticles: %"ISYM" (before: ThisGrid = %"ISYM").\n", + TotalNumberOfParticles, NumberOfParticles); /* Allocate space for the particles. */ diff --git a/src/enzo/RebuildHierarchy.C b/src/enzo/RebuildHierarchy.C index ce3b98eef..4a89b0c90 100644 --- a/src/enzo/RebuildHierarchy.C +++ b/src/enzo/RebuildHierarchy.C @@ -55,7 +55,8 @@ int CommunicationTransferSubgridParticles(LevelHierarchyEntry *LevelArray[], int CommunicationTransferParticles(grid *GridPointer[], int NumberOfGrids); int CommunicationTransferStars(grid *GridPointer[], int NumberOfGrids); int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], int level, - bool ParticlesAreLocal, int CollectMode); + bool ParticlesAreLocal, + bool SyncNumberOfParticles, int CollectMode); int CommunicationSyncNumberOfParticles(HierarchyEntry *GridHierarchyPointer[], int NumberOfGrids); int FastSiblingLocatorInitialize(ChainingMeshStructure *Mesh, int Rank, @@ -94,7 +95,7 @@ int RebuildHierarchy(TopGridData *MetaData, if (debug) printf("RebuildHierarchy: level = %"ISYM"\n", level); ReportMemoryUsage("Rebuild pos 1"); - bool ParticlesAreLocal; + bool ParticlesAreLocal, SyncNumberOfParticles = true; int i, j, k, grids, grids2, subgrids, MoveParticles; int TotalFlaggedCells, FlaggedGrids; FLOAT ZeroVector[MAX_DIMENSION]; @@ -200,13 +201,14 @@ int RebuildHierarchy(TopGridData *MetaData, tt0 = ReturnWallTime(); if (level > MaximumStaticSubgridLevel) { ParticlesAreLocal = false; + SyncNumberOfParticles = false; CommunicationCollectParticles(LevelArray, level, ParticlesAreLocal, - SIBLINGS_ONLY); + SyncNumberOfParticles, SIBLINGS_ONLY); ParticlesAreLocal = true; + SyncNumberOfParticles = true; } tt1 = ReturnWallTime(); RHperf[2] += tt1-tt0; - /* --------------------------------------------------------------------- */ /* if this is level 0 then transfer particles between grids. */ @@ -226,6 +228,18 @@ int RebuildHierarchy(TopGridData *MetaData, CommunicationTransferParticles(GridPointer, grids); CommunicationTransferStars(GridPointer, grids); + /* We need to collect particles again */ + + if (level > MaximumStaticSubgridLevel) { + ParticlesAreLocal = false; + SyncNumberOfParticles = true; + CommunicationCollectParticles(LevelArray, level, ParticlesAreLocal, + SyncNumberOfParticles, SIBLINGS_ONLY); + ParticlesAreLocal = true; + SyncNumberOfParticles = true; + } + + } // ENDIF level 0 tt1 = ReturnWallTime(); RHperf[1] += tt1-tt0; @@ -376,7 +390,7 @@ int RebuildHierarchy(TopGridData *MetaData, tt0 = ReturnWallTime(); CommunicationCollectParticles(LevelArray, i, ParticlesAreLocal, - SUBGRIDS_LOCAL); + SyncNumberOfParticles, SUBGRIDS_LOCAL); tt1 = ReturnWallTime(); RHperf[7] += tt1-tt0; @@ -480,7 +494,7 @@ int RebuildHierarchy(TopGridData *MetaData, for (j = level; j <= MaximumStaticSubgridLevel+1; j++) if (LevelArray[j] != NULL) CommunicationCollectParticles(LevelArray, j, ParticlesAreLocal, - SIBLINGS_ONLY); + SyncNumberOfParticles, SIBLINGS_ONLY); tt1 = ReturnWallTime(); RHperf[14] += tt1-tt0; From 114a1bd084d97d1841d896c35e848b49c274fb09 Mon Sep 17 00:00:00 2001 From: John Wise Date: Fri, 13 Nov 2009 14:27:56 -0500 Subject: [PATCH 19/19] Forgot to reset NumberOfParticles=0 if we're not sync'ing it. --HG-- branch : week-of-code --- src/enzo/CommunicationCollectParticles.C | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/enzo/CommunicationCollectParticles.C b/src/enzo/CommunicationCollectParticles.C index 7d65dc120..ed59ad6e6 100644 --- a/src/enzo/CommunicationCollectParticles.C +++ b/src/enzo/CommunicationCollectParticles.C @@ -456,6 +456,14 @@ int CommunicationCollectParticles(LevelHierarchyEntry *LevelArray[], if (SyncNumberOfParticles) CommunicationSyncNumberOfParticles(GridHierarchyPointer, NumberOfGrids); + else { + for (i = StartGrid; i < EndGrid; i++) + if (MyProcessorNumber != + GridHierarchyPointer[i]->GridData->ReturnProcessorNumber()) { + GridHierarchyPointer[i]->GridData->SetNumberOfParticles(0); + GridHierarchyPointer[i]->GridData->SetNumberOfStars(0); + } + } #ifdef DEBUG_CCP for (i = StartGrid; i < EndGrid; i++)