Smartstar

doc/manual/source/parameters/index.rst

@@ -2598,12 +2598,12 @@ SmartStar Feedback
     Default: 0
 
 ``SmartStarStellarRadiativeFeedback`` (external)
-    This parameter controls whether stellar feedback is activated or not. For feedback from PopIII or SMSs then this needs to be on.
-    The stellar radiative feedback is divided up into 5 energy bins. The energy bins have energies of 2.0 eV, 12.8 eV, 14.0 eV, 25.0 eV
-    and 200 eV. The fraction of energy assigned to each bin is determined using the PopIII tables from Schaerer et. al 2002 Table 4.
-    The spectrum for a PopIII star and SMS are different. For a PopIII star a spectrum for a 40 Msolar star is assumed and
-    weighted accordingly. For a SMS a 1000 Msolar star is assumed and weighted accordingly.
-    Future improvements to the SEDs employed here are under active investigation. 
+    This parameter controls whether stellar feedback is activated or not.
+    For feedback from PopII, PopIII or SMSs then this needs to be on.
+    The stellar radiative feedback is divided up into 5 energy bins. See `DetermineSEDParameters.C` for details. 
+    For PopII stars the energy bins are: 2.0 eV, 12.8 eV, 21.62 eV, 30.0 eV and 60 eV.
+    For PopIII stars the energy bins are: 2.0 eV, 12.8 eV, 28.0 eV, 30.0 eV and 58 eV. 
+    For SMS stars the energy bins are: 2.0 eV and 12.8 eV
     Default: 0
 
 ``SmartStarBHFeedback`` (external)
@@ -2637,13 +2637,9 @@ SmartStar Feedback
 
 ``SmartStarSpin`` (external)
     The dimensionless spin of the SmartStar particle. This is a very unconstrained parameter and cannot be readily computed on the fly. This parameter
-    should be set if you want to have jet feedback. Setting this is zero and turning on jet feedback wouldn't make sense. The default is set to be 0.7 and
-    this is probably reasonable. 
-    Default: 0.7
-
-``SmartStarSMSLifetime`` (external)
-    This is the lifetime for a supermassive star in years. After this time has elapsed a SmartStar particle which is behaving like a SMS will collapse
-    directly into a black hole with no supernova event. Default: 1e6
+    should be set if you want to have jet feedback. Setting this is zero and turning on jet feedback wouldn't make sense. The default is set to be 0.674351 and
+    this leads to a radiative efficiency of 1.
+    Default: 0.674351
 
 ``SmartStarJetVelocity`` (external)
     The velocity that the jets are ejected at. Typically jets are observed to travel at a substantial fraction of the speed of light -

doc/manual/source/physics/active_particles.rst

@@ -60,7 +60,8 @@ ______________________________
 
 The ``SmartStar`` particle is built on top of the ``AccretingParticle`` type with additional feedback and accretion protocols attached.
 The ``SmartStar`` particle was designed to be a single particle type that could adjust to the environment
-in which it finds itself. Currently it can represent a PopIII star,
+in which it finds itself. ``SmartStar`` is designed to model resolved star formation.
+Currently it can represent a PopII star, PopIII star,
 a super-massive star or a black hole. However, there is no inherent limit to the physical object it can represent. In that sense
 it may be suitable to augment the ``SmartStar`` particle with your required feature rather than implmenting a new feature. You
 will also be able to build on ``SmartStar`` tests and documentation too rather than starting from scratch. 

run/Hydro/Hydro-3D/BurkertBodenheimerActiveParticles/BBCollapseTestSmartStarParticle.enzo

@@ -0,0 +1,104 @@
+#
+#  AMR PROBLEM DEFINITION FILE: Non-cosmological Collapse test
+#  Description: a sphere collapses until becoming pressure supported.
+#
+#  define problem
+#
+ProblemType                = 27         // Collapse test
+TopGridRank                = 3
+TopGridDimensions          = 64 64 64
+SelfGravity                = 1          // gravity on
+TopGridGravityBoundary     = 0          // periodic
+LeftFaceBoundaryCondition  = 3 3 3      // periodic
+RightFaceBoundaryCondition = 3 3 3
+
+#
+# problem parameters
+#
+
+CollapseTestRefineAtStart   = 1         // check refinement before running
+CollapseTestNumberOfSpheres = 1
+CollapseTestUseParticles    = 0
+CollapseTestInitialTemperature = 1e3    // temperature of the background gas
+CollapseTestSpherePosition[0]   = 0.5 0.5 0.5
+CollapseTestSphereVelocity[0]   = 0.0 0.0 0.0
+CollapseTestSphereRadius[0]     = 0.15625
+CollapseTestSphereCoreRadius[0] = 0.05  // only used with sphere type 5
+CollapseTestSphereDensity[0]    = 1e2   //100   // sphere density, the background density is 1
+CollapseTestSphereTemperature[0] = 10    // put sphere in pressure equilibrium (rho * T is constant)
+CollapseTestSphereRotationPeriod[0] = 8.72734
+CollapseTestFracKeplerianRot[0] = 0.7
+CollapseTestSphereType[0]       = 11     // constant density
+                                        // 1: uniform
+					// 2: r^-2 power-law
+					// 3: NFW
+					// 4: Gaussian
+					// 5: r^-2 power-law with a core
+					// 11: Burkert & Bodenheimer setup
+#EOSType = 
+#  no cosmology for this run
+#
+ComovingCoordinates   = 0              // Expansion OFF
+#
+#  units
+#
+DensityUnits          = 3.82e-20      // 10^4 cm^-3
+LengthUnits           = 3.2e17        #5e17          // 1 pc in cm
+TimeUnits             = 1.0e12        // 10^4 yrs
+GravitationalConstant = 0.0320327594   // 4*pi*G_{cgs}*DensityUnits*TimeUnits^2
+#
+#  set I/O and stop/start parameters
+#
+StopTime          = 1.45
+dtDataDump        = 0.01
+#CycleSkipDataDump = 1
+DataDumpDir       = DD
+DataDumpName      = DD
+OutputTemperature = 1                  // Output temperature field.
+#
+#  set hydro parameters
+#
+Gamma                       = 1.001
+PPMDiffusionParameter       = 0        // diffusion off
+DualEnergyFormalism         = 1        // use total & internal energy
+InterpolationMethod         = 1        // SecondOrderA
+CourantSafetyNumber         = 0.3
+FluxCorrection              = 1
+ConservativeInterpolation   = 0
+RiemannSolver = 4
+HydroMethod                 = 0        // PPM
+SmallRho                    = 1e-8
+SmallE			    = 1e-8
+#
+#  chemistry/cooling
+#
+MultiSpecies                = 0        // chemistry off
+RadiativeCooling            = 0        // cooling off
+Mu  			    = 3.0
+#
+#  set grid refinement parameters
+#
+StaticHierarchy           = 0          // dynamic hierarchy
+MaximumRefinementLevel    = 4         // use up to 4 levels
+RefineBy                  = 2          // refinement factor
+CellFlaggingMethod        = 2 6         // refine on baryon mass and Jeans Length
+MinimumEfficiency         = 0.3
+#OutputFirstTimeAtLevel    = 4         // output when level 4, 5, 6, etc reached (commented out for now)
+#StopFirstTimeAtLevel      = 10         // stop if/when level 10 reached
+RefineByJeansLengthSafetyFactor = 16    // resolve Jeans length by 16 cells (used with CellFlaggingMethod 6)
+
+
+#
+#  set some global parameters
+#
+GreensFunctionMaxNumber   = 10         // # of greens function at any one time
+PotentialIterations	  = 100
+ComputePotential          = 1
+WritePotential            = 0
+#
+#  active particle parameters
+#
+AppendActiveParticleType = SmartStar
+ActiveParticleDensityThreshold = 1e+30   //Jeans density criteria + others
+SmartStarFeedback              = 0
+

run/Hydro/Hydro-3D/BurkertBodenheimerActiveParticles/notes.txt

@@ -0,0 +1,33 @@
+Burkert-Bodenheimer Test
+#This is based off the initial conditions given in Fedderath et al. 2010. It shows the collapse of a
+#rotating sphere. Particles should form in the bar at the centre. The test is not part of the test
+#infrastrure because it is, at the moment, difficult to automate
+
+
+
+# To run the test
+# There is a bit of a palaver to get this running as the mechanism for
+# running enzo with the gravitational potential (GP) turned on is a little
+# difficult - in particular I don't think there is a currently a way to get
+# enzo to run with the GP calculation from the start. So here is how I do it:
+
+1. ./enzo.exe -d BBCollapseTestSmartStarParticle.enzo
+2. This will die after a couple of seconds - that's ok
+3. ./enzo.exe -d -g DD0000/DD0000    #this turns on the gravitational potential
+4. mpirun -np X ./enzo.exe -d -r GravPotential/DD0000/DD0000
+
+This will run relatively slowly i.e. tens of minutes on around 32 cores. This is another reason
+why this is currently not part of the test infrastructure. With the current setup the first particles
+should appears in output DD0133. 
+
+
+#Results
+The test should form 2 particles initially at the ends of a bar that forms in the middle of the
+collapsing sphere. Several particles will subsequently form along the bar as the bar fragments. 
+
+#Plots
+Plots can be made of the collapse using sinkslices.py
+
+
+
+

run/Hydro/Hydro-3D/BurkertBodenheimerActiveParticles/sinkslices.py

@@ -0,0 +1,141 @@
+import matplotlib
+matplotlib.use('Agg')
+import yt
+import numpy as np
+import sys
+import glob
+import matplotlib.pyplot as plt
+from yt.units import G, kboltz, mh
+from yt.units.yt_array import YTQuantity, YTArray
+yt.enable_parallelism()
+Mu = 3.0
+CENTRE = [0.5, 0.5, 0.5]
+#CENTRE = [0.50231934,  0.49035645,  0.49865723]
+AccretionRadius = 4
+BASE = "./GravPotential/"
+fns = glob.glob(BASE + "DD014?/DD????.hierarchy")
+fns.sort()
+#print fns
+WIDTH = 8000
+
+ActiveParticle = "SmartStar"
+
+
+for f in fns:
+    ff = f
+    ds = yt.load(ff)
+    print("Stats = ", ds.index.get_smallest_dx().in_units("au"))
+    dx = ds.index.get_smallest_dx().in_units("au")
+    print("dx = ", dx.d)
+    dd = ds.all_data()
+    centre = CENTRE
+    print("Time = ", ds.current_time.in_units("kyr"))
+    #Now create a slice of the density field through the sink particle
+    #slc = yt.SlicePlot(ds, "x", "density", center='c', width=(0.1, 'pc'))
+    #slc.set_zlim('density', 1e-17, 1e-20)
+    #slc.save("Slice_x_%s.png" % (ds))
+    NumAPs = 0
+    flist = dir(ds.fields)
+    ap = flist[0]
+    try:
+        NumAPs = len(dd[ap, "level"])
+        print("Number of %s active particles = %d" % (flist[0], NumAPs))
+    except:
+        print("No active particles found")
+
+    #v,c = ds.find_max("velocity_magnitude")
+    #print "Max Vel = ", v.in_units("km/s")
+    #print "Loc = ", c
+    #centre = CENTRE
+    #v,c = ds.find_min("velocity_magnitude")
+    #print "Min Vel = ", v.in_units("km/s")
+    #print "Loc = ", c
+    #prj = yt.ProjectionPlot(ds, "x", 'density', center='c',width=(WIDTH, 'au'),
+    #                    weight_field=None)
+    #prj.set_zlim('density', 1e-1, 1e2)
+    #prj.annotate_scale(corner='upper_right')
+    #prj.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
+
+    #prj.save("Proj_x_%s.png" % (ds))
+
+    #slc = yt.SlicePlot(ds, "y", "density", center='c', width=(0.25, 'pc'))
+    #slc.set_zlim('density', 1e-17, 1e-20)
+    #slc.save("Slice_y_%s.png" % (ds))
+    #prj = yt.ProjectionPlot(ds, "y", 'density', center='c',width=(WIDTH, 'au'),
+    #                    weight_field=None)
+    #prj.set_zlim('density', 1e-1, 1e2)
+    #prj.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
+    #prj.annotate_scale(corner='upper_right')
+    #prj.save("Proj_y_%s.png" % (ds))
+
+
+    slc = yt.SlicePlot(ds, "z", "density", center=centre, width=(WIDTH, 'au'))
+    #slc.set_zlim('density', 1e-13, 1e-19)
+    if(NumAPs > 0):
+        for i in range(NumAPs):
+            pos = dd[ActiveParticle, "particle_position"][i]
+            slc.annotate_sphere(pos, radius=(dx.d*AccretionRadius, 'au'),
+                                circle_args={'color':'black', 'fill':True})
+    slc.save("Slice_z_%s.png" % (ds))
+
+    prj = yt.ProjectionPlot(ds, "z", 'number_density', center=centre, width=(WIDTH, 'au'),
+                            weight_field='density')
+    #prj.set_zlim('density', 1e-1, 1e2)
+    prj.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
+    #prj.annotate_scale(corner='upper_right')
+    if(NumAPs > 0):
+        for i in range(NumAPs):
+            pos = dd[ActiveParticle, "particle_position"][i]
+            prj.annotate_sphere(pos, radius=(dx.d*AccretionRadius,'au'),
+                                circle_args={'color':'black', 'fill':True})
+    prj.save("Proj_z_%s.png" % (ds))
+
+    slc = yt.SlicePlot(ds, "z", "velocity_magnitude", center=centre, width=(WIDTH, 'au'))
+    slc.set_unit("velocity_magnitude", "km/s")
+    #slc.set_zlim('velocity_magnitude', 0.01, 100)
+    if(NumAPs > 0):
+        for i in range(NumAPs):
+            pos = dd[ActiveParticle, "particle_position"][i]
+            slc.annotate_sphere(pos, radius=(dx.d*AccretionRadius, 'au'),
+                                circle_args={'color':'black', 'fill':True})
+    slc.save("Slice_z_Vel_%s.png" % (ds))
+
+    prj = yt.ProjectionPlot(ds, "z", 'velocity_magnitude', center=centre,width=(WIDTH, 'au'),
+                            weight_field="velocity_magnitude")
+    prj.set_unit("velocity_magnitude", "km/s")
+    #prj.set_zlim('velocity_magnitude', 0.9, 2)
+    prj.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
+    #prj.annotate_scale(corner='upper_right')
+    if(NumAPs > 0):
+        for i in range(NumAPs):
+            pos = dd[ActiveParticle, "particle_position"][i]
+            prj.annotate_sphere(pos, radius = (dx.d*AccretionRadius, 'au'),
+                                circle_args={'color':'black', 'fill':True})
+    prj.save("Proj_z_Vel_%s.png" % (ds))
+    #print dir(plt)
+    #
+    # Save the image.
+    # Optionally, give a string as an argument
+    # to name files with a keyword.
+    #slc = yt.SlicePlot(ds, "z", "PotentialField", center='c', width=(WIDTH, 'au'))
+    #print "Min/Max = ", np.nanmin(dd["PotentialField"]), np.nanmax(dd["PotentialField"])
+    #slc.set_log('PotentialField', False)
+    #if(NumAPs > 0):
+    #    for i in range(NumAPs):
+    #        pos = dd["SmartStar", "particle_position"][i]
+    #        slc.annotate_sphere(pos, radius=(5, 'au'),
+    #                            circle_args={'color':'black', 'fill':True})
+    #slc.save("PotentialSlice_z_%s.png" % (ds))
+
+
+    #prj = yt.ProjectionPlot(ds, "z", 'PotentialField', center='c',width=(WIDTH, 'au'),
+    #                    weight_field=None)
+    #prj.set_zlim('density', 1e-1, 1e2)
+    #prj.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
+    #prj.annotate_scale(corner='upper_right')
+    #if(NumAPs > 0):
+    #    for i in range(NumAPs):
+    #        pos = dd["SmartStar", "particle_position"][i]
+    #        prj.annotate_sphere(pos, radius=(5, 'au'),
+     #                           circle_args={'color':'black', 'fill':True})
+    #prj.save("Proj_z_%s.png" % (ds))