Pass pre-defined matrices to FCN subroutine of IVPAG

[febypa]

Hey all,

As the title describes I would like to know if it is possible to pass a matrix to the FCN subroutine of IVPAG. I am trying to solve a stiff ODE whose YPRIME definition contains components from another matrix (of predefined dimension) that is defined in the program unit which calls FCN.

I tried to implement with common blocks. But the array dimension in the FCN subroutine gets initialized with junk values (very large negative dimension). And the output aborts.

I have read here that arrays in common blocks is not recommended.

Could someone point me in the right direction?

Thank You
Feby

[mecej4]

If by "pass" you mean pass an additional argument to the subroutine, the answer is 'No', since the IVPAG routine has no provision to do so.

In lieu of Common blocks, you may declare and define shared variables, such as the matrix that you speak of, in a module, and use that module in FCN.

For example, consider the example in the IVPAG section of the FNL7.1 manual. You can remove the Common blocks by (i) adding a module with the lines:

     MODULE GlobVars
         REAL       HX
     END MODULE

and, (ii) in every subprogram that uses the Common block at present, insert 'USE GlobVars' after the Subroutine, Program or Function line.

[febypa]

Thanks @mecej4 for the suggestion. I will look into it.

A quick question, is MODULE feature available in Fortran77? I am writing a UMAT for Abaqus. So I have to do it in Fortran77.


- Feby

[mecej4]

There are no current Fortran compilers that cannot digest more recent versions of Fortran than Fortran 77.

I am not an Abaqus user, but I am pretty sure that there is no requirement that you use Fortran 77. Abaqus merely calls your UMAT routines in the DLL. Whether the DLL came out of compiling Fortran, Pascal, C or whatever else, including assembler code, as long as the calling convention matches the one Abaqus expects, it should work.

[febypa]

Thanks once again. The MODULE command seems to solve the global scope issue. But now the code fails during runtime. Using a write out command, I see that during execution it enters the IVPAG function. But does not exit successfully. The cmd output window displays the following:

forrt1: severe <32>: invalid logical unit number, unit -858993460, file unknown

And the following Abort window from Visual C++ Runtime Library pops up.

Capture.jpg

I am also attaching the output windows screen-grab.
Capture1.jpg

I have attached the bare-bones code that I am using. Could someone please suggest where I am going wrong?

      
      MODULE GLOBVARS
      REAL ARR0(3,3) , ARR1(3,3) , DELTIME
      END MODULE GLOBVARS
      
      PROGRAM MAIN
      INCLUDE 'link_fnl_static.h'
      
      USE GLOBVARS
      USE IVPAG_INT
      USE UMACH_INT !Routine UMACH sets or retrieves the input, output, or error output device unit numbers.
      USE LINRG_INT
      
      IMPLICIT NONE
      
      INTEGER MXPARM, N
      
      PARAMETER (MXPARM=50, N=4)
      
! SPECIFICATIONS FOR IVPAG PARAMETERS
      INTEGER MABSE, MBDF, MSOLVE
      PARAMETER (MABSE=1, MBDF=2, MSOLVE=2)
      
! SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER IDO, ISTEP, NOUT
      REAL A(1,1), PARAM(MXPARM), T, TEND, TOL, Y(N)
!
      REAL D0(3,3), D1(3,3),DT
!
! SPECIFICATIONS FOR FUNCTIONS
      EXTERNAL FCN, FCNJ
!     
      DATA D0/1.985,0.0,0.0,0.0,0.725,0.0,0.0,0.0,0.7/, 
     1D1/2.0,0.0,0.0,0.0,0.718,0.0,0.0,0.0,0.695/, DT/0.1/  
!     
      ARR0 = D0
      ARR1 = D1
      DELTIME = DT
!
!Routine UMACH sets or retrieves the input, output, or error output device unit numbers.
      CALL UMACH (2, NOUT) !2 - Retrieves and stores the output unit number in NOUT    
! Set initial conditions
      T = 0.0
      Y(1) = 1.0
      Y(2) = 1.0
      Y(3) = 1.0
      Y(4) = 1.0
! Set TEND
      TEND = 10.0
      
! Set error tolerance
      TOL = 0.000001
      
! Set PARAM to defaults
      PARAM = 0.0E0
!
      PARAM(10) = MABSE
! Select BDF method
      PARAM(12) = MBDF
! Select chord method and a divided difference Jacobian.
      PARAM(13) = MSOLVE
      
! Solve & Print
      WRITE (NOUT,99998)
99998 FORMAT (5X, 'Time', 9X, 'Y1', 11X, 'Y2', 11X,'Y3', 
     111X, 'Y4')
      IDO = 1
10    CONTINUE
      WRITE(NOUT,'(A)') 'CALLING IVPAG'
      CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y, TOL=TOL, PARAM=PARAM)
      WRITE(NOUT,'(A)') 'EXITING IVPAG'
      IF (T .LE. TEND .AND. IDO .NE. 3) THEN
          WRITE (NOUT,'(5F12.6)') T, Y          
          IF (T .EQ. TEND) THEN 
              IDO = 3
              GO TO 10
          END IF
      END IF
      
! Show number of function calls.
      WRITE (NOUT,99999) PARAM(35)
99999 FORMAT (4X, 'Number of fcn calls with IVPAG =', I6)
      END
      
      SUBROUTINE FCN (N, T, Y, YPRIME)
      USE GLOBVARS
! SPECIFICATIONS FOR ARGUMENTS
      INTEGER N
      REAL T, Y(N), YPRIME(N), A(3,3), B(3,0), INV(3,3), X(3,0), 
     1FDOT(3,3), FINV(3,3), A1(3,3)
!
      WRITE(NOUT,'(A)') 'ENTERING SUBROUTINE FCN'
!      
      FDOT = (ARR1-ARR0)/DELTIME
      A1 = ARR1
! 
      CALL LINRG (A1, AINV)
      FINV = AINV
!
      YPRIME(1) = (FDOT(1,1)*FINV(1,1)*Y(1) + FDOT(1,2)*FINV(2,1)*Y(1) 
     1+FDOT(1,2)*FINV(2,2)*Y(4)+FDOT(1,3)*FINV(3,1)*Y(1))
     2+0.5*((3/((Y(2)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))
     3+(Y(1)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))+(1/Y(3)))) - Y(1))
!        
      YPRIME(2) = (FDOT(2,1)*FINV(1,1)*Y(4) + FDOT(2,1)*FINV(1,2)*Y(2)
     1+FDOT(2,3)*FINV(3,1)*Y(4) + FDOT(2,3)*FINV(3,2)*Y(2))
     2+0.5*((3/((Y(2)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))
     3+(Y(1)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))+(1/Y(3)))) - Y(2))
!      
      YPRIME(3) = (FDOT(3,1)*FINV(1,3)*Y(3) + FDOT(3,2)*FINV(2,3)*Y(3)
     1+FDOT(3,3)*FINV(3,3)*Y(3))
     2+0.5*(3/((Y(2)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4)))
     3+(Y(1)*Y(3)/(Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4))) + (1/Y(3))) -Y(3))
!
      YPRIME(4) = (FDOT(1,1)*FINV(1,1)*Y(4) + FDOT(1,1)*FINV(1,2)*Y(2)
     1+FDOT(1,2)*FINV(2,1)*Y(4) + FDOT(1,2)*FINV(2,2)*Y(2)
     2+FDOT(1,3)*FINV(3,1)*Y(4) + FDOT(1,3)*FINV(3,2)*Y(2))
      RETURN
      END
      
      SUBROUTINE FCNJ (N, T, Y, DYPDY)
! SPECIFICATIONS FOR ARGUMENTS
      INTEGER N
      REAL T, Y(N), DYPDY(N,*)
!
      RETURN    
      END

Thank You

- Feby

[mecej4]

You have to debug the code yourself.

You have at least two errors in FCN. Check all the declarations, and check (by comparison with hand calculated results for at least one value of Y) whether FCN returns the correct result in YPRIME for some Y that you select, before attempting to use the ODE solver.

You divide by Y(1)*Y(2)*Y(3) - Y(3)*Y(4)*Y(4) in a couple of places. At the initial point, all the Y are equal to 1.0, so this causes division by zero.