CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > System Analysis

Crank-Nicolson FTN95 code

Register Blogs Community New Posts Updated Threads Search

Like Tree2Likes
  • 1 Post By tas38
  • 1 Post By beer

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 5, 2015, 15:15
Exclamation Crank-Nicolson FTN95 code
  #1
New Member
 
Join Date: Jan 2015
Posts: 24
Rep Power: 11
Athos1387 is on a distinguished road
Hi,
I must solve the question below using crank-nicolson method and Thomas algorithm by writing a code in fortran.



I've written a code for FTN95 as below. It works without a problem and gives me the answers, the problem is that the answers are wrong. I've solved it with FTCS method and analytically,and I know what the right answers are. I've checked the code 100 times and I don't understand what's wrong with it. can anyone pleeeeeas help me?!!!

************************************************** **
program HW36u
implicit none
real , dimension(19) :: u , uold , x , d , b , dprime
real , dimension(2:19) :: a
real , dimension(18) :: cprime , c
real , dimension( 19 , 19 ) :: K
real :: t , dt , dx , r , u0 , u20
integer :: j , i , n , m
write (*,*)'type your desirable time'
read (*,*) t
write (*,*)'type your desirable time step'
read (*,*) dt
n = t / dt
dx = 0.05
r = dt / (dx**2)
u0 = 0 !boundary condition

!Initial Condition
do i=1,19
x(i) = 0 + (i)*dx
if ( x(i) < 0.75 ) then
uold(i)= 0
else if ( x(i) > 0.75 ) then
uold(i) = 1
else
uold(i)= 0.5
end if
end do

!constructing the diagonals of the coeficents matrix A
do i=1,19
b(i) = (-2) * (r+1)
end do
do i=1,18
c(i) = r
end do
do i=2,19
a(i) = r
end do


!constructing the coeficents of knwon u, matrix K
do i=1,19
do j=1,19
if ( i == j ) then
K(i,j) = (2) * (r-1)
else if ( j-i == 1 .or. i-j == 1 ) then
K(i,j) = (-1) * r
else
K(i,j) = 0
end if
end do
end do

!computing the value of heat eq. for the desired time

!Thomas Algorithm
do m = 1,n

!constructing the knwon matrix d
d = matmul (K , uold)

!constructing c'
do i=1,18
if (i==1) then
cprime(i) = c(i) / b(i)
else
cprime(i) = c(i) / ( b(i) - cprime(i-1)*a(i) )
end if
end do

!constructing d'
do i=1,19
if (i==1) then
dprime(i) = d(i) / b(i)
else
dprime(i) = ( d(i) - dprime(i-1)*a(i) ) / ( b(i) - cprime(i-1)*a(i) )
end if
end do

!calculating u
do i=19,1,-1
if ( i == 19 ) then
u(i) = dprime(i)
else
u(i) = dprime(i) - ( cprime(i)*u(i+1) )
end if
end do

uold = u

end do

!Neuman boundary condition
u20 = ( (real(4)/3)*u(19) ) - ( (real(1)/3)*u(18) )

!output on screen
write (*,*) 'The data will be outputed to a file'

!output data to a file
open (1, file='hw3-6.txt' , status='new')
write(1,*) u0
do i=1,19
write(1,*) u(i)
end do
write(1,*) u20
close(1)

end program HW36u
************************************************** ***
Athos1387 is offline   Reply With Quote

Old   January 6, 2015, 20:26
Default
  #2
Senior Member
 
Troy Snyder
Join Date: Jul 2009
Location: Akron, OH
Posts: 220
Rep Power: 19
tas38 is on a distinguished road
I have a code that solves the 1d heat eqn. via Crank-Nicolson, I will compare tomorrow.
nbarman likes this.
tas38 is offline   Reply With Quote

Old   January 8, 2015, 11:52
Default
  #3
New Member
 
Join Date: Jan 2015
Posts: 24
Rep Power: 11
Athos1387 is on a distinguished road
thanks a lot!
Athos1387 is offline   Reply With Quote

Old   January 22, 2015, 11:14
Default
  #4
Member
 
Join Date: Dec 2012
Posts: 92
Rep Power: 14
beer is on a distinguished road
Do you still need help? I found the code quite interesting using the Thomas algorithm, so I've taken a look.

I think your coefficients are not quite right. Think about the equation again, separate the time steps n+1 and n and then build the coefficient matrix. At the Neuman boundary condition I used the shortcut that zero gradient means zero flux, so you don't need a dummy point outside the domain to reconstruct the gradient. Other than that your code seems all right and after fixing the coefficients the solution looks ok I guess.

Code:

program HW36u
implicit none
integer, parameter:: mag = 21
real , dimension(mag) :: u , uold , x , d , b , dprime
real , dimension(2:mag) :: a
real , dimension(mag-1) :: cprime , c
real , dimension( mag , mag ) :: K
real :: t , dt , dx , r , u0 , u20
integer :: j , i , n , m
write (*,*)'type your desirable time'
read (*,*) t
write (*,*)'type your desirable time step'
read (*,*) dt
n = t / dt
dx = 1. / (mag-1)
r = dt / (dx**2)
u0 = 0 !boundary condition

!Initial Condition
do i=1,mag
x(i) = 0 + (i-1)*dx
if ( x(i) < 0.75 ) then
u(i)= 0
else if ( x(i) > 0.75 ) then
u(i) = 1
else
u(i)= 0.5
end if
end do

!! Build coefficient matrix for steady state solution
!constructing the diagonals of the coeficents matrix A
do i=1,mag
b(i) = 2*r
end do
b(mag) = r
do i=1,mag-1
c(i) = -r
end do
c(1) = 0
do i=2,mag
a(i) = -r
end do
!constructing the coeficents of knwon u, matrix K
do i=1,mag
do j=1,mag
if ( i == j ) then
K(i,j) = -b(i)
else if ( j-i == 1) then
K(i,j) = -c(i)
else if ( i-j == 1 ) then
K(i,j) = -a(i)
else
K(i,j) = 0
end if
end do
end do

!! Now add transient term on both sides
do i=1,mag
b(i) = b(i) + 2
K(i,i) = K(i,i) + 2
end do

!computing the value of heat eq. for the desired time

!Thomas Algorithm
do m = 1,n

uold = u

!constructing the knwon matrix d
d = matmul (K , uold)

!constructing c'
do i=1,mag-1
if (i==1) then
cprime(i) = c(i) / b(i)
else
cprime(i) = c(i) / ( b(i) - cprime(i-1)*a(i) )
end if
end do

!constructing d'
do i=1,mag
if (i==1) then
dprime(i) = d(i) / b(i)
else
dprime(i) = ( d(i) - dprime(i-1)*a(i) ) / ( b(i) - cprime(i-1)*a(i) )
end if
end do

!calculating u
do i=mag,1,-1
if ( i == mag ) then
u(i) = dprime(i)
else
u(i) = dprime(i) - ( cprime(i)*u(i+1) )
end if
end do
end do


!output on screen
write (*,*) 'The data will be outputed to a file'

!output data to a file
open (1, file='hw3-6.txt')
do i=1,mag
write(1,100) x(i), ' ', u(i)
end do
100 format (F6.3,A3,F6.3)
close(1)

end program HW36u



Solution:
x u
0.000 0.000
0.050 0.021
0.100 0.042
0.150 0.063
0.200 0.086
0.250 0.110
0.300 0.135
0.350 0.162
0.400 0.190
0.450 0.219
0.500 0.248
0.550 0.279
0.600 0.308
0.650 0.337
0.700 0.365
0.750 0.390
0.800 0.412
0.850 0.431
0.900 0.446
0.950 0.456
1.000 0.461
Athos1387 likes this.

Last edited by beer; January 23, 2015 at 08:23.
beer is offline   Reply With Quote

Old   January 27, 2015, 13:37
Default
  #5
New Member
 
Join Date: Jan 2015
Posts: 24
Rep Power: 11
Athos1387 is on a distinguished road
Thanks a lot for the help! it was great!
Athos1387 is offline   Reply With Quote

Reply

Tags
crank nicolson, crank-nicolson, fortran code, ftn95, heat equation


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
The FOAM Documentation Project - SHUT-DOWN holger_marschall OpenFOAM 242 March 7, 2013 13:30
How to make code run in parallel? cwang5 OpenFOAM Programming & Development 1 May 30, 2011 05:47
Open Source Vs Commercial Software MechE OpenFOAM 28 May 16, 2011 12:02
Small 3-D code Zdravko Stojanovic Main CFD Forum 2 July 19, 2010 11:11
crank nicolson....emergency morteza08 Main CFD Forum 5 June 1, 2010 09:21


All times are GMT -4. The time now is 20:57.