Problem with periodic boundary condition:

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

Problem with periodic boundary condition:

Postby debsankar » Sat Apr 27, 2013 9:47 am

I was trying to solve Diffution Eq in finite difference implicit method with periodic boundary condition.... Thats actually solving Ax=B in a loop with substituting B=x after each calculation.
I used DGESV but it does not give me right answer.Here the matrix A(NxN) is a sparse matrix which has A(1,N) and A(N,1) non-zero due to PBC.
Here is my code:

[code][/code]
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! CODE TO SOLVE DIFFUSION EQUATION USING IMPLICIT METHOD
! Df= diffusion co-eff, delt and delx are grid spacing of time
! and space in discretized equations.
! the files:'initial.txt' contains the initial distribution and
! 'diffusion.txt' contains time evoluted distribution.


program diffusion
implicit none
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 x(1:100),A(1:100,1:100),p(1:100)
real*8 Df,delx,delt,time
integer i,j,k,L,T,pivot(100),OK,nhrs

Df=1.d0
L=100
delx=1.d0/L
delt=4.d0*delx*delx/Df
time=0.005d0
T=int(time/delt)
nhrs=1


open(unit=01,file='initial.txt')
open(unit=02,file='diffusion.txt')
!cccccccccccccccccccccccccccccccccccccccccccccccc
! Initialisation
!cccccccccccccccccccccccccccccccccccccccccccccccc

do i=1,L
x(i)=0.d0
if(i*delx.ge.0.4) then
if(i*delx.le.0.6) then
x(i)=1.d0
endif
endif
enddo
! x(L/2)=1.d0

do i=1,L
do j=1,L
A(i,j)=0.d0
enddo
enddo

do i=1,L
write(01,*) i*delx,x(i)
enddo

!cccccccccccccccccccccccccccccccccccccccccccccccc
! Constructing the sparse matrix
!cccccccccccccccccccccccccccccccccccccccccccccccc

do i=2,L-1
A(i,i)=1.d0+((2.d0*Df*delt)/(delx*delx))
A(i,i-1)=-(Df*delt)/(delx*delx)
A(i,i+1)=-(Df*delt)/(delx*delx)
! if(i.eq.1) A(i,L)=A(i,i-1)
! if(i.eq.L) A(i,1)=A(i,i+1)
! write(02,*) A(i,i-1),A(i,i),A(i,i+1)

enddo
A(1,1)=1.d0+((2.d0*Df*delt)/(delx*delx))
A(1,2)=-(Df*delt)/(delx*delx)
A(1,L)=-(Df*delt)/(delx*delx)

A(L,L)=1.d0+((2.d0*Df*delt)/(delx*delx))
A(L,1)=-(Df*delt)/(delx*delx)
A(L,L-1)=-(Df*delt)/(delx*delx)

! write(*,*) A(1,L),A(L,1)



!cccccccccccccccccccccccccccccccccccccccccccccccc
! Time Loop starts:
!cccccccccccccccccccccccccccccccccccccccccccccccc

do k=1,T
do i=1,L
p(i)=x(i)
! write(*,*) p(i)
enddo


call DGESV( L, nhrs, A, L, pivot, p, L, OK )


do i=1,L
x(i)=p(i)
!write(*,*) p(i)
enddo


enddo

write(*,*) OK

do i=1,L
write(02,*) i*delx,x(i)
enddo



stop
end program diffusion
************************************************************************************
Its really frustating if I can not figure out where it goes wrong.... Please Help.
Regards,
Deb.
debsankar
 
Posts: 1
Joined: Sat Apr 27, 2013 9:24 am

Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 3 guests

cron