I'm trying to use some existing code for integrating delay differential equations in a python project (radar5 from this page ). I thought I had successfully wrapped the fortran code with f2py because I was able to call a fortran program and see output from the integrator.
But, all python code after the call to radar5 does not execute.
The python code:
import numpy
import radar5_f2pytest
print(radar5_f2pytest.__doc__)
i=[2]
x=[0]
rpar=numpy.array([1.34, 1.6E9, 8.0E3, 4.0E7, 1, 1, 6E-2, 6E-2, 15E-2],order='F')
ND=2
NRDENS=1
NGRID=1
NLAGS=1
NJACL=2
MXST=4000
LWORK=11
LIWORK=16
IJAC=1
MLJAC=ND
IMAS=0
IOUT=1
X=0
Y=numpy.array([1E-10,1E-5],order='F')
TAU=rpar[8]
XEND=100.5
ITOL=0
RTOL=1E-9
ATOL=RTOL*1E-9
H=1E-6
IWORK=numpy.zeros((16), order='F', dtype=numpy.int64)
WORK=numpy.array([1E-16, 0.9, 0.001, min(0.03, RTOL**0.5), 1, 1.2, XEND-X, 0.2, 8, 0, 5], order='F')
IWORK[1]=1000000
IWORK[4]=ND
IWORK[2]=7
IWORK[8]=1
IWORK[10]=2
IWORK[11]=MXST
IWORK[12]=NGRID
IWORK[13]=1
IWORK[14]=NRDENS
IPAST=numpy.zeros((NRDENS+1), order='F')
IPAST[0]=2
GRID=numpy.array([TAU], order='F')
ipar=numpy.zeros((1))
#print((WORK))
#print(IWORK)
print("WHAT")
radar5_f2pytest.radar5_f2py(X, Y, XEND, H,
RTOL, ATOL, ITOL, IJAC,
MLJAC, NLAGS, NJACL,
IMAS,ND,0, IOUT, WORK,
IWORK,GRID, IPAST, rpar, ipar)
print("WHAT")
The fortran code. Note that there are external dependencies, and other subroutines, but for brevity, I will only include the subroutine that I am wrapping with f2py.
! -*- f90 -*-
SUBROUTINE radar5_f2py(N, X, Y, XEND, H, RTOL, ATOL, ITOL,&
& IJAC, MLJAC, MUJAC, NLAGS, NJACL, &
& IMAS, MLMAS, MUMAS, IOUT, WORK, IWORK, &
& GRID, IPAST, RPAR, IPAR, IDID, RESULT)
IMPLICIT NONE
INTEGER, PARAMETER :: DP=kind(1D0)
INTEGER, intent(in) :: N, NLAGS, NJACL
REAL(kind=8), intent(inout) :: X
REAL(kind=8), intent(in) :: XEND
REAL(kind=8), intent(in) :: H
LOGICAL, intent(in) :: ITOL, IJAC, IMAS, IOUT
INTEGER, intent(in) :: MLJAC, MLMAS, MUMAS
INTEGER, intent(in) :: MUJAC
!f2py integer, optional, intent(in) :: MUJAC
REAL(kind=8), dimension(N), intent(inout) :: Y
REAL(kind=8), intent(in), dimension(11) :: WORK
REAL(kind=8), intent(in) :: ATOL,RTOL
INTEGER, intent(in), dimension(16) :: IWORK
REAL(kind=8), intent(inout) :: GRID
INTEGER, intent(inout) :: IPAST
REAL(kind=8), intent(in), dimension(9) :: RPAR
INTEGER, intent(in) :: IPAR
INTEGER, intent(out) :: IDID
REAL*8, intent(out), dimension(1,N) :: RESULT
REAL(kind=8), EXTERNAL :: PHI
REAL(kind=8), EXTERNAL :: ARGLAG
EXTERNAL :: FCN, JFCN, JACLAG, SOLOUT, DUMMY
!PRINT *, X, Y, WORK, ATOL, RTOL, IWORK, GRID, IPAST, RPAR, IPAR
CALL RADAR5(N,FCN,PHI,ARGLAG,X,Y,XEND,H,&
& RTOL,ATOL,ITOL, &
& JFCN,IJAC,MLJAC,MUJAC, &
& JACLAG,NLAGS,NJACL, &
& IMAS,SOLOUT,IOUT, &
& WORK,IWORK,RPAR,IPAR,IDID,&
& GRID,IPAST,DUMMY,MLMAS,MUMAS,RESULT)
RESULT(1,1)=1
!PRINT *, X, IDID
RETURN
END SUBROUTINE
The f2py module was compiled with this command:
flang -c radar5.f contr5.f dc_decdel.f decsol.f dontr5.f
================================================================
start f2py...
f2py -c -m radar5_f2pytest radar5_template.f90 radar5.o contr5.o dc_decdel.o decsol.o dontr5.o only: radar5_f2py :
I first compile the fortran source files and then feed them to f2py with a source file containing the subroutine above.
To reiterate, I can see output from fortran indicating that the call to RADAR5 is working. And the results match examples given by the authors of radar5.
I have tried debugging the python code, but that merely confirmed that nothing after the second call to the f2py module executed.
EDIT:
Sorry for the delay and the lack of a reproducible example. I have been trying to make an example that shows this issue without breaking something else, but radar5 is quite large and complicated, and I am inexperienced with fortran. The issue is related to the call to radar5, since if I do not call radar5 at all, everything works as expected. In the meantime, here is the output:
This module 'radar5_f2pytest' is auto-generated with f2py (version:1.21.5).
Functions:
idid,result = radar5_f2py(x,y,xend,h,rtol,atol,itol,ijac,mljac,nlags,njacl,imas,mlmas,mumas,iout,work,iwork,grid,ipast,rpar,ipar,n=len(y),mujac=)
.
WHAT
STARTING INTEGRATION...
NUMBER OF PRESCRIBED GRID POINTS: 1
NUMBER OF DELAYED COMPONENTS: 1
INTEGRATION...
0 COMPUTED BREAKING POINTS:
Here, I know that radar5 is being called and is executing because of the last line. I can also print out results (from fortran) like IDID or the final point of the integration. These results have been as expected.
EDIT 2:
Running the python code in an interactive window causes a crash:
Version details:
Python 3.9.13 (main, Aug 25 2022, 23:51:50) [MSC v.1916 64 bit (AMD64)]
Type 'copyright', 'credits' or 'license' for more information
IPython 7.31.1 -- An enhanced Interactive Python. Type '?' for help.
Output in interactive window:
Canceled future for execute_request message before replies were done
The Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click here for more info. View Jupyter log for further details.
