Hi Matt, I will try to explain what my goal is. I have a 1D profile of a star that is physically at the equilibrium. Regarding the problem initialization, I set up the density and the total specific energy fields according to the model, velocities equal 0. However, the star is not numerically stable so I have to relax it. Therefore, I divide the velocity by 2 at each timestep, and I let Enzo run for a few dynamical times. Then, I read the last data dump and create an relaxed_model.dat: ########################################################## region = pf.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) output = 'relaxed_model.dat' f = open(output,'w') print >>f, '# Grid Index rho pressure' density = region["Density"] pressure = region["Pressure"] size = NP.size(density) for i in range(0,size,1): dtmp = density[i] ptmp = pressure[i] print >>f, "0", i, '%06e' % dtmp, '%06e' % ptmp ########################################################## so I end up with a file like this: ########################################################## # Restart file for only one grid # Grid Index rho pressure 0 0 5.099011e-10 2.568020e+05 0 1 5.099011e-10 2.568020e+05 0 2 5.099011e-10 2.568020e+05 ... ########################################################## Finally, I read index, density, pressure from that previous file and do: ########################################################## for (k = 0; k < GridDimension[2]; k++) for (j = 0; j < GridDimension[1]; j++) for (i = 0; i < GridDimension[0]; i++) { index = i + GridDimension[0]*(j + GridDimension[1]*k); density /= DensityUnits; pressure /= PressureUnits; BaryonField[0][index] = density; BaryonField[1][index] = pressure / ((Gamma - 1.0) * density); ########################################################## The problem is GridDimensions contain the ghost zones so the variable index does not match with the actual index read in the file relaxed_model.dat. That is why I wanted to have that file values for the ghost zones as well. Does it make sense ? Do you have any suggestion ? Thanks for your help, JC Matthew Turk a écrit :

Hi Jean-Claude,

There are a couple aspects to this. The first is that Enzo doesn't output the ghost zones -- so any ghost zones handled inside yt are generated by yt. Were Enzo to output the ghost zones, we would probably be able to handle this, but it does not.

Derived fields can depend on the generation of ghost zones, but keep in mind that these are ghost zones generated by yt. These ghost zones are constructed in a similar manner to how Enzo generates them, but there may be minor differences. You can manually inspect ghost zones on a *grid* by calling retrieve_ghost_zones on that grid.

If you could tell us a bit more about your goal, maybe we could help out a bit more?

-Matt

On Wed, Jun 2, 2010 at 9:32 AM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

...

Dear all,

in order to set up my simulation, I need to access the values of certain fields in the ghost zones. If I do something like:

pf = load(data) region = pf.h.region(...) x = region[field]

x contains the values of field for the physical grid only. Is there a way I can get the same 1D-array but with the ghost zones included as well ?

Thanks for your help,

Jean-Claude

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

Or if you just want to read the data directly back in, why not do a restart ("enzo -r")? Greg On Jun 2, 2010, at 3:36 PM, David Collins wrote:

Hi, JC--

Would it be easier to skip YT all together, and just have your problem generator read the HDF5 file produced by the 1d run directly? This would streamline your processes, and you wouldn't have the noise incurred by both the ASCII file format and the two units conversions.

I would probably approach your task by looping over only the active zones in enzo, and incrementing an extra counter.

int counter = 0; for(int field=0;field<NumberOfBaryonFields;field++) for(int k=GridStartIndex[2];k<=GridEndIndex[2]; k++) for(int j=GridStartIndex[1];j<=GridEndIndex[1]; j++) for(int i=GridStartIndex[0];i<=GridEndIndex[0]; i++){ index = i + GridDimension[0]*(j + GridDimension[1]*k); BaryonField[ DensNum ][ i ] = DensityFromFile[ counter ]; counter++ }

This code snippit came from the wiki page on Baryon access, for future reference: http://lca.ucsd.edu/projects/enzo/wiki/Tutorials/BaryonFieldAccess

Does that help? d.

On Wed, Jun 2, 2010 at 12:07 PM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

...

Hi Matt,

I will try to explain what my goal is. I have a 1D profile of a star that is physically at the equilibrium. Regarding the problem initialization, I set up the density and the total specific energy fields according to the model, velocities equal 0. However, the star is not numerically stable so I have to relax it. Therefore, I divide the velocity by 2 at each timestep, and I let Enzo run for a few dynamical times. Then, I read the last data dump and create an relaxed_model.dat:

########################################################## region = pf.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) output = 'relaxed_model.dat' f = open(output,'w') print >>f, '# Grid Index rho pressure' density = region["Density"] pressure = region["Pressure"] size = NP.size(density) for i in range(0,size,1): dtmp = density[i] ptmp = pressure[i] print >>f, "0", i, '%06e' % dtmp, '%06e' % ptmp ##########################################################

so I end up with a file like this:

########################################################## # Restart file for only one grid # Grid Index rho pressure 0 0 5.099011e-10 2.568020e+05 0 1 5.099011e-10 2.568020e+05 0 2 5.099011e-10 2.568020e+05 ... ########################################################## Finally, I read index, density, pressure from that previous file and do:

########################################################## for (k = 0; k < GridDimension[2]; k++) for (j = 0; j < GridDimension[1]; j++) for (i = 0; i < GridDimension[0]; i++) { index = i + GridDimension[0]*(j + GridDimension[1]*k); density /= DensityUnits; pressure /= PressureUnits; BaryonField[0][index] = density; BaryonField[1][index] = pressure / ((Gamma - 1.0) * density); ########################################################## The problem is GridDimensions contain the ghost zones so the variable index does not match with the actual index read in the file relaxed_model.dat. That is why I wanted to have that file values for the ghost zones as well.

Does it make sense ? Do you have any suggestion ?

Thanks for your help,

JC

Matthew Turk a écrit :

Hi Jean-Claude,

There are a couple aspects to this. The first is that Enzo doesn't output the ghost zones -- so any ghost zones handled inside yt are generated by yt. Were Enzo to output the ghost zones, we would probably be able to handle this, but it does not.

Derived fields can depend on the generation of ghost zones, but keep in mind that these are ghost zones generated by yt. These ghost zones are constructed in a similar manner to how Enzo generates them, but there may be minor differences. You can manually inspect ghost zones on a *grid* by calling retrieve_ghost_zones on that grid.

If you could tell us a bit more about your goal, maybe we could help out a bit more?

-Matt

On Wed, Jun 2, 2010 at 9:32 AM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

Dear all,

in order to set up my simulation, I need to access the values of certain fields in the ghost zones. If I do something like:

pf = load(data) region = pf.h.region(...) x = region[field]

x contains the values of field for the physical grid only. Is there a way I can get the same 1D-array but with the ghost zones included as well ?

Thanks for your help,

Jean-Claude

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

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

-- Sent from my Stone Tablet and carried by my Pterodactyl. _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

The issue is I can't trust any simulation done with restarting Enzo. Look at

I wouldn't be so sure that the solution you propose isn't going to suffer from the same problems that you see with restarting. The restarting machinery doesn't do anything fancy, there aren't many places for errors to creep in there. So baring a few things, this is potentially a deeper issue. One thing to check, which you may have done already, is to ensure that you're writing with the same precision that you're computing with. Ensure that the top 5 things in make show-config are all 64 (or 32, if you must) CONFIG_PRECISION: 64 CONFIG_PARTICLES: 64 CONFIG_INTEGERS: 64 CONFIG_INITS: 64 CONFIG_IO: 64 (Inits doesn't matter for you, but the last one, CONFIG_IO, is important) Have you checked that restarting without changing anything at all gives different answers? What compiler optimizations are you using? Try with -O1 or -O0, sometimes more aggressive compiler options can change the answer. d.

those two plots: http://drop.io/0etjfmc. I don't want to bother you with too many details but they represent the separation between my two particles - the core of the RGB star and the compact companion.

The first run was with dumping every 80 cycles - dumping based on CycleSkip. It stopped at dumping #196. I restarted it and got data #197, #198,...Then, I wanted to check that restarting did not mess things up, so I run the same sim again but with a dumping frequency twice smaller. So data #98 and #99 correspond to the same cycle as file #196 and #198 respectively. Files #98 and #196 are exactly the same (which is expected) but #99 and #198 are not. They are slightly different and lead at the end to a unphysical behavior - the companion bouncing outwards instead of sinking deeper.

I don't what is causing that, may be the new structure of particles I am using. Has anyone already had that kind of problem ? I will inspect the data at some point in the future but in the meantime, I cannot restart Enzo.

Cheers,

JC

Greg Bryan a écrit :

Or if you just want to read the data directly back in, why not do a restart ("enzo -r")?

Greg

On Jun 2, 2010, at 3:36 PM, David Collins wrote:

Hi, JC--

Would it be easier to skip YT all together, and just have your problem generator read the HDF5 file produced by the 1d run directly? This would streamline your processes, and you wouldn't have the noise incurred by both the ASCII file format and the two units conversions.

I would probably approach your task by looping over only the active zones in enzo, and incrementing an extra counter.

int counter = 0; for(int field=0;field<NumberOfBaryonFields;field++) for(int k=GridStartIndex[2];k<=GridEndIndex[2]; k++) for(int j=GridStartIndex[1];j<=GridEndIndex[1]; j++) for(int i=GridStartIndex[0];i<=GridEndIndex[0]; i++){ index = i + GridDimension[0]*(j + GridDimension[1]*k); BaryonField[ DensNum ][ i ] = DensityFromFile[ counter ]; counter++ }

This code snippit came from the wiki page on Baryon access, for future reference: http://lca.ucsd.edu/projects/enzo/wiki/Tutorials/BaryonFieldAccess

Does that help? d.

On Wed, Jun 2, 2010 at 12:07 PM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

Hi Matt,

I will try to explain what my goal is. I have a 1D profile of a star that is physically at the equilibrium. Regarding the problem initialization, I set up the density and the total specific energy fields according to the model, velocities equal 0. However, the star is not numerically stable so I have to relax it. Therefore, I divide the velocity by 2 at each timestep, and I let Enzo run for a few dynamical times. Then, I read the last data dump and create an relaxed_model.dat:

########################################################## region = pf.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) output = 'relaxed_model.dat' f = open(output,'w') print >>f, '# Grid Index rho pressure' density = region["Density"] pressure = region["Pressure"] size = NP.size(density) for i in range(0,size,1): dtmp = density[i] ptmp = pressure[i] print >>f, "0", i, '%06e' % dtmp, '%06e' % ptmp ##########################################################

so I end up with a file like this:

########################################################## # Restart file for only one grid # Grid Index rho pressure 0 0 5.099011e-10 2.568020e+05 0 1 5.099011e-10 2.568020e+05 0 2 5.099011e-10 2.568020e+05 ... ########################################################## Finally, I read index, density, pressure from that previous file and do:

########################################################## for (k = 0; k < GridDimension[2]; k++) for (j = 0; j < GridDimension[1]; j++) for (i = 0; i < GridDimension[0]; i++) { index = i + GridDimension[0]*(j + GridDimension[1]*k); density /= DensityUnits; pressure /= PressureUnits; BaryonField[0][index] = density; BaryonField[1][index] = pressure / ((Gamma - 1.0) * density); ########################################################## The problem is GridDimensions contain the ghost zones so the variable index does not match with the actual index read in the file relaxed_model.dat. That is why I wanted to have that file values for the ghost zones as well.

Does it make sense ? Do you have any suggestion ?

Thanks for your help,

JC

Matthew Turk a écrit :

Hi Jean-Claude,

There are a couple aspects to this. The first is that Enzo doesn't output the ghost zones -- so any ghost zones handled inside yt are generated by yt. Were Enzo to output the ghost zones, we would probably be able to handle this, but it does not.

Derived fields can depend on the generation of ghost zones, but keep in mind that these are ghost zones generated by yt. These ghost zones are constructed in a similar manner to how Enzo generates them, but there may be minor differences. You can manually inspect ghost zones on a *grid* by calling retrieve_ghost_zones on that grid.

If you could tell us a bit more about your goal, maybe we could help out a bit more?

-Matt

On Wed, Jun 2, 2010 at 9:32 AM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

Dear all,

in order to set up my simulation, I need to access the values of certain fields in the ghost zones. If I do something like:

pf = load(data) region = pf.h.region(...) x = region[field]

x contains the values of field for the physical grid only. Is there a way I can get the same 1D-array but with the ghost zones included as well ?

Thanks for your help,

Jean-Claude

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

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

-- Sent from my Stone Tablet and carried by my Pterodactyl. _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

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

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

-- Sent from my Stone Tablet and carried by my Pterodactyl.

your guess was totally right. I was computing with 64 bits and writing with 32 bits only. I fixed it and restarting Enzo worked well. Great !

Great, I'm glad that worked. Entire fields of science have been created because of that very issue. It's subtle, but a killer.

    - change some global parameters like StopCycle, CycleDataDump, GlobalDir,... Would it work if I only edit the files DDnnnn/CommonEnvelopennnn ? Wouldn't it create a conflict somewhere ?

Changing cycle based outputs won't change the answer. Changing timestep based outputs (like DtDataDump, Redshift) might introduce a little diffusion, because there will usually be a short timestep to make the output time right. As long as you replace GlobalDir everywhere, you're fine. ls -1 |grep -v grid | sed -i 's/GlobalDir.*/GlobalDir = $A' should do it, but double check. (that's ls -One, not ls -Ell)(On some platforms, like the native sed on OSX and AIX, sed doesn't have a -i option, so you have to do this through some temp file.)

    - keep all BaryonFields the same except the velocity that I need to set up to 0. Is there an easy way to do that ?

This is pretty easy with Python and h5py. Basically, loop over all the grid files, open them, and re-write the file into a new directory,but make V = 0 as you go. Then copy all the files that aren't *.grid* files to your new directory. (I don't think there's a way to kill a dataset directly, so you need to copy everything. ) Something like: for grid in glob.glob("*.grid"): file1 = h5py.File(grid,'r') file2 = h5py.File(other_dir + grid, 'w') for group in file1.listitems(): open the group. for field in group: if field is not velocity: file2.create_datasete(field, shape, data=file1[group][field) else: file2.create_dataset( field, shape, data=numpy.zeros( file1[group][field].shape) ) with the appropriate syntax improvements. d.

Thanks for your help,

JC

David Collins a écrit :

The issue is I can't trust any simulation done with restarting Enzo. Look at

I wouldn't be so sure that the solution you propose isn't going to suffer from the same problems that you see with restarting. The restarting machinery doesn't do anything fancy, there aren't many places for errors to creep in there. So baring a few things, this is potentially a deeper issue.

One thing to check, which you may have done already, is to ensure that you're writing with the same precision that you're computing with. Ensure that the top 5 things in make show-config are all 64 (or 32, if you must) CONFIG_PRECISION: 64 CONFIG_PARTICLES: 64 CONFIG_INTEGERS: 64 CONFIG_INITS: 64 CONFIG_IO: 64 (Inits doesn't matter for you, but the last one, CONFIG_IO, is important)

Have you checked that restarting without changing anything at all gives different answers? What compiler optimizations are you using? Try with -O1 or -O0, sometimes more aggressive compiler options can change the answer.

d.

those two plots: http://drop.io/0etjfmc. I don't want to bother you with too many details but they represent the separation between my two particles - the core of the RGB star and the compact companion.

The first run was with dumping every 80 cycles - dumping based on CycleSkip. It stopped at dumping #196. I restarted it and got data #197, #198,...Then, I wanted to check that restarting did not mess things up, so I run the same sim again but with a dumping frequency twice smaller. So data #98 and #99 correspond to the same cycle as file #196 and #198 respectively. Files #98 and #196 are exactly the same (which is expected) but #99 and #198 are not. They are slightly different and lead at the end to a unphysical behavior - the companion bouncing outwards instead of sinking deeper.

I don't what is causing that, may be the new structure of particles I am using. Has anyone already had that kind of problem ? I will inspect the data at some point in the future but in the meantime, I cannot restart Enzo.

Cheers,

JC

Greg Bryan a écrit :

Or if you just want to read the data directly back in, why not do a restart ("enzo -r")?

Greg

On Jun 2, 2010, at 3:36 PM, David Collins wrote:

Hi, JC--

Would it be easier to skip YT all together, and just have your problem generator read the HDF5 file produced by the 1d run directly? This would streamline your processes, and you wouldn't have the noise incurred by both the ASCII file format and the two units conversions.

I would probably approach your task by looping over only the active zones in enzo, and incrementing an extra counter.

int counter = 0; for(int field=0;field<NumberOfBaryonFields;field++) for(int k=GridStartIndex[2];k<=GridEndIndex[2]; k++) for(int j=GridStartIndex[1];j<=GridEndIndex[1]; j++) for(int i=GridStartIndex[0];i<=GridEndIndex[0]; i++){ index = i + GridDimension[0]*(j + GridDimension[1]*k); BaryonField[ DensNum ][ i ] = DensityFromFile[ counter ]; counter++ }

This code snippit came from the wiki page on Baryon access, for future reference: http://lca.ucsd.edu/projects/enzo/wiki/Tutorials/BaryonFieldAccess

Does that help? d.

On Wed, Jun 2, 2010 at 12:07 PM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

Hi Matt,

I will try to explain what my goal is. I have a 1D profile of a star that is physically at the equilibrium. Regarding the problem initialization, I set up the density and the total specific energy fields according to the model, velocities equal 0. However, the star is not numerically stable so I have to relax it. Therefore, I divide the velocity by 2 at each timestep, and I let Enzo run for a few dynamical times. Then, I read the last data dump and create an relaxed_model.dat:

########################################################## region = pf.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]) output = 'relaxed_model.dat' f = open(output,'w') print >>f, '# Grid Index rho pressure' density = region["Density"] pressure = region["Pressure"] size = NP.size(density) for i in range(0,size,1): dtmp = density[i] ptmp = pressure[i] print >>f, "0", i, '%06e' % dtmp, '%06e' % ptmp ##########################################################

so I end up with a file like this:

########################################################## # Restart file for only one grid # Grid Index rho pressure 0 0 5.099011e-10 2.568020e+05 0 1 5.099011e-10 2.568020e+05 0 2 5.099011e-10 2.568020e+05 ... ########################################################## Finally, I read index, density, pressure from that previous file and do:

########################################################## for (k = 0; k < GridDimension[2]; k++) for (j = 0; j < GridDimension[1]; j++) for (i = 0; i < GridDimension[0]; i++) { index = i + GridDimension[0]*(j + GridDimension[1]*k); density /= DensityUnits; pressure /= PressureUnits; BaryonField[0][index] = density; BaryonField[1][index] = pressure / ((Gamma - 1.0) * density); ########################################################## The problem is GridDimensions contain the ghost zones so the variable index does not match with the actual index read in the file relaxed_model.dat. That is why I wanted to have that file values for the ghost zones as well.

Does it make sense ? Do you have any suggestion ?

Thanks for your help,

JC

Matthew Turk a écrit :

Hi Jean-Claude,

There are a couple aspects to this. The first is that Enzo doesn't output the ghost zones -- so any ghost zones handled inside yt are generated by yt. Were Enzo to output the ghost zones, we would probably be able to handle this, but it does not.

Derived fields can depend on the generation of ghost zones, but keep in mind that these are ghost zones generated by yt. These ghost zones are constructed in a similar manner to how Enzo generates them, but there may be minor differences. You can manually inspect ghost zones on a *grid* by calling retrieve_ghost_zones on that grid.

If you could tell us a bit more about your goal, maybe we could help out a bit more?

-Matt

On Wed, Jun 2, 2010 at 9:32 AM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

Dear all,

in order to set up my simulation, I need to access the values of certain fields in the ghost zones. If I do something like:

pf = load(data) region = pf.h.region(...) x = region[field]

x contains the values of field for the physical grid only. Is there a way I can get the same 1D-array but with the ghost zones included as well ?

Thanks for your help,

Jean-Claude

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

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

-- Sent from my Stone Tablet and carried by my Pterodactyl. _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

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

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

-- Sent from my Stone Tablet and carried by my Pterodactyl.

There is no more *.grid* file. That was for Enzo 1.0 if I am not wrong. You probably meant the *.cpu* files, right ?

That's correct.

To go further, is there a way I can:

    - restart Enzo with WritePotential = 1 using a file that had WritePotential = 0 ?     - restart Enzo and use Tracer Particles with a file that did not have any of them ?

Both are possible, but you have to rewrite the datasets as well a the supporting meta data (boundary and hierarchy). I don't know the details about writing tracer particles, but very recently I wrote a script to add the potential field to a dataset. Unfortunately, the script was for my version of Enzo, which doesn't use packed AMR and has extra fields for my MHD implementation. I've attached a draft of a code that will do the same thing for the more recent Enzo 1.5 with the PackedAMR version, though I don't have time right now to do much more. Basically, you need to have a GravPotential field in the hdf5 files, it needs to have a label in the parameter file, it needs to have it's field_type listed and be counted in NumberOfBaryonFields in the hierarchy file and the boundary file, and (if you're not using OutOfCore Boundary) needs to have space in the hdf5 boundary file. There are a lot of things to get right. Alternatively, you could start a simulation with the desired size and processor layout, and use h5py to write the density, energy, and velocity fields into that. That will ensure the hierarchy and boundary are constructed properly. Easiest still would be to just start the initial relax sims with the gravitational potential written, that will have everything set up right. If anyone out there has an easier solution, that would be sweet. d.

Thanks again for your help,

JC

-- Sent from my Stone Tablet and carried by my Pterodactyl.

Hi Jean-Claude, David, I think this is better suited to enzo-users-l. Best, Matt On Tue, Jun 15, 2010 at 7:11 AM, Jean-Claude Passy <jcpassy@gmail.com> wrote:

Hi David,

I was able modify the outputs as I wanted. Again, you were right: it took me a little bit of time but it was totally worth it.

Now, let me ask you another quick question. I want to add a particle to the data. Typically, I consider the relaxed data of my primary model, set up the velocity fields to zero, add a particle describing the companion and start my "real" simulation. Therefore, I need to modify the fields 'particle...' accordingly and it should be pretty easy now that I know how to use h5py.

My only problem is that I have created my own particles type and I need to specify it. In macros_and_parameters.h:

#define PARTICLE_TYPE_POINT_MASS   5

Where/how do I specify the particle type ? I can't find those parameters in the output data...

Thanks a lot,

Jean-Claude

David Collins a écrit :

your guess was totally right. I was computing with 64 bits and writing with 32 bits only. I fixed it and restarting Enzo worked well. Great !

Great, I'm glad that worked. Entire fields of science have been created because of that very issue. It's subtle, but a killer.

- change some global parameters like StopCycle, CycleDataDump, GlobalDir,... Would it work if I only edit the files DDnnnn/CommonEnvelopennnn ? Wouldn't it create a conflict somewhere ?

Changing cycle based outputs won't change the answer. Changing timestep based outputs (like DtDataDump, Redshift) might introduce a little diffusion, because there will usually be a short timestep to make the output time right.

As long as you replace GlobalDir everywhere, you're fine.

ls -1 |grep -v grid | sed -i 's/GlobalDir.*/GlobalDir = $A'

should do it, but double check. (that's ls -One, not ls -Ell)(On some platforms, like the native sed on OSX and AIX, sed doesn't have a -i option, so you have to do this through some temp file.)

- keep all BaryonFields the same except the velocity that I need to set up to 0. Is there an easy way to do that ?

This is pretty easy with Python and h5py. Basically, loop over all the grid files, open them, and re-write the file into a new directory,but make V = 0 as you go. Then copy all the files that aren't *.grid* files to your new directory. (I don't think there's a way to kill a dataset directly, so you need to copy everything. )

Something like:

for grid in glob.glob("*.grid"): file1 = h5py.File(grid,'r') file2 = h5py.File(other_dir + grid, 'w') for group in file1.listitems(): open the group. for field in group: if field is not velocity: file2.create_datasete(field, shape, data=file1[group][field) else: file2.create_dataset( field, shape, data=numpy.zeros( file1[group][field].shape) )

with the appropriate syntax improvements.

d.

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