Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

STARSS Star Formation and Feedback #198

Merged
merged 82 commits into from
Jun 30, 2022
Merged

Conversation

WillHicks96
Copy link
Contributor

@WillHicks96 WillHicks96 commented Jun 10, 2022

STARSS = Scale-intelligent Terminal-momentum Algorithm for Realistic Stellar Sources
The STARSS algorithm is being actively developed by Azton Wells. Here's a link to his open PR in enzo-dev that has some more info -> enzo-project/enzo-dev#145

This PR adds two new methods EnzoMethodStarMakerSTARSS and EnzoMethodFeedbackSTARSS, along with some changes and bugfixes to the EnzoMethodStarMaker class. This depends on PR #143.

This is the final form of the EnzoMethodFire2<StarMaker/Feedback> methods that Andrew Emerick began working on, and I was able to reuse some of his code. There are still some miscellaneous TODO`s in the code (e.g. giving parameters better names and adding more tests), but I think that this is ready for review.

I've so far added three files to the input directory:

  1. input/STARSS/SF_FB.incl - This outlines all of the relevant parameters, in the configuration that I've found works best in the cosmology tests I've done so far
  2. input/STARSS/method_feedback_starss.in - This is an ideal test that sets off a single supernova on the corner of a block. This is meant to test that the supernova mass/energy/momentum remains consistent with expectation when the initial deposition crosses block boundaries
  3. input/test_cosmo-dd-SF_FB.in - This is an example that shows how you would call SF+FB in a cosmology simulation. This is a pretty small test, so it won't form its first star particle until around redshift 7.5ish (~cycle 270) with the current parameters. This will take 15-20 minutes on 8 processors. A cosmology test that forms stars more quickly would probably be necessary at some point -- maybe we could have one that restarts from a checkpoint at a redshift where star formation is about to happen.

Most of my testing has been within the context of cosmology simulations. Here is a snapshot from a recent 256^3 simulation that shows SF+FB doing things:

with_stars_Projection_z_density_density
ENZOE_N256_SF_FB_lev4-009600_Projection_z_metallicity_density

William Hicks and others added 30 commits January 25, 2022 13:35
  - seems to work if deposition doesn't leak into ghost zones (need to add refresh)
…deposition still slightly wrong for particles that explode right on the edge of a block.
…I think)

Long-term evolution still looks slightly asymmetric though
Code currently hangs up if feedback isn't the last method in method:list (with refresh enabled).
Ideal test not conserving momentum across block boundaries (might not be feedback related)
Other than that, seems to be fully working in cosmo sims
…derdeposition of mass and metals

Added tool for generating .block_list file with new_output
…unding cells

Changed to be consistent with Enzo version, which also made this change recently
…bug where I could only call feedback as the last method.

Storing MMW as double instead of int....................................................................
…gy/potential

back to saving particle masses as mass
Added debug statement
Moved initialization of deposit fields as `tiny_number` to compute_() to fix bug where if multiple supernovae/stellar wind events occur near a block boundary in a cycle, ghost zone values would get overwritten each time an event occurs for that cycle instead of just adding to a running total of energy, momentum, etc. to deposit.
…kSTARSS

moved compute_done() back to compute(). Cosmology simulations are hanging when it's in p_method_feedback_starss_end() for some reason
Copy link
Contributor

@stefanarridge stefanarridge left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is my last day, I'll approve this PR so as not to hold it up. You can decide whether or not to implement the changes I suggested.

@tumlinson
Copy link
Contributor

Hi @WillHicks96 -- the isolated SN test is working great. I have scaled it up to 512^3 and 256 proc on Pleiades with no problems. The cosmology test also starts and has run to z= 8. However, when I try to open the outputs I get a mysterious yt error. I so am uncertain whether this relates to the code of the PR itself or whether its independent. If you recommend that this is separate, I'm ready to approve the PR. Thanks.

In [1]: import yt
/home1/jtumlins/anaconda3/lib/python3.9/site-packages/yt/utilities/logger.py:4: VisibleDeprecationWarning: The configuration file /home1/jtumlins/.config/yt/ytrc is deprecated in favor of /home1/jtumlins/.config/yt/yt.toml. Currently, both are present. Please manually remove the deprecated one to silence this warning.
Deprecated since v4.0.0. This feature will be removed in v4.1.0
from yt.config import ytcfg

In [2]: ds = yt.load('Dir_COSMO_SF_FB_0100/Dir_COSMO_SF_FB_0100.block_list')
yt : [INFO ] 2022-06-24 06:55:35,070 Parameters: current_time = 7.555469551567898
yt : [INFO ] 2022-06-24 06:55:35,070 Parameters: domain_dimensions = [32 32 32]
yt : [INFO ] 2022-06-24 06:55:35,071 Parameters: domain_left_edge = [0. 0. 0.]
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: domain_right_edge = [1. 1. 1.]
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: cosmological_simulation = 1
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: current_redshift = 21.687504086446236
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: omega_lambda = 0.7
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: omega_matter = 0.3
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: omega_radiation = 0.0
yt : [INFO ] 2022-06-24 06:55:35,072 Parameters: hubble_constant = 0.7

In [3]: p = yt.ProjectionPlot(ds, 'x', 'density')
Parsing Hierarchy: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 73/73 [00:00<00:00, 4725.50it/s]

RuntimeError Traceback (most recent call last)
in
----> 1 p = yt.ProjectionPlot(ds, 'x', 'density')

~/anaconda3/lib/python3.9/site-packages/yt/visualization/plot_window.py in init(self, ds, axis, fields, center, width, axes_unit, weight_field, max_level, origin, right_handed, fontsize, field_parameters, data_source, method, proj_style, window_size, buff_size, aspect)
1862 if weight_field is None and method == "integrate":
1863 self.projected = True
-> 1864 (bounds, center, display_center) = get_window_parameters(
1865 axis, center, width, ds
1866 )

~/anaconda3/lib/python3.9/site-packages/yt/visualization/plot_window.py in get_window_parameters(axis, center, width, ds)
76
77 def get_window_parameters(axis, center, width, ds):
---> 78 width = ds.coordinates.sanitize_width(axis, width, None)
79 center, display_center = ds.coordinates.sanitize_center(center, axis)
80 xax = ds.coordinates.x_axis[axis]

~/anaconda3/lib/python3.9/site-packages/yt/geometry/coordinates/coordinate_handler.py in sanitize_width(self, axis, width, depth)
215 if width is None:
216 # initialize the index if it is not already initialized
--> 217 self.ds.index
218 # Default to code units
219 if not is_sequence(axis):

~/anaconda3/lib/python3.9/site-packages/yt/data_objects/static_output.py in index(self)
537 oldsettings = np.geterr()
538 np.seterr(all="ignore")
--> 539 self.create_field_info()
540 np.seterr(**oldsettings)
541 return self._instantiated_index

~/anaconda3/lib/python3.9/site-packages/yt/data_objects/static_output.py in create_field_info(self)
590 mylog.debug("Creating Particle Union 'all'")
591 pu = ParticleUnion("all", list(self.particle_types_raw))
--> 592 nfields = self.add_particle_union(pu)
593 if nfields == 0:
594 mylog.debug("zero common fields: skipping particle union 'all'")

~/anaconda3/lib/python3.9/site-packages/yt/data_objects/static_output.py in add_particle_union(self, union)
733 # ...if we can't find them, we set them up as defaults.
734 new_fields = self._setup_particle_types([union.name])
--> 735 self.field_info.find_dependencies(new_fields)
736 return len(new_fields)
737

~/anaconda3/lib/python3.9/site-packages/yt/fields/field_info_container.py in find_dependencies(self, loaded)
414
415 def find_dependencies(self, loaded):
--> 416 deps, unavailable = self.check_derived_fields(loaded)
417 self.ds.field_dependencies.update(deps)
418 # Note we may have duplicated

~/anaconda3/lib/python3.9/site-packages/yt/fields/field_info_container.py in check_derived_fields(self, fields_to_check)
611 try:
612 # fd: field detector
--> 613 fd = fi.get_dependencies(ds=self.ds)
614 except blacklist as err:
615 print(f"{err.class} raised for field {field}")

~/anaconda3/lib/python3.9/site-packages/yt/fields/derived_field.py in get_dependencies(self, *args, **kwargs)
251 """
252 e = FieldDetector(*args, **kwargs)
--> 253 e[self.name]
254 return e
255

~/anaconda3/lib/python3.9/site-packages/yt/fields/field_detector.py in missing(self, item)
123 vv = finfo(self)
124 if not permute_params:
--> 125 vv = finfo(self)
126 except NeedsGridType as exc:
127 ngz = exc.ghost_zones

~/anaconda3/lib/python3.9/site-packages/yt/fields/derived_field.py in call(self, data)
292 )
293 with self.unit_registry(data):
--> 294 dd = self._function(self, data)
295 for field_name in data.keys():
296 if field_name not in original_fields:

~/anaconda3/lib/python3.9/site-packages/yt/fields/particle_fields.py in _deposit_field(field, data)
162 pos = data[ptype, "particle_position"]
163 # Get back into density
--> 164 pden = data[ptype, "particle_mass"]
165 top = data.deposit(pos, [pden * data[(ptype, fname)]], method=method)
166 bottom = data.deposit(pos, [pden], method=method)

~/anaconda3/lib/python3.9/site-packages/yt/fields/field_detector.py in missing(self, item)
123 vv = finfo(self)
124 if not permute_params:
--> 125 vv = finfo(self)
126 except NeedsGridType as exc:
127 ngz = exc.ghost_zones

~/anaconda3/lib/python3.9/site-packages/yt/fields/derived_field.py in call(self, data)
292 )
293 with self.unit_registry(data):
--> 294 dd = self._function(self, data)
295 for field_name in data.keys():
296 if field_name not in original_fields:

~/anaconda3/lib/python3.9/site-packages/yt/fields/particle_fields.py in _cat_field(field, data)
891
892 def _cat_field(field, data):
--> 893 return uconcatenate(
894 [data[dep_type, field_name] for dep_type in data.ds.particle_types_raw]
895 )

~/anaconda3/lib/python3.9/site-packages/unyt/array.py in uconcatenate(arrs, axis)
2099 """
2100 v = np.concatenate(arrs, axis=axis)
-> 2101 v = _validate_numpy_wrapper_units(v, arrs)
2102 return v
2103

~/anaconda3/lib/python3.9/site-packages/unyt/array.py in _validate_numpy_wrapper_units(v, arrs)
2077 a1 = arrs[0]
2078 if not all(a.units == a1.units for a in arrs[1:]):
-> 2079 raise RuntimeError("Your arrays must have identical units.")
2080 v.units = a1.units
2081 return v

RuntimeError: Your arrays must have identical units.

@WillHicks96
Copy link
Contributor Author

Thanks, @stefanarridge! I think I'm almost there. There is still some more to do over the next few days before this is ready to merge, mainly finishing the documentation and some more general cleaning.

@WillHicks96
Copy link
Contributor Author

WillHicks96 commented Jun 24, 2022

@tumlinson Thanks for trying out the tests! I've seen that yt error before, and it seems to be fixed in the latest dev version of yt. So if you update your yt, you should hopefully be able to access the data. I notice you're using a version of yt installed with conda.

The error happens when yt is trying to create the ("all", "particle_mass") union field, so I'm guessing the issue popped up in yt when PR #89 merged.

EDIT: Looks like the issue actually wasn't fixed in yt, and that I just commented out the relevant line in yt when I was frantically preparing for the Enzo-E workshop and then immediately forgot about it. I'll post an issue on the yt github. From what I can tell though, this is a yt issue, and not specifically an issue with this PR

@tumlinson
Copy link
Contributor

tumlinson commented Jun 24, 2022

I am on yt 4.1.dev0 and still getting the error, but since this is a yt thing and the Enzo-E tests run I can go ahead and approve.

@tumlinson tumlinson assigned tumlinson and unassigned tumlinson Jun 24, 2022
No longer need to specify them in a parameter file
@tumlinson tumlinson self-requested a review June 24, 2022 23:45
Copy link
Contributor

@tumlinson tumlinson left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

approved based on test simulations running successfully on Pleiades.

@jwise77
Copy link
Contributor

jwise77 commented Jun 29, 2022

The code is failing to compile now.

Removed duplicated code in SolveHydroEquations
@WillHicks96
Copy link
Contributor Author

WillHicks96 commented Jun 29, 2022

@jwise77 Think I found the issue. std::max didn't like being fed a value of type enzo_float, which caused the failure in single precision. It was compiling fine for me in double precision though.

William Hicks and others added 5 commits June 29, 2022 18:00
…stricted_sn=false, it will only set off a Type IA supernova if it's been determined that no Type II supernovae will go off. How it was written before, no Type IA supernovae would go off if unrestricted_sn=false, which is not the intended behavior.
…nergy formalism.

Mainly because it hasn't been tested without dual energy formalism at all
@peeples
Copy link
Contributor

peeples commented Jun 30, 2022

@WillHicks96 this looks ready to merge; please do so if there isn't anything left to edit!

@peeples peeples merged commit 0d20e20 into enzo-project:main Jun 30, 2022
@WillHicks96 WillHicks96 deleted the starss-merge branch April 12, 2023 19:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants