yt-project / yt

Main yt repository
http://yt-project.org
Other
461 stars 276 forks source link

changing units manually in FLASH datasets #907

Closed yt-fido closed 7 years ago

yt-fido commented 9 years ago

Originally reported by: Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown)


Hello everyone, I am a new yt user; working with FLASH datafiles. FLASH files work smoothly in general, but the front-end seems to lack the possibility to change the unit system in reading, on the contrary of some other dataset (e.g. Athena) where one can specify units at load time.

I need to specify a custom unit system; either at load time or later. I mean, I want yt to read the exact numbers on the file, but attach different physical units to them.

My files may be a bit special under this aspect, but it's true that the FLASH front-end assumes cgs only. Under http://yt.readthedocs.org/en/latest/examining/supported_frontends_data.html, the user is warned that this will happen, but there's no solution or workaround proposed...

I'd please like to know whether there's a way to change the units of a dataset, or to include the possibility of doing this manually in the front-end at load time.

Thank you!


yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


Resolved by PR # 1236

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


Resolved by PR #1236

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


I see the problem here--you need to specify units in your field definition:

#!python

def _density_squared(field, data):
   return data["dens"]**2

ds.add_field("density_squared", function=_density_squared, units="g**2/cm**6")

Or something like that--not sure which is appropriate for your case, since you're changing the units. You could also do something like this, which would be unambiguous:

#!python

def _density_squared(field, data):
   return (data["dens"]**2).in_cgs()

ds.add_field("density_squared", function=_density_squared, units="g**2/cm**6")

Where we added an in_cgs to make sure the densities get into cgs first. If you need to do something else and need help, let me know.

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


Ok, if you don't get the error, it's very likely that I have done something wrong. The scripts follow almost perfectly the photon_simulator tutorial (I had to add the yt namespace in a few spots). You may get the file here:

https://www.dropbox.com/s/yti4xy98naij0hb/s200m5P42b640l10_hdf5_plt_cnt_0000?dl=0

#!python
import yt

# My data
uo = {"length_unit":(1.0,"Mpc"), "time_unit":(13770.0,"Myr"), "mass_unit":(5.23e12,"Msun")}
ds = yt.load("s200m5P42b640l10_hdf5_plt_cnt_0000", units_override=uo)

# The Athena data
#ds = yt.load("MHDSloshing/virgo_low_res.0054.vtk")

def _density_squared(field, data):
   return data["dens"]**2

ds.add_field("density_squared", function=_density_squared)

slc = yt.SlicePlot(ds, "x", ["density_squared"])
# The error shows up here
slc.set_cmap("density_squared", "gray_r")
slc.save()

I couldn't change the units from cgs in athena file: yt says that the old way is deprecated, the units override maybe wanted a different syntax... I don't know. But the slice plot should work anyway.... I also tried to add the new datafield to yt before loading the file; same error. Maybe I should update my yt version, too... is it right if I clone again the one you suggested earlier?

Thanks again!

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


Hey @sacielo, I can't reproduce the error with the GasSloshing dataset. Can you provide scripts that are failing, both for your dataset and for that one?

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


Thank you!

I'll have a look at the files in the PR, too.

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


I'll try to have a look at your points tonight or tomorrow.

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


Apparently I can't add datafields even to the GasSloshing dataset now... I get the same error (unit-related). So, maybe it is due to the units override... ?

I could have a look at it my self, but I don't know where to look for.

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


Dear John,

I have seen the PR for the various data front-end, and I think it's really a useful one!

I was testing the changes for the FLASH case, but I couldn't do much... I could print on screen and plot the data fileds on the original file and some derived one.

The units and numbers are correct, but nothing else seems to work... I think it is not something related to the units override (if I don't use it I have exactly the same errors), but I can't say much.

Most derived data fields (everyone requiring some spacial derivative) give errors like:

File "selection_routines.pyx", line 729, in yt.geometry.selection_routines.RegionSelector.__init__ (yt/geometry/selection_routines.c:20594)
RuntimeError: Error: bad Region in non-periodic domain along dimension 0. Region left edge = -0.41 code_length, Region right edge = -0.15 code_lengthDataset left edge = -0.32 code_length, Dataset right edge = 0.32 code_length

which I believe is just due to the FLASH mesh more than to the units.

I also can't add any derived field myself. I tried to experiment with the tutorial of the photon_simulator module (the one I'm mostly interested in using)

http://yt-project.org/doc/analyzing/analysis_modules/photon_simulator.html (I know you know it)

and I couldn't even add the density_squared data field. I keep getting an error like:


 File "/hydra/u/sacielo/software/yt/yt-3.x/yt/units/yt_array.py", line 402, in _unit_repr_check_same
    self.units, self.units.dimensions, units, units.dimensions)
yt.utilities.exceptions.YTUnitConversionError: Unit dimensionalities do not match. Tried to convert between code_mass**2/code_length**6 (dim (mass)**2/(length)**6) and dimensionless (dim 1).

regardless whether I use the units override or not.

I know I didn't simply type something wrong, because if then I try

g = ds.index.grids[1]
g["density_squared"]

it prints the correct value on screen (but I can't plot it in any way).

I went on with the photon_simulator tutorial, but got a nearly-blank image (but can be due to many factors).

So, I couldn't really test much, sorry...

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


Dear John,

thank you very much! I have been away for a couple of days, but now I tested your changes. By my first tests, all seems to work flawlessly. I have successfully read some of my data files (e.g. https://www.dropbox.com/s/1rv5d8vm7aurhla/sample_FLASH_custom_units_hdf5_plt_cnt_0000?dl=0) with the correct unit system, and yt could print the right value of the density, both in FLY units and in cgs (which is still the default output system).

I may have to make some other test (e.g. with derived data fields), but I have no reason to believe that there will be problems. I do not make use of magnetic fields, so I can't test any MHD quantity, but again I don't see why it should not work.

I will keep you updated on my tests, but I'm really happy of it so far!

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


Salvatore, can you test out some of the changes I made? This should be working for FLASH datasets already, at least in a bare-bones way. You can use the following steps:

#!bash

hg clone http://bitbucket.org/jzuhone/yt-3.x
cd yt-3.x
python setup.py develop

Then, using the GasSloshing dataset from the data page, try out this script:

#!python

import yt
units_override={"length_unit":(1.0,"kpc"), "time_unit":(1.0,"Myr"), "mass_unit":(1.0,"Msun")}
ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0100", units_override=units_override)
sp = ds.sphere("c", (0.1,"unitary"))
print sp["density"]

It won't make any sense, because this dataset is in CGS, but then you can try applying it to yours. We'll have to iterate on this a few times but this is the basic form it will take.

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


I think that probably the best way to resolve this issue is to take the machinery present in the Athena frontend and bring it up into the Dataset class. It will probably take me a few weeks to set this up (currently swamped with a lot of things), but it's definitely something I think we should do.

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


FLASH4.0

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


Got it! One other quick question before I formulate some thoughts on this--FLASH2 or FLASH3/4?

yt-fido commented 9 years ago

Original comment by Salvatore Cielo (Bitbucket: sacielo, GitHub: Unknown):


Hi John.

Thanks for your quick reply; sorry if this post is a bit long.

About the version of yt I'm using, I think I got a bit confused, my bad. I am actually using the 3.1-dev; I can downgrade to the stable release if you advise me to. I also see the difference with version 2.x.

Your analysis of the problem is correct: usually the calculations in FLASH are in cgs; but I have changed that in my setup. I wanted to use a unit system that is only partially implemented in FLASH, the FLY system (astrophysical friendly). In order to do that, I manually changed the values of the physical constants and other quantities (e.g. the radiative cooling, provided in the supplementary units). I have to say that it was a real struggle, especially it took long to test. But I checked that it actually prevents some numerical overflow in my setup; with these units one may work in single precision with astronomical numbers with little worry. I never advertised this, because I myself think it's a "dirty" solution; FLASH actually believes to work in cgs, but it does not. I think my troubles come from here. So, FLASH thinks to output numbers in centimeters, but in reality they are... Megaparsec!

I am actually very impressed by the way in which yt deals with units; it's... just the way it should be! I think that once I specify the correct units, yt would be able to use the correct values of the physical constants for the analysis (am I right?). And I plan, in my analysis, to use some astrophysical package (another very appealing feature of yt) for which units matter; so I have to have them right!

This is why I said that my files are "special", but I think that changing units is a freedom one wants to have... isn't it?

Thanks again for reading through this very long post! :)

Best

yt-fido commented 9 years ago

Original comment by John ZuHone (Bitbucket: jzuhone, GitHub: jzuhone):


Hi Salvatore,

First, welcome to yt!

Second, a quick question--which version of yt are you using? The docs link you pointed to is from version 2.x, which is not the most recent version. Support for units is much better in yt 3.0: http://yt-project.org/docs/3.0/analyzing/units/index.html#units

The updated version of the docs for FLASH data is here: http://yt-project.org/docs/3.0/examining/loading_data.html#flash-data

Now, to your actual question. The reason why the FLASH frontend assumes CGS is that almost all FLASH calculations are technically done in CGS. It is true that you can set up length, time, etc. to be whatever you want, but the internal constants (G, kB, etc.) will all be in CGS. The only exception I am aware of for this is MKS, but I don't know if anyone is using it (maybe you are?), and whether or not it even works correctly.

However, there is definitely room for more flexibility on this on the yt side and probably also the FLASH side (for which I also develop from time to time). Could you give more details on your exact use case, in particular how the units are set up in your simulation? We should be able to work something out.