Results 1 to 6 of 6

Thread: Pass pre-defined matrices to FCN subroutine of IVPAG

  1. #1
    Junior Member
    Join Date
    Sep 2016
    Location
    Chennai, India
    Posts
    3

    Pass pre-defined matrices to FCN subroutine of IVPAG

    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

  2. #2
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    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:
    Code:
         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.
    Last edited by mecej4; 10-08-2016 at 10:01 AM.

  3. #3
    Junior Member
    Join Date
    Sep 2016
    Location
    Chennai, India
    Posts
    3
    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

  4. #4
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    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.

  5. #5
    Junior Member
    Join Date
    Sep 2016
    Location
    Chennai, India
    Posts
    3
    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:
    Code:
    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?

    Code:
          
          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

  6. #6
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    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.

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •