CalculateCTFs Subroutine

private subroutine CalculateCTFs(nrf, SolutionDimensions)

Arguments

Type	Intent	Optional	Attributes		Name
integer,	intent(out)			::	nrf
integer,	intent(in)			::	SolutionDimensions
Called By

Help
Source Code

CalculateCTFs
Source Code

SUBROUTINE CalculateCTFs(nrf,SolutionDimensions)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Russ Taylor
          !       DATE WRITTEN   June 1990
          !       MODIFIED       Apr 96, RKS, cosmetic, algorithm neutral changes
          !       RE-ENGINEERED  July 1996, RKS; Nov 1999, LKL (allocatable arrays)

          ! PURPOSE OF THIS SUBROUTINE:
          ! Subprogram to calculate the Sj and ej coefficients of Seem's
          ! dissertation.  Follows Seem's technique to compute the coefficients
          ! in order with minimum storage requirements.

          ! METHODOLOGY EMPLOYED:
          ! Combine the results of the ExponMatrix, InvertMatrix, and
          ! CalculateGammas routines together to arrive at the temperature
          ! coefficients (s, s0) and the heat flux history coefficients (e) of
          ! the CTFs.  The outline of this subroutine is based on step 5 of
          ! Seem's suggested implementation of the state space method found on
          ! pages 26+27 of his dissertation.

          ! REFERENCES:
          ! The state space method of calculating CTFs is outlined in the
          ! doctoral dissertation of John Seem, "Modeling of Heat Transfer in
          ! Buildings", Department of Mechanical Engineering, University of
          ! Wisconsin-Madison, 1987.  In particular, the equations used for
          ! these calculations are equations (2.1.24) through (2.1.26) in Seem's
          ! dissertation.

          ! USE STATEMENTS:
          ! none

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(OUT) :: nrf                ! Number of response factor terms
  INTEGER, INTENT(IN)  :: SolutionDimensions ! Integer relating whether a 1- or 2-D solution is required

          ! SUBROUTINE PARAMETER DEFINITIONS:
  REAL(r64), PARAMETER :: ConvrgLim = 1.0d-13   ! Convergence limit (ratio) for cutting off the calculation of further
                                                       ! CTFs.  This value was found to give suitable accuracy in IBLAST.

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE DUMMY VARIABLE DECLARATIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

  REAL(r64)                                :: avg       ! Intermediate calculation variable (average)
  LOGICAL                                              :: CTFConvrg ! Set after CTFs are calculated, based on whether there are
                                                                    ! too many CTFs terms
  INTEGER                                              :: i         ! Loop counter
  INTEGER                                              :: ic        ! Loop counter
  INTEGER                                              :: inum      ! Loop counter
  INTEGER                                              :: ir        ! Loop counter
  INTEGER                                              :: is        ! Loop counter
  INTEGER                                              :: is2       ! Loop counter
  INTEGER                                              :: j         ! Loop counter
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)        :: PhiR0     ! Product of Phi( = AExp) and R0 matrices from the state
                                                                    ! space method
  REAL(r64)                                :: rat       ! Intermediate calculation variable (ratio of flux history
                                                                    ! terms)
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)        :: Rnew      ! Current R matrix
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)        :: Rold      ! R matrix from the last iteration
  INTEGER                                              :: SurfNode  ! Loop counter (for nodes at a surface)
  REAL(r64)                                :: SurfNodeFac ! Multiplying factor applied to various surface nodes
  REAL(r64)                                :: trace     ! Trace of the product of Phi( = AExp) and R0

          ! FLOW:

          ! Subroutine initializations
  ALLOCATE(PhiR0(rcmax,rcmax))
  ALLOCATE(Rnew(rcmax,rcmax))
  ALLOCATE(Rold(rcmax,rcmax))
  PhiR0 = 0.0D0
  Rold  = 0.0D0

  s0   = 0.0D0
  s    = 0.0D0
  e    = 0.0D0
  Rnew = IdenMatrix ! Rnew initialized to the identity matrix

          ! Calculate Gamma1-Gamma2.  Gamma1 is not used by itself in the
          ! equations, only Gamma1-Gamma2.  Thus, reset Gamma1 to:
          ! Gamma1-Gamma2
  DO i = 1,rcmax
    DO j = 1,3
      Gamma1(i,j) = Gamma1(i,j)-Gamma2(i,j)
    END DO
  END DO

          ! Compute s0.  See Seem's thesis equation (2.1.24) which states that:
          ! s0  =  (CMat*R0*Gamma2) + (DMat)
          ! Note that for a two-dimensional solution, there is more than one
          ! node at the surface and the effect of each of these must be added
          ! together.
  IF (SolutionDimensions == 1) THEN
    s0(1,1) = CMat(1)*Gamma2(1,1) + DMat(1)
    s0(1,2) = CMat(1)*Gamma2(1,2)
    s0(1,3) = CMat(1)*Gamma2(1,3)
    s0(2,1) = CMat(2)*Gamma2(rcmax,1)
    s0(2,2) = CMat(2)*Gamma2(rcmax,2) + DMat(2)
    s0(2,3) = CMat(2)*Gamma2(rcmax,3)
  ELSE  ! SolutionDimensions = 2
    DO SurfNode = 1, NumOfPerpendNodes
      IF ( (SurfNode == 1) .OR. (SurfNode == NumOfPerpendNodes) ) THEN
        SurfNodeFac = 0.5d0
      ELSE
        SurfNodeFac = 1.0d0
      END IF
      s0(1,1) = s0(1,1) + SurfNodeFac*CMat(1)*Gamma2(SurfNode,1)
      s0(1,2) = s0(1,2) + SurfNodeFac*CMat(1)*Gamma2(SurfNode,2)
      s0(1,3) = s0(1,3) + SurfNodeFac*CMat(1)*Gamma2(SurfNode,3)
      s0(2,1) = s0(2,1) + SurfNodeFac*CMat(2)*Gamma2(rcmax+1-SurfNode,1)
      s0(2,2) = s0(2,2) + SurfNodeFac*CMat(2)*Gamma2(rcmax+1-SurfNode,2)
      s0(2,3) = s0(2,3) + SurfNodeFac*CMat(2)*Gamma2(rcmax+1-SurfNode,3)
    END DO
    s0(1,1) = s0(1,1) + REAL(NumOfPerpendNodes-1,r64)*DMat(1)
    s0(2,2) = s0(2,2) + REAL(NumOfPerpendNodes-1,r64)*DMat(2)
  END IF

  IF (NodeSource > 0) THEN
    s0(3,1) = Gamma2(NodeSource,1)
    s0(3,2) = Gamma2(NodeSource,2)
    s0(3,3) = Gamma2(NodeSource,3)
  END IF
  IF (NodeUserTemp > 0) THEN
    s0(4,1) = Gamma2(NodeUserTemp,1)
    s0(4,2) = Gamma2(NodeUserTemp,2)
    s0(4,3) = Gamma2(NodeUserTemp,3)
  END IF

          ! Check for and enforce symmetry in the cross term (Y)
  IF (ABS(s0(1,2)) /= ABS(s0(2,1))) THEN
    avg = (ABS(s0(1,2))+ABS(s0(2,1))) * 0.5D0
    s0(1,2) = avg*s0(1,2)/ABS(s0(1,2))
    s0(2,1) = avg*s0(2,1)/ABS(s0(2,1))
  END IF

          ! Compute S's and e's from 1 to n-1.  See equations (2.1.25) and
          ! (2.1.26) and Appendix C.
  inum = 1            ! Set history term counter
  CTFConvrg = .FALSE. ! Set the convergence logical to false

          ! The following DO WHILE loop calculates each successive set of time
          ! history terms until there are rcmax number of history terms or the
          ! latest flux history term is negligibly small compared to the first
          ! flux history term.
  DO WHILE ((.NOT.CTFConvrg) .AND. (inum < rcmax))   ! Begin CTF calculation loop ...

          ! Compute e(inum) based on Appendix C (Seem's dissertation). First,
          ! compute the new PhiR0 and its trace.

    trace = 0.0D0

    DO ir = 1,rcmax

      DO ic = 1,rcmax
        PhiR0(ir,ic) = 0.0D0
        DO is = 1,rcmax
          ! Make sure the next term won't cause an underflow.  If it will end up being
          ! so small as to go below TinyLimit, then ignore it since it won't add anything
          ! to PhiR0 anyway.
          IF (ABS(Rnew(is,ic)) > TinyLimit) THEN
            IF (ABS(AExp(ir,is)) > ABS(TinyLimit/Rnew(is,ic))) &
              PhiR0(ir,ic) = PhiR0(ir,ic)+AExp(ir,is)*Rnew(is,ic)
          END IF
        END DO
      END DO

      trace = trace + PhiR0(ir,ir)

    END DO

          ! Now calculate ej from the trace.  According to Appendix C:
          ! e(j) = -Trace[AExp*R(j-1)]/j

    e(inum) = -trace/REAL(inum,r64)

          ! Update Rold and compute Rnew.  Note:  PhiR0 = AExp*R(j-1) here.
          ! According to Appendix C:  R(j) = AExp*R(j-1) + e(j-1)

    DO ir = 1,rcmax
      DO ic = 1,rcmax
        Rold(ir,ic) = Rnew(ir,ic)
        Rnew(ir,ic) = PhiR0(ir,ic)
      END DO
      Rnew(ir,ir) = Rnew(ir,ir)+e(inum)
    END DO

          ! Compute S(inum) based on eq.(2.1.25) which states:
          ! S(j)  =  CMat*[R(j-1)*(Gamma1-Gamma2)+R(j)*Gamma2]
          !          + e(j)*DMat
    IF (SolutionDimensions == 1) THEN
      DO j = 1, 3
        DO is2 = 1,rcmax
          s(inum,1,j) = s(inum,1,j) + CMat(1)*( Rold(1,is2)*Gamma1(is2,j)  &
                                               +Rnew(1,is2)*Gamma2(is2,j) )
          s(inum,2,j) = s(inum,2,j) + CMat(2)*( Rold(rcmax,is2)*Gamma1(is2,j)  &
                                               +Rnew(rcmax,is2)*Gamma2(is2,j) )
          IF (NodeSource > 0) THEN
            s(inum,3,j) = s(inum,3,j) + ( Rold(NodeSource,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeSource,is2)*Gamma2(is2,j) )
          END IF
          IF (NodeUserTemp > 0) THEN
            s(inum,4,j) = s(inum,4,j) + ( Rold(NodeUserTemp,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeUserTemp,is2)*Gamma2(is2,j) )
          END IF

        END DO
        IF (j /= 3) s(inum,j,j) = s(inum,j,j) + e(inum)*DMat(j)
      END DO
    ELSE ! SolutionDimensions = 2
      DO j = 1, 3
        DO is2 = 1,rcmax
          DO SurfNode = 1, NumOfPerpendNodes
            IF ( (SurfNode == 1) .OR. (SurfNode == NumOfPerpendNodes) ) THEN
              SurfNodeFac = 0.5d0
            ELSE
              SurfNodeFac = 1.0d0
            END IF
            s(inum,1,j) = s(inum,1,j) + SurfNodeFac*CMat(1)*( Rold(SurfNode,is2)*Gamma1(is2,j)  &
                                                             +Rnew(SurfNode,is2)*Gamma2(is2,j) )
            s(inum,2,j) = s(inum,2,j) + SurfNodeFac*CMat(2)*( Rold(rcmax+1-SurfNode,is2)*Gamma1(is2,j)  &
                                                             +Rnew(rcmax+1-SurfNode,is2)*Gamma2(is2,j) )
          END DO
          IF (NodeSource > 0) THEN
            s(inum,3,j) = s(inum,3,j) + ( Rold(NodeSource,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeSource,is2)*Gamma2(is2,j) )
          END IF
          IF (NodeUserTemp > 0) THEN
            s(inum,4,j) = s(inum,4,j) + ( Rold(NodeUserTemp,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeUserTemp,is2)*Gamma2(is2,j) )
          END IF
        END DO
      END DO
      s(inum,1,1) = s(inum,1,1) + e(inum)*DMat(1)*REAL(NumOfPerpendNodes-1,r64)
      s(inum,2,2) = s(inum,2,2) + e(inum)*DMat(2)*REAL(NumOfPerpendNodes-1,r64)
    END IF

          ! Check for and enforce symmetry in the cross term (Y)
    IF (ABS(s(inum,1,2)) /= ABS(s(inum,2,1))) THEN
      avg = (ABS(s(inum,1,2))+ABS(s(inum,2,1))) * 0.5D0
      s(inum,1,2) = avg*s(inum,1,2)/ABS(s(inum,1,2))
      s(inum,2,1) = avg*s(inum,2,1)/ABS(s(inum,2,1))
    END IF

          ! Check for convergence of the CTFs.
    IF (e(1) == 0.0D0) THEN

      nrf = 1                   ! e(1) is zero, so there are no history terms.
      CTFConvrg = .TRUE.        ! CTF calculations have converged--set logical.

    ELSE
          ! e(1) is non-zero -- Calculate and compare the ratio of the flux
          ! terms to the convergence limit.
      rat = ABS(e(inum)/e(1))

      IF (rat < ConvrgLim) THEN

          ! If the ratio is less than the convergence limit, then any other
          ! terms would have a neglible impact on the CTF-based energy balances.
        nrf = inum
        CTFConvrg = .TRUE.    ! CTF calculations have converged--set logical.

      END IF
    END IF  ! ... end of convergence check block.

    inum = inum+1

  END DO    ! ... end of CTF calculation loop.
            ! Continue to the next coefficient if the solution has not converged
                    !
  IF (.NOT.CTFConvrg) THEN  ! Compute last e and S, if still unconverged.

          ! Compute e(inum) based on Appendix C (Seem's dissertation) or see
          ! equation above.  First compute the new PhiR0 and its trace.

    trace = 0.0D0

    DO ir = 1,rcmax
      DO is = 1,rcmax
        trace = trace+AExp(ir,is)*Rnew(is,ir)
      END DO
    END DO

    e(rcmax) = -trace/REAL(rcmax,r64)     ! Now calculate ej from the trace.

          ! Compute S(inum) based on eq.(2.1.25) which states:
          !   S(last) = CMat*R(last-1)*(Gamma1-Gamma2)+e(last)*DMat

    IF (SolutionDimensions == 1) THEN
      DO j = 1, 3
        DO is2 = 1,rcmax
          s(rcmax,1,j) = s(rcmax,1,j) + CMat(1)*Rnew(1,is2)*Gamma1(is2,j)
          s(rcmax,2,j) = s(rcmax,2,j) + CMat(2)*Rnew(rcmax,is2)*Gamma1(is2,j)
          IF (NodeSource > 0) THEN
            s(rcmax,3,j) = s(rcmax,3,j) + Rnew(NodeSource,is2)*Gamma1(is2,j)
          END IF
          IF (NodeUserTemp > 0) THEN
            s(rcmax,4,j) = s(rcmax,4,j) + Rnew(NodeUserTemp,is2)*Gamma1(is2,j)
          END IF
        END DO
      END DO
      s(rcmax,1,1) = s(rcmax,1,1)+e(rcmax)*DMat(1)
      s(rcmax,2,2) = s(rcmax,2,2)+e(rcmax)*DMat(2)
      nrf = rcmax
    ELSE ! SolutionDimensions = 2
      DO j = 1, 3
        DO is2 = 1,rcmax
          DO SurfNode = 1, NumOfPerpendNodes
            IF ( (SurfNode == 1) .OR. (SurfNode == NumOfPerpendNodes) ) THEN
              SurfNodeFac = 0.5d0
            ELSE
              SurfNodeFac = 1.0d0
            END IF
            s(rcmax,1,j) = s(rcmax,1,j) + SurfNodeFac*CMat(1)*Rnew(SurfNode,is2)*Gamma1(is2,j)
            s(rcmax,2,j) = s(rcmax,2,j) + SurfNodeFac*CMat(2)*Rnew(rcmax+1-SurfNode,is2)*Gamma1(is2,j)
          END DO
          IF (NodeSource > 0) THEN
            s(rcmax,3,j) = s(rcmax,3,j) + Rnew(NodeSource,is2)*Gamma1(is2,j)
          END IF
          IF (NodeUserTemp > 0) THEN
            s(rcmax,4,j) = s(rcmax,4,j) + Rnew(NodeUserTemp,is2)*Gamma1(is2,j)
          END IF
        END DO
      END DO
      s(rcmax,1,1) = s(rcmax,1,1) + e(rcmax)*DMat(1)*REAL(NumOfPerpendNodes-1,r64)
      s(rcmax,2,2) = s(rcmax,2,2) + e(rcmax)*DMat(2)*REAL(NumOfPerpendNodes-1,r64)
    END IF

          ! Check for and enforce symmetry in the cross term (Y)

    IF (ABS(s(rcmax,1,2)) /= ABS(s(rcmax,2,1))) THEN
      avg = (ABS(s(rcmax,1,2))+ABS(s(rcmax,2,1)))*0.5D0
      s(rcmax,1,2) = avg*s(rcmax,1,2)/ABS(s(rcmax,1,2))
      s(rcmax,2,1) = avg*s(rcmax,2,1)/ABS(s(rcmax,2,1))
    END IF

  END IF    ! ... end of IF block for calculation of last e and S.

  DEALLOCATE(PhiR0)
  DEALLOCATE(Rnew)
  DEALLOCATE(Rold)

  RETURN    ! The array e and the matrices s and s0 now contain the conduction
            ! transfer functions for this construction.

END SUBROUTINE CalculateCTFs
All Procedures
