[Users] 3D Carpet component with one grid point along the z direction

[Erik Schnetter]

On Wed, Apr 14, 2021 at 2:06 PM Gabriele Bozzola
<bozzola.gabriele at 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 at 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 at 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 at einsteintoolkit.org
>> > http://lists.einsteintoolkit.org/mailman/listinfo/users
>>
>>
>>
>> --
>> Erik Schnetter <schnetter at gmail.com>
>> http://www.perimeterinstitute.ca/personal/eschnetter/



-- 
Erik Schnetter <schnetter at gmail.com>
http://www.perimeterinstitute.ca/personal/eschnetter/