Thanks Britton, it works perfectly. Cheers, JC On 5/1/12 1:04 PM, Britton Smith wrote:

Hi JC,

This is happening because the generation of the field happens grid by grid, so you're comparing the field for the entire region to the field for just one grid.

A very similar question was posed on yt-dev about a week ago, but for some reason I can't find a record of it in the archive. I've found it in my email so I'll paste it here. Sorry for it being in reverse order and a bit difficult to read.

Britton

---------- Forwarded message ---------- From: *Geoffrey So* mailto:gsiisg@gmail.com> Date: Thu, Apr 26, 2012 at 4:11 PM Subject: Re: [yt-dev] derived field from two datasets To: yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org

ok changing the line to data2 = pf2.h.grids[data.id http://data.id/ - data._id_offset] worked

Thanks so much:-)

From G.S.

On Thu, Apr 26, 2012 at 1:04 PM, Matthew Turk mailto:matthewturk@gmail.com> wrote:

Ah, don't ask for [field], just get .grids[ ... ].

-Matt

On Thu, Apr 26, 2012 at 4:02 PM, Geoffrey So mailto:gsiisg@gmail.com> wrote: > The new add field line > pf1.field_info.add_field("CompareIonization", function = _NewField, > validators = [ValidateGridType()]) > but that got me a new set of errors > KeyError: 0x1065cbb90> > > full trace > Traceback pasted to http://paste.yt-project.org/show/sIW43mg7AdRaKJELo0U1 > > From > G.S. > > On Thu, Apr 26, 2012 at 12:54 PM, Matthew Turk mailto:matthewturk@gmail.com> > wrote: >> >> Hi Geoffrey, >> >> AH! Rats. I forgot one thing. In the call to add_field, include a >> validator to ensure you're getting grids, not regions: >> >> validators = [ValidateGridType()] >> >> This script, by the way, will *only* work with unigrid where the grid >> IDs match between the two datasets. (Doing it with AMR is possible >> but will be a hassle.) What that line does is cross-correlate grids >> between the two parameter files. >> >> -Matt >> >> On Thu, Apr 26, 2012 at 3:50 PM, Geoffrey So mailto:gsiisg@gmail.com> wrote: >> > Yeah I think that should do exactly what I want. >> > >> > The script is at >> > http://paste.yt-project.org/show/2324/ >> > >> > gave me these error: >> > http://paste.yt-project.org/show/5XFq1xU2iVxiA9MTtT3v >> > >> > pf_compare was ran without error, the error occurs when I tried to call >> > on dd["CompareIonization"] >> > >> > Can you also explain what this line does with the data.id http://data.id/ and >> > _id_offset? >> > data2 = pf2.h.grids[data.id http://data.id/ - data._id_offset][field] >> > >> > From >> > G.S. >> > >> > On Thu, Apr 26, 2012 at 12:15 PM, Matthew Turk mailto:matthewturk@gmail.com> >> > wrote: >> >> >> >> Hi Geoffrey, >> >> >> >> On Thu, Apr 26, 2012 at 3:09 PM, Geoffrey So mailto:gsiisg@gmail.com> wrote: >> >> > This is very close to what I was trying to do, except I wanted >> >> > something >> >> > slightly more general where I can compare two distinct fields from >> >> > two >> >> > datasets. (what I am trying to do is to use >> >> > HI_Density/(HI_Density+HII_Density)) where the HI and HII are stored >> >> > separately, so I'll need pf1 with field1 and pf2 with field2.) >> >> > >> >> > Rick recently showed me a rendering where people at Argonne made a >> >> > rendering >> >> > of the fractional difference in HI or HII from two simulations, one >> >> > with >> >> > radiation transport one without, to compare the effects. The ability >> >> > to >> >> > call two different fields from two different simulation parameters in >> >> > creating a derived field would be essential to doing something like >> >> > that. >> >> >> >> Here's a video of it, for everyone curious: >> >> >> >> http://www.youtube.com/watch?v=hDh4NB31e2s >> >> >> >> > >> >> > Forgive my yt skills. I tried to first use the script you provided >> >> > as >> >> > an >> >> > example, but I wasn't sure how to get values from pf_compare, I tried >> >> > >> >> > dd1=pf1.h.all_data() >> >> > dd1['NewField'] >> >> > >> >> > but it gives: >> >> > >> >> > TypeError: _NewField() takes exactly 4 arguments (2 given) >> >> >> >> Okay, I should have been more clear. What you want to do is call: >> >> >> >> pf_compare(pf1, pf2, some_field) >> >> >> >> And then pf1 should have the field accessible. It doesn't return the >> >> field; what it does is make the field accessible to yt functionality. >> >> >> >> Here's a closer example, based on what you have asked for: >> >> >> >> def pf_compare(pf1, pf2): >> >> def _NewField(field, data): >> >> data2 = pf2.h.grids[data.id http://data.id/ - data._id_offset][field] >> >> return >> >> data["HI_Density"]/(data2["HI_Density"]+data2["HII_Density"])) >> >> pf1.field_info.add_field("CompareIonization", function = _NewField) >> >> >> >> pf1 = load( ... ) >> >> pf2 = load( ... ) >> >> pf1.h >> >> pf_compare(pf1, pf2) >> >> >> >> # Now you can do: >> >> >> >> dd = pf1.h.all_data() >> >> dd["CompareIonization"] >> >> >> >> >> >> -Matt >> >> >> >> > >> >> > >> >> > >> >> > -------------------------------------------------------------------------------------- >> >> > >> >> > What I was hoping to be able to do is something like: >> >> > >> >> > def pf_compare(pf1, field1, pf2, field2): >> >> > def _NewField(field1, data1, field2, data2): >> >> > return data1[field1] / data2[field2] # or any calculation >> >> > involving data1,2 and field1,2 >> >> > pf1.field_info.add_field("NewField", function = _NewField) >> >> > >> >> > pf1 = load('DD###1') >> >> > pf2 = load('DD###2') >> >> > >> >> > and call >> >> > >> >> > pf1.h.all_data()['NewField'] >> >> > >> >> > Where "NewField" involves two different fields from two different >> >> > dataset >> >> > >> >> > This seems very possible from the example you wrote Matt, but I guess >> >> > I'm >> >> > stuck on how to call the values or make a plot of "NewField" once it >> >> > is >> >> > defined. >> >> > >> >> > From >> >> > G.S. >> >> > >> >> > >> >> > On Thu, Apr 26, 2012 at 5:02 AM, Matthew Turk mailto:matthewturk@gmail.com> >> >> > wrote: >> >> >> >> >> >> Hi Geoffrey, >> >> >> >> >> >> There's no example per se, but I'll try to give you one. This >> >> >> sketch >> >> >> should give you an idea. There's a way to do this for AMR grids >> >> >> which >> >> >> would be similar (you'd use a covering grid) but I think you're >> >> >> using >> >> >> static refinement so that makes it a bit easier. We want to wrap >> >> >> this >> >> >> in a function: >> >> >> >> >> >> def pf_compare(pf1, pf2, field): >> >> >> def _NewField(field, data): >> >> >> return data[field] / pf2.h.grids[data.id http://data.id/ - >> >> >> data._id_offset][field] >> >> >> pf1.field_info.add_field("Compare%s" % field, function = >> >> >> _NewField) >> >> >> >> >> >> pf1 = load( ... ) >> >> >> pf2 = load( ... ) >> >> >> pf1.h >> >> >> pf_compare(pf1, pf2, "Density") >> >> >> >> >> >> Now you should have CompareDensity. >> >> >> >> >> >> Let me know if that works. In the discussion document for yt 3.0 >> >> >> Britton brought up that we should have a grammar for this type of >> >> >> situation. >> >> >> >> >> >> -Matt >> >> >> >> >> >> On Wed, Apr 25, 2012 at 3:39 PM, Geoffrey So mailto:gsiisg@gmail.com> >> >> >> wrote: >> >> >> > Hi, I was wondering if there's an example of how to create a >> >> >> > derived >> >> >> > field >> >> >> > from two different datasets. >> >> >> > >> >> >> > From >> >> >> > G.S. >> >> >> > >> >> >> > _______________________________________________ >> >> >> > yt-dev mailing list >> >> >> > yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org >> >> >> > http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org >> >> >> > >> >> >> _______________________________________________ >> >> >> yt-dev mailing list >> >> >> yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org >> >> >> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org >> >> > >> >> > >> >> > >> >> > _______________________________________________ >> >> > yt-dev mailing list >> >> > yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org >> >> > http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org >> >> > >> >> _______________________________________________ >> >> yt-dev mailing list >> >> yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org >> >> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org >> > >> > >> > >> > _______________________________________________ >> > yt-dev mailing list >> > yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org >> > http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org >> > >> _______________________________________________ >> yt-dev mailing list >> yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org >> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org > > > > _______________________________________________ > yt-dev mailing list > yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org > http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org > _______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org

_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org mailto:yt-dev@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org

On Tue, May 1, 2012 at 12:27 PM, Jean-Claude Passy mailto:jcpassy@gmail.com> wrote:

Hi all,

I am trying to compare the energy between two datasets. I execute the following script that uses the derived field 'RelativeEnergy':

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

def _RelativeEnergy(field, data):

return (data["TotalEnergy"]-data.pf.parameters["TotalEnergy_ref"])/data.pf.parameters["TotalEnergy_ref"]

add_field("RelativeEnergy", function=_RelativeEnergy)

filen1 = 'helix/DD0000/CE0000'

pf1 = load(filen1)

region1 = pf1.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0])

energy1 = region1["TotalEnergy"]

filen2 = 'parallel/DD0000/CE0000'

pf2 = load(filen2)

region2 = pf2.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0])

energy2 = region2["TotalEnergy"]

# Adding reference energy into pf2

pf2.parameters["TotalEnergy_ref"] = energy1

delta_energy = region2["RelativeEnergy"]

pc = PlotCollection(pf2)

pc.add_slice("RelativeEnergy",2,center=[0.5, 0.5, 0.5])

pc.set_lim((0.4, 0.6, 0.4, 0.6))

pc.save('DeltaEnergy')

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

I get the error pasted here: http://paste.yt-project.org/show/DPgMjV44JktpXExCljUj/ The variable delta_energy looks good (same type and size and energy1), but the problem arises when I try to take a slice of it.

Any help much appreciated. Thanks a lot,

JC

_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org mailto:yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

Hi JC, Are the hierarchies identical between pf1 and pf2? If not, I think that could be a problem. It could be that with the slice you passed only through root grids where the hierarchy is identical. Doing a projection that acts upon the entire domain would expose this issue. If this is the case, you may need to create a constant resolution data container for use in the comparison. Britton On Tue, May 1, 2012 at 3:58 PM, Jean-Claude Passy wrote:

Oops, I spoke too fast. I can make a slice but if I try a projection, the same kind of error occurs: http://paste.yt-project.org/show/7QiwjVoYXKYxCn3Iu2u9/

This is the script I am using:

------------------------------------------------------------------------- def pf_compare(pf1, pf2): def _NewField(field, data): data2 = pf2.h.grids[data.id - data._id_offset] return abs(data2["TotalEnergy"]-data["TotalEnergy"])/data["TotalEnergy"] pf1.field_info.add_field("CompareEnergy", function = _NewField, validators = [ValidateGridType()]) pf1 = load('helix/DD0000/CE0000') pf2 = load('parallel/DD0000/CE0000') # Comparison pf_compare(pf1, pf2) region = pf1.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) compare_energy = region["CompareEnergy"] pc = PlotCollection(pf1) pc.add_slice("CompareEnergy",2,center=[0.5, 0.5, 0.5]) pc.add_projection("CompareEnergy",2) pc.set_lim((0.4, 0.6, 0.4, 0.6)) pc.save() ------------------------------------------------------------------------- Any idea? Cheers,

JC

On 5/1/12 3:35 PM, Jean-Claude Passy wrote:

Thanks Britton, it works perfectly.

Cheers,

JC

_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org