19.9.3. Junction Box Example 3: Timestep Control

This Junction Box subroutine should be called at the end of each timestep in a transient run. It is designed to dynamically control the timestep size during the transient run. The algorithm is similar to one of the adaptive timestepping algorithms available in the solver, but serves to illustrate how a junction box can be used to control the timestep size. The algorithm for modifying the timestep is as follows: if the previous timestep needed more iterations than MAXLOOPS, the next timestep size is decreased by a factor FACDEC. If the previous timestep needed fewer iterations than MINLOOPS, the next timestep size is increased by a factor FACINC. The junction box calls USER_GET_TRANS_INFO to retrieve current timestep data and creates a data area called ‘NEWDT’ to store the next timestep size. The procedure requires a User CEL function called Retrieve_deltat to use the timestep calculated by the junction box. Retrieve_deltat requires an additional argument (passed in through the argument list) called ‘Timesteps’, which will be the initial timestep size. The routines t_step_control.F and Retrieve_deltat.F are supplied in the <CFXROOT>/examples/UserFortran/ directory of your CFX installation.

The following is a listing of t_step_control.F:

#include "cfx5ext.h"
dllexport(tstepcontrol)
C================================================================
      SUBROUTINE tstepcontrol(CZ,DZ,IZ,LZ,RZ)
CC
CD Recalculates timestep size to fit into a maximum and minimum
CD coefficient loops.
CC
CC --------------------
CC        Modified
CC --------------------
CC
CV Stacks
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC --------------------
CC        Output
CC --------------------
CC
CC --------------------
CC Modification History
CC --------------------
CC
CC
CC --------------------
CC Variable dictionary
CC --------------------
CC
CC --------------------
CC Local variables & parameters
CC --------------------
CC
CC   CCLOOP,LCLOOP,NCLOOP    : Coefficient loop name,length, and number
CC   CSTEP,LTSTEP,NTSTEP     : Coefficient step name,length, and number
CC   DT                      : Timestep size
CC   MAXDT,MINDT             : Maximum and minimum timestep size
CC   FACINC,FACDEC           : Increment and decrement factors
CC
CC======================================================================
CC ------------------------------
C     Preprocessor includes
C ------------------------------
C
#include "MMS.h"
#include "stack_point.h"
#include "cfd_constants.h"
C
C ------------------------------
C        Argument List
C ------------------------------
C
      CHARACTER CZ(*)*(1)
C
      DOUBLE PRECISION DZ(*)
C
      INTEGER IZ(*)
C
      LOGICAL LZ(*)
C
      REAL RZ(*)
C
C ------------------------------
C        External routines
C ------------------------------
C     
      LOGICAL LFIL      
      EXTERNAL CFROMR,CFROMI,LFIL
C
C ------------------------------
C        Local Variables
C ------------------------------
C
      CHARACTER*4 CRESLT
      CHARACTER*(10) CFROMI
      CHARACTER*(20) CCLOOP,CTSTEP,CFROMR,FNAME
C
      INTEGER LCLOOP,NCLOOP,LTSTEP,NTSTEP,MINLOOPS,MAXLOOPS,IUNIT,ICALL
C      
      REAL DT,MAXDT,MINDT,TIME,FACINC,FACDEC
C
C ------------------------------
C        Local Parameters
C ------------------------------
C
      PARAMETER(IUNIT=91,FNAME='dt_ncloop.dat')
C
C--Saved data.
C      
      SAVE ICALL,FACINC,FACDEC,MINDT,MAXDT,MINLOOPS,MAXLOOPS
C
C ------------------------------
C        Stack pointers
C ------------------------------
C
      __stack_point__ pOINT
C
C ---------------------------
C    Executable Statements
C ---------------------------
C
      IF (ICALL.EQ.0) THEN
        ICALL=1
C
C---- Set defaults:
C     MINDT is used to enforce a positive timestep size
C
        MINDT = 0.
        MAXDT = 1.E10
        MINLOOPS=3
        MAXLOOPS=8
C
C---- Negative means: increase by 10/decrease by 0.1 each given 
C     iteration number, positive means: factor itself
C
        FACINC=-40.
        FACDEC=-10.
C
C---- Get from USER area if there:
C
        CALL USER_PEEKR('MINDT',1,MINDT,'SKIP',CRESLT,CZ)
        CALL USER_PEEKR('MAXDT',1,MAXDT,'SKIP',CRESLT,CZ)
        CALL USER_PEEKI('MINLOOPS',1,MINLOOPS,'SKIP',CRESLT,CZ)
        CALL USER_PEEKI('MAXLOOPS',1,MAXLOOPS,'SKIP',CRESLT,CZ)
        CALL USER_PEEKR('FACINC',1,FACINC,'SKIP',CRESLT,CZ)
        CALL USER_PEEKR('FACDEC',1,FACDEC,'SKIP',CRESLT,CZ)
C
C---- Increment by factor 10 each ## time steps
C
        IF (FACINC.LE.0) FACINC=10.**(1./ABS(FACINC))
C
C---- Decrement by factor 0.1 each ## time steps
C
        IF (FACDEC.LE.0) FACDEC=0.1**(1./ABS(FACDEC))
C
C---- Output parameters
C
        CALL MESAGE('WRITE','MINDT is '//CFROMR(MINDT))
        CALL MESAGE('WRITE','MAXDT is '//CFROMR(MAXDT))
        CALL MESAGE('WRITE','MINLOOPS is '//CFROMI(MINLOOPS))
        CALL MESAGE('WRITE','MAXLOOPS is '//CFROMI(MAXLOOPS))
        CALL MESAGE('WRITE','FACINC is '//CFROMR(FACINC)//
     &              ' ('//CFROMI(INT( 1./LOG10(FACINC)+0.5))//')')
        CALL MESAGE('WRITE','FACDEC is '//CFROMR(FACDEC)//
     &              ' ('//CFROMI(INT(-1./LOG10(FACDEC)+0.5))//')')
        IF (IUNIT>0) OPEN(UNIT=IUNIT,FILE=FNAME)
      END IF
C
C---- Get transient info
C
      CALL USER_GET_TRANS_INFO('TSTEPCONTROL','GET',' ',
     &                         CZ,DZ,IZ,LZ,RZ)
C
C---- Get time step
C
      CALL PEEKR('/USER/TRANS_INFO/DT',1,DT,'STOP',CRESLT,RZ)
C
C---- Write into file
C
      IF (IUNIT>0) THEN
C
C---- Get time and current coefficient loop
C
        CALL PEEKR('/USER/TRANS_INFO/ATIME',1,TIME,'STOP',CRESLT,RZ)      
        CALL PEEKI('/USER/TRANS_INFO/NCLOOP',1,NCLOOP,'STOP',CRESLT,IZ)
        WRITE(IUNIT,*) DT,NCLOOP,TIME
C
      END IF
C
C---- Control structure
C
      CALL MESAGE('WRITE','Number of coeff loops was: '//CFROMI(NCLOOP))
C
      IF (NCLOOP.LT.MINLOOPS) THEN
        DT=MIN(FACINC*DT,MAXDT)
        CALL MESAGE('WRITE','Time step increased to: '//CFROMR(DT))
      ELSE IF (NCLOOP.GT.MAXLOOPS) THEN
        DT=FACDEC*DT
        IF (MINDT.LE.0.AND.DT.LE.-MINDT) THEN
          IF (IUNIT.GT.0) CLOSE(IUNIT)
          CALL MESAGE('WRITE','Stop program due to bad convergence.')
          OPEN(UNIT=91,FILE='stp')
          WRITE(91,*) ' '
          CLOSE(91)
        END IF
          DT=MAX(DT,MINDT)
          CALL MESAGE('WRITE','Time step decreased to: '//CFROMR(DT))
      ELSE
        CALL MESAGE('WRITE','Time step stayed at: '//CFROMR(DT))
      END IF
C
C---- Set new time step
C
      IF (.NOT.LFIL('/USER/NEWDT')) THEN
        CALL PSHDIR ('/USER','STOP',CRESLT)
          CALL MAKDAT('NEWDT','REAL','STOP',1,pOINT,CRESLT)
        CALL POPDIR ('STOP',CRESLT)          
      ENDIF
      CALL POKER('/USER/NEWDT',1,DT,'STOP',CRESLT,RZ)
C
C---- Release transient data
C
      CALL USER_GET_TRANS_INFO('TSTEPCONTROL','RELEASE',' ',
     &                         CZ,DZ,IZ,LZ,RZ)
C     
      END

The following is a listing of Retrieve_deltat.F:

#include "cfx5ext.h"
dllexport(get_timestep)
      SUBROUTINE get_timestep (
     &  NLOC,NRET,NARG,RET,ARGS,CRESLT,CZ,DZ,IZ,LZ,RZ)
CC
CD User routine: retrieves delta t stored in /USER directory
CC
CC --------------------
CC        Input
CC --------------------
CC
CC  NLOC   - size of current locale
CC  NRET   - number of components in result
CC  NARG   - number of arguments in call
CC  ARGS() - (NLOC,NARG) argument values
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC  Stacks possibly.
CC
CC --------------------
CC        Output
CC --------------------
CC
CC  RET()  - (NLOC,NRET) return values
CC  CRESLT - 'GOOD' for success
CC
CC --------------------
CC       Details
CC --------------------
CC
CC======================================================================
C
C ------------------------------
C     Preprocessor includes
C ------------------------------
C
#include "MMS.h"
#include "cfd_constants.h"
C
C ------------------------------
C        Global Parameters
C ------------------------------
C
C
C ------------------------------
C        Argument list
C ------------------------------
C
      CHARACTER CRESLT*(*)
      CHARACTER CZ(*)*(1)   
C
      DOUBLE PRECISION DZ(*)
C
      INTEGER IZ(*)
      INTEGER NLOC,NARG,NRET
C
      LOGICAL LZ(*)
C
      REAL RZ(*)
      REAL ARGS(NLOC,NARG), RET(NLOC,NRET)
C
C
C ------------------------------
C        External routines
C ------------------------------
C
C
C ------------------------------
C        Local Parameters
C ------------------------------
C
C
C ------------------------------
C        Local Variables
C ------------------------------
C
      CHARACTER*(4) CRES
C
      REAL  DT
C
C ---------------------------
C    Executable Statements
C ---------------------------
C
C
C Set success flag.
C
      CRES = 'GOOD'
C
C---- Retrieves DT value. If there is none, it uses default the value
C     passed as argument.
C      
      CALL PEEKR('/USER/NEWDT',IONE,DT,'SKIP',CRES,RZ)
      IF (CRES .NE. 'GOOD') DT = ARGS(1,1)
      RET(1,1) = DT
C
      END