Generating Metallicity Field

[Ryan Tanner]

I am just learning how to use trident. I am loading Athena++ data to use with Trident. The Athena hdf5 files don't include metallicity as a field, unlike the enzo example files. How would I generate a simple solar metallicity field needed for make_simple_ray?

Thanks,

Ryan Tanner

[Cameron Hummels]

Hi Ryan,

Welcome to Trident!  

One thing that you could do is just create a dummy metallicity field with all the same metallicity value for your dataset.  This could be done by following the instructions for "Creating a derived field" in the yt documentation: https://yt-project.org/docs/dev/developing/creating_derived_fields.html .

Specifically, I think you could just add this to the front of your code before you load your dataset:

def _metallicity(field, data):
    factor = 0.3
    return (
        data.ds.arr(data["index", "ones"].d, "Zsun") * factor
    )

yt.add_field(
    ("gas", "metallicity"),
    function=_metallicity,
    sampling_type="local",
    units="Zsun",
)

That will add a metallicity field if one does not exist for any dataset loaded in the same python session.  By default, I've set "factor = 0.3" to set the metallicity to 0.3 solar everywhere, but you could change it to whatever you want depending on the nature of your problem.  To make sure this is working, you could just print out the metallicity field after you've loaded the datset with:

ds = yt.load('your_dataset')
print(ds.r['gas', 'metallicity'])

Let me know if this works for you.  

Cameron

--
You received this message because you are subscribed to the Google Groups "trident-project-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to trident-project-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/trident-project-users/a5ae2b2a-8a5f-4125-8535-7f9a5bc02625n%40googlegroups.com.

--

Cameron Hummels

Computational Astrophysicist

California Institute of Technology

http://chummels.org

[Ryan Tanner]

Thanks, that worked! But then I immediately ran into another problem. I realized that my version of yt wasn't up to date so I updated yt (which created its own mess). Now with yt version 4.0.1 and trident version 1.2.3 I try importing trident and I get:

>>> import trident

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/trident/__init__.py", line 41, in <module>

    from trident.instrument import \

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/trident/instrument.py", line 17, in <module>

    from trident.absorption_spectrum.absorption_spectrum import \

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/trident/absorption_spectrum/absorption_spectrum.py", line 26, in <module>

    from yt.convenience import load

ModuleNotFoundError: No module named 'yt.convenience'

[Cameron Hummels]

Currently there is an issue where the stable version of Trident (v1.2.3) won't work with the stable version of yt because the `yt.load()` command moved from convenience.py to loaders.py.  This has been the case for the last month or so since the new yt release of yt4.0.  I need to release a new stable version to address this.  I'll try to get to this ASAP, but for now, your best bet is to just use the development version of Trident: https://trident.readthedocs.io/en/latest/installation.html#step-2-install-trident

I hope this works for you!  Sorry about all of the problems.  I saw you were having issues with the new cmaps stuff in yt too.  I hope Clement got that fixed for you.

Cameron

To view this discussion on the web visit https://groups.google.com/d/msgid/trident-project-users/6a40df53-8adb-40c4-a097-ea2e3ed97726n%40googlegroups.com.

[Ryan Tanner]

I got the previous problem resolved and instantly hit another one. Sorry, I feel like I am hitting every single obscure problem possible. 

I can run the example with the enzo data and it works fine, and I can run the same thing with my athena++ data, but only with hydrogen lines. If I include any other atoms I get this error after it finishes calculating the hydrogen lines:

-----------------

Traceback (most recent call last):

  File "<ipython-input-2-5a28162496c0>", line 1, in <module>

    runfile('/mnt/c/Users/rtanner2/OneDrive - NASA/research/python/trident/single_ray_test.py', wdir='/mnt/c/Users/rtanner2/OneDrive - NASA/research/python/trident')

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py", line 827, in runfile

    execfile(filename, namespace)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py", line 110, in execfile

    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/mnt/c/Users/rtanner2/OneDrive - NASA/research/python/trident/single_ray_test.py", line 48, in <module>

    sg.make_spectrum(ray, lines=line_list)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/trident/trident/spectrum_generator.py", line 425, in make_spectrum

    min_tau=min_tau, njobs=njobs)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/trident/trident/absorption_spectrum/absorption_spectrum.py", line 493, in make_spectrum

    min_tau=min_tau, njobs=njobs)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/trident/trident/absorption_spectrum/absorption_spectrum.py", line 709, in _add_lines_to_spectrum

    column_density = field_data[line['field_name']] * field_data['dl']

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/data_containers.py", line 258, in __getitem__

    self.get_data(f)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/selection_objects/data_selection_objects.py", line 211, in get_data

    self._generate_fields(fields_to_generate)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/selection_objects/data_selection_objects.py", line 232, in _generate_fields

    fd = self._generate_field(field)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/data_containers.py", line 299, in _generate_field

    tr = self._generate_fluid_field(field)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/data_containers.py", line 318, in _generate_fluid_field

    rv = finfo(gen_obj)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/fields/derived_field.py", line 297, in __call__

    dd = self._function(self, data)

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/fields/geometric_fields.py", line 104, in _ones

    tmp = np.ones(data.ires.shape, dtype="float64")

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/selection_objects/data_selection_objects.py", line 424, in ires

    return self._current_chunk.ires

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/geometry/geometry_handler.py", line 255, in cached_func

    tr = self._accumulate_values(n[1:])

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/geometry/geometry_handler.py", line 292, in _accumulate_values

    arrs.append(f(self.dobj))

  File "/home/rtanner2/anaconda3/lib/python3.7/site-packages/yt/data_objects/index_subobjects/particle_container.py", line 17, in _func_non_indexed

    raise YTNonIndexedDataContainer(self)

YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type.  Operations such as ires, icoords, fcoords and fwidth will not work on it.

-----------------

Drummond Fielding recommended Trident to me and from what I have seen so far with it it looks great. I'm sorry I'm having such a hard time with it.

[Cameron Hummels]

Hi Ryan,

Thanks for being persistent with all these hurdles--clearly I've got some work to do to clean up the code! :)

Would you mind sharing your script?  I *think* I know the issue, but I'm a bit surprised that an Athena++ dataset would have this issue.  But I think it's that the metallicity field definition I gave you assumed that all the fluid elements were grid-based and it looks like they might be particle-based?  But it'll be easier to debug if I can just see the script you're using.  I can try to reproduce it with one of the sample Athena++ datasets here: https://yt-project.org/data/ .

Cameron

To view this discussion on the web visit https://groups.google.com/d/msgid/trident-project-users/952a608a-488b-4050-83a9-f16f692ecf06n%40googlegroups.com.

[Ryan Tanner]

Here is my script:

----START----

import yt

import trident

units_override = {"length_unit":(1.0,"pc"),

                  "time_unit":(1.0,"Myr"),

                  "mass_unit":(0.0279311,"Msun")}

fn = '../../data/SB0_AGN1/Starburst.out2.00010.athdf'

def _metallicity(field, data):

    factor = 0.3

    return (

        data.ds.arr(data["index", "ones"].d, "Zsun") * factor

    )

yt.add_field(

    ("gas", "metallicity"),

    function=_metallicity,

    sampling_type="local",

    units="Zsun",

)

ds = yt.load(fn, units_override=units_override, unit_system='cgs')

ray_start = ds.domain_left_edge

ray_end = ds.domain_right_edge

line_list = ['H', 'C', 'N', 'O', 'Mg']

ray = trident.make_simple_ray(ds,

                              start_position=ray_start,

                              end_position=ray_end,

                              data_filename="ray.h5",

                              lines=line_list,

                              ftype='gas')

p = yt.ProjectionPlot(ds, 'x', 'density')

p.annotate_ray(ray, arrow=True)

p.save('projection.png')

sg = trident.SpectrumGenerator('COS-G130M')

sg.make_spectrum(ray, lines=line_list)

sg.save_spectrum('spec_raw.txt')

sg.plot_spectrum('spec_raw.png')

sg.add_qso_spectrum()

sg.add_milky_way_foreground()

sg.apply_lsf()

sg.add_gaussian_noise(30)

sg.save_spectrum('spec_final.txt')

sg.plot_spectrum('spec_final.png')

-----END-----

If you want here is a file from my simulation:

https://drive.google.com/file/d/1rud36qINYtdk7NzHI8KimWUVUgbB84wg/view?usp=sharing

It's 2.5 GB.

[Cameron Hummels]

Yeah, my hunch was correct.  I guess I had incorrectly assumed that all Athena++ datasets were grid-based, but it seems like the fluid elements in this simulation are being recognized as particles according to the yt infrastructure.  If you swap out the code I gave you before for the dummy metallicity field with the code below, it works.  For reference, all we're doing is creating a metallicity field by using the same type as the gas density field and setting it to be 0.3 Zsun everywhere.  

OK, so hopefully this addresses all of your trident troubles! But let us know if you have more.  Also, the slack workspace is also an option for shorter-latency contact if you're interested, but the email list works too.

Cameron

import numpy as np

def _metallicity(field, data):
    factor = 0.3
    return (

        data.ds.arr(np.ones_like(data["gas", "density"]), 'Zsun') * factor

    )

yt.add_field(
    ("gas", "metallicity"),
    function=_metallicity,
    sampling_type="local",
    units="Zsun",
)

To view this discussion on the web visit https://groups.google.com/d/msgid/trident-project-users/969fcd1e-39cb-47d5-911e-6e51875f4e2fn%40googlegroups.com.

[Ryan Tanner]

Thanks! That worked great!

Now I can work on doing some science. 

But, um, one more thing. The link to join the trident slack channel doesn't work. On your trident website the link to slack just gives:

This link is no longer active

To join this workspace, you’ll need to ask the person who originally invited you for a new link.

I swear this is the last bug I'll find.

RT

[Cameron Hummels]

Great! I'm glad that everything (but slack) worked!

Here is a new slack link: https://join.slack.com/t/trident-project/shared_invite/zt-42h0uuwy-fBggZbeymnq2cB9ivtWloA .  I'll update our docs.

Cameron

To view this discussion on the web visit https://groups.google.com/d/msgid/trident-project-users/317de229-a533-418f-8b30-0bda5abc7e01n%40googlegroups.com.