On Wed, Apr 14, 2021 at 2:06 PM Gabriele Bozzola bozzola.gabriele@gmail.com wrote:

Hi Erik,

thank you.

...

To correct the actual bug you noticed, one could apply a correction that is local to OutputSlice.cc. The code there has access to the grid spacing, which is the same for all components, and which could therefore be calculated in the same place where coord_lower and coord_upper are calculated (function "GetCoordinates"), and could then be passed to the routine that outputs the "delta" attribute.

Is there a way to do this cleanly? In my simulation I assume dx=dy=dz, so I can check if coord_delta[d] = 0 and set coord_delta to the correct value as read from the other directions (in the function "GetCoordinates"). But this would only be a workaround. If in the future I forget about this and run with different grid spacings, this would lead to a small error that would be pretty much impossible to find.

Also, when I look at the function I see:

rvect const cctk_delta_space = delta_space.at(m) * rvect(mglevelfact);

what is delta_space? Does the 0 come from here?

No; this is where the good data live. The problem is that this routine doesn't return the "delta" value, only "lower" and "upper", and "delta" is later re-calculated from these. Which is not possible if "lower" and "upper" are the same.

I created a Carpet branch "eschnett/delta" where I jotted down my idea. I did not test it. Could you? In particular, there could be an error in calculating the grid spacing, and it might now be wrong, which should be easy to detect in a small test run.

-erik

Thanks, Gabriele

On Tue, Apr 13, 2021 at 11:45 AM Erik Schnetter schnetter@gmail.com wrote:

...

Gabriele

This is indeed a funny component shape. To my knowledge, certain thorns will not handle it correctly, because they assume that the width of the domain is at least the same as the number of ghost or boundary points. (This assumes you're using more than one ghost and boundary points.)

There is special code in Carpet to handle the case where a component is only one point wide, and which sets the width to 0 in this case. The reason is that Cactus does not store the grid spacing internally, but only the lower and upper bound of the grid. Both are identical in this case, so it is not possible to determine the grid spacing.

In your case, the proper solution would be to ensure that no such component is created, as it is quite inefficient. Unfortunately, the domain decomposition algorithm already has to satisfy many constraints, and modifying it is thus difficult. One way would be to set "Carpet::granularity" to e.g. 4, which forces all grids to have a size that is a multiple of 4. This would also restrict the total domain sizes you can have, so using this is not straightforward.

To correct the actual bug you noticed, one could apply a correction that is local to OutputSlice.cc. The code there has access to the grid spacing, which is the same for all components, and which could therefore be calculated in the same place where coord_lower and coord_upper are calculated (function "GetCoordinates"), and could then be passed to the routine that outputs the "delta" attribute.

-erik

On Tue, Apr 13, 2021 at 11:04 AM Gabriele Bozzola bozzola.gabriele@gmail.com wrote:

...

Hello,

one of my evolutions produced the attached file.

What is interesting about this file is that refinement level 10 has a component that has shape (1, 53, 105). So, we have only one point along the z direction. I am not entirely sure how this came to be, but the problem is that the HDF5 attribute 'delta' reports a cell width of 0 along the z direction:

In [12]: with h5py.File("alp.xyz.file_0.h5") as f: ...: print(f['ADMBASE::alp it=55296 tl=0 rl=10 c=0'].shape) (1, 53, 105)

In [13]: with h5py.File("alp.xyz.file_0.h5") as f: ...: print(f['ADMBASE::alp it=55296 tl=0 rl=10 c=0'].attrs['delta']) [0.072 0.072 0. ]

I would have expected the value of 0.072, even if the cell has only one grid point.

Is this a bug?

Thanks, Gabriele _______________________________________________ Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users

-- Erik Schnetter schnetter@gmail.com http://www.perimeterinstitute.ca/personal/eschnetter/