1. 2.
    0
    c---|
    c..A POINT/LINE ITERATIVE SOLVER FOR ELLIPTIC PDEs |
    c Course: AE305 | |
    c---|

    program ELLIPTIC

    parameter (imx=101, jmx=201)
    common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
    common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
    common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
    real error, error_tol, error0
    integer i,j,k
    c..Read the input data, generate the grid data and initialize the solution
    call INIT
    write(*,*) 'enter error_tol'
    read(*,*) error_tol
    c..Start the iterative solution loop

    k = 0
    m=alog10(error/error0)
    DO WHILE (m.gt.error_tol)
    k = k + 1
    c..Point iterative solutions
    call ITERATIVE
    c..Line iterative solution
    c call ADI(... )

    c..Update phik array and evaluate residual
    rl2=0.
    do j = 1,jmax
    do i = 1,imax
    rl2 = rl2 + (phikp(i,j) - phik(i,j))**2
    phik(i,j) = phikp(i,j)
    enddo
    enddo
    rl2 = sqrt(rl2/(imax*jmax))
    print*, 'Residual@k =',k,rl2
    c... error definition
    error=0
    do j=2,jmax
    do i=2,imax
    error=error+(phik(i,j)-phik(i,j))**2
    enddo
    enddo
    error=error/float(imax*jmax)
    error=sqrt(error)
    c...
    call sSTRESS
    c..Output intermediate solutions
    if( mod(k,iostep).eq.0 .and. k.ne.kmx) call IO(k)
    ENDDO
    call IO(k)

    stop
    end

    c---
    subroutine INIT
    parameter (imx=101, jmx=201)
    common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
    common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
    common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
    xl=0.05
    yl=0.10
    Torque = 100
    c..set the second moment of area, J and the rhs
    sma = xl*yl*(xl**2+yl**2)/3.
    rhs = -2*Torque/sma

    c..Read input parameters like imax, jmax
    print*, 'Enter Imax and Jmax :'
    read(*,*) imax, jmax
    dx = xl/float(imax-1)
    dy = yl/float(jmax-1)
    beta2 = (dx/dy)**2

    c..Grid generation
    do j=1,jmax
    do i=1,imax
    x(i,j)= dx*float(i-1)
    y(i,j)= dy*float(j-1)
    enddo
    enddo

    c..Initialize the phi arrays
    do j = 1,jmax
    do i = 1,imax
    phik(i,j) = 0.
    tauzx(i,j)=0.
    tauzy(i,j)=0.
    enddo
    enddo
    error0=1.
    error=error0
    return
    end

    c---
    subroutine ITERATIVE
    parameter (imx=101, jmx=201)
    common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
    common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)

    c..Solve for phi^k+1
    do j = 2,jmax-1
    do i = 2,imax-1
    c..Implement, Jacobi, Gauss-Seidel and SOR methods
    phikp(i,j) =(1./2.*(1+beta2))*(-rhs+phik(i-1,j)+phik(i+1,j)+beta2*
    &(phik(i,j-1)+phik(i,j+1)))
    enddo
    enddo
    return
    end
    c---
    subroutine sSTRESS
    parameter (imx=101, jmx=201)
    common /var/beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
    common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
    common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
    do j=2,jmax-1
    do i=2,imax-1
    phik(i,j)=phikp(i,j)
    tauzx(i,j)=-(phik(i+1,j)-phik(i-1,j))/2.*dx
    tauzy(i,j)=(phik(i,j+1)-phik(i,j-1))/2.*dy
    enddo
    enddo
    return
    end

    c---
    subroutine IO(k)
    c..Output solution in tecplot format
    parameter (imx=101, jmx=201)
    common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
    common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
    common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
    character fname*32,string*8,ext*5

    write(string,'(f8.5)') float(k)/100000
    read(string,'(3x,a5)') ext
    fname = 'phi-'//ext//'.tec'
    open(1,file=fname, form='formatted')
    write(1,*) ' variables="x","y","phi","tauzx","tauzy" '
    write(1,*) ' zone i=',imax, ', j=',jmax
    do j = 1,jmax
    do i = 1,imax
    write(1,*) x(i,j),y(i,j),phik(i,j),tauzx(i,j),tauzy(i,j)
    enddo
    enddo
    close(1)

    return
    end

    c---
    c subroutine THOMAS(AA,BB,CC,FF,n)
    ---
    ! Solution of a tridiagonal system of n equations of the form
    ! aa(i)*x(i-1) + bb(i)*x(i) + cc(i)*x(i+1) = ff(k) for i=1,n
    ! the solution x(k) is stored in ff(k)
    ! aa(1) and cc(n) is not used.
    ! aa,bb,cc,ff are arrays to be filled by the caller program
    !______________________________________________________________
    c real aa(1), bb(1), cc(1), ff(1)
    c
    c n2 = 2
    c n1n = n+1
    c..forward sweep
    c
    c bb(1) = 1./bb(1)
    c aa(1) = ff(1)*bb(1)
    c do k = n2,n
    c k1 = k-1
    c cc(k1) = cc(k1)*bb(k1)
    c bb(k) = 1./(bb(k)-aa(k)*cc(k1))
    c aa(k) = (ff(k)-aa(k)*aa(k1))*bb(k)
    c enddo
    c
    c..back substitution
    c
    c ff(n) = aa(n)
    c do k1 = n2,n
    c k = n1n-k1
    c ff(k) = aa(k)-cc(k)*ff(k+1)
    c enddo
    c
    c return
    c end
    Tümünü Göster
    ···
  2. 1.
    0
    reserved
    ···