search

ivan-pi / pdecheb

Chebyshev Polynomial Software for Elliptic-Parabolic Systems of PDEs
2 stars 0 forks source link

Aliasing in UVINIT #7

Open ivan-pi opened 1 year ago

ivan-pi commented 1 year ago

The V argument is passed as an array of size 1.

https://github.com/ivan-pi/pdecheb/blob/f4294351471375835c79cbfbe508eb8e2f1b39b0/pdecheb/cset.f#L1-L11

...

https://github.com/ivan-pi/pdecheb/blob/f4294351471375835c79cbfbe508eb8e2f1b39b0/pdecheb/cset.f#L122

This can have some surprising effects, if NV = 0, but the solution array would be declared as Y(NPDE*NPTS + 1) in the calling program. Here's the initialization routine that exposes the issue:

SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV, V)
IMPLICIT NONE
INTEGER :: NPDE, NPTS, NV
DOUBLE PRECISION :: X(NPTS), U(NPDE,NPTS), TIME, V(NV)
U(1,:) = 0
U(1,NPTS) = 1
V(1) = 0       ! <-- eliminates effect of previous line
END SUBROUTINE

In practice this should be flagged as an out-of-bounds violation (or produce a segfault).

ivan-pi commented 1 year ago

The problem is more obvious in INICHB where CSET gets called:

      NVST = NPDE*NPTS + 1
C ...
      IV = NPDE*NPTS
      IF (NV.GT.0) THEN
         IV = NVST
C ...
      END IF
      CALL CSET(NPDE,NPTS,U,WK(I6),WK,WK(I2),WK(I5),NEL,NPTL,WK(I4),
     *          WK(I12),XBK,IBK,WK(I3),U(IV),NV)

If NV = 0, the element U(IV) will be the alias for U(NPDE*NPTS) or U(1,NPTS) in the UVINIT routine.