commit/enzo-dev: 15 new changesets
1 view
Skip to first unread message
commits...@bitbucket.org
Feb 4, 2019, 4:11:00 PM2/4/19
to enzo...@googlegroups.com
15 new commits in enzo-dev:
https://bitbucket.org/enzo/enzo-dev/commits/67819df51893/
Changeset: 67819df51893
Branch: week-of-code
User: jwise77
Date: 2018-04-10 13:42:59+00:00
Summary: Adding an option to make split particles into must-refine ones
Affected #: 5 files
diff -r 149dd1d6513b -r 67819df51893 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -185,6 +185,13 @@
NumberOfNewParticles, MyProcessorNumber);
#endif
+ /* If specified, set particle type to must-refine */
+
+ if (ParticleSplitterMustRefine == TRUE) {
+ for (i = 0; i < NumberOfNewParticles; i++)
+ tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ }
+
/* If not set in the above routine, then set the metal fraction to zero. */
if (MetallicityField == FALSE && NumberOfParticleAttributes > 1)
diff -r 149dd1d6513b -r 67819df51893 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1262,6 +1262,8 @@
&ParticleSplitterRandomSeed);
ret += sscanf(line, "ParticleSplitterChildrenParticleSeparation = %"FSYM,
&ParticleSplitterChildrenParticleSeparation);
+ ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
+ &ParticleSplitterMustRefine);
ret += sscanf(line, "ResetMagneticField = %"ISYM,
&ResetMagneticField);
ret += sscanf(line, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM,
diff -r 149dd1d6513b -r 67819df51893 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -962,6 +962,7 @@
ParticleSplitterIterations = FALSE;
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
+ ParticleSplitterMustRefine = FALSE;
/* Magnetic Field Resetter */
diff -r 149dd1d6513b -r 67819df51893 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -624,6 +624,8 @@
ParticleSplitterRandomSeed);
fprintf(fptr, "ParticleSplitterChildrenParticleSeparation = %"FSYM"\n",
ParticleSplitterChildrenParticleSeparation);
+ fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
+ ParticleSplitterMustRefine);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
ResetMagneticField);
fprintf(fptr, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM"\n",
diff -r 149dd1d6513b -r 67819df51893 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -984,6 +984,7 @@
EXTERN int ParticleSplitterIterations;
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
+EXTERN int ParticleSplitterMustRefine;
/* Magnetic Field Resetter */
https://bitbucket.org/enzo/enzo-dev/commits/cfe6395292a7/
Changeset: cfe6395292a7
Branch: week-of-code
User: jwise77
Date: 2018-04-13 16:51:38+00:00
Summary: Adding additional particle splitting functionality from John Regan's
enzo-3.0 fork.
Affected #: 9 files
diff -r 67819df51893 -r cfe6395292a7 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2409,13 +2409,14 @@
/* Particle splitter routine. */
- int ParticleSplitter(int level);
+ int ParticleSplitter(int level, int iter);
int CreateChildParticles(float dx, int NumberOfParticles, float *ParticleMass,
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles);
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles);
/* Magnetic field resetting routine. */
diff -r 67819df51893 -r cfe6395292a7 src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -33,6 +33,8 @@
#include "fortran.def"
#include "CosmologyParameters.h"
+#define DEBUG_PS
+
void mt_init(unsigned_int seed);
unsigned_long_int mt_random();
@@ -41,13 +43,15 @@
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles)
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles)
{
int partnum = 0, i = 0, child = 0, m = 0, numpart = 0, innerchild = 0;
int xindex=0, yindex=0, zindex=0;
float rad = 0.0, sin60 = 0.0;
FLOAT NewPos[3][CHILDRENPERPARENT];
+ FLOAT LeftEdge[3], RightEdge[3];
int ChildrenPerParent = CHILDRENPERPARENT, total_children = 0;
float alpha[3];
float l11 = 0.0, l12 = 0.0, l13 = 0.0;
@@ -67,6 +71,63 @@
//Reproducible random seed
mt_init(((unsigned_int) ParticleSplitterRandomSeed));
+ FLOAT sep[3], midpoint[3], newsep[3];
+
+ /*
+ * The following three options determine over what region
+ * particle splitting takes place.
+ * By default the splitting occurs over the entire most refined
+ * region. At times this can be too large so two further options
+ * exist controlled via the parameter file. The first just selects
+ * a fraction of that region centred on the centre of the most
+ * refined region.
+ * The second allows the user to specify the centre of the splitting
+ * and the size of the region.
+ */
+
+ /* Default Option */
+ for (i = 0; i < 3; i++) {
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ /* Change particle refinement region in basic way */
+ if(ParticleSplitterFraction[iter] != 1.0)
+ {
+ fprintf(stderr, "Setting Particle Refinement using fractions.\n");
+ for (i = 0; i < 3; i++) {
+ sep[i] = RefineRegionRightEdge[i] - RefineRegionLeftEdge[i];
+ midpoint[i] = sep[i]/2.0 + RefineRegionLeftEdge[i];
+ newsep[i] = sep[i]*ParticleSplitterFraction[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ }
+ }
+ /* Centre particle refinement region around a specific point */
+ if(ParticleSplitterCenter[0] > 0.0 && ParticleSplitterCenterRegion[iter] > 0.0) {
+ fprintf(stderr, "Setting Particle Refinement around point.\n");
+ for (i = 0; i < 3; i++) {
+ midpoint[i] = ParticleSplitterCenter[i];
+ newsep[i] = ParticleSplitterCenterRegion[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ if(LeftEdge[i] < RefineRegionLeftEdge[i])
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ if(RightEdge[i] > RefineRegionRightEdge[i])
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ }
+
+ #ifdef DEBUG_PS
+ fprintf(stdout, "%s: Iteration %d: midpoint = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, midpoint[0], midpoint[1], midpoint[2]);
+ fprintf(stdout, "%s: Iteration %d: newsep = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, newsep[0], newsep[1], newsep[2]);
+ fprintf(stdout, "%s: Iteration %d: Refine Particles between (%lf, %lf, %lf) and (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, LeftEdge[0], LeftEdge[1], LeftEdge[2],
+ RightEdge[0], RightEdge[1], RightEdge[2]);
+ #endif
+
+
/*
* Loop over existing (parent) particles; It implicitly assumes that
* only DM and conventional star particles get splitted. Other
@@ -84,12 +145,12 @@
* blatantly discriminate against non-refined types around
* here.
*/
- if(ParticlePosition[0][partnum] < RefineRegionLeftEdge[0] ||
- ParticlePosition[0][partnum] > RefineRegionRightEdge[0] ||
- ParticlePosition[1][partnum] < RefineRegionLeftEdge[1] ||
- ParticlePosition[1][partnum] > RefineRegionRightEdge[1] ||
- ParticlePosition[2][partnum] < RefineRegionLeftEdge[2] ||
- ParticlePosition[2][partnum] > RefineRegionRightEdge[2])
+ if(ParticlePosition[0][partnum] < LeftEdge[0] ||
+ ParticlePosition[0][partnum] > RightEdge[0] ||
+ ParticlePosition[1][partnum] < LeftEdge[1] ||
+ ParticlePosition[1][partnum] > RightEdge[1] ||
+ ParticlePosition[2][partnum] < LeftEdge[2] ||
+ ParticlePosition[2][partnum] > RightEdge[2])
{
//fprintf(stdout, "grid::PS: Particle outside RR - ignoring\n");
continue;
@@ -302,8 +363,9 @@
}
}
#ifdef DEBUG_PS
- fprintf(stdout,"%d new child particles created in this grid from %d parents",
- total_children, NumberOfParticles);
+ if(total_children > 0)
+ fprintf(stdout,"Iteration %d: %d new child particles created in this grid from %d parents\n",
+ iter, total_children, NumberOfParticles);
#endif
*NumberOfNewParticles = total_children;
diff -r 67819df51893 -r cfe6395292a7 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -40,22 +40,7 @@
float *VelocityUnits, FLOAT Time);
int FindField(int field, int farray[], int numfields);
-extern "C" void FORTRAN_NAME(particle_splitter)(int *nx, int *ny, int *nz,
- int *idual, int *imetal, hydro_method *imethod, float *dt,
- float *r, float *dx, FLOAT *t, float *z,
- float *d1, float *x1, float *v1, float *t1,
- FLOAT *xstart, FLOAT *ystart, FLOAT *zstart, int *ibuff,
- int *npart,
- FLOAT *xpold, FLOAT *ypold, FLOAT *zpold, float *upold, float *vpold, float *wpold,
- float *mpold, float *tdpold, float *tcpold, float *metalfold, int *typeold,
- int *nmax, int *npartnew, int *children, int *level,
- FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
- float *mp, float *tdp, float *tcp, float *metalf, int *type,
- int *iterations, float *separation, int *ran1_init,
- FLOAT *rr_leftedge, FLOAT *rr_rightedge);
-
-
-int grid::ParticleSplitter(int level)
+int grid::ParticleSplitter(int level, int iteration)
{
if (ParticleSplitterIterations == 0)
@@ -131,7 +116,7 @@
/* Allocate space for new particles. */
- int ChildrenPerParent = 12; //12+1 = 13 will be the final number of particles per parent
+ int ChildrenPerParent = CHILDRENPERPARENT; //12+1 = 13 will be the final number of particles per parent
int MaximumNumberOfNewParticles = ChildrenPerParent*NumberOfParticles+1;
tg->AllocateNewParticles(MaximumNumberOfNewParticles);
@@ -168,7 +153,8 @@
if(!tg->CreateChildParticles(CellWidthTemp, NumberOfParticles, ParticleMass,
ParticleType, ParticlePosition, ParticleVelocity,
ParticleAttribute, CellLeftEdge, GridDimension,
- MaximumNumberOfNewParticles, &NumberOfNewParticles))
+ MaximumNumberOfNewParticles, iteration,
+ &NumberOfNewParticles))
{
fprintf(stdout, "Failed to create child particles in grid %d\n", this->GetGridID());
return FAIL;
diff -r 67819df51893 -r cfe6395292a7 src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -38,7 +38,7 @@
#include "CosmologyParameters.h"
#include "CommunicationUtilities.h"
-#define NO_DEBUG_PS
+#define NO_DEBUG_PS
int RebuildHierarchy(TopGridData *MetaData,
LevelHierarchyEntry *LevelArray[], int level);
@@ -70,6 +70,13 @@
ParticleSplitterChildrenParticleSeparation <=0)
return SUCCESS;
+ if(ParticleSplitterIterations > MAX_SPLIT_ITERATIONS) {
+ fprintf(stderr, "WARNING: Splitting iterations exceeds maximum allowed\n" \
+ "Resetting to %d\n", MAX_SPLIT_ITERATIONS);
+ ParticleSplitterIterations = MAX_SPLIT_ITERATIONS;
+ }
+
+
/* Find total NumberOfParticles in all grids; this is needed in
CommunicationUpdateStarParticleCount below */
@@ -78,9 +85,6 @@
/* Initialize all star particles if this is a restart */
- if (ParticleSplitterIterations > 1)
- fprintf(stderr, "WARNING: ParticleSplitterIterations > 1 is not properly tested yet.\n");
-
for (i = 0; i < ParticleSplitterIterations; i++) {
for (level = 0; level < MAX_DEPTH_OF_HIERARCHY-1; level++) {
@@ -98,7 +102,6 @@
RecordTotalStarParticleCount(Grids, NumberOfGrids,
TotalStarParticleCountPrevious);
-
// for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
// fprintf(stdout, "TotalStarParticleCountPrevious[grid=%d] = %d\n", grid1,
// TotalStarParticleCountPrevious[grid1]);
@@ -110,7 +113,7 @@
Grids[grid1]->GridData->ReturnNumberOfParticles());
#endif
- if (Grids[grid1]->GridData->ParticleSplitter(level) == FAIL) {
+ if (Grids[grid1]->GridData->ParticleSplitter(level, i) == FAIL) {
ENZO_FAIL("Error in grid::ParticleSplitter.\n");
}
diff -r 67819df51893 -r cfe6395292a7 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1264,6 +1264,14 @@
&ParticleSplitterChildrenParticleSeparation);
ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
&ParticleSplitterMustRefine);
+ ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
+ ParticleSplitterFraction+3);
+ ret += sscanf(line, "ParticleSplitterCenter = %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterCenter+0, ParticleSplitterCenter+1, ParticleSplitterCenter+2);
+ ret += sscanf(line, "ParticleSplitterCenterRegion = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterCenterRegion+0, ParticleSplitterCenterRegion+1,
+ ParticleSplitterCenterRegion+2, ParticleSplitterCenterRegion+3);
ret += sscanf(line, "ResetMagneticField = %"ISYM,
&ResetMagneticField);
ret += sscanf(line, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM,
diff -r 67819df51893 -r cfe6395292a7 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -963,6 +963,12 @@
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
ParticleSplitterMustRefine = FALSE;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterFraction[i] = 1.0;
+ for(int i = 0; i < MAX_DIMENSION; i++)
+ ParticleSplitterCenter[i] = -1.0;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterCenterRegion[i] = -1.0;
/* Magnetic Field Resetter */
diff -r 67819df51893 -r cfe6395292a7 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -624,6 +624,12 @@
ParticleSplitterRandomSeed);
fprintf(fptr, "ParticleSplitterChildrenParticleSeparation = %"FSYM"\n",
ParticleSplitterChildrenParticleSeparation);
+ fprintf(fptr, "ParticleSplitterFraction = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterFraction);
+ fprintf(fptr, "ParticleSplitterCenter = ");
+ WriteListOfFloats(fptr, 3, ParticleSplitterCenter);
+ fprintf(fptr, "ParticleSplitterCenterRegion = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterCenterRegion);
fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
ParticleSplitterMustRefine);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
diff -r 67819df51893 -r cfe6395292a7 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -985,6 +985,9 @@
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
EXTERN int ParticleSplitterMustRefine;
+EXTERN float ParticleSplitterFraction[MAX_SPLIT_ITERATIONS];
+EXTERN FLOAT ParticleSplitterCenter[MAX_DIMENSION];
+EXTERN float ParticleSplitterCenterRegion[MAX_SPLIT_ITERATIONS];
/* Magnetic Field Resetter */
diff -r 67819df51893 -r cfe6395292a7 src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -478,6 +478,11 @@
#define PARTICLE_TYPE_SIMPLE_SOURCE 10
#define CHILDRENPERPARENT 12
+
+/* You probably don't want to do this too often */
+
+#define MAX_SPLIT_ITERATIONS 4
+
/* Ways to deposit particles from a subgrid. */
#define CIC_DEPOSIT 0
https://bitbucket.org/enzo/enzo-dev/commits/fbcce23071d3/
Changeset: fbcce23071d3
Branch: week-of-code
User: jwise77
Date: 2018-04-16 11:15:16+00:00
Summary: Correcting float macro for particle splitter region center
Affected #: 1 file
diff -r cfe6395292a7 -r fbcce23071d3 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1267,7 +1267,7 @@
ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
ParticleSplitterFraction+3);
- ret += sscanf(line, "ParticleSplitterCenter = %"FSYM" %"FSYM" %"FSYM"",
+ ret += sscanf(line, "ParticleSplitterCenter = %"PSYM" %"PSYM" %"PSYM"",
ParticleSplitterCenter+0, ParticleSplitterCenter+1, ParticleSplitterCenter+2);
ret += sscanf(line, "ParticleSplitterCenterRegion = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
ParticleSplitterCenterRegion+0, ParticleSplitterCenterRegion+1,
https://bitbucket.org/enzo/enzo-dev/commits/ff55ab979b6b/
Changeset: ff55ab979b6b
Branch: week-of-code
User: jwise77
Date: 2018-04-18 15:22:30+00:00
Summary: Additional MRP check to not include DM particles that were created
from old star particles and merging routines.
Affected #: 1 file
diff -r fbcce23071d3 -r ff55ab979b6b src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -92,13 +92,23 @@
// Flag particles as must refine particles
int *IsParticleMustRefine;
IsParticleMustRefine = new int[NumberOfParticles];
+ bool OriginalParticle;
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
+ if (NumberOfParticleAttributes > 0)
+ OriginalParticle = (ParticleAttribute[0][i] < 0);
+ else
+ OriginalParticle = true;
- // check particle type and uniform mass
+ // check particle type and uniform mass. Also check particle
+ // creation time for DM particles that are positive, indicating
+ // that they are either inert stellar remnants or "leftovers" from
+ // particle merging
+
rules[0] = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
ParticleType[i] == PARTICLE_TYPE_MBH) ||
- (ParticleMass[i] < UniformParticleMass);
+ (ParticleMass[i] < UniformParticleMass &&
+ ParticleType[i] == PARTICLE_TYPE_DARK_MATTER && OriginalParticle);
// check particle mass greater than minimum mass
rules[1] = (ParticleMass[i] > MustRefineParticlesMinimumMass);
@@ -110,6 +120,7 @@
// set flag for this particle
for (int j = 0; j < NumberOfRules; j++)
IsParticleMustRefine[i] *= rules[j];
+
}
PFORTRAN_NAME(cic_flag)(IsParticleMustRefine,
https://bitbucket.org/enzo/enzo-dev/commits/2d5a79c5d3a5/
Changeset: 2d5a79c5d3a5
Branch: week-of-code
User: jwise77
Date: 2018-04-19 18:31:32+00:00
Summary: Adding capability to read particle IDs to mark as must-refine during
particle splitting.
Affected #: 8 files
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2409,7 +2409,8 @@
/* Particle splitter routine. */
- int ParticleSplitter(int level, int iter);
+ int ParticleSplitter(int level, int iter, int NumberOfIDs,
+ long *MustRefineIDs);
int CreateChildParticles(float dx, int NumberOfParticles, float *ParticleMass,
int *ParticleType, FLOAT *ParticlePosition[],
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -134,10 +134,14 @@
* particles - which usually become Star class particles - doesn't
* seem to have any reason to be splitted. (as of Oct.2009)
*/
+
+ /* (Apr 2018) Include must-refine particles */
for(partnum = 0; partnum < NumberOfParticles; partnum++)
{
- if(ParticleMass[partnum] > 0.0 && ParticleType[partnum] <= 2)
+ if(ParticleMass[partnum] > 0.0 &&
+ (ParticleType[partnum] <= 2 ||
+ ParticleType[partnum] == PARTICLE_TYPE_MUST_REFINE))
{
/*
* Check that particle is within the most refined region.
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -40,7 +40,8 @@
float *VelocityUnits, FLOAT Time);
int FindField(int field, int farray[], int numfields);
-int grid::ParticleSplitter(int level, int iteration)
+int grid::ParticleSplitter(int level, int iteration, int NumberOfIDs,
+ long *MustRefineIDs)
{
if (ParticleSplitterIterations == 0)
@@ -130,7 +131,7 @@
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles before splitting = %d, MyProcessorNumber = %d\n",
NumberOfParticles, MyProcessorNumber);
#endif
-
+
if (NumberOfParticles > 0) {
#define NO_PARTICLE_IN_GRID_CHECK
@@ -150,6 +151,24 @@
xindex, yindex, zindex, level);
}
#endif
+
+ /* Before splitting, set particle type to must-refine if given in
+ the provided list. Only do this for the first iteration because
+ in the later iterations, the particle types are propagated to
+ the split particles. */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs != NULL &&
+ iteration == 0) {
+ for (i = 0; i < NumberOfParticles; i++) {
+ for (j = 0; j < NumberOfIDs; j++) {
+ if (this->ParticleNumber[i] == MustRefineIDs[j]) {
+ this->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ break;
+ }
+ } // ENDFOR j
+ } // ENDFOR i
+ }
+
if(!tg->CreateChildParticles(CellWidthTemp, NumberOfParticles, ParticleMass,
ParticleType, ParticlePosition, ParticleVelocity,
ParticleAttribute, CellLeftEdge, GridDimension,
@@ -171,9 +190,10 @@
NumberOfNewParticles, MyProcessorNumber);
#endif
- /* If specified, set particle type to must-refine */
-
- if (ParticleSplitterMustRefine == TRUE) {
+ /* If specified and no ID list given, set all particle types to
+ must-refine */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs == NULL) {
for (i = 0; i < NumberOfNewParticles; i++)
tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
}
@@ -204,7 +224,7 @@
#ifdef DEBUG_PS
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles(New) = %d, MyProcessorNumber = %d\n",
- NumberOfParticles, MyProcessorNumber);
+ NumberOfParticles, MyProcessorNumber);
#endif
} // end: if (NumberOfNewParticles > 0)
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -24,6 +24,7 @@
#include <stdlib.h>
#include <stdio.h>
+#include <hdf5.h>
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "typedefs.h"
@@ -83,6 +84,34 @@
MetaData->NumberOfParticles = FindTotalNumberOfParticles(LevelArray);
NumberOfOtherParticles = MetaData->NumberOfParticles - NumberOfStarParticles;
+ /* Read must-refine particle positions if specified */
+
+ long *MustRefineIDs = NULL;
+ int NumberOfIDs = 0;
+
+ if (ParticleSplitterMustRefine &&
+ ParticleSplitterMustRefineIDFile != NULL) {
+
+ int NumberOfMustRefineIDs;
+ hid_t file_id, dataset_id, dataspace_id;
+ herr_t status;
+ hsize_t dims[1], maxdims[1];
+
+ file_id = H5Fopen(ParticleSplitterMustRefineIDFile,
+ H5F_ACC_RDONLY, H5P_DEFAULT);
+ dataset_id = H5Dopen2(file_id, "/particle_identifier", H5P_DEFAULT);
+ dataspace_id = H5Dget_space(dataset_id);
+ // Get number of particles and allocate variable
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, maxdims);
+ MustRefineIDs = new long[dims[0]];
+ NumberOfIDs = dims[0];
+ status = H5Dread(dataset_id, H5T_NATIVE_LLONG, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ MustRefineIDs);
+ status = H5Dclose(dataset_id);
+ status = H5Fclose(file_id);
+
+ }
+
/* Initialize all star particles if this is a restart */
for (i = 0; i < ParticleSplitterIterations; i++) {
@@ -113,7 +142,8 @@
Grids[grid1]->GridData->ReturnNumberOfParticles());
#endif
- if (Grids[grid1]->GridData->ParticleSplitter(level, i) == FAIL) {
+ if (Grids[grid1]->GridData->ParticleSplitter(level, i, NumberOfIDs,
+ MustRefineIDs) == FAIL) {
ENZO_FAIL("Error in grid::ParticleSplitter.\n");
}
@@ -151,6 +181,7 @@
/* Set splitter parameter zero; otherwise it will split particles again at next restart */
+ delete [] MustRefineIDs;
ParticleSplitterIterations = 0;
return SUCCESS;
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1264,6 +1264,8 @@
&ParticleSplitterChildrenParticleSeparation);
ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
&ParticleSplitterMustRefine);
+ if (sscanf(line, "ParticleSplitterMustRefineIDFile = %s", dummy) == 1)
+ ParticleSplitterMustRefineIDFile = dummy;
ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
ParticleSplitterFraction+3);
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -963,6 +963,7 @@
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
ParticleSplitterMustRefine = FALSE;
+ ParticleSplitterMustRefineIDFile = NULL;
for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
ParticleSplitterFraction[i] = 1.0;
for(int i = 0; i < MAX_DIMENSION; i++)
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -632,6 +632,8 @@
WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterCenterRegion);
fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
ParticleSplitterMustRefine);
+ fprintf(fptr, "ParticleSplitterMustRefineIDFile = %s\n",
+ ParticleSplitterMustRefineIDFile);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
ResetMagneticField);
fprintf(fptr, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM"\n",
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -985,6 +985,7 @@
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
EXTERN int ParticleSplitterMustRefine;
+EXTERN char *ParticleSplitterMustRefineIDFile;
EXTERN float ParticleSplitterFraction[MAX_SPLIT_ITERATIONS];
EXTERN FLOAT ParticleSplitterCenter[MAX_DIMENSION];
EXTERN float ParticleSplitterCenterRegion[MAX_SPLIT_ITERATIONS];
https://bitbucket.org/enzo/enzo-dev/commits/5409f66b0ab2/
Changeset: 5409f66b0ab2
Branch: week-of-code
User: jwise77
Date: 2018-04-20 14:24:19+00:00
Summary: Don't rebuild hierarchy before particle splitting. Will lose information if using must-refine particles.
Affected #: 2 files
diff -r 2d5a79c5d3a5 -r 5409f66b0ab2 src/enzo/EvolveHierarchy.C
--- a/src/enzo/EvolveHierarchy.C
+++ b/src/enzo/EvolveHierarchy.C
@@ -270,20 +270,26 @@
*/
+ /* Particle Splitter. Split particles into 13 (=1+12) child
+ particles. The hierarchy is rebuilt inside this routine. */
+
+ if (MetaData.FirstTimestepAfterRestart == TRUE &&
+ ParticleSplitterIterations > 0) {
+
+ ParticleSplitter(LevelArray, 0, &MetaData);
+
+ } else {
+
/* Do the first grid regeneration. */
- if(CheckpointRestart == FALSE) {
- RebuildHierarchy(&MetaData, LevelArray, 0);
- }
+ if(CheckpointRestart == FALSE) {
+ RebuildHierarchy(&MetaData, LevelArray, 0);
+ }
+
+ } // ENDELSE particle splitting
PrintMemoryUsage("1st rebuild");
- /* Particle Splitter. Split particles into 13 (=1+12) child particles */
-
- if (MetaData.FirstTimestepAfterRestart == TRUE &&
- ParticleSplitterIterations > 0)
- ParticleSplitter(LevelArray, 0, &MetaData);
-
/* Reset magnetic fields if requested. */
if (MetaData.FirstTimestepAfterRestart == TRUE &&
diff -r 2d5a79c5d3a5 -r 5409f66b0ab2 src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -61,15 +61,13 @@
HierarchyEntry **Grids;
int NumberOfGrids;
- /* First, rebuild the hierarchy */
-
- RebuildHierarchy(MetaData, LevelArray, 0);
-
- /* Return if this does not concern us */
+ /* Rebuild hierarchy and return if this does not concern us */
if (ParticleSplitterIterations <= 0 ||
- ParticleSplitterChildrenParticleSeparation <=0)
+ ParticleSplitterChildrenParticleSeparation <=0) {
+ RebuildHierarchy(MetaData, LevelArray, 0);
return SUCCESS;
+ }
if(ParticleSplitterIterations > MAX_SPLIT_ITERATIONS) {
fprintf(stderr, "WARNING: Splitting iterations exceeds maximum allowed\n" \
https://bitbucket.org/enzo/enzo-dev/commits/828e88f05f0f/
Changeset: 828e88f05f0f
Branch: week-of-code
User: jwise77
Date: 2018-04-23 14:18:14+00:00
Summary: Additional restrictions on split must-refine particles
Affected #: 2 files
diff -r 5409f66b0ab2 -r 828e88f05f0f src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -93,10 +93,11 @@
int *IsParticleMustRefine;
IsParticleMustRefine = new int[NumberOfParticles];
bool OriginalParticle;
+ const float OriginalCrit = FLOAT_UNDEFINED+1;
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
if (NumberOfParticleAttributes > 0)
- OriginalParticle = (ParticleAttribute[0][i] < 0);
+ OriginalParticle = (ParticleAttribute[0][i] < OriginalCrit);
else
OriginalParticle = true;
diff -r 5409f66b0ab2 -r 828e88f05f0f src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -197,6 +197,12 @@
for (i = 0; i < NumberOfNewParticles; i++)
tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
}
+
+ if (NumberOfParticleAttributes > 0)
+ for (i = 0; i < NumberOfNewParticles; i++) {
+ if (tg->ParticleType[i] == PARTICLE_TYPE_DARK_MATTER)
+ tg->ParticleAttribute[0][i] = -1.0;
+ }
/* If not set in the above routine, then set the metal fraction to zero. */
https://bitbucket.org/enzo/enzo-dev/commits/68d91a7109f9/
Changeset: 68d91a7109f9
Branch: week-of-code
User: jwise77
Date: 2018-05-04 15:44:42+00:00
Summary: More debugging of split must-refine particles
Affected #: 3 files
diff -r 828e88f05f0f -r 68d91a7109f9 src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -47,7 +47,7 @@
int *NumberOfNewParticles)
{
- int partnum = 0, i = 0, child = 0, m = 0, numpart = 0, innerchild = 0;
+ int partnum = 0, i = 0, child = 0, m = 0, innerchild = 0;
int xindex=0, yindex=0, zindex=0;
float rad = 0.0, sin60 = 0.0;
FLOAT NewPos[3][CHILDRENPERPARENT];
@@ -350,8 +350,13 @@
for(i = 0; i < 3; i++)
this->ParticleVelocity[i][child] = ParticleVelocity[i][partnum];
this->ParticleType[child] = ParticleType[partnum];
+ // Flag that a DM particle was split (originally -99999 or 0)
+ if (NumberOfParticleAttributes > 0 &&
+ ParticleType[partnum] == PARTICLE_TYPE_DARK_MATTER)
+ if (ParticleAttribute[0][partnum] <= 0)
+ ParticleAttribute[0][partnum] = tiny_number;
for(i = 0; i < NumberOfParticleAttributes; i++)
- this->ParticleAttribute[i][child] = ParticleAttribute[i][numpart];
+ this->ParticleAttribute[i][child] = ParticleAttribute[i][partnum];
}
/* Loop forward by CHILDRENPERPARENT each time. */
diff -r 828e88f05f0f -r 68d91a7109f9 src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -93,11 +93,10 @@
int *IsParticleMustRefine;
IsParticleMustRefine = new int[NumberOfParticles];
bool OriginalParticle;
- const float OriginalCrit = FLOAT_UNDEFINED+1;
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
if (NumberOfParticleAttributes > 0)
- OriginalParticle = (ParticleAttribute[0][i] < OriginalCrit);
+ OriginalParticle = (ParticleAttribute[0][i] <= 0.0);
else
OriginalParticle = true;
@@ -115,7 +114,6 @@
rules[1] = (ParticleMass[i] > MustRefineParticlesMinimumMass);
// add more rules here
-
//
// set flag for this particle
diff -r 828e88f05f0f -r 68d91a7109f9 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -198,12 +198,6 @@
tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
}
- if (NumberOfParticleAttributes > 0)
- for (i = 0; i < NumberOfNewParticles; i++) {
- if (tg->ParticleType[i] == PARTICLE_TYPE_DARK_MATTER)
- tg->ParticleAttribute[0][i] = -1.0;
- }
-
/* If not set in the above routine, then set the metal fraction to zero. */
if (MetallicityField == FALSE && NumberOfParticleAttributes > 1)
https://bitbucket.org/enzo/enzo-dev/commits/3e4499464d7b/
Changeset: 3e4499464d7b
Branch: week-of-code
User: jwise77
Date: 2018-10-09 15:43:08+00:00
Summary: Adding must-refine particle splitting parameters to docs
Affected #: 1 file
diff -r 68d91a7109f9 -r 3e4499464d7b doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1865,10 +1865,44 @@
Currently it implicitly assumes that only DM (type=1) and
conventional star particles (type=2) inside the ``RefineRegion`` get
split. Other particles, which usually become Star class objects,
- seem to have no reason to be split. Default: 0
+ seem to have no reason to be split. It is reset to zero after a
+ restart to avoid resplitting in subsequent datasets. It can be set
+ to a maximum of 4. Default: 0
+``ParticleSplitterRandomSeed`` (external)
+ Random seed used when randomly rotating the hexagonal close
+ packed array on whose vertices the split particles are
+ placed. Default: 131180
+``ParticleSplitterMustRefine`` (external)
+ Set to 1 to mark the split particles as must-refine particles. The
+ user must also set associated must-refine parameters to enable its
+ machinery. Default: 0
+``ParticleSplitterMustRefineIDFile`` (external)
+ Filename for the HDF5 file that has a dataset containing the
+ particle IDs that should be marked as must-refine. All other
+ particles within the region marked for splitting will retain their
+ original types. This must be used in conjunction with
+ ``ParticleSplitterMustRefine = 1``. The dataset must be named
+ ``particle_identifier`` in the base group. Default: (null)
+``ParticleSplitterFraction`` (external)
+ An array of four values that represent the width of the splitting
+ region in units of the original refine region set by
+ ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. Each
+ successive value represents the next nested split region. Valid
+ up to ``ParticleSplitterIterations`` times. Default: 1.0 (all 4
+ values)
+``ParticleSplitterCenter`` (external)
+ The center of split region in code units. Specify if the split
+ region does not correspond to the center of the refine region.
+ Not used if negative. Default: -1.0 -1.0 -1.0
+``ParticleSplitterCenterRegion`` (external)
+ The width of the split region in code units. Must be used in
+ conjunction with ``ParticleSplitterCenter``. Each successive
+ value represents the next nested split region. Valid up to
+ ``ParticleSplitterIterations`` times. Not used if
+ negative. Default: -1.0 (all 4 values)
``ParticleSplitterChildrenParticleSeparation`` (external)
This is the spacing between the child particles placed on a
- hexagonal close-packed (HCP) array. In the unit of a cell size
+ hexagonal close-packed (HCP) array. In units of a cell size
which the parent particle resides in. Default: 1.0
.. _starparticleparameters:
https://bitbucket.org/enzo/enzo-dev/commits/40e4a6cb3d8a/
Changeset: 40e4a6cb3d8a
Branch: week-of-code
User: jwise77
Date: 2018-10-16 15:49:24+00:00
Summary: Adding an example yt script to generate particle splitting ID HDF5 file
Affected #: 1 file
diff -r 3e4499464d7b -r 40e4a6cb3d8a doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1882,7 +1882,22 @@
particles within the region marked for splitting will retain their
original types. This must be used in conjunction with
``ParticleSplitterMustRefine = 1``. The dataset must be named
- ``particle_identifier`` in the base group. Default: (null)
+ ``particle_identifier`` in the base group. An example yt script is
+ provided below, selecting the particles in a sphere centered at
+ [0.5, 0.5, 0.5] with a radius 0.05 in code length units, using yt::
+
+ import yt
+ import h5py as h5
+
+ ds = yt.load('DD0040/DD0040')
+ center = ds.arr([0.5, 0.5, 0.5], 'code_length')
+ radius = ds.quan(0.05, 'code_length')
+ sp = ds.sphere(center, radius)
+ fp = h5.File('particle-ids.h5', 'w')
+ fp['particle_identifier'] = sp['particle_index'].astype('int')
+ fp.close()
+
+ Default: (null)
``ParticleSplitterFraction`` (external)
An array of four values that represent the width of the splitting
region in units of the original refine region set by
https://bitbucket.org/enzo/enzo-dev/commits/93c2de54c3fc/
Changeset: 93c2de54c3fc
Branch: week-of-code
User: jwise77
Date: 2018-10-16 16:29:06+00:00
Summary: Correct syntax highlighting for sample code
Affected #: 1 file
diff -r 40e4a6cb3d8a -r 93c2de54c3fc doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1882,9 +1882,13 @@
particles within the region marked for splitting will retain their
original types. This must be used in conjunction with
``ParticleSplitterMustRefine = 1``. The dataset must be named
- ``particle_identifier`` in the base group. An example yt script is
- provided below, selecting the particles in a sphere centered at
- [0.5, 0.5, 0.5] with a radius 0.05 in code length units, using yt::
+ ``particle_identifier`` in the base group. Default: (null).
+
+ An example yt script is provided below, selecting the particles in
+ a sphere centered at [0.5, 0.5, 0.5] with a radius 0.05 in code
+ length units.
+
+ .. code-block:: python
import yt
import h5py as h5
@@ -1897,7 +1901,6 @@
fp['particle_identifier'] = sp['particle_index'].astype('int')
fp.close()
- Default: (null)
``ParticleSplitterFraction`` (external)
An array of four values that represent the width of the splitting
region in units of the original refine region set by
https://bitbucket.org/enzo/enzo-dev/commits/ff98226315b7/
Changeset: ff98226315b7
Branch: week-of-code
User: jwise77
Date: 2019-01-09 05:01:55+00:00
Summary: Merging with latest mainline updates
Affected #: 39 files
diff -r 93c2de54c3fc -r ff98226315b7 doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -904,7 +904,21 @@
constrained by static grids. Instead, refinement around dark matter
particles is allowed down to the level of a particle's generation level.
Refinement beyond this level is allowed around particles within the MUSIC
- ellipsoidal making region.
+ ellipsoidal masking region. Note, dark matter particles corresponding to
+ a generation level N are guaranteed to be refined to at least level N, but
+ may also exist on levels N > 1 if in the vicinity of an N > 1 dark matter
+ particle or a tagged must-refine particle.
+ 4. Similar to setting 3, except dark matter particles corresponding to a
+ generation level N are refined only to level N and no further. If two
+ dark matter particles from different levels occupy the same cell, that
+ cell will be refined to the coarser level. Tagged must-refine particles
+ near coarse dark matter particles will be similarly de-refined. Compared
+ to option 3, this can be used to prevent unnecessary mesh refinement in
+ regions where coarser particles enter into higher resolution regions,
+ slowing down the simulation. Note, with this setting, coarse boundary
+ particles entering into a high resolution region will eventually lead
+ to total derefinement of the region of interest.
+
``MustRefineParticlesRefineToLevel`` (external)
The maximum level on which ``MustRefineParticles`` are required to
diff -r 93c2de54c3fc -r ff98226315b7 doc/manual/source/physics/star_particles.rst
--- a/doc/manual/source/physics/star_particles.rst
+++ b/doc/manual/source/physics/star_particles.rst
@@ -552,29 +552,36 @@
----------------------------
*Source: hydro_rk/SuperNovaSeedField.C*
-Select this method by setting ``UseSupernovaSeedFieldSourceTerms = 1``
-(Default = 0) and
-specifying the following parameters:
+Select this method by setting ``UseMagneticSupernovaFeedback = 1``
+(Default = 0) and specifying the parameters below. If
+``UseMagneticSupernovaFeedback == 2``, the parameters
+``MagneticSupernovaRadius`` and ``MagneticSupernovaDuration`` will
+be calculated to be the minimum allowed values (see below)
+at runtime based on the cell width and time step of the most-refined grid.
-``SupernovaSeedFieldTotalEnergy`` (in units of ergs) is the total amount
-of magnetic energy to be injected by a single supernova event. Defualt = 0.0.
-``SupernovaSeedFieldRadius`` (in units of parsecs) gives the scale over
+``MagneticSupernovaEnergy`` (in units of ergs) is the total amount
+of magnetic energy to be injected by a single supernova event. Defualt = 1e51
+ergs.
+
+``MagneticSupernovaRadius`` (in units of parsecs) gives the scale over
which to inject supernova energy. The injection mechanism normalizes the
spatial exponential decay of the injected supernova energy so that all of the
energy is contained within the specified radius. For this reason, the
-``SupernovaSeedFieldRadius`` should be at least 3 times the minimum cell width of
-the simulation. Default = 0.0.
+``MagneticSupernovaRadius`` should be at least 1.5 times the minimum cell width of
+the simulation (in pc). Default = 300 pc.
-``SupernovaSeedFieldDuration`` (in units of Myr) gives the duration of the
+``MagneticSupernovaDuration`` (in units of years) gives the duration of the
supernova magnetic energy injection. The injection mechanism is normalized so
-that all of the ``SupernovaSeedFieldTotalEnergy`` is injected over this
-time scale. In order to inject the correct amount of energy, ``SupernovaSeedFieldDuration`` should be set to at least 4
-times the minimum time step of the simulation. Default = 0.0.
+that all of the ``MagneticSupernovaEnergy`` is injected over this
+time scale. In order to inject the correct amount of energy, ``MagneticSupernovaDuration`` should be set to at least 5
+times the minimum time step of the simulation. Default = 50,000 years.
-The following applies to Methods 0 (Cen & Ostriker) and 1 (+
-stochastic star formation). The magnetic feedback method is described fully in `Butsky et al. (2017)
+The following applies to all star formation methods that produce a
+PARTICLE_TYPE_STAR object. Methods 0 (Cen & Ostriker) and 1 (+
+stochastic star formation) have been tested extensively.
+The magnetic feedback method is described fully in `Butsky et al. (2017)
<http://adsabs.harvard.edu/abs/2017ApJ...843..113B>`_.
When a star cluster particle reaches the end of its lifetime, we inject a
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/CosmologyReadParameters.C
--- a/src/enzo/CosmologyReadParameters.C
+++ b/src/enzo/CosmologyReadParameters.C
@@ -45,6 +45,7 @@
FinalRedshift = 0;
CosmologyTableNumberOfBins = 1000;
CosmologyTableLogt = NULL;
+ CosmologyTableLoga = NULL;
CosmologyTableLogaInitial = -6.0;
CosmologyTableLogaFinal = 0.0;
CosmologyTableLogtIndex = 0;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/EnzoTiming.h
--- a/src/enzo/EnzoTiming.h
+++ b/src/enzo/EnzoTiming.h
@@ -29,8 +29,15 @@
#include <cstring>
#include <map>
-#define min(A,B) ((A) < (B) ? (A) : (B))
-#define max(A,B) ((A) > (B) ? (A) : (B))
+template <typename T>
+bool min(const T& A, const T& B) {
+ return A < B ? A : B;
+}
+
+template <typename T>
+bool max(const T& A, const T& B) {
+ return A > B ? A : B;
+}
double ReturnWallTime(void);
void Reduce_Times(double time, double *time_array);
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/EvolveLevel.C
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -540,6 +540,10 @@
Grids[grid1]->GridData->CopyBaryonFieldToOldBaryonField();
+ // Find recently-supernova stars to add them the MagneticSupernovaList
+ if ((UseMagneticSupernovaFeedback) && (level == MaximumRefinementLevel))
+ Grids[grid1]->GridData->AddMagneticSupernovaeToList();
+
/* Call hydro solver and save fluxes around subgrids.
* HD_RK and MHD_RK are the 2nd order Runge-Kutta integrations, which
* require two steps (*_1stStep and *_2ndStep)
@@ -547,7 +551,7 @@
* All others (PPM, Zeus, MHD_Li/CT) are called from SolveHydroEquations
*/
-
+
if( HydroMethod != HD_RK && HydroMethod != MHD_RK ){
Grids[grid1]->GridData->SolveHydroEquations(LevelCycleCount[level],
NumberOfSubgrids[grid1], SubgridFluxesEstimate[grid1], level);
@@ -713,7 +717,10 @@
if (ComovingCoordinates)
Grids[grid1]->GridData->ComovingExpansionTerms();
- } // end loop over grids
+ if (UseMagneticSupernovaFeedback)
+ Grids[grid1]->GridData->MagneticSupernovaList.clear();
+
+ } // end loop over grids
/* Finalize (accretion, feedback, etc.) star particles */
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FinalizeFluxes.C
--- a/src/enzo/FinalizeFluxes.C
+++ b/src/enzo/FinalizeFluxes.C
@@ -48,10 +48,11 @@
DeleteFluxes(SubgridFluxesEstimate[grid1][subgrid]);
delete SubgridFluxesEstimate[grid1][subgrid];
}
- delete [] SubgridFluxesEstimate[grid1];
-
}
-
+
+ delete [] SubgridFluxesEstimate[grid1];
+ SubgridFluxesEstimate[grid1] = NULL;
+
} // end of loop over grids
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FindSubgrids.C
--- a/src/enzo/FindSubgrids.C
+++ b/src/enzo/FindSubgrids.C
@@ -85,7 +85,8 @@
}
/* Set the static (permanent) regions. */
- if (MustRefineParticlesCreateParticles != 3){
+ if (MustRefineParticlesCreateParticles != 3 &&
+ MustRefineParticlesCreateParticles != 4) {
if (CurrentGrid->SetFlaggingFieldStaticRegions(level, NumberOfFlaggedCells)
== FAIL) {
ENZO_FAIL("Error in grid->SetFlaggingFieldStaticRegions.");
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FreeRealMem.C
--- /dev/null
+++ b/src/enzo/FreeRealMem.C
@@ -0,0 +1,199 @@
+#ifdef TASKMAP
+
+#ifdef USE_LL
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <sys/types.h>
+#include "llapi.h"
+
+#include "macros_and_parameters.h"
+
+
+Eint64 FreeRealMem( char *node )
+{
+
+ LL_element *queryObject;
+ LL_element *job, *step;
+
+ Eint32 rc, obj_count, err_code, ok;
+
+ Eint64 free;
+ Eint64 freemem;
+ Eint64 def_mem;
+
+ def_mem = (Eint64)(12.0); /* minimum memory on DataStar */
+
+ Eint32 node_number;
+
+#ifdef SEABORG
+ sscanf(node, "s%5d", &node_number);
+ fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
+
+ if ( 701 <= node_number && node_number <= 713 ) {
+ def_mem = (Eint64)(56.0);
+ } else {
+
+ if ( 3701 <= node_number && node_number <= 5213 ) {
+ def_mem = (Eint64)(28.0);
+ } else {
+ def_mem = (Eint64)(14.0);
+ }
+
+ }
+#endif
+
+#ifdef DATASTAR
+ sscanf(node, "ds%3d", &node_number);
+ fprintf(stderr, "DataStar Node %s number %d\n", node, node_number);
+
+ if ( 4 <= node_number && node_number <= 11 ) {
+ def_mem = (Eint64)(120.0);
+ } else {
+
+ if ( 300 <= node_number && node_number <= 395 ) {
+ def_mem = (Eint64)(27.0);
+ } else {
+ def_mem = (Eint64)(12.0);
+ }
+
+ }
+#endif
+
+ def_mem = def_mem * 1024.0;
+ fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
+
+ ok = 0;
+
+ /* Initialize the query */
+
+ queryObject = ll_query(MACHINES);
+
+ if (!queryObject) {
+ fprintf(stderr, "Query MACHINES: ll_query() returns NULL.\n");
+ freemem = def_mem;
+ return(freemem);
+ }
+
+ char **nodelist;
+ nodelist = (char **) malloc(2*sizeof(char *));
+ char nodename[6];
+ strcpy(nodename, node);
+
+ nodelist[0] = nodename;
+ nodelist[1] = NULL;
+
+ rc = ll_set_request(queryObject, QUERY_HOST, nodelist, ALL_DATA);
+
+ if (rc) {
+ fprintf(stderr, "Query MACHINES: ll_set_request() return code is non-zero.\n");
+ freemem = def_mem;
+ return(freemem);
+ }
+
+ job = ll_get_objs(queryObject, LL_CM, NULL, &obj_count, &err_code);
+
+ /* Process the job objects */
+
+ rc = ll_get_data(job, LL_MachineFreeRealMemory64, &free);
+
+ if (rc) {
+ fprintf(stderr, "Get data: ll_get_data() return code is non-zero: %d\n", rc);
+ freemem = def_mem;
+ return(freemem);
+ }
+
+/*
+ fprintf(stderr, "Free: %lld MBytes\n", free);
+*/
+
+ ll_free_objs(queryObject);
+ ll_deallocate(queryObject);
+
+ freemem = (Eint64)(free);
+
+ return(freemem);
+}
+
+#else
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "macros_and_parameters.h"
+
+Eint64 FreeRealMem( char *node )
+{
+
+ Eint64 freemem;
+ Eint64 def_mem;
+
+ def_mem = (Eint64)(12.5); /* minimum memory of DataStar */
+
+ Eint32 node_number;
+
+
+#ifdef SEABORG
+ sscanf(node, "s%5d", &node_number);
+ fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
+
+ if ( 701 <= node_number && node_number <= 713 ) {
+ def_mem = (Eint64)(56.0);
+ } else {
+
+ if ( 3701 <= node_number && node_number <= 5213 ) {
+ def_mem = (Eint64)(28.0);
+ } else {
+ def_mem = (Eint64)(14.0);
+ }
+
+ }
+#endif
+
+#ifdef DATASTAR
+ sscanf(node, "ds%3d", &node_number);
+ fprintf(stderr, "DataStarNode %s number %d\n", node, node_number);
+
+ if ( 4 <= node_number && node_number <= 11 ) {
+ def_mem = (Eint64)(120.0);
+ } else {
+
+ if ( 300 <= node_number && node_number <= 395 ) {
+ def_mem = (Eint64)(27.0);
+ } else {
+ def_mem = (Eint64)(12.0);
+ }
+
+ }
+#endif
+
+ def_mem = def_mem * 1024.0;
+ fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
+
+ freemem = def_mem;
+
+ return(freemem);
+}
+
+#endif
+
+#else
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "macros_and_parameters.h"
+
+Eint64 FreeRealMem( char *node )
+{
+
+ Eint64 freemem;
+
+ fprintf(stderr, "Error in FreeRealMem - the TASKMAP macro is not defined! Default to 12.0 GBytes.\n");
+
+ freemem = (Eint64)(12.0*1024.0); /* minimum memory on DataStar */
+
+ return(freemem);
+}
+
+#endif
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FreeRealMem.c
--- a/src/enzo/FreeRealMem.c
+++ /dev/null
@@ -1,199 +0,0 @@
-#ifdef TASKMAP
-
-#ifdef USE_LL
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include "llapi.h"
-
-#include "macros_and_parameters.h"
-
-
-Eint64 FreeRealMem( char *node )
-{
-
- LL_element *queryObject;
- LL_element *job, *step;
-
- Eint32 rc, obj_count, err_code, ok;
-
- Eint64 free;
- Eint64 freemem;
- Eint64 def_mem;
-
- def_mem = (Eint64)(12.0); /* minimum memory on DataStar */
-
- Eint32 node_number;
-
-#ifdef SEABORG
- sscanf(node, "s%5d", &node_number);
- fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
-
- if ( 701 <= node_number && node_number <= 713 ) {
- def_mem = (Eint64)(56.0);
- } else {
-
- if ( 3701 <= node_number && node_number <= 5213 ) {
- def_mem = (Eint64)(28.0);
- } else {
- def_mem = (Eint64)(14.0);
- }
-
- }
-#endif
-
-#ifdef DATASTAR
- sscanf(node, "ds%3d", &node_number);
- fprintf(stderr, "DataStar Node %s number %d\n", node, node_number);
-
- if ( 4 <= node_number && node_number <= 11 ) {
- def_mem = (Eint64)(120.0);
- } else {
-
- if ( 300 <= node_number && node_number <= 395 ) {
- def_mem = (Eint64)(27.0);
- } else {
- def_mem = (Eint64)(12.0);
- }
-
- }
-#endif
-
- def_mem = def_mem * 1024.0;
- fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
-
- ok = 0;
-
- /* Initialize the query */
-
- queryObject = ll_query(MACHINES);
-
- if (!queryObject) {
- fprintf(stderr, "Query MACHINES: ll_query() returns NULL.\n");
- freemem = def_mem;
- return(freemem);
- }
-
- char **nodelist;
- nodelist = (char **) malloc(2*sizeof(char *));
- char nodename[6];
- strcpy(nodename, node);
-
- nodelist[0] = nodename;
- nodelist[1] = NULL;
-
- rc = ll_set_request(queryObject, QUERY_HOST, nodelist, ALL_DATA);
-
- if (rc) {
- fprintf(stderr, "Query MACHINES: ll_set_request() return code is non-zero.\n");
- freemem = def_mem;
- return(freemem);
- }
-
- job = ll_get_objs(queryObject, LL_CM, NULL, &obj_count, &err_code);
-
- /* Process the job objects */
-
- rc = ll_get_data(job, LL_MachineFreeRealMemory64, &free);
-
- if (rc) {
- fprintf(stderr, "Get data: ll_get_data() return code is non-zero: %d\n", rc);
- freemem = def_mem;
- return(freemem);
- }
-
-/*
- fprintf(stderr, "Free: %lld MBytes\n", free);
-*/
-
- ll_free_objs(queryObject);
- ll_deallocate(queryObject);
-
- freemem = (Eint64)(free);
-
- return(freemem);
-}
-
-#else
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "macros_and_parameters.h"
-
-Eint64 FreeRealMem( char *node )
-{
-
- Eint64 freemem;
- Eint64 def_mem;
-
- def_mem = (Eint64)(12.5); /* minimum memory of DataStar */
-
- Eint32 node_number;
-
-
-#ifdef SEABORG
- sscanf(node, "s%5d", &node_number);
- fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
-
- if ( 701 <= node_number && node_number <= 713 ) {
- def_mem = (Eint64)(56.0);
- } else {
-
- if ( 3701 <= node_number && node_number <= 5213 ) {
- def_mem = (Eint64)(28.0);
- } else {
- def_mem = (Eint64)(14.0);
- }
-
- }
-#endif
-
-#ifdef DATASTAR
- sscanf(node, "ds%3d", &node_number);
- fprintf(stderr, "DataStarNode %s number %d\n", node, node_number);
-
- if ( 4 <= node_number && node_number <= 11 ) {
- def_mem = (Eint64)(120.0);
- } else {
-
- if ( 300 <= node_number && node_number <= 395 ) {
- def_mem = (Eint64)(27.0);
- } else {
- def_mem = (Eint64)(12.0);
- }
-
- }
-#endif
-
- def_mem = def_mem * 1024.0;
- fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
-
- freemem = def_mem;
-
- return(freemem);
-}
-
-#endif
-
-#else
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "macros_and_parameters.h"
-
-Eint64 FreeRealMem( char *node )
-{
-
- Eint64 freemem;
-
- fprintf(stderr, "Error in FreeRealMem - the TASKMAP macro is not defined! Default to 12.0 GBytes.\n");
-
- freemem = (Eint64)(12.0*1024.0); /* minimum memory on DataStar */
-
- return(freemem);
-}
-
-#endif
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -11,6 +11,7 @@
************************************************************************/
#ifndef GRID_DEFINED__
#define GRID_DEFINED__
+#include <vector>
#include "ProtoSubgrid.h"
#include "ListOfParticles.h"
#include "region.h"
@@ -20,7 +21,6 @@
#include "Star.h"
#include "FOF_allvars.h"
#include "MemoryPool.h"
-#include "list.h"
#include "hydro_rk/SuperNova.h"
#ifdef ECUDA
#include "hydro_rk/CudaMHD.h"
@@ -3013,11 +3013,10 @@
int MHDCT_ConvertEnergyToConservedS();
int MHDCT_ConvertEnergyToSpecificS();
- //List of SuperNova objects that each grid needs to keep track of
-
- List<SuperNova> SuperNovaList;
-
-
+ // used if UseMagneticSupernovaFeedback is turned on
+ int AddMagneticSupernovaeToList();
+ //List of SuperNova objects that each grid needs to keep track of
+ std::vector<SuperNova> MagneticSupernovaList;
};
// inline int grid::ReadRandomForcingFields (FILE *main_file_pointer);
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_AddFeedbackSphere.C
--- a/src/enzo/Grid_AddFeedbackSphere.C
+++ b/src/enzo/Grid_AddFeedbackSphere.C
@@ -43,8 +43,6 @@
int FindField(int field, int farray[], int numfields);
-void mt_init(unsigned_int seed);
-
int grid::AddFeedbackSphere(Star *cstar, int level, float radius, float DensityUnits,
float LengthUnits, float VelocityUnits,
float TemperatureUnits, float TimeUnits, double EjectaDensity,
@@ -990,48 +988,7 @@
/* Now it's done, unmark. */
//cstar->FeedbackFlag = NO_FEEDBACK;
-
-
- // Create a randomly-oriented SuperNova object and add it to the SuperNova Grid list
- if (cstar->FeedbackFlag == SUPERNOVA_SEEDFIELD) {
- if(UseSupernovaSeedFieldSourceTerms){
-
- SuperNova *P = new SuperNova();
- int mt_seed = 100;
- mt_init(mt_seed);
- float random_u = (float)mt_seed/100.0; // random variable from 0 to 1
- float random_v = (float)mt_seed/100.0;
- float random_phi = 2*M_PI*random_u; // 0 to 2pi
- float random_theta = acos(2*random_v-1); // 0 to pi
- // Setting up randomly oriented magnetic feedback of supernova
- float phi_x = sin(random_theta)*cos(random_phi);
- float phi_y = sin(random_theta)*sin(random_phi);
- float phi_z = cos(random_theta);
-
- // Birthtime of supernova is at the end of a star particle's life
- float sn_birthtime = cstar->BirthTime + cstar->LifeTime;
-
- // Convert units to system units
- // Converting time from years to seconds, then internal units
- float sn_duration = SupernovaSeedFieldDuration * 3.1556952e7 / TimeUnits;
- // Converting radius from parsecs to cm, then internal units
- float sn_radius = SupernovaSeedFieldRadius * 3.0856775714e18 / LengthUnits;
- // Converting energy from ergs to internal units
- float MassUnits = DensityUnits * POW(LengthUnits, 3);
- float sn_energy = SupernovaSeedFieldEnergy / (MassUnits*VelocityUnits*VelocityUnits);
-
- // Creates a supernova with magnetic feedback set by user-defined parameters and
- // adds it to supernova list
- if((Time > sn_birthtime) && (Time < sn_birthtime + SupernovaSeedFieldDuration)){
- P->setValues(phi_x, phi_y, phi_z, cstar->pos[0], cstar->pos[1], cstar->pos[2],
- sn_radius, sn_birthtime, sn_duration, sn_energy);
-
- this->SuperNovaList.append(P);
- printf("Added supernova to list \n");
- }
- }
- }
-
+
return SUCCESS;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_AddMagneticSupernovaeToList.C
--- /dev/null
+++ b/src/enzo/Grid_AddMagneticSupernovaeToList.C
@@ -0,0 +1,112 @@
+/***********************************************************************
+/
+/ FIND DEAD STAR PARTICLES AND ADD THEM TO MAGNETIC SUPERNOVA LIST
+/
+/ written by: Iryna Butsky
+/ date: June 2018
+/
+/ PURPOSE: To apply magnetic supernova feedback, we need to first add
+/ all stars that have recently gone supernova to the global list.
+/
+************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "phys_constants.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "TopGridData.h"
+#include "Grid.h"
+#include "hydro_rk/SuperNova.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+ float *TemperatureUnits, float *TimeUnits,
+ float *VelocityUnits, double *MassUnits, FLOAT Time);
+
+void mt_init(unsigned_int seed);
+unsigned_long_int mt_random();
+
+int grid::AddMagneticSupernovaeToList()
+{
+
+ if (ProcessorNumber != MyProcessorNumber) {
+ return SUCCESS;
+ }
+
+ float random_u, random_v, random_phi, random_theta, phi_x, phi_y, phi_z;
+ float sn_birthtime, sn_duration, sn_radius, sn_energy, star_birthtime, star_lifetime;
+
+ float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, VelocityUnits, MassUnits, EnergyUnits;
+ if (GetUnits(&DensityUnits, &LengthUnits,&TemperatureUnits, &TimeUnits,
+ &VelocityUnits, &MassUnits, Time) == FAIL){
+ fprintf(stderr, "Error in GetUnits.\n");
+ return FAIL;
+ }
+
+
+ // if UseMagneticSupernovaFeedback > 1, then we set the magnetic feedback radius and duration
+ // below based on the resolution of the grid at the highest refinement level
+ if (UseMagneticSupernovaFeedback > 1){
+ sn_duration = 5.0 * dtFixed;
+ MagneticSupernovaDuration = sn_duration * TimeUnits / SecondsPerYear;
+ sn_radius = 1.5 * this->CellWidth[0][0];
+ MagneticSupernovaRadius = sn_radius * LengthUnits / pc;
+ }
+ else {
+ // Converting time from years to seconds, then internal units
+ sn_duration = MagneticSupernovaDuration * SecondsPerYear / TimeUnits;
+ if(sn_duration < 5.0 * dtFixed) {
+ printf("WARNING: Magnetic supernova feedback duration is less than the recommended minimum of 5 x dtFixed\n");
+ printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / SecondsPerYear);
+ }
+ // Converting radius from parsecs to cm, then internal units
+ sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
+ if(sn_radius < 1.5 * this->CellWidth[0][0]){
+ printf("WARNING: Magnetic supernova feedback radius is less than the recommended minimum of 1.5 x CellWidth\n");
+ printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc);
+ }
+
+ }
+ // Converting energy from ergs to internal units
+ sn_energy = MagneticSupernovaEnergy / (MassUnits*VelocityUnits*VelocityUnits);
+
+ for (int star_id = 0; star_id < this->NumberOfParticles; star_id++){
+ if (this->ParticleType[star_id] == PARTICLE_TYPE_STAR){
+
+ star_birthtime = this->ParticleAttribute[0][star_id];
+ star_lifetime = this->ParticleAttribute[1][star_id];
+ sn_birthtime = star_birthtime + star_lifetime;
+
+ // Creates a supernova with magnetic feedback set by user-defined parameters and
+ // adds it to supernova list
+ if((this->Time > sn_birthtime) && (this->Time < sn_birthtime + sn_duration)){
+ SuperNova P = SuperNova();
+ mt_init((unsigned_int) this->ParticleNumber[star_id]);
+
+ random_u = (float)(mt_random()%32768)/32768.0; // random variable from 0 to 1
+ random_v = (float)(mt_random()%32768)/32768.0;
+ random_phi = 2*M_PI*random_u; // 0 to 2pi
+ random_theta = acos(2*random_v-1); // 0 to pi
+
+ // Setting up randomly oriented magnetic feedback of supernova
+ phi_x = sin(random_theta)*cos(random_phi);
+ phi_y = sin(random_theta)*sin(random_phi);
+ phi_z = cos(random_theta);
+
+ P.setValues(phi_x, phi_y, phi_z,
+ ParticlePosition[0][star_id], ParticlePosition[1][star_id], ParticlePosition[2][star_id],
+ sn_radius, sn_birthtime, sn_duration, sn_energy);
+ this->MagneticSupernovaList.push_back(P);
+ }
+ }
+ }
+ return SUCCESS;
+}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -80,19 +80,42 @@
FlaggingField[i] = 0;
float UniformParticleMass = 0.0;
- if (ProblemType == 30 && MustRefineParticlesCreateParticles == 3)
+ if (ProblemType == 30 &&
+ (MustRefineParticlesCreateParticles == 3 ||
+ MustRefineParticlesCreateParticles == 4))
UniformParticleMass = OmegaDarkMatterNow / OmegaMatterNow;
/* Loop over all particles, marking wich ones are must refine
To add rules, modify number of rules here and add to loop below */
- bool *rules;
+ bool *rules, rule0;
const int NumberOfRules = 2;
rules = new bool[NumberOfRules];
+ // Rules to prevent refinement, cancelling out the above rules.
+ bool *antirules;
+ int *AntiFlaggingField;
+ int NumberOfAntiRules = 0;
+ antirules = new bool[NumberOfAntiRules];
+
+ // Add an antirule to unflag over-refined dark matter particles.
+ if (MustRefineParticlesCreateParticles == 4) {
+ NumberOfAntiRules++;
+ }
+
+ if (NumberOfAntiRules > 0) {
+ antirules = new bool[NumberOfAntiRules];
+ AntiFlaggingField = new int[size];
+ for (i = 0; i < size; i++)
+ AntiFlaggingField[i] = 0;
+ }
+
// Flag particles as must refine particles
- int *IsParticleMustRefine;
+ int *IsParticleMustRefine, *IsParticleNotMustRefine;
+ bool OriginalParticle;
IsParticleMustRefine = new int[NumberOfParticles];
- bool OriginalParticle;
+ if (NumberOfAntiRules > 0) {
+ IsParticleNotMustRefine = new int[NumberOfParticles];
+ }
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
if (NumberOfParticleAttributes > 0)
@@ -105,21 +128,34 @@
// that they are either inert stellar remnants or "leftovers" from
// particle merging
- rules[0] = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
+ // check particle type and uniform mass
+ rule0 = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
ParticleType[i] == PARTICLE_TYPE_MBH) ||
(ParticleMass[i] < UniformParticleMass &&
ParticleType[i] == PARTICLE_TYPE_DARK_MATTER && OriginalParticle);
+ rules[0] = rule0;
// check particle mass greater than minimum mass
rules[1] = (ParticleMass[i] > MustRefineParticlesMinimumMass);
// add more rules here
+
//
// set flag for this particle
for (int j = 0; j < NumberOfRules; j++)
IsParticleMustRefine[i] *= rules[j];
+ // anti-rules
+ if (NumberOfAntiRules > 0) {
+ IsParticleNotMustRefine[i] = 1;
+ // check for over-refined dark matter particles
+ antirules[0] = !rule0;
+ }
+
+ // set antiflag for this particle
+ for (int j = 0; j < NumberOfAntiRules; j++)
+ IsParticleNotMustRefine[i] *= antirules[j];
}
PFORTRAN_NAME(cic_flag)(IsParticleMustRefine,
@@ -128,6 +164,18 @@
LeftEdge, GridDimension, GridDimension+1, GridDimension+2,
&CellSize, &ParticleBufferSize);
+ if (NumberOfAntiRules > 0) {
+ PFORTRAN_NAME(cic_flag)(IsParticleNotMustRefine,
+ ParticlePosition[0], ParticlePosition[1], ParticlePosition[2],
+ &GridRank, &NumberOfParticles, AntiFlaggingField,
+ LeftEdge, GridDimension, GridDimension+1, GridDimension+2,
+ &CellSize, &ParticleBufferSize);
+
+ for (i = 0;i < size;i++) {
+ FlaggingField[i] *= !(AntiFlaggingField[i]);
+ }
+ }
+
/* Increase particle mass flagging field for definite refinement */
float MustRefineMass =
@@ -141,7 +189,9 @@
particles, and don't use the particle mass field. */
int NumberOfFlaggedCells = 0;
- if (!(ProblemType == 30 && MustRefineParticlesCreateParticles == 3 &&
+ if (!(ProblemType == 30 &&
+ (MustRefineParticlesCreateParticles == 3 ||
+ MustRefineParticlesCreateParticles == 4) &&
level == MustRefineParticlesRefineToLevel)) {
for (i = 0; i < size; i++)
if (FlaggingField[i]) {
@@ -170,6 +220,12 @@
delete [] IsParticleMustRefine;
delete [] rules;
+ if (NumberOfAntiRules > 0) {
+ delete [] AntiFlaggingField;
+ delete [] IsParticleNotMustRefine;
+ delete [] antirules;
+ }
+
return NumberOfFlaggedCells;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_FindAllStarParticles.C
--- a/src/enzo/Grid_FindAllStarParticles.C
+++ b/src/enzo/Grid_FindAllStarParticles.C
@@ -50,7 +50,6 @@
//StarType = abs(ParticleType[i]);
if (ParticleType[i] == PARTICLE_TYPE_SINGLE_STAR ||
ParticleType[i] == PARTICLE_TYPE_BLACK_HOLE ||
- (ParticleType[i] == PARTICLE_TYPE_STAR && UseSupernovaSeedFieldSourceTerms) ||
ParticleType[i] == PARTICLE_TYPE_CLUSTER ||
ParticleType[i] == PARTICLE_TYPE_COLOR_STAR ||
ParticleType[i] == PARTICLE_TYPE_MBH ||
@@ -58,7 +57,7 @@
(StarParticleRadiativeFeedback == TRUE &&
ParticleType[i] == PARTICLE_TYPE_STAR)) {
- if (ParticleType[i] == PARTICLE_TYPE_STAR)
+ if (ParticleType[i] == PARTICLE_TYPE_STAR && UseMagneticSupernovaFeedback == 0)
LifetimeFactor = 12.0;
else
LifetimeFactor = 1.0;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_FindNewStarParticles.C
--- a/src/enzo/Grid_FindNewStarParticles.C
+++ b/src/enzo/Grid_FindNewStarParticles.C
@@ -39,9 +39,6 @@
for (i = 0; i < NumberOfParticles; i++)
if (ParticleType[i] == -PARTICLE_TYPE_SINGLE_STAR ||
- (ParticleType[i] == -PARTICLE_TYPE_STAR && UseSupernovaSeedFieldSourceTerms &&
- (this->Time >= ParticleAttribute[0][i] &&
- this->Time <= ParticleAttribute[0][i]+ParticleAttribute[1][i])) ||
ParticleType[i] == -PARTICLE_TYPE_BLACK_HOLE ||
ParticleType[i] == -PARTICLE_TYPE_CLUSTER ||
ParticleType[i] == -PARTICLE_TYPE_COLOR_STAR ||
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_WalkPhotonPackage.C
--- a/src/enzo/Grid_WalkPhotonPackage.C
+++ b/src/enzo/Grid_WalkPhotonPackage.C
@@ -245,15 +245,15 @@
// solid angle associated with package (= 4 Pi/N_package[on this level])
uint64_t Hlevel = (*PP)->level, res = 2L, BRP = 12L;
- uint64_t Nlevel = (BRP * (res << (res*Hlevel-1)));
+ uint64_t Nlevel = BRP * (1L << (res*Hlevel));
double n_on_this_level = Nlevel;
double omega_package=4*M_PI/(n_on_this_level);
double dtheta = sqrt(omega_package);
if(n_on_this_level <= 0.0 || omega_package == INFINITY) {
fprintf(stdout, "%s: level = %llu\n", __FUNCTION__, level);
fprintf(stdout, "%s: 2*(*PP)->level-1 = %llu\n", __FUNCTION__, res*(*PP)->level-1);
- fprintf(stdout, "%s: 12 * (2 << (2*(*PP)->level-1)) = %llu\n", __FUNCTION__,
- BRP * (res << (res*(*PP)->level-1)));
+ fprintf(stdout, "%s: 12 * (1 << (2*(*PP)->level)) = %llu\n", __FUNCTION__,
+ BRP * (1L << (res*(*PP)->level)));
fprintf(stdout, "%s: Hlevel = %llu\n", __FUNCTION__, Hlevel);
fprintf(stdout, "%s: Nlevel = %llu\n", __FUNCTION__, Nlevel);
fprintf(stdout, "%s: n_on_this_level = %lf\n", __FUNCTION__, n_on_this_level);
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_constructor.C
--- a/src/enzo/Grid_constructor.C
+++ b/src/enzo/Grid_constructor.C
@@ -200,12 +200,11 @@
MagneticDims[field][dim] = -100;
}}
-
/* For once-per-rootgrid-timestep star formation, the following flag
determines whether SF is about to occur or not. It's currently
(April 2012) only implemented for H2REG_STAR and completely
ignored for all other star makers. */
MakeStars = 0;
- if (UseSupernovaSeedFieldSourceTerms == 1) List<SuperNova> SuperNovaList;
+ if (UseMagneticSupernovaFeedback) std::vector<SuperNova> MagneticSupernovaList;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Make.config.objects
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -307,6 +307,7 @@
Grid_AddFeedbackSphere.o \
Grid_AddFieldMassToMassFlaggingField.o \
Grid_AddFields.o \
+ Grid_AddMagneticSupernovaeToList.o \
Grid_AddOneParticleFromList.o \
Grid_AddOverlappingParticleMassField.o \
Grid_AddParticlesFromList.o \
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Make.mach.ncsa-bluewaters-gnu
--- a/src/enzo/Make.mach.ncsa-bluewaters-gnu
+++ b/src/enzo/Make.mach.ncsa-bluewaters-gnu
@@ -69,7 +69,7 @@
MACH_CXXFLAGS =
MACH_FFLAGS = -fno-second-underscore -m64
MACH_F90FLAGS = -fno-second-underscore -m64
-MACH_LDFLAGS = -Bdynamic
+MACH_LDFLAGS = -Bstatic
#-----------------------------------------------------------------------
# Optimization flags
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/PhotonGrid_Methods.h
--- a/src/enzo/PhotonGrid_Methods.h
+++ b/src/enzo/PhotonGrid_Methods.h
@@ -102,7 +102,7 @@
int RadiativeTransferXRays(PhotonPackageEntry **PP, FLOAT *dPi, int cellindex,
int species, FLOAT ddr, FLOAT tau, float geo_correction,
- float photonrate, float *excessrate, float *ion_factor2,
+ FLOAT photonrate, FLOAT *excessrate, float *ion_factor2,
float heat_factor, const int *kphNum, int gammaNum);
int RadiativeTransferComptonHeating(PhotonPackageEntry **PP, FLOAT *dPi, int cellindex,
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/RadiativeTransferInitialize.C
--- a/src/enzo/RadiativeTransferInitialize.C
+++ b/src/enzo/RadiativeTransferInitialize.C
@@ -78,7 +78,7 @@
/* Check for radiation fields and delete them */
- NumberOfObsoleteFields = 10;
+ NumberOfObsoleteFields = 13;
ObsoleteFields[0] = kphHI;
ObsoleteFields[1] = PhotoGamma;
ObsoleteFields[2] = kphHeI;
@@ -89,6 +89,9 @@
ObsoleteFields[7] = kphHM;
ObsoleteFields[8] = kdissH2II;
ObsoleteFields[9] = RaySegments;
+ ObsoleteFields[10] = RadPressure0;
+ ObsoleteFields[11] = RadPressure1;
+ ObsoleteFields[12] = RadPressure2;
for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel)
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/RadiativeTransferXRays.C
--- a/src/enzo/RadiativeTransferXRays.C
+++ b/src/enzo/RadiativeTransferXRays.C
@@ -25,7 +25,7 @@
#include "Grid.h"
#include "CosmologyParameters.h"
-int grid::RadiativeTransferXRays(PhotonPackageEntry **PP, float *dPXray, int cellindex,
+int grid::RadiativeTransferXRays(PhotonPackageEntry **PP, FLOAT *dPXray, int cellindex,
int species, FLOAT ddr, FLOAT tau, float geo_correction,
FLOAT photonrate, FLOAT *excessrate, float *ion2_factor,
float heat_factor, const int *kphNum, int gammaNum)
@@ -57,7 +57,7 @@
#define COMPTON 3
-int grid::RadiativeTransferComptonHeating(PhotonPackageEntry **PP, float *dPXray, int cellindex,
+int grid::RadiativeTransferComptonHeating(PhotonPackageEntry **PP, FLOAT *dPXray, int cellindex,
float LengthUnits, float photonrate,
int TemperatureField, FLOAT ddr, double dN,
float geo_correction, int gammaNum)
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1285,10 +1285,10 @@
ret += sscanf(line, "GasDragCoefficient = %"GSYM, &GasDragCoefficient);
// Parameters for magnetic feedback from supernovae
- ret += sscanf(line, "UseSupernovaSeedFieldSourceTerms = %"ISYM, &UseSupernovaSeedFieldSourceTerms);
- ret += sscanf(line,"SupernovaSeedFieldRadius = %"FSYM, &SupernovaSeedFieldRadius);
- ret += sscanf(line,"SupernovaSeedFieldEnergy = %"FSYM, &SupernovaSeedFieldEnergy);
- ret += sscanf(line,"SupernovaSeedFieldDuration = %"FSYM, &SupernovaSeedFieldDuration);
+ ret += sscanf(line, "UseMagneticSupernovaFeedback = %"ISYM, &UseMagneticSupernovaFeedback);
+ ret += sscanf(line,"MagneticSupernovaRadius = %"FSYM, &MagneticSupernovaRadius);
+ ret += sscanf(line,"MagneticSupernovaEnergy = %"FSYM, &MagneticSupernovaEnergy);
+ ret += sscanf(line,"MagneticSupernovaDuration = %"FSYM, &MagneticSupernovaDuration);
/* If the dummy char space was used, then make another. */
@@ -1657,6 +1657,14 @@
if (initialize_chemistry_data(&grackle_units) == FAIL) {
ENZO_FAIL("Error in Grackle initialize_chemistry_data.\n");
}
+
+ // Need to set these after initialize_chemistry_data since
+ // that function sets them automatically based on the tables.
+ if (FinalRedshift < grackle_data->UVbackground_redshift_off) {
+ grackle_data->UVbackground_redshift_off = FinalRedshift;
+ grackle_data->UVbackground_redshift_drop = FinalRedshift;
+ }
+
} // if (grackle_data->use_grackle == TRUE)
else {
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/RestartPhotons.C
--- a/src/enzo/RestartPhotons.C
+++ b/src/enzo/RestartPhotons.C
@@ -147,6 +147,9 @@
/* Optically thin Lyman-Werner (H2) radiation field */
+ if (AllStars == NULL)
+ return SUCCESS;
+
int NumberOfSources = 0;
Star *cstar = AllStars->NextStar;
while (cstar != NULL) {
@@ -157,12 +160,10 @@
if (RadiativeTransferOpticallyThinH2)
for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel) {
- if (Temp->GridData->InitializeTemperatureFieldForH2Shield() == FAIL) {
- ENZO_FAIL("Error in InitializeTemperatureFieldForH2Shield.\n");
- }
+ Temp->GridData->InitializeTemperatureFieldForH2Shield();
Temp->GridData->AddH2Dissociation(AllStars, NumberOfSources);
}
-
+
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -986,12 +986,12 @@
UseGasDrag = 0;
GasDragCoefficient = 0.;
- /* Supernova magnetic seed field */
+ /* Supernova magnetic feedback */
/* Default == 0 -> no magnetic field contribution */
- UseSupernovaSeedFieldSourceTerms = 0;
- SupernovaSeedFieldRadius = 0.0;
- SupernovaSeedFieldDuration = 0.0;
- SupernovaSeedFieldEnergy = 0.0;
+ UseMagneticSupernovaFeedback = 0; // 0 = off; 1+ = on; 2 = Radius and duration calculated during runtime
+ MagneticSupernovaRadius = 300.0; // Total injection radius of magnetic field in parsecs
+ MagneticSupernovaDuration = 5e4 ; // Total duration of magnetic feedback in years
+ MagneticSupernovaEnergy = 1.0e51; // Total energy (ergs) injected per star particle (supernova)
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/StarParticleAddFeedback.C
--- a/src/enzo/StarParticleAddFeedback.C
+++ b/src/enzo/StarParticleAddFeedback.C
@@ -88,10 +88,7 @@
&TimeUnits, &VelocityUnits, Time);
count = 0;
- // clear list of Supernovae at each timestep to avoid adding duplicates in Grid_AddFeedbackSphere
- if(UseSupernovaSeedFieldSourceTerms){
- LevelArray[level]->GridData->SuperNovaList.clear();
- }
+
for (cstar = AllStars; cstar; cstar = cstar->NextStar, count++) {
AddedFeedback[count] = false;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Star_Accrete.C
--- a/src/enzo/Star_Accrete.C
+++ b/src/enzo/Star_Accrete.C
@@ -29,15 +29,10 @@
int Star::Accrete(void)
{
- if (UseSupernovaSeedFieldSourceTerms == 1){
- if (this->CurrentGrid == NULL ||(this->naccretions == 0))
- return SUCCESS;
- }
- else {
- if (this->CurrentGrid == NULL ||
- (this->naccretions == 0 && fabs(this->DeltaMass) < tiny_number))
- return SUCCESS;
- }
+ if (this->CurrentGrid == NULL ||
+ (this->naccretions == 0 && fabs(this->DeltaMass) < tiny_number))
+ return SUCCESS;
+
const double Msun = 1.989e33, yr = 3.1557e7;
int dim, i, n, count;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Star_HitEndpoint.C
--- a/src/enzo/Star_HitEndpoint.C
+++ b/src/enzo/Star_HitEndpoint.C
@@ -101,9 +101,6 @@
case PopIII_CF:
break;
- case NormalStar:
- break;
-
} // ENDSWITCH
return result;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Star_SetFeedbackFlag.C
--- a/src/enzo/Star_SetFeedbackFlag.C
+++ b/src/enzo/Star_SetFeedbackFlag.C
@@ -135,12 +135,7 @@
//this->type = abs_type;
- case PARTICLE_TYPE_STAR:
- if(UseSupernovaSeedFieldSourceTerms)
- this->FeedbackFlag = SUPERNOVA_SEEDFIELD;
- break;
-}
-
+ }
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -1210,10 +1210,10 @@
fprintf(fptr, "CorrectParentBoundaryFlux = %d\n", CorrectParentBoundaryFlux);
/* Supernova magnetic seed field */
- fprintf(fptr, "UseSupernovaSeedFieldSourceTerms = %"ISYM"\n", UseSupernovaSeedFieldSourceTerms);
- fprintf(fptr, "SupernovaSeedFieldRadius = %"GSYM"\n", SupernovaSeedFieldRadius);
- fprintf(fptr,"SupernovaSeedFieldEnergy = %"GSYM"\n",SupernovaSeedFieldEnergy);
- fprintf(fptr,"SupernovaSeedFieldDuration = %"GSYM"\n",SupernovaSeedFieldDuration);
+ fprintf(fptr, "UseMagneticSupernovaFeedback = %"ISYM"\n", UseMagneticSupernovaFeedback);
+ fprintf(fptr, "MagneticSupernovaRadius = %"GSYM"\n", MagneticSupernovaRadius);
+ fprintf(fptr,"MagneticSupernovaEnergy = %"GSYM"\n",MagneticSupernovaEnergy);
+ fprintf(fptr,"MagneticSupernovaDuration = %"GSYM"\n",MagneticSupernovaDuration);
/* Output current time */
time_t ID;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -1114,9 +1114,9 @@
EXTERN float GalaxySimulationPreWindVelocity[MAX_DIMENSION];
/* Supernova magnetic seed field */
-EXTERN int UseSupernovaSeedFieldSourceTerms;
-EXTERN float SupernovaSeedFieldRadius;
-EXTERN float SupernovaSeedFieldDuration;
-EXTERN float SupernovaSeedFieldEnergy;
+EXTERN int UseMagneticSupernovaFeedback;
+EXTERN float MagneticSupernovaRadius;
+EXTERN float MagneticSupernovaDuration;
+EXTERN float MagneticSupernovaEnergy;
#endif
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/hydro_rk/Grid_MHDSourceTerms.C
--- a/src/enzo/hydro_rk/Grid_MHDSourceTerms.C
+++ b/src/enzo/hydro_rk/Grid_MHDSourceTerms.C
@@ -19,13 +19,14 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-#include "list.h"
+#include "phys_constants.h"
#include "Fluxes.h"
#include "GridList.h"
#include "ExternalBoundary.h"
#include "TopGridData.h"
#include "Grid.h"
#include "EOS.h"
+#include "hydro_rk/SuperNova.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -33,6 +34,8 @@
int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
int FindField(int field, int farray[], int numfields);
+void mt_init(unsigned_int seed);
+unsigned_long_int mt_random();
int grid::MHDSourceTerms(float **dU)
{
@@ -357,49 +360,73 @@
}
}
-if((UseSupernovaSeedFieldSourceTerms == 1)) {
+ if(UseMagneticSupernovaFeedback) {
+
+ int n, active_x, active_y, center_i, center_j, center_k, num_sn_cells_x, num_sn_cells_y, num_sn_cells_z;
+ snsf_source_terms S;
+ float dx, dy, dz, dist_to_sn, magnetic_energy_density;
+ float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, VelocityUnits;
+
+ if (GetUnits(&DensityUnits, &LengthUnits,&TemperatureUnits, &TimeUnits,
+ &VelocityUnits, Time) == FAIL){
+ fprintf(stderr, "Error in GetUnits.\n");
+ return FAIL;
+ }
+
+ // Converting radius from parsecs to cm, then internal units
+ float sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
- int n = 0, igrid;
- int iden=FindField(Density, FieldType, NumberOfBaryonFields);
- snsf_source_terms S;
- List<SuperNova>::Iterator *P = this->SuperNovaList.begin();
- FLOAT cell_center[3];
- FLOAT dx, dy, dz, dist_to_sn;
- int temp =1;
- int entered = 0;
+ active_x = GridDimension[0] - 2*NumberOfGhostZones;
+ active_y = GridDimension[1] - 2*NumberOfGhostZones;
- 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++, n++) {
- igrid = i+(j+k*GridDimension[1])*GridDimension[0];
+ // the number of additional cells that receive magnetic supernova feedback.
+ // assuming the supernova is in the center of one cell
+ num_sn_cells_x = (int) ((sn_radius - 0.5*CellWidth[0][0])/ CellWidth[0][0]);
+ num_sn_cells_y = (int) ((sn_radius - 0.5*CellWidth[1][0])/ CellWidth[1][0]);
+ num_sn_cells_z = (int) ((sn_radius - 0.5*CellWidth[2][0])/ CellWidth[2][0]);
+
+ for (std::vector<SuperNova>::iterator current_sn = this->MagneticSupernovaList.begin();
+ current_sn != this->MagneticSupernovaList.end(); current_sn++){
- cell_center[0] = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
- cell_center[1] = CellLeftEdge[1][j] + 0.5*CellWidth[1][j];
- cell_center[2] = CellLeftEdge[2][k] + 0.5*CellWidth[2][k];
+ // find index of the cell nearest to the supernova center
+ // assuming that supernova in the center of that cell
+ center_i = (int)((current_sn->getPosition()[0] - GridLeftEdge[0]) / CellWidth[0][0]);
+ center_j = (int)((current_sn->getPosition()[1] - GridLeftEdge[1]) / CellWidth[1][0]);
+ center_k = (int)((current_sn->getPosition()[2] - GridLeftEdge[2]) / CellWidth[2][0]);
- while(P != this->SuperNovaList.end()){
- dx = P->get()->getPosition()[0] - cell_center[0];
- dy = P->get()->getPosition()[1] - cell_center[1];
- dz = P->get()->getPosition()[2] - cell_center[2];
-
- dist_to_sn = sqrt(dx*dx + dy*dy + dz*dz);
- if (dist_to_sn < 1.1*SupernovaSeedFieldRadius){
- S = P->get()->getSourceTerms(dx, dy, dz, Time);
- double rho = BaryonField[DensNum][igrid];
+ for(int k = center_k - num_sn_cells_z; k <= center_k + num_sn_cells_z; k++){
+ for(int j = center_j - num_sn_cells_y; j <= center_j + num_sn_cells_y; j++){
+ for(int i = center_i - num_sn_cells_x; i <= center_i + num_sn_cells_x; i++){
+
+ // only add magnetic feedback on the active grid cells
+ if ((k >= GridStartIndex[2]) && (k <= GridEndIndex[2]) &&
+ (j >= GridStartIndex[1]) && (j <= GridEndIndex[1]) &&
+ (i >= GridStartIndex[0]) && (i <= GridEndIndex[0])){
+
+ dx = CellWidth[0][0] * (float)(i-center_i);
+ dy = CellWidth[1][0] * (float)(j-center_j);
+ dz = CellWidth[2][0] * (float)(k-center_k);
+
+ dist_to_sn = sqrt(dx*dx + dy*dy + dz*dz);
+ S = current_sn->getSourceTerms(dx, dy, dz, Time);
+
+ // solving for index n
+ // analogous to how igrid is calculated, but taking into acount Ghost Zones
+ n = (i - GridStartIndex[0])+((j-GridStartIndex[1]) + (k-GridStartIndex[2])*active_y) * active_x;
+
+ dU[iBx][n] += S.dbx*dtFixed;
+ dU[iBy][n] += S.dby*dtFixed;
+ dU[iBz][n] += S.dbz*dtFixed;
- dU[iBx][n] += S.dbx*dtFixed;
- dU[iBy][n] += S.dby*dtFixed;
- dU[iBz][n] += S.dbz*dtFixed;
- dU[iEtot][n] += S.dUtot*dtFixed;
+ dU[iEtot][n] += S.dUtot * dtFixed;
+
+ }
- }
- P = P->next();
- }// End of SuperNovaList iteration
- } // End of k for-loop
- } // End of j for-loop
- } // End of i for-loop
-
-} // End of UseSuperNovaSeedFieldSourceTerms scope
+ } // End of k for-loop
+ } // End of j for-loop
+ } // End of i for-loop
+ } // End of MagneticSupernovaList loop
+ } // End of UseMagneticSupernovaFeedback scope
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/hydro_rk/LLF_Zero_MHD.C
--- a/src/enzo/hydro_rk/LLF_Zero_MHD.C
+++ b/src/enzo/hydro_rk/LLF_Zero_MHD.C
@@ -34,7 +34,7 @@
// compute priml and primr
int iprim;
- for (int field = 0; field < NEQ_MHD; field++) {
+ for (int field = 0; field < NEQ_MHD-DualEnergyFormalism; field++) {
for (int i = 0; i < ActiveSize+1; i++) {
iprim = i + NumberOfGhostZones - 1;
priml[field][i] = prim[field][iprim];
@@ -53,9 +53,9 @@
iprim = i+NumberOfGhostZones-1;
for (int field = 0; field < NSpecies; field++) {
if (FluxLine[iD][i] >= 0) {
- species[field][i] = prim[field+5][iprim ];
+ species[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism][iprim ];
} else {
- species[field][i] = prim[field+5][iprim+1];
+ species[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism][iprim+1];
}
}
sum = 0;
@@ -76,9 +76,9 @@
iprim = i+NumberOfGhostZones-1;
for (int field = 0; field < NColor; field++) {
if (FluxLine[iD][i] >= 0) {
- colors[field][i] = prim[field+5+NSpecies][iprim ];
+ colors[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism+NSpecies][iprim ];
} else {
- colors[field][i] = prim[field+5+NSpecies][iprim+1];
+ colors[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism+NSpecies][iprim+1];
}
}
for (int field = 0; field < NColor; field++) {
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/hydro_rk/SuperNovaSeedField.C
--- a/src/enzo/hydro_rk/SuperNovaSeedField.C
+++ b/src/enzo/hydro_rk/SuperNovaSeedField.C
@@ -127,14 +127,14 @@
snsf_source_terms SuperNova::getSourceTerms(double dx, double dy, double dz, double enzoTime){
// ------------------------------------------------------------------- //
- // Inputs: (1) distances to current cell, dx, dy, dz //
+ // Inputs: (1) distances to current cell, dx, dy, dz //
// (2) current time of the simulation //
// //
- // Output: (1) A struct of source terms corresponding to the given //
- // supernova //
- // (i) x,y,z components of time derivative of B-field //
- // (ii) time derivative of magnetic energy density //
- // (iii) time derivative of total energy density //
+ // Output: (1) A struct of source terms corresponding to the given //
+ // supernova //
+ // (i) x,y,z components of time derivative of B-field //
+ // (ii) time derivative of magnetic energy density //
+ // (iii) time derivative of total energy density //
// ------------------------------------------------------------------- //
FLOAT phi[3], zhat_cell[3];
@@ -180,6 +180,9 @@
FLOAT r_scale = (r_cyl/characteristicLength)*\
exp(- r_s*r_s /(characteristicLength*characteristicLength));
+
+ FLOAT db_r_scale = (r_cyl/characteristicLength)*\
+ exp(- r_s*r_s /(2.0*characteristicLength*characteristicLength));
FLOAT t_exp = exp(-t_s/characteristicTime);
FLOAT db_t_exp = t_exp / characteristicTime;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -1,5 +1,12 @@
#ifndef __macros_and_parameters_h_
#define __macros_and_parameters_h_
+#include "string.h"
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <algorithm>
+#include <map>
+#include <cassert>
/***********************************************************************
/
/ MACRO DEFINITIONS AND PARAMETERS
@@ -522,7 +529,6 @@
#define MBH_THERMAL 7
#define MBH_JETS 8
#define COLOR_FIELD 9
-#define SUPERNOVA_SEEDFIELD 11
/* Sink particle accretion modes */
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/phys_constants.h
--- a/src/enzo/phys_constants.h
+++ b/src/enzo/phys_constants.h
@@ -60,4 +60,7 @@
#define kpc 3.0857e21
#define pc 3.0857e18
+/* Time units */
+#define SecondsPerYear 3.1556952e7
+
#endif
https://bitbucket.org/enzo/enzo-dev/commits/a74f32a0094e/
Changeset: a74f32a0094e
Branch: week-of-code
User: jwise77
Date: 2019-01-11 17:04:53+00:00
Summary: More descriptive comment for MAX_SPLIT_ITERATIONS
Affected #: 1 file
diff -r ff98226315b7 -r a74f32a0094e src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -486,7 +486,8 @@
#define CHILDRENPERPARENT 12
-/* You probably don't want to do this too often */
+/* Splitting particles more than 4 times can induce numerical
+ artifacts */
#define MAX_SPLIT_ITERATIONS 4
https://bitbucket.org/enzo/enzo-dev/commits/bb230291c6f2/
Changeset: bb230291c6f2
Branch: week-of-code
User: jwise77
Date: 2019-02-01 17:29:40+00:00
Summary: Updating particle splitting docs to clarify some details.
Affected #: 1 file
diff -r a74f32a0094e -r bb230291c6f2 doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1889,14 +1889,17 @@
``ParticleSplitterMustRefine`` (external)
Set to 1 to mark the split particles as must-refine particles. The
user must also set associated must-refine parameters to enable its
- machinery. Default: 0
+ machinery that can be used to restrict AMR only to the must-refine
+ particles. Default: 0
``ParticleSplitterMustRefineIDFile`` (external)
Filename for the HDF5 file that has a dataset containing the
particle IDs that should be marked as must-refine. All other
particles within the region marked for splitting will retain their
- original types. This must be used in conjunction with
- ``ParticleSplitterMustRefine = 1``. The dataset must be named
- ``particle_identifier`` in the base group. Default: (null).
+ original types. If not set, all particles within the must-refine
+ region will be must-refine particles. This must be used in
+ conjunction with ``ParticleSplitterMustRefine = 1``. The dataset
+ must be named ``particle_identifier`` in the base group. Default:
+ (null).
An example yt script is provided below, selecting the particles in
a sphere centered at [0.5, 0.5, 0.5] with a radius 0.05 in code
@@ -1918,10 +1921,11 @@
``ParticleSplitterFraction`` (external)
An array of four values that represent the width of the splitting
region in units of the original refine region set by
- ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. Each
+ ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. The
+ splitting region is centered on the refine region center. Each
successive value represents the next nested split region. Valid
- up to ``ParticleSplitterIterations`` times. Default: 1.0 (all 4
- values)
+ up to ``ParticleSplitterIterations`` times. Cannot be used with
+ ``ParticleSplitterCenterRegion``. Default: 1.0 (all 4 values)
``ParticleSplitterCenter`` (external)
The center of split region in code units. Specify if the split
region does not correspond to the center of the refine region.
@@ -1929,13 +1933,14 @@
``ParticleSplitterCenterRegion`` (external)
The width of the split region in code units. Must be used in
conjunction with ``ParticleSplitterCenter``. Each successive
- value represents the next nested split region. Valid up to
+ value represents the next nested split region. Cannot be used
+ with ``ParticleSplitterFraction``. Valid up to
``ParticleSplitterIterations`` times. Not used if
negative. Default: -1.0 (all 4 values)
``ParticleSplitterChildrenParticleSeparation`` (external)
This is the spacing between the child particles placed on a
- hexagonal close-packed (HCP) array. In units of a cell size
- which the parent particle resides in. Default: 1.0
+ hexagonal close-packed (HCP) array. In units of a cell size which
+ the parent particle resides in. Default: 1.0
.. _starparticleparameters:
https://bitbucket.org/enzo/enzo-dev/commits/e58dbfd476cb/
Changeset: e58dbfd476cb
Branch: week-of-code
User: gbryan
Date: 2019-02-04 21:10:50+00:00
Summary: Merged in jwise77/clean-copy (pull request #426)
Updates to particle splitter
Approved-by: John Regan <johnanth...@gmail.com>
Approved-by: Greg Bryan <gbr...@astro.columbia.edu>
Approved-by: Britton Smith <britto...@gmail.com>
Affected #: 12 files
diff -r 603a8a701fe4 -r e58dbfd476cb doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1879,11 +1879,68 @@
Currently it implicitly assumes that only DM (type=1) and
conventional star particles (type=2) inside the ``RefineRegion`` get
split. Other particles, which usually become Star class objects,
- seem to have no reason to be split. Default: 0
+ seem to have no reason to be split. It is reset to zero after a
+ restart to avoid resplitting in subsequent datasets. It can be set
+ to a maximum of 4. Default: 0
+``ParticleSplitterRandomSeed`` (external)
+ Random seed used when randomly rotating the hexagonal close
+ packed array on whose vertices the split particles are
+ placed. Default: 131180
+``ParticleSplitterMustRefine`` (external)
+ Set to 1 to mark the split particles as must-refine particles. The
+ user must also set associated must-refine parameters to enable its
+ machinery that can be used to restrict AMR only to the must-refine
+ particles. Default: 0
+``ParticleSplitterMustRefineIDFile`` (external)
+ Filename for the HDF5 file that has a dataset containing the
+ particle IDs that should be marked as must-refine. All other
+ particles within the region marked for splitting will retain their
+ original types. If not set, all particles within the must-refine
+ region will be must-refine particles. This must be used in
+ conjunction with ``ParticleSplitterMustRefine = 1``. The dataset
+ must be named ``particle_identifier`` in the base group. Default:
+ (null).
+
+ An example yt script is provided below, selecting the particles in
+ a sphere centered at [0.5, 0.5, 0.5] with a radius 0.05 in code
+ length units.
+
+ .. code-block:: python
+
+ import yt
+ import h5py as h5
+
+ ds = yt.load('DD0040/DD0040')
+ center = ds.arr([0.5, 0.5, 0.5], 'code_length')
+ radius = ds.quan(0.05, 'code_length')
+ sp = ds.sphere(center, radius)
+ fp = h5.File('particle-ids.h5', 'w')
+ fp['particle_identifier'] = sp['particle_index'].astype('int')
+ fp.close()
+
+``ParticleSplitterFraction`` (external)
+ An array of four values that represent the width of the splitting
+ region in units of the original refine region set by
+ ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. The
+ splitting region is centered on the refine region center. Each
+ successive value represents the next nested split region. Valid
+ up to ``ParticleSplitterIterations`` times. Cannot be used with
+ ``ParticleSplitterCenterRegion``. Default: 1.0 (all 4 values)
+``ParticleSplitterCenter`` (external)
+ The center of split region in code units. Specify if the split
+ region does not correspond to the center of the refine region.
+ Not used if negative. Default: -1.0 -1.0 -1.0
+``ParticleSplitterCenterRegion`` (external)
+ The width of the split region in code units. Must be used in
+ conjunction with ``ParticleSplitterCenter``. Each successive
+ value represents the next nested split region. Cannot be used
+ with ``ParticleSplitterFraction``. Valid up to
+ ``ParticleSplitterIterations`` times. Not used if
+ negative. Default: -1.0 (all 4 values)
``ParticleSplitterChildrenParticleSeparation`` (external)
This is the spacing between the child particles placed on a
- hexagonal close-packed (HCP) array. In the unit of a cell size
- which the parent particle resides in. Default: 1.0
+ hexagonal close-packed (HCP) array. In units of a cell size which
+ the parent particle resides in. Default: 1.0
.. _starparticleparameters:
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/EvolveHierarchy.C
--- a/src/enzo/EvolveHierarchy.C
+++ b/src/enzo/EvolveHierarchy.C
@@ -270,20 +270,26 @@
*/
+ /* Particle Splitter. Split particles into 13 (=1+12) child
+ particles. The hierarchy is rebuilt inside this routine. */
+
+ if (MetaData.FirstTimestepAfterRestart == TRUE &&
+ ParticleSplitterIterations > 0) {
+
+ ParticleSplitter(LevelArray, 0, &MetaData);
+
+ } else {
+
/* Do the first grid regeneration. */
- if(CheckpointRestart == FALSE) {
- RebuildHierarchy(&MetaData, LevelArray, 0);
- }
+ if(CheckpointRestart == FALSE) {
+ RebuildHierarchy(&MetaData, LevelArray, 0);
+ }
+
+ } // ENDELSE particle splitting
PrintMemoryUsage("1st rebuild");
- /* Particle Splitter. Split particles into 13 (=1+12) child particles */
-
- if (MetaData.FirstTimestepAfterRestart == TRUE &&
- ParticleSplitterIterations > 0)
- ParticleSplitter(LevelArray, 0, &MetaData);
-
/* Reset magnetic fields if requested. */
if (MetaData.FirstTimestepAfterRestart == TRUE &&
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2409,13 +2409,15 @@
/* Particle splitter routine. */
- int ParticleSplitter(int level);
+ int ParticleSplitter(int level, int iter, int NumberOfIDs,
+ long *MustRefineIDs);
int CreateChildParticles(float dx, int NumberOfParticles, float *ParticleMass,
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles);
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles);
/* Magnetic field resetting routine. */
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -33,6 +33,8 @@
#include "fortran.def"
#include "CosmologyParameters.h"
+#define DEBUG_PS
+
void mt_init(unsigned_int seed);
unsigned_long_int mt_random();
@@ -41,13 +43,15 @@
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles)
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles)
{
- int partnum = 0, i = 0, child = 0, m = 0, numpart = 0, innerchild = 0;
+ int partnum = 0, i = 0, child = 0, m = 0, innerchild = 0;
int xindex=0, yindex=0, zindex=0;
float rad = 0.0, sin60 = 0.0;
FLOAT NewPos[3][CHILDRENPERPARENT];
+ FLOAT LeftEdge[3], RightEdge[3];
int ChildrenPerParent = CHILDRENPERPARENT, total_children = 0;
float alpha[3];
float l11 = 0.0, l12 = 0.0, l13 = 0.0;
@@ -67,16 +71,77 @@
//Reproducible random seed
mt_init(((unsigned_int) ParticleSplitterRandomSeed));
+ FLOAT sep[3], midpoint[3], newsep[3];
+
+ /*
+ * The following three options determine over what region
+ * particle splitting takes place.
+ * By default the splitting occurs over the entire most refined
+ * region. At times this can be too large so two further options
+ * exist controlled via the parameter file. The first just selects
+ * a fraction of that region centred on the centre of the most
+ * refined region.
+ * The second allows the user to specify the centre of the splitting
+ * and the size of the region.
+ */
+
+ /* Default Option */
+ for (i = 0; i < 3; i++) {
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ /* Change particle refinement region in basic way */
+ if(ParticleSplitterFraction[iter] != 1.0)
+ {
+ fprintf(stderr, "Setting Particle Refinement using fractions.\n");
+ for (i = 0; i < 3; i++) {
+ sep[i] = RefineRegionRightEdge[i] - RefineRegionLeftEdge[i];
+ midpoint[i] = sep[i]/2.0 + RefineRegionLeftEdge[i];
+ newsep[i] = sep[i]*ParticleSplitterFraction[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ }
+ }
+ /* Centre particle refinement region around a specific point */
+ if(ParticleSplitterCenter[0] > 0.0 && ParticleSplitterCenterRegion[iter] > 0.0) {
+ fprintf(stderr, "Setting Particle Refinement around point.\n");
+ for (i = 0; i < 3; i++) {
+ midpoint[i] = ParticleSplitterCenter[i];
+ newsep[i] = ParticleSplitterCenterRegion[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ if(LeftEdge[i] < RefineRegionLeftEdge[i])
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ if(RightEdge[i] > RefineRegionRightEdge[i])
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ }
+
+ #ifdef DEBUG_PS
+ fprintf(stdout, "%s: Iteration %d: midpoint = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, midpoint[0], midpoint[1], midpoint[2]);
+ fprintf(stdout, "%s: Iteration %d: newsep = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, newsep[0], newsep[1], newsep[2]);
+ fprintf(stdout, "%s: Iteration %d: Refine Particles between (%lf, %lf, %lf) and (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, LeftEdge[0], LeftEdge[1], LeftEdge[2],
+ RightEdge[0], RightEdge[1], RightEdge[2]);
+ #endif
+
+
/*
* Loop over existing (parent) particles; It implicitly assumes that
* only DM and conventional star particles get splitted. Other
* particles - which usually become Star class particles - doesn't
* seem to have any reason to be splitted. (as of Oct.2009)
*/
+
+ /* (Apr 2018) Include must-refine particles */
for(partnum = 0; partnum < NumberOfParticles; partnum++)
{
- if(ParticleMass[partnum] > 0.0 && ParticleType[partnum] <= 2)
+ if(ParticleMass[partnum] > 0.0 &&
+ (ParticleType[partnum] <= 2 ||
+ ParticleType[partnum] == PARTICLE_TYPE_MUST_REFINE))
{
/*
* Check that particle is within the most refined region.
@@ -84,12 +149,12 @@
* blatantly discriminate against non-refined types around
* here.
*/
- if(ParticlePosition[0][partnum] < RefineRegionLeftEdge[0] ||
- ParticlePosition[0][partnum] > RefineRegionRightEdge[0] ||
- ParticlePosition[1][partnum] < RefineRegionLeftEdge[1] ||
- ParticlePosition[1][partnum] > RefineRegionRightEdge[1] ||
- ParticlePosition[2][partnum] < RefineRegionLeftEdge[2] ||
- ParticlePosition[2][partnum] > RefineRegionRightEdge[2])
+ if(ParticlePosition[0][partnum] < LeftEdge[0] ||
+ ParticlePosition[0][partnum] > RightEdge[0] ||
+ ParticlePosition[1][partnum] < LeftEdge[1] ||
+ ParticlePosition[1][partnum] > RightEdge[1] ||
+ ParticlePosition[2][partnum] < LeftEdge[2] ||
+ ParticlePosition[2][partnum] > RightEdge[2])
{
//fprintf(stdout, "grid::PS: Particle outside RR - ignoring\n");
continue;
@@ -285,8 +350,13 @@
for(i = 0; i < 3; i++)
this->ParticleVelocity[i][child] = ParticleVelocity[i][partnum];
this->ParticleType[child] = ParticleType[partnum];
+ // Flag that a DM particle was split (originally -99999 or 0)
+ if (NumberOfParticleAttributes > 0 &&
+ ParticleType[partnum] == PARTICLE_TYPE_DARK_MATTER)
+ if (ParticleAttribute[0][partnum] <= 0)
+ ParticleAttribute[0][partnum] = tiny_number;
for(i = 0; i < NumberOfParticleAttributes; i++)
- this->ParticleAttribute[i][child] = ParticleAttribute[i][numpart];
+ this->ParticleAttribute[i][child] = ParticleAttribute[i][partnum];
}
/* Loop forward by CHILDRENPERPARENT each time. */
@@ -302,8 +372,9 @@
}
}
#ifdef DEBUG_PS
- fprintf(stdout,"%d new child particles created in this grid from %d parents",
- total_children, NumberOfParticles);
+ if(total_children > 0)
+ fprintf(stdout,"Iteration %d: %d new child particles created in this grid from %d parents\n",
+ iter, total_children, NumberOfParticles);
#endif
*NumberOfNewParticles = total_children;
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -111,17 +111,28 @@
// Flag particles as must refine particles
int *IsParticleMustRefine, *IsParticleNotMustRefine;
+ bool OriginalParticle;
IsParticleMustRefine = new int[NumberOfParticles];
if (NumberOfAntiRules > 0) {
IsParticleNotMustRefine = new int[NumberOfParticles];
}
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
+ if (NumberOfParticleAttributes > 0)
+ OriginalParticle = (ParticleAttribute[0][i] <= 0.0);
+ else
+ OriginalParticle = true;
+
+ // check particle type and uniform mass. Also check particle
+ // creation time for DM particles that are positive, indicating
+ // that they are either inert stellar remnants or "leftovers" from
+ // particle merging
// check particle type and uniform mass
rule0 = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
ParticleType[i] == PARTICLE_TYPE_MBH) ||
- (ParticleMass[i] < UniformParticleMass);
+ (ParticleMass[i] < UniformParticleMass &&
+ ParticleType[i] == PARTICLE_TYPE_DARK_MATTER && OriginalParticle);
rules[0] = rule0;
// check particle mass greater than minimum mass
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -40,22 +40,8 @@
float *VelocityUnits, FLOAT Time);
int FindField(int field, int farray[], int numfields);
-extern "C" void FORTRAN_NAME(particle_splitter)(int *nx, int *ny, int *nz,
- int *idual, int *imetal, hydro_method *imethod, float *dt,
- float *r, float *dx, FLOAT *t, float *z,
- float *d1, float *x1, float *v1, float *t1,
- FLOAT *xstart, FLOAT *ystart, FLOAT *zstart, int *ibuff,
- int *npart,
- FLOAT *xpold, FLOAT *ypold, FLOAT *zpold, float *upold, float *vpold, float *wpold,
- float *mpold, float *tdpold, float *tcpold, float *metalfold, int *typeold,
- int *nmax, int *npartnew, int *children, int *level,
- FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
- float *mp, float *tdp, float *tcp, float *metalf, int *type,
- int *iterations, float *separation, int *ran1_init,
- FLOAT *rr_leftedge, FLOAT *rr_rightedge);
-
-
-int grid::ParticleSplitter(int level)
+int grid::ParticleSplitter(int level, int iteration, int NumberOfIDs,
+ long *MustRefineIDs)
{
if (ParticleSplitterIterations == 0)
@@ -131,7 +117,7 @@
/* Allocate space for new particles. */
- int ChildrenPerParent = 12; //12+1 = 13 will be the final number of particles per parent
+ int ChildrenPerParent = CHILDRENPERPARENT; //12+1 = 13 will be the final number of particles per parent
int MaximumNumberOfNewParticles = ChildrenPerParent*NumberOfParticles+1;
tg->AllocateNewParticles(MaximumNumberOfNewParticles);
@@ -145,7 +131,7 @@
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles before splitting = %d, MyProcessorNumber = %d\n",
NumberOfParticles, MyProcessorNumber);
#endif
-
+
if (NumberOfParticles > 0) {
#define NO_PARTICLE_IN_GRID_CHECK
@@ -165,10 +151,29 @@
xindex, yindex, zindex, level);
}
#endif
+
+ /* Before splitting, set particle type to must-refine if given in
+ the provided list. Only do this for the first iteration because
+ in the later iterations, the particle types are propagated to
+ the split particles. */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs != NULL &&
+ iteration == 0) {
+ for (i = 0; i < NumberOfParticles; i++) {
+ for (j = 0; j < NumberOfIDs; j++) {
+ if (this->ParticleNumber[i] == MustRefineIDs[j]) {
+ this->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ break;
+ }
+ } // ENDFOR j
+ } // ENDFOR i
+ }
+
if(!tg->CreateChildParticles(CellWidthTemp, NumberOfParticles, ParticleMass,
ParticleType, ParticlePosition, ParticleVelocity,
ParticleAttribute, CellLeftEdge, GridDimension,
- MaximumNumberOfNewParticles, &NumberOfNewParticles))
+ MaximumNumberOfNewParticles, iteration,
+ &NumberOfNewParticles))
{
fprintf(stdout, "Failed to create child particles in grid %d\n", this->GetGridID());
return FAIL;
@@ -185,6 +190,14 @@
NumberOfNewParticles, MyProcessorNumber);
#endif
+ /* If specified and no ID list given, set all particle types to
+ must-refine */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs == NULL) {
+ for (i = 0; i < NumberOfNewParticles; i++)
+ tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ }
+
/* If not set in the above routine, then set the metal fraction to zero. */
if (MetallicityField == FALSE && NumberOfParticleAttributes > 1)
@@ -211,7 +224,7 @@
#ifdef DEBUG_PS
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles(New) = %d, MyProcessorNumber = %d\n",
- NumberOfParticles, MyProcessorNumber);
+ NumberOfParticles, MyProcessorNumber);
#endif
} // end: if (NumberOfNewParticles > 0)
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -24,6 +24,7 @@
#include <stdlib.h>
#include <stdio.h>
+#include <hdf5.h>
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "typedefs.h"
@@ -38,7 +39,7 @@
#include "CosmologyParameters.h"
#include "CommunicationUtilities.h"
-#define NO_DEBUG_PS
+#define NO_DEBUG_PS
int RebuildHierarchy(TopGridData *MetaData,
LevelHierarchyEntry *LevelArray[], int level);
@@ -60,15 +61,20 @@
HierarchyEntry **Grids;
int NumberOfGrids;
- /* First, rebuild the hierarchy */
-
- RebuildHierarchy(MetaData, LevelArray, 0);
-
- /* Return if this does not concern us */
+ /* Rebuild hierarchy and return if this does not concern us */
if (ParticleSplitterIterations <= 0 ||
- ParticleSplitterChildrenParticleSeparation <=0)
+ ParticleSplitterChildrenParticleSeparation <=0) {
+ RebuildHierarchy(MetaData, LevelArray, 0);
return SUCCESS;
+ }
+
+ if(ParticleSplitterIterations > MAX_SPLIT_ITERATIONS) {
+ fprintf(stderr, "WARNING: Splitting iterations exceeds maximum allowed\n" \
+ "Resetting to %d\n", MAX_SPLIT_ITERATIONS);
+ ParticleSplitterIterations = MAX_SPLIT_ITERATIONS;
+ }
+
/* Find total NumberOfParticles in all grids; this is needed in
CommunicationUpdateStarParticleCount below */
@@ -76,10 +82,35 @@
MetaData->NumberOfParticles = FindTotalNumberOfParticles(LevelArray);
NumberOfOtherParticles = MetaData->NumberOfParticles - NumberOfStarParticles;
- /* Initialize all star particles if this is a restart */
+ /* Read must-refine particle positions if specified */
+
+ long *MustRefineIDs = NULL;
+ int NumberOfIDs = 0;
+
+ if (ParticleSplitterMustRefine &&
+ ParticleSplitterMustRefineIDFile != NULL) {
- if (ParticleSplitterIterations > 1)
- fprintf(stderr, "WARNING: ParticleSplitterIterations > 1 is not properly tested yet.\n");
+ int NumberOfMustRefineIDs;
+ hid_t file_id, dataset_id, dataspace_id;
+ herr_t status;
+ hsize_t dims[1], maxdims[1];
+
+ file_id = H5Fopen(ParticleSplitterMustRefineIDFile,
+ H5F_ACC_RDONLY, H5P_DEFAULT);
+ dataset_id = H5Dopen2(file_id, "/particle_identifier", H5P_DEFAULT);
+ dataspace_id = H5Dget_space(dataset_id);
+ // Get number of particles and allocate variable
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, maxdims);
+ MustRefineIDs = new long[dims[0]];
+ NumberOfIDs = dims[0];
+ status = H5Dread(dataset_id, H5T_NATIVE_LLONG, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ MustRefineIDs);
+ status = H5Dclose(dataset_id);
+ status = H5Fclose(file_id);
+
+ }
+
+ /* Initialize all star particles if this is a restart */
for (i = 0; i < ParticleSplitterIterations; i++) {
@@ -98,7 +129,6 @@
RecordTotalStarParticleCount(Grids, NumberOfGrids,
TotalStarParticleCountPrevious);
-
// for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
// fprintf(stdout, "TotalStarParticleCountPrevious[grid=%d] = %d\n", grid1,
// TotalStarParticleCountPrevious[grid1]);
@@ -110,7 +140,8 @@
Grids[grid1]->GridData->ReturnNumberOfParticles());
#endif
- if (Grids[grid1]->GridData->ParticleSplitter(level) == FAIL) {
+ if (Grids[grid1]->GridData->ParticleSplitter(level, i, NumberOfIDs,
+ MustRefineIDs) == FAIL) {
ENZO_FAIL("Error in grid::ParticleSplitter.\n");
}
@@ -148,6 +179,7 @@
/* Set splitter parameter zero; otherwise it will split particles again at next restart */
+ delete [] MustRefineIDs;
ParticleSplitterIterations = 0;
return SUCCESS;
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1262,6 +1262,18 @@
&ParticleSplitterRandomSeed);
ret += sscanf(line, "ParticleSplitterChildrenParticleSeparation = %"FSYM,
&ParticleSplitterChildrenParticleSeparation);
+ ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
+ &ParticleSplitterMustRefine);
+ if (sscanf(line, "ParticleSplitterMustRefineIDFile = %s", dummy) == 1)
+ ParticleSplitterMustRefineIDFile = dummy;
+ ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
+ ParticleSplitterFraction+3);
+ ret += sscanf(line, "ParticleSplitterCenter = %"PSYM" %"PSYM" %"PSYM"",
+ ParticleSplitterCenter+0, ParticleSplitterCenter+1, ParticleSplitterCenter+2);
+ ret += sscanf(line, "ParticleSplitterCenterRegion = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterCenterRegion+0, ParticleSplitterCenterRegion+1,
+ ParticleSplitterCenterRegion+2, ParticleSplitterCenterRegion+3);
ret += sscanf(line, "ResetMagneticField = %"ISYM,
&ResetMagneticField);
ret += sscanf(line, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM,
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -962,6 +962,14 @@
ParticleSplitterIterations = FALSE;
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
+ ParticleSplitterMustRefine = FALSE;
+ ParticleSplitterMustRefineIDFile = NULL;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterFraction[i] = 1.0;
+ for(int i = 0; i < MAX_DIMENSION; i++)
+ ParticleSplitterCenter[i] = -1.0;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterCenterRegion[i] = -1.0;
/* Magnetic Field Resetter */
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -624,6 +624,16 @@
ParticleSplitterRandomSeed);
fprintf(fptr, "ParticleSplitterChildrenParticleSeparation = %"FSYM"\n",
ParticleSplitterChildrenParticleSeparation);
+ fprintf(fptr, "ParticleSplitterFraction = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterFraction);
+ fprintf(fptr, "ParticleSplitterCenter = ");
+ WriteListOfFloats(fptr, 3, ParticleSplitterCenter);
+ fprintf(fptr, "ParticleSplitterCenterRegion = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterCenterRegion);
+ fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
+ ParticleSplitterMustRefine);
+ fprintf(fptr, "ParticleSplitterMustRefineIDFile = %s\n",
+ ParticleSplitterMustRefineIDFile);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
ResetMagneticField);
fprintf(fptr, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM"\n",
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -984,6 +984,11 @@
EXTERN int ParticleSplitterIterations;
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
+EXTERN int ParticleSplitterMustRefine;
+EXTERN char *ParticleSplitterMustRefineIDFile;
+EXTERN float ParticleSplitterFraction[MAX_SPLIT_ITERATIONS];
+EXTERN FLOAT ParticleSplitterCenter[MAX_DIMENSION];
+EXTERN float ParticleSplitterCenterRegion[MAX_SPLIT_ITERATIONS];
/* Magnetic Field Resetter */
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -485,6 +485,12 @@
#define PARTICLE_TYPE_SIMPLE_SOURCE 10
#define CHILDRENPERPARENT 12
+
+/* Splitting particles more than 4 times can induce numerical
+ artifacts */
+
+#define MAX_SPLIT_ITERATIONS 4
+
/* Ways to deposit particles from a subgrid. */
#define CIC_DEPOSIT 0
Repository URL: https://bitbucket.org/enzo/enzo-dev/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
https://bitbucket.org/enzo/enzo-dev/commits/67819df51893/
Changeset: 67819df51893
Branch: week-of-code
User: jwise77
Date: 2018-04-10 13:42:59+00:00
Summary: Adding an option to make split particles into must-refine ones
Affected #: 5 files
diff -r 149dd1d6513b -r 67819df51893 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -185,6 +185,13 @@
NumberOfNewParticles, MyProcessorNumber);
#endif
+ /* If specified, set particle type to must-refine */
+
+ if (ParticleSplitterMustRefine == TRUE) {
+ for (i = 0; i < NumberOfNewParticles; i++)
+ tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ }
+
/* If not set in the above routine, then set the metal fraction to zero. */
if (MetallicityField == FALSE && NumberOfParticleAttributes > 1)
diff -r 149dd1d6513b -r 67819df51893 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1262,6 +1262,8 @@
&ParticleSplitterRandomSeed);
ret += sscanf(line, "ParticleSplitterChildrenParticleSeparation = %"FSYM,
&ParticleSplitterChildrenParticleSeparation);
+ ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
+ &ParticleSplitterMustRefine);
ret += sscanf(line, "ResetMagneticField = %"ISYM,
&ResetMagneticField);
ret += sscanf(line, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM,
diff -r 149dd1d6513b -r 67819df51893 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -962,6 +962,7 @@
ParticleSplitterIterations = FALSE;
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
+ ParticleSplitterMustRefine = FALSE;
/* Magnetic Field Resetter */
diff -r 149dd1d6513b -r 67819df51893 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -624,6 +624,8 @@
ParticleSplitterRandomSeed);
fprintf(fptr, "ParticleSplitterChildrenParticleSeparation = %"FSYM"\n",
ParticleSplitterChildrenParticleSeparation);
+ fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
+ ParticleSplitterMustRefine);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
ResetMagneticField);
fprintf(fptr, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM"\n",
diff -r 149dd1d6513b -r 67819df51893 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -984,6 +984,7 @@
EXTERN int ParticleSplitterIterations;
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
+EXTERN int ParticleSplitterMustRefine;
/* Magnetic Field Resetter */
https://bitbucket.org/enzo/enzo-dev/commits/cfe6395292a7/
Changeset: cfe6395292a7
Branch: week-of-code
User: jwise77
Date: 2018-04-13 16:51:38+00:00
Summary: Adding additional particle splitting functionality from John Regan's
enzo-3.0 fork.
Affected #: 9 files
diff -r 67819df51893 -r cfe6395292a7 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2409,13 +2409,14 @@
/* Particle splitter routine. */
- int ParticleSplitter(int level);
+ int ParticleSplitter(int level, int iter);
int CreateChildParticles(float dx, int NumberOfParticles, float *ParticleMass,
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles);
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles);
/* Magnetic field resetting routine. */
diff -r 67819df51893 -r cfe6395292a7 src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -33,6 +33,8 @@
#include "fortran.def"
#include "CosmologyParameters.h"
+#define DEBUG_PS
+
void mt_init(unsigned_int seed);
unsigned_long_int mt_random();
@@ -41,13 +43,15 @@
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles)
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles)
{
int partnum = 0, i = 0, child = 0, m = 0, numpart = 0, innerchild = 0;
int xindex=0, yindex=0, zindex=0;
float rad = 0.0, sin60 = 0.0;
FLOAT NewPos[3][CHILDRENPERPARENT];
+ FLOAT LeftEdge[3], RightEdge[3];
int ChildrenPerParent = CHILDRENPERPARENT, total_children = 0;
float alpha[3];
float l11 = 0.0, l12 = 0.0, l13 = 0.0;
@@ -67,6 +71,63 @@
//Reproducible random seed
mt_init(((unsigned_int) ParticleSplitterRandomSeed));
+ FLOAT sep[3], midpoint[3], newsep[3];
+
+ /*
+ * The following three options determine over what region
+ * particle splitting takes place.
+ * By default the splitting occurs over the entire most refined
+ * region. At times this can be too large so two further options
+ * exist controlled via the parameter file. The first just selects
+ * a fraction of that region centred on the centre of the most
+ * refined region.
+ * The second allows the user to specify the centre of the splitting
+ * and the size of the region.
+ */
+
+ /* Default Option */
+ for (i = 0; i < 3; i++) {
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ /* Change particle refinement region in basic way */
+ if(ParticleSplitterFraction[iter] != 1.0)
+ {
+ fprintf(stderr, "Setting Particle Refinement using fractions.\n");
+ for (i = 0; i < 3; i++) {
+ sep[i] = RefineRegionRightEdge[i] - RefineRegionLeftEdge[i];
+ midpoint[i] = sep[i]/2.0 + RefineRegionLeftEdge[i];
+ newsep[i] = sep[i]*ParticleSplitterFraction[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ }
+ }
+ /* Centre particle refinement region around a specific point */
+ if(ParticleSplitterCenter[0] > 0.0 && ParticleSplitterCenterRegion[iter] > 0.0) {
+ fprintf(stderr, "Setting Particle Refinement around point.\n");
+ for (i = 0; i < 3; i++) {
+ midpoint[i] = ParticleSplitterCenter[i];
+ newsep[i] = ParticleSplitterCenterRegion[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ if(LeftEdge[i] < RefineRegionLeftEdge[i])
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ if(RightEdge[i] > RefineRegionRightEdge[i])
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ }
+
+ #ifdef DEBUG_PS
+ fprintf(stdout, "%s: Iteration %d: midpoint = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, midpoint[0], midpoint[1], midpoint[2]);
+ fprintf(stdout, "%s: Iteration %d: newsep = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, newsep[0], newsep[1], newsep[2]);
+ fprintf(stdout, "%s: Iteration %d: Refine Particles between (%lf, %lf, %lf) and (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, LeftEdge[0], LeftEdge[1], LeftEdge[2],
+ RightEdge[0], RightEdge[1], RightEdge[2]);
+ #endif
+
+
/*
* Loop over existing (parent) particles; It implicitly assumes that
* only DM and conventional star particles get splitted. Other
@@ -84,12 +145,12 @@
* blatantly discriminate against non-refined types around
* here.
*/
- if(ParticlePosition[0][partnum] < RefineRegionLeftEdge[0] ||
- ParticlePosition[0][partnum] > RefineRegionRightEdge[0] ||
- ParticlePosition[1][partnum] < RefineRegionLeftEdge[1] ||
- ParticlePosition[1][partnum] > RefineRegionRightEdge[1] ||
- ParticlePosition[2][partnum] < RefineRegionLeftEdge[2] ||
- ParticlePosition[2][partnum] > RefineRegionRightEdge[2])
+ if(ParticlePosition[0][partnum] < LeftEdge[0] ||
+ ParticlePosition[0][partnum] > RightEdge[0] ||
+ ParticlePosition[1][partnum] < LeftEdge[1] ||
+ ParticlePosition[1][partnum] > RightEdge[1] ||
+ ParticlePosition[2][partnum] < LeftEdge[2] ||
+ ParticlePosition[2][partnum] > RightEdge[2])
{
//fprintf(stdout, "grid::PS: Particle outside RR - ignoring\n");
continue;
@@ -302,8 +363,9 @@
}
}
#ifdef DEBUG_PS
- fprintf(stdout,"%d new child particles created in this grid from %d parents",
- total_children, NumberOfParticles);
+ if(total_children > 0)
+ fprintf(stdout,"Iteration %d: %d new child particles created in this grid from %d parents\n",
+ iter, total_children, NumberOfParticles);
#endif
*NumberOfNewParticles = total_children;
diff -r 67819df51893 -r cfe6395292a7 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -40,22 +40,7 @@
float *VelocityUnits, FLOAT Time);
int FindField(int field, int farray[], int numfields);
-extern "C" void FORTRAN_NAME(particle_splitter)(int *nx, int *ny, int *nz,
- int *idual, int *imetal, hydro_method *imethod, float *dt,
- float *r, float *dx, FLOAT *t, float *z,
- float *d1, float *x1, float *v1, float *t1,
- FLOAT *xstart, FLOAT *ystart, FLOAT *zstart, int *ibuff,
- int *npart,
- FLOAT *xpold, FLOAT *ypold, FLOAT *zpold, float *upold, float *vpold, float *wpold,
- float *mpold, float *tdpold, float *tcpold, float *metalfold, int *typeold,
- int *nmax, int *npartnew, int *children, int *level,
- FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
- float *mp, float *tdp, float *tcp, float *metalf, int *type,
- int *iterations, float *separation, int *ran1_init,
- FLOAT *rr_leftedge, FLOAT *rr_rightedge);
-
-
-int grid::ParticleSplitter(int level)
+int grid::ParticleSplitter(int level, int iteration)
{
if (ParticleSplitterIterations == 0)
@@ -131,7 +116,7 @@
/* Allocate space for new particles. */
- int ChildrenPerParent = 12; //12+1 = 13 will be the final number of particles per parent
+ int ChildrenPerParent = CHILDRENPERPARENT; //12+1 = 13 will be the final number of particles per parent
int MaximumNumberOfNewParticles = ChildrenPerParent*NumberOfParticles+1;
tg->AllocateNewParticles(MaximumNumberOfNewParticles);
@@ -168,7 +153,8 @@
if(!tg->CreateChildParticles(CellWidthTemp, NumberOfParticles, ParticleMass,
ParticleType, ParticlePosition, ParticleVelocity,
ParticleAttribute, CellLeftEdge, GridDimension,
- MaximumNumberOfNewParticles, &NumberOfNewParticles))
+ MaximumNumberOfNewParticles, iteration,
+ &NumberOfNewParticles))
{
fprintf(stdout, "Failed to create child particles in grid %d\n", this->GetGridID());
return FAIL;
diff -r 67819df51893 -r cfe6395292a7 src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -38,7 +38,7 @@
#include "CosmologyParameters.h"
#include "CommunicationUtilities.h"
-#define NO_DEBUG_PS
+#define NO_DEBUG_PS
int RebuildHierarchy(TopGridData *MetaData,
LevelHierarchyEntry *LevelArray[], int level);
@@ -70,6 +70,13 @@
ParticleSplitterChildrenParticleSeparation <=0)
return SUCCESS;
+ if(ParticleSplitterIterations > MAX_SPLIT_ITERATIONS) {
+ fprintf(stderr, "WARNING: Splitting iterations exceeds maximum allowed\n" \
+ "Resetting to %d\n", MAX_SPLIT_ITERATIONS);
+ ParticleSplitterIterations = MAX_SPLIT_ITERATIONS;
+ }
+
+
/* Find total NumberOfParticles in all grids; this is needed in
CommunicationUpdateStarParticleCount below */
@@ -78,9 +85,6 @@
/* Initialize all star particles if this is a restart */
- if (ParticleSplitterIterations > 1)
- fprintf(stderr, "WARNING: ParticleSplitterIterations > 1 is not properly tested yet.\n");
-
for (i = 0; i < ParticleSplitterIterations; i++) {
for (level = 0; level < MAX_DEPTH_OF_HIERARCHY-1; level++) {
@@ -98,7 +102,6 @@
RecordTotalStarParticleCount(Grids, NumberOfGrids,
TotalStarParticleCountPrevious);
-
// for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
// fprintf(stdout, "TotalStarParticleCountPrevious[grid=%d] = %d\n", grid1,
// TotalStarParticleCountPrevious[grid1]);
@@ -110,7 +113,7 @@
Grids[grid1]->GridData->ReturnNumberOfParticles());
#endif
- if (Grids[grid1]->GridData->ParticleSplitter(level) == FAIL) {
+ if (Grids[grid1]->GridData->ParticleSplitter(level, i) == FAIL) {
ENZO_FAIL("Error in grid::ParticleSplitter.\n");
}
diff -r 67819df51893 -r cfe6395292a7 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1264,6 +1264,14 @@
&ParticleSplitterChildrenParticleSeparation);
ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
&ParticleSplitterMustRefine);
+ ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
+ ParticleSplitterFraction+3);
+ ret += sscanf(line, "ParticleSplitterCenter = %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterCenter+0, ParticleSplitterCenter+1, ParticleSplitterCenter+2);
+ ret += sscanf(line, "ParticleSplitterCenterRegion = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterCenterRegion+0, ParticleSplitterCenterRegion+1,
+ ParticleSplitterCenterRegion+2, ParticleSplitterCenterRegion+3);
ret += sscanf(line, "ResetMagneticField = %"ISYM,
&ResetMagneticField);
ret += sscanf(line, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM,
diff -r 67819df51893 -r cfe6395292a7 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -963,6 +963,12 @@
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
ParticleSplitterMustRefine = FALSE;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterFraction[i] = 1.0;
+ for(int i = 0; i < MAX_DIMENSION; i++)
+ ParticleSplitterCenter[i] = -1.0;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterCenterRegion[i] = -1.0;
/* Magnetic Field Resetter */
diff -r 67819df51893 -r cfe6395292a7 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -624,6 +624,12 @@
ParticleSplitterRandomSeed);
fprintf(fptr, "ParticleSplitterChildrenParticleSeparation = %"FSYM"\n",
ParticleSplitterChildrenParticleSeparation);
+ fprintf(fptr, "ParticleSplitterFraction = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterFraction);
+ fprintf(fptr, "ParticleSplitterCenter = ");
+ WriteListOfFloats(fptr, 3, ParticleSplitterCenter);
+ fprintf(fptr, "ParticleSplitterCenterRegion = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterCenterRegion);
fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
ParticleSplitterMustRefine);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
diff -r 67819df51893 -r cfe6395292a7 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -985,6 +985,9 @@
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
EXTERN int ParticleSplitterMustRefine;
+EXTERN float ParticleSplitterFraction[MAX_SPLIT_ITERATIONS];
+EXTERN FLOAT ParticleSplitterCenter[MAX_DIMENSION];
+EXTERN float ParticleSplitterCenterRegion[MAX_SPLIT_ITERATIONS];
/* Magnetic Field Resetter */
diff -r 67819df51893 -r cfe6395292a7 src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -478,6 +478,11 @@
#define PARTICLE_TYPE_SIMPLE_SOURCE 10
#define CHILDRENPERPARENT 12
+
+/* You probably don't want to do this too often */
+
+#define MAX_SPLIT_ITERATIONS 4
+
/* Ways to deposit particles from a subgrid. */
#define CIC_DEPOSIT 0
https://bitbucket.org/enzo/enzo-dev/commits/fbcce23071d3/
Changeset: fbcce23071d3
Branch: week-of-code
User: jwise77
Date: 2018-04-16 11:15:16+00:00
Summary: Correcting float macro for particle splitter region center
Affected #: 1 file
diff -r cfe6395292a7 -r fbcce23071d3 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1267,7 +1267,7 @@
ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
ParticleSplitterFraction+3);
- ret += sscanf(line, "ParticleSplitterCenter = %"FSYM" %"FSYM" %"FSYM"",
+ ret += sscanf(line, "ParticleSplitterCenter = %"PSYM" %"PSYM" %"PSYM"",
ParticleSplitterCenter+0, ParticleSplitterCenter+1, ParticleSplitterCenter+2);
ret += sscanf(line, "ParticleSplitterCenterRegion = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
ParticleSplitterCenterRegion+0, ParticleSplitterCenterRegion+1,
https://bitbucket.org/enzo/enzo-dev/commits/ff55ab979b6b/
Changeset: ff55ab979b6b
Branch: week-of-code
User: jwise77
Date: 2018-04-18 15:22:30+00:00
Summary: Additional MRP check to not include DM particles that were created
from old star particles and merging routines.
Affected #: 1 file
diff -r fbcce23071d3 -r ff55ab979b6b src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -92,13 +92,23 @@
// Flag particles as must refine particles
int *IsParticleMustRefine;
IsParticleMustRefine = new int[NumberOfParticles];
+ bool OriginalParticle;
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
+ if (NumberOfParticleAttributes > 0)
+ OriginalParticle = (ParticleAttribute[0][i] < 0);
+ else
+ OriginalParticle = true;
- // check particle type and uniform mass
+ // check particle type and uniform mass. Also check particle
+ // creation time for DM particles that are positive, indicating
+ // that they are either inert stellar remnants or "leftovers" from
+ // particle merging
+
rules[0] = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
ParticleType[i] == PARTICLE_TYPE_MBH) ||
- (ParticleMass[i] < UniformParticleMass);
+ (ParticleMass[i] < UniformParticleMass &&
+ ParticleType[i] == PARTICLE_TYPE_DARK_MATTER && OriginalParticle);
// check particle mass greater than minimum mass
rules[1] = (ParticleMass[i] > MustRefineParticlesMinimumMass);
@@ -110,6 +120,7 @@
// set flag for this particle
for (int j = 0; j < NumberOfRules; j++)
IsParticleMustRefine[i] *= rules[j];
+
}
PFORTRAN_NAME(cic_flag)(IsParticleMustRefine,
https://bitbucket.org/enzo/enzo-dev/commits/2d5a79c5d3a5/
Changeset: 2d5a79c5d3a5
Branch: week-of-code
User: jwise77
Date: 2018-04-19 18:31:32+00:00
Summary: Adding capability to read particle IDs to mark as must-refine during
particle splitting.
Affected #: 8 files
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2409,7 +2409,8 @@
/* Particle splitter routine. */
- int ParticleSplitter(int level, int iter);
+ int ParticleSplitter(int level, int iter, int NumberOfIDs,
+ long *MustRefineIDs);
int CreateChildParticles(float dx, int NumberOfParticles, float *ParticleMass,
int *ParticleType, FLOAT *ParticlePosition[],
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -134,10 +134,14 @@
* particles - which usually become Star class particles - doesn't
* seem to have any reason to be splitted. (as of Oct.2009)
*/
+
+ /* (Apr 2018) Include must-refine particles */
for(partnum = 0; partnum < NumberOfParticles; partnum++)
{
- if(ParticleMass[partnum] > 0.0 && ParticleType[partnum] <= 2)
+ if(ParticleMass[partnum] > 0.0 &&
+ (ParticleType[partnum] <= 2 ||
+ ParticleType[partnum] == PARTICLE_TYPE_MUST_REFINE))
{
/*
* Check that particle is within the most refined region.
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -40,7 +40,8 @@
float *VelocityUnits, FLOAT Time);
int FindField(int field, int farray[], int numfields);
-int grid::ParticleSplitter(int level, int iteration)
+int grid::ParticleSplitter(int level, int iteration, int NumberOfIDs,
+ long *MustRefineIDs)
{
if (ParticleSplitterIterations == 0)
@@ -130,7 +131,7 @@
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles before splitting = %d, MyProcessorNumber = %d\n",
NumberOfParticles, MyProcessorNumber);
#endif
-
+
if (NumberOfParticles > 0) {
#define NO_PARTICLE_IN_GRID_CHECK
@@ -150,6 +151,24 @@
xindex, yindex, zindex, level);
}
#endif
+
+ /* Before splitting, set particle type to must-refine if given in
+ the provided list. Only do this for the first iteration because
+ in the later iterations, the particle types are propagated to
+ the split particles. */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs != NULL &&
+ iteration == 0) {
+ for (i = 0; i < NumberOfParticles; i++) {
+ for (j = 0; j < NumberOfIDs; j++) {
+ if (this->ParticleNumber[i] == MustRefineIDs[j]) {
+ this->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ break;
+ }
+ } // ENDFOR j
+ } // ENDFOR i
+ }
+
if(!tg->CreateChildParticles(CellWidthTemp, NumberOfParticles, ParticleMass,
ParticleType, ParticlePosition, ParticleVelocity,
ParticleAttribute, CellLeftEdge, GridDimension,
@@ -171,9 +190,10 @@
NumberOfNewParticles, MyProcessorNumber);
#endif
- /* If specified, set particle type to must-refine */
-
- if (ParticleSplitterMustRefine == TRUE) {
+ /* If specified and no ID list given, set all particle types to
+ must-refine */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs == NULL) {
for (i = 0; i < NumberOfNewParticles; i++)
tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
}
@@ -204,7 +224,7 @@
#ifdef DEBUG_PS
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles(New) = %d, MyProcessorNumber = %d\n",
- NumberOfParticles, MyProcessorNumber);
+ NumberOfParticles, MyProcessorNumber);
#endif
} // end: if (NumberOfNewParticles > 0)
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -24,6 +24,7 @@
#include <stdlib.h>
#include <stdio.h>
+#include <hdf5.h>
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "typedefs.h"
@@ -83,6 +84,34 @@
MetaData->NumberOfParticles = FindTotalNumberOfParticles(LevelArray);
NumberOfOtherParticles = MetaData->NumberOfParticles - NumberOfStarParticles;
+ /* Read must-refine particle positions if specified */
+
+ long *MustRefineIDs = NULL;
+ int NumberOfIDs = 0;
+
+ if (ParticleSplitterMustRefine &&
+ ParticleSplitterMustRefineIDFile != NULL) {
+
+ int NumberOfMustRefineIDs;
+ hid_t file_id, dataset_id, dataspace_id;
+ herr_t status;
+ hsize_t dims[1], maxdims[1];
+
+ file_id = H5Fopen(ParticleSplitterMustRefineIDFile,
+ H5F_ACC_RDONLY, H5P_DEFAULT);
+ dataset_id = H5Dopen2(file_id, "/particle_identifier", H5P_DEFAULT);
+ dataspace_id = H5Dget_space(dataset_id);
+ // Get number of particles and allocate variable
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, maxdims);
+ MustRefineIDs = new long[dims[0]];
+ NumberOfIDs = dims[0];
+ status = H5Dread(dataset_id, H5T_NATIVE_LLONG, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ MustRefineIDs);
+ status = H5Dclose(dataset_id);
+ status = H5Fclose(file_id);
+
+ }
+
/* Initialize all star particles if this is a restart */
for (i = 0; i < ParticleSplitterIterations; i++) {
@@ -113,7 +142,8 @@
Grids[grid1]->GridData->ReturnNumberOfParticles());
#endif
- if (Grids[grid1]->GridData->ParticleSplitter(level, i) == FAIL) {
+ if (Grids[grid1]->GridData->ParticleSplitter(level, i, NumberOfIDs,
+ MustRefineIDs) == FAIL) {
ENZO_FAIL("Error in grid::ParticleSplitter.\n");
}
@@ -151,6 +181,7 @@
/* Set splitter parameter zero; otherwise it will split particles again at next restart */
+ delete [] MustRefineIDs;
ParticleSplitterIterations = 0;
return SUCCESS;
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1264,6 +1264,8 @@
&ParticleSplitterChildrenParticleSeparation);
ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
&ParticleSplitterMustRefine);
+ if (sscanf(line, "ParticleSplitterMustRefineIDFile = %s", dummy) == 1)
+ ParticleSplitterMustRefineIDFile = dummy;
ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
ParticleSplitterFraction+3);
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -963,6 +963,7 @@
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
ParticleSplitterMustRefine = FALSE;
+ ParticleSplitterMustRefineIDFile = NULL;
for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
ParticleSplitterFraction[i] = 1.0;
for(int i = 0; i < MAX_DIMENSION; i++)
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -632,6 +632,8 @@
WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterCenterRegion);
fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
ParticleSplitterMustRefine);
+ fprintf(fptr, "ParticleSplitterMustRefineIDFile = %s\n",
+ ParticleSplitterMustRefineIDFile);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
ResetMagneticField);
fprintf(fptr, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM"\n",
diff -r ff55ab979b6b -r 2d5a79c5d3a5 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -985,6 +985,7 @@
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
EXTERN int ParticleSplitterMustRefine;
+EXTERN char *ParticleSplitterMustRefineIDFile;
EXTERN float ParticleSplitterFraction[MAX_SPLIT_ITERATIONS];
EXTERN FLOAT ParticleSplitterCenter[MAX_DIMENSION];
EXTERN float ParticleSplitterCenterRegion[MAX_SPLIT_ITERATIONS];
https://bitbucket.org/enzo/enzo-dev/commits/5409f66b0ab2/
Changeset: 5409f66b0ab2
Branch: week-of-code
User: jwise77
Date: 2018-04-20 14:24:19+00:00
Summary: Don't rebuild hierarchy before particle splitting. Will lose information if using must-refine particles.
Affected #: 2 files
diff -r 2d5a79c5d3a5 -r 5409f66b0ab2 src/enzo/EvolveHierarchy.C
--- a/src/enzo/EvolveHierarchy.C
+++ b/src/enzo/EvolveHierarchy.C
@@ -270,20 +270,26 @@
*/
+ /* Particle Splitter. Split particles into 13 (=1+12) child
+ particles. The hierarchy is rebuilt inside this routine. */
+
+ if (MetaData.FirstTimestepAfterRestart == TRUE &&
+ ParticleSplitterIterations > 0) {
+
+ ParticleSplitter(LevelArray, 0, &MetaData);
+
+ } else {
+
/* Do the first grid regeneration. */
- if(CheckpointRestart == FALSE) {
- RebuildHierarchy(&MetaData, LevelArray, 0);
- }
+ if(CheckpointRestart == FALSE) {
+ RebuildHierarchy(&MetaData, LevelArray, 0);
+ }
+
+ } // ENDELSE particle splitting
PrintMemoryUsage("1st rebuild");
- /* Particle Splitter. Split particles into 13 (=1+12) child particles */
-
- if (MetaData.FirstTimestepAfterRestart == TRUE &&
- ParticleSplitterIterations > 0)
- ParticleSplitter(LevelArray, 0, &MetaData);
-
/* Reset magnetic fields if requested. */
if (MetaData.FirstTimestepAfterRestart == TRUE &&
diff -r 2d5a79c5d3a5 -r 5409f66b0ab2 src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -61,15 +61,13 @@
HierarchyEntry **Grids;
int NumberOfGrids;
- /* First, rebuild the hierarchy */
-
- RebuildHierarchy(MetaData, LevelArray, 0);
-
- /* Return if this does not concern us */
+ /* Rebuild hierarchy and return if this does not concern us */
if (ParticleSplitterIterations <= 0 ||
- ParticleSplitterChildrenParticleSeparation <=0)
+ ParticleSplitterChildrenParticleSeparation <=0) {
+ RebuildHierarchy(MetaData, LevelArray, 0);
return SUCCESS;
+ }
if(ParticleSplitterIterations > MAX_SPLIT_ITERATIONS) {
fprintf(stderr, "WARNING: Splitting iterations exceeds maximum allowed\n" \
https://bitbucket.org/enzo/enzo-dev/commits/828e88f05f0f/
Changeset: 828e88f05f0f
Branch: week-of-code
User: jwise77
Date: 2018-04-23 14:18:14+00:00
Summary: Additional restrictions on split must-refine particles
Affected #: 2 files
diff -r 5409f66b0ab2 -r 828e88f05f0f src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -93,10 +93,11 @@
int *IsParticleMustRefine;
IsParticleMustRefine = new int[NumberOfParticles];
bool OriginalParticle;
+ const float OriginalCrit = FLOAT_UNDEFINED+1;
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
if (NumberOfParticleAttributes > 0)
- OriginalParticle = (ParticleAttribute[0][i] < 0);
+ OriginalParticle = (ParticleAttribute[0][i] < OriginalCrit);
else
OriginalParticle = true;
diff -r 5409f66b0ab2 -r 828e88f05f0f src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -197,6 +197,12 @@
for (i = 0; i < NumberOfNewParticles; i++)
tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
}
+
+ if (NumberOfParticleAttributes > 0)
+ for (i = 0; i < NumberOfNewParticles; i++) {
+ if (tg->ParticleType[i] == PARTICLE_TYPE_DARK_MATTER)
+ tg->ParticleAttribute[0][i] = -1.0;
+ }
/* If not set in the above routine, then set the metal fraction to zero. */
https://bitbucket.org/enzo/enzo-dev/commits/68d91a7109f9/
Changeset: 68d91a7109f9
Branch: week-of-code
User: jwise77
Date: 2018-05-04 15:44:42+00:00
Summary: More debugging of split must-refine particles
Affected #: 3 files
diff -r 828e88f05f0f -r 68d91a7109f9 src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -47,7 +47,7 @@
int *NumberOfNewParticles)
{
- int partnum = 0, i = 0, child = 0, m = 0, numpart = 0, innerchild = 0;
+ int partnum = 0, i = 0, child = 0, m = 0, innerchild = 0;
int xindex=0, yindex=0, zindex=0;
float rad = 0.0, sin60 = 0.0;
FLOAT NewPos[3][CHILDRENPERPARENT];
@@ -350,8 +350,13 @@
for(i = 0; i < 3; i++)
this->ParticleVelocity[i][child] = ParticleVelocity[i][partnum];
this->ParticleType[child] = ParticleType[partnum];
+ // Flag that a DM particle was split (originally -99999 or 0)
+ if (NumberOfParticleAttributes > 0 &&
+ ParticleType[partnum] == PARTICLE_TYPE_DARK_MATTER)
+ if (ParticleAttribute[0][partnum] <= 0)
+ ParticleAttribute[0][partnum] = tiny_number;
for(i = 0; i < NumberOfParticleAttributes; i++)
- this->ParticleAttribute[i][child] = ParticleAttribute[i][numpart];
+ this->ParticleAttribute[i][child] = ParticleAttribute[i][partnum];
}
/* Loop forward by CHILDRENPERPARENT each time. */
diff -r 828e88f05f0f -r 68d91a7109f9 src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -93,11 +93,10 @@
int *IsParticleMustRefine;
IsParticleMustRefine = new int[NumberOfParticles];
bool OriginalParticle;
- const float OriginalCrit = FLOAT_UNDEFINED+1;
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
if (NumberOfParticleAttributes > 0)
- OriginalParticle = (ParticleAttribute[0][i] < OriginalCrit);
+ OriginalParticle = (ParticleAttribute[0][i] <= 0.0);
else
OriginalParticle = true;
@@ -115,7 +114,6 @@
rules[1] = (ParticleMass[i] > MustRefineParticlesMinimumMass);
// add more rules here
-
//
// set flag for this particle
diff -r 828e88f05f0f -r 68d91a7109f9 src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -198,12 +198,6 @@
tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
}
- if (NumberOfParticleAttributes > 0)
- for (i = 0; i < NumberOfNewParticles; i++) {
- if (tg->ParticleType[i] == PARTICLE_TYPE_DARK_MATTER)
- tg->ParticleAttribute[0][i] = -1.0;
- }
-
/* If not set in the above routine, then set the metal fraction to zero. */
if (MetallicityField == FALSE && NumberOfParticleAttributes > 1)
https://bitbucket.org/enzo/enzo-dev/commits/3e4499464d7b/
Changeset: 3e4499464d7b
Branch: week-of-code
User: jwise77
Date: 2018-10-09 15:43:08+00:00
Summary: Adding must-refine particle splitting parameters to docs
Affected #: 1 file
diff -r 68d91a7109f9 -r 3e4499464d7b doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1865,10 +1865,44 @@
Currently it implicitly assumes that only DM (type=1) and
conventional star particles (type=2) inside the ``RefineRegion`` get
split. Other particles, which usually become Star class objects,
- seem to have no reason to be split. Default: 0
+ seem to have no reason to be split. It is reset to zero after a
+ restart to avoid resplitting in subsequent datasets. It can be set
+ to a maximum of 4. Default: 0
+``ParticleSplitterRandomSeed`` (external)
+ Random seed used when randomly rotating the hexagonal close
+ packed array on whose vertices the split particles are
+ placed. Default: 131180
+``ParticleSplitterMustRefine`` (external)
+ Set to 1 to mark the split particles as must-refine particles. The
+ user must also set associated must-refine parameters to enable its
+ machinery. Default: 0
+``ParticleSplitterMustRefineIDFile`` (external)
+ Filename for the HDF5 file that has a dataset containing the
+ particle IDs that should be marked as must-refine. All other
+ particles within the region marked for splitting will retain their
+ original types. This must be used in conjunction with
+ ``ParticleSplitterMustRefine = 1``. The dataset must be named
+ ``particle_identifier`` in the base group. Default: (null)
+``ParticleSplitterFraction`` (external)
+ An array of four values that represent the width of the splitting
+ region in units of the original refine region set by
+ ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. Each
+ successive value represents the next nested split region. Valid
+ up to ``ParticleSplitterIterations`` times. Default: 1.0 (all 4
+ values)
+``ParticleSplitterCenter`` (external)
+ The center of split region in code units. Specify if the split
+ region does not correspond to the center of the refine region.
+ Not used if negative. Default: -1.0 -1.0 -1.0
+``ParticleSplitterCenterRegion`` (external)
+ The width of the split region in code units. Must be used in
+ conjunction with ``ParticleSplitterCenter``. Each successive
+ value represents the next nested split region. Valid up to
+ ``ParticleSplitterIterations`` times. Not used if
+ negative. Default: -1.0 (all 4 values)
``ParticleSplitterChildrenParticleSeparation`` (external)
This is the spacing between the child particles placed on a
- hexagonal close-packed (HCP) array. In the unit of a cell size
+ hexagonal close-packed (HCP) array. In units of a cell size
which the parent particle resides in. Default: 1.0
.. _starparticleparameters:
https://bitbucket.org/enzo/enzo-dev/commits/40e4a6cb3d8a/
Changeset: 40e4a6cb3d8a
Branch: week-of-code
User: jwise77
Date: 2018-10-16 15:49:24+00:00
Summary: Adding an example yt script to generate particle splitting ID HDF5 file
Affected #: 1 file
diff -r 3e4499464d7b -r 40e4a6cb3d8a doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1882,7 +1882,22 @@
particles within the region marked for splitting will retain their
original types. This must be used in conjunction with
``ParticleSplitterMustRefine = 1``. The dataset must be named
- ``particle_identifier`` in the base group. Default: (null)
+ ``particle_identifier`` in the base group. An example yt script is
+ provided below, selecting the particles in a sphere centered at
+ [0.5, 0.5, 0.5] with a radius 0.05 in code length units, using yt::
+
+ import yt
+ import h5py as h5
+
+ ds = yt.load('DD0040/DD0040')
+ center = ds.arr([0.5, 0.5, 0.5], 'code_length')
+ radius = ds.quan(0.05, 'code_length')
+ sp = ds.sphere(center, radius)
+ fp = h5.File('particle-ids.h5', 'w')
+ fp['particle_identifier'] = sp['particle_index'].astype('int')
+ fp.close()
+
+ Default: (null)
``ParticleSplitterFraction`` (external)
An array of four values that represent the width of the splitting
region in units of the original refine region set by
https://bitbucket.org/enzo/enzo-dev/commits/93c2de54c3fc/
Changeset: 93c2de54c3fc
Branch: week-of-code
User: jwise77
Date: 2018-10-16 16:29:06+00:00
Summary: Correct syntax highlighting for sample code
Affected #: 1 file
diff -r 40e4a6cb3d8a -r 93c2de54c3fc doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1882,9 +1882,13 @@
particles within the region marked for splitting will retain their
original types. This must be used in conjunction with
``ParticleSplitterMustRefine = 1``. The dataset must be named
- ``particle_identifier`` in the base group. An example yt script is
- provided below, selecting the particles in a sphere centered at
- [0.5, 0.5, 0.5] with a radius 0.05 in code length units, using yt::
+ ``particle_identifier`` in the base group. Default: (null).
+
+ An example yt script is provided below, selecting the particles in
+ a sphere centered at [0.5, 0.5, 0.5] with a radius 0.05 in code
+ length units.
+
+ .. code-block:: python
import yt
import h5py as h5
@@ -1897,7 +1901,6 @@
fp['particle_identifier'] = sp['particle_index'].astype('int')
fp.close()
- Default: (null)
``ParticleSplitterFraction`` (external)
An array of four values that represent the width of the splitting
region in units of the original refine region set by
https://bitbucket.org/enzo/enzo-dev/commits/ff98226315b7/
Changeset: ff98226315b7
Branch: week-of-code
User: jwise77
Date: 2019-01-09 05:01:55+00:00
Summary: Merging with latest mainline updates
Affected #: 39 files
diff -r 93c2de54c3fc -r ff98226315b7 doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -904,7 +904,21 @@
constrained by static grids. Instead, refinement around dark matter
particles is allowed down to the level of a particle's generation level.
Refinement beyond this level is allowed around particles within the MUSIC
- ellipsoidal making region.
+ ellipsoidal masking region. Note, dark matter particles corresponding to
+ a generation level N are guaranteed to be refined to at least level N, but
+ may also exist on levels N > 1 if in the vicinity of an N > 1 dark matter
+ particle or a tagged must-refine particle.
+ 4. Similar to setting 3, except dark matter particles corresponding to a
+ generation level N are refined only to level N and no further. If two
+ dark matter particles from different levels occupy the same cell, that
+ cell will be refined to the coarser level. Tagged must-refine particles
+ near coarse dark matter particles will be similarly de-refined. Compared
+ to option 3, this can be used to prevent unnecessary mesh refinement in
+ regions where coarser particles enter into higher resolution regions,
+ slowing down the simulation. Note, with this setting, coarse boundary
+ particles entering into a high resolution region will eventually lead
+ to total derefinement of the region of interest.
+
``MustRefineParticlesRefineToLevel`` (external)
The maximum level on which ``MustRefineParticles`` are required to
diff -r 93c2de54c3fc -r ff98226315b7 doc/manual/source/physics/star_particles.rst
--- a/doc/manual/source/physics/star_particles.rst
+++ b/doc/manual/source/physics/star_particles.rst
@@ -552,29 +552,36 @@
----------------------------
*Source: hydro_rk/SuperNovaSeedField.C*
-Select this method by setting ``UseSupernovaSeedFieldSourceTerms = 1``
-(Default = 0) and
-specifying the following parameters:
+Select this method by setting ``UseMagneticSupernovaFeedback = 1``
+(Default = 0) and specifying the parameters below. If
+``UseMagneticSupernovaFeedback == 2``, the parameters
+``MagneticSupernovaRadius`` and ``MagneticSupernovaDuration`` will
+be calculated to be the minimum allowed values (see below)
+at runtime based on the cell width and time step of the most-refined grid.
-``SupernovaSeedFieldTotalEnergy`` (in units of ergs) is the total amount
-of magnetic energy to be injected by a single supernova event. Defualt = 0.0.
-``SupernovaSeedFieldRadius`` (in units of parsecs) gives the scale over
+``MagneticSupernovaEnergy`` (in units of ergs) is the total amount
+of magnetic energy to be injected by a single supernova event. Defualt = 1e51
+ergs.
+
+``MagneticSupernovaRadius`` (in units of parsecs) gives the scale over
which to inject supernova energy. The injection mechanism normalizes the
spatial exponential decay of the injected supernova energy so that all of the
energy is contained within the specified radius. For this reason, the
-``SupernovaSeedFieldRadius`` should be at least 3 times the minimum cell width of
-the simulation. Default = 0.0.
+``MagneticSupernovaRadius`` should be at least 1.5 times the minimum cell width of
+the simulation (in pc). Default = 300 pc.
-``SupernovaSeedFieldDuration`` (in units of Myr) gives the duration of the
+``MagneticSupernovaDuration`` (in units of years) gives the duration of the
supernova magnetic energy injection. The injection mechanism is normalized so
-that all of the ``SupernovaSeedFieldTotalEnergy`` is injected over this
-time scale. In order to inject the correct amount of energy, ``SupernovaSeedFieldDuration`` should be set to at least 4
-times the minimum time step of the simulation. Default = 0.0.
+that all of the ``MagneticSupernovaEnergy`` is injected over this
+time scale. In order to inject the correct amount of energy, ``MagneticSupernovaDuration`` should be set to at least 5
+times the minimum time step of the simulation. Default = 50,000 years.
-The following applies to Methods 0 (Cen & Ostriker) and 1 (+
-stochastic star formation). The magnetic feedback method is described fully in `Butsky et al. (2017)
+The following applies to all star formation methods that produce a
+PARTICLE_TYPE_STAR object. Methods 0 (Cen & Ostriker) and 1 (+
+stochastic star formation) have been tested extensively.
+The magnetic feedback method is described fully in `Butsky et al. (2017)
<http://adsabs.harvard.edu/abs/2017ApJ...843..113B>`_.
When a star cluster particle reaches the end of its lifetime, we inject a
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/CosmologyReadParameters.C
--- a/src/enzo/CosmologyReadParameters.C
+++ b/src/enzo/CosmologyReadParameters.C
@@ -45,6 +45,7 @@
FinalRedshift = 0;
CosmologyTableNumberOfBins = 1000;
CosmologyTableLogt = NULL;
+ CosmologyTableLoga = NULL;
CosmologyTableLogaInitial = -6.0;
CosmologyTableLogaFinal = 0.0;
CosmologyTableLogtIndex = 0;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/EnzoTiming.h
--- a/src/enzo/EnzoTiming.h
+++ b/src/enzo/EnzoTiming.h
@@ -29,8 +29,15 @@
#include <cstring>
#include <map>
-#define min(A,B) ((A) < (B) ? (A) : (B))
-#define max(A,B) ((A) > (B) ? (A) : (B))
+template <typename T>
+bool min(const T& A, const T& B) {
+ return A < B ? A : B;
+}
+
+template <typename T>
+bool max(const T& A, const T& B) {
+ return A > B ? A : B;
+}
double ReturnWallTime(void);
void Reduce_Times(double time, double *time_array);
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/EvolveLevel.C
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -540,6 +540,10 @@
Grids[grid1]->GridData->CopyBaryonFieldToOldBaryonField();
+ // Find recently-supernova stars to add them the MagneticSupernovaList
+ if ((UseMagneticSupernovaFeedback) && (level == MaximumRefinementLevel))
+ Grids[grid1]->GridData->AddMagneticSupernovaeToList();
+
/* Call hydro solver and save fluxes around subgrids.
* HD_RK and MHD_RK are the 2nd order Runge-Kutta integrations, which
* require two steps (*_1stStep and *_2ndStep)
@@ -547,7 +551,7 @@
* All others (PPM, Zeus, MHD_Li/CT) are called from SolveHydroEquations
*/
-
+
if( HydroMethod != HD_RK && HydroMethod != MHD_RK ){
Grids[grid1]->GridData->SolveHydroEquations(LevelCycleCount[level],
NumberOfSubgrids[grid1], SubgridFluxesEstimate[grid1], level);
@@ -713,7 +717,10 @@
if (ComovingCoordinates)
Grids[grid1]->GridData->ComovingExpansionTerms();
- } // end loop over grids
+ if (UseMagneticSupernovaFeedback)
+ Grids[grid1]->GridData->MagneticSupernovaList.clear();
+
+ } // end loop over grids
/* Finalize (accretion, feedback, etc.) star particles */
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FinalizeFluxes.C
--- a/src/enzo/FinalizeFluxes.C
+++ b/src/enzo/FinalizeFluxes.C
@@ -48,10 +48,11 @@
DeleteFluxes(SubgridFluxesEstimate[grid1][subgrid]);
delete SubgridFluxesEstimate[grid1][subgrid];
}
- delete [] SubgridFluxesEstimate[grid1];
-
}
-
+
+ delete [] SubgridFluxesEstimate[grid1];
+ SubgridFluxesEstimate[grid1] = NULL;
+
} // end of loop over grids
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FindSubgrids.C
--- a/src/enzo/FindSubgrids.C
+++ b/src/enzo/FindSubgrids.C
@@ -85,7 +85,8 @@
}
/* Set the static (permanent) regions. */
- if (MustRefineParticlesCreateParticles != 3){
+ if (MustRefineParticlesCreateParticles != 3 &&
+ MustRefineParticlesCreateParticles != 4) {
if (CurrentGrid->SetFlaggingFieldStaticRegions(level, NumberOfFlaggedCells)
== FAIL) {
ENZO_FAIL("Error in grid->SetFlaggingFieldStaticRegions.");
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FreeRealMem.C
--- /dev/null
+++ b/src/enzo/FreeRealMem.C
@@ -0,0 +1,199 @@
+#ifdef TASKMAP
+
+#ifdef USE_LL
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <sys/types.h>
+#include "llapi.h"
+
+#include "macros_and_parameters.h"
+
+
+Eint64 FreeRealMem( char *node )
+{
+
+ LL_element *queryObject;
+ LL_element *job, *step;
+
+ Eint32 rc, obj_count, err_code, ok;
+
+ Eint64 free;
+ Eint64 freemem;
+ Eint64 def_mem;
+
+ def_mem = (Eint64)(12.0); /* minimum memory on DataStar */
+
+ Eint32 node_number;
+
+#ifdef SEABORG
+ sscanf(node, "s%5d", &node_number);
+ fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
+
+ if ( 701 <= node_number && node_number <= 713 ) {
+ def_mem = (Eint64)(56.0);
+ } else {
+
+ if ( 3701 <= node_number && node_number <= 5213 ) {
+ def_mem = (Eint64)(28.0);
+ } else {
+ def_mem = (Eint64)(14.0);
+ }
+
+ }
+#endif
+
+#ifdef DATASTAR
+ sscanf(node, "ds%3d", &node_number);
+ fprintf(stderr, "DataStar Node %s number %d\n", node, node_number);
+
+ if ( 4 <= node_number && node_number <= 11 ) {
+ def_mem = (Eint64)(120.0);
+ } else {
+
+ if ( 300 <= node_number && node_number <= 395 ) {
+ def_mem = (Eint64)(27.0);
+ } else {
+ def_mem = (Eint64)(12.0);
+ }
+
+ }
+#endif
+
+ def_mem = def_mem * 1024.0;
+ fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
+
+ ok = 0;
+
+ /* Initialize the query */
+
+ queryObject = ll_query(MACHINES);
+
+ if (!queryObject) {
+ fprintf(stderr, "Query MACHINES: ll_query() returns NULL.\n");
+ freemem = def_mem;
+ return(freemem);
+ }
+
+ char **nodelist;
+ nodelist = (char **) malloc(2*sizeof(char *));
+ char nodename[6];
+ strcpy(nodename, node);
+
+ nodelist[0] = nodename;
+ nodelist[1] = NULL;
+
+ rc = ll_set_request(queryObject, QUERY_HOST, nodelist, ALL_DATA);
+
+ if (rc) {
+ fprintf(stderr, "Query MACHINES: ll_set_request() return code is non-zero.\n");
+ freemem = def_mem;
+ return(freemem);
+ }
+
+ job = ll_get_objs(queryObject, LL_CM, NULL, &obj_count, &err_code);
+
+ /* Process the job objects */
+
+ rc = ll_get_data(job, LL_MachineFreeRealMemory64, &free);
+
+ if (rc) {
+ fprintf(stderr, "Get data: ll_get_data() return code is non-zero: %d\n", rc);
+ freemem = def_mem;
+ return(freemem);
+ }
+
+/*
+ fprintf(stderr, "Free: %lld MBytes\n", free);
+*/
+
+ ll_free_objs(queryObject);
+ ll_deallocate(queryObject);
+
+ freemem = (Eint64)(free);
+
+ return(freemem);
+}
+
+#else
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "macros_and_parameters.h"
+
+Eint64 FreeRealMem( char *node )
+{
+
+ Eint64 freemem;
+ Eint64 def_mem;
+
+ def_mem = (Eint64)(12.5); /* minimum memory of DataStar */
+
+ Eint32 node_number;
+
+
+#ifdef SEABORG
+ sscanf(node, "s%5d", &node_number);
+ fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
+
+ if ( 701 <= node_number && node_number <= 713 ) {
+ def_mem = (Eint64)(56.0);
+ } else {
+
+ if ( 3701 <= node_number && node_number <= 5213 ) {
+ def_mem = (Eint64)(28.0);
+ } else {
+ def_mem = (Eint64)(14.0);
+ }
+
+ }
+#endif
+
+#ifdef DATASTAR
+ sscanf(node, "ds%3d", &node_number);
+ fprintf(stderr, "DataStarNode %s number %d\n", node, node_number);
+
+ if ( 4 <= node_number && node_number <= 11 ) {
+ def_mem = (Eint64)(120.0);
+ } else {
+
+ if ( 300 <= node_number && node_number <= 395 ) {
+ def_mem = (Eint64)(27.0);
+ } else {
+ def_mem = (Eint64)(12.0);
+ }
+
+ }
+#endif
+
+ def_mem = def_mem * 1024.0;
+ fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
+
+ freemem = def_mem;
+
+ return(freemem);
+}
+
+#endif
+
+#else
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "macros_and_parameters.h"
+
+Eint64 FreeRealMem( char *node )
+{
+
+ Eint64 freemem;
+
+ fprintf(stderr, "Error in FreeRealMem - the TASKMAP macro is not defined! Default to 12.0 GBytes.\n");
+
+ freemem = (Eint64)(12.0*1024.0); /* minimum memory on DataStar */
+
+ return(freemem);
+}
+
+#endif
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/FreeRealMem.c
--- a/src/enzo/FreeRealMem.c
+++ /dev/null
@@ -1,199 +0,0 @@
-#ifdef TASKMAP
-
-#ifdef USE_LL
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include "llapi.h"
-
-#include "macros_and_parameters.h"
-
-
-Eint64 FreeRealMem( char *node )
-{
-
- LL_element *queryObject;
- LL_element *job, *step;
-
- Eint32 rc, obj_count, err_code, ok;
-
- Eint64 free;
- Eint64 freemem;
- Eint64 def_mem;
-
- def_mem = (Eint64)(12.0); /* minimum memory on DataStar */
-
- Eint32 node_number;
-
-#ifdef SEABORG
- sscanf(node, "s%5d", &node_number);
- fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
-
- if ( 701 <= node_number && node_number <= 713 ) {
- def_mem = (Eint64)(56.0);
- } else {
-
- if ( 3701 <= node_number && node_number <= 5213 ) {
- def_mem = (Eint64)(28.0);
- } else {
- def_mem = (Eint64)(14.0);
- }
-
- }
-#endif
-
-#ifdef DATASTAR
- sscanf(node, "ds%3d", &node_number);
- fprintf(stderr, "DataStar Node %s number %d\n", node, node_number);
-
- if ( 4 <= node_number && node_number <= 11 ) {
- def_mem = (Eint64)(120.0);
- } else {
-
- if ( 300 <= node_number && node_number <= 395 ) {
- def_mem = (Eint64)(27.0);
- } else {
- def_mem = (Eint64)(12.0);
- }
-
- }
-#endif
-
- def_mem = def_mem * 1024.0;
- fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
-
- ok = 0;
-
- /* Initialize the query */
-
- queryObject = ll_query(MACHINES);
-
- if (!queryObject) {
- fprintf(stderr, "Query MACHINES: ll_query() returns NULL.\n");
- freemem = def_mem;
- return(freemem);
- }
-
- char **nodelist;
- nodelist = (char **) malloc(2*sizeof(char *));
- char nodename[6];
- strcpy(nodename, node);
-
- nodelist[0] = nodename;
- nodelist[1] = NULL;
-
- rc = ll_set_request(queryObject, QUERY_HOST, nodelist, ALL_DATA);
-
- if (rc) {
- fprintf(stderr, "Query MACHINES: ll_set_request() return code is non-zero.\n");
- freemem = def_mem;
- return(freemem);
- }
-
- job = ll_get_objs(queryObject, LL_CM, NULL, &obj_count, &err_code);
-
- /* Process the job objects */
-
- rc = ll_get_data(job, LL_MachineFreeRealMemory64, &free);
-
- if (rc) {
- fprintf(stderr, "Get data: ll_get_data() return code is non-zero: %d\n", rc);
- freemem = def_mem;
- return(freemem);
- }
-
-/*
- fprintf(stderr, "Free: %lld MBytes\n", free);
-*/
-
- ll_free_objs(queryObject);
- ll_deallocate(queryObject);
-
- freemem = (Eint64)(free);
-
- return(freemem);
-}
-
-#else
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "macros_and_parameters.h"
-
-Eint64 FreeRealMem( char *node )
-{
-
- Eint64 freemem;
- Eint64 def_mem;
-
- def_mem = (Eint64)(12.5); /* minimum memory of DataStar */
-
- Eint32 node_number;
-
-
-#ifdef SEABORG
- sscanf(node, "s%5d", &node_number);
- fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
-
- if ( 701 <= node_number && node_number <= 713 ) {
- def_mem = (Eint64)(56.0);
- } else {
-
- if ( 3701 <= node_number && node_number <= 5213 ) {
- def_mem = (Eint64)(28.0);
- } else {
- def_mem = (Eint64)(14.0);
- }
-
- }
-#endif
-
-#ifdef DATASTAR
- sscanf(node, "ds%3d", &node_number);
- fprintf(stderr, "DataStarNode %s number %d\n", node, node_number);
-
- if ( 4 <= node_number && node_number <= 11 ) {
- def_mem = (Eint64)(120.0);
- } else {
-
- if ( 300 <= node_number && node_number <= 395 ) {
- def_mem = (Eint64)(27.0);
- } else {
- def_mem = (Eint64)(12.0);
- }
-
- }
-#endif
-
- def_mem = def_mem * 1024.0;
- fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
-
- freemem = def_mem;
-
- return(freemem);
-}
-
-#endif
-
-#else
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "macros_and_parameters.h"
-
-Eint64 FreeRealMem( char *node )
-{
-
- Eint64 freemem;
-
- fprintf(stderr, "Error in FreeRealMem - the TASKMAP macro is not defined! Default to 12.0 GBytes.\n");
-
- freemem = (Eint64)(12.0*1024.0); /* minimum memory on DataStar */
-
- return(freemem);
-}
-
-#endif
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -11,6 +11,7 @@
************************************************************************/
#ifndef GRID_DEFINED__
#define GRID_DEFINED__
+#include <vector>
#include "ProtoSubgrid.h"
#include "ListOfParticles.h"
#include "region.h"
@@ -20,7 +21,6 @@
#include "Star.h"
#include "FOF_allvars.h"
#include "MemoryPool.h"
-#include "list.h"
#include "hydro_rk/SuperNova.h"
#ifdef ECUDA
#include "hydro_rk/CudaMHD.h"
@@ -3013,11 +3013,10 @@
int MHDCT_ConvertEnergyToConservedS();
int MHDCT_ConvertEnergyToSpecificS();
- //List of SuperNova objects that each grid needs to keep track of
-
- List<SuperNova> SuperNovaList;
-
-
+ // used if UseMagneticSupernovaFeedback is turned on
+ int AddMagneticSupernovaeToList();
+ //List of SuperNova objects that each grid needs to keep track of
+ std::vector<SuperNova> MagneticSupernovaList;
};
// inline int grid::ReadRandomForcingFields (FILE *main_file_pointer);
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_AddFeedbackSphere.C
--- a/src/enzo/Grid_AddFeedbackSphere.C
+++ b/src/enzo/Grid_AddFeedbackSphere.C
@@ -43,8 +43,6 @@
int FindField(int field, int farray[], int numfields);
-void mt_init(unsigned_int seed);
-
int grid::AddFeedbackSphere(Star *cstar, int level, float radius, float DensityUnits,
float LengthUnits, float VelocityUnits,
float TemperatureUnits, float TimeUnits, double EjectaDensity,
@@ -990,48 +988,7 @@
/* Now it's done, unmark. */
//cstar->FeedbackFlag = NO_FEEDBACK;
-
-
- // Create a randomly-oriented SuperNova object and add it to the SuperNova Grid list
- if (cstar->FeedbackFlag == SUPERNOVA_SEEDFIELD) {
- if(UseSupernovaSeedFieldSourceTerms){
-
- SuperNova *P = new SuperNova();
- int mt_seed = 100;
- mt_init(mt_seed);
- float random_u = (float)mt_seed/100.0; // random variable from 0 to 1
- float random_v = (float)mt_seed/100.0;
- float random_phi = 2*M_PI*random_u; // 0 to 2pi
- float random_theta = acos(2*random_v-1); // 0 to pi
- // Setting up randomly oriented magnetic feedback of supernova
- float phi_x = sin(random_theta)*cos(random_phi);
- float phi_y = sin(random_theta)*sin(random_phi);
- float phi_z = cos(random_theta);
-
- // Birthtime of supernova is at the end of a star particle's life
- float sn_birthtime = cstar->BirthTime + cstar->LifeTime;
-
- // Convert units to system units
- // Converting time from years to seconds, then internal units
- float sn_duration = SupernovaSeedFieldDuration * 3.1556952e7 / TimeUnits;
- // Converting radius from parsecs to cm, then internal units
- float sn_radius = SupernovaSeedFieldRadius * 3.0856775714e18 / LengthUnits;
- // Converting energy from ergs to internal units
- float MassUnits = DensityUnits * POW(LengthUnits, 3);
- float sn_energy = SupernovaSeedFieldEnergy / (MassUnits*VelocityUnits*VelocityUnits);
-
- // Creates a supernova with magnetic feedback set by user-defined parameters and
- // adds it to supernova list
- if((Time > sn_birthtime) && (Time < sn_birthtime + SupernovaSeedFieldDuration)){
- P->setValues(phi_x, phi_y, phi_z, cstar->pos[0], cstar->pos[1], cstar->pos[2],
- sn_radius, sn_birthtime, sn_duration, sn_energy);
-
- this->SuperNovaList.append(P);
- printf("Added supernova to list \n");
- }
- }
- }
-
+
return SUCCESS;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_AddMagneticSupernovaeToList.C
--- /dev/null
+++ b/src/enzo/Grid_AddMagneticSupernovaeToList.C
@@ -0,0 +1,112 @@
+/***********************************************************************
+/
+/ FIND DEAD STAR PARTICLES AND ADD THEM TO MAGNETIC SUPERNOVA LIST
+/
+/ written by: Iryna Butsky
+/ date: June 2018
+/
+/ PURPOSE: To apply magnetic supernova feedback, we need to first add
+/ all stars that have recently gone supernova to the global list.
+/
+************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "phys_constants.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "TopGridData.h"
+#include "Grid.h"
+#include "hydro_rk/SuperNova.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+ float *TemperatureUnits, float *TimeUnits,
+ float *VelocityUnits, double *MassUnits, FLOAT Time);
+
+void mt_init(unsigned_int seed);
+unsigned_long_int mt_random();
+
+int grid::AddMagneticSupernovaeToList()
+{
+
+ if (ProcessorNumber != MyProcessorNumber) {
+ return SUCCESS;
+ }
+
+ float random_u, random_v, random_phi, random_theta, phi_x, phi_y, phi_z;
+ float sn_birthtime, sn_duration, sn_radius, sn_energy, star_birthtime, star_lifetime;
+
+ float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, VelocityUnits, MassUnits, EnergyUnits;
+ if (GetUnits(&DensityUnits, &LengthUnits,&TemperatureUnits, &TimeUnits,
+ &VelocityUnits, &MassUnits, Time) == FAIL){
+ fprintf(stderr, "Error in GetUnits.\n");
+ return FAIL;
+ }
+
+
+ // if UseMagneticSupernovaFeedback > 1, then we set the magnetic feedback radius and duration
+ // below based on the resolution of the grid at the highest refinement level
+ if (UseMagneticSupernovaFeedback > 1){
+ sn_duration = 5.0 * dtFixed;
+ MagneticSupernovaDuration = sn_duration * TimeUnits / SecondsPerYear;
+ sn_radius = 1.5 * this->CellWidth[0][0];
+ MagneticSupernovaRadius = sn_radius * LengthUnits / pc;
+ }
+ else {
+ // Converting time from years to seconds, then internal units
+ sn_duration = MagneticSupernovaDuration * SecondsPerYear / TimeUnits;
+ if(sn_duration < 5.0 * dtFixed) {
+ printf("WARNING: Magnetic supernova feedback duration is less than the recommended minimum of 5 x dtFixed\n");
+ printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / SecondsPerYear);
+ }
+ // Converting radius from parsecs to cm, then internal units
+ sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
+ if(sn_radius < 1.5 * this->CellWidth[0][0]){
+ printf("WARNING: Magnetic supernova feedback radius is less than the recommended minimum of 1.5 x CellWidth\n");
+ printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc);
+ }
+
+ }
+ // Converting energy from ergs to internal units
+ sn_energy = MagneticSupernovaEnergy / (MassUnits*VelocityUnits*VelocityUnits);
+
+ for (int star_id = 0; star_id < this->NumberOfParticles; star_id++){
+ if (this->ParticleType[star_id] == PARTICLE_TYPE_STAR){
+
+ star_birthtime = this->ParticleAttribute[0][star_id];
+ star_lifetime = this->ParticleAttribute[1][star_id];
+ sn_birthtime = star_birthtime + star_lifetime;
+
+ // Creates a supernova with magnetic feedback set by user-defined parameters and
+ // adds it to supernova list
+ if((this->Time > sn_birthtime) && (this->Time < sn_birthtime + sn_duration)){
+ SuperNova P = SuperNova();
+ mt_init((unsigned_int) this->ParticleNumber[star_id]);
+
+ random_u = (float)(mt_random()%32768)/32768.0; // random variable from 0 to 1
+ random_v = (float)(mt_random()%32768)/32768.0;
+ random_phi = 2*M_PI*random_u; // 0 to 2pi
+ random_theta = acos(2*random_v-1); // 0 to pi
+
+ // Setting up randomly oriented magnetic feedback of supernova
+ phi_x = sin(random_theta)*cos(random_phi);
+ phi_y = sin(random_theta)*sin(random_phi);
+ phi_z = cos(random_theta);
+
+ P.setValues(phi_x, phi_y, phi_z,
+ ParticlePosition[0][star_id], ParticlePosition[1][star_id], ParticlePosition[2][star_id],
+ sn_radius, sn_birthtime, sn_duration, sn_energy);
+ this->MagneticSupernovaList.push_back(P);
+ }
+ }
+ }
+ return SUCCESS;
+}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -80,19 +80,42 @@
FlaggingField[i] = 0;
float UniformParticleMass = 0.0;
- if (ProblemType == 30 && MustRefineParticlesCreateParticles == 3)
+ if (ProblemType == 30 &&
+ (MustRefineParticlesCreateParticles == 3 ||
+ MustRefineParticlesCreateParticles == 4))
UniformParticleMass = OmegaDarkMatterNow / OmegaMatterNow;
/* Loop over all particles, marking wich ones are must refine
To add rules, modify number of rules here and add to loop below */
- bool *rules;
+ bool *rules, rule0;
const int NumberOfRules = 2;
rules = new bool[NumberOfRules];
+ // Rules to prevent refinement, cancelling out the above rules.
+ bool *antirules;
+ int *AntiFlaggingField;
+ int NumberOfAntiRules = 0;
+ antirules = new bool[NumberOfAntiRules];
+
+ // Add an antirule to unflag over-refined dark matter particles.
+ if (MustRefineParticlesCreateParticles == 4) {
+ NumberOfAntiRules++;
+ }
+
+ if (NumberOfAntiRules > 0) {
+ antirules = new bool[NumberOfAntiRules];
+ AntiFlaggingField = new int[size];
+ for (i = 0; i < size; i++)
+ AntiFlaggingField[i] = 0;
+ }
+
// Flag particles as must refine particles
- int *IsParticleMustRefine;
+ int *IsParticleMustRefine, *IsParticleNotMustRefine;
+ bool OriginalParticle;
IsParticleMustRefine = new int[NumberOfParticles];
- bool OriginalParticle;
+ if (NumberOfAntiRules > 0) {
+ IsParticleNotMustRefine = new int[NumberOfParticles];
+ }
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
if (NumberOfParticleAttributes > 0)
@@ -105,21 +128,34 @@
// that they are either inert stellar remnants or "leftovers" from
// particle merging
- rules[0] = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
+ // check particle type and uniform mass
+ rule0 = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
ParticleType[i] == PARTICLE_TYPE_MBH) ||
(ParticleMass[i] < UniformParticleMass &&
ParticleType[i] == PARTICLE_TYPE_DARK_MATTER && OriginalParticle);
+ rules[0] = rule0;
// check particle mass greater than minimum mass
rules[1] = (ParticleMass[i] > MustRefineParticlesMinimumMass);
// add more rules here
+
//
// set flag for this particle
for (int j = 0; j < NumberOfRules; j++)
IsParticleMustRefine[i] *= rules[j];
+ // anti-rules
+ if (NumberOfAntiRules > 0) {
+ IsParticleNotMustRefine[i] = 1;
+ // check for over-refined dark matter particles
+ antirules[0] = !rule0;
+ }
+
+ // set antiflag for this particle
+ for (int j = 0; j < NumberOfAntiRules; j++)
+ IsParticleNotMustRefine[i] *= antirules[j];
}
PFORTRAN_NAME(cic_flag)(IsParticleMustRefine,
@@ -128,6 +164,18 @@
LeftEdge, GridDimension, GridDimension+1, GridDimension+2,
&CellSize, &ParticleBufferSize);
+ if (NumberOfAntiRules > 0) {
+ PFORTRAN_NAME(cic_flag)(IsParticleNotMustRefine,
+ ParticlePosition[0], ParticlePosition[1], ParticlePosition[2],
+ &GridRank, &NumberOfParticles, AntiFlaggingField,
+ LeftEdge, GridDimension, GridDimension+1, GridDimension+2,
+ &CellSize, &ParticleBufferSize);
+
+ for (i = 0;i < size;i++) {
+ FlaggingField[i] *= !(AntiFlaggingField[i]);
+ }
+ }
+
/* Increase particle mass flagging field for definite refinement */
float MustRefineMass =
@@ -141,7 +189,9 @@
particles, and don't use the particle mass field. */
int NumberOfFlaggedCells = 0;
- if (!(ProblemType == 30 && MustRefineParticlesCreateParticles == 3 &&
+ if (!(ProblemType == 30 &&
+ (MustRefineParticlesCreateParticles == 3 ||
+ MustRefineParticlesCreateParticles == 4) &&
level == MustRefineParticlesRefineToLevel)) {
for (i = 0; i < size; i++)
if (FlaggingField[i]) {
@@ -170,6 +220,12 @@
delete [] IsParticleMustRefine;
delete [] rules;
+ if (NumberOfAntiRules > 0) {
+ delete [] AntiFlaggingField;
+ delete [] IsParticleNotMustRefine;
+ delete [] antirules;
+ }
+
return NumberOfFlaggedCells;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_FindAllStarParticles.C
--- a/src/enzo/Grid_FindAllStarParticles.C
+++ b/src/enzo/Grid_FindAllStarParticles.C
@@ -50,7 +50,6 @@
//StarType = abs(ParticleType[i]);
if (ParticleType[i] == PARTICLE_TYPE_SINGLE_STAR ||
ParticleType[i] == PARTICLE_TYPE_BLACK_HOLE ||
- (ParticleType[i] == PARTICLE_TYPE_STAR && UseSupernovaSeedFieldSourceTerms) ||
ParticleType[i] == PARTICLE_TYPE_CLUSTER ||
ParticleType[i] == PARTICLE_TYPE_COLOR_STAR ||
ParticleType[i] == PARTICLE_TYPE_MBH ||
@@ -58,7 +57,7 @@
(StarParticleRadiativeFeedback == TRUE &&
ParticleType[i] == PARTICLE_TYPE_STAR)) {
- if (ParticleType[i] == PARTICLE_TYPE_STAR)
+ if (ParticleType[i] == PARTICLE_TYPE_STAR && UseMagneticSupernovaFeedback == 0)
LifetimeFactor = 12.0;
else
LifetimeFactor = 1.0;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_FindNewStarParticles.C
--- a/src/enzo/Grid_FindNewStarParticles.C
+++ b/src/enzo/Grid_FindNewStarParticles.C
@@ -39,9 +39,6 @@
for (i = 0; i < NumberOfParticles; i++)
if (ParticleType[i] == -PARTICLE_TYPE_SINGLE_STAR ||
- (ParticleType[i] == -PARTICLE_TYPE_STAR && UseSupernovaSeedFieldSourceTerms &&
- (this->Time >= ParticleAttribute[0][i] &&
- this->Time <= ParticleAttribute[0][i]+ParticleAttribute[1][i])) ||
ParticleType[i] == -PARTICLE_TYPE_BLACK_HOLE ||
ParticleType[i] == -PARTICLE_TYPE_CLUSTER ||
ParticleType[i] == -PARTICLE_TYPE_COLOR_STAR ||
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_WalkPhotonPackage.C
--- a/src/enzo/Grid_WalkPhotonPackage.C
+++ b/src/enzo/Grid_WalkPhotonPackage.C
@@ -245,15 +245,15 @@
// solid angle associated with package (= 4 Pi/N_package[on this level])
uint64_t Hlevel = (*PP)->level, res = 2L, BRP = 12L;
- uint64_t Nlevel = (BRP * (res << (res*Hlevel-1)));
+ uint64_t Nlevel = BRP * (1L << (res*Hlevel));
double n_on_this_level = Nlevel;
double omega_package=4*M_PI/(n_on_this_level);
double dtheta = sqrt(omega_package);
if(n_on_this_level <= 0.0 || omega_package == INFINITY) {
fprintf(stdout, "%s: level = %llu\n", __FUNCTION__, level);
fprintf(stdout, "%s: 2*(*PP)->level-1 = %llu\n", __FUNCTION__, res*(*PP)->level-1);
- fprintf(stdout, "%s: 12 * (2 << (2*(*PP)->level-1)) = %llu\n", __FUNCTION__,
- BRP * (res << (res*(*PP)->level-1)));
+ fprintf(stdout, "%s: 12 * (1 << (2*(*PP)->level)) = %llu\n", __FUNCTION__,
+ BRP * (1L << (res*(*PP)->level)));
fprintf(stdout, "%s: Hlevel = %llu\n", __FUNCTION__, Hlevel);
fprintf(stdout, "%s: Nlevel = %llu\n", __FUNCTION__, Nlevel);
fprintf(stdout, "%s: n_on_this_level = %lf\n", __FUNCTION__, n_on_this_level);
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Grid_constructor.C
--- a/src/enzo/Grid_constructor.C
+++ b/src/enzo/Grid_constructor.C
@@ -200,12 +200,11 @@
MagneticDims[field][dim] = -100;
}}
-
/* For once-per-rootgrid-timestep star formation, the following flag
determines whether SF is about to occur or not. It's currently
(April 2012) only implemented for H2REG_STAR and completely
ignored for all other star makers. */
MakeStars = 0;
- if (UseSupernovaSeedFieldSourceTerms == 1) List<SuperNova> SuperNovaList;
+ if (UseMagneticSupernovaFeedback) std::vector<SuperNova> MagneticSupernovaList;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Make.config.objects
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -307,6 +307,7 @@
Grid_AddFeedbackSphere.o \
Grid_AddFieldMassToMassFlaggingField.o \
Grid_AddFields.o \
+ Grid_AddMagneticSupernovaeToList.o \
Grid_AddOneParticleFromList.o \
Grid_AddOverlappingParticleMassField.o \
Grid_AddParticlesFromList.o \
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Make.mach.ncsa-bluewaters-gnu
--- a/src/enzo/Make.mach.ncsa-bluewaters-gnu
+++ b/src/enzo/Make.mach.ncsa-bluewaters-gnu
@@ -69,7 +69,7 @@
MACH_CXXFLAGS =
MACH_FFLAGS = -fno-second-underscore -m64
MACH_F90FLAGS = -fno-second-underscore -m64
-MACH_LDFLAGS = -Bdynamic
+MACH_LDFLAGS = -Bstatic
#-----------------------------------------------------------------------
# Optimization flags
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/PhotonGrid_Methods.h
--- a/src/enzo/PhotonGrid_Methods.h
+++ b/src/enzo/PhotonGrid_Methods.h
@@ -102,7 +102,7 @@
int RadiativeTransferXRays(PhotonPackageEntry **PP, FLOAT *dPi, int cellindex,
int species, FLOAT ddr, FLOAT tau, float geo_correction,
- float photonrate, float *excessrate, float *ion_factor2,
+ FLOAT photonrate, FLOAT *excessrate, float *ion_factor2,
float heat_factor, const int *kphNum, int gammaNum);
int RadiativeTransferComptonHeating(PhotonPackageEntry **PP, FLOAT *dPi, int cellindex,
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/RadiativeTransferInitialize.C
--- a/src/enzo/RadiativeTransferInitialize.C
+++ b/src/enzo/RadiativeTransferInitialize.C
@@ -78,7 +78,7 @@
/* Check for radiation fields and delete them */
- NumberOfObsoleteFields = 10;
+ NumberOfObsoleteFields = 13;
ObsoleteFields[0] = kphHI;
ObsoleteFields[1] = PhotoGamma;
ObsoleteFields[2] = kphHeI;
@@ -89,6 +89,9 @@
ObsoleteFields[7] = kphHM;
ObsoleteFields[8] = kdissH2II;
ObsoleteFields[9] = RaySegments;
+ ObsoleteFields[10] = RadPressure0;
+ ObsoleteFields[11] = RadPressure1;
+ ObsoleteFields[12] = RadPressure2;
for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel)
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/RadiativeTransferXRays.C
--- a/src/enzo/RadiativeTransferXRays.C
+++ b/src/enzo/RadiativeTransferXRays.C
@@ -25,7 +25,7 @@
#include "Grid.h"
#include "CosmologyParameters.h"
-int grid::RadiativeTransferXRays(PhotonPackageEntry **PP, float *dPXray, int cellindex,
+int grid::RadiativeTransferXRays(PhotonPackageEntry **PP, FLOAT *dPXray, int cellindex,
int species, FLOAT ddr, FLOAT tau, float geo_correction,
FLOAT photonrate, FLOAT *excessrate, float *ion2_factor,
float heat_factor, const int *kphNum, int gammaNum)
@@ -57,7 +57,7 @@
#define COMPTON 3
-int grid::RadiativeTransferComptonHeating(PhotonPackageEntry **PP, float *dPXray, int cellindex,
+int grid::RadiativeTransferComptonHeating(PhotonPackageEntry **PP, FLOAT *dPXray, int cellindex,
float LengthUnits, float photonrate,
int TemperatureField, FLOAT ddr, double dN,
float geo_correction, int gammaNum)
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1285,10 +1285,10 @@
ret += sscanf(line, "GasDragCoefficient = %"GSYM, &GasDragCoefficient);
// Parameters for magnetic feedback from supernovae
- ret += sscanf(line, "UseSupernovaSeedFieldSourceTerms = %"ISYM, &UseSupernovaSeedFieldSourceTerms);
- ret += sscanf(line,"SupernovaSeedFieldRadius = %"FSYM, &SupernovaSeedFieldRadius);
- ret += sscanf(line,"SupernovaSeedFieldEnergy = %"FSYM, &SupernovaSeedFieldEnergy);
- ret += sscanf(line,"SupernovaSeedFieldDuration = %"FSYM, &SupernovaSeedFieldDuration);
+ ret += sscanf(line, "UseMagneticSupernovaFeedback = %"ISYM, &UseMagneticSupernovaFeedback);
+ ret += sscanf(line,"MagneticSupernovaRadius = %"FSYM, &MagneticSupernovaRadius);
+ ret += sscanf(line,"MagneticSupernovaEnergy = %"FSYM, &MagneticSupernovaEnergy);
+ ret += sscanf(line,"MagneticSupernovaDuration = %"FSYM, &MagneticSupernovaDuration);
/* If the dummy char space was used, then make another. */
@@ -1657,6 +1657,14 @@
if (initialize_chemistry_data(&grackle_units) == FAIL) {
ENZO_FAIL("Error in Grackle initialize_chemistry_data.\n");
}
+
+ // Need to set these after initialize_chemistry_data since
+ // that function sets them automatically based on the tables.
+ if (FinalRedshift < grackle_data->UVbackground_redshift_off) {
+ grackle_data->UVbackground_redshift_off = FinalRedshift;
+ grackle_data->UVbackground_redshift_drop = FinalRedshift;
+ }
+
} // if (grackle_data->use_grackle == TRUE)
else {
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/RestartPhotons.C
--- a/src/enzo/RestartPhotons.C
+++ b/src/enzo/RestartPhotons.C
@@ -147,6 +147,9 @@
/* Optically thin Lyman-Werner (H2) radiation field */
+ if (AllStars == NULL)
+ return SUCCESS;
+
int NumberOfSources = 0;
Star *cstar = AllStars->NextStar;
while (cstar != NULL) {
@@ -157,12 +160,10 @@
if (RadiativeTransferOpticallyThinH2)
for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel) {
- if (Temp->GridData->InitializeTemperatureFieldForH2Shield() == FAIL) {
- ENZO_FAIL("Error in InitializeTemperatureFieldForH2Shield.\n");
- }
+ Temp->GridData->InitializeTemperatureFieldForH2Shield();
Temp->GridData->AddH2Dissociation(AllStars, NumberOfSources);
}
-
+
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -986,12 +986,12 @@
UseGasDrag = 0;
GasDragCoefficient = 0.;
- /* Supernova magnetic seed field */
+ /* Supernova magnetic feedback */
/* Default == 0 -> no magnetic field contribution */
- UseSupernovaSeedFieldSourceTerms = 0;
- SupernovaSeedFieldRadius = 0.0;
- SupernovaSeedFieldDuration = 0.0;
- SupernovaSeedFieldEnergy = 0.0;
+ UseMagneticSupernovaFeedback = 0; // 0 = off; 1+ = on; 2 = Radius and duration calculated during runtime
+ MagneticSupernovaRadius = 300.0; // Total injection radius of magnetic field in parsecs
+ MagneticSupernovaDuration = 5e4 ; // Total duration of magnetic feedback in years
+ MagneticSupernovaEnergy = 1.0e51; // Total energy (ergs) injected per star particle (supernova)
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/StarParticleAddFeedback.C
--- a/src/enzo/StarParticleAddFeedback.C
+++ b/src/enzo/StarParticleAddFeedback.C
@@ -88,10 +88,7 @@
&TimeUnits, &VelocityUnits, Time);
count = 0;
- // clear list of Supernovae at each timestep to avoid adding duplicates in Grid_AddFeedbackSphere
- if(UseSupernovaSeedFieldSourceTerms){
- LevelArray[level]->GridData->SuperNovaList.clear();
- }
+
for (cstar = AllStars; cstar; cstar = cstar->NextStar, count++) {
AddedFeedback[count] = false;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Star_Accrete.C
--- a/src/enzo/Star_Accrete.C
+++ b/src/enzo/Star_Accrete.C
@@ -29,15 +29,10 @@
int Star::Accrete(void)
{
- if (UseSupernovaSeedFieldSourceTerms == 1){
- if (this->CurrentGrid == NULL ||(this->naccretions == 0))
- return SUCCESS;
- }
- else {
- if (this->CurrentGrid == NULL ||
- (this->naccretions == 0 && fabs(this->DeltaMass) < tiny_number))
- return SUCCESS;
- }
+ if (this->CurrentGrid == NULL ||
+ (this->naccretions == 0 && fabs(this->DeltaMass) < tiny_number))
+ return SUCCESS;
+
const double Msun = 1.989e33, yr = 3.1557e7;
int dim, i, n, count;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Star_HitEndpoint.C
--- a/src/enzo/Star_HitEndpoint.C
+++ b/src/enzo/Star_HitEndpoint.C
@@ -101,9 +101,6 @@
case PopIII_CF:
break;
- case NormalStar:
- break;
-
} // ENDSWITCH
return result;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/Star_SetFeedbackFlag.C
--- a/src/enzo/Star_SetFeedbackFlag.C
+++ b/src/enzo/Star_SetFeedbackFlag.C
@@ -135,12 +135,7 @@
//this->type = abs_type;
- case PARTICLE_TYPE_STAR:
- if(UseSupernovaSeedFieldSourceTerms)
- this->FeedbackFlag = SUPERNOVA_SEEDFIELD;
- break;
-}
-
+ }
return SUCCESS;
}
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -1210,10 +1210,10 @@
fprintf(fptr, "CorrectParentBoundaryFlux = %d\n", CorrectParentBoundaryFlux);
/* Supernova magnetic seed field */
- fprintf(fptr, "UseSupernovaSeedFieldSourceTerms = %"ISYM"\n", UseSupernovaSeedFieldSourceTerms);
- fprintf(fptr, "SupernovaSeedFieldRadius = %"GSYM"\n", SupernovaSeedFieldRadius);
- fprintf(fptr,"SupernovaSeedFieldEnergy = %"GSYM"\n",SupernovaSeedFieldEnergy);
- fprintf(fptr,"SupernovaSeedFieldDuration = %"GSYM"\n",SupernovaSeedFieldDuration);
+ fprintf(fptr, "UseMagneticSupernovaFeedback = %"ISYM"\n", UseMagneticSupernovaFeedback);
+ fprintf(fptr, "MagneticSupernovaRadius = %"GSYM"\n", MagneticSupernovaRadius);
+ fprintf(fptr,"MagneticSupernovaEnergy = %"GSYM"\n",MagneticSupernovaEnergy);
+ fprintf(fptr,"MagneticSupernovaDuration = %"GSYM"\n",MagneticSupernovaDuration);
/* Output current time */
time_t ID;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -1114,9 +1114,9 @@
EXTERN float GalaxySimulationPreWindVelocity[MAX_DIMENSION];
/* Supernova magnetic seed field */
-EXTERN int UseSupernovaSeedFieldSourceTerms;
-EXTERN float SupernovaSeedFieldRadius;
-EXTERN float SupernovaSeedFieldDuration;
-EXTERN float SupernovaSeedFieldEnergy;
+EXTERN int UseMagneticSupernovaFeedback;
+EXTERN float MagneticSupernovaRadius;
+EXTERN float MagneticSupernovaDuration;
+EXTERN float MagneticSupernovaEnergy;
#endif
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/hydro_rk/Grid_MHDSourceTerms.C
--- a/src/enzo/hydro_rk/Grid_MHDSourceTerms.C
+++ b/src/enzo/hydro_rk/Grid_MHDSourceTerms.C
@@ -19,13 +19,14 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-#include "list.h"
+#include "phys_constants.h"
#include "Fluxes.h"
#include "GridList.h"
#include "ExternalBoundary.h"
#include "TopGridData.h"
#include "Grid.h"
#include "EOS.h"
+#include "hydro_rk/SuperNova.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -33,6 +34,8 @@
int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
int FindField(int field, int farray[], int numfields);
+void mt_init(unsigned_int seed);
+unsigned_long_int mt_random();
int grid::MHDSourceTerms(float **dU)
{
@@ -357,49 +360,73 @@
}
}
-if((UseSupernovaSeedFieldSourceTerms == 1)) {
+ if(UseMagneticSupernovaFeedback) {
+
+ int n, active_x, active_y, center_i, center_j, center_k, num_sn_cells_x, num_sn_cells_y, num_sn_cells_z;
+ snsf_source_terms S;
+ float dx, dy, dz, dist_to_sn, magnetic_energy_density;
+ float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, VelocityUnits;
+
+ if (GetUnits(&DensityUnits, &LengthUnits,&TemperatureUnits, &TimeUnits,
+ &VelocityUnits, Time) == FAIL){
+ fprintf(stderr, "Error in GetUnits.\n");
+ return FAIL;
+ }
+
+ // Converting radius from parsecs to cm, then internal units
+ float sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
- int n = 0, igrid;
- int iden=FindField(Density, FieldType, NumberOfBaryonFields);
- snsf_source_terms S;
- List<SuperNova>::Iterator *P = this->SuperNovaList.begin();
- FLOAT cell_center[3];
- FLOAT dx, dy, dz, dist_to_sn;
- int temp =1;
- int entered = 0;
+ active_x = GridDimension[0] - 2*NumberOfGhostZones;
+ active_y = GridDimension[1] - 2*NumberOfGhostZones;
- 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++, n++) {
- igrid = i+(j+k*GridDimension[1])*GridDimension[0];
+ // the number of additional cells that receive magnetic supernova feedback.
+ // assuming the supernova is in the center of one cell
+ num_sn_cells_x = (int) ((sn_radius - 0.5*CellWidth[0][0])/ CellWidth[0][0]);
+ num_sn_cells_y = (int) ((sn_radius - 0.5*CellWidth[1][0])/ CellWidth[1][0]);
+ num_sn_cells_z = (int) ((sn_radius - 0.5*CellWidth[2][0])/ CellWidth[2][0]);
+
+ for (std::vector<SuperNova>::iterator current_sn = this->MagneticSupernovaList.begin();
+ current_sn != this->MagneticSupernovaList.end(); current_sn++){
- cell_center[0] = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
- cell_center[1] = CellLeftEdge[1][j] + 0.5*CellWidth[1][j];
- cell_center[2] = CellLeftEdge[2][k] + 0.5*CellWidth[2][k];
+ // find index of the cell nearest to the supernova center
+ // assuming that supernova in the center of that cell
+ center_i = (int)((current_sn->getPosition()[0] - GridLeftEdge[0]) / CellWidth[0][0]);
+ center_j = (int)((current_sn->getPosition()[1] - GridLeftEdge[1]) / CellWidth[1][0]);
+ center_k = (int)((current_sn->getPosition()[2] - GridLeftEdge[2]) / CellWidth[2][0]);
- while(P != this->SuperNovaList.end()){
- dx = P->get()->getPosition()[0] - cell_center[0];
- dy = P->get()->getPosition()[1] - cell_center[1];
- dz = P->get()->getPosition()[2] - cell_center[2];
-
- dist_to_sn = sqrt(dx*dx + dy*dy + dz*dz);
- if (dist_to_sn < 1.1*SupernovaSeedFieldRadius){
- S = P->get()->getSourceTerms(dx, dy, dz, Time);
- double rho = BaryonField[DensNum][igrid];
+ for(int k = center_k - num_sn_cells_z; k <= center_k + num_sn_cells_z; k++){
+ for(int j = center_j - num_sn_cells_y; j <= center_j + num_sn_cells_y; j++){
+ for(int i = center_i - num_sn_cells_x; i <= center_i + num_sn_cells_x; i++){
+
+ // only add magnetic feedback on the active grid cells
+ if ((k >= GridStartIndex[2]) && (k <= GridEndIndex[2]) &&
+ (j >= GridStartIndex[1]) && (j <= GridEndIndex[1]) &&
+ (i >= GridStartIndex[0]) && (i <= GridEndIndex[0])){
+
+ dx = CellWidth[0][0] * (float)(i-center_i);
+ dy = CellWidth[1][0] * (float)(j-center_j);
+ dz = CellWidth[2][0] * (float)(k-center_k);
+
+ dist_to_sn = sqrt(dx*dx + dy*dy + dz*dz);
+ S = current_sn->getSourceTerms(dx, dy, dz, Time);
+
+ // solving for index n
+ // analogous to how igrid is calculated, but taking into acount Ghost Zones
+ n = (i - GridStartIndex[0])+((j-GridStartIndex[1]) + (k-GridStartIndex[2])*active_y) * active_x;
+
+ dU[iBx][n] += S.dbx*dtFixed;
+ dU[iBy][n] += S.dby*dtFixed;
+ dU[iBz][n] += S.dbz*dtFixed;
- dU[iBx][n] += S.dbx*dtFixed;
- dU[iBy][n] += S.dby*dtFixed;
- dU[iBz][n] += S.dbz*dtFixed;
- dU[iEtot][n] += S.dUtot*dtFixed;
+ dU[iEtot][n] += S.dUtot * dtFixed;
+
+ }
- }
- P = P->next();
- }// End of SuperNovaList iteration
- } // End of k for-loop
- } // End of j for-loop
- } // End of i for-loop
-
-} // End of UseSuperNovaSeedFieldSourceTerms scope
+ } // End of k for-loop
+ } // End of j for-loop
+ } // End of i for-loop
+ } // End of MagneticSupernovaList loop
+ } // End of UseMagneticSupernovaFeedback scope
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/hydro_rk/LLF_Zero_MHD.C
--- a/src/enzo/hydro_rk/LLF_Zero_MHD.C
+++ b/src/enzo/hydro_rk/LLF_Zero_MHD.C
@@ -34,7 +34,7 @@
// compute priml and primr
int iprim;
- for (int field = 0; field < NEQ_MHD; field++) {
+ for (int field = 0; field < NEQ_MHD-DualEnergyFormalism; field++) {
for (int i = 0; i < ActiveSize+1; i++) {
iprim = i + NumberOfGhostZones - 1;
priml[field][i] = prim[field][iprim];
@@ -53,9 +53,9 @@
iprim = i+NumberOfGhostZones-1;
for (int field = 0; field < NSpecies; field++) {
if (FluxLine[iD][i] >= 0) {
- species[field][i] = prim[field+5][iprim ];
+ species[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism][iprim ];
} else {
- species[field][i] = prim[field+5][iprim+1];
+ species[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism][iprim+1];
}
}
sum = 0;
@@ -76,9 +76,9 @@
iprim = i+NumberOfGhostZones-1;
for (int field = 0; field < NColor; field++) {
if (FluxLine[iD][i] >= 0) {
- colors[field][i] = prim[field+5+NSpecies][iprim ];
+ colors[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism+NSpecies][iprim ];
} else {
- colors[field][i] = prim[field+5+NSpecies][iprim+1];
+ colors[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism+NSpecies][iprim+1];
}
}
for (int field = 0; field < NColor; field++) {
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/hydro_rk/SuperNovaSeedField.C
--- a/src/enzo/hydro_rk/SuperNovaSeedField.C
+++ b/src/enzo/hydro_rk/SuperNovaSeedField.C
@@ -127,14 +127,14 @@
snsf_source_terms SuperNova::getSourceTerms(double dx, double dy, double dz, double enzoTime){
// ------------------------------------------------------------------- //
- // Inputs: (1) distances to current cell, dx, dy, dz //
+ // Inputs: (1) distances to current cell, dx, dy, dz //
// (2) current time of the simulation //
// //
- // Output: (1) A struct of source terms corresponding to the given //
- // supernova //
- // (i) x,y,z components of time derivative of B-field //
- // (ii) time derivative of magnetic energy density //
- // (iii) time derivative of total energy density //
+ // Output: (1) A struct of source terms corresponding to the given //
+ // supernova //
+ // (i) x,y,z components of time derivative of B-field //
+ // (ii) time derivative of magnetic energy density //
+ // (iii) time derivative of total energy density //
// ------------------------------------------------------------------- //
FLOAT phi[3], zhat_cell[3];
@@ -180,6 +180,9 @@
FLOAT r_scale = (r_cyl/characteristicLength)*\
exp(- r_s*r_s /(characteristicLength*characteristicLength));
+
+ FLOAT db_r_scale = (r_cyl/characteristicLength)*\
+ exp(- r_s*r_s /(2.0*characteristicLength*characteristicLength));
FLOAT t_exp = exp(-t_s/characteristicTime);
FLOAT db_t_exp = t_exp / characteristicTime;
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -1,5 +1,12 @@
#ifndef __macros_and_parameters_h_
#define __macros_and_parameters_h_
+#include "string.h"
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <algorithm>
+#include <map>
+#include <cassert>
/***********************************************************************
/
/ MACRO DEFINITIONS AND PARAMETERS
@@ -522,7 +529,6 @@
#define MBH_THERMAL 7
#define MBH_JETS 8
#define COLOR_FIELD 9
-#define SUPERNOVA_SEEDFIELD 11
/* Sink particle accretion modes */
diff -r 93c2de54c3fc -r ff98226315b7 src/enzo/phys_constants.h
--- a/src/enzo/phys_constants.h
+++ b/src/enzo/phys_constants.h
@@ -60,4 +60,7 @@
#define kpc 3.0857e21
#define pc 3.0857e18
+/* Time units */
+#define SecondsPerYear 3.1556952e7
+
#endif
https://bitbucket.org/enzo/enzo-dev/commits/a74f32a0094e/
Changeset: a74f32a0094e
Branch: week-of-code
User: jwise77
Date: 2019-01-11 17:04:53+00:00
Summary: More descriptive comment for MAX_SPLIT_ITERATIONS
Affected #: 1 file
diff -r ff98226315b7 -r a74f32a0094e src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -486,7 +486,8 @@
#define CHILDRENPERPARENT 12
-/* You probably don't want to do this too often */
+/* Splitting particles more than 4 times can induce numerical
+ artifacts */
#define MAX_SPLIT_ITERATIONS 4
https://bitbucket.org/enzo/enzo-dev/commits/bb230291c6f2/
Changeset: bb230291c6f2
Branch: week-of-code
User: jwise77
Date: 2019-02-01 17:29:40+00:00
Summary: Updating particle splitting docs to clarify some details.
Affected #: 1 file
diff -r a74f32a0094e -r bb230291c6f2 doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1889,14 +1889,17 @@
``ParticleSplitterMustRefine`` (external)
Set to 1 to mark the split particles as must-refine particles. The
user must also set associated must-refine parameters to enable its
- machinery. Default: 0
+ machinery that can be used to restrict AMR only to the must-refine
+ particles. Default: 0
``ParticleSplitterMustRefineIDFile`` (external)
Filename for the HDF5 file that has a dataset containing the
particle IDs that should be marked as must-refine. All other
particles within the region marked for splitting will retain their
- original types. This must be used in conjunction with
- ``ParticleSplitterMustRefine = 1``. The dataset must be named
- ``particle_identifier`` in the base group. Default: (null).
+ original types. If not set, all particles within the must-refine
+ region will be must-refine particles. This must be used in
+ conjunction with ``ParticleSplitterMustRefine = 1``. The dataset
+ must be named ``particle_identifier`` in the base group. Default:
+ (null).
An example yt script is provided below, selecting the particles in
a sphere centered at [0.5, 0.5, 0.5] with a radius 0.05 in code
@@ -1918,10 +1921,11 @@
``ParticleSplitterFraction`` (external)
An array of four values that represent the width of the splitting
region in units of the original refine region set by
- ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. Each
+ ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. The
+ splitting region is centered on the refine region center. Each
successive value represents the next nested split region. Valid
- up to ``ParticleSplitterIterations`` times. Default: 1.0 (all 4
- values)
+ up to ``ParticleSplitterIterations`` times. Cannot be used with
+ ``ParticleSplitterCenterRegion``. Default: 1.0 (all 4 values)
``ParticleSplitterCenter`` (external)
The center of split region in code units. Specify if the split
region does not correspond to the center of the refine region.
@@ -1929,13 +1933,14 @@
``ParticleSplitterCenterRegion`` (external)
The width of the split region in code units. Must be used in
conjunction with ``ParticleSplitterCenter``. Each successive
- value represents the next nested split region. Valid up to
+ value represents the next nested split region. Cannot be used
+ with ``ParticleSplitterFraction``. Valid up to
``ParticleSplitterIterations`` times. Not used if
negative. Default: -1.0 (all 4 values)
``ParticleSplitterChildrenParticleSeparation`` (external)
This is the spacing between the child particles placed on a
- hexagonal close-packed (HCP) array. In units of a cell size
- which the parent particle resides in. Default: 1.0
+ hexagonal close-packed (HCP) array. In units of a cell size which
+ the parent particle resides in. Default: 1.0
.. _starparticleparameters:
https://bitbucket.org/enzo/enzo-dev/commits/e58dbfd476cb/
Changeset: e58dbfd476cb
Branch: week-of-code
User: gbryan
Date: 2019-02-04 21:10:50+00:00
Summary: Merged in jwise77/clean-copy (pull request #426)
Updates to particle splitter
Approved-by: John Regan <johnanth...@gmail.com>
Approved-by: Greg Bryan <gbr...@astro.columbia.edu>
Approved-by: Britton Smith <britto...@gmail.com>
Affected #: 12 files
diff -r 603a8a701fe4 -r e58dbfd476cb doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -1879,11 +1879,68 @@
Currently it implicitly assumes that only DM (type=1) and
conventional star particles (type=2) inside the ``RefineRegion`` get
split. Other particles, which usually become Star class objects,
- seem to have no reason to be split. Default: 0
+ seem to have no reason to be split. It is reset to zero after a
+ restart to avoid resplitting in subsequent datasets. It can be set
+ to a maximum of 4. Default: 0
+``ParticleSplitterRandomSeed`` (external)
+ Random seed used when randomly rotating the hexagonal close
+ packed array on whose vertices the split particles are
+ placed. Default: 131180
+``ParticleSplitterMustRefine`` (external)
+ Set to 1 to mark the split particles as must-refine particles. The
+ user must also set associated must-refine parameters to enable its
+ machinery that can be used to restrict AMR only to the must-refine
+ particles. Default: 0
+``ParticleSplitterMustRefineIDFile`` (external)
+ Filename for the HDF5 file that has a dataset containing the
+ particle IDs that should be marked as must-refine. All other
+ particles within the region marked for splitting will retain their
+ original types. If not set, all particles within the must-refine
+ region will be must-refine particles. This must be used in
+ conjunction with ``ParticleSplitterMustRefine = 1``. The dataset
+ must be named ``particle_identifier`` in the base group. Default:
+ (null).
+
+ An example yt script is provided below, selecting the particles in
+ a sphere centered at [0.5, 0.5, 0.5] with a radius 0.05 in code
+ length units.
+
+ .. code-block:: python
+
+ import yt
+ import h5py as h5
+
+ ds = yt.load('DD0040/DD0040')
+ center = ds.arr([0.5, 0.5, 0.5], 'code_length')
+ radius = ds.quan(0.05, 'code_length')
+ sp = ds.sphere(center, radius)
+ fp = h5.File('particle-ids.h5', 'w')
+ fp['particle_identifier'] = sp['particle_index'].astype('int')
+ fp.close()
+
+``ParticleSplitterFraction`` (external)
+ An array of four values that represent the width of the splitting
+ region in units of the original refine region set by
+ ``RefineRegionLeftEdge`` and ``RefineRegionRightEdge``. The
+ splitting region is centered on the refine region center. Each
+ successive value represents the next nested split region. Valid
+ up to ``ParticleSplitterIterations`` times. Cannot be used with
+ ``ParticleSplitterCenterRegion``. Default: 1.0 (all 4 values)
+``ParticleSplitterCenter`` (external)
+ The center of split region in code units. Specify if the split
+ region does not correspond to the center of the refine region.
+ Not used if negative. Default: -1.0 -1.0 -1.0
+``ParticleSplitterCenterRegion`` (external)
+ The width of the split region in code units. Must be used in
+ conjunction with ``ParticleSplitterCenter``. Each successive
+ value represents the next nested split region. Cannot be used
+ with ``ParticleSplitterFraction``. Valid up to
+ ``ParticleSplitterIterations`` times. Not used if
+ negative. Default: -1.0 (all 4 values)
``ParticleSplitterChildrenParticleSeparation`` (external)
This is the spacing between the child particles placed on a
- hexagonal close-packed (HCP) array. In the unit of a cell size
- which the parent particle resides in. Default: 1.0
+ hexagonal close-packed (HCP) array. In units of a cell size which
+ the parent particle resides in. Default: 1.0
.. _starparticleparameters:
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/EvolveHierarchy.C
--- a/src/enzo/EvolveHierarchy.C
+++ b/src/enzo/EvolveHierarchy.C
@@ -270,20 +270,26 @@
*/
+ /* Particle Splitter. Split particles into 13 (=1+12) child
+ particles. The hierarchy is rebuilt inside this routine. */
+
+ if (MetaData.FirstTimestepAfterRestart == TRUE &&
+ ParticleSplitterIterations > 0) {
+
+ ParticleSplitter(LevelArray, 0, &MetaData);
+
+ } else {
+
/* Do the first grid regeneration. */
- if(CheckpointRestart == FALSE) {
- RebuildHierarchy(&MetaData, LevelArray, 0);
- }
+ if(CheckpointRestart == FALSE) {
+ RebuildHierarchy(&MetaData, LevelArray, 0);
+ }
+
+ } // ENDELSE particle splitting
PrintMemoryUsage("1st rebuild");
- /* Particle Splitter. Split particles into 13 (=1+12) child particles */
-
- if (MetaData.FirstTimestepAfterRestart == TRUE &&
- ParticleSplitterIterations > 0)
- ParticleSplitter(LevelArray, 0, &MetaData);
-
/* Reset magnetic fields if requested. */
if (MetaData.FirstTimestepAfterRestart == TRUE &&
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2409,13 +2409,15 @@
/* Particle splitter routine. */
- int ParticleSplitter(int level);
+ int ParticleSplitter(int level, int iter, int NumberOfIDs,
+ long *MustRefineIDs);
int CreateChildParticles(float dx, int NumberOfParticles, float *ParticleMass,
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles);
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles);
/* Magnetic field resetting routine. */
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid_CreateChildParticles.C
--- a/src/enzo/Grid_CreateChildParticles.C
+++ b/src/enzo/Grid_CreateChildParticles.C
@@ -33,6 +33,8 @@
#include "fortran.def"
#include "CosmologyParameters.h"
+#define DEBUG_PS
+
void mt_init(unsigned_int seed);
unsigned_long_int mt_random();
@@ -41,13 +43,15 @@
int *ParticleType, FLOAT *ParticlePosition[],
float *ParticleVelocity[], float *ParticleAttribute[],
FLOAT *CellLeftEdge[], int *GridDimension,
- int MaximumNumberOfNewParticles, int *NumberOfNewParticles)
+ int MaximumNumberOfNewParticles, int iter,
+ int *NumberOfNewParticles)
{
- int partnum = 0, i = 0, child = 0, m = 0, numpart = 0, innerchild = 0;
+ int partnum = 0, i = 0, child = 0, m = 0, innerchild = 0;
int xindex=0, yindex=0, zindex=0;
float rad = 0.0, sin60 = 0.0;
FLOAT NewPos[3][CHILDRENPERPARENT];
+ FLOAT LeftEdge[3], RightEdge[3];
int ChildrenPerParent = CHILDRENPERPARENT, total_children = 0;
float alpha[3];
float l11 = 0.0, l12 = 0.0, l13 = 0.0;
@@ -67,16 +71,77 @@
//Reproducible random seed
mt_init(((unsigned_int) ParticleSplitterRandomSeed));
+ FLOAT sep[3], midpoint[3], newsep[3];
+
+ /*
+ * The following three options determine over what region
+ * particle splitting takes place.
+ * By default the splitting occurs over the entire most refined
+ * region. At times this can be too large so two further options
+ * exist controlled via the parameter file. The first just selects
+ * a fraction of that region centred on the centre of the most
+ * refined region.
+ * The second allows the user to specify the centre of the splitting
+ * and the size of the region.
+ */
+
+ /* Default Option */
+ for (i = 0; i < 3; i++) {
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ /* Change particle refinement region in basic way */
+ if(ParticleSplitterFraction[iter] != 1.0)
+ {
+ fprintf(stderr, "Setting Particle Refinement using fractions.\n");
+ for (i = 0; i < 3; i++) {
+ sep[i] = RefineRegionRightEdge[i] - RefineRegionLeftEdge[i];
+ midpoint[i] = sep[i]/2.0 + RefineRegionLeftEdge[i];
+ newsep[i] = sep[i]*ParticleSplitterFraction[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ }
+ }
+ /* Centre particle refinement region around a specific point */
+ if(ParticleSplitterCenter[0] > 0.0 && ParticleSplitterCenterRegion[iter] > 0.0) {
+ fprintf(stderr, "Setting Particle Refinement around point.\n");
+ for (i = 0; i < 3; i++) {
+ midpoint[i] = ParticleSplitterCenter[i];
+ newsep[i] = ParticleSplitterCenterRegion[iter]/2.0;
+ LeftEdge[i] = midpoint[i] - newsep[i];
+ RightEdge[i] = midpoint[i] + newsep[i];
+ if(LeftEdge[i] < RefineRegionLeftEdge[i])
+ LeftEdge[i] = RefineRegionLeftEdge[i];
+ if(RightEdge[i] > RefineRegionRightEdge[i])
+ RightEdge[i] = RefineRegionRightEdge[i];
+ }
+ }
+
+ #ifdef DEBUG_PS
+ fprintf(stdout, "%s: Iteration %d: midpoint = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, midpoint[0], midpoint[1], midpoint[2]);
+ fprintf(stdout, "%s: Iteration %d: newsep = (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, newsep[0], newsep[1], newsep[2]);
+ fprintf(stdout, "%s: Iteration %d: Refine Particles between (%lf, %lf, %lf) and (%lf, %lf, %lf)\n",
+ __FUNCTION__, iter, LeftEdge[0], LeftEdge[1], LeftEdge[2],
+ RightEdge[0], RightEdge[1], RightEdge[2]);
+ #endif
+
+
/*
* Loop over existing (parent) particles; It implicitly assumes that
* only DM and conventional star particles get splitted. Other
* particles - which usually become Star class particles - doesn't
* seem to have any reason to be splitted. (as of Oct.2009)
*/
+
+ /* (Apr 2018) Include must-refine particles */
for(partnum = 0; partnum < NumberOfParticles; partnum++)
{
- if(ParticleMass[partnum] > 0.0 && ParticleType[partnum] <= 2)
+ if(ParticleMass[partnum] > 0.0 &&
+ (ParticleType[partnum] <= 2 ||
+ ParticleType[partnum] == PARTICLE_TYPE_MUST_REFINE))
{
/*
* Check that particle is within the most refined region.
@@ -84,12 +149,12 @@
* blatantly discriminate against non-refined types around
* here.
*/
- if(ParticlePosition[0][partnum] < RefineRegionLeftEdge[0] ||
- ParticlePosition[0][partnum] > RefineRegionRightEdge[0] ||
- ParticlePosition[1][partnum] < RefineRegionLeftEdge[1] ||
- ParticlePosition[1][partnum] > RefineRegionRightEdge[1] ||
- ParticlePosition[2][partnum] < RefineRegionLeftEdge[2] ||
- ParticlePosition[2][partnum] > RefineRegionRightEdge[2])
+ if(ParticlePosition[0][partnum] < LeftEdge[0] ||
+ ParticlePosition[0][partnum] > RightEdge[0] ||
+ ParticlePosition[1][partnum] < LeftEdge[1] ||
+ ParticlePosition[1][partnum] > RightEdge[1] ||
+ ParticlePosition[2][partnum] < LeftEdge[2] ||
+ ParticlePosition[2][partnum] > RightEdge[2])
{
//fprintf(stdout, "grid::PS: Particle outside RR - ignoring\n");
continue;
@@ -285,8 +350,13 @@
for(i = 0; i < 3; i++)
this->ParticleVelocity[i][child] = ParticleVelocity[i][partnum];
this->ParticleType[child] = ParticleType[partnum];
+ // Flag that a DM particle was split (originally -99999 or 0)
+ if (NumberOfParticleAttributes > 0 &&
+ ParticleType[partnum] == PARTICLE_TYPE_DARK_MATTER)
+ if (ParticleAttribute[0][partnum] <= 0)
+ ParticleAttribute[0][partnum] = tiny_number;
for(i = 0; i < NumberOfParticleAttributes; i++)
- this->ParticleAttribute[i][child] = ParticleAttribute[i][numpart];
+ this->ParticleAttribute[i][child] = ParticleAttribute[i][partnum];
}
/* Loop forward by CHILDRENPERPARENT each time. */
@@ -302,8 +372,9 @@
}
}
#ifdef DEBUG_PS
- fprintf(stdout,"%d new child particles created in this grid from %d parents",
- total_children, NumberOfParticles);
+ if(total_children > 0)
+ fprintf(stdout,"Iteration %d: %d new child particles created in this grid from %d parents\n",
+ iter, total_children, NumberOfParticles);
#endif
*NumberOfNewParticles = total_children;
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid_DepositMustRefineParticles.C
--- a/src/enzo/Grid_DepositMustRefineParticles.C
+++ b/src/enzo/Grid_DepositMustRefineParticles.C
@@ -111,17 +111,28 @@
// Flag particles as must refine particles
int *IsParticleMustRefine, *IsParticleNotMustRefine;
+ bool OriginalParticle;
IsParticleMustRefine = new int[NumberOfParticles];
if (NumberOfAntiRules > 0) {
IsParticleNotMustRefine = new int[NumberOfParticles];
}
for (i = 0; i < NumberOfParticles; i ++){
IsParticleMustRefine[i] = 1;
+ if (NumberOfParticleAttributes > 0)
+ OriginalParticle = (ParticleAttribute[0][i] <= 0.0);
+ else
+ OriginalParticle = true;
+
+ // check particle type and uniform mass. Also check particle
+ // creation time for DM particles that are positive, indicating
+ // that they are either inert stellar remnants or "leftovers" from
+ // particle merging
// check particle type and uniform mass
rule0 = (ParticleType[i] == PARTICLE_TYPE_MUST_REFINE ||
ParticleType[i] == PARTICLE_TYPE_MBH) ||
- (ParticleMass[i] < UniformParticleMass);
+ (ParticleMass[i] < UniformParticleMass &&
+ ParticleType[i] == PARTICLE_TYPE_DARK_MATTER && OriginalParticle);
rules[0] = rule0;
// check particle mass greater than minimum mass
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/Grid_ParticleSplitter.C
--- a/src/enzo/Grid_ParticleSplitter.C
+++ b/src/enzo/Grid_ParticleSplitter.C
@@ -40,22 +40,8 @@
float *VelocityUnits, FLOAT Time);
int FindField(int field, int farray[], int numfields);
-extern "C" void FORTRAN_NAME(particle_splitter)(int *nx, int *ny, int *nz,
- int *idual, int *imetal, hydro_method *imethod, float *dt,
- float *r, float *dx, FLOAT *t, float *z,
- float *d1, float *x1, float *v1, float *t1,
- FLOAT *xstart, FLOAT *ystart, FLOAT *zstart, int *ibuff,
- int *npart,
- FLOAT *xpold, FLOAT *ypold, FLOAT *zpold, float *upold, float *vpold, float *wpold,
- float *mpold, float *tdpold, float *tcpold, float *metalfold, int *typeold,
- int *nmax, int *npartnew, int *children, int *level,
- FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
- float *mp, float *tdp, float *tcp, float *metalf, int *type,
- int *iterations, float *separation, int *ran1_init,
- FLOAT *rr_leftedge, FLOAT *rr_rightedge);
-
-
-int grid::ParticleSplitter(int level)
+int grid::ParticleSplitter(int level, int iteration, int NumberOfIDs,
+ long *MustRefineIDs)
{
if (ParticleSplitterIterations == 0)
@@ -131,7 +117,7 @@
/* Allocate space for new particles. */
- int ChildrenPerParent = 12; //12+1 = 13 will be the final number of particles per parent
+ int ChildrenPerParent = CHILDRENPERPARENT; //12+1 = 13 will be the final number of particles per parent
int MaximumNumberOfNewParticles = ChildrenPerParent*NumberOfParticles+1;
tg->AllocateNewParticles(MaximumNumberOfNewParticles);
@@ -145,7 +131,7 @@
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles before splitting = %d, MyProcessorNumber = %d\n",
NumberOfParticles, MyProcessorNumber);
#endif
-
+
if (NumberOfParticles > 0) {
#define NO_PARTICLE_IN_GRID_CHECK
@@ -165,10 +151,29 @@
xindex, yindex, zindex, level);
}
#endif
+
+ /* Before splitting, set particle type to must-refine if given in
+ the provided list. Only do this for the first iteration because
+ in the later iterations, the particle types are propagated to
+ the split particles. */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs != NULL &&
+ iteration == 0) {
+ for (i = 0; i < NumberOfParticles; i++) {
+ for (j = 0; j < NumberOfIDs; j++) {
+ if (this->ParticleNumber[i] == MustRefineIDs[j]) {
+ this->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ break;
+ }
+ } // ENDFOR j
+ } // ENDFOR i
+ }
+
if(!tg->CreateChildParticles(CellWidthTemp, NumberOfParticles, ParticleMass,
ParticleType, ParticlePosition, ParticleVelocity,
ParticleAttribute, CellLeftEdge, GridDimension,
- MaximumNumberOfNewParticles, &NumberOfNewParticles))
+ MaximumNumberOfNewParticles, iteration,
+ &NumberOfNewParticles))
{
fprintf(stdout, "Failed to create child particles in grid %d\n", this->GetGridID());
return FAIL;
@@ -185,6 +190,14 @@
NumberOfNewParticles, MyProcessorNumber);
#endif
+ /* If specified and no ID list given, set all particle types to
+ must-refine */
+
+ if (ParticleSplitterMustRefine == TRUE && MustRefineIDs == NULL) {
+ for (i = 0; i < NumberOfNewParticles; i++)
+ tg->ParticleType[i] = PARTICLE_TYPE_MUST_REFINE;
+ }
+
/* If not set in the above routine, then set the metal fraction to zero. */
if (MetallicityField == FALSE && NumberOfParticleAttributes > 1)
@@ -211,7 +224,7 @@
#ifdef DEBUG_PS
fprintf(stdout, "grid::ParticleSplitter: NumberOfParticles(New) = %d, MyProcessorNumber = %d\n",
- NumberOfParticles, MyProcessorNumber);
+ NumberOfParticles, MyProcessorNumber);
#endif
} // end: if (NumberOfNewParticles > 0)
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/ParticleSplitter.C
--- a/src/enzo/ParticleSplitter.C
+++ b/src/enzo/ParticleSplitter.C
@@ -24,6 +24,7 @@
#include <stdlib.h>
#include <stdio.h>
+#include <hdf5.h>
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "typedefs.h"
@@ -38,7 +39,7 @@
#include "CosmologyParameters.h"
#include "CommunicationUtilities.h"
-#define NO_DEBUG_PS
+#define NO_DEBUG_PS
int RebuildHierarchy(TopGridData *MetaData,
LevelHierarchyEntry *LevelArray[], int level);
@@ -60,15 +61,20 @@
HierarchyEntry **Grids;
int NumberOfGrids;
- /* First, rebuild the hierarchy */
-
- RebuildHierarchy(MetaData, LevelArray, 0);
-
- /* Return if this does not concern us */
+ /* Rebuild hierarchy and return if this does not concern us */
if (ParticleSplitterIterations <= 0 ||
- ParticleSplitterChildrenParticleSeparation <=0)
+ ParticleSplitterChildrenParticleSeparation <=0) {
+ RebuildHierarchy(MetaData, LevelArray, 0);
return SUCCESS;
+ }
+
+ if(ParticleSplitterIterations > MAX_SPLIT_ITERATIONS) {
+ fprintf(stderr, "WARNING: Splitting iterations exceeds maximum allowed\n" \
+ "Resetting to %d\n", MAX_SPLIT_ITERATIONS);
+ ParticleSplitterIterations = MAX_SPLIT_ITERATIONS;
+ }
+
/* Find total NumberOfParticles in all grids; this is needed in
CommunicationUpdateStarParticleCount below */
@@ -76,10 +82,35 @@
MetaData->NumberOfParticles = FindTotalNumberOfParticles(LevelArray);
NumberOfOtherParticles = MetaData->NumberOfParticles - NumberOfStarParticles;
- /* Initialize all star particles if this is a restart */
+ /* Read must-refine particle positions if specified */
+
+ long *MustRefineIDs = NULL;
+ int NumberOfIDs = 0;
+
+ if (ParticleSplitterMustRefine &&
+ ParticleSplitterMustRefineIDFile != NULL) {
- if (ParticleSplitterIterations > 1)
- fprintf(stderr, "WARNING: ParticleSplitterIterations > 1 is not properly tested yet.\n");
+ int NumberOfMustRefineIDs;
+ hid_t file_id, dataset_id, dataspace_id;
+ herr_t status;
+ hsize_t dims[1], maxdims[1];
+
+ file_id = H5Fopen(ParticleSplitterMustRefineIDFile,
+ H5F_ACC_RDONLY, H5P_DEFAULT);
+ dataset_id = H5Dopen2(file_id, "/particle_identifier", H5P_DEFAULT);
+ dataspace_id = H5Dget_space(dataset_id);
+ // Get number of particles and allocate variable
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, maxdims);
+ MustRefineIDs = new long[dims[0]];
+ NumberOfIDs = dims[0];
+ status = H5Dread(dataset_id, H5T_NATIVE_LLONG, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ MustRefineIDs);
+ status = H5Dclose(dataset_id);
+ status = H5Fclose(file_id);
+
+ }
+
+ /* Initialize all star particles if this is a restart */
for (i = 0; i < ParticleSplitterIterations; i++) {
@@ -98,7 +129,6 @@
RecordTotalStarParticleCount(Grids, NumberOfGrids,
TotalStarParticleCountPrevious);
-
// for (grid1 = 0; grid1 < NumberOfGrids; grid1++)
// fprintf(stdout, "TotalStarParticleCountPrevious[grid=%d] = %d\n", grid1,
// TotalStarParticleCountPrevious[grid1]);
@@ -110,7 +140,8 @@
Grids[grid1]->GridData->ReturnNumberOfParticles());
#endif
- if (Grids[grid1]->GridData->ParticleSplitter(level) == FAIL) {
+ if (Grids[grid1]->GridData->ParticleSplitter(level, i, NumberOfIDs,
+ MustRefineIDs) == FAIL) {
ENZO_FAIL("Error in grid::ParticleSplitter.\n");
}
@@ -148,6 +179,7 @@
/* Set splitter parameter zero; otherwise it will split particles again at next restart */
+ delete [] MustRefineIDs;
ParticleSplitterIterations = 0;
return SUCCESS;
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1262,6 +1262,18 @@
&ParticleSplitterRandomSeed);
ret += sscanf(line, "ParticleSplitterChildrenParticleSeparation = %"FSYM,
&ParticleSplitterChildrenParticleSeparation);
+ ret += sscanf(line, "ParticleSplitterMustRefine = %"ISYM,
+ &ParticleSplitterMustRefine);
+ if (sscanf(line, "ParticleSplitterMustRefineIDFile = %s", dummy) == 1)
+ ParticleSplitterMustRefineIDFile = dummy;
+ ret += sscanf(line, "ParticleSplitterFraction = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterFraction+0, ParticleSplitterFraction+1, ParticleSplitterFraction+2,
+ ParticleSplitterFraction+3);
+ ret += sscanf(line, "ParticleSplitterCenter = %"PSYM" %"PSYM" %"PSYM"",
+ ParticleSplitterCenter+0, ParticleSplitterCenter+1, ParticleSplitterCenter+2);
+ ret += sscanf(line, "ParticleSplitterCenterRegion = %"FSYM" %"FSYM" %"FSYM" %"FSYM"",
+ ParticleSplitterCenterRegion+0, ParticleSplitterCenterRegion+1,
+ ParticleSplitterCenterRegion+2, ParticleSplitterCenterRegion+3);
ret += sscanf(line, "ResetMagneticField = %"ISYM,
&ResetMagneticField);
ret += sscanf(line, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM,
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -962,6 +962,14 @@
ParticleSplitterIterations = FALSE;
ParticleSplitterChildrenParticleSeparation = 1.0;
ParticleSplitterRandomSeed = 131180;
+ ParticleSplitterMustRefine = FALSE;
+ ParticleSplitterMustRefineIDFile = NULL;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterFraction[i] = 1.0;
+ for(int i = 0; i < MAX_DIMENSION; i++)
+ ParticleSplitterCenter[i] = -1.0;
+ for(int i = 0; i < MAX_SPLIT_ITERATIONS; i++)
+ ParticleSplitterCenterRegion[i] = -1.0;
/* Magnetic Field Resetter */
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -624,6 +624,16 @@
ParticleSplitterRandomSeed);
fprintf(fptr, "ParticleSplitterChildrenParticleSeparation = %"FSYM"\n",
ParticleSplitterChildrenParticleSeparation);
+ fprintf(fptr, "ParticleSplitterFraction = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterFraction);
+ fprintf(fptr, "ParticleSplitterCenter = ");
+ WriteListOfFloats(fptr, 3, ParticleSplitterCenter);
+ fprintf(fptr, "ParticleSplitterCenterRegion = ");
+ WriteListOfFloats(fptr, MAX_SPLIT_ITERATIONS, ParticleSplitterCenterRegion);
+ fprintf(fptr, "ParticleSplitterMustRefine = %"ISYM"\n",
+ ParticleSplitterMustRefine);
+ fprintf(fptr, "ParticleSplitterMustRefineIDFile = %s\n",
+ ParticleSplitterMustRefineIDFile);
fprintf(fptr, "ResetMagneticField = %"ISYM"\n",
ResetMagneticField);
fprintf(fptr, "ResetMagneticFieldAmplitude = %"GSYM" %"GSYM" %"GSYM"\n",
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -984,6 +984,11 @@
EXTERN int ParticleSplitterIterations;
EXTERN float ParticleSplitterChildrenParticleSeparation;
EXTERN int ParticleSplitterRandomSeed;
+EXTERN int ParticleSplitterMustRefine;
+EXTERN char *ParticleSplitterMustRefineIDFile;
+EXTERN float ParticleSplitterFraction[MAX_SPLIT_ITERATIONS];
+EXTERN FLOAT ParticleSplitterCenter[MAX_DIMENSION];
+EXTERN float ParticleSplitterCenterRegion[MAX_SPLIT_ITERATIONS];
/* Magnetic Field Resetter */
diff -r 603a8a701fe4 -r e58dbfd476cb src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -485,6 +485,12 @@
#define PARTICLE_TYPE_SIMPLE_SOURCE 10
#define CHILDRENPERPARENT 12
+
+/* Splitting particles more than 4 times can induce numerical
+ artifacts */
+
+#define MAX_SPLIT_ITERATIONS 4
+
/* Ways to deposit particles from a subgrid. */
#define CIC_DEPOSIT 0
Repository URL: https://bitbucket.org/enzo/enzo-dev/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
Reply all
Reply to author
Forward
0 new messages