[yt-users] Problem with making SlicePlot and ProjectionPlot from new derived field
Bili Dong - Gmail
qobilidop at gmail.com
Thu Oct 5 16:40:44 PDT 2017
Hi Tazkera,
Now I understand where the error comes from. It's not from the unit
system but from the difference between 'Gas' and 'gas' fields. Let me show
you the code that works and explain later. The corrected code looks like
this:
```
dd = ds.all_data()
# derive a new domain width with conversion factor
(6Mpc*current_time/hubble_constant)
def _NewDomainWidth(field, data):
return ds.domain_width*((3.086e18*6e6*0.909090893796)/(0.7))
ds.add_field(('gas', 'PhysicalDomainWidth'), function=_NewDomainWidth,
units="code_length")
# derive a new domain center with conversion factor
(6Mpc*current_time/hubble_constant)
def _NewDomainCenter(field, data):
return ds.domain_center*((3.086e18*6e6*0.909090893796)/(0.7))
ds.add_field(('gas', 'PhysicalDomainCenter'), function=_NewDomainCenter,
units="code_length")
# derive a new electron pressure field
def _ElectronPressure(field, data):
return data['Gas','Density']*(5.49e-7)*data['Gas','Temperature']
ds.add_field(('Gas', 'ElectronPressure'), function=_ElectronPressure,
units="(code_mass*K)/(code_length**3)",
particle_type=True)
f_smooth = ds.add_smoothed_particle_field(('Gas', 'ElectronPressure'))
# f_smooth equals ('deposit', 'Gas_smoothed_ElectronPressure')
# integration of pressure field along line of sight
yt.ProjectionPlot(ds, 'z', f_smooth, width = dd['gas',
'PhysicalDomainWidth'],
center = dd['gas', 'PhysicalDomainCenter'] )
```
First, let's go back to the error message. Where do the magic numbers (39678,)
(54608,) come from? Print out dd['Gas', 'Density'].shape and dd['gas',
'density'].shape and you'll see. The issue with your original code is that
you tried to create a 'gas' field from a 'Gas' field but their dimension
doesn't match. Here 'Gas' field is a particle field retrieved from your
snapshot directly, but 'gas' field is an SPH field calculated from the
particle field. They have different dimensions. More detailed documentation
could be found here:
http://yt-project.org/docs/dev/analyzing/fields.html#sph-fields
That's for solving the error. As a side note, I do not recommend defining
fields like your dd['gas', 'PhysicalDomainCenter'], though it works (to my
surprise). Because usually we assume the same field type means the same
dimension (as in the first number of shape). Your case is more suitable to
define a new variable.
Hope it helps.
Best,
Bili
On Thu, Oct 5, 2017 at 3:54 PM, Bili Dong - Gmail <qobilidop at gmail.com>
wrote:
>
> On Thu, Oct 5, 2017 at 3:48 PM, Bili Dong - Gmail <qobilidop at gmail.com>
> wrote:
>
>> Hi Tazkera,
>>
>> I think yt should be able to understand your code units and do the unit
>> conversion for you, or easily configured to interpret your code units. See
>> an example here
>> <http://yt-project.org/docs/dev/examining/loading_data.html#tipsy-data>.
>> This I feel is the better solution.
>>
>> But if you want to the conversion by hand, a quick fix to your code is to
>> define new variables for domain width and domain center:
>>
>> ```
>> new_domain_width = ds.domain_width*((3.086e18*6e6*0.909090893796)/(0.7))
>> new_domain_center = ds.domain_center*((3.086e18*
>> 6e6*0.909090893796)/(0.7))
>> yt.ProjectionPlot(ds, 'z', ('gas', 'ElectronPressure'), width=
>> new_domain_width[0],
>> center=new_domain_center)
>> ```
>>
>> The issue with your code is that ds.domain_width and ds.domain_center
>> are not fields themselves. If you try to add new fields for those constant
>> values, you'll get a long list of those values as you could read from the
>> error message. You could print out dd['gas', 'PhysicalDomainWidth'] or dd['gas',
>> 'PhysicalDomainCenter'] in your code to see what they look like. That
>> might help you understand what was wrong.
>>
>
> Sorry I was wrong. This is not true.
>
>>
>> Best,
>> Bili
>>
>>
>> On Thu, Oct 5, 2017 at 2:59 PM, tazkera haque <h.tazkera at gmail.com>
>> wrote:
>>
>>> Dear yt people,
>>>
>>> I am working on a tipsy output file from a n-body simulation. My file
>>> contains some scaled units for fields like density, coordinates, mass etc.
>>> To make a physical sense of the data, I had to convert those fields units
>>> to some physical units (c.g.s) and I derived new fields having cgs unit.
>>> for example:
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> *dd = ds.all_data()# derive a new domain width with conversion factor
>>> (6Mpc*current_time/hubble_constant) def _NewDomainWidth(field, data):
>>> return
>>> ds.domain_width*((3.086e18*6e6*0.909090893796)/(0.7))ds.add_field(('gas',
>>> 'PhysicalDomainWidth'), function=_NewDomainWidth, units="code_length")*
>>>
>>>
>>>
>>>
>>>
>>> *# derive a new domain center with conversion factor
>>> (6Mpc*current_time/hubble_constant) def _NewDomainCenter(field, data):
>>> return
>>> ds.domain_center*((3.086e18*6e6*0.909090893796)/(0.7))ds.add_field(('gas',
>>> 'PhysicalDomainCenter'), function=_NewDomainCenter, units="code_length")*
>>>
>>> note that here "code_length" is actually in "cm". but I could not just
>>> do dd["Gas",'Coordinates'].in_units('cm'), because this conversion does
>>> not take into account some constant corrections and returns the wrong
>>> length magnitude.
>>>
>>> My goal was to make sliceplot and projectionplot of a newly derived
>>> field :
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> *# derive a new electron pressure fielddef _ElectronPressure(field,
>>> data): return
>>> data['Gas','Density']*(5.49e-7)*data['Gas','Temperature']ds.add_field(('gas',
>>> 'ElectronPressure'), function=_ElectronPressure,
>>> units="(code_mass*K)/(code_length**3)")*
>>>
>>> *# integration of pressure field along line of sight*
>>>
>>> *yt.ProjectionPlot(ds, 'z', ('gas', 'ElectronPressure'), width =
>>> dd['gas', 'PhysicalDomainWidth'], center = dd['gas',
>>> 'PhysicalDomainCenter'] )*
>>>
>>>
>>> *and it is returning me *
>>>
>>> ValueError Traceback (most recent call last)<ipython-input-8-fa351b6c09d1> in <module>() 1 yt.ProjectionPlot(ds, 'z', ('gas', 'ElectronPressure'), width = dd['gas', 'PhysicalDomainWidth'], ----> 2 center = dd['gas', 'PhysicalDomainCenter'] )
>>>
>>> ValueError: operands could not be broadcast together with shapes (39678,) (54608,)
>>>
>>> Now I cannot use width =ds.domain_width and center = ds.domain_center for the same unit error. I was wondering how I can make the projectionplot for my newly derived field and newly derived domain width and domain center.
>>>
>>>
>>>
>>> Your suggestion is most welcome
>>>
>>>
>>> Best
>>>
>>> Tazkera
>>>
>>>
>>> _______________________________________________
>>> yt-users mailing list
>>> yt-users at lists.spacepope.org
>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20171005/e9527c4f/attachment-0001.html>
More information about the yt-users
mailing list