Incompact3d icon indicating copy to clipboard operation
Incompact3d copied to clipboard

IBM - Undefined values in ypraf in case of periodicity

Open mathrack opened this issue 1 year ago • 0 comments

The size of the array ypraf is nyraf = nym nraf. In case of periodicity, nym=ny. Otherwise nym=ny-1.

The loop used to define ypraf will populate the array from 1 to (ny-1) nraf.

https://github.com/xcompact3d/Incompact3d/blob/46e874e236265c2f2671f254a9e3f4415d1e40d0/src/genepsi3d.f90#L228

In case of periodicity, the array will contain nraf undefined values at the end. This is an issue because it can lead to floating point exceptions when computing yepsi

https://github.com/xcompact3d/Incompact3d/blob/46e874e236265c2f2671f254a9e3f4415d1e40d0/src/genepsi3d.f90#L233

When the mesh is regular, one can simply defined ypraf using

    do j = 1, nyraf
          ypraf(j) = (j-1) * dyraf
    enddo

In case of stretching, I think one would have to use a subroutine similar to stretching

https://github.com/xcompact3d/Incompact3d/blob/46e874e236265c2f2671f254a9e3f4415d1e40d0/src/tools.f90#L1215

As a temporary workaround, I suggest to implement the simple loop to fix the issue in case of periodicity without stretching. In case of stretching combined with periodicity, we could keep the current code, print a warning and populate the end of the array ypraf with yly + (j-(ny-1)*nraf)*dyraf or something similar.

mathrack avatar Sep 11 '24 15:09 mathrack