diff --git a/config/expanse_icc.py b/config/expanse_icc.py deleted file mode 100755 index 0412f43de5..0000000000 --- a/config/expanse_icc.py +++ /dev/null @@ -1,75 +0,0 @@ -import os -from _common_search_paths import charm_path_search, grackle_path_search - -is_arch_valid = 1 -use_gfortran = 0 -smp = 0 - -flags_arch = '-Wall -O3 -g' -#flags_arch = '-fprofile-arcs -ftest-coverage' -#flags_arch = '-Wall -g' -#flags_arch = '-Wall -g -fsanitize=address -fno-omit-frame-pointer' -#flags_arch = '-Wall -O3 -pg' - -# rdynamic required for backtraces -#flags_link_charm = '-rdynamic' -#flags_link_charm = '-memory paranoid' - -intel_dir = '/cm/shared/apps/spack/cpu/opt/spack/linux-centos8-zen/gcc-8.3.1/intel-19.1.1.217-4d42ptjd6wsnh5bgbzcv6lp44vxpjwut/compilers_and_libraries_2020.1.217/linux/bin/intel64' -cc = intel_dir + '/icc' -f90 = intel_dir + '/ifort' - -flags_prec_single = '' -flags_prec_double = '-r8' - -libpath_fortran = '/cm/shared/apps/spack/cpu/opt/spack/linux-centos8-zen/gcc-8.3.1/intel-19.1.1.217-4d42ptjd6wsnh5bgbzcv6lp44vxpjwut/compilers_and_libraries_2020.1.217/linux/compiler/lib/intel64' -libs_fortran = ['ifcore', 'ifport'] - - -#USE GFORTRAN INSTEAD OF IFORT -if use_gfortran: - f90 = 'gfortran' - libpath_fortran = '/cm/shared/apps/spack/cpu/opt/spack/linux-centos8-zen/gcc-8.3.1/gcc-10.2.0-n7su7jf54rc7l2ozegds5xksy6qhrjin/lib64' - libs_fortran = ['gfortran'] - flags_arch_fortran = '-ffixed-line-length-132' - flags_prec_double = '-fdefault-real-8 -fdefault-double-8' - flags_arch = '-O3 -Wall' - flags_fc = '' - -############################# - -home = os.getenv('HOME') - -charm_path = charm_path_search(home) - -use_papi=0 -papi_inc = '/usr/local/include' -papi_lib = '/usr/local/lib' - -boost_path = os.getenv('BOOST_HOME') -boost_inc = boost_path + '/include' -boost_lib = boost_path + '/lib' - -hdf5_inc = os.getenv('HDF5HOME') + '/include' -if hdf5_inc is None: - if os.path.exists('/usr/include/hdf5.h'): - hdf5_inc = '/usr/include' - elif os.path.exists('/usr/include/hdf5/serial/hdf5.h'): - hdf5_inc = '/usr/include/hdf5/serial' - else: - raise Exception('HDF5 include file was not found. Try setting the HDF5_INC environment variable such that $HDF5_INC/hdf5.h exists.') - -hdf5_lib = os.getenv('HDF5HOME') + '/lib' -if hdf5_lib is None: - if os.path.exists('/usr/lib/libhdf5.a'): - hdf5_lib = '/usr/lib' - elif os.path.exists('/usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5.a'): - hdf5_lib = '/usr/lib/x86_64-linux-gnu/hdf5/serial' - else: - raise Exception('HDF5 lib file was not found. Try setting the HDF5_LIB environment variable such that $HDF5_LIB/libhdf5.a exists.') - -png_path = os.getenv('LIBPNG_HOME') -if png_path is None: - png_path = '/lib/x86_64-linux-gnu' - -grackle_path = os.getenv('GRACKLE_HOME') diff --git a/doc/source/param/method.incl b/doc/source/param/method.incl index 8462b2f1a6..f384663cc3 100644 --- a/doc/source/param/method.incl +++ b/doc/source/param/method.incl @@ -706,6 +706,16 @@ most up-to-date description of Grackle parameters, see the `Grackle parameters < :Todo: :o:`write` :Status: **Not accessed** +star_maker +---------- + +..include:: method_star_maker.incl + +feedback +-------- + +..include:: method_feedback.incl + heat ---- diff --git a/doc/source/param/method_feedback.incl b/doc/source/param/method_feedback.incl new file mode 100644 index 0000000000..59cb329903 --- /dev/null +++ b/doc/source/param/method_feedback.incl @@ -0,0 +1,63 @@ +:p: `Method:feedback` parameters. + +---- + +Parameter: :p:`Method` : :p:`feedback` : :p:`single_sn` +:Summary: :s:`Whether to turn on supernova feedback in EnzoMethodFeedbackSTARSS` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Whether to turn on supernova feedback in EnzoMethodFeedbackSTARSS` + +---- + +Parameter: :p:`Method` : :p:`feedback` : :p:`unrestricted_sn` +:Summary: :s:`Whether to turn on supernova feedback in EnzoMethodFeedbackSTARSS` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Allow for > 1 supernova event per star particle in a timestep.` + +---- + +Parameter: :p:`Method` : :p:`feedback` : :p:`stellar_winds` +:Summary: :s:`Whether to turn on stellar winds in EnzoMethodFeedbackSTARSS` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Whether to turn on stellar wind feedback in EnzoMethodFeedbackSTARSS` + +---- + +Parameter: :p:`Method` : :p:`feedback` : :p:`analytic_SNR_shell_mass` +:Summary: :s:`Calculates supernova remnant shell mass via analytic formulae` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Calculates supernova remnant shell mass via analytic formulae` + +---- + +Parameter: :p:`Method` : :p:`feedback` : :p:`fade_SNR` +:Summary: :s:`Allow coupling to fading phase of SNe in EnzoMethodFeedbackSTARSS` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Allow coupling to fading phase of SNe in EnzoMethodFeedbackSTARSS if cell width is greater than fading radius.` + +---- + +Parameter: :p:`Method` : :p:`feedback` : :p:`NEvents` +:Summary: :s:`Manually set off supernovae in test problems` +:Type: :t:`integer` +:Default: :d:`-1` +:Scope: :z:`Enzo` + +:e:`If -1, will probabilistically model supernova. If "NEvents" > 1, will set off N-supernovae per particle, (1 per timestep). Mostly for testing purposes` + +---- \ No newline at end of file diff --git a/doc/source/param/method_star_maker.incl b/doc/source/param/method_star_maker.incl new file mode 100644 index 0000000000..3d0773839f --- /dev/null +++ b/doc/source/param/method_star_maker.incl @@ -0,0 +1,191 @@ +:p: `Method:star_maker` parameters. + +---- + +:Parameter: :p:`Method` : :p:`star_maker` : :p:`flavor` +:Summary: :s:`Which star_maker method to use` +:Type: :t:`string` +:Default: :d:`STARSS` +:Scope: :z:`Enzo` + +:e:`Options: "STARSS", "stochastic"` + +---- + +:Parameter: :p:`Method` : :p:`star_maker` : :p:`use_density_threshold` +:Summary: :s:`Use number density threshold for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Flag to enable number density threshold for star formation.` + +---- + +:Parameter: :p:`Method` : :p:`star_maker` : :p:`number_density_threshold` +:Summary: :s:`Use number density threshold for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Set number density threshold for star formation. Requires "use_density_threshold"=true.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_overdensity_threshold` +:Summary: :s:`Use overdensity threshold` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Flag to enable overdensity threshold for star formation.` + +---- + +:Parameter: :p:`Method` : :p:`star_maker` : :p:`overdensity_threshold` +:Summary: :s:`Use overdensity threshold for star formation` +:Type: :t:`float` +:Default: :d:`0.0` +:Scope: :z:`Enzo` + +:e:`Set overdensity threshold for star formation. Requires "use_overdensity_threshold"=true.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_velocity_divergence` +:Summary: :s:`Use converging flow criterion for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Flag to check whether div(V) < 0` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_cooling_time` +:Summary: :s:`Check if cooling_time < dynamical_time for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Flag to check if cooling_time < dynamical_time` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_temperature_threshold` +:Summary: :s:`Use temperature threshold for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Flag to enable temperature threshold check for star formation` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`temperature_threshold` +:Summary: :s:`Temperature threshold for star formation` +:Type: :t:`float` +:Default: :d:`1e4` +:Scope: :z:`Enzo` + +:e:`Set temperature threshold required for star formation. Requires "use_temperature_threshold"=true.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_self_gravitating` +:Summary: :s:`Use FIRE2 virial parameter criterion for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Checks that alpha < 1, where alpha is the virial parameter calculated using the FIRE2 prescription.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_altAlpha` +:Summary: :s:`Use alternate virial parameter criterion for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Checks that alpha < 1, where alpha is the virial parameter calculated as "potential/total_energy". Currently only accessed by EnzoMethodStarMakerSTARSS. Requires DEBUG_COPY_POTENTIAL flag to be set in EnzoMethodGravity.cpp` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_h2_self_shielding` +:Summary: :s:`Use H2 self-shielding criterion for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Checks that f_shield < 0, f_shield is calculated using fits from Krumholz & Gnedin (2011).' + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`use_jeans_mass` +:Summary: :s:`Use Jeans mass criterion for star formation` +:Type: :t:`logical` +:Default: :d:`false` +:Scope: :z:`Enzo` + +:e:`Checks that cell_mass > jeans_mass in a cell' + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`critical_metallicity` +:Summary: :s:`Metallicity threshold for star formation` +:Type: :t:`float` +:Default: :d:`0.0` +:Scope: :z:`Enzo` + +:e:`Set metallicity threshold required for star formation` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`maximum_mass_fraction` +:Summary: :s:`Max fraction of gas in a cell that can be converted into a star particle per formation event.` +:Type: :t:`float` +:Default: :d:`0.5` +:Scope: :z:`Enzo` + +:e:`Max fraction of gas in a cell that can be converted into a star particle per formation event.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`min_level` +:Summary: :s:`Minimum AMR level required for star formation.` +:Type: :t:`integer` +:Default: :d:`0` +:Scope: :z:`Enzo` + +:e:`Set minimum AMR level required for star formation.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`minimum_star_mass` +:Summary: :s:`Minimum star particle mass` +:Type: :t:`float` +:Default: :d:`0.0` +:Scope: :z:`Enzo` + +:e:`Set minimum star particle mass.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`maximum_star_mass` +:Summary: :s:`Maximum star particle mass` +:Type: :t:`float` +:Default: :d:`1e4` +:Scope: :z:`Enzo` + +:e:`Set maximum star particle mass. For no limit, set "maximum_star_mass" < 0.` + +---- + +Parameter: :p:`Method` : :p:`star_maker` : :p:`turn_off_probability` +:Summary: :s:`Turn off probablistic elements of EnzoMethodStarMakerSTARSS.` +:Type: :t:`logical` +:Default: :d:`fapse` +:Scope: :z:`Enzo` + +:e:`Turn off probablistic elements of EnzoMethodStarMakerSTARSS. Mostly meant for debugging.` \ No newline at end of file diff --git a/input/STARSS/SF_FB.incl b/input/STARSS/SF_FB.incl new file mode 100644 index 0000000000..54d06921e8 --- /dev/null +++ b/input/STARSS/SF_FB.incl @@ -0,0 +1,78 @@ +Particle { + list += ["star"]; + star { + attributes = [ "x", "default", + "y", "default", + "z", "default", + "vx", "default", + "vy", "default", + "vz", "default", + "ax", "default", + "ay", "default", + "az", "default", + "mass", "default", + "creation_time", "default", + "creation_level", "default", + "lifetime", "default", + "metal_fraction", "default", + "is_copy", "int64", + "number_of_sn", "int64", + "id", "int64" ]; + position = [ "x", "y", "z" ]; + velocity = [ "vx", "vy", "vz" ]; + group_list = ["is_gravitating"]; + } +} + +Method { + star_maker { + flavor = "STARSS"; # pick star formation method + + min_level = 4; # minimum level for star formation + + use_density_threshold = false; # check number density for SF + use_overdensity_threshold= false; # check overdensity + + use_velocity_divergence = true; # converging flow criterion + use_cooling_time = false; # check that cooling_time < dynamical_time + + use_self_gravitating = true; # check that virial parameter < 1 + use_altAlpha = false; # replace virial parameter calculation with total_energy/potential + # requires DEBUG_COPY_POTENTIAL to be set in EnzoMethodGravity + use_h2_self_shielding = false; + use_temperature_threshold= false; + use_critical_metallicity = false; + use_jeans_mass = true; + + overdensity_threshold = 512.0; + minimum_star_mass = 500.0; + maximum_star_mass = 1000.0; + number_density_threshold = 100.0; + temperature_threshold = 1e4; + critical_metallicity = 4.1e-8 / 0.012; + turn_off_probability = false; # turn off probablistic part of SF for testing + maximum_mass_fraction = 0.05; # maximum fraction of a cell that can be converted into stars + }; + + feedback { + flavor = "STARSS"; + single_sn = 1; # whether to have SNe at all + unrestricted_sn = 1; # Allow > 1 SN per particle in a timestep + stellar_winds = 1; # whether to have stellar winds + analytic_SNR_shell_mass = 1; # calculate shell mass and subtract this from blast interior + # during feedback deposition + + fade_SNR = 1; # Allow coupling of FB out to the fading phase, where expansion of the SNR + # is comparable to the local sound speed. + + NEvents = -1; # Use this to manually set off SN events in ideal tests. + # if -1, SN rates are calculated for the particles to determine when + # they blow up + }; + + output { + particle_list += ["star"]; + } + + +} diff --git a/input/STARSS/method_feedback_starss.in b/input/STARSS/method_feedback_starss.in new file mode 100644 index 0000000000..9038f3c518 --- /dev/null +++ b/input/STARSS/method_feedback_starss.in @@ -0,0 +1,236 @@ +# +# +# Very simple feedback test problem +# +# Drops a star particle of desired mass +# into a unigrid, uniform domain. Commented-out +# Adapt blocks below can be used to add in AMR +# (static, nested grids or adaptive). +# + +Adapt { + list = []; + mask { + type = "mask"; + value = [1.0, y<=0.5, 0.0]; + } + max_level = 0; + } + +Boundary { + type = "periodic"; +} + +Domain { + lower = [ 0.0, 0.0, 0.0]; + upper = [ 1.0, 1.0, 1.0]; +} + +Mesh { + root_blocks = [4,4,4]; + root_rank = 3; + root_size = [64,64,64]; # given length units, res = 1024 pc / 64 = 16 pc + + +} + +Field { + alignment = 8; + gamma = 1.40; + ghost_depth = 4; + + list = ["density", "internal_energy", "total_energy", + "velocity_x", "velocity_y", "velocity_z", + "pressure", "temperature"]; + list += ["HI_density","HII_density","HeI_density","HeII_density","HeIII_density","e_density","metal_density"]; +} + +Group { + list = ["color","derived","conserved","make_field_conservative"]; + + color { + field_list = ["metal_density"]; + field_list += ["HI_density","HII_density","HeI_density","HeII_density","HeIII_density","e_density"]; + } + derived { + field_list = ["temperature","pressure"]; + } + conserved { + field_list = ["density", "HI_density", "HII_density", "HeI_density", "HeII_density", "HeIII_density", "e_density", "metal_density"]; + } + make_field_conservative { + field_list = ["velocity_x","velocity_y","velocity_z","total_energy","internal_energy"]; + } + +} + +Method { + list = ["ppm", "grackle", "feedback"]; + + null { + dt = 10.0; # force a minimum dt - code units + }; + + flux_correct {enable=false;} + + feedback { + flavor = "STARSS"; + single_sn = 1; + unrestricted_sn = 1; + stellar_winds = 1; + analytic_SNR_shell_mass = 1; + fade_SNR = 1; + NEvents = 1; + }; + + ppm { + diffusion = true; + riemann_solver = "two_shock"; + dual_energy = true; + flattening = 3; + steepening = true; + mol_weight = 1.2; + courant = 0.04; + density_floor = 1.0E-20; + number_density_floor = 1.0E-20; + pressure_floor = 1.0E-20; + temperature_floor = 1.0E-20; + }; + + grackle { + data_file = "input/CloudyData_UVB=HM2012.h5"; + with_radiative_cooling = 1; + metal_cooling = 1; + primordial_chemistry = 1; + UVbackground = 0; + } ; +} + +Particle { + list = ["star"]; + + star { + attributes = [ "x", "double", + "y", "double", + "z", "double", + "vx", "double", + "vy", "double", + "vz", "double", + "ax", "double", + "ay", "double", + "az", "double", + "mass", "double", + "creation_time", "double", + "lifetime", "double", + "number_of_sn", "int64", + "metal_fraction", "double", + "is_copy", "int64", + "id", "int64"]; + position = [ "x", "y", "z" ]; + velocity = [ "vx", "vy", "vz" ]; + group_list = ["is_gravitating"]; + } +} +Units { + length = 64.0 * 10.0 *3.0866E18; # middle number = cell resolution in pc + time = 3.15576E13; # 1 Myr + density = 1.2E-24; + } + +Initial { + + list = ["feedback_test"]; # name of IC problem + + feedback_test { + density = 4.0*1.2E-24; + HI_density = 0.7 * 4.0*1.2E-24; + HII_density = 1e-10 * 4.0*1.2E-24; + HeI_density = 0.3 * 4.0*1.2E-24; + HeII_density = 1e-10 * 4.0*1.2E-24; + HeIII_density = 1e-10 * 4.0*1.2E-24; + e_density = 1e-10 * 4.0*1.2E-24; # uniform mass density for ICs + metal_fraction = 1e-2*0.012; + + position = [0.5,0.5,0.5]; # unigrid corner (or one block) + + temperature = 100.0; # in K + # particles are still local to each block and are kicked away from + # boundaries if NxN feedback zone crosses a block boundary + star_mass = 1000.0; # particle mass in Msun + } + +} + +Output { + list = ["de","te","vx"]; + de { + dir = [ "FEEDBACK_TEST_%04d", "cycle" ]; + field_list = [ "density" ]; + image_ghost = false; + image_size = [ 512, 512 ]; + image_type = "data"; + name = [ "de-%04d.png", "cycle" ]; + schedule { + start = 0; #1; + step = 1; + var = "cycle"; + }; + type = "image"; + }; + te { + dir = [ "FEEDBACK_TEST_%04d", "cycle" ]; + field_list = [ "total_energy" ]; + image_ghost = false; + image_size = [ 512, 512 ]; + image_type = "data"; + name = [ "te-%04d.png", "cycle" ]; + schedule { + start = 0; #1; + step = 1; + var = "cycle"; + }; + type = "image"; + }; + vx { + dir = [ "FEEDBACK_TEST_%04d", "cycle" ]; + field_list = [ "velocity_x" ]; + image_ghost = false; + image_size = [ 512, 512 ]; + image_type = "data"; + name = [ "vx-%04d.png", "cycle" ]; + schedule { + start = 0; #1; + step = 1; + var = "cycle"; + }; + type = "image"; + }; + + data { + field_list = [ "density" , "metal_density","velocity_x","velocity_y","velocity_z", + "pressure", "temperature","total_energy","internal_energy"]; + + field_list += ["HI_density","HII_density","HeI_density","HeII_density","HeIII_density","e_density"]; + + + particle_list = ["star"]; + + dir = ["FEEDBACK_TEST_%04d","cycle"]; + name = [ "feedback-test-data-%04d-%03d.h5", "cycle", "proc" ]; + type = "data"; + + # schedule the output based on 'var' and 'step' + # step refers to the interval if 'var' to output on + schedule { + var = "cycle"; + start = 0; + step = 1; # time in code units (Myr) + } + }; + + +} + +Stopping { + cycle = 401; +} diff --git a/input/bb_unigrid_SF_FB.in b/input/bb_unigrid_SF_FB.in index 62bca7da14..3870010b24 100644 --- a/input/bb_unigrid_SF_FB.in +++ b/input/bb_unigrid_SF_FB.in @@ -131,7 +131,7 @@ Particle { }; list = [ "null", "ppm", - "grackle", "star_maker","feedback", + "grackle", #"star_maker","feedback", "pm_deposit", "gravity", "pm_update"]; null { @@ -139,7 +139,7 @@ Particle { }; grackle { - data_file = "CloudyData_noUVB.h5"; + data_file = "input/CloudyData_noUVB.h5"; with_radiative_cooling = 1; primordial_chemistry = 0; metal_cooling = 1; @@ -202,7 +202,7 @@ Particle { }; type = "data"; }; - list = [ "data" ]; + #list = [ "data" ]; } Solver { diff --git a/input/test_cosmo-dd-SF_FB.in b/input/test_cosmo-dd-SF_FB.in new file mode 100644 index 0000000000..fd2195f3ab --- /dev/null +++ b/input/test_cosmo-dd-SF_FB.in @@ -0,0 +1,39 @@ +include "input/test_cosmo-dd-multispecies.in" +include "input/STARSS/SF_FB.incl" + +# This is a cosmology problem that is set up to use six-species chemistry, star formation, and feedback. +# Since this is a very small problem (32^3 root grid), the first star particle won't form until around redshift 8. +# To get to this point, it should take around 15-20 minutes on 8 processors. + +Method { + list = [ "star_maker", "feedback", "pm_deposit", "gravity", "ppm", "grackle", "pm_update", "comoving_expansion"]; + output { particle_list += ["star"];} +} + +Output { + list = ["hdf5","de","star"]; + hdf5 { + dir = ["Dir_COSMO_SF_FB_%04d","cycle"]; + particle_list += ["star"]; + } + de { dir = [ "Dir_COSMO_SF_FB_%04d", "cycle" ]; } + hdf5 { dir = [ "Dir_COSMO_SF_FB_%04d", "cycle" ]; } + + star { + colormap = [ 0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000, 1.000000000000000, 0.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000 ]; + dir = [ "Dir_COSMO_SF_FB_%04d", "cycle" ]; + image_ghost = false; + image_size = [ 1024, 1024 ]; + image_type = "data"; + name = [ "star-%02d.png", "count" ]; + particle_list = [ "star" ]; + schedule { + start = 0; + step = 100; + var = "cycle"; + }; + type = "image"; + }; +} + +Stopping {cycle = 1000;} diff --git a/input/test_cosmo-dd-fc0.in b/input/test_cosmo-dd-fc0.in index d6ffba10ad..ac1f6977f0 100644 --- a/input/test_cosmo-dd-fc0.in +++ b/input/test_cosmo-dd-fc0.in @@ -7,6 +7,8 @@ Adapt { min_level = -2; } }; } +Stopping { cycle = 300; } + Solver { list = [ "dd", "dd_root", "dd_domain", "dd_smooth", "root_coarse", "root_pre", "root_post" ]; diff --git a/input/test_cosmo-dd-multispecies.in b/input/test_cosmo-dd-multispecies.in index 604620a8f4..0a1dbe5a52 100755 --- a/input/test_cosmo-dd-multispecies.in +++ b/input/test_cosmo-dd-multispecies.in @@ -72,7 +72,7 @@ Method { Output { - list = ["hdf5"]; + list = ["hdf5","de"]; hdf5 { dir = ["Dir_COSMO_MULTI_%04d","cycle"]; field_list += ["HI_density", @@ -86,16 +86,6 @@ Output { "temperature"]; } de { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - depa { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - ax { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - ay { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - az { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - dark { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - mesh { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - po { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - hdf5 { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - dep { dir = [ "Dir_COSMO_MULTI_%04d", "cycle" ]; } - check { dir = [ "Dir_COSMO_MULTI_%04d-checkpoint", "count" ]; } } Stopping {cycle = 160;} diff --git a/input/test_cosmo.incl b/input/test_cosmo.incl index 2240e0f35d..5cefb9df17 100644 --- a/input/test_cosmo.incl +++ b/input/test_cosmo.incl @@ -133,7 +133,7 @@ Adapt { order = 2; solver = "cg"; }; - list = ["pm_deposit", "gravity", "pm_update", "comoving_expansion" ]; + list = ["ppm", "pm_deposit", "gravity", "pm_update", "comoving_expansion" ]; ppm { courant = 0.5000000000000000; diffusion = false; @@ -142,7 +142,7 @@ Adapt { } Output { - list = [ "de", "depa", "ax", "dark", "mesh", "po", "check" ]; + list = [ "de", "depa", "ax", "dark", "mesh", "po", "check", "hdf5" ]; de { dir = [ "COSMO_CG_%04d", "cycle" ]; field_list = [ "density" ]; diff --git a/src/Enzo/_enzo.hpp b/src/Enzo/_enzo.hpp index 526d51d395..47907a9614 100644 --- a/src/Enzo/_enzo.hpp +++ b/src/Enzo/_enzo.hpp @@ -236,6 +236,7 @@ extern "C" { #include "enzo_EnzoMethodCosmology.hpp" #include "enzo_EnzoMethodFeedback.hpp" #include "enzo_EnzoMethodDistributedFeedback.hpp" +#include "enzo_EnzoMethodFeedbackSTARSS.hpp" #include "enzo_EnzoMethodGrackle.hpp" #include "enzo_EnzoMethodGravity.hpp" #include "enzo_EnzoMethodHeat.hpp" @@ -247,6 +248,7 @@ extern "C" { #include "enzo_EnzoMethodTurbulence.hpp" #include "enzo_EnzoMethodStarMaker.hpp" #include "enzo_EnzoMethodStarMakerStochasticSF.hpp" +#include "enzo_EnzoMethodStarMakerSTARSS.hpp" #include "enzo_EnzoMethodMHDVlct.hpp" #include "enzo_EnzoMethodMergeSinks.hpp" #include "enzo_EnzoMethodAccretion.hpp" diff --git a/src/Enzo/cic_deposit.F b/src/Enzo/cic_deposit.F index a3d66f04a0..74419aec6b 100644 --- a/src/Enzo/cic_deposit.F +++ b/src/Enzo/cic_deposit.F @@ -69,7 +69,7 @@ subroutine cic_deposit(posx, posy, posz, ndim, npositions, ! write(0,*) npositions, leftedge, dim1, dim2, dim3, cellsize - print*, 'cic_deposit.F mass(1)=',mass(1) +c print*, 'cic_deposit.F mass(1)=',mass(1) fact = 1._PKIND/cellsize refine = cellsize/cloudsize half = 0.5001_PKIND/refine @@ -155,7 +155,7 @@ subroutine cic_deposit(posx, posy, posz, ndim, npositions, do n=1, npositions c c Compute the position of the central cell -c +c xpos = min(max((posx(n) - leftedge(1))*fact, half),edge1) ypos = min(max((posy(n) - leftedge(2))*fact, half),edge2) zpos = min(max((posz(n) - leftedge(3))*fact, half),edge3) @@ -175,7 +175,7 @@ subroutine cic_deposit(posx, posy, posz, ndim, npositions, c c Interpolate from field into sumfield c -C print*, 'DEBUG_CIC',dx,dy,dz +c print*, 'DEBUG_CIC ',dx,dy,dz,i1,j1,k1,field(i1,j1,k1) field(i1 ,j1 ,k1 ) = field(i1 ,j1 ,k1 ) + & mass(n)* dx * dy * dz field(i1+1,j1 ,k1 ) = field(i1+1,j1 ,k1 ) + @@ -195,7 +195,6 @@ subroutine cic_deposit(posx, posy, posz, ndim, npositions, c enddo - endif c return diff --git a/src/Enzo/enzo.ci b/src/Enzo/enzo.ci index aa06716637..47f18367d1 100644 --- a/src/Enzo/enzo.ci +++ b/src/Enzo/enzo.ci @@ -87,6 +87,7 @@ module enzo { PUPable EnzoMethodBackgroundAcceleration; PUPable EnzoMethodFeedback; PUPable EnzoMethodDistributedFeedback; + PUPable EnzoMethodFeedbackSTARSS; PUPable EnzoMethodComovingExpansion; PUPable EnzoMethodCosmology; @@ -115,6 +116,7 @@ module enzo { PUPable EnzoMethodStarMaker; PUPable EnzoMethodStarMakerStochasticSF; + PUPable EnzoMethodStarMakerSTARSS; PUPable EnzoBfieldMethodCT; PUPable EnzoPhysicsCosmology; @@ -161,7 +163,10 @@ module enzo { entry EnzoBlock (MsgRefine * msg); #endif entry EnzoBlock(); - + + // EnzoMethodFeedbackSTARSS synchronization entry methods + entry void p_method_feedback_starss_end(); + // EnzoMethodTurbulence synchronization entry methods entry void r_method_turbulence_end(CkReductionMsg *msg); diff --git a/src/Enzo/enzo_EnzoBlock.hpp b/src/Enzo/enzo_EnzoBlock.hpp index b74bb5bc3c..00216e2bb1 100644 --- a/src/Enzo/enzo_EnzoBlock.hpp +++ b/src/Enzo/enzo_EnzoBlock.hpp @@ -331,6 +331,9 @@ class EnzoBlock : public CBase_EnzoBlock void solver_mg0_prolong_recv(FieldMsg * msg); void p_solver_mg0_restrict_recv(FieldMsg * msg); + // EnzoMethodFeedbackSTARSS + + void p_method_feedback_starss_end(); void print() { Block::print(); diff --git a/src/Enzo/enzo_EnzoConfig.cpp b/src/Enzo/enzo_EnzoConfig.cpp index cb921df785..bbe6d9b06d 100644 --- a/src/Enzo/enzo_EnzoConfig.cpp +++ b/src/Enzo/enzo_EnzoConfig.cpp @@ -72,6 +72,12 @@ EnzoConfig::EnzoConfig() throw () initial_collapse_temperature(0.0), // EnzoInitialFeedbackTest initial_feedback_test_density(), + initial_feedback_test_HI_density(), + initial_feedback_test_HII_density(), + initial_feedback_test_HeI_density(), + initial_feedback_test_HeII_density(), + initial_feedback_test_HeIII_density(), + initial_feedback_test_e_density(), initial_feedback_test_star_mass(), initial_feedback_test_temperature(), initial_feedback_test_from_file(), @@ -211,6 +217,7 @@ EnzoConfig::EnzoConfig() throw () method_hydro_reconstruct_positive(0), method_hydro_riemann_solver(""), // EnzoMethodFeedback, + method_feedback_flavor(""), method_feedback_ejecta_mass(0.0), method_feedback_supernova_energy(1.0), method_feedback_ejecta_metal_fraction(0.0), @@ -220,19 +227,37 @@ EnzoConfig::EnzoConfig() throw () method_feedback_ke_fraction(0.0), method_feedback_use_ionization_feedback(false), method_feedback_time_first_sn(-1), // in Myr + // EnzoMethodFeedbackSTARSS, + method_feedback_single_sn(0), + method_feedback_unrestricted_sn(0), + method_feedback_stellar_winds(0), + method_feedback_gas_return_fraction(0.0), + method_feedback_min_level(0), + method_feedback_analytic_SNR_shell_mass(0), + method_feedback_fade_SNR(0), + method_feedback_NEvents(0), // EnzoMethodStarMaker, method_star_maker_flavor(""), // star maker type to use - method_star_maker_use_density_threshold(true), // check above density threshold before SF - method_star_maker_use_velocity_divergence(true), // check for converging flow before SF - method_star_maker_use_dynamical_time(true), // compute t_ff / t_dyn. Otherwise take as 1.0 + method_star_maker_use_altAlpha(false), + method_star_maker_use_density_threshold(false), // check above density threshold before SF + method_star_maker_use_velocity_divergence(false), // check for converging flow before SF + method_star_maker_use_dynamical_time(false), // compute t_ff / t_dyn. Otherwise take as 1.0 + method_star_maker_use_cooling_time(false), // check if t_cool < t_dyn + method_star_maker_use_overdensity_threshold(false), + method_star_maker_use_temperature_threshold(false), method_star_maker_use_self_gravitating(false), // method_star_maker_use_h2_self_shielding(false), method_star_maker_use_jeans_mass(false), method_star_maker_number_density_threshold(0.0), // Number density threshold in cgs + method_star_maker_overdensity_threshold(0.0), + method_star_maker_critical_metallicity(0.0), + method_star_maker_temperature_threshold(1.0E4), method_star_maker_maximum_mass_fraction(0.5), // maximum cell mass fraction to convert to stars method_star_maker_efficiency(0.01), // star maker efficiency per free fall time method_star_maker_minimum_star_mass(1.0E4), // minimum star particle mass in solar masses method_star_maker_maximum_star_mass(1.0E4), // maximum star particle mass in solar masses + method_star_maker_min_level(0), // minimum AMR level for star formation + method_star_maker_turn_off_probability(false), // EnzoMethodTurbulence method_turbulence_edot(0.0), method_turbulence_mach_number(0.0), @@ -492,6 +517,12 @@ void EnzoConfig::pup (PUP::er &p) PUParray(p, initial_feedback_test_position,3); p | initial_feedback_test_density; + p | initial_feedback_test_HI_density; + p | initial_feedback_test_HII_density; + p | initial_feedback_test_HeI_density; + p | initial_feedback_test_HeII_density; + p | initial_feedback_test_HeIII_density; + p | initial_feedback_test_e_density; p | initial_feedback_test_star_mass; p | initial_feedback_test_temperature; p | initial_feedback_test_from_file; @@ -558,6 +589,7 @@ void EnzoConfig::pup (PUP::er &p) p | method_hydro_reconstruct_positive; p | method_hydro_riemann_solver; + p | method_feedback_flavor; p | method_feedback_ejecta_mass; p | method_feedback_supernova_energy; p | method_feedback_ejecta_metal_fraction; @@ -568,18 +600,37 @@ void EnzoConfig::pup (PUP::er &p) p | method_feedback_use_ionization_feedback; p | method_feedback_time_first_sn; + p | method_feedback_single_sn; + p | method_feedback_unrestricted_sn; + p | method_feedback_stellar_winds; + p | method_feedback_gas_return_fraction; + p | method_feedback_min_level; + p | method_feedback_analytic_SNR_shell_mass; + p | method_feedback_fade_SNR; + p | method_feedback_NEvents; + p | method_star_maker_flavor; + p | method_star_maker_use_altAlpha; p | method_star_maker_use_density_threshold; + p | method_star_maker_use_overdensity_threshold; + p | method_star_maker_use_temperature_threshold; + p | method_star_maker_use_critical_metallicity; p | method_star_maker_use_velocity_divergence; p | method_star_maker_use_dynamical_time; + p | method_star_maker_use_cooling_time; p | method_star_maker_use_self_gravitating; p | method_star_maker_use_h2_self_shielding; p | method_star_maker_use_jeans_mass; p | method_star_maker_number_density_threshold; + p | method_star_maker_overdensity_threshold; + p | method_star_maker_critical_metallicity; + p | method_star_maker_temperature_threshold; p | method_star_maker_maximum_mass_fraction; p | method_star_maker_efficiency; p | method_star_maker_minimum_star_mass; p | method_star_maker_maximum_star_mass; + p | method_star_maker_min_level; + p | method_star_maker_turn_off_probability; p | method_turbulence_edot; @@ -1231,6 +1282,24 @@ void EnzoConfig::read_initial_feedback_test_(Parameters * p) initial_feedback_test_density = p->value_float ("Initial:feedback_test:density", 1.0E-24); + initial_feedback_test_HI_density = p->value_float + ("Initial:feedback_test:HI_density", 1.0E-24); + + initial_feedback_test_HII_density = p->value_float + ("Initial:feedback_test:HII_density", 1.0E-100); + + initial_feedback_test_HeI_density = p->value_float + ("Initial:feedback_test:HeI_density", 1.0E-100); + + initial_feedback_test_HeII_density = p->value_float + ("Initial:feedback_test:HeII_density", 1.0E-100); + + initial_feedback_test_HeIII_density = p->value_float + ("Initial:feedback_test:HeIII_density", 1.0E-100); + + initial_feedback_test_e_density = p->value_float + ("Initial:feedback_test:e_density", 1.0E-100); + initial_feedback_test_star_mass = p->value_float ("Initial:feedback_test:star_mass", 1000.0); @@ -1313,7 +1382,6 @@ void EnzoConfig::read_method_grackle_(Parameters * p) method_grackle_metallicity_floor = p-> value_float ("Method:grackle:metallicity_floor", 0.0); - // Set Grackle parameters from parameter file method_grackle_chemistry->with_radiative_cooling = p->value_integer ("Method:grackle:with_radiative_cooling", @@ -1339,6 +1407,14 @@ void EnzoConfig::read_method_grackle_(Parameters * p) ("Method:grackle:cmb_temperature_floor", method_grackle_chemistry->cmb_temperature_floor); + method_grackle_chemistry->h2_charge_exchange_rate = p->value_integer + ("Method:grackle:h2_charge_exchange_rate", + method_grackle_chemistry->h2_charge_exchange_rate); + + method_grackle_chemistry->h2_h_cooling_rate = p->value_integer + ("Method:grackle:h2_h_cooling_rate", + method_grackle_chemistry->h2_h_cooling_rate); + std::string grackle_data_file_ = p->value_string ("Method:grackle:data_file", ""); ASSERT("EnzoConfig::read", @@ -1443,6 +1519,9 @@ void EnzoConfig::read_method_grackle_(Parameters * p) void EnzoConfig::read_method_feedback_(Parameters * p) { + method_feedback_flavor = p->value_string + ("Method:feedback:flavor","distributed"); + method_feedback_ejecta_mass = p->value_float ("Method:feedback:ejecta_mass",0.0); @@ -1469,6 +1548,31 @@ void EnzoConfig::read_method_feedback_(Parameters * p) method_feedback_use_ionization_feedback = p->value_logical ("Method:feedback:use_ionization_feedback", false); + + // MethodFeedbackSTARSS parameters + method_feedback_single_sn = p->value_integer + ("Method:feedback:single_sn",0); + + method_feedback_unrestricted_sn = p->value_integer + ("Method:feedback:unrestricted_sn",0); + + method_feedback_stellar_winds = p->value_integer + ("Method:feedback:stellar_winds",0); + + method_feedback_gas_return_fraction = p->value_float + ("Method:feedback:gas_return_fraction",0.0); + + method_feedback_min_level = p->value_integer + ("Method:feedback:min_level",0); + + method_feedback_analytic_SNR_shell_mass = p->value_integer + ("Method:feedback:analytic_SNR_shell_mass",0); + + method_feedback_fade_SNR = p->value_integer + ("Method:feedback:fade_SNR",0); + + method_feedback_NEvents = p->value_integer + ("Method:feedback:NEvents",-1); } //---------------------------------------------------------------------- @@ -1479,14 +1583,23 @@ void EnzoConfig::read_method_star_maker_(Parameters * p) method_star_maker_flavor = p->value_string ("Method:star_maker:flavor","stochastic"); + method_star_maker_use_altAlpha = p->value_logical + ("Method:star_maker:use_altAlpha",false); + method_star_maker_use_density_threshold = p->value_logical - ("Method:star_maker:use_density_threshold",true); + ("Method:star_maker:use_density_threshold",false); + + method_star_maker_use_overdensity_threshold = p->value_logical + ("Method:star_maker:use_overdensity_threshold",false); method_star_maker_use_velocity_divergence = p->value_logical - ("Method:star_maker:use_velocity_divergence",true); + ("Method:star_maker:use_velocity_divergence",false); method_star_maker_use_dynamical_time = p->value_logical - ("Method:star_maker:use_dynamical_time",true); + ("Method:star_maker:use_dynamical_time",false); + + method_star_maker_use_cooling_time = p->value_logical + ("Method:star_maker:use_cooling_time",false); method_star_maker_use_self_gravitating = p->value_logical ("Method:star_maker:use_self_gravitating", false); @@ -1497,9 +1610,24 @@ void EnzoConfig::read_method_star_maker_(Parameters * p) method_star_maker_use_jeans_mass = p->value_logical ("Method:star_maker:use_jeans_mass", false); + method_star_maker_use_temperature_threshold = p->value_logical + ("Method:star_maker:use_temperature_threshold",false); + + method_star_maker_use_critical_metallicity = p->value_logical + ("Method:star_maker:use_critical_metallicity",false); + method_star_maker_number_density_threshold = p->value_float ("Method:star_maker:number_density_threshold",0.0); + method_star_maker_overdensity_threshold = p->value_float + ("Method:star_maker:overdensity_threshold",0.0); + + method_star_maker_temperature_threshold = p->value_float + ("Method:star_maker:temperature_threshold",1.0E4); + + method_star_maker_critical_metallicity = p->value_float + ("Method:star_maker:critical_metallicity",0.0); + method_star_maker_maximum_mass_fraction = p->value_float ("Method:star_maker:maximum_mass_fraction",0.5); @@ -1511,6 +1639,12 @@ void EnzoConfig::read_method_star_maker_(Parameters * p) method_star_maker_maximum_star_mass = p->value_float ("Method:star_maker:maximum_star_mass",1.0E4); + + method_star_maker_min_level = p->value_integer + ("Method:star_maker:min_level",0); + + method_star_maker_turn_off_probability = p->value_logical + ("Method:star_maker:turn_off_probability",false); } diff --git a/src/Enzo/enzo_EnzoConfig.hpp b/src/Enzo/enzo_EnzoConfig.hpp index 2d5c36b026..5cad96c630 100644 --- a/src/Enzo/enzo_EnzoConfig.hpp +++ b/src/Enzo/enzo_EnzoConfig.hpp @@ -39,6 +39,8 @@ inline void operator|(PUP::er &p, chemistry_data &c){ } p | c.cmb_temperature_floor; + p | c.h2_charge_exchange_rate; + p | c.h2_h_cooling_rate; p | c.Gamma; p | c.h2_on_dust; p | c.use_dust_density_field; @@ -181,6 +183,12 @@ class EnzoConfig : public Config { #endif /* CONFIG_USE_GRACKLE */ // EnzoInitialFeedbackTest initial_feedback_test_density(), + initial_feedback_test_HI_density(), + initial_feedback_test_HII_density(), + initial_feedback_test_HeI_density(), + initial_feedback_test_HeII_density(), + initial_feedback_test_HeIII_density(), + initial_feedback_test_e_density(), initial_feedback_test_star_mass(), initial_feedback_test_temperature(), initial_feedback_test_from_file(), @@ -313,6 +321,7 @@ class EnzoConfig : public Config { method_hydro_reconstruct_positive(false), method_hydro_riemann_solver(""), /// EnzoMethodFeedback + method_feedback_flavor(""), method_feedback_ejecta_mass(0.0), method_feedback_supernova_energy(1.0), method_feedback_ejecta_metal_fraction(0.0), @@ -322,19 +331,37 @@ class EnzoConfig : public Config { method_feedback_ke_fraction(0.0), method_feedback_use_ionization_feedback(false), method_feedback_time_first_sn(-1.0), // in Myr + /// EnzoMethodFeedbackSTARSS + method_feedback_single_sn(0), + method_feedback_unrestricted_sn(0), + method_feedback_stellar_winds(0), + method_feedback_gas_return_fraction(0.0), + method_feedback_min_level(0), + method_feedback_analytic_SNR_shell_mass(0), + method_feedback_fade_SNR(0), + method_feedback_NEvents(0), /// EnzoMethodStarMaker method_star_maker_flavor(""), - method_star_maker_use_density_threshold(true), // check above density threshold before SF - method_star_maker_use_velocity_divergence(true), // check for converging flow before SF - method_star_maker_use_dynamical_time(true), // + method_star_maker_use_density_threshold(false), // check above density threshold before SF + method_star_maker_use_velocity_divergence(false), // check for converging flow before SF + method_star_maker_use_dynamical_time(false), // + method_star_maker_use_altAlpha(false), // alternate virial parameter calculation + method_star_maker_use_cooling_time(false), method_star_maker_use_self_gravitating(false), // method_star_maker_use_h2_self_shielding(false), method_star_maker_use_jeans_mass(false), + method_star_maker_use_overdensity_threshold(false), + method_star_maker_use_critical_metallicity(false), + method_star_maker_use_temperature_threshold(false), + method_star_maker_critical_metallicity(0.0), + method_star_maker_temperature_threshold(1.0E4), method_star_maker_number_density_threshold(0.0), // Number density threshold in cgs method_star_maker_maximum_mass_fraction(0.5), // maximum cell mass fraction to convert to stars method_star_maker_efficiency(0.01), // star maker efficiency method_star_maker_minimum_star_mass(1.0E4), // minium star particle mass in solar masses method_star_maker_maximum_star_mass(1.0E4), // maximum star particle mass in solar masses + method_star_maker_min_level(0), // minimum refinement level for star formation + method_star_maker_turn_off_probability(false), // EnzoMethodTurbulence method_turbulence_edot(0.0), method_turbulence_mach_number(0.0), @@ -672,6 +699,12 @@ class EnzoConfig : public Config { double initial_feedback_test_position[3]; double initial_feedback_test_density; + double initial_feedback_test_HI_density; + double initial_feedback_test_HII_density; + double initial_feedback_test_HeI_density; + double initial_feedback_test_HeII_density; + double initial_feedback_test_HeIII_density; + double initial_feedback_test_e_density; double initial_feedback_test_star_mass; double initial_feedback_test_temperature; bool initial_feedback_test_from_file; @@ -732,6 +765,7 @@ class EnzoConfig : public Config { /// EnzoMethodFeedback + std::string method_feedback_flavor; double method_feedback_ejecta_mass; double method_feedback_supernova_energy; double method_feedback_ejecta_metal_fraction; @@ -742,20 +776,40 @@ class EnzoConfig : public Config { bool method_feedback_shift_cell_center; bool method_feedback_use_ionization_feedback; + /// EnzoMethodFeedbackSTARSS + + int method_feedback_single_sn; + int method_feedback_unrestricted_sn; + int method_feedback_stellar_winds; + double method_feedback_gas_return_fraction; + int method_feedback_min_level; + int method_feedback_analytic_SNR_shell_mass; + int method_feedback_fade_SNR; + int method_feedback_NEvents; /// EnzoMethodStarMaker std::string method_star_maker_flavor; + bool method_star_maker_use_altAlpha; bool method_star_maker_use_density_threshold; + bool method_star_maker_use_overdensity_threshold; + bool method_star_maker_use_temperature_threshold; + bool method_star_maker_use_critical_metallicity; bool method_star_maker_use_velocity_divergence; + bool method_star_maker_use_cooling_time; bool method_star_maker_use_dynamical_time; bool method_star_maker_use_h2_self_shielding; bool method_star_maker_use_jeans_mass; bool method_star_maker_use_self_gravitating; double method_star_maker_number_density_threshold; + double method_star_maker_overdensity_threshold; + double method_star_maker_temperature_threshold; + double method_star_maker_critical_metallicity; double method_star_maker_maximum_mass_fraction; double method_star_maker_efficiency; double method_star_maker_minimum_star_mass; double method_star_maker_maximum_star_mass; + int method_star_maker_min_level; + bool method_star_maker_turn_off_probability; /// EnzoMethodTurbulence double method_turbulence_edot; diff --git a/src/Enzo/enzo_EnzoInitialBurkertBodenheimer.cpp b/src/Enzo/enzo_EnzoInitialBurkertBodenheimer.cpp index 953a81b6e8..d903460aa1 100644 --- a/src/Enzo/enzo_EnzoInitialBurkertBodenheimer.cpp +++ b/src/Enzo/enzo_EnzoInitialBurkertBodenheimer.cpp @@ -255,11 +255,11 @@ void EnzoInitialBurkertBodenheimer::enforce_block if (metal) metal[i] = d[i]*inner_metal_fraction; - if(i == 20) { - CkPrintf("Density = %e\n", d[i]); - CkPrintf("Angular Velocity = %e rad/s", sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i])/radius / enzo_units->time()); + //if(i == 20) { + //CkPrintf("Density = %e\n", d[i]); + //CkPrintf("Angular Velocity = %e rad/s", sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i])/radius / enzo_units->time()); //CkExit(-99); - } + //} } } } diff --git a/src/Enzo/enzo_EnzoInitialFeedbackTest.cpp b/src/Enzo/enzo_EnzoInitialFeedbackTest.cpp index 1472c80800..32bcccfae8 100644 --- a/src/Enzo/enzo_EnzoInitialFeedbackTest.cpp +++ b/src/Enzo/enzo_EnzoInitialFeedbackTest.cpp @@ -108,6 +108,14 @@ void EnzoInitialFeedbackTest::enforce_block enzo_float * metal = (enzo_float *) field.values("metal_density"); + enzo_float * d_HI = (enzo_float *) field.values("HI_density"); + enzo_float * d_HII = (enzo_float *) field.values("HII_density"); + enzo_float * d_HeI = (enzo_float *) field.values("HeI_density"); + enzo_float * d_HeII = (enzo_float *) field.values("HeII_density"); + enzo_float * d_HeIII = (enzo_float *) field.values("HeIII_density"); + enzo_float * d_electron = (enzo_float *) field.values("e_density"); + enzo_float * temperature = (enzo_float *) field.values("temperature"); + // Block size (excluding ghosts) int nx,ny,nz; field.size(&nx,&ny,&nz); @@ -129,14 +137,22 @@ void EnzoInitialFeedbackTest::enforce_block int ngx = nx + 2*gx; int ngy = ny + 2*gy; int ngz = nz + 2*gz; - + for (int iz = 0; iz < ngz; iz++){ for (int iy = 0; iy < ngy; iy++){ for (int ix = 0; ix < ngx; ix++){ int i = INDEX(ix,iy,iz,ngx,ngy); - d[i] = enzo_config->initial_feedback_test_density / enzo_units->density(); + // values are specified in CGS in the parameter file + d[i] = enzo_config->initial_feedback_test_density / enzo_units->density(); + d_HI[i] = enzo_config->initial_feedback_test_HI_density / enzo_units->density(); + d_HII[i] = enzo_config->initial_feedback_test_HII_density / enzo_units->density(); + d_HeI[i] = enzo_config->initial_feedback_test_HeI_density / enzo_units->density(); + d_HeII[i] = enzo_config->initial_feedback_test_HeII_density / enzo_units->density(); + d_HeIII[i] = enzo_config->initial_feedback_test_HeIII_density / enzo_units->density(); + d_electron[i] = enzo_config->initial_feedback_test_e_density / enzo_units->density(); + temperature[i] = enzo_config->initial_feedback_test_temperature / enzo_units->kelvin_per_energy_units(); for (int dim = 0; dim < 3; dim++) v3[dim][i] = 0.0; @@ -156,8 +172,6 @@ void EnzoInitialFeedbackTest::enforce_block // drop in a particle - if (enzo_block->level() != 0) return; // don't do particles below root grid - ParticleDescr * particle_descr = cello::particle_descr(); Particle particle = block->data()->particle(); @@ -170,7 +184,7 @@ void EnzoInitialFeedbackTest::enforce_block int ia_vx = particle.attribute_index (it, "vx"); int ia_vy = particle.attribute_index (it, "vy"); int ia_vz = particle.attribute_index (it, "vz"); - int ia_loc = particle.attribute_index (it, "is_local"); + int ia_cop = particle.attribute_index (it, "is_copy"); int ia_id = particle.attribute_index (it, "id"); int ia_to = particle.has_attribute(it,"creation_time") ? @@ -197,7 +211,7 @@ void EnzoInitialFeedbackTest::enforce_block enzo_float * pmetal = 0; enzo_float * plifetime = 0; enzo_float * pform = 0; - int64_t * is_local = 0; + int64_t * is_copy = 0; int64_t * id = 0; @@ -226,7 +240,7 @@ void EnzoInitialFeedbackTest::enforce_block pmetal = (enzo_float *) particle.attribute_array(it, ia_metal, ib); plifetime = (enzo_float *) particle.attribute_array(it, ia_l, ib); pform = (enzo_float *) particle.attribute_array(it, ia_to, ib); - is_local = (int64_t *) particle.attribute_array(it, ia_loc, ib); + is_copy = (int64_t *) particle.attribute_array(it, ia_cop, ib); ipp = 0; for (int i = 0; i < this->num_particles; i++){ @@ -251,7 +265,7 @@ void EnzoInitialFeedbackTest::enforce_block plifetime[ipp] = 1.00E9* enzo_constants::yr_s / enzo_units->time(); pform[ipp] = 1.0E-10 * enzo_constants::yr_s / enzo_units->time(); // really just needs to be non-zero - is_local[ipp] = 1; + is_copy[ipp] = 1; ipp++; } diff --git a/src/Enzo/enzo_EnzoMethodFeedbackSTARSS.cpp b/src/Enzo/enzo_EnzoMethodFeedbackSTARSS.cpp new file mode 100644 index 0000000000..d3db331858 --- /dev/null +++ b/src/Enzo/enzo_EnzoMethodFeedbackSTARSS.cpp @@ -0,0 +1,1527 @@ + +/// See LICENSE_CELLO file for license and copyright information + +/// @file enzo_EnzoMethodFeedbackSTARSS.cpp +/// @author Will Hicks (whicks@ucsd.edu) +// Andrew Emerick (aemerick11@gmail.com) +// +/// @date +/// @brief Implements the STARSS model for stellar feedback +/// as implemented into Enzo as the MechStars model +/// by Azton Wells. This is a direct copy/paste of the +/// MechStars methods in Enzo, adapted for Enzo-E. + + +#include "cello.hpp" +#include "enzo.hpp" + +#include + +//#ifdef NOTDEFINED // for now... since not done coding + +#define DEBUG_FEEDBACK_STARSS +//#define DEBUG_FEEDBACK_STARSS_SN + +// ============================================================================= +// splice these off to a different file (later) +// TODO: Maybe create EnzoStarParticle class and add rate-calculating functions there? + +int EnzoMethodFeedbackSTARSS::determineSN(double age_Myr, int* nSNII, int* nSNIA, + double pmass_Msun, double tunit, float dt){ + + const EnzoConfig * enzo_config = enzo::config(); + std::mt19937 mt(std::time(nullptr)); + + if (NEvents > 0){ + *nSNII = 1; + NEvents -= 1; + return 1; + } + /* else, calculate SN rate, probability and determine number of events */ + *nSNII = 0; + *nSNIA = 0; + double RII=0, RIA=0, PII=0, PIA=0, random = 0; + if (enzo_config->method_feedback_single_sn && NEvents < 0) + { + /* age-dependent rates */ + if (age_Myr < 3.401) + { + RII = 0.0; + RIA = 0.0; + } + if (3.401 <= age_Myr && age_Myr < 10.37) + { + RII = 5.408e-4; + RIA = 0.0; + } + if (10.37 <= age_Myr && age_Myr < 37.53) + { + RII = 2.516e-4; + RIA = 0.0; + } + if (37.53 <= age_Myr) + { + RII = 0.0; + RIA = 5.3e-8+1.6e-5*exp(-0.5*pow((age_Myr-50.0)/10.0, 2)); + } + /* rates -> probabilities */ + if (RII > 0){ + PII = RII * pmass_Msun / enzo_constants::Myr_s *tunit*dt; + double random = double(mt())/double(mt.max()); + if (PII > 1.0){ + if (enzo_config->method_feedback_unrestricted_sn) { + int round = (int)PII; + *nSNII = round; // number of Type II SNe that will definitely go off + PII -= round; // probability for setting off one more supernova + } + else { + // Restrict particles to at most one supernova per timestep + // if unrestricted_sn = false + *nSNII = 1; + PII = 0; // set probability for setting off another supernova to zero + } + } + + if (random <= PII){ + *nSNII += 1; + } + } + + #ifdef DEBUG_FEEDBACK_STARSS_SN + CkPrintf("MethodFeedbackSTARSS::determineSN() -- pmass_Msun = %f; age = %f; RII = %f; RIA = %f\n", + pmass_Msun, age_Myr, RII, RIA); + #endif + + if (RIA > 0){ + PIA = RIA*pmass_Msun / enzo_constants::Myr_s *tunit*dt; + float random = float(rand())/float(RAND_MAX); + + if (PIA > 1.0) + { + if (enzo_config->method_feedback_unrestricted_sn) { + int round = int(PIA); + *nSNIA = round; // number of Type Ia SNe that will definitely go off + PIA -= round; // probability for setting off one more supernova + } + else { + *nSNIA = 1; + PIA = 0; + } + } + + if (!enzo_config->method_feedback_unrestricted_sn && *nSNII > 0) { + // for unrestricted_sn = false, only set off a Type IA supernova + // if nSNII == 0 + *nSNIA = 0; + PIA = 0; + } + + if (random < PIA){ + *nSNIA += 1; + } + } + + #ifdef DEBUG_FEEDBACK_STARSS_SN + CkPrintf("MethodFeedbackSTARSS::determineSN() -- PII = %f; PIA = %f\n", PII, PIA); + #endif + } + + return 1; +} + +int EnzoMethodFeedbackSTARSS::determineWinds(double age_Myr, double * eWinds, double * mWinds, double * zWinds, + double pmass_Msun, double metallicity_Zsun, double tunit, double dt) + { + // age in Myr + + const EnzoConfig * enzo_config = enzo::config(); + bool oldEnough = (age_Myr < 0.0001)?(false):(true); + double windE = 0, wind_mass_solar = 0, windZ = 0.0; + double wind_factor = 0.0; + double e_factor = 0.0; + + if (pmass_Msun > 11 && oldEnough){ + + if (0.001 < age_Myr && age_Myr < 1.0){ + wind_factor = 4.763 * std::min((0.01 + metallicity_Zsun), 1.0); + } + if (1 <= age_Myr && age_Myr < 3.5){ + wind_factor = 4.763*std::min(0.01+metallicity_Zsun, 1.0)* + pow(age_Myr, 1.45+0.08*std::min(log(metallicity_Zsun), 1.0)); + } + if (3.5 <= age_Myr && age_Myr < 100){ + wind_factor = 29.4*pow(age_Myr/3.5, -3.25)+0.0042; + } + if (age_Myr < 100){ + double d = powl(1+age_Myr/2.5, 1.4); + double a50 = powl(double(age_Myr)/10.0, 5.0); + e_factor = 5.94e4 / d + a50 +4.83; + } + if (100 <= age_Myr){ + e_factor = 4.83; + wind_factor = 0.42*pow(age_Myr/1000, -1.1)/(19.81/log(age_Myr)); + } + wind_mass_solar = pmass_Msun * wind_factor; // Msun/Gyr + wind_mass_solar = wind_mass_solar*dt*tunit/(1e3 * enzo_constants::Myr_s); // Msun + + if (wind_mass_solar > pmass_Msun){ + CkPrintf("Winds too large Mw = %e, Mp = %e age=%f, Z = %e\n", + wind_mass_solar, pmass_Msun, age_Myr, metallicity_Zsun); + wind_mass_solar = 0.125*pmass_Msun; // limit loss to huge if necessary. + } + windZ = std::max(enzo_constants::metallicity_solar, 0.016+0.0041*std::max(metallicity_Zsun, 1.65)+0.0118)*wind_mass_solar; + windE = e_factor * 1e12 * wind_mass_solar; + + + *mWinds = wind_mass_solar; + *zWinds = windZ; + *eWinds = windE; + +} + + return 1; + +} + +void EnzoMethodFeedbackSTARSS::transformComovingWithStar(enzo_float * density, + enzo_float * velocity_x, enzo_float * velocity_y, enzo_float * velocity_z, + const enzo_float up, const enzo_float vp, const enzo_float wp, + const int mx, const int my, const int mz, int direction) const throw() +{ + int size = mx*my*mz; + if (direction > 0) + { + // to comoving with star + // NOTE: This transforms the velocity field into a momentum density + // field for the sake of depositing momentum easily + + for (int ind = 0; indrefresh_set_name(ir_post_,name()); + Refresh * refresh = cello::refresh(ir_post_); + refresh->add_all_fields(); + + sf_minimum_level_ = enzo_config->method_star_maker_min_level; + single_sn_ = enzo_config->method_feedback_single_sn; + + // Initialize temporary fields + i_d_dep = cello::field_descr()->insert_temporary(); + i_te_dep = cello::field_descr()->insert_temporary(); + i_ge_dep = cello::field_descr()->insert_temporary(); + i_mf_dep = cello::field_descr()->insert_temporary(); + i_vx_dep = cello::field_descr()->insert_temporary(); + i_vy_dep = cello::field_descr()->insert_temporary(); + i_vz_dep = cello::field_descr()->insert_temporary(); + i_d_shell= cello::field_descr()->insert_temporary(); + + i_d_dep_a = cello::field_descr()->insert_temporary(); + i_te_dep_a = cello::field_descr()->insert_temporary(); + i_ge_dep_a = cello::field_descr()->insert_temporary(); + i_mf_dep_a = cello::field_descr()->insert_temporary(); + i_vx_dep_a = cello::field_descr()->insert_temporary(); + i_vy_dep_a = cello::field_descr()->insert_temporary(); + i_vz_dep_a = cello::field_descr()->insert_temporary(); + i_d_shell_a= cello::field_descr()->insert_temporary(); + + // Deposition across grid boundaries is handled using refresh with set_accumulate=true. + // The set_accumulate flag tells Cello to include ghost zones in the refresh operation, + // and adds the ghost zone values from the "src" field to the corresponding active zone + // values in the "dst" field. + + // I'm currently using a set of two initially empty fields for each field that SNe directly affect. + // CIC deposition directly modifies the *_deposit fields (including ghost zones). The ghost zone + // values are then sent to the *deposit_accumulate fields during the refresh operation. Values are then + // copied back to the original field. + + ir_feedback_ = add_refresh_(); + cello::simulation()->refresh_set_name(ir_feedback_,name()+":add"); + Refresh * refresh_fb = cello::refresh(ir_feedback_); + + refresh_fb->set_accumulate(true); + + refresh_fb->add_field_src_dst + (i_d_dep, i_d_dep_a); + refresh_fb->add_field_src_dst + (i_te_dep, i_te_dep_a); + refresh_fb->add_field_src_dst + (i_ge_dep, i_ge_dep_a); + refresh_fb->add_field_src_dst + (i_mf_dep, i_mf_dep_a); + refresh_fb->add_field_src_dst + (i_vx_dep, i_vx_dep_a); + refresh_fb->add_field_src_dst + (i_vy_dep, i_vy_dep_a); + refresh_fb->add_field_src_dst + (i_vz_dep, i_vz_dep_a); + refresh_fb->add_field_src_dst + (i_d_shell, i_d_shell_a); + + refresh_fb->set_callback(CkIndex_EnzoBlock::p_method_feedback_starss_end()); + + std::mt19937 mt(std::time(nullptr)); // initialize random function + + // initialize NEvents parameter (mainly for testing). Sets off 'NEvents' supernovae, + // with at most one supernova per star particle per cycle. + this->NEvents = enzo_config->method_feedback_NEvents; + return; +} + +void EnzoMethodFeedbackSTARSS::pup (PUP::er &p) +{ + /// NOTE: Change this function whenever attributes change + + TRACEPUP; + + Method::pup(p); + + p | sf_minimum_level_; + p | single_sn_; + p | NEvents; + p | ir_feedback_; + + p | i_d_dep; + p | i_te_dep; + p | i_ge_dep; + p | i_mf_dep; + p | i_vx_dep; + p | i_vy_dep; + p | i_vz_dep; + p | i_d_shell; + + p | i_d_dep_a; + p | i_te_dep_a; + p | i_ge_dep_a; + p | i_mf_dep_a; + p | i_vx_dep_a; + p | i_vy_dep_a; + p | i_vz_dep_a; + p | i_d_shell_a; + + return; +} + +void EnzoMethodFeedbackSTARSS::compute (Block * block) throw() +{ + + if (block->is_leaf()){ + this->compute_(block); + } + + else { + block->compute_done(); + } + return; +} + +void EnzoMethodFeedbackSTARSS::add_accumulate_fields(EnzoBlock * enzo_block) throw() +{ + int mx, my, mz, gx, gy, gz, nx, ny, nz; + + double xm, ym, zm, xp, yp, zp, hx, hy, hz; + Field field = enzo_block->data()->field(); + field.size(&nx,&ny,&nz); + field.ghost_depth(0,&gx,&gy,&gz); + enzo_block->data()->lower(&xm,&ym,&zm); + enzo_block->data()->upper(&xp,&yp,&zp); + field.cell_width(xm,xp,&hx,ym,yp,&hy,zm,zp,&hz); + + mx = nx + 2*gx; + my = ny + 2*gy; + mz = nz + 2*gz; + + // add accumulated values over and reset them to zero + + enzo_float * d = (enzo_float *) field.values("density"); + enzo_float * te = (enzo_float *) field.values("total_energy"); + enzo_float * ge = (enzo_float *) field.values("internal_energy"); + enzo_float * mf = (enzo_float *) field.values("metal_density"); + enzo_float * vx = (enzo_float *) field.values("velocity_x"); + enzo_float * vy = (enzo_float *) field.values("velocity_y"); + enzo_float * vz = (enzo_float *) field.values("velocity_z"); + + enzo_float * d_dep = (enzo_float *) field.values(i_d_dep); + enzo_float * te_dep = (enzo_float *) field.values(i_te_dep); + enzo_float * ge_dep = (enzo_float *) field.values(i_ge_dep); + enzo_float * mf_dep = (enzo_float *) field.values(i_mf_dep); + enzo_float * vx_dep = (enzo_float *) field.values(i_vx_dep); + enzo_float * vy_dep = (enzo_float *) field.values(i_vy_dep); + enzo_float * vz_dep = (enzo_float *) field.values(i_vz_dep); + + enzo_float * d_dep_a = (enzo_float *) field.values(i_d_dep_a); + enzo_float * te_dep_a = (enzo_float *) field.values(i_te_dep_a); + enzo_float * ge_dep_a = (enzo_float *) field.values(i_ge_dep_a); + enzo_float * mf_dep_a = (enzo_float *) field.values(i_mf_dep_a); + enzo_float * vx_dep_a = (enzo_float *) field.values(i_vx_dep_a); + enzo_float * vy_dep_a = (enzo_float *) field.values(i_vy_dep_a); + enzo_float * vz_dep_a = (enzo_float *) field.values(i_vz_dep_a); + + enzo_float * d_shell = (enzo_float *) field.values(i_d_shell); + enzo_float * d_shell_a = (enzo_float *) field.values(i_d_shell_a); + + EnzoUnits * enzo_units = enzo::units(); + double cell_volume_code = hx*hy*hz; + double cell_volume_cgs = cell_volume_code * enzo_units->volume(); + double rhounit = enzo_units->density(); + + double maxEvacFraction = 0.75; // TODO: make this a parameter + + // multiply by this value to convert cell density (in code units) + // to cell mass in Msun + double rho_to_m = rhounit*cell_volume_cgs / enzo_constants::mass_solar; + + for (int iz=gz; iz (this->method()); + method->add_accumulate_fields(this); + + compute_done(); + return; +} + +void EnzoMethodFeedbackSTARSS::compute_ (Block * block) +{ + + //---------------------------------------------------- + // some constants here that might get moved to parameters or something else + const float SNII_ejecta_mass_Msun = 10.5; + const float SNIa_ejecta_mass_Msun = 1.4; + const float z_solar = enzo_constants::metallicity_solar; //Solar metal fraction (0.012) + //----------------------------------------------------- + + EnzoBlock * enzo_block = enzo::block(block); + const EnzoConfig * enzo_config = enzo::config(); + Particle particle = enzo_block->data()->particle(); + EnzoUnits * enzo_units = enzo::units(); + + double munit = enzo_units->mass(); + double tunit = enzo_units->time(); + + double current_time = block->time(); + + Field field = enzo_block->data()->field(); + + enzo_float * d = (enzo_float *) field.values("density"); + enzo_float * te = (enzo_float *) field.values("total_energy"); + enzo_float * ge = (enzo_float *) field.values("internal_energy"); + enzo_float * mf = (enzo_float *) field.values("metal_density"); + + // Obtain grid sizes and ghost sizes + int mx, my, mz, gx, gy, gz, nx, ny, nz; + double xm, ym, zm, xp, yp, zp, hx, hy, hz; + field.size(&nx,&ny,&nz); + field.ghost_depth(0,&gx,&gy,&gz); + block->data()->lower(&xm,&ym,&zm); + block->data()->upper(&xp,&yp,&zp); + field.cell_width(xm,xp,&hx,ym,yp,&hy,zm,zp,&hz); + + mx = nx + 2*gx; + my = ny + 2*gy; + mz = nz + 2*gz; + + const int rank = cello::rank(); + + // apply feedback depending on particle type + // for now, just do this for all star particles + + int it = particle.type_index("star"); + + // AE: Copying this from MechStars_FeedbackRoutine + + int numSN = 0; // counter of SN events + int count = 0; // counter of particles + + allocate_temporary_(enzo_block); + + // initialize temporary fields as zero + enzo_float * d_dep = (enzo_float *) field.values(i_d_dep); + enzo_float * te_dep = (enzo_float *) field.values(i_te_dep); + enzo_float * ge_dep = (enzo_float *) field.values(i_ge_dep); + enzo_float * mf_dep = (enzo_float *) field.values(i_mf_dep); + enzo_float * vx_dep = (enzo_float *) field.values(i_vx_dep); + enzo_float * vy_dep = (enzo_float *) field.values(i_vy_dep); + enzo_float * vz_dep = (enzo_float *) field.values(i_vz_dep); + + enzo_float * d_dep_a = (enzo_float *) field.values(i_d_dep_a); + enzo_float * te_dep_a = (enzo_float *) field.values(i_te_dep_a); + enzo_float * ge_dep_a = (enzo_float *) field.values(i_ge_dep_a); + enzo_float * mf_dep_a = (enzo_float *) field.values(i_mf_dep_a); + enzo_float * vx_dep_a = (enzo_float *) field.values(i_vx_dep_a); + enzo_float * vy_dep_a = (enzo_float *) field.values(i_vy_dep_a); + enzo_float * vz_dep_a = (enzo_float *) field.values(i_vz_dep_a); + + enzo_float * d_shell = (enzo_float *) field.values(i_d_shell); + enzo_float * d_shell_a = (enzo_float *) field.values(i_d_shell_a); + + + for (int i=0; i 0.0 && plifetime[ipdl] > 0.0){ + const double age = (current_time - pcreation[ipdc]) * enzo_units->time() / enzo_constants::Myr_s; + count++; // increment particles examined here + + // compute coordinates of central feedback cell + // this must account for ghost zones + + int ix = (int) floor((px[ipdp] - xm) / hx + gx); + int iy = (int) floor((py[ipdp] - ym) / hy + gy); + int iz = (int) floor((pz[ipdp] - zm) / hz + gz); + + int index = INDEX(ix,iy,iz,mx,my); // 1D cell index of star position + + int nSNII = 0, nSNIa = 0; + double SNMassEjected = 0.0, SNMetalEjected = 0.0; + + if (single_sn_){ + + /* Determine number of SN events from rates (currently taken from Hopkins 2018) */ + + determineSN(age, &nSNII, &nSNIa, pmass_solar, + tunit, block->dt()); + + numSN += nSNII + nSNIa; + + #ifdef DEBUG_FEEDBACK_STARSS + if (nSNII + nSNIa > 0){ + CkPrintf("MethodFeedbackSTARSS SNe %d %d level = %d age = %f\n", nSNII, nSNIa, block->level(), age); + } + #endif + + if (nSNII+nSNIa > 0){ + /* set feedback properties based on number and types of SN */ + double energySN = (nSNII+nSNIa) * 1.0e51; + + SNMassEjected = SNII_ejecta_mass_Msun * nSNII + + SNIa_ejecta_mass_Msun * nSNIa; // AE: split this in two channels for yields + + const double starZ = pmetal[ipdmf] / z_solar; + + /* Fixed mass ejecta */ + + this->deposit_feedback( block, energySN, SNMassEjected, SNMetalEjected, + pvx[ipdv],pvy[ipdv],pvz[ipdv], + px[ipdp],py[ipdp],pz[ipdp], + ix, iy, iz, 0, nSNII, nSNIa, starZ); // removed P3 + + // add counter for number of SNe for the particle + psncounter[ipsn] += nSNII+nSNIa; + + } // if nSNII or nSNIa > 0 + } // if single SN + + // Now do stellar winds + double windMass=0.0, windMetals=0.0, windEnergy=0.0; + if (enzo_config->method_feedback_stellar_winds){ + + const double starZ = pmetal[ipdmf] / z_solar; + + determineWinds(age, &windEnergy, &windMass, &windMetals, + pmass_solar, + starZ, tunit, block->dt()); + + if (windMass > 0){ + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Adding stellar winds...\n"); + #endif + this->deposit_feedback( block, windEnergy, windMass, windMetals, + pvx[ipdv],pvy[ipdv],pvz[ipdv], + px[ipdp],py[ipdp],pz[ipdp], + ix, iy, iz, 1, 0, 0, starZ); // removed P3 + + } // if wind mass > 0 + } // if winds + + // windMass and SNMassEjected are in units of Msun + pmass[ipdm] -= std::max(0.0, + (windMass + SNMassEjected) / + (munit/enzo_constants::mass_solar)); + + + } // if mass and lifetime > 0 + + } // end particle loop + } // end batch loop + + + if (count > 0){ + CkPrintf("FeedbackSTARSS: Num FB particles = %d Events = %d FeedbackTime %e\n", + count, numSN, 0.00); + } + + // refresh + cello::refresh(ir_feedback_)->set_active(enzo_block->is_leaf()); + enzo_block->refresh_start(ir_feedback_, CkIndex_EnzoBlock::p_method_feedback_starss_end()); + return; +} + + + + +// ---------------------------------------------------------------------------- + + +void EnzoMethodFeedbackSTARSS::deposit_feedback (Block * block, + double ejectaEnergy, + double ejectaMass, + double ejectaMetals, + const enzo_float up, const enzo_float vp, const enzo_float wp, + const enzo_float xp, const enzo_float yp, const enzo_float zp, + const int ix, const int iy, const int iz, + const int winds, const int nSNII, + const int nSNIA, + const double starZ) + + const throw(){ + /* + This routine will create a cube of coupling particles, where we determine + the feedback quantities. The vertices of the cube are ~coupled particles + Each vertex particle will then be CIC deposited to the grid! + */ + EnzoBlock * enzo_block = enzo::block(block); + Field field = enzo_block->data()->field(); + EnzoUnits * enzo_units = enzo::units(); + double vunit = enzo_units->velocity(); + double rhounit = enzo_units->density(); + double lunit = enzo_units->length(); + double munit = enzo_units->mass(); + double tunit = enzo_units->time(); + double eunit = munit*vunit*vunit; // energy units (NOTE: not specific energy) + + const EnzoConfig * enzo_config = enzo::config(); + + bool AnalyticSNRShellMass = enzo_config->method_feedback_analytic_SNR_shell_mass; + + // Obtain grid sizes and ghost sizes + int mx, my, mz, gx, gy, gz, nx, ny, nz; + double xm, ym, zm, xp_, yp_, zp_, hx, hy, hz; + field.size(&nx,&ny,&nz); + field.ghost_depth(0,&gx,&gy,&gz); + enzo_block->data()->lower(&xm,&ym,&zm); + enzo_block->data()->upper(&xp_,&yp_,&zp_); + field.cell_width(xm,xp_,&hx,ym,yp_,&hy,zm,zp_,&hz); + + mx = nx + 2*gx; + my = ny + 2*gy; + mz = nz + 2*gz; + + int size = mx*my*mz; + + const int rank = cello::rank(); + + double cell_volume_code = hx*hy*hz; + double cell_volume = cell_volume_code*enzo_units->volume(); + + // conversion from code_density to mass in Msun + double rho_to_m = rhounit*cell_volume / enzo_constants::mass_solar; + + enzo_float * d = (enzo_float *) field.values("density"); + enzo_float * te = (enzo_float *) field.values("total_energy"); + enzo_float * ge = (enzo_float *) field.values("internal_energy"); + + enzo_float * vx = (enzo_float *) field.values("velocity_x"); + enzo_float * vy = (enzo_float *) field.values("velocity_y"); + enzo_float * vz = (enzo_float *) field.values("velocity_z"); + + enzo_float * mf = (enzo_float *) field.values("metal_density"); + + enzo_float * temperature = (enzo_float *) field.values("temperature"); + + enzo_float * dHI = field.is_field("HI_density") ? + (enzo_float*) field.values("HI_density") : NULL; + enzo_float * dHII = field.is_field("HII_density") ? + (enzo_float*) field.values("HII_density") : NULL; + enzo_float * dHeI = field.is_field("HeI_density") ? + (enzo_float*) field.values("HeI_density") : NULL; + enzo_float * dHeII = field.is_field("HeII_density") ? + (enzo_float*) field.values("HeII_density") : NULL; + enzo_float * dHeIII = field.is_field("HeIII_density") ? + (enzo_float*) field.values("HeIII_density") : NULL; + enzo_float * d_el = field.is_field("e_density") ? + (enzo_float*) field.values("e_density") : NULL; + + enzo_float * dH2I = field.is_field("H2I_density") ? + (enzo_float*) field.values("H2I_density") : NULL; + enzo_float * dH2II = field.is_field("H2II_density") ? + (enzo_float*) field.values("H2II_density") : NULL; + enzo_float * dHM = field.is_field("HM_density") ? + (enzo_float*) field.values("HM_density") : NULL; + + enzo_float * dDI = field.is_field("DI_density") ? + (enzo_float *) field.values("DI_density") : NULL; + enzo_float * dDII = field.is_field("DII_density") ? + (enzo_float *) field.values("DII_density") : NULL; + enzo_float * dHDI = field.is_field("HDI_density") ? + (enzo_float *) field.values("HDI_density") : NULL; + + // "deposit" fields that hold running total for all exploding particles in this block + enzo_float * d_dep_tot = (enzo_float *) field.values(i_d_dep); + enzo_float * te_dep_tot = (enzo_float *) field.values(i_te_dep); + enzo_float * ge_dep_tot = (enzo_float *) field.values(i_ge_dep); + enzo_float * mf_dep_tot = (enzo_float *) field.values(i_mf_dep); + enzo_float * vx_dep_tot = (enzo_float *) field.values(i_vx_dep); + enzo_float * vy_dep_tot = (enzo_float *) field.values(i_vy_dep); + enzo_float * vz_dep_tot = (enzo_float *) field.values(i_vz_dep); + + // holds just shell densities (used for refresh+accumulate) + enzo_float * d_shell = (enzo_float *) field.values(i_d_shell); + + // allocate another set of temporary deposit fields for this event + enzo_float * d_dep = new enzo_float[size]; + enzo_float * te_dep = new enzo_float[size]; + enzo_float * ge_dep = new enzo_float[size]; + enzo_float * mf_dep = new enzo_float[size]; + enzo_float * vx_dep = new enzo_float[size]; + enzo_float * vy_dep = new enzo_float[size]; + enzo_float * vz_dep = new enzo_float[size]; + + // initialize temporary deposit fields as zero + // Not doing so gives in screwy results + for (int i=0; itransformComovingWithStar(d,vx,vy,vz,up,vp,wp,mx,my,mz, 1); + this->transformComovingWithStar(d_shell,vx_dep_tot,vy_dep_tot,vz_dep_tot,up,vp,wp,mx,my,mz, 1); + + /* + Use averaged quantities across multiple cells so that deposition is stable. + vmean is used to determine whether the supernova shell calculation should proceed: + M_shell > 0 iff v_shell > v_gas + */ + double Z_mean=0, d_mean=0, n_mean=0, v_mean=0, mu_mean=0; + double mu = enzo_config->ppm_mol_weight; + + for (int ix_ = ix-1; ix_ < ix+2; ix_++) { + for (int iy_ = iy-1; iy_ < iy+2; iy_++) { + for (int iz_ = iz-1; iz_ < iz+2; iz_++) { + int ind = INDEX(ix_,iy_,iz_, mx,my); + Z_mean += mf[ind] / d[ind]; + // TODO: make EnzoComputeMolecularWeight, and access mu_field here? + + #ifdef CONFIG_USE_GRACKLE + int primordial_chemistry = (enzo::config()->method_grackle_chemistry)->primordial_chemistry; + if (primordial_chemistry > 0) { + mu = d_el[ind] + dHI[ind] + dHII[ind] + 0.25*(dHeI[ind]+dHeII[ind]+dHeIII[ind]); + + if (primordial_chemistry > 1) { + mu += dHM[ind] + 0.5*(dH2I[ind]+dH2II[ind]); + } + if (primordial_chemistry > 2) { + mu += 0.5*(dDI[ind] + dDII[ind]) + dHDI[ind]/3.0; + } + } + mu /= d[ind]; + #endif + + mu_mean += mu; + d_mean += d[ind]; + } + } + } + + v_mean *= vunit / 27; // cm/s + Z_mean *= 1/(27 * enzo_constants::metallicity_solar); // Zsun + mu_mean /= 27; + d_mean *= rhounit/27; // g/cm^3 + n_mean = d_mean / (enzo_constants::mass_hydrogen/mu_mean); + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Zmean = %e Dmean = %e (%e) mu_mean = %e ", Z_mean, d_mean, d_mean / rhounit, mu_mean); + CkPrintf("Nmean = %f vmean = %f\n", n_mean, v_mean/1e5); + #endif + + double Z_Zsun = Z_mean; + double fz = std::min(2.0, pow(std::max(Z_Zsun,0.01),-0.14)); + + // Cooling radius as in Hopkins, but as an average over cells + double CoolingRadius = 28.4 * pow(std::max(0.1, n_mean), -3.0 / 7.0) * pow(ejectaEnergy / 1.0e51, 2.0 / 7.0) * fz; // pc + double coupledEnergy = ejectaEnergy; + double cellwidth_pc = hx*lunit / enzo_constants::pc_cm; // pc + double dxRatio = cellwidth_pc / CoolingRadius; + + /* + We want to couple one of four phases: free expansion, Sedov-taylor, shell formation, or terminal + The first three phases are take forms from Kim & Ostriker 2015, the last from Cioffi 1988 + */ + + + // if we resolve free expansion, all energy is thermally coupled + double p_free = sqrt(2*ejectaMass*enzo_constants::mass_solar*ejectaEnergy) / enzo_constants::mass_solar/1e5; + double r_free = 2.75*pow(ejectaMass / 3 / n_mean, 1.0/3); // free expansion radius eq. 2 + + double t3_sedov = pow(std::max(r_free, std::max(cellwidth_pc, r_free)) * enzo_constants::pc_cm / + (5.0 * enzo_constants::pc_cm * pow(ejectaEnergy / 1e51 / n_mean, 1.0 / 5.0)), 5./ 2.); + + double p_sedov = 2.21e4 * pow(ejectaEnergy / 1e51, 4. / 5.) * + pow(n_mean, 1. / 5.) * pow(t3_sedov, 3. / 5.); // eq 16 + + // shell formation radius eq. 8 + double r_shellform = 22.6*pow(ejectaEnergy/1e51,0.29) * pow(n_mean,-0.42); + double p_shellform = 3.1e5 * pow(ejectaEnergy / 1e51, 0.94) * pow(n_mean, -0.13); // p_sf = m_sf*v_sf eq 9,11 + + /* + terminal momentum + We picked the analytic form from Thornton as it matched high-resolution SN we ran the best. + */ + + double pTerminal; + if (Z_Zsun > 0.01) { + // Thornton, 1998, M_s * V_s, eqs 22, 23, 33, 34 + pTerminal = 1.6272e5 * pow(ejectaEnergy/1e51, 13./14.) * pow(n_mean, -0.25) * pow(Z_Zsun, -0.36); + } + + else { + pTerminal = 8.3619e5 * pow(ejectaEnergy/1e51, 13./14.) * pow(n_mean, -0.25); + } + + // compute the temperature (returns temperature in Kelvin) + EnzoComputeTemperature compute_temperature + (enzo_config->ppm_density_floor, + enzo_config->ppm_temperature_floor, + enzo_config->ppm_mol_weight, + enzo_config->physics_cosmology); + + compute_temperature.compute(enzo_block); + + double T = temperature[index]; + + double cSound = sqrt(enzo_constants::kboltz*T/(mu_mean*enzo_constants::mass_hydrogen)) / 1e5; // km/s + + // fading radius of a supernova, using gas energy of the host cell and ideal gas approximations + double r_fade = std::max(66.0*pow(ejectaEnergy/1e51, 0.32)*pow(n_mean, -0.37)*pow(cSound/10, -2.0/5.0), + CoolingRadius * 1.5); + double fadeRatio = cellwidth_pc/r_fade; + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Fading: T = %e; Cs = %e; R_f = %e; fadeR = %f\n", T, cSound, r_fade, fadeRatio); + #endif + + double coupledMomenta = 0.0; + double eKinetic = 0.0; + + double cw_eff = cellwidth_pc; + + double dx_eff = cw_eff / CoolingRadius; // our resolution with respect to the cooling radius + double fader = cw_eff / r_fade; + + double shellVelocity = 413.0 * pow(n_mean, 1.0 / 7.0) * pow(Z_Zsun, 3.0 / 14.0) + * pow(coupledEnergy / 1e51, 1.0 / 14.0); + + double ratio_pds = cw_eff/r_shellform; + + shellVelocity *= ratio_pds > 1 ? pow(dx_eff, -7.0 / 3.0) : 1; //km/s + + float beta = std::max( 1.0, shellVelocity / std::max(1.0, cSound)); + + // critical density to skip snowplough (remnant combines with ISM before radiative phase); eq 4.9 cioffi 1988 + float nCritical = 0.0038 *(pow(n_mean * T / 1e4, 7.0/9.0) + * pow(beta, 14.0/9.0))/(pow(ejectaEnergy/1e51, 1.0/9.0) + * pow(std::max(0.0001, Z_Zsun), 1.0/3.0)); // n/cc + + float r_merge = std::max(151.0 * pow((ejectaEnergy/1e51) / + beta / beta / n_mean / T * 1e4, 1./3.), 1.5 * r_fade); // pc + float merger = cw_eff / r_merge; + bool faded = fader > 1; + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf( + "STARSS_FB: RADII: cell = %e, free = %e, shellform = %e, cooling = %e, fade = %e t_3=%e, R_m=%e\n", + cellwidth_pc, r_free, r_shellform, CoolingRadius, r_fade, t3_sedov, r_merge); + #endif + + if (! winds) + { // this calculation for SNe only + + // free-expansion phase + if (cw_eff < r_free){ + coupledMomenta = 0.0;// Thermal coupling only at free expansion limit. + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: modifying free phase: p = %e\n", coupledMomenta); + #endif + } + + // Sedov phase + if (r_free < cw_eff && dx_eff <= 1){ + coupledMomenta = std::min(p_sedov, pTerminal*dx_eff); + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Coupling Sedov-Terminal phase: p = %e (ps = %e, pt = %e, dxe = %e)\n", + coupledMomenta, p_sedov, pTerminal, dx_eff); + #endif + } + + // terminal phase + if (dx_eff > 1){ + coupledMomenta = pTerminal/ sqrt(std::min(1.5, dx_eff)); + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Coupling Terminal phase: p = %e; dx_eff = %e\n", coupledMomenta, dx_eff); + #endif + } + + // fading phase + if (fader > 1 && enzo_config->method_feedback_fade_SNR){ + coupledMomenta = pTerminal * (1.0 - tanh(pow(fader * merger, 2.5))); + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Coupling Fading phase: p = %e\n", coupledMomenta); + #endif + } + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Checking critical density metric..." + "(n_mean = %e; N_Crit = %e; factors: %e %e %e; beta = %e/%e == %e; rmerge = %e)\n", + n_mean, nCritical, pow(n_mean * T / 1e4, 7.0/9.0), + pow(ejectaEnergy/1e51, 1.0/9.0), pow(fz, 1.0/3.0), + shellVelocity , cSound, beta, r_merge); + #endif + + if (n_mean <= 10.0 * nCritical){ // in high-pressure, low nb, p_t doesn't hold since there is essentially no radiative phase. + // thermal energy dominates the evolution (Tang, 2005, doi 10.1086/430875 ) + // We inject 100% thermal energy to simulate this recombining with the ISM + // and rely on the hydro and the thermal radiation to arrive at the right solution + coupledMomenta = coupledMomenta * (1.0-tanh(pow(1.45*nCritical/n_mean, 6.5))); + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Adjusting for high-pressure low-n phase (thermal coupling: Nc = %e): p = %e\n", + nCritical, coupledMomenta); + #endif + } + + if (T > 1e6 && coupledMomenta > 1e5) { + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Coupling high momenta to very hot gas!! (p= %e, T= %e, n_c = %e)\n", coupledMomenta, T, nCritical); + #endif + } + } // endif no winds + + else { // if winds + coupledMomenta = sqrt(ejectaMass*enzo_constants::mass_solar * 0.5 * ejectaEnergy)/enzo_constants::mass_solar/1e5; + } + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Calculated p = %e (sq_fact = %e; p_f = %e; p_t = %e; mcell = %e; mcpl = %e)\n", + coupledMomenta, (d_mean * cell_volume/enzo_constants::mass_solar) / ejectaMass * 63, p_free, + pTerminal, d_mean * cell_volume /enzo_constants::mass_solar, ejectaMass/27.0); + #endif + + /* + If resolution is in a range comparable to Rcool and + Analytic SNR shell mass is on, adjust the shell mass + upper range of applicability for shell mass is determined by + local average gas velocity (v_shell ~ c_s = no shell) + */ + + double centralMass = 0.0; + double centralMetals = 0.0; + double maxEvacFraction = 0.75; // TODO: make this a parameter + double shellMass = 0.0; + + if (coupledEnergy > 0 && AnalyticSNRShellMass && !winds) + { + if (dx_eff < 1) + shellMass = 4*cello::pi/3 * d_mean * pow(cw_eff*enzo_constants::pc_cm,3) / enzo_constants::mass_solar; + + else if (dx_eff > 1) { + if (Z_Zsun > 0.01) + shellMass = 1.41e4*pow(ejectaEnergy/1e51, 6./7.) * pow(d_mean, -0.24) * pow(Z_Zsun, -0.27); + else + shellMass = 4.89e4*pow(ejectaEnergy/1e51, 6./7.) * pow(d_mean, -0.24); + } + + if (shellMass > 0) + { + // get total mass in cells affected by SN + for (int ix_ = ix-1; ix_ <= ix+1; ix_++) { + for (int iy_ = iy-1; iy_ <= iy+1; iy_++) { + for (int iz_ = iz-1; iz_ <= iz+1; iz_++) { + int flat = INDEX(ix_,iy_,iz_,mx,my); + centralMass += d[flat]; + centralMetals += mf[flat]; + } + } + } // endfor ix_ + + centralMass *= rho_to_m; // Mass in Msun + + // Can't let host cells evacuate completely! + // Enforce maximum shellMass + if (shellMass > maxEvacFraction * centralMass) { + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS: Shell mass too high for host cells: Rescaling %e -> %e\n", + shellMass, maxEvacFraction * centralMass); + #endif + shellMass = maxEvacFraction * centralMass; + } + } // endif shellMass > 0 + + } // endif coupledEnergy >0 && AnalyticSNRShellMass && !winds + + double shellMetals = std::min(maxEvacFraction*centralMetals * rho_to_m, + Z_Zsun * enzo_constants::metallicity_solar * shellMass); //Msun + +#ifdef DEBUG_FEEDBACK_STARSS + if (AnalyticSNRShellMass) { + CkPrintf("STARSS_FB: Shell_m = %e Shell_z = %e shell_V= %e P = %e M_C = %e\n", + shellMass, shellMetals, shellVelocity, coupledMomenta, centralMass); + } +#endif + + double coupledMass = shellMass + ejectaMass; + + eKinetic = coupledMomenta * coupledMomenta / (2.0 *(coupledMass)) * enzo_constants::mass_solar * 1e10; + if (eKinetic > coupledEnergy && !winds){ + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Rescaling high kinetic energy %e -> ", eKinetic); + #endif + + coupledMomenta = sqrt(2.0 * (coupledMass*enzo_constants::mass_solar) * ejectaEnergy)/enzo_constants::mass_solar/1e5; + eKinetic = coupledMomenta * coupledMomenta / (2.0 *(coupledMass) * enzo_constants::mass_solar) * enzo_constants::mass_solar * 1e10; + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: %e; new p = %e\n", eKinetic, coupledMomenta); + #endif + } + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Ekinetic = %e Mass = %e\n", + eKinetic, d_mean * pow(lunit * hx, 3) / enzo_constants::mass_solar); + #endif + + if (eKinetic > 1e60) + { + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: winds = %d, Ekinetic = %e Mass = %e\n", + winds, eKinetic, d_mean * pow(lunit * hx, 3) / enzo_constants::mass_solar); + #endif + ERROR("EnzoMethodFeedbackSTARSS::deposit_feedback()","Ekinetic > 1e60 erg!\n"); + } + + double coupledGasEnergy = std::max(ejectaEnergy - eKinetic, 0.0); + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Coupled Gas Energy = %e\n", coupledGasEnergy); + #endif + + if (dxRatio > 1.0 && !winds){ + // Fade internal energy deposited if our resolution is worse + // than a cooling radius. + // + // Don't need to do this for stellar winds because winds are + // a lot weaker than SNe, so it makes essentially no difference. + coupledGasEnergy = (coupledGasEnergy * pow(dxRatio, -6.5)); + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Reducing gas energy... GE = %e\n", coupledGasEnergy); + #endif + } + + double coupledMetals = 0.0;//, SNIAmetals = 0.0, SNIImetals = 0.0, P3metals = 0.0; + if (winds) coupledMetals = ejectaMetals; // winds only couple to metals + if (AnalyticSNRShellMass && !winds) coupledMetals += shellMetals; + + ejectaMetals += nSNIA * 1.4; // add ejecta from Type Ia SNe + ejectaMetals += nSNII * (1.91 + 0.0479 * std::max(starZ, 1.65)); // add ejecta from Type II + + coupledMetals += ejectaMetals; + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Coupled Metals: %e %e %e\n", ejectaMetals, shellMetals, coupledMetals); + #endif + + if (!winds) coupledEnergy = std::min(coupledEnergy, eKinetic); + + // Subtract shell mass from the central cells + double minusRho=0, minusZ=0, msubtracted=0; + double remainMass=shellMass/rho_to_m, remainZ = shellMetals/rho_to_m; + if (shellMass > 0 && AnalyticSNRShellMass) + { + double zsubtracted=0; + msubtracted=0; + + for (int ix_ = ix-1; ix_ <= ix+1; ix_++) { + for (int iy_ = iy-1; iy_ <= iy+1; iy_++) { + for (int iz_ = iz-1; iz_ <= iz+1; iz_++) { + int flat = INDEX(ix_,iy_,iz_,mx,my); + + // cell left edges for CiC (starting from ghost zones) + double xcell = xm + (ix_+0.5 - gx)*hx; + double ycell = ym + (iy_+0.5 - gy)*hy; + double zcell = zm + (iz_+0.5 - gz)*hz; + + double window = Window(xp-xcell, yp-ycell, zp-zcell, hx); // CiC fraction + if (window <= 0) continue; + + double dpre = d[flat]; + double zpre = mf[flat]; + double pre_z_frac = zpre / dpre; + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS: Baryon Prior: d_Msun = %e, mc = %e, ms = %e, m_z = %e, z = %e\n", + d[flat] * rho_to_m, centralMass, shellMass, shellMetals, pre_z_frac); + #endif + + // subtract values from the "deposit" fields + d_dep[flat] -= std::min(window * remainMass, maxEvacFraction*dpre); + + minusRho += -1*d_dep[flat]; + msubtracted += -1*d_dep[flat]; + + mf_dep[flat] -= std::min(window * remainZ, maxEvacFraction*zpre); + + minusZ += -1*mf_dep[flat]; + zsubtracted += -1*mf_dep[flat]; + } // endfor iz_ + } // endfor iy_ + } // endfor ix_ + remainMass -= msubtracted; + remainZ -= zsubtracted; + } // endif shellMass > 0 && AnalyticSNRShellMass + + minusRho *= rho_to_m; //Msun + minusZ *= rho_to_m; //Msun + + double minusM = minusRho; + if (minusM != coupledMass - ejectaMass && shellMass > 0 && AnalyticSNRShellMass) + { + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: Of %e, only subtracted %e; rescaling the coupling mass\n", shellMass, + minusRho); + #endif + + coupledMass = minusM + ejectaMass; + coupledMetals = minusZ + ejectaMetals; // metal_density in Msun/code_volume + + // if we're in here, we found an inconsistent place where the mass within + // cannot support the expected shell. + // the only real choice, to keep things from exploding in a bad way, + // is to rescale the momentum accordingly. + // If this is triggered, dont expect the terminal momenta relations to hold up. + + if (coupledMomenta*coupledMomenta / + (2.0*coupledMass*enzo_constants::mass_solar)*enzo_constants::mass_solar * 1e10 > 1e51) { + coupledMomenta = std::min(coupledMomenta, sqrt(2.0 * (coupledMass*enzo_constants::mass_solar) + * ejectaEnergy)/enzo_constants::mass_solar/1e5); + + #ifdef DEBUG_FEEDBACK_STARSS + CkPrintf("STARSS_FB: rescaled momentum to %e (est KE = %e)\n", + coupledMomenta, coupledMomenta*coupledMomenta / (2*coupledMass) + * enzo_constants::mass_solar * 1e10); + #endif + } + + } // endif minusM != coupledMass - ejectaMass + + coupledEnergy += coupledGasEnergy; + + // put everything back into code units before CIC + // these values correspond to TOTAL (energy, mass, momentum) + // over all the coupling particles. + + coupledEnergy /= (eunit * cell_volume_code); // energy -> energy density + coupledGasEnergy /= (eunit * cell_volume_code); + coupledMass /= rho_to_m; // mass -> density + coupledMetals /= rho_to_m; + coupledMomenta /= (rho_to_m * vunit / 1e5); // momentum -> momentum_density + + + // Create arrays to pass into CiC routine + enzo_float coupledMass_list[nCouple], coupledMetals_list[nCouple]; + enzo_float coupledMomenta_x[nCouple], coupledMomenta_y[nCouple], coupledMomenta_z[nCouple]; + enzo_float coupledEnergy_list[nCouple], coupledGasEnergy_list[nCouple]; + + for (int n=0; ntransformComovingWithStar(d,vx,vy,vz,up,vp,wp,mx,my,mz, -1); + this->transformComovingWithStar(d_shell,vx_dep_tot,vy_dep_tot,vz_dep_tot,up,vp,wp,mx,my,mz, -1); + + // clear temporary fields + delete [] d_dep; + delete [] mf_dep; + delete [] te_dep; + delete [] ge_dep; + delete [] vx_dep; + delete [] vy_dep; + delete [] vz_dep; +} + + + +double EnzoMethodFeedbackSTARSS::timestep (Block * block) throw() +{ + // In general this is not needed, but could imagine putting timestep + // limiters in situations where, for example, one would want + // dt < star_lifetime (or something like that), especially if + // important things happen throughout the star's lifetime. + EnzoUnits * enzo_units = enzo::units(); + const EnzoConfig * enzo_config = enzo::config(); + + double dtStar = std::numeric_limits::max(); + if (block->level() >= sf_minimum_level_){ + const double pSNmax = 0.0005408 * enzo_config->method_star_maker_minimum_star_mass * + block->dt() * enzo_units->time() / enzo_constants::Myr_s * 1.25; + if (pSNmax > 1.0) dtStar = block->dt() * 1.0 / pSNmax; + } + + return dtStar; +} + + +// #endif diff --git a/src/Enzo/enzo_EnzoMethodFeedbackSTARSS.hpp b/src/Enzo/enzo_EnzoMethodFeedbackSTARSS.hpp new file mode 100644 index 0000000000..25b32bea73 --- /dev/null +++ b/src/Enzo/enzo_EnzoMethodFeedbackSTARSS.hpp @@ -0,0 +1,144 @@ +/// See LICENSE_CELLO file for license and copyright information + +/// @file enzo_EnzoMethodFeedbackSTARSS.hpp +/// @author Will Hicks (whicks@ucsd.edu) +/// @date +/// @brief Implements the STARSS model for stellar Feedback +/// as implemented into Enzo from Azton Wells +/// + + +//#ifdef NOTDEFINED //uncomment when ready to compile + +#ifndef ENZO_ENZO_METHOD_FEEDBACK_STARSS +#define ENZO_ENZO_METHOD_FEEDBACK_STARSS + +class EnzoMethodFeedbackSTARSS : public Method { + + /// @class EnzoMethodFeedbackSTARSS + /// @ingroup Enzo + /// @btief [\ref Enzo] Encapsulate Feedback Routines + +public: + + EnzoMethodFeedbackSTARSS(); + + /// Destructor + virtual ~EnzoMethodFeedbackSTARSS() throw() {}; + + /// Charm++ Pup::able declarations + PUPable_decl(EnzoMethodFeedbackSTARSS); + + /// Charm++ Pup::able migration Constructor + EnzoMethodFeedbackSTARSS (CkMigrateMessage *m) + : Method (m) + , ir_feedback_(-1) + { } + + /// Charm++ Pack / Unpack function + void pup(PUP::er &p); + + /// Apply the method + virtual void compute (Block * block) throw(); + + void compute_ (Block * block); + + /// name + virtual std::string name() throw() + { return "feedback"; } + + // Compute the maximum timestep for this method + virtual double timestep (Block * block) throw(); + + int determineSN (double age_Myr, int * nSNII, int * nSNIA, + double mass_Msun, double tunit, float dt); + + int determineWinds(double age_Myr, double * eWinds, double * mWinds, double * zWinds, + double mass_Msun, double metallicity_Zsun, double tunit, double dt); + + // this can raise errors -- remove const throw() ??? + void deposit_feedback (Block * block, + double ejectaEnergy, double ejectaMass, double ejectaMetals, + const enzo_float up, const enzo_float vp, const enzo_float wp, + const enzo_float xp, const enzo_float yp, const enzo_float zp, + const int ix, const int iy, const int iz, + const int winds, const int nSNII, const int nSNIa, + const double starZ) const throw(); + + void transformComovingWithStar(enzo_float * density, + enzo_float * velocity_x, enzo_float * velocity_y, enzo_float * velocity_z, + const enzo_float up, const enzo_float vp, const enzo_float wp, + const int mx, const int my, const int mz, int direction) const throw(); + + void add_accumulate_fields(EnzoBlock * enzo_block) throw(); + + // window function for calculating CiC fractions. May make more sense to put this in Cello somewhere + double Window(double xd, double yd, double zd, double width) const throw(); + +protected: // methods + void allocate_temporary_(EnzoBlock * enzo_block) + { + Field field = enzo_block->data()->field(); + field.allocate_temporary(i_d_dep); + field.allocate_temporary(i_te_dep); + field.allocate_temporary(i_ge_dep); + field.allocate_temporary(i_mf_dep); + field.allocate_temporary(i_vx_dep); + field.allocate_temporary(i_vy_dep); + field.allocate_temporary(i_vz_dep); + field.allocate_temporary(i_d_shell); + + field.allocate_temporary(i_d_dep_a); + field.allocate_temporary(i_te_dep_a); + field.allocate_temporary(i_ge_dep_a); + field.allocate_temporary(i_mf_dep_a); + field.allocate_temporary(i_vx_dep_a); + field.allocate_temporary(i_vy_dep_a); + field.allocate_temporary(i_vz_dep_a); + field.allocate_temporary(i_d_shell_a); + } + + void deallocate_temporary_(EnzoBlock * enzo_block) + { + Field field = enzo_block->data()->field(); + field.deallocate_temporary(i_d_dep); + field.deallocate_temporary(i_te_dep); + field.deallocate_temporary(i_ge_dep); + field.deallocate_temporary(i_mf_dep); + field.deallocate_temporary(i_vx_dep); + field.deallocate_temporary(i_vy_dep); + field.deallocate_temporary(i_vz_dep); + field.deallocate_temporary(i_d_shell); + + field.deallocate_temporary(i_d_dep_a); + field.deallocate_temporary(i_te_dep_a); + field.deallocate_temporary(i_ge_dep_a); + field.deallocate_temporary(i_mf_dep_a); + field.deallocate_temporary(i_vx_dep_a); + field.deallocate_temporary(i_vy_dep_a); + field.deallocate_temporary(i_vz_dep_a); + field.deallocate_temporary(i_d_shell_a); + } + +protected: + + int sf_minimum_level_; + int single_sn_; + int NEvents; + + // Refresh ID + int ir_feedback_; + + // deposit field id's + int i_d_dep , i_d_dep_a; + int i_te_dep, i_te_dep_a; + int i_ge_dep, i_ge_dep_a; + int i_mf_dep, i_mf_dep_a; + int i_vx_dep, i_vx_dep_a; + int i_vy_dep, i_vy_dep_a; + int i_vz_dep, i_vz_dep_a; + int i_d_shell, i_d_shell_a; + +}; + +#endif diff --git a/src/Enzo/enzo_EnzoMethodFire2Feedback.cpp b/src/Enzo/enzo_EnzoMethodFire2Feedback.cpp deleted file mode 100644 index 94349fab31..0000000000 --- a/src/Enzo/enzo_EnzoMethodFire2Feedback.cpp +++ /dev/null @@ -1,652 +0,0 @@ - -/// See LICENSE_CELLO file for license and copyright information - -/// @file enzo_EnzoMethodFire2Feedback.cpp -/// @author Andrew Emerick (aemerick11@gmail.com) -/// @date -/// @brief Implements the FIRE2 model for stellar feedback -/// as implemented into Enzo as the MechStars model -/// by Azton Wells. This is intended to be a direct copy-over of the -/// MechStars methods. -/// TODO: List differences / changes here -/// - - -#include "cello.hpp" -#include "enzo.hpp" - -#include - -#ifdef NOTDEFINED // for now... since not done coding - -// ============================================================================= -// splice these off to a different file (later) - -int determineSN(float age, int* nSNII, int* nSNIA, - float massmass_solar, float TimeUnits, float dt){ - - const EnzoConfig * enzo_config = enzo::config(); - - if (NEvents > 0){ - *nSNII = 1; - NEvents -= 1; - return SUCCESS; - } - /* else, calculate SN rate, probability and determine number of events */ - int seed = time(NULL): - *nSNII = 0; - *nSNIA = 0; - float RII=0, RIA=0, PII=0, PIA=0, random = 0; - if (enzo_config->method_feedback_fire2_single_sn && NEvents < 0) - { - // printf("Calculating rates\n"); - /* age-dependent rates */ - if (age < 3.401) - { - RII = 0.0; - RIA = 0.0; - } - if (3.401 <= age && age< 10.37) - { - RII = 5.408e-4; - RIA = 0.0; - } - if (10.37 <= age && age < 37.53) - { - RII = 2.516e-4; - RIA = 0.0; - } - if (37.53 <= age) - { - RII = 0.0; - RIA = 5.3e-8+1.6e-5*exp(-0.5*pow((age-50.0)/10.0, 2)); - } - // fprintf(stdout, "Rates: %f %f %f\n", age, RII, RIA); - /* rates -> probabilities */ - if (RII > 0){ - srand(seed); - // printf("Zcpl = %e", zCouple); - PII = RII * massmass_solar / enzo_constants::Myr_s *TimeUnits*dt; - // printf("PII =%f\n %f %e %f\n", PII, RII, massmass_solar, age); - random = float(rand())/float(RAND_MAX); - if (PII > 1.0 && enzo_config->method_feedback_fire2_unrestricted_SN){ - int round = (int)PII; - *nSNII = round; - PII -= round; - } - if (PII > 1.0 && !enzo_config->method_feedback_fire2_unrestricted_SN){ - ERROR("determineSN: PII too large!"); - } - int psn = *nSNII; - if (random < PII){ - *nSNII = psn+1; - } - } - // printf("RANDOM = %f\n", random); - // printf("N SNII=%d\n",*nSNII); - - if (RIA > 0){ - srand(seed); - PIA = RIA*massmass_solar / cello:Myr_s *TimeUnits*dt; - float random = float(rand())/float(RAND_MAX); - - if (PIA > 1.0 && enzo_config->method_feedback_fire2_unrestricted_SN) - { - int round = int(PIA); - *nSNIA = round; - PIA -= round; - } - int psn = *nSNIA; - - if (random < PIA) - *nSNIA = psn+1; - // if (*nSNIA > 0) - // fprintf(stdout, "PIA = %f\n", PIA); - } - } - return 1; -} - -int determineWinds(float age, float* eWinds, float* mWinds, float* zWinds, - float massMsun, float zZsun, float TimeUnits, float dt){ - - // age in Myr - - const EnzoConfig * enzo_config = enzo::config(); - - float windE = 0, windM = 0, windZ = 0.0; - float wind_factor = 0.0; - float e_factor = 0.0; - - - if (enzo_config->method_feedback_fire2_stellarwinds){ - - if(zZsun > 3.0) zZsun = 3.0; - if(zZsun < 0.01) zZsun = 0.01; - - if (age <= 1.0){ - wind_factor = 11.6846; - } else { - if (age <= 3.5){ - wind_factor = 11.6846 * zZsun * - pow(10.0, 1.838*(0.79+log10(zZsun))*(log10(age))); - } else { - if (age <= 100.0){ - wind_factor = 72.1215*pow(0.001*age/0.0035,-3.25)+0.0103; - } else { - wind_factor = 1.03*pow(0.001*age,-1.1)/(12.9-log(age*0.001)); - } // age < 100 - } // age < 3.5 - - } // age < 1.0 - - if (age <= 100.0){ - e_factor = 0.0013+16.0/(1.0+pow(age*0.001/0.0025,1.4)+pow(age*0.001/0.01,5.0)); - } else{ - e_factor = 0.0013; - } - - wind_factor *= (enzo_config->method_feedback_fire2_gasreturnfraction *\ - 0.001 * (dt * (TimeUnits / enzo_constants::Myr_s)); // 0.001 converts from Msun/Gyr to Msun/Myr - wind_factor = 1.0 - exp(-wind_factor); - wind_factor = std::min(1.0, wind_factor); - wind_factor *= 1.4 * 0.291175; - - double n_wind = (double)floor(wind_factor/rfrac); wind_factor -= n_wind *rfrac; - windM = n_wind * rfrac; - if( float(rand())/float(RAND_MAX) < wind_factor / rfrac){windM += rfrac}; - - windM *= massMsun; // convert to total mas from fraction - windZ = (0.016+0.0041+0.0118) * windM; // C+N+O yields with no Z scaling (total masss of metals) - windE = e_factor * 1e12 * windM - } - - - *mWinds = windM; - *zWinds = windZ; - *eWinds = windE; - - return 1; -/* - - if (enzo_config->method_feedback_fire2_stellarwinds && - oldEnough && massMsun > 11){ - - if (0.001 < age && age < 1.0){ - wind_factor =4.763 * min((0.01 + zZsun), 1.0) ; - } - if (1 <= age && age < 3.5){ - wind_factor = 4.763*min(0.01+zZsun, 1.0)* - pow(age, 1.45+0.08*min(log(zZsun), 1.0)); - } - if (3.5 <= age && age < 100){ - wind_factor = 29.4*pow(age/3.5, -3.25)+0.0042; - - } - if (age < 100){ - float d = powl(1+age/2.5, 1.4); - float a50 = powl(double(age)/10.0, 5.0); - e_factor = 5.94e4 / d + a50 +4.83; - - } - if (100 <= age){ - e_factor = 4.83; - wind_factor = 0.42*pow(age/1000, -1.1)/(19.81/log(age)); - } - windM = massMsun * wind_factor; //Msun/Gyr - windM = windM*dtFixed*TimeUnits/3.1557e16; //Msun - // if (debug) printf("First winds mass = %e\nFrom wf = %f, dt=%f Z = %e\n", windM, wind_factor, dtFixed, zZsun); - //printf("eFactor = %f age = %f\n", e_factor, age); - if (windM > massMsun){ - printf("Winds too large Mw = %e, Mp = %e age=%f, Z = %e\n", - windM, massMsun, age, zZsun); - windM = 0.125*massMsun; // limit loss to huge if necessary. - } - windZ = max(0.02, 0.016+0.0041*max(zZsun, 1.65)+0.0118)*windM; - windE = e_factor * 1e12 * windM; - //fprintf(stdout, "Age = %e Ewinds = %e Mwinds = %e Zwinds = %e Zsun = %e\n", - // age, windE, windM, windZ, zZsun); - *mWinds = windM; - *zWinds = windZ; - *eWinds = windE; - } - - return SUCCESS; - } -*/ -} - -// ============================================================================= -// ============================================================================= -// ============================================================================= - -EnzoMethodFire2Feedback::EnzoMethodFire2Feedback -() - : Method() -{ - cello::particle_descr()->check_particle_attribute("star","mass"); - FieldDescr * field_descr = cello::field_descr(); - const EnzoConfig * enzo_config = enzo::config(); - EnzoUnits * enzo_units = enzo::units(); - - // required fields - cello::define_field ("density"); - cello::define_field ("pressure"); - cello::define_field ("internal_energy"); - cello::define_field ("total_energy"); - - cello::simulation()->refresh_set_name(ir_post_,name()); - Refresh * refresh = cello::refresh(ir_post_); - refresh->add_all_fields(); - - sf_minimum_level_ = enzo_config->method_feedback_fire2_minimum_level; - single_sn_ = enzo_config->method_feedback_fire2_single_sn; - - return; -} - -void EnzoMethodFire2Feedback::pup (PUP::er &p) -{ - /// NOTE: Change this function whenever attributes change - - TRACEPUP; - - Method::pup(p); - - p | sf_minimum_level_; - p | single_sn_; - - return; -} - -void EnzoMethodFire2Feedback::compute (Block * block) throw() -{ - - if (block->is_leaf()){ - this->compute_(block); - } - - block->compute_done(); - - return; -} - -void EnzoMethodFire2Feedback::compute_ (Block * block) -{ - - //---------------------------------------------------- - // some constants here that might get moved to parameters or something else - const float SNII_ejecta_mass = 10.5; - const float SNIa_ejecta_mass = 1.4; - const float z_solar = 0.02; // Solar metal fraction (fire2) - //----------------------------------------------------- - - EnzoBlock * enzo_block = enzo::block(block); - const EnzoConfig * enzo_config = enzo::config(); - Particle particle = enzo_block->data()->particle(); - EnzoUnits * enzo_units = enzo::units(); - - double current_time = block->time(); - - Field field = enzo_block->data()->field(); - - enzo_float * d = (enzo_float *) field.values("density"); - enzo_float * te = (enzo_float *) field.values("total_energy"); - enzo_float * ge = (enzo_float *) field.values("internal_energy"); - enzo_float * mf = (enzo_float *) field.values("metal_density"); - - // Obtain grid sizes and ghost sizes - int mx, my, mz, gx, gy, gz, nx, ny, nz; - double xm, ym, zm, xp, yp, zp, hx, hy, hz; - field.size(&nx,&ny,&nz); - field.ghost_depth(0,&gx,&gy,&gz); - block->data()->lower(&xm,&ym,&zm); - block->data()->upper(&xp,&yp,&zp); - field.cell_width(xm,xp,&hx,ym,yp,&hy,zm,zp,&hz); - - mx = nx + 2*gx; - my = ny + 2*gy; - mz = nz + 2*gz; - - EnzoPhysicsCosmology * cosmology = enzo::cosmology(); - enzo_float cosmo_a = 1.0; - - const int rank = cello::rank(); - - double current_time = block->time(); - if (cosmology) { - enzo_float cosmo_dadt = 0.0; - double dt = block->dt(); - cosmology->compute_expansion_factor(&cosmo_a,&cosmo_dadt,current_time+0.5*dt); - if (rank >= 1) hx *= cosmo_a; - if (rank >= 2) hy *= cosmo_a; - if (rank >= 3) hz *= cosmo_a; - } - - // apply feedback depending on particle type - // for now, just do this for all star particles - - int it = particle.type_index("star"); - int count = 0; - - // AE: Copying this from MechStars_FeedbackRoutine - - int numSN = 0; // counter of SN events - int count = 0; // counter of particles - - if (particle.num_particles(it) <= 0) return; // nothing to do here - - const int ia_m = particle.attribute_index (it, "mass"); - const int ia_x = particle.attribute_index (it, "x"); - const int ia_y = particle.attribute_index (it, "y"); - const int ia_z = particle.attribute_index (it, "z"); - const int ia_vx = particle.attribute_index (it, "vx"); - const int ia_vy = particle.attribute_index (it, "vy"); - const int ia_vz = particle.attribute_index (it, "vz"); - const int ia_l = particle.attribute_index (it, "lifetime"); - const int ia_c = particle.attribute_index (it, "creation_time"); - const int ia_mf = particle.attribute_index (it, "metal_fraction"); - const int ia_sn = particle.attribute_index (it, "number_of_sn"); // name change? - - const int dm = particle.stride(it, ia_m); - const int dp = particle.stride(it, ia_x); - const int dv = particle.stride(it, ia_vx); - const int dl = particle.stride(it, ia_l); - const int dc = particle.stride(it, ia_c); - const int dmf = particle.stride(it, ia_mf); - const int dsn = particle.stride(it, ia_sn); - - const int nb = particle.num_batches(it); - - for (int ib=0; ib 0.0 && plifetime[ipdl] > 0.0){ - const double age = (current_time - pcreation[ipdc]) * enzo_units->time() / enzo_constants::Myr_s; - count++; // increment particles examined here - - // compute coordinates of central feedback cell - // this must account for ghost zones - double xcell = (xpos - xm) / hx + gx - 0.5; - double ycell = (ypos - ym) / hy + gy - 0.5; - double zcell = (zpos - zm) / hz + gz - 0.5; - - int ix = ((int) floor(xcell + 0.5)); - int iy = ((int) floor(ycell + 0.5)); - int iz = ((int) floor(zcell + 0.5)); - - // do cell shifting here if we need to... assume not for now ... - - int index = INDEX(ix,iy,iz,mx,my); // 1D cell index of star position - - //if(pmass[ipdm] * enzo_units->mass_units() / enzo_constants::mass_solar - // < enzo_config->method_star_maker_minimum_star_mass){ - // - //} - // skipping continual formation routines here (May need later) - - int nSNII = 0, nSNIa = 0; - double SNMassEjected = 0.0, SNMetalEjected = 0.0; - - /* determine how many supernova events */ - if (single_sn_){ - - /* Determine SN events from rates (currently taken from Hopkins 2018) */ - - determineSN(age, &nSNII, &nSNIa, pmass[ipdm] * enzo_units->mass_units() / enzo_constants::mass_solar, - enzo_units->time_units(), block->dt()); - - numSN += nSNII + nSNIa; - -#ifdef FEEDBACK_DEBUG - if (numSN){ - CkPrintf("Fire2Feedback SNe %d %d level = %d age = %f\n", nSNII, nSNIa, level, age); - } -#endif - - /* AJE: Can I just change this to a single call to deposit feedback for - the total ejecta of SNe and winds ? */ - - if (numSN > 0){ - /* set feedback properties based on number and types of SN */ - double energySN = (numSN) * 1.0e51; - - SNMassEjected = SNII_ejecta_mass * nSNII + - SNIa_ejecta_mass * nSNIa; // AE: split this in two channels for yields - - const double starZ = pmetal[ipdmf] / z_solar; - - // compute metal mass ejected here - - SNMetalEjected = ; - - /* Fixed mass ejecta */ - - this->deposit_feedback( energySN, SNMassEjected, SNMetalEjected, - pvx[ipdv],pvy[ipdv],pvz[ipdv], - px[ipdp],py[ipdp],pz[ipdp], - ix, iy, iz, mx, my, mz, 0, nSNII, nSNIa, starZ); // removed P3 - - // add counter for number of SNe for the particle - psncounter[ipsn] += numSN; - - - - } // if nSNII or nSNIa > 0 - } // if single SN - - // Now do stellar winds - if (enzo_config->method_feedback_fire2_stellarwinds){ - - const double starZ = pmetal[ipdmf] / z_solar; - - determineWinds(age, &windEnergy, &windMass, &windMetals, - pmass[ipdm] * pmass[ipdm] * enzo_units->mass_units() / enzo_constants::mass_solar, - starZ, enzo_units->time_units(), block->dt()); - - /* Error messages go here */ - - if (windMass > 0){ - this->deposit_feedback( windEnergy, windMass, windMetals, - pvx[ipdv],pvy[ipdv],pvz[ipdv], - px[ipdp],py[ipdp],pz[ipdp], - ix, iy, iz, mx, my, mz, 1, 0, 0, 0.0); // removed P3 - - } // if wind mass > 0 - } // if winds - - pmass[ipdm] -= std::max(0.0, - (windMass + SNMassEjected) *enzo_units->mass_units()/enzo_constants::mass_solar); - - // - - - } // if mass and lifetime > 0 - - } // end particle loop - } // end batch loop - - - if (count > 0){ - CkPrintf("Fire2Feedback: Num FB particles = %d Events = %d FeedbackTime %e\n", - count, numSN, 0.00); - } - - return; -} - - - - -// ---------------------------------------------------------------------------- - - -void EnzoMethodFire2Feedback::deposit_feedbak (EnzoBlock * enzo_block, - const float ejectaEnergy, - const float ejectaMass, - const float ejectaMetals, - const enzo_float up, const enzo_float vp, const enzo_float wp, - const enzo_float xp, const enzo_float yp, const enzo_float zp, - const int ix, const int iy, const int iz, - const int wind_mode, const int nSNII, - const int nSNIa, - enzo_float starZ); - -) const throw(){ - /* - This routine will create an isocahedron of coupling particles, where we determine - the feedback quantities. The vertices of the isocahedron are ~coupled particles - and all have radius dx from the source particle. - Each vertex particle will then be CIC deposited to the grid! - */ - - Field field = enzo_block->data()->field(); - // Obtain grid sizes and ghost sizes - int mx, my, mz, gx, gy, gz, nx, ny, nz; - double xm, ym, zm, xp, yp, zp, hx, hy, hz; - field.size(&nx,&ny,&nz); - field.ghost_depth(0,&gx,&gy,&gz); - enzo_block->data()->lower(&xm,&ym,&zm); - enzo_block->data()->upper(&xp,&yp,&zp); - field.cell_width(xm,xp,&hx,ym,yp,&hy,zm,zp,&hz); - - mx = nx + 2*gx; - my = ny + 2*gy; - mz = nz + 2*gz; - - enzo_float * d = (enzo_float *) field.values("density"); - enzo_float * te = (enzo_float *) field.values("total_energy"); - enzo_float * ge = (enzo_float *) field.values("internal_energy"); - - enzo_float * vx = (enzo_float *) field.values("velocity_x"); - enzo_float * vy = (enzo_float *) field.values("velocity_y"); - enzo_float * vz = (enzo_float *) field.values("velocity_z"); - - enzo_float * mf = (enzo_float *) field.values("metal_density"); - - const int index = INDEX(ix,iy,iz,mx,my); - - int stretch_factor = 1.0; //1.5/sin(M_PI/10.0); // How far should cloud particles be from their host - // in units of dx. Since the cloud forms a sphere shell, stretchFactor > 1 is not recommended - - const float phi = (1.0 + sqrt(5.0)) * 0.5 ; // golden ratio - const float iphi = 1.0 / phi; - - /* A DODECAHEDRON+ISOCAHEDRON */ - const int nCouple = 26; - const float A = stretchFactor * dx; - - /* - Make a cloud of coupling particles; - each is on a grid with dx spacing - from the host particles. - */ - - enzo_float CloudParticleVectorX[] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1}; - // {1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0, - // iphi, iphi, -iphi, -iphi, phi, phi, -phi, -phi, - // 0, 0, 0, 0, 1, 1, -1, -1, phi,-phi, phi,-phi}; - enzo_float CloudParticleVectorY[] = {1, 1, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, 0, -1, -1, -1, 1, 1, 1, 0, 0, -1, -1, -1, 0}; - // {1,1,-1,-1, 1, 1, -1, -1, iphi, iphi, -iphi, -iphi, - // phi, -phi, phi,-phi, 0, 0, 0, 0,1, 1, -1, -1, - // phi, -phi, -phi, phi, 0, 0, 0, 0}; - enzo_float CloudParticleVectorZ[] = {1, 0, -1, 1, -1, 1, 0, -1, 0, 1, 0, -1, 1, -1, 1, 0, -1, 1, 0, -1, 1, -1, 1, 0, -1, 0}; - // {1,-1, 1,-1, 1,-1, 1,-1, phi,-phi, phi,-phi, - // 0, 0, 0, 0, iphi, -iphi, iphi, -iphi, - // phi, -phi, -phi, phi, 0, 0, 0, 0, 1, 1, -1, -1}; - enzo_float weightsVector[nCouple]; - /* Set position of feedback cloud particles */ - - enzo_float CloudParticlePositionX[nCouple]; - enzo_float CloudParticlePositionY[nCouple]; - enzo_float CloudParticlePositionZ[nCouple]; - - /*all possible values of x,y,z with origin particle at x=y=z=0.0 */ - for (int cpInd = 0; cpInd < nCouple; cpInd++) - { - float norm = sqrt(CloudParticleVectorX[cpInd] * CloudParticleVectorX[cpInd] + - CloudParticleVectorY[cpInd] * CloudParticleVectorY[cpInd] + - CloudParticleVectorZ[cpInd] * CloudParticleVectorZ[cpInd]); - float inv_norm = 1.0 / norm; - float xbaMag = A * A * norm * norm; - /* in this cloud, take the coupling particle position as 0.5, 0.5, 0.5 */ - CloudParticlePositionX[cpInd] = xp - CloudParticleVectorX[cpInd]*inv_norm * A; - CloudParticlePositionY[cpInd] = yp - CloudParticleVectorY[cpInd]*inv_norm * A; - CloudParticlePositionZ[cpInd] = zp - CloudParticleVectorZ[cpInd]*inv_norm * A; - weightsVector[cpInd] = 0.5 * (1. - 1. / (1. + 1. / 4. / cello::pi / 26. / xbaMag / cello::pi)); - /* turn the vectors back into unit-vectors */ - CloudParticleVectorZ[cpInd] *= inv_norm; - CloudParticleVectorY[cpInd] *= inv_norm; - CloudParticleVectorX[cpInd] *= inv_norm; - } - float weightsSum = 0.0; - for (int wind = 0; wind < nCouple; wind++) - { - weightsSum += weightsVector[wind]; - } - for (int wind = 0; wind < nCouple; wind++) - { - weightsVector[wind] /= weightsSum; - if (weightsVector[wind] == 0 || isnan(weightsVector[wind])) - { - ERROR("EnzoMethodFire2Feedback: deposit_feedback: NaN weight Vector!") - } - } - -/* AJE LEFT OFF HERE */ - - - return; -} - - - -double EnzoMethodFire2Feedback::timestep (Block * block) throw() -{ - // In general this is not needed, but could imagine putting timestep - // limiters in situations where, for example, one would want - // dt < star_lifetime (or something like that), especially if - // important things happen throughout the star's lifetime. - EnzoUnits * enzo_units = enzo::units(); - const EnzoConfig * enzo_config = enzo::config(); - - double dtStar = std::numeric_limits::max(); - if (block->level() >= sf_minimum_level_){ - const double pSNmax = 0.0005408 * enzo_config->method_star_maker_minimum_star_mass * - block->dt() * enzo_units->time() / enzo_constants::Myr_s * 1.25; - if (pSNmax > 1.0) dtStar = dt * 1.0 / pSNmax; - } - - return dtStar; -} - - -#endif diff --git a/src/Enzo/enzo_EnzoMethodFire2Feedback.hpp b/src/Enzo/enzo_EnzoMethodFire2Feedback.hpp deleted file mode 100644 index 2886191ed1..0000000000 --- a/src/Enzo/enzo_EnzoMethodFire2Feedback.hpp +++ /dev/null @@ -1,59 +0,0 @@ -/// See LICENSE_CELLO file for license and copyright information - -/// @file enzo_EnzoMethodFire2Feedback.hpp -/// @author Andrew Emerick (aemerick11@gmail.com) -/// @date -/// @brief Implements the FIRE2 model for stellar Feedback -/// as implemented into Enzo from Azton Wells -/// - - -#ifdef NOTDEFINED - -#ifndef ENZO_ENZO_METHOD_FIRE2_FEEDBACK -#define ENZO_ENZO_METHOD_FIRE2_FEEDBACK - -class EnzoMethodFire2Feedback : public Method { - - /// @class EnzoMethodFire2Feedback - /// @ingroup Enzo - /// @btief [\ref Enzo] Encapsulate Feedback Routines - -public: - - EnzoMethodFire2Feedback(); - - /// Destructor - virtual ~EnzoMethodFire2Feedback() throw() {}; - - /// Charm++ Pup::able declarations - PUPable_decl(EnzoMethodFire2Feedback); - - /// Charm++ Pup::able migration Constructor - EnzoMethodFire2Feedback (CkMigrateMessage *m) - : Method (m) - { } - - /// Charm++ Pack / Unpack function - void pup(PUP::er &p); - - /// Apply the method - virtual void compute (Block * block) throw(); - - void compute_ (Block * block); - - /// name - virtual std::string name() throw() - { return "feedback"; } - - // Compute the maximum timestep for this method - virtual double timestep (Block * block) throw(); - -protected: - - int sf_minimum_level_; - int single_sn_; - - -} -#endif diff --git a/src/Enzo/enzo_EnzoMethodGrackle.cpp b/src/Enzo/enzo_EnzoMethodGrackle.cpp index ac874c4717..be4d7eee7c 100644 --- a/src/Enzo/enzo_EnzoMethodGrackle.cpp +++ b/src/Enzo/enzo_EnzoMethodGrackle.cpp @@ -43,7 +43,6 @@ EnzoMethodGrackle::EnzoMethodGrackle /// Define Grackle's internal data structures time_grackle_data_initialized_ = ENZO_FLOAT_UNDEFINED; initialize_grackle_chemistry_data(time); - #endif /* CONFIG_USE_GRACKLE */ } @@ -403,6 +402,7 @@ void EnzoMethodGrackle::update_grackle_density_fields( chemistry_data * grackle_chemistry = enzo::config()->method_grackle_chemistry; + double metallicity_floor_ = enzo::config()->method_grackle_metallicity_floor; for (int iz = 0; izadd_field("acceleration_x"); refresh->add_field("acceleration_y"); refresh->add_field("acceleration_z"); - // Accumulate is used when particles are deposited into density_total if (accumulate) { @@ -123,7 +122,6 @@ void EnzoMethodGravity::compute(Block * block) throw() // Initialize the linear system Field field = block->data()->field(); - /// access problem-defining fields for eventual RHS and solution const int ib = field.field_id ("B"); const int id = field.field_id("density"); @@ -140,7 +138,6 @@ void EnzoMethodGravity::compute(Block * block) throw() field.ghost_depth(0,&gx,&gy,&gz); const int m = mx*my*mz; - enzo_float * B = (enzo_float*) field.values (ib); #ifdef DEBUG_COPY_B const int ib_copy = field.field_id ("B_copy"); @@ -162,7 +159,6 @@ void EnzoMethodGravity::compute(Block * block) throw() if (block->is_leaf()) { if (cosmology) { - int gx,gy,gz; field.ghost_depth(0,&gx,&gy,&gz); gx=gy=gz=0; @@ -197,7 +193,6 @@ void EnzoMethodGravity::compute(Block * block) throw() const int ix = field.field_id ("potential"); std::shared_ptr A (std::make_shared(order_)); - solver->set_field_x(ix); solver->set_field_b(ib); #ifdef DEBUG_COPY_B @@ -208,7 +203,6 @@ void EnzoMethodGravity::compute(Block * block) throw() if (DT_copy) for (int i=0; iapply (A, block); } diff --git a/src/Enzo/enzo_EnzoMethodStarMaker.cpp b/src/Enzo/enzo_EnzoMethodStarMaker.cpp index 9816bd71ea..c7bea2b008 100644 --- a/src/Enzo/enzo_EnzoMethodStarMaker.cpp +++ b/src/Enzo/enzo_EnzoMethodStarMaker.cpp @@ -15,7 +15,7 @@ #include "cello.hpp" #include "enzo.hpp" -// #define DEBUG_SF +//#define DEBUG_SF //------------------------------------------------------------------- @@ -41,8 +41,8 @@ EnzoMethodStarMaker::EnzoMethodStarMaker // Copy over parameters from config to local names here for convenience use_density_threshold_ = enzo_config->method_star_maker_use_density_threshold; use_velocity_divergence_ = enzo_config->method_star_maker_use_velocity_divergence; - use_dynamical_time_ = enzo_config->method_star_maker_use_dynamical_time; use_self_gravitating_ = enzo_config->method_star_maker_use_self_gravitating; + use_altAlpha_ = enzo_config->method_star_maker_use_altAlpha; use_h2_self_shielding_ = enzo_config->method_star_maker_use_h2_self_shielding; use_jeans_mass_ = enzo_config->method_star_maker_use_jeans_mass; number_density_threshold_ = enzo_config->method_star_maker_number_density_threshold; @@ -50,6 +50,15 @@ EnzoMethodStarMaker::EnzoMethodStarMaker maximum_star_fraction_ = enzo_config->method_star_maker_maximum_mass_fraction; star_particle_min_mass_ = enzo_config->method_star_maker_minimum_star_mass; star_particle_max_mass_ = enzo_config->method_star_maker_maximum_star_mass; + use_dynamical_time_ = enzo_config->method_star_maker_use_dynamical_time; + + use_overdensity_threshold_ = enzo_config->method_star_maker_use_overdensity_threshold; + overdensity_threshold_ = enzo_config->method_star_maker_overdensity_threshold; + use_critical_metallicity_ = enzo_config->method_star_maker_use_critical_metallicity; + critical_metallicity_ = enzo_config->method_star_maker_critical_metallicity; + use_cooling_time_ = enzo_config->method_star_maker_use_cooling_time; + use_temperature_threshold_ = enzo_config->method_star_maker_use_temperature_threshold; + temperature_threshold_ = enzo_config->method_star_maker_temperature_threshold; } //------------------------------------------------------------------- @@ -64,7 +73,6 @@ void EnzoMethodStarMaker::pup (PUP::er &p) p | use_density_threshold_; p | use_velocity_divergence_; - p | use_dynamical_time_; p | number_density_threshold_; p | efficiency_; p | maximum_star_fraction_; @@ -73,6 +81,9 @@ void EnzoMethodStarMaker::pup (PUP::er &p) p | use_self_gravitating_; p | use_h2_self_shielding_; p | use_jeans_mass_; + p | use_cooling_time_; + p | use_overdensity_threshold_; + p | critical_metallicity_; return; } @@ -118,8 +129,7 @@ void EnzoMethodStarMaker::rescale_densities(EnzoBlock * enzo_block, for (int ic = 0; ic < nc; ic++){ enzo_float * cfield = (enzo_float *) field.values(field_groups->item("color",ic)); - - cfield[index] *= density_ratio; + cfield[index] *= density_ratio; } @@ -192,10 +202,34 @@ int EnzoMethodStarMaker::check_number_density_threshold( (d >= this->number_density_threshold_); } +int EnzoMethodStarMaker::check_overdensity_threshold(const double &rho) +{ + //#ifdef DEBUG_SF + // CkPrintf("MethodStarMaker -- overdensity = %f\n", rho); + //#endif + return !(this->use_overdensity_threshold_) + + (rho >= this->overdensity_threshold_); +} + + +int EnzoMethodStarMaker::check_self_gravitating_new(const double total_energy, const double potential) +{ + + if (!this->use_self_gravitating_) + return 1; + + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- alpha = %f\n",total_energy/potential); + #endif + + return (total_energy/potential < 1); + +} + int EnzoMethodStarMaker::check_self_gravitating( - const double mean_particle_mass, const double rho_cgs, const enzo_float temperature, + const double mean_particle_mass, const double density, const enzo_float temperature, enzo_float *vx, enzo_float *vy, enzo_float *vz, - const double lunit, const double vunit, + const double lunit, const double vunit, const double rhounit, const int &index, const int &dix, const int &diy, const int &diz, const double dx, const double dy, const double dz) { @@ -213,20 +247,24 @@ int EnzoMethodStarMaker::check_self_gravitating( // Frobenius norm of the velocity gradient tensor div_v_norm2 = (pow(vx[index+dix] - vx[index-dix], 2) + pow(vy[index+dix] - vy[index-dix], 2) + - pow(vz[index+dix] - vz[index-dix], 2)) / dx2 + + pow(vz[index+dix] - vz[index-dix], 2)) / (4*dx2) + (pow(vx[index+diy] - vx[index-diy], 2) + pow(vy[index+diy] - vy[index-diy], 2) + - pow(vz[index+diy] - vz[index-diy], 2)) / dy2 + + pow(vz[index+diy] - vz[index-diy], 2)) / (4*dy2) + (pow(vx[index+diz] - vx[index-diz], 2) + pow(vy[index+diz] - vy[index-diz], 2) + - pow(vz[index+diz] - vz[index-diz], 2)) / dz2; - div_v_norm2 *= (vunit*vunit); + pow(vz[index+diz] - vz[index-diz], 2)) / (4*dz2); + div_v_norm2 *= (vunit*vunit); // constant for testing. TODO: change to variable const double gamma = 5.0 / 3.0; cs2 = (gamma * enzo_constants::kboltz * temperature) / mean_particle_mass; - alpha = (div_v_norm2 + cs2/dx2) / (8 * cello::pi * enzo_constants::grav_constant * rho_cgs); + alpha = (div_v_norm2 + cs2/dx2) / (8 * cello::pi * enzo_constants::grav_constant * density*rhounit); + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- alpha = %f\n",alpha); + #endif + return (alpha < 1); } @@ -248,21 +286,25 @@ double EnzoMethodStarMaker::h2_self_shielding_factor( const double rho_cgs = rho[index] * dunit; - grad_rho = sqrt(pow((rho[index+dix] - rho[index-dix]) / dx, 2) + - pow((rho[index+diy] - rho[index-diy]) / dy, 2) + - pow((rho[index+diz] - rho[index-diz]) / dz, 2)); + grad_rho = sqrt(pow((rho[index+dix] - rho[index-dix]) / (2*dx), 2) + + pow((rho[index+diy] - rho[index-diy]) / (2*dy), 2) + + pow((rho[index+diz] - rho[index-diz]) / (2*dz), 2)); grad_rho *= dunit / lunit; - tau = 434.8 * rho_cgs * (dx + rho_cgs / grad_rho); // 434.8 cm^2 / g + tau = 434.8 * rho_cgs * (dx*lunit + rho_cgs / grad_rho); // 434.8 cm^2 / g phi = 0.756 * pow(1.0 + 3.1 * metallicity, 0.365); psi = (0.6 * tau * (0.01 + metallicity)) / (log(1.0 + 0.6*phi + 0.01*phi*phi)); f_shield = 1.0 - 3.0 / (1.0 + 4.0*psi); + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- f_shield = %f\n",f_shield); + #endif return f_shield; } int EnzoMethodStarMaker::check_jeans_mass( const double temperature, const double mean_particle_mass, - const double rho_cgs, const double mass + const double density, const double mass, + const double munit, const double rhounit ) { if (!use_jeans_mass_) @@ -271,15 +313,21 @@ int EnzoMethodStarMaker::check_jeans_mass( const double gamma = 5.0 / 3.0; const double minimum_jeans_mass = 1000 * enzo_constants::mass_solar; double cs2 = (gamma * enzo_constants::kboltz * temperature) / mean_particle_mass; - double m_jeans = (cello::pi/6) * pow(cs2, 1.5) / (pow(enzo_constants::grav_constant, 1.5) * sqrt(rho_cgs)); + + double m_jeans = (cello::pi/6) * pow(cs2, 1.5) / + (pow(enzo_constants::grav_constant, 1.5) * sqrt(density*rhounit)); + double m_jcrit = MAX(minimum_jeans_mass, m_jeans); - return (mass < m_jcrit); + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- jeans_mass = %f\n",m_jeans); + #endif + return (mass*munit < m_jcrit); } int EnzoMethodStarMaker::check_velocity_divergence( enzo_float *vx, enzo_float *vy, enzo_float *vz, const int &index, const int &dix, const int &diy, - const int &diz){ + const int &diz, const double dx, const double dy, const double dz) { /// Apply the criteria that the divergence of the velocity /// be negative, if so desired by user (use_velocity_divergence). @@ -288,21 +336,16 @@ int EnzoMethodStarMaker::check_velocity_divergence( return 1.0; } - int result = 0; + double div = 0.0; + if (vx) div += 0.5 * (vx[index+dix] - vx[index-dix]) / dx; // in units of dx + if (vy) div += 0.5 * (vy[index+diy] - vy[index-diy]) / dy; // in units of dy + if (vz) div += 0.5 * (vz[index+diz] - vz[index-diz]) / dz; // in units of dz - if(vx){ - result = (vx[index+dix] - vx[index-dix] < 0) ? 1 : 0; - } + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- velocity_divergence = %f\n",div); + #endif - if(vy && result){ - result = (vy[index+diy] - vy[index-diy] < 0) ? 1 : 0; - } - - if(vz && result){ - result = (vz[index+diz] - vz[index-diz] < 0) ? 1 : 0; - } - - return result; + return div < 0; } int EnzoMethodStarMaker::check_mass(const double &m){ @@ -315,3 +358,38 @@ int EnzoMethodStarMaker::check_mass(const double &m){ return minlimit; } + +int EnzoMethodStarMaker::check_cooling_time(const double &cooling_time,const double &total_density, + const double tunit, const double rhounit) +{ + // Check whether cooling_time < dynamical_time + if (!(this->use_cooling_time_)) { + return 1; + } + + double dynamical_time = pow(3.0*cello::pi/32.0/enzo_constants::grav_constant/(total_density*rhounit),0.5); //s + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- cooling_time = %f, dynamical_time = %f\n",cooling_time*tunit, dynamical_time); + #endif + return cooling_time*tunit < dynamical_time; + +} + +int EnzoMethodStarMaker::check_metallicity(const double &Z) +{ + // Enforce a critical metallicity for star formation + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- Z = %f, Zcrit = %f\n",Z, critical_metallicity_); + #endif + return !(this->use_critical_metallicity_) + + (Z >= critical_metallicity_); +} + +int EnzoMethodStarMaker::check_temperature(const double &T) +{ + #ifdef DEBUG_SF + CkPrintf("MethodStarMaker -- T = %f, Tcrit = %f\n", T, temperature_threshold_); + #endif + return !(this->use_temperature_threshold_) + + (T < temperature_threshold_); +} diff --git a/src/Enzo/enzo_EnzoMethodStarMaker.hpp b/src/Enzo/enzo_EnzoMethodStarMaker.hpp index 634fff1ea1..0eacb47ed8 100644 --- a/src/Enzo/enzo_EnzoMethodStarMaker.hpp +++ b/src/Enzo/enzo_EnzoMethodStarMaker.hpp @@ -68,18 +68,24 @@ class EnzoMethodStarMaker : public Method { // Routine functions for checking certain conditions // int check_number_density_threshold(const double &d); + int check_overdensity_threshold(const double &rho); int check_velocity_divergence( enzo_float *vx, enzo_float *vy, enzo_float *vz, const int &index, const int &dix, const int &diy, - const int &diz); + const int &diz, + const double dx, const double dy, const double dz); + int check_mass(const double &m); + int check_self_gravitating( - const double mean_particle_mass, const double rho_cgs, + const double mean_particle_mass, const double density, const enzo_float temperature, enzo_float *vx, enzo_float *vy, enzo_float *vz, - const double lunit, const double vunit, + const double lunit, const double vunit, const double rhounit, const int &index, const int &dix, const int &diy, const int &diz, const double dx, const double dy, const double dz); + + int check_self_gravitating_new(const double total_energy, const double potential); double h2_self_shielding_factor( enzo_float *rho, const double metallicity, const double dunit, const double lunit, @@ -87,22 +93,35 @@ class EnzoMethodStarMaker : public Method { const double dx, const double dy, const double dz); int check_jeans_mass( const double temperature, const double mean_particle_mass, - const double rho_cgs, const double mass); + const double density, const double mass, + const double munit, const double rhounit); + int check_cooling_time(const double &cooling_time, const double &total_density, + const double rhounit, const double tunit); + int check_metallicity(const double &Z); + int check_temperature(const double &T); protected: // attributes bool use_density_threshold_; bool use_velocity_divergence_; - bool use_dynamical_time_; bool use_self_gravitating_; + bool use_altAlpha_; bool use_h2_self_shielding_; bool use_jeans_mass_; + bool use_overdensity_threshold_; + bool use_critical_metallicity_; + bool use_temperature_threshold_; + bool use_cooling_time_; + bool use_dynamical_time_; + double overdensity_threshold_; + double critical_metallicity_; double number_density_threshold_; double efficiency_; double maximum_star_fraction_; double star_particle_min_mass_; double star_particle_max_mass_; + double temperature_threshold_; // variables to be passsed here }; diff --git a/src/Enzo/enzo_EnzoMethodStarMakerFIRE2.cpp b/src/Enzo/enzo_EnzoMethodStarMakerFIRE2.cpp deleted file mode 100644 index 5a4982c63b..0000000000 --- a/src/Enzo/enzo_EnzoMethodStarMakerFIRE2.cpp +++ /dev/null @@ -1,282 +0,0 @@ -/// See LICENSE_CELLO file for license and copyright information - -/// @file enzo_EnzoMethodStarMakerFIRE2.cpp -/// @author John Wise (jwise@physics.gatech.edu) -/// @date -/// @brief Implements the star maker class using FIRE-2 -/// -/// Derived star maker class that actually makes stars. This is -/// adapted from the algorithm described in Hopkins et al. -/// (2018, MNRAS, 480, 800) - -#include "cello.hpp" -#include "enzo.hpp" -#include - -// TODO (JHW, 27 Feb 2020) -// Copied from EnzoMethodStarMakerStochasticSF with no changes whatsoever. -// Plan is to make this a separate class that inherits from the stochastic -// algorithm. - -//------------------------------------------------------------------- - -EnzoMethodStarMakerStochasticSF::EnzoMethodStarMakerStochasticSF -() - : EnzoMethodStarMaker() -{ - // To Do: Make the seed an input parameter - //srand(time(NULL)); // need randum number generator for later - return; -} - -//------------------------------------------------------------------- - -void EnzoMethodStarMakerStochasticSF::pup (PUP::er &p) -{ - // NOTE: Change this function whenever attributes change - - TRACEPUP; - - EnzoMethodStarMaker::pup(p); // call parent class pup - - return; -} - -//------------------------------------------------------------------ - -void EnzoMethodStarMakerStochasticSF::compute ( Block *block) throw() -{ - - // Loop through the grid and check star formation criteria - // stochastically form stars if zone meets these criteria - - int count = 0; - - // Are we at the highest level? - if (! block->is_leaf()){ - block->compute_done(); - return; - } - - EnzoBlock * enzo_block = enzo::block(block); - const EnzoConfig * enzo_config = enzo::config(); - EnzoUnits * enzo_units = enzo::units(); - - - Particle particle = enzo_block->data()->particle(); - Field field = enzo_block->data()->field(); - - double dx, dy, dz; - block->cell_width(&dx, &dy, &dz); - - double lx, ly, lz; - block->lower(&lx,&ly,&lz); - - // declare particle position arrays - // default particle type is "star", but this will default - // to subclass particle_type - const int it = particle.type_index (this->particle_type()); - - const int ia_m = particle.attribute_index (it, "mass"); - const int ia_x = particle.attribute_index (it, "x"); - const int ia_y = particle.attribute_index (it, "y"); - const int ia_z = particle.attribute_index (it, "z"); - const int ia_vx = particle.attribute_index (it, "vx"); - const int ia_vy = particle.attribute_index (it, "vy"); - const int ia_vz = particle.attribute_index (it, "vz"); - - // additional particle attributes - const int ia_metal = particle.attribute_index (it, "metal_fraction"); - const int ia_to = particle.attribute_index (it, "creation_time"); - const int ia_l = particle.attribute_index (it, "lifetime"); - - int ib = 0; // batch counter - int ipp = 0; // counter - - /// pointers for particle attribut arrays (later) - enzo_float * pmass = 0; - enzo_float * px = 0; - enzo_float * py = 0; - enzo_float * pz = 0; - enzo_float * pvx = 0; - enzo_float * pvy = 0; - enzo_float * pvz = 0; - /// - enzo_float * pmetal = 0; - enzo_float * pform = 0; - enzo_float * plifetime = 0; - - // obtain the particle stride length - const int ps = particle.stride(it, ia_m); - - int rank = cello::rank(); - - int gx,gy,gz; - field.ghost_depth (0, &gx, &gy, &gz); - - int nx, ny, nz; - field.size ( &nx, &ny, &nz); - - int mx = nx + 2*gx; - int my = ny + 2*gy; - int mz = nz + 2*gz; - - - enzo_float * density = (enzo_float *) field.values("density"); - enzo_float * temperature = (enzo_float *) field.values("temperature"); - - enzo_float * velocity_x = (rank >= 1) ? - (enzo_float *)field.values("velocity_x") : NULL; - enzo_float * velocity_y = (rank >= 2) ? - (enzo_float *)field.values("velocity_y") : NULL; - enzo_float * velocity_z = (rank >= 3) ? - (enzo_float *)field.values("velocity_z") : NULL; - - enzo_float * metal = field.is_field("metal_density") ? - (enzo_float *) field.values("metal_density") : NULL; - - const double Zsolar = 0.02; // TODO: Update to more accurate value - - // Idea for multi-metal species - group these using 'group' - // class in IC parameter file and in SF / Feedback routines simply - // check if this group exists, and if it does, loop over all of these - // fields to assign particle chemical tags and deposit yields - - // compute the temperature (we need it here) - EnzoComputeTemperature compute_temperature - (enzo_config->ppm_density_floor, - enzo_config->ppm_temperature_floor, - enzo_config->ppm_mol_weight, - enzo_config->physics_cosmology); - - compute_temperature.compute(enzo_block); - - // iterate over all cells (not including ghost zones) - // - // To Do: Allow for multi-zone star formation by adding mass in - // surrounding cells if needed to accumulte enough mass - // to hit target star particle mass () - for (int iz=gz; izdensity(); - double mean_particle_mass = enzo_config->ppm_mol_weight * enzo_constants::mass_hydrogen; - double ndens = rho_cgs / mean_particle_mass; - - double mass = density[i] *dx*dy*dz * enzo_units->mass() / enzo_constants::mass_solar; - double metallicity = (metal) ? metal[i]/density[i]/Zsolar : 0.0; - - // - // Apply the criteria for star formation - // - if (! this->check_number_density_threshold(ndens)) continue; - // if (! this->check_self_gravitating(mean_particle_mass, rho_cgs, temperature[i], - // enzo_units->length(), enzo_units->density(), - // velocity_x, velocity_y, velocity_z, - // i, 1, my, my*mz, dx)) continue; - // if (! this->check_h2_self_shielding(density, metallicity, i, 1, my, my*mz, dx)) continue; - if (! this->check_velocity_divergence(velocity_x, velocity_y, - velocity_z, i, - 1, my, my*mz)) continue; - // Check whether mass in [min_mass, max_range] range and if specified, Jeans unstable - if (! this->check_mass(mass)) continue; - - double tdyn = sqrt(3.0 * cello::pi / 32.0 / enzo_constants::grav_constant / - (density[i] * enzo_units->density())); - - // - // compute fraction that can / will be converted to stars this step - // (just set to efficiency if dynamical time is ignored) - // - double star_fraction = this->use_dynamical_time_ ? - std::min(this->efficiency_ * enzo_block->dt * enzo_units->time() / tdyn, 1.0) : - this->efficiency_ ; - - // if this is less than the mass of a single particle, - // use a random number draw to generate the particle - if ( star_fraction * mass < this->star_particle_min_mass_){ - // get a random number - double rnum = (double(rand())) / (double(RAND_MAX)); - double probability = this->efficiency_ * mass / this->star_particle_min_mass_; - if (rnum > probability){ - continue; // do not form stars - } else{ - star_fraction = this->star_particle_min_mass_ / mass; - } - } else { - // (else, leave star fraction alone ) - star_fraction = this->star_particle_min_mass_ / mass; - } - - count++; // - - // now create a star particle - // insert_particles( particle_type, number_of_particles ) - int my_particle = particle.insert_particles(it, 1); - - // For the inserted particle, obtain the batch number (ib) - // and the particle index (ipp) - particle.index(my_particle, &ib, &ipp); - - int io = ipp; // ipp*ps - // pointer to mass array in block - pmass = (enzo_float *) particle.attribute_array(it, ia_m, ib); - - pmass[io] = star_fraction * (density[i] * dx * dy * dz); - px = (enzo_float *) particle.attribute_array(it, ia_x, ib); - py = (enzo_float *) particle.attribute_array(it, ia_y, ib); - pz = (enzo_float *) particle.attribute_array(it, ia_z, ib); - - // need to double check that these are correctly handling ghost zones - // I believe lx is lower coordinates of active region, but - // ix is integer index of whole grid (active + ghost) - // - px[io] = lx + (ix - gx + 0.5) * dx; - py[io] = ly + (iy - gy + 0.5) * dy; - pz[io] = lz + (iz - gz + 0.5) * dz; - - pvx = (enzo_float *) particle.attribute_array(it, ia_vx, ib); - pvy = (enzo_float *) particle.attribute_array(it, ia_vy, ib); - pvz = (enzo_float *) particle.attribute_array(it, ia_vz, ib); - - pvx[io] = velocity_x[i]; - if (velocity_y) pvy[io] = velocity_y[i]; - if (velocity_z) pvz[io] = velocity_z[i]; - - // finalize attributes - plifetime = (enzo_float *) particle.attribute_array(it, ia_l, ib); - pform = (enzo_float *) particle.attribute_array(it, ia_to, ib); - - pform[io] = enzo_block->time(); // formation time - plifetime[io] = 10.0 * enzo_constants::Myr_s / enzo_units->time() ; // lifetime - - if (metal){ - pmetal = (enzo_float *) particle.attribute_array(it, ia_metal, ib); - pmetal[io] = metal[i] / density[i]; - } - - // Remove mass from grid and rescale fraction fields - density[i] = (1.0 - star_fraction) * density[i]; - double scale = (1.0 - star_fraction) / 1.0; - - // rescale tracer fields to maintain constant mass fraction - // with the corresponding new density... - // scale = new_density / old_density - rescale_densities(enzo_block, i, scale); - } - } - } // end loop iz - - if (count > 0){ - std::cout << "Number of particles formed: " << count << "\n"; - } - - - block->compute_done(); - - return; -} diff --git a/src/Enzo/enzo_EnzoMethodStarMakerSTARSS.cpp b/src/Enzo/enzo_EnzoMethodStarMakerSTARSS.cpp new file mode 100644 index 0000000000..f809993e54 --- /dev/null +++ b/src/Enzo/enzo_EnzoMethodStarMakerSTARSS.cpp @@ -0,0 +1,565 @@ +/// See LICENSE_CELLO file for license and copyright information + +/// @file enzo_EnzoMethodStarMakerSTARSS.cpp +/// @author Will Hicks (whicks@ucsd.edu) +/// @date +/// @brief Implements the star maker class using STARSS algorithm developed +// by Azton Wells (TODO: link paper when it's published) +/// +/// Derived star maker class that actually makes stars. This is +/// adapted from the algorithm described in Hopkins et al. +/// (2018, MNRAS, 480, 800) + +#include "cello.hpp" +#include "enzo.hpp" +#include + +// #define DEBUG_SF_CRITERIA +// #define DEBUG_SF_CRITERIA_EXTRA +// #define DEBUG_STORE_INITIAL_PROPERTIES +//------------------------------------------------------------------- + +EnzoMethodStarMakerSTARSS::EnzoMethodStarMakerSTARSS +() + : EnzoMethodStarMaker() +{ + cello::simulation()->refresh_set_name(ir_post_,name()); + Refresh * refresh = cello::refresh(ir_post_); + refresh->add_all_fields(); + refresh->add_all_particles(); + + ParticleDescr * particle_descr = cello::particle_descr(); + + return; +} + +//------------------------------------------------------------------- + +void EnzoMethodStarMakerSTARSS::pup (PUP::er &p) +{ + // NOTE: Change this function whenever attributes change + + TRACEPUP; + + EnzoMethodStarMaker::pup(p); // call parent class pup + + return; +} + +//------------------------------------------------------------------ + +void EnzoMethodStarMakerSTARSS::compute ( Block *block) throw() +{ + + // Loop through the grid and check star formation criteria + // stochastically form stars if zone meets these criteria + + // TODO: Make random seed an input parameter + std::mt19937 mt(std::time(nullptr)); + int count = 0; + + const EnzoConfig * enzo_config = enzo::config(); + + // Are we at the highest level? + // Can we form stars at this level? + if ( (! block->is_leaf() ) || (block->level() < enzo_config->method_star_maker_min_level) ){ + block->compute_done(); + return; + } + + + EnzoBlock * enzo_block = enzo::block(block); + EnzoUnits * enzo_units = enzo::units(); + double lunit = enzo_units->length(); + double tunit = enzo_units->time(); + double vunit = enzo_units->velocity(); + double rhounit = enzo_units->density(); + double munit = enzo_units->mass(); + double munit_solar = munit / enzo_constants::mass_solar; + double eunit = vunit*vunit; // specific energy units + + Particle particle = enzo_block->data()->particle(); + Field field = enzo_block->data()->field(); + + double dx, dy, dz; + block->cell_width(&dx, &dy, &dz); + + const int rank = cello::rank(); + + double dt = enzo_block->dt; // timestep in code units + enzo_float cosmo_a = 1.0; + enzo_float cosmo_dadt = 0.0; + EnzoPhysicsCosmology * cosmology = enzo::cosmology(); + + double cell_volume = dx*dy*dz; + double lx, ly, lz; + block->lower(&lx,&ly,&lz); + + // declare particle position arrays + const int it = particle.type_index (this->particle_type()); + + const int ia_m = particle.attribute_index (it, "mass"); + const int ia_x = particle.attribute_index (it, "x"); + const int ia_y = particle.attribute_index (it, "y"); + const int ia_z = particle.attribute_index (it, "z"); + const int ia_vx = particle.attribute_index (it, "vx"); + const int ia_vy = particle.attribute_index (it, "vy"); + const int ia_vz = particle.attribute_index (it, "vz"); + + #ifdef DEBUG_STORE_INITIAL_PROPERTIES + const int ia_m_0 = particle.attribute_index (it, "mass_0"); + const int ia_x_0 = particle.attribute_index (it, "x_0"); + const int ia_y_0 = particle.attribute_index (it, "y_0"); + const int ia_z_0 = particle.attribute_index (it, "z_0"); + const int ia_vx_0 = particle.attribute_index (it, "vx_0"); + const int ia_vy_0 = particle.attribute_index (it, "vy_0"); + const int ia_vz_0 = particle.attribute_index (it, "vz_0"); + #endif + + // additional particle attributes + const int ia_metal = particle.attribute_index (it, "metal_fraction"); + const int ia_to = particle.attribute_index (it, "creation_time"); + const int ia_l = particle.attribute_index (it, "lifetime"); + const int ia_lev = particle.attribute_index (it, "creation_level"); + + int ib = 0; // batch counter + int ipp = 0; // counter + + /// pointers for particle attribute arrays (later) + enzo_float * pmass = 0; + enzo_float * px = 0; + enzo_float * py = 0; + enzo_float * pz = 0; + enzo_float * pvx = 0; + enzo_float * pvy = 0; + enzo_float * pvz = 0; + /// + enzo_float * pmetal = 0; + enzo_float * pform = 0; + enzo_float * plifetime = 0; + enzo_float * plevel = 0; + + int gx,gy,gz; + field.ghost_depth (0, &gx, &gy, &gz); + + int nx, ny, nz; + field.size ( &nx, &ny, &nz); + + int mx = nx + 2*gx; + int my = ny + 2*gy; + int mz = nz + 2*gz; + + // array increments + const int idx = 1; + const int idy = mx; + const int idz = mx*my; + + + enzo_float * density = (enzo_float *) field.values("density"); + enzo_float * temperature = (enzo_float *) field.values("temperature"); + + enzo_float * total_energy = (enzo_float *) field.values("total_energy"); + + enzo_float * potential = field.is_field("potential_copy") ? + (enzo_float *) field.values("potential_copy") : NULL; + + enzo_float * velocity_x = (rank >= 1) ? + (enzo_float *)field.values("velocity_x") : NULL; + enzo_float * velocity_y = (rank >= 2) ? + (enzo_float *)field.values("velocity_y") : NULL; + enzo_float * velocity_z = (rank >= 3) ? + (enzo_float *)field.values("velocity_z") : NULL; + + enzo_float * metal = field.is_field("metal_density") ? + (enzo_float *) field.values("metal_density") : NULL; + + enzo_float * cooling_time = field.is_field("cooling_time") ? + (enzo_float *) field.values("cooling_time") : NULL; + + enzo_float * density_particle_accumulate = field.is_field("density_particle_accumulate") ? + (enzo_float *) field.values("density_particle_accumulate") : NULL; + + enzo_float * dHI = field.is_field("HI_density") ? + (enzo_float*) field.values("HI_density") : NULL; + enzo_float * dHII = field.is_field("HII_density") ? + (enzo_float*) field.values("HII_density") : NULL; + enzo_float * dHeI = field.is_field("HeI_density") ? + (enzo_float*) field.values("HeI_density") : NULL; + enzo_float * dHeII = field.is_field("HeII_density") ? + (enzo_float*) field.values("HeII_density") : NULL; + enzo_float * dHeIII = field.is_field("HeIII_density") ? + (enzo_float*) field.values("HeIII_density") : NULL; + enzo_float * d_el = field.is_field("e_density") ? + (enzo_float*) field.values("e_density") : NULL; + + enzo_float * dH2I = field.is_field("H2I_density") ? + (enzo_float*) field.values("H2I_density") : NULL; + enzo_float * dH2II = field.is_field("H2II_density") ? + (enzo_float*) field.values("H2II_density") : NULL; + enzo_float * dHM = field.is_field("HM_density") ? + (enzo_float*) field.values("HM_density") : NULL; + + enzo_float * dDI = field.is_field("DI_density") ? + (enzo_float *) field.values("DI_density") : NULL; + enzo_float * dDII = field.is_field("DII_density") ? + (enzo_float *) field.values("DII_density") : NULL; + enzo_float * dHDI = field.is_field("HDI_density") ? + (enzo_float *) field.values("HDI_density") : NULL; + +/* + bool use_altAlpha = this->use_altAlpha_; + #ifndef DEBUG_COPY_POTENTIAL + if (use_altAlpha) { + CkPrintf("MethodStarMaker -- WARNING: use_altAlpha = true, " + " but DEBUG_COPY_POTENTIAL flag is not " + " set in EnzoMethodGravity.\n" + " Defaulting to use_altAlpha = false\n"); + use_altAlpha = false; + } + #endif +*/ + if (use_altAlpha_) { + // "potential_copy" field holds potential in proper coordinates, + // need to convert back to comoving. + if (cosmology) { + enzo_float cosmo_a = 1.0; + enzo_float cosmo_dadt = 0.0; + double dt = enzo_block->timestep(); + double time = enzo_block->time(); + cosmology-> compute_expansion_factor (&cosmo_a,&cosmo_dadt,time+0.5*dt); + for (int i=0; ippm_density_floor, + enzo_config->ppm_temperature_floor, + enzo_config->ppm_mol_weight, + enzo_config->physics_cosmology); + + compute_temperature.compute(enzo_block); + + + double mu = enzo_config->ppm_mol_weight; + // iterate over all cells (not including ghost zones) + for (int iz=gz; izmethod_grackle_chemistry)->primordial_chemistry; + if (primordial_chemistry > 0) { + mu = d_el[i] + dHI[i] + dHII[i] + 0.25*(dHeI[i]+dHeII[i]+dHeIII[i]); + + if (primordial_chemistry > 1) { + mu += dHM[i] + 0.5*(dH2I[i]+dH2II[i]); + } + if (primordial_chemistry > 2) { + mu += 0.5*(dDI[i] + dDII[i]) + dHDI[i]/3.0; + } + } + mu /= density[i]; + mu = 1.0/mu; + #endif + + double rho_cgs = density[i] * enzo_units->density(); + double mean_particle_mass = mu * enzo_constants::mass_hydrogen; + double ndens = rho_cgs / mean_particle_mass; + + double cell_mass = density[i] * cell_volume; + double metallicity = (metal) ? metal[i]/density[i]/enzo_constants::metallicity_solar : 0.0; + + // + // Apply the criteria for star formation + // + + // check number density threshold + if (! this->check_number_density_threshold(ndens) ) continue; + + // check overdensity threshold + // In cosmology, units are scaled such that mean(density) = 1, + // so density IS overdensity in these units + if (! this->check_overdensity_threshold(density[i]) ) continue; + + #ifdef DEBUG_SF_CRITERIA + CkPrintf("MethodFeedbackSTARSS -- overdensity threshold passed! rho=%f\n",dmean); + #endif + + // check velocity divergence < 0 + if (! this->check_velocity_divergence(velocity_x, velocity_y, + velocity_z, i, + idx, idy, idz, dx, dy, dz)) continue; + + #ifdef DEBUG_SF_CRITERIA_EXTRA + CkPrintf("MethodStarMakerSTARSS -- div(v) < 0 in cell %d\n", i); + #endif + + // check that alpha < 1 + if (use_altAlpha_) { + // NOTE: this requires the DEBUG_COPY_POTENTIAL flag to be set in EnzoMethodGravity + // TODO: Make saving potential an input parameter instead of debug flag? + if (! this->check_self_gravitating_new(total_energy[i], potential[i])) continue; + } + + else { + if (! this->check_self_gravitating(mean_particle_mass, density[i], temperature[i], + velocity_x, velocity_y, velocity_z, + lunit, vunit, rhounit, + i, idx, idy, idz, dx, dy, dz)) continue; + } + + #ifdef DEBUG_SF_CRITERIA_EXTRA + CkPrintf("MethodStarMakerSTARSS -- alpha < 1 in cell %d\n", i); + #endif + + // check that (T Tcrit + if (! this->check_temperature(temperature[i])) continue; + + if (cooling_time){ // if we are evolving a "cooling_time" field + if (! this->check_cooling_time(cooling_time[i], total_density, tunit, rhounit)) continue; + } + + // check that M > Mjeans + if (! check_jeans_mass(temperature[i], mean_particle_mass, density[i], cell_mass, + munit,rhounit )) continue; + + #ifdef DEBUG_SF_CRITERIA_EXTRA + CkPrintf("MethodStarMakerSTARSS -- M > M_jeans in cell %d\n", i); + #endif + + // check that H2 self shielded fraction f_shield > 0 + double f_shield = this->h2_self_shielding_factor(density,metallicity, + rhounit,lunit,i,idx,idy,idz,dx,dy,dz); + if (f_shield < 0) continue; + + // check that Z > Z_crit + if (! check_metallicity(metallicity)) continue; + +//-----------------------------------CREATION ROUTINE------------------------------- + + #ifdef DEBUG_SF_CRITERIA + CkPrintf("MethodStarMakerSTARSS -- SF criteria passed in cell %d\n", i); + #endif + + //free fall time in code units + double tff = sqrt(3*cello::pi/(32*enzo_constants::grav_constant*density[i]*rhounit))/tunit; + /* Determine Mass of new particle + WARNING: this removes the mass of the formed particle from the + host cell. If your simulation has very small (>15 Msun) baryon mass + per cell, it will break your sims! - AIW + */ + double divisor = std::max(1.0, tff * tunit/enzo_constants::Myr_s); + double maximum_star_mass = enzo_config->method_star_maker_maximum_star_mass; + double minimum_star_mass = enzo_config->method_star_maker_minimum_star_mass; + + if (maximum_star_mass < 0){ + maximum_star_mass = this->maximum_star_fraction_ * cell_mass * munit_solar; //Msun + } + + double bulk_SFR = f_shield * this->maximum_star_fraction_ * cell_mass*munit_solar/divisor; + + // Probability has the last word + // FIRE-2 uses p = 1 - exp (-MassShouldForm*dt / M_gas_particle) to convert a whole particle to star particle + // We convert a fixed portion of the baryon mass (or the calculated amount) + + double p_form = 1.0 - std::exp(-bulk_SFR*dt*(tunit/enzo_constants::Myr_s)/ + (this->maximum_star_fraction_*cell_mass*munit_solar)); + + if (enzo_config->method_star_maker_turn_off_probability) p_form = 1.0; + + double random = double(mt()) / double(mt.max()); + + /* New star is mass_should_form up to f_shield*maximum_star_fraction_ * baryon mass of the cell, + but at least 15 msun */ + double new_mass = std::min(f_shield*maximum_star_fraction_*cell_mass, maximum_star_mass/munit_solar); + + #ifdef DEBUG_SF_CRITERIA + CkPrintf("MethodStarMakerSTARSS -- new_mass = %f; p_form = %f\n", new_mass*munit_solar, p_form); + CkPrintf("MethodStarMakerSTARSS -- cell_mass = %f Msun; divisor = %f\n", cell_mass*munit_solar,divisor); + CkPrintf("MethodStarMakerSTARSS -- (ix, iy, iz) = (%d, %d, %d)\n", ix,iy,iz); + #endif + + if ( + (new_mass * munit_solar < minimum_star_mass) // too small + || (random > p_form) // too unlikely + || (new_mass > cell_mass) // too big compared to cell + ) + { + #ifdef DEBUG_SF_CRITERIA + CkPrintf("MethodStarMakerSTARSS -- star mass is either too big, too small, or failed the dice roll\n"); + #endif + continue; + } + + int n_newStars = 1; // track how many stars to form + double mFirmed = 0.0; // track total mass formed thus far + double massPerStar = new_mass; + double max_massPerStar = 5e3; // TODO: either make this a new parameter or replace 'maximum_star_mass' + double mass_split = 1e3; + if (new_mass * munit_solar > max_massPerStar) { // + // for large particles, split them into several 1e3-ish Msun particles + // that have slightly different birth times, + // spread over three dynamical times. This is to prevent huge particles suddenly dumping a HUGE + // amount of ionizing radiation at once. + // TODO: Enforce that if dt < 3 dynamical times, new star particles form in the next timestep + // Can do this by creating the particle now so we can still access it at the next timestep, + // but not assigning it any attributes until it's actually supposed to form + + n_newStars = std::floor(new_mass * munit_solar / mass_split); + massPerStar = new_mass / n_newStars; + #ifdef DEBUG_SF_CRITERIA + CkPrintf("MethodStarMakerSTARSS -- Predicted cluster mass %1.3e Msun > %1.3e Msun;\n" + " splitting into %d particles with mass %1.3e Msun\n", + new_mass*munit_solar, max_massPerStar, n_newStars, massPerStar*munit_solar); + #endif + } + for (int n=0; ntime(); + if (n > 0) { + double mod = n * 3.0 * tff/tunit/n_newStars; + ctime += mod; + } + count++; // time to form a star! + + #ifdef DEBUG_SF_CRITERIA + CkPrintf("MethodStarMakerSTARSS -- Forming star in gas with number density %f cm^-3\n", ndens); + #endif + + // now create a star particle + int my_particle = particle.insert_particles(it, 1); + + // For the inserted particle, obtain the batch number (ib) + // and the particle index (ipp) + particle.index(my_particle, &ib, &ipp); + + int io = ipp; // ipp*ps + // pointer to mass array in block + pmass = (enzo_float *) particle.attribute_array(it, ia_m, ib); + + pmass[io] = massPerStar; + px = (enzo_float *) particle.attribute_array(it, ia_x, ib); + py = (enzo_float *) particle.attribute_array(it, ia_y, ib); + pz = (enzo_float *) particle.attribute_array(it, ia_z, ib); + + // give it position at center of host cell + // TODO: Calculate CM instead? + px[io] = lx + (ix - gx + 0.5) * dx; + py[io] = ly + (iy - gy + 0.5) * dy; + pz[io] = lz + (iz - gz + 0.5) * dz; + + pvx = (enzo_float *) particle.attribute_array(it, ia_vx, ib); + pvy = (enzo_float *) particle.attribute_array(it, ia_vy, ib); + pvz = (enzo_float *) particle.attribute_array(it, ia_vz, ib); + + // average particle velocity over many cells to prevent runaway + double rhosum = 0.0; + double vx = 0.0; + double vy = 0.0; + double vz = 0.0; + for (int ix_=std::max(0,ix-3); ix_ <= std::min(ix+3,mx); ix_++) { + for (int iy_=std::max(0,iy-3); iy_ <= std::min(iy+3,my); iy_++) { + for (int iz_=std::max(0,iz-3); iz_ <= std::min(iz+3,mz); iz_++) { + int i_ = INDEX(ix_,iy_,iz_,mx,my); + vx += velocity_x[i_]*density[i_]; + vy += velocity_y[i_]*density[i_]; + vz += velocity_z[i_]*density[i_]; + rhosum += density[i_]; + } + } + } + vx /= rhosum; + vy /= rhosum; + vz /= rhosum; + + // TODO: Make this an input parameter + double max_velocity = 150e5/vunit; + if (std::abs(vx) > max_velocity) vx = vx/std::abs(vx) * max_velocity; + if (std::abs(vy) > max_velocity) vy = vy/std::abs(vy) * max_velocity; + if (std::abs(vz) > max_velocity) vz = vz/std::abs(vz) * max_velocity; + + pvx[io] = vx; + pvy[io] = vy; + pvz[io] = vz; + + // finalize attributes + plifetime = (enzo_float *) particle.attribute_array(it, ia_l, ib); + pform = (enzo_float *) particle.attribute_array(it, ia_to, ib); + plevel = (enzo_float *) particle.attribute_array(it, ia_lev, ib); + + pform[io] = ctime; // formation time + + //TODO: Need to have some way of calculating lifetime based on particle mass + plifetime[io] = 25.0 * enzo_constants::Myr_s / enzo_units->time() ; // lifetime (not accessed for STARSS FB) + + plevel[io] = enzo_block->level(); // formation level + + if (metal){ + pmetal = (enzo_float *) particle.attribute_array(it, ia_metal, ib); + pmetal[io] = metal[i] / density[i]; // in ABSOLUTE units + } + + // Remove mass from grid and rescale fraction fields + // TODO: If particle position is updated to CM instead of being cell-centered, will have to + // remove mass using CiC, which could complicate things because that CiC cloud could + // leak into the ghost zones. Would have to use same refresh+accumulate machinery + // as MethodFeedbackSTARSS to account for this. + double scale = (1.0 - pmass[io] / cell_mass); + density[i] *= scale; + // rescale color fields too + this->rescale_densities(enzo_block, i, scale); + + #ifdef DEBUG_STORE_INITIAL_PROPERTIES + enzo_float * pmass0 = (enzo_float *) particle.attribute_array(it, ia_m_0 , ib); + enzo_float * px0 = (enzo_float *) particle.attribute_array(it, ia_x_0 , ib); + enzo_float * py0 = (enzo_float *) particle.attribute_array(it, ia_y_0 , ib); + enzo_float * pz0 = (enzo_float *) particle.attribute_array(it, ia_z_0 , ib); + enzo_float * pvx0 = (enzo_float *) particle.attribute_array(it, ia_vx_0, ib); + enzo_float * pvy0 = (enzo_float *) particle.attribute_array(it, ia_vy_0, ib); + enzo_float * pvz0 = (enzo_float *) particle.attribute_array(it, ia_vz_0, ib); + + pmass0[io] = pmass[io]; + px0 [io] = px[io]; + py0 [io] = py[io]; + pz0 [io] = pz[io]; + pvx0[io] = pvx[io]; + pvy0[io] = pvy[io]; + pvz0[io] = pvz[io]; + #endif + + } // end loop through particles created in this cell + + + } + } + } // end loop iz + + if (count > 0){ + CkPrintf("MethodStarMakerSTARSS -- Number of particles formed: %d\n", count); + } + +/* TODO: Add this part in once Pop III SF/FB is implemented + + if (gridShouldFormStars && MechStarsSeedField && (nCreated == 0)){ + // set off a p3 supernova at at the last cell that could + // host star formation in the grid if the + // grid can host star formation but has no appreciable metals + fprintf(stdout, "\n\n\n[%d] %d %d %d Creating seed field!\n\n\n\n", + ID,seedIndex[0], seedIndex[1], seedIndex[2]) ; + MechStars_SeedSupernova(&totalMetal[0], Temperature, seedIndex); + + } +*/ + + block->compute_done(); + + return; +} diff --git a/src/Enzo/enzo_EnzoMethodStarMakerSTARSS.hpp b/src/Enzo/enzo_EnzoMethodStarMakerSTARSS.hpp new file mode 100644 index 0000000000..5241608659 --- /dev/null +++ b/src/Enzo/enzo_EnzoMethodStarMakerSTARSS.hpp @@ -0,0 +1,49 @@ +/// + +/// +/// +/// +/// + +#ifndef ENZO_ENZO_METHOD_STARMAKER_STARSS +#define ENZO_ENZO_METHOD_STARMAKER_STARSS + +/// +/// +/// +/// + +class EnzoMethodStarMakerSTARSS : public EnzoMethodStarMaker { + +public: + // Create new EnzoStarMakerSTARSS object + EnzoMethodStarMakerSTARSS(); + + // Charm++ PUP::able declarations + PUPable_decl(EnzoMethodStarMakerSTARSS); + + // Charm++ PUP::able declarations + EnzoMethodStarMakerSTARSS (CkMigrateMessage *m) + : EnzoMethodStarMaker (m) + { } + + /// Charm++ Pack / Unpack function + void pup (PUP::er &p); + + /// Apply method + virtual void compute ( Block * block) throw(); + + virtual std::string particle_type () throw() + { return "star";} + + /// Name + virtual std::string name () throw() + { return "star_maker";} + + virtual ~EnzoMethodStarMakerSTARSS() throw() {}; + +protected: + +}; + +#endif /* EnzoMethodStarMakerSTARSS */ diff --git a/src/Enzo/enzo_EnzoMethodStarMakerStochasticSF.cpp b/src/Enzo/enzo_EnzoMethodStarMakerStochasticSF.cpp index c33b2e50f6..c265df91f3 100644 --- a/src/Enzo/enzo_EnzoMethodStarMakerStochasticSF.cpp +++ b/src/Enzo/enzo_EnzoMethodStarMakerStochasticSF.cpp @@ -177,8 +177,9 @@ void EnzoMethodStarMakerStochasticSF::compute ( Block *block) throw() if (! this->check_number_density_threshold(ndens)) continue; if (! this->check_self_gravitating( mean_particle_mass, rho_cgs, temperature[i], velocity_x, velocity_y, velocity_z, - enzo_units->length(), enzo_units->density(), - i, 1, my, my*mz, dx, dy, dz)) continue; + enzo_units->length(), enzo_units->velocity(), + enzo_units->density(), + i, 1, mx, mx*my, dx, dy, dz)) continue; // AJE: TO DO --- // If Grackle is used, check for this and use the H2 @@ -190,13 +191,13 @@ void EnzoMethodStarMakerStochasticSF::compute ( Block *block) throw() metallicity, enzo_units->density(), enzo_units->length(), - i, 1, my, my*mz, + i, 1, mx, mx*my, dx, dy, dz); mass *= f_h2; // apply correction (f_h2 = 1 if not used) if (! this->check_velocity_divergence(velocity_x, velocity_y, velocity_z, i, - 1, my, my*mz)) continue; + 1, mx, mx*my, dx, dy, dz)) continue; // Check whether mass in [min_mass, max_range] range and if specified, Jeans unstable if (! this->check_mass(mass)) continue; diff --git a/src/Enzo/enzo_EnzoProblem.cpp b/src/Enzo/enzo_EnzoProblem.cpp index 1068addb84..4816e060b5 100644 --- a/src/Enzo/enzo_EnzoProblem.cpp +++ b/src/Enzo/enzo_EnzoProblem.cpp @@ -702,6 +702,9 @@ Method * EnzoProblem::create_method_ // should generalize this to enable multiple maker types if (enzo_config->method_star_maker_flavor == "stochastic"){ method = new EnzoMethodStarMakerStochasticSF(); + } else if (enzo_config->method_star_maker_flavor == "STARSS" || + enzo_config->method_star_maker_flavor == "starss") { + method = new EnzoMethodStarMakerSTARSS(); } else{ // does not do anything method = new EnzoMethodStarMaker(); } @@ -709,7 +712,14 @@ Method * EnzoProblem::create_method_ } else if (name == "feedback") { // need a similar type swtich as in star maker - method = new EnzoMethodDistributedFeedback(); + if (enzo_config->method_feedback_flavor == "distributed"){ + method = new EnzoMethodDistributedFeedback(); + } else if (enzo_config->method_feedback_flavor == "STARSS" || + enzo_config->method_feedback_flavor == "starss") { + method = new EnzoMethodFeedbackSTARSS(); + } else { // does not do anything + method = new EnzoMethodFeedback(); + } } else if (name == "merge_sinks") { diff --git a/src/Enzo/enzo_EnzoSolverMg0.cpp b/src/Enzo/enzo_EnzoSolverMg0.cpp index c36d6ed176..cb5afb2130 100644 --- a/src/Enzo/enzo_EnzoSolverMg0.cpp +++ b/src/Enzo/enzo_EnzoSolverMg0.cpp @@ -231,7 +231,6 @@ void EnzoSolverMg0::apply ( std::shared_ptr A, Block * block) throw() Sync * sync_prolong = psync_prolong(block); sync_prolong->set_stop(1 + 1); // self and parent - enter_solver_ (enzo_block); } @@ -281,7 +280,6 @@ void EnzoSolverMg0::enter_solver_ (EnzoBlock * enzo_block) throw() enzo_block->contribute(2*sizeof(long double), &reduce, sum_long_double_2_type, callback); - } else { SOLVER_CONTROL(enzo_block,"min","max","2 calling begin_solve_2 (no shift)"); @@ -316,7 +314,6 @@ void EnzoSolverMg0::compute_shift_ void EnzoBlock::r_solver_mg0_begin_solve(CkReductionMsg* msg) { performance_start_(perf_compute,__FILE__,__LINE__); - static_cast (solver())->begin_solve(this,msg); performance_stop_(perf_compute,__FILE__,__LINE__);