MODULE SetEqn USE IntrType USE Data CONTAINS SUBROUTINE SetCoeffs ! ! Set coefficients in the finite volume conduction ! ! Programmed by John Mahaffy 1/01 ! IMPLICIT NONE REAL(sdk) :: cols(5), coeffs(5), cxm, cxp, cyp, cym, fth INTEGER(sik) :: i,j, ieq, ieqxp, ieqxm, ieqym, ieqyp, ii, ic, ieqt INTEGER(sik) :: nCoef ! ! Count the number of non-zero coefficients in each row ! as(1:n)%nCoef = 5 as(1:nx)%nCoef = as(1:nx)%nCoef -1 ! Bottom Boundary as(n-nx+1:n)%nCoef = as(n-nx+1:n)%nCoef -1 ! Top Boundary as(1:n:nx)%nCoef = as(1:n:nx)%nCoef -1 ! Left Boundary as(ny:n:nx)%nCoef = as(ny:n:nx)%nCoef - 1 ! Right Boundary ! DO i = 1,n ii = as(i)%nCoef ALLOCATE(as(i)%a(ii), as(i)%iCol(ii), as(i)%iDim(ii)) ENDDO as(1:n)%nCoef = 0 ieq = 0 ! First calculate matrix contribution from heat flux ! across interior faces DO j = 1,ny ! y Coefficients IF(j.GT.1) THEN cym = d/dy**2 ELSE cym = 0.0_sdk ENDIF IF(j.LT.ny) THEN cyp = d/dy**2 ELSE cyp = 0.0_sdk ENDIF DO i=1,nx ieq = ieq+1 ieqxm = ieq - 1 ieqxp = ieq + 1 ieqym = ieq - nx ieqyp = ieq + nx IF(i.GT.1) THEN cxm = d/dx**2 ! x Coefficients ELSE cxm = 0.0_sdk ENDIF IF(i.LT.nx) THEN cxp = d/dx**2 ELSE cxp = 0.0_sdk ENDIF ! List of columns with nonzero coefficients cols=(/ieqym, ieqxm, ieq, ieqxp, ieqyp/) ! List of those nonzero coefficients coeffs = (/-cym, -cxm, cym+cxm+cxp+cyp+dton, -cxp, -cyp/) ! DO ii = 1,5 ic = cols(ii) IF(ic.lt.1.OR.ic.GT.n) CYCLE a(ieq,ic) = coeffs(ii) ENDDO ENDDO ENDDO ! ! Adjustments for the y face boundary conditions ! DO ieq=1,nx a(ieq,ieq) = a(ieq,ieq)+2*d*hy/(dy*(hy*dy+2*cnd)) !Bottom ieqt = ieq + n - nx a(ieqt,ieqt) = a(ieqt,ieqt)+2*d*hy/(dy*(hy*dy+2*cnd)) !Top ENDDO ! ! Adjustments for the x face boundary conditions ! DO ieq=1,n,nx a(ieq,ieq) = a(ieq,ieq)+2*d*hx/(dx*(hx*dx+2*cnd)) !Left ieqt = ieq + nx - 1 a(ieqt,ieqt) = a(ieqt,ieqt)+2*d*hx/(dx*(hx*dx+2*cnd)) !Right ENDDO ! ! Load the sparse version of the matrix ! DO ieq = 1,n DO ic = 1,n IF(a(ieq,ic).EQ.0.0_sdk) CYCLE nCoef = as(ieq)%nCoef + 1 as(ieq)%nCoef = nCoef as(ieq)%a(nCoef) = a(ieq,ic) as(ieq)%iCol(nCoef) = ic IF(ieq.EQ.ic) THEN as(ieq)%iDim(nCoef) = 0 as(ieq)%iDiag = nCoef ELSE IF(ic.LE.ieq-nx.OR.ic.GE.ieq+nx) THEN as(ieq)%iDim(nCoef) = 2 ELSE as(ieq)%iDim(nCoef) = 1 ENDIF ENDDO ENDDO END SUBROUTINE SetCoeffs SUBROUTINE SetRHS ! Set the right hand side of the equation. Note that when dton is 1.0 ! you get the right hand side of a transient equation. When dton is 0.0 ! you get the right hand side of a steady state equation. ! ! Programmed by John Mahaffy 1/01 IMPLICIT NONE INTEGER(sik) :: j, i, k, kp ! ! Contributions common to all volumes ! DO j=1,ny DO i=1,nx k = (j-1)*nx + i b(k) = p + dton*t(i,j) ! (dton allows use of coding for transients ENDDO ENDDO ! ! Adjustments for the y face boundary conditions ! DO k=1,nx b(k) = b(k)+2*d*hy*tf/(dy*(hy*dy+2*cnd)) kp = k + n - nx b(kp) = b(kp)+2*d*hy*tf/(dy*(hy*dy+2*cnd)) ENDDO ! Adjustments for the x face boundary conditions DO k=1,n,nx b(k) = b(k)+2*d*hx*tf/(dx*(hx*dx+2*cnd)) kp = k + nx - 1 b(kp) = b(kp)+2*d*hx*tf/(dx*(hx*dx+2*cnd)) ENDDO RETURN END SUBROUTINE SetRHS END MODULE