[yt-users] fields in ghost zones

David Collins dcollins at physics.ucsd.edu
Sun Jun 6 15:59:58 PDT 2010


>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 at 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 at 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 at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
>
>
> _______________________________________________
> yt-users mailing list
> yt-users at 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 at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
>
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
>
>
>



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



More information about the yt-users mailing list