|
[Sponsors] |
January 5, 2015, 15:15 |
Crank-Nicolson FTN95 code
|
#1 |
New Member
Join Date: Jan 2015
Posts: 24
Rep Power: 11 |
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 ************************************************** *** |
|
January 6, 2015, 20:26 |
|
#2 |
Senior Member
Troy Snyder
Join Date: Jul 2009
Location: Akron, OH
Posts: 220
Rep Power: 19 |
I have a code that solves the 1d heat eqn. via Crank-Nicolson, I will compare tomorrow.
|
|
January 8, 2015, 11:52 |
|
#3 |
New Member
Join Date: Jan 2015
Posts: 24
Rep Power: 11 |
thanks a lot!
|
|
January 22, 2015, 11:14 |
|
#4 |
Member
Join Date: Dec 2012
Posts: 92
Rep Power: 14 |
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 Last edited by beer; January 23, 2015 at 08:23. |
|
January 27, 2015, 13:37 |
|
#5 |
New Member
Join Date: Jan 2015
Posts: 24
Rep Power: 11 |
Thanks a lot for the help! it was great!
|
|
Tags |
crank nicolson, crank-nicolson, fortran code, ftn95, heat equation |
|
|
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 |