Bug: STRIDE ZVODE integration produces NaN due to stack overflow in `ode_run`

[d-burg]

Branch: slayer_growthrate
File: ode.F

Description

STRIDE failed silently with ZVODE reporting NaN values and a negative istate error code, producing garbage output. The root cause was stack memory corruption from a large automatic (stack-allocated) array inside ode_run.

Root Cause

identityMat was declared as a fixed-size automatic array:

COMPLEX(r8), DIMENSION(2*mpert,2*mpert) :: identityMat

With mpert=34 (a typical value), this allocates a 68×68 complex(r8) array (~72 KB) on the stack every time ode_run is called. Combined with other large automatic arrays already on the stack in that scope, the total stack usage exceeded the system stack limit, which silently corrupted memory. ZVODE subsequently operated on corrupted data, which produced NaNs.

The corruption manifested as identityMat having a max element value of 2.22775E+08 immediately after initialization.

Fix

Changed identityMat from a stack-allocated automatic array to a heap-allocated ALLOCATABLE array:

! Before
COMPLEX(r8), DIMENSION(2*mpert,2*mpert) :: identityMat

! After
COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: identityMat
...
ALLOCATE(identityMat(2*mpert,2*mpert))
! ... use identityMat ...
DEALLOCATE(identityMat)

Result

After the fix, STRIDE produces correct results:

plasma =  -8.812E-01
vacuum =   2.138E+00
total  =   1.257E+00

Now for my question: should all other vulnerable automatic arrays be converted?

I can see other large automatic arrays in STRIDE whose sizes depend on runtime parameters (mpert, msing, mthsurf, etc.). The ones inside routines called during ODE integration or per-rational-surface may be at the highest risk. Should these be systematically converted to ALLOCATABLE?

ode.F — highest risk (called inside ZVODE integration loop):

Variable(s)	Type/Dimension	Approx. size (mpert=34)
uFM, duFM	COMPLEX(r8)(2*mpert,2*mpert)	~72 KB each
uFMInv	COMPLEX(r8)(2*mpert,2*mpert)	~72 KB
u_i, tM (×2 locations)	COMPLEX(r8)(2*mpert,2*mpert)	~72 KB each
u_tot, A, tM, umat	COMPLEX(r8)(2*mpert,2*mpert)	~72 KB each
X0, eV	COMPLEX(r8)(mpert,mpert)	~18 KB each
ua	COMPLEX(r8)(mpert,2*mpert,2)	~36 KB

riccati.f — high risk (called per Riccati integration step):

Variable(s)	Type/Dimension	Approx. size
ac (×2 subroutines)	COMPLEX(r8)(2*mpert,2*mpert)	~72 KB each
pmat, pmat1, pmat1t	COMPLEX(r8)(mpert,mpert)	~18 KB each

free.f — moderate risk (called per output step):

Variable(s)	Type/Dimension	Approx. size
temp, tempUIn, tempU, tempVt	COMPLEX(r8)(mpert,mpert)	~18 KB each
u, du	COMPLEX(r8)(mpert,mpert,2)	~36 KB each
xin, xis	COMPLEX(r8)(mpert,mpert)	~18 KB each
jqvec	COMPLEX(r8)(mpert,mpert,3)	~55 KB
bvec0, bvec	COMPLEX(r8)(mpert,0:mthsurf,3)	mthsurf-dependent

sing.F — moderate risk:

Variable(s)	Type/Dimension	Approx. size
v	COMPLEX(r8)(mpert,2*mpert,2)	~36 KB
idMat	COMPLEX(r8)(mpert,mpert)	~18 KB

The cumulative stack pressure from just the ode.F inner-loop scope alone can exceed 600 KB at mpert=34. A global audit and conversion to ALLOCATABLE (or SAVEd module-level arrays) would make the code robust against stack overflows at higher mpert values.

Other notes from Claude:

• This class of bug is difficult to detect because neither the compiler nor the runtime necessarily raises an error — the array is simply allocated over whatever happens to be adjacent on the stack.
• The default stack size on macOS is 8 MB; Linux HPC systems often set it lower or impose per-thread limits.
• Any large array whose size depends on a runtime parameter inside a frequently-called subroutine should be declared ALLOCATABLE.

[logan-nc]

@d-burg are you going to knock out the high risk ones in #173 then? Or do you want a separate branch just to address this issue?

It's pretty safe to allocate as long as you are sure to deallocate appropriately - so go for it.