CFD Review  
Serving the CFD Community with News, Articles, and Discussion
 
CFD Review

User Preferences
Site Sponsorship
Headline Feeds
Mobile Edition
Privacy Policy
Terms of Service
twitter

Submit a CFD Story

Site Sponsors
Siemens PLM Software
Pointwise: Reliable CFD meshing
Software Cradle

Tell a Friend
Help this site to grow by sending a friend an invitation to visit this site.

CFD News by Email
Did you know that you can get today's CFD Review headlines mailed to your inbox? Just log in and select Email Headlines Each Night on your User Preferences page.

 
Help With Sample CFD Program
Posted Tue August 08, 2006 @04:58PM
Print version Email story Tweet story
Help Desk Caprisun writes "Can someone explain this program to me please? I'm trying to understand it for a project and I'm kinda new to fortran. My lecturer gave it to me as an example :S..."

Sponsor CFD Review

PROGRAM QUESTION3
INTEGER I,NI,J,NJ,N,NP,NMAX,II,IL,IR
PARAMETER (NI=10,NJ=10)
REAL DX,DY,XTOTAL,YTOTAL,D,TINY, dt, tfin, t
PARAMETER(NP=100,N=100,NMAX=100)
DIMENSION A(NP,NP),INDEX(N),B(N)
REAL XNI,YNJ

        XNI=10.0 
        YNJ=10.0
        
        XTOTAL=1.
        YTOTAL=1.
              
        DX=XTOTAL/(XNI-1)
        DY=YTOTAL/(YNJ-1)
		

!INITIALISE MATRIX AND RHS

DO I=1,N
        DO J=1,N
           A(I,J)=0.
        ENDDO
        B(I)=0.
ENDDO  
!COMPLETED

!FILL IN MATRIX WITH COEFFICIENTS FOR INTERNAL NODES

DO I=2,NI-1
        DO J=2,NJ-1
                IL=NI*(J-2)+I
                II=NI*(J-1)+I
                IR=NI*(J)+I
                A(II,IL)=1./(DY*DY)
                A(II,II-1)=1./(DX*DX)
                A(II,II)=-2.*(1./(DX*DX)+1./(DY*DY))
                A(II,II+1)=1./(DX*DX)
                A(II,IR)= 1./(DY*DY)
        ENDDO
ENDDO
!COMPLETED  

 ! LEFT SIDE
   I=1
DO J=1,NJ-1
        II=NI*(J-1)+I
        A(II,II)=1.
        B(II)=300
ENDDO
! CORNER NODE 
        I=1
        J=NJ
        II=NI*(J-1)+I
        A(II,II)=1.
        B(II)=350
!COMPLETED

!TOP SIDE
J=NJ
DO I=2,NI
        II=NI*(J-1)+I
        A(II,II)=1.
        B(II)=400
ENDDO
!COMPLETED

!BOTTOM SIDE                       
   J=2
DO I=2,NI-1
        II=I
        IR=NI*(J-1)+I
        A(II,II)=1./(DY)
        A(II,IR)=-1./(DY)
ENDDO
!COMPLETED

!CORNER NODE
A(10,10)=1. ! modify here
A(10,19)=-1. ! modify here
!COMPLETED

!RIGHT SIDE                
DO J=2,NJ-1
        II=NI*(J-1)+NI
        IL=NI*(J-1)+(NI-1)
        A(II,IL)=-1./(DX)
        A(II,II)=1./(DX)
ENDDO
!COMPLETED 

                   
  !INITIALISE MATRIX AND RHS
       
        CALL LUDCMP(A,N,NP,INDEX,D)
        
        CALL LUBKSB(A,N,NP,INDEX,B)    
    
        DO I=1,N
                WRITE(*,*) B(I)
        ENDDO
        
         OPEN(1,FILE='results.txt') 
         DO I=1,NI
                 DO J=1,NJ
               II=NI*(J-1)+I
                 WRITE(1,*)B(II)
                 ENDDO
         ENDDO
        
END PROGRAM               
      
!------------------------------------------------

SUBROUTINE LUDCMP( A, N, NP, INDEX, D )

  INTEGER N, NP, INDEX(N), NMAX
  REAL D,A(NP,NP),TINY
  PARAMETER (NMAX=100,TINY=1.0E-20)
  INTEGER I,IMAX,J,K
  REAL  AAMAX,DUM,SUM,VV(NMAX)

   D=1.

   DO I=1,N
     AAMAX=0.
     DO J=1,N

     IF (ABS(A(I,J)).GT.AAMAX)AAMAX=ABS(A(I,J))

     ENDDO

     IF (AAMAX.EQ.0.) PAUSE 'SINGULAR MATRIX ON LUDCOMP'

     VV(I)=1./AAMAX

  ENDDO

  DO J=1,N
    DO I=1,J-1
       SUM=A(I,J)

       DO K=1,I-1
         SUM=SUM-A(I,K)*A(K,J)
       ENDDO

      A(I,J)=SUM
    ENDDO

    AAMAX=0.

    DO I=J,N
      SUM=A(I,J)

      DO K=1,J-1
        SUM=SUM-A(I,K)*A(K,J)
      ENDDO

      A(I,J)=SUM

      DUM=VV(I)*ABS(SUM)

      IF(DUM.GE.AAMAX)THEN
            IMAX=I
            AAMAX=DUM
      ENDIF
    ENDDO

    IF (J.NE.IMAX) THEN
      DO K=1,N
              DUM=A(IMAX,K)
              A(IMAX,K)=A(J,K)
              A(J,K)=DUM
      ENDDO
      D=-D
      VV(IMAX)=VV(J)
    ENDIF

    INDEX(J)=IMAX

    IF(A(J,J).EQ.0.)A(J,J)=TINY

    IF(J.NE.N)THEN
            DUM=1./A(J,J)
            DO I=J+1,N
            A(I,J)=A(I,J)*DUM
            ENDDO
    ENDIF

  ENDDO

  RETURN
END
!------------------------------------------------
!------------------------------------------------
SUBROUTINE LUBKSB( A, N, NP, INDEX, B )

  INTEGER N,NP,INDEX(N)
  REAL  A(NP,NP),B(N)
  INTEGER I,II,J,LL
  REAL  SUM

  II=0

  DO I=1,N
  LL=INDEX(I)
  SUM=B(LL)
  B(LL)=B(I)

      IF(II.NE.0.)THEN
        DO J=II,I-1
          SUM=SUM-A(I,J)*B(J)
        ENDDO
      ELSEIF(SUM.NE.0.) THEN
         II=I
      END IF

  B(I)=SUM
  ENDDO

  DO I = N, 1, -1
  SUM=B(I)

  DO J=I+1,N
          SUM=SUM-A(I,J)*B(J)
  END DO
  B(I)=SUM/A(I,I)
  ENDDO

  RETURN
END
!------------------------------------------------

[ Post Comment ]

NASA Ames Purchases CEI's EnSight Gold | Tecplot "Enjoy the View Day"  >

 

 
CFD Review Login
User name:

Password:

Create an Account

Related Links
  • More on Help Desk
  • Also by nwyman
  • This discussion has been archived. No new comments can be posted.

    Your lover will never wish to leave you. All content except comments
    ©2017, Viable Computing.

    [ home | submit story | search | polls | faq | preferences | privacy | terms of service | rss  ]