MODULE Output USE IntrType USE Data CHARACTER*1 :: tab = ACHAR(9) ! easy way use tabs in output ! ! Subroutines associated with output ! CONTAINS ! SUBROUTINE edit ! ! Edit key data ! Calculate surface conditions, and energy balance ! Output is to UNIT 11 ! ! Programmed by John Mahaffy 1/01 ! IMPLICIT NONE INTEGER :: i,j REAL(sdk) :: qsurf, qvol, tmax, tmin LOGICAL :: isOpen ! ! OPEN the output file INQUIRE (UNIT=11,OPENED=isOpen) IF(.NOT.isopen) OPEN(11,FILE='trans.txt') ! ! Surface temperatures and heat loss at the surface ! txs(1:ny,1) = (3*hx*tf+9*cnd/dx*t(1,1:ny)-cnd/dx*t(2,1:ny))/(3*hx+8*cnd/dx) txs(1:ny,2) = (3*hx*tf+9*cnd/dx*t(nx,1:ny)-cnd/dx*t(nx-1,1:ny))/(3*hx+8*cnd/dx) tys(1:nx,1) = (3*hy*tf+9*cnd/dy*t(1:nx,1)-cnd/dy*t(1:nx,2))/(3*hy+8*cnd/dy) tys(1:nx,2) = (3*hy*tf+9*cnd/dy*t(1:nx,ny)-cnd/dy*t(1:nx,ny-1))/(3*hy+8*cnd/dy) qsurf = SUM(hx*(txs(1:ny,1)-tf)*dy) + & SUM(hx*(txs(1:ny,2)-tf)*dy) + & SUM(hy*(tys(1:nx,1)-tf)*dx) + & SUM(hy*(tys(1:nx,2)-tf)*dx) qvol = q*xsize*ysize ! Total heat generation ! ! Check for a match of total predicted surface heat flux with ! total heat produced in the region. ! IF (numStepsTried .EQ. 0) THEN WRITE(11,2000) qsurf, qvol 2000 FORMAT(/,5x, 'Surface Energy Loss', 3x,'Energy Generated', /, & 1x, 1P, 2E22.12) RETURN END IF WRITE(11,20) time, tab, MAXVAL(tss-t), tab, MAXVAL(t), tab,MINVAL(t), & tab, qsurf, tab, qvol ! ! Tabs are added to data to ease pickup by spreadsheets ! 20 FORMAT(1x, F7.3, 3(A, F18.12), 2(A, 1P, E22.12)) ! RETURN END SUBROUTINE edit SUBROUTINE TimeRichardson ! ! This applies Richardson based error analysis to a series of runs with ! time steps increasing from a base value with a fixed ratio. Note that ! the results in tLast(:,1) correspond to the run with the smallest time step. ! ! Programmed by John Mahaffy 2/05 IMPLICIT NONE INTEGER(sik) :: i, iloc(1), j REAL(sdk) :: avgErr, maxErr, minErr, avgOrder, maxOrder, minOrder ! ! Loop over all time step sets tried. If more than three step sizes ! were chosen, this gives the opportunity to check convergence of the ! estimated order of accuracy. ! IF (numStepsTried .LT. 3) RETURN ! WRITE (*,'(/a,i3,a,i3)')'Richardson Time Analysis for nx =', nx, ' and ny = ', ny ! DO i = 1, numStepsTried-2 timeOrder(:) = (tLast(:,i+2)-tlast(:,i+1))/(tLast(:,i+1)-tlast(:,i)) ! ! Evaluate Richardson based estimates of order of accuracy and error. ! For all points ! WHERE (timeOrder .GT. 0) timeOrder = log(timeOrder)/log(rTime) timeErr = (tLast(:,i+1)-tLast(:,i))/(rTime**timeOrder - 1.0_sdk) ELSEWHERE timeOrder = -1 timeErr = 0 END WHERE avgErr = SUM(ABS(timeErr),timeOrder.GT.0)/MAX(1,COUNT(timeOrder.GT.0)) iloc = MAXLOC(ABS(timeErr)) maxErr = timeErr(iloc(1)) iloc = MINLOC(ABS(timeErr),timeOrder.GT.0) minErr = timeErr(iloc(1)) avgOrder = SUM(timeOrder,timeOrder.GT.0)/MAX(1,COUNT(timeOrder.GT.0)) iloc = MAXLOC(timeOrder) maxOrder = timeOrder(iloc(1)) iloc = MINLOC(timeOrder,timeOrder.GT.0) minOrder = timeOrder(iloc(1)) ! WRITE (*, '(1p,3(a,e9.2))') ' Time Step Set = ', timeSteps(i),', ', & & timeSteps(i+1),',',timeSteps(i+2) WRITE (*,'(1p,3(a,e10.2))') ' For smallest step mean error = ', avgErr, & & ', Minimum error = ',minErr, ', Maximum error = ',maxErr ! WRITE (*,'(1p,3(a,e10.2))') ' Average order = ', avgOrder, & & ', Minimum Order = ',minOrder, ', Maximum Order = ',maxOrder END DO END SUBROUTINE TimeRichardson SUBROUTINE SpatialRichardson ! ! Calculate estimate of spatial order of accuracy and error estimates ! for each mesh. ! This requires evaluation on at least 3 spatial meshes, and monotonic ! error change. If error is not monotonic the value returned for p is 0 ! ! Programmed by John Mahaffy PSU 1/05 ! IMPLICIT NONE INTEGER(sik) :: i REAL(sdk) :: pC, pWx, pWy, errC, errWx, errWy, exactC, exactWx, exactWy REAL(sdk) :: error ! DO i = 3, numMeshTried ! ! Center point ! pC = (tCenter(i-2)-tCenter(i-1))/(tCenter(i-1)-tCenter(i)) IF (pC .GT. 0.0_sdk) THEN pC = log(pC)/log(rSpace) exactC = tCenter(i)-(tCenter(i-1)-tCenter(i))/(rSpace**pC-1) errC = tCenter(i)-exactC ELSE pC = 0.0_sdk exactC = 999 errC = 999 ENDIF ! ! Center of a surface facing in the negative x direction ! pWx = (txwc(i-2)-txwc(i-1))/(txwc(i-1)-txwc(i)) IF (pWx .GT. 0.0_sdk) THEN pWx = log(pWx)/log(rSpace) exactWx = txwc(i)-(txwc(i-1)-txwc(i))/(rSpace**pWx-1) errWx = txwc(i)-exactWx ELSE pWx = 0.0_sdk exactWx = 999 errWx = 999 ENDIF ! ! Center of a surface facing in the negative y direction ! pWy = (tywc(i-2)-tywc(i-1))/(tywc(i-1)-tywc(i)) IF (pWy .GT. 0.0_sdk) THEN pWy = log(pWy)/log(rSpace) exactWy = tywc(i)-(tywc(i-1)-tywc(i))/(rSpace**pWy-1) errWy = tywc(i)-exactWy ELSE pWy = 0.0_sdk exactWy = 999 errWy = 999 ENDIF IF (i .EQ. 3) THEN WRITE (*,*) 'Richardson Based Analysis for Mesh' WRITE (*,'(/,a,i3,a)')' For ',nCells(1),' mesh cells in each direction' error = tCenter(1)-exactC WRITE (*,'(1p,3(a,e15.7))')' Center temperature ', tCenter(1), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pC, ', "exact =', exactC error = txwc(1)-exactWx WRITE (*,'(1p,3(a,e15.7))')'X Surface temperature ', txwc(1), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pWx, ', "exact =', exactWx error = tywc(1)-exactWy WRITE (*,'(1p,3(a,e15.7))')'Y Surface temperature ', tywc(1), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pWy, ', "exact =', exactWy ! WRITE (*,'(/,a,i3,a)')' For ',nCells(2),' mesh cells in each direction' error = tCenter(2)-exactC WRITE (*,'(1p,3(a,e15.7))')' Center temperature ', tCenter(2), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pC, ', "exact =', exactC error = txwc(2)-exactWx WRITE (*,'(1p,3(a,e15.7))')'X Surface temperature ', txwc(2), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pWx, ', "exact =', exactWx error = tywc(2)-exactWy WRITE (*,'(1p,3(a,e15.7))')'Y Surface temperature ', tywc(2), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pWy, ', "exact =', exactWy WRITE (22,*) 'Richardson Based Analysis for Mesh' WRITE (22,'(/,a,i3,a)')' For ',nCells(1),' mesh cells in each direction' error = tCenter(1)-exactC WRITE (22,'(1p,3(a,e15.7))')' Center temperature ', tCenter(1), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pC, ', "exact =', exactC error = txwc(1)-exactWx WRITE (22,'(1p,3(a,e15.7))')'X Surface temperature ', txwc(1), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pWx, ', "exact =', exactWx error = tywc(1)-exactWy WRITE (22,'(1p,3(a,e15.7))')'Y Surface temperature ', tywc(1), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pWy, ', "exact =', exactWy ! WRITE (22,'(/,a,i3,a)')' For ',nCells(2),' mesh cells in each direction' error = tCenter(2)-exactC WRITE (22,'(1p,3(a,e15.7))')' Center temperature ', tCenter(2), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pC, ', "exact =', exactC error = txwc(2)-exactWx WRITE (22,'(1p,3(a,e15.7))')'X Surface temperature ', txwc(2), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pWx, ', "exact =', exactWx error = tywc(2)-exactWy WRITE (22,'(1p,3(a,e15.7))')'Y Surface temperature ', tywc(2), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pWy, ', "exact =', exactWy END IF ! WRITE (*,*) 'Richardson Based Analysis for Mesh' WRITE (*,'(/,a,i3,a)')' For ',nCells(i),' mesh cells in each direction' error = tCenter(i)-exactC WRITE (*,'(1p,3(a,e15.7))')' Center temperature ', tCenter(i), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pC, ', "exact =', exactC error = txwc(i)-exactWx WRITE (*,'(1p,3(a,e15.7))')'X Surface temperature ', txwc(i), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pWx, ', "exact =', exactWx error = tywc(i)-exactWy WRITE (*,'(1p,3(a,e15.7))')'Y Surface temperature ', tywc(i), ', error ', & error WRITE (*,'(1p,3(a,e15.7))')' Order ', pWy, ', "exact =', exactWy WRITE (22,*) 'Richardson Based Analysis for Mesh' WRITE (22,'(/,a,i3,a)')' For ',nCells(i),' mesh cells in each direction' error = tCenter(i)-exactC WRITE (22,'(1p,3(a,e15.7))')' Center temperature ', tCenter(i), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pC, ', "exact =', exactC error = txwc(i)-exactWx WRITE (22,'(1p,3(a,e15.7))')'X Surface temperature ', txwc(i), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pWx, ', "exact =', exactWx error = tywc(i)-exactWy WRITE (22,'(1p,3(a,e15.7))')'Y Surface temperature ', tywc(i), ', error ', & error WRITE (22,'(1p,3(a,e15.7))')' Order ', pWy, ', "exact =', exactWy ! END DO END SUBROUTINE SpatialRichardson END MODULE Output