Sorry, the previous code is a kind of mixed up… @derived_field(name = "CylCellVolume") def CylCellVolume(field, data): return data["dx"]*data["dy"]*2*math.pi*data["x"] @derived_field(name = "CylCellMass") def CylCellMass(field, data): return data["CylCellVolume"]*data["Density"] Best wishes. Suoqing On Jun 18, 2013, at 1:00 PM, Suoqing Ji wrote:

Hi Roman,

I assume the quantity 'CylindricalCellMass' is a derived field defined by yourself? I usually define the 2d cylindrical cell volume as followings:

@derived_field(name = "CylCellVolume") # CellVolume in yt is wrong def realVolume(field, data): return data["dx"]*data["dy"]*2*math.pi*data["x"]

@derived_field(name = "CylCellMass") def realCellMass(field, data): return data["CylCellVolume"]*data["Density"]

And for my data set, the output given by YT (use dd.quantities['TotalQuantity']('CylCellMass')[0]) and FLASH subroutine (IO_writeIntegralQuantities) match down to the machine accuracy.

Best wishes, Suoqing

On Jun 18, 2013, at 12:46 PM, Roman Yurchak wrote:

...

Hello,

I'm interested in getting some integrated quantities from a cylindrical 2d FLASH simulation.

Just to be sure, I started by comparing the total simulation mass computed by yt ( with pf.h.all_data().quantities['TotalQuantity'] etc.) and the one provided in FLASH logs, and the results are rather different [see attachment]: * _Values_: difference of up 8% between the two methods * _Time evolution_: total mass is decreasing in time according to FLASH (which is ok with open boundary conditions) but it is increasing with the yt's integration. * _Sensitivity to resolution_: there are jumps in total density given by yt for the time steps where the maximum refinement level is manually decreased.

Would anyone know what might be the reason for this behaviour?

Just in case, I tried to do a second order integration instead of simply summing the values in all the cells, but if didn't make an error, the correction was negligible.

I can provide data and the yt script I'm using by PM if some could look into this.

Thank you, Regards, -- Roman Yurchak <integrated_quantities_FLASH_cyl.txt>_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

Hi Roman, On Tue, Jun 18, 2013 at 6:25 PM, Roman Yurchak <rth@crans.org> wrote:

Hi Suoqing,

On 18/06/13 19:00, Suoqing Ji wrote:

...

I assume the quantity 'CylindricalCellMass' is a derived field defined by yourself? I usually define the 2d cylindrical cell volume as followings: [...] yes, that's pretty much the expression I was using.

...

And for my data set, the output given by YT (use dd.quantities['TotalQuantity']('CylCellMass')[0]) and FLASH subroutine (IO_writeIntegralQuantities) match down to the machine accuracy.

ahh, thank you for your response! It made me realize two things:

(a) I made a stupid mistake while sorting different time steps, and for a normal 2d cylindrical FLASH simulation I do have convergence to machine precision between the two methods.

Ah! Good, that is great to hear.

(b) I was post-processing a simulation that has a few coarse blocks removed and used as new boundary conditions in order to get a non-rectangular domain (see source/Simulation/SimulationMain/WindTunnel/Simulation_defineDomain.F90 in FLASH source, or the illustration here: http://perso.crans.org/yurchak/i/illustations_yt.pdf ). This seems to be the origin of the issue I've been experiencing (and that's still there).

I think I see what's happening here. There are a few diagnostics to run to see if we can understand and then correct yt's understanding of the mesh; basically, yt will construct a full mesh based on what it sees in the file, and we need to make sure that the mesh matches what is being evolved in the simulation. For instance, it's not clear to me completely whether or not the block that has been removed exists in the mesh described in the FLASH output file, and if so, if we can simply remove it. We may be able to do this directly via the FLASH frontend. I do not think so, but there may be assumptions in the FLASH frontend that the system resembles the typical refinement pattern for FLASH, which this may not match. In the past we've been able to do funny things to the hierarchy via the Stream frontend -- interposing IO and whatnot, mocking up new hierarchies -- which I would be happy to make work with this output if we can't make the FLASH frontend work exactly as desired.

I looks like yt, treats this configuration partly ok as for instance none of leaf blocks has for parent a block that was removed. On the other side though: * there this issue with d.quantities['TotalQuantity'] * a slice shows that there is some plasma propagating through the region that's supposed to be treated as outside of the simulation (with removed blocks) and there is nothing there in VisIt.

So, anyway, thanks for narrowing it down.

Would you have any suggestions how to deal with this issue?

So let's start out by taking a look at a few things. These are some common diagnostics to figure out whether or not a hierarchy is "well-posed" in yt. One thing I am wondering is if the non-regular cutout is causing irregular problems -- specifically, I am thinking that it may not be "masked" in the parent grids with which it overlaps. So those operations that utilize "masking" information will not know about the missing block, but operations that *don't* will not. I wrote up some simple diagnostics that I *think* will be sufficient to detect problems like this. They work properly for the yt sample FLASH dataset, and I think they should draw out problems like the ones you describe. Here is a script: http://paste.yt-project.org/show/3619/ Can you run this on your dataset and see if it works? As a note, there have been issues recently with hierarchies getting serialized; I'd recommend removing and .yt files scattered about and also setting this in ~/.yt/config : [yt] serialize = False Just to be on the safe side. -Matt

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