Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations Mike Lewis on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

heat equation using crank-nicolsan scheme in fortran 1

Status
Not open for further replies.
watto8

watto8

Programmer
Feb 3, 2014
29
GB
Program crank_nicolson
call thomas(a,b,c,d,JI)
Real, allocatable :: x:)),u:)),a:)),b:)),c:)),d:))
Real:: m,dx,dt,Tmax
Integer:: j,NI,JI
Print*, 'Enter the total number of time steps'
read*, NI
Print*, 'Enter the final time'
read*, Tmax

dt=Tmax/NI !The size of timestep
Print*, 'This gives stepsize of time dt=',dt
dx= 0.1 !delta x =0.1
allocate (x(0:NI),u(0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))
m = dt/(2*dx**2)
JI= int(2/dx)-1
open(10,file='crank_nicolsan.m')

!Initial Condition
do j=1,JI-1
x(j)= -1+j*dx
u(j)=(1+x(j))*(1-x(j))**2 !x=-1+j*dx
end do

do n=1 ,NI-1 !BC
u(n)=0 !b.c. at j=0
u(JI)=0 !b.c. at j=JI

do j = 1, JI
a(j) = -m
b(j) = 1+2*m
c(j) = -m
d(j) = m*u(j+1)+(1-2*m)*u(j)+m*u(j-1)
end do
end do
do j=1, JI-1
u(j)=d(j)
end do

!Print out the Approximate solution in matlab file
write(10,*) 'ApproximateSolution =[',x(0),u(0)
do j =1, JI
u(j) = initialfun(j,dx)
end do
write(10,*)x(JI),(JI),']'

end Program crank_nicolson

subroutine thomas (a,b,c,d,JI)
implicit none
real, intent(inout) :: a(*),b(*),c(*),d(*)
integer, intent(in) :: JI
integer j

do j = 2,JI !combined decomposition and forward substitution
a(j) = a(j)/b(j-1)
b(j) = b(j)-a(j)*c(j-1)
d(j) = d(j)-a(j)*d(j-1)
end do

!back substitution
d(j) = d(j)/b(j)
do j = JI-1,1,-1
d(j) = (d(j)-c(j)*d(j+1))/b(j)
end do
return
end subroutine thomas



Error: Statement ordering error - REAL cannot appear after executable statements in this line Real, allocatable :: x:)),u:)),a:)),b:)),c:)),d:))
There might be some other errors and missing code too. Can help me and check my code please. Thanks in advance














 
Sort by date Sort by votes
the error message tells it all, let alone that you should know better than that, if you have typed this much code yourself...just my suspicion.

so, the second line in the program is out of place...where does it really go?
 
Upvote 0 Downvote
The error statement is pretty clear. These lines:

Code: 
Real, allocatable :: x(:),u(:),a(:),b(:),c(:),d(:)
Real:: m,dx,dt,Tmax
Integer:: j,NI,JI

have to appear before this line:

Code: 
call thomas(a,b,c,d,JI)

I'm not sure why you're calling the Thomas solver before you even assemble a, b, c, and d. You need to define JI before you allocate the arrays. Lastly, the function initialfun(.,.) is not defined anywhere.
 
Upvote 0 Downvote
Code: 
Program crank_nicolson
implicit none
Real, allocatable :: x(:),u(:,:),a(:),b(:),c(:),d(:)
 Real::  m,dx,dt,Tmax
Integer:: n,j,NI,JI
Print*, 'Enter the total number of time steps'
read*, NI
Print*, 'Enter the final time'
read*, Tmax
dt=Tmax/NI !The size of timestep
Print*, 'This gives stepsize of time dt=',dt 
dx= 0.1 !delta x =0.1
allocate (x(0:JI),u(0:NI,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))
m = dt/(2*dx**2)
JI= 20
open(10,file='crank_nicolsan.m')

!Initial Condition
do j=1,JI-1
  x(j)= -1+j*dx
u(0,j)=(1+x(j))*(1-x(j))**2  !x=-1+j*dx
end do
x(0)= -1
x(JI)= 1

!Boundary condition
do n=0,NI
u(n,0)=0      !b.c. at j=0
u(n,JI)=0     !b.c. at j=JI
end do

do n =0 , NI !loop
do j = 1, JI  ! a,b,c,d are the coefficients of C-N scheme 
          a(j) = -m
        b(j) = 1+2*m
        c(j) = -m
       d(j) = m*u(n,j+1)+(1-2*m)*u(n,j)+m*u(n,j-1)    
end do
end do
   
do j=1, JI-1
     do n=1, NI-1
       u(n,j)=d(j)
       end do
       end do
      
call thomas(a,b,c,d,JI)

!Print out the Approximate solution in matlab file
write(10,*)  'ApproximateSolution =[',x(0),u(0,0)
do j =1, JI
write(10,*) x(j),u(n,j)
 end do
write(10,*)x(JI),u(NI,JI),']'
 
end Program crank_nicolson
 
subroutine thomas (a,b,c,d,JI)
implicit none
real, intent(inout) :: a(*),b(*),c(*),d(*)
   integer, intent(in) :: JI 
    integer j
  
  do j = 2,JI                   !combined decomposition and forward substitution
      a(j) = a(j)/b(j-1)
      b(j) = b(j)-a(j)*c(j-1)
      d(j) = d(j)-a(j)*d(j-1)
   end do     
     
!back substitution
   d(j) = d(j)/b(j)
   do j = JI-1,1,-1
      d(j) = (d(j)-c(j)*d(j+1))/b(j)
   end do
  return
   end subroutine thomas

After doing some changes, im still getting the another Runtime-Error: Undefined variable, element or arrays of functions in this line: allocate (x(0:JI),u(0:NI,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI)).
There were no complier errors this time.
Please help me. Thank You
 
Upvote 0 Downvote
After moving JI=20 up, it works but still getting another run-time errors in this line: d(j) = m*u(n,j+1)+(1-2.0*m)*u(n,j)+m*u(n,j-1).
Run-time errors saying: Undefined variablse, arrays element or functions reference to that line.
There were no complier errors.

Code: 
 Program crank_nicolson
implicit none
Real, allocatable :: x(:),u(:,:),a(:),b(:),c(:),d(:)
Real:: m,dx,dt,Tmax
Integer:: n,j,NI,JI
Print*, 'Enter the total number of time steps'
read*, NI
Print*, 'Enter the final time'
read*, Tmax
dt=Tmax/NI !The size of timestep
Print*, 'This gives stepsize of time dt=',dt 
dx= 0.1 !delta x =0.1
JI= 20
allocate (x(0:JI),u(0:NI,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))
m = dt/(2*dx**2)
open(10,file='crank_nicolsan.m')
!Initial Condition
do j=1,JI-1
x(j)= -1+j*dx
u(0,j)=(1+x(j))*(1-x(j))**2 !x=-1+j*dx
end do
!x(0)= -1
!x(JI)= 1
!Boundary condition
do n=0,NI
u(n,0)=0 !b.c. at j=0
u(n,JI)=0 !b.c. at j=JI
end do
do n =0 , NI !loop
do j = 1, JI-1 ! a,b,c,d are the coefficients of C-N scheme 
a(j) = -m
b(j) = 1+2.0*m
c(j) = -m
d(j) = m*u(n,j+1)+(1-2.0*m)*u(n,j)+m*u(n,j-1) 
end do
end do

do j=1, JI-1
do n=1, NI-1
u(n,j)=d(j)
end do
end do

call thomas(a,b,c,d,JI)
!Print out the Approximate solution in matlab file
write(10,*) 'ApproximateSolution =[',x(0),u(0,0)
do j =1, JI
write(10,*) x(j),u(n,j)
end do
write(10,*)x(JI),u(NI,JI),']'

end Program crank_nicolson

subroutine thomas (a,b,c,d,JI)
implicit none
real, intent(inout) :: a(*),b(*),c(*),d(*)
integer, intent(in) :: JI 
integer j

do j = 2,JI !combined decomposition and forward substitution
a(j) = a(j)/b(j-1)
b(j) = b(j)-a(j)*c(j-1)
d(j) = d(j)-a(j)*d(j-1)
end do 

!back substitution
d(j) = d(j)/b(j)
do j = JI-1,1,-1
d(j) = (d(j)-c(j)*d(j+1))/b(j)
end do
return
end subroutine thomas
 
Upvote 0 Downvote
I am also getting some others run-time errors in these lines:

call thomas(a,b,c,d,JI)

a(j) = a(j)/b(j-1)

Run-time errors saying: Undefined variablse, arrays element or functions reference to those lines.

There were no complier errors.
Please help me.
 
Upvote 0 Downvote
Just after a rapid glance, I already detect two troubles in your THOMAS routine :

- it expects arrays from 1 to JI but you pass arrays from 1 to JI+1 (because declared 0:JI in the calling program). I suggest to call that routine with JI+1 as last argument.

- after the comment !back substitution, the next instruction is incorrect : the variable j is "out of bound" because that instruction is out of any loop over j. Put JI instead => d(JI) = d(JI)/b(JI)

Last advice : run your program under a debugger please, execute it instruction by instruction and check at each step whether this program really computes what you want!

François Jacq
 
Upvote 0 Downvote
Program crank_nicolson
implicit none
Real, allocatable :: x:)),u:),:),a:)),b:)),c:)),d:))
Real:: m,dx,dt,Tmax
Integer:: n,j,NI,JI
Print*, 'Enter the total number of time steps'
read*, NI
Print*, 'Enter the final time'
read*, Tmax
dt=Tmax/NI !The size of timestep
Print*, 'This gives stepsize of time dt=',dt
dx= 0.1 !delta x =0.1
JI= 20
allocate (x(0:JI),u(0:NI,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))
m = dt/(2*dx**2)
open(10,file='crank_nicolsan.m')

!Initial Condition
do j=1,JI-1
x(j)= -1+j*dx
u(0,j)=(1+x(j))*(1-x(j))**2 !x=-1+j*dx
end do
x(0)= -1
x(JI)= 1

!Boundary condition
do n=0,NI
u(n,0)=0 !b.c. at j=0
u(n,JI)=0 !b.c. at j=JI
end do

do n =0 , NI !loop
do j = 1, JI-1 ! a,b,c are the coefficients of C-N scheme and d
a(j) = -m
b(j) = 1+2.0*m
c(j) = -m
d(j) = m*u(n,j+1)+(1-2.0*m)*u(n,j)+m*u(n,j-1)
call thomas(a,b,c,d,JI)
u(n,j)=d(j)
end do
end do

!Print out the Approximate solution in matlab file
write(10,*) 'ApproximateSolution =[',x(0),u(0,0)
do j =1, JI
write(10,*) x(j),u(n,j)
end do
write(10,*)x(JI),u(NI,JI),']'
close(10)
end Program crank_nicolson

subroutine thomas (a,b,c,d,JI)
implicit none
real, intent(inout) :: a(*),b(*),c(*),d(*)
integer, intent(in) :: JI
integer j

do j = 2,JI-1 !combined decomposition and forward substitution
a(j) = a(j)/b(j-1)
b(j) = b(j)-a(j)*c(j-1)
d(j) = d(j)-a(j)*d(j-1)
end do

!back substitution
d(JI) = d(JI)/b(JI)
do j = JI-1,1,-1
d(j) = (d(j)-c(j)*d(j+1))/b(j)
end do
return
end subroutine thomas



(I DID SOME CHANGES WITH THE HELP OF FRANCOIS JACQ. Where do i need to write JI+1 instead JI. IS IT IN allocate (x(0:JI),u(0:NI,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI)) OF u,a,b,c and d)? THANKS FRANCOIS.

BUT,I am still getting run-time errors in these lines:

THOMAS - in file crank_nicolson.f95 at line a(j) = a(j)/b(j-1)
main - in file crank_nicolson.f95 at line call thomas(a,b,c,d,JI)

Run-time errors saying: Undefined variablse, array element or functions result reference to those lines.

There were no complier errors.
















 
Upvote 0 Downvote
I have added JI Plus 1 in this line call thomas(a,b,c,d,JI+1).

BUT,I am still getting run-time errors in these lines:

THOMAS - in file crank_nicolson.f95 at line a(j) = a(j)/b(j-1)
main - in file crank_nicolson.f95 at line call thomas(a,b,c,d,JI)

Run-time errors saying: Undefined variablse, array element or functions result reference to those lines.

And if i add JI PLUS 1 in this line subroutine thomas (a,b,c,d,JI) and also in this line call thomas(a,b,c,d,JI+1).

i got error 471 - Invalid item(s) in argument list in this line: subroutine thomas (a,b,c,d,JI)
 
Upvote 0 Downvote
Give us the exact and full error messages without paraphrasing.

If you have a runtime error with instructions like a(j) = a(j)/b(j-1), this is probably because b(j-1) is equal to zero. But "runtime error" is a too general message. The exact message might be SEGFAULT, SIGABRT or directly "division by zero".

"call thomas(a,b,c,d,JI+1)" is obviously what I suggested...

But you should use a debugging tool. Do you have one ?

In reading your main program, I noticed that the call to the subroutine thomas is not placed correctly. You wrote :

Code: 
do n =0 , NI !loop
  do j = 1, JI-1 ! a,b,c are the coefficients of C-N scheme and d 
    a(j) = -m
    b(j) = 1+2.0*m
    c(j) = -m
    d(j) = m*u(n,j+1)+(1-2.0*m)*u(n,j)+m*u(n,j-1) 
    call thomas(a,b,c,d,JI) 
    u(n,j)=d(j)
  end do
end do

I am almost sure that it should be something like :
Code: 
do n =0 , NI !loop
  do j = 1, JI-1 ! a,b,c are the coefficients of C-N scheme and d 
    a(j) = -m
    b(j) = 1+2.0*m
    c(j) = -m
    d(j) = m*u(n,j+1)+(1-2.0*m)*u(n,j)+m*u(n,j-1) 
  end do
  call thomas(a,b,c,d,JI) 
  u(n+1,:)=d ! u(n,:) is doubtful 
end do

This is why, for instance, b is not initialize correctly and you probably divide by zero or use an invalid number.



François Jacq
 
Upvote 0 Downvote
Code: 
 Program crank_nicolson
implicit none
Real, allocatable :: x(:),u(:,:),a(:),b(:),c(:),d(:)
 Real::  m,dx,dt,Tmax
Integer:: n,j,NI,JI
Print*, 'Enter the total number of time steps'
read*, NI
Print*, 'Enter the final time'
read*, Tmax
dt=Tmax/NI !The size of timestep
Print*, 'This gives stepsize of time dt=',dt 
dx= 0.1 !delta x =0.1
JI= 20
allocate (x(0:JI),u(0:NI+1,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))
m = dt/(2*dx**2)
open(10,file='crank_nicolsan.m')

!Initial Condition
do j=1,JI-1
  x(j)= -1+j*dx
u(0,j)=(1+x(j))*(1-x(j))**2  !x=-1+j*dx
end do
x(0)= -1
x(JI)= 1

!Boundary condition
do n=0,NI
u(n,0)=0      !b.c. at j=0
u(n,JI)=0     !b.c. at j=JI
end do
do n =0 , NI !loop
do j = 1, JI-1  ! a,b,c are the coefficients of C-N scheme and d is the right part
          a(j) = -m
        b(j) = 1+2.0*m
        c(j) = -m
       d(j) = m*u(n,j+1)+(1-2.0*m)*u(n,j)+m*u(n,j-1)  
       end do  
   call thomas_algorithm(a,b,c,d,JI+1)     
u(n+1,j)=d(j)
end do

!Print out the Approximate solution in matlab file
write(10,*)  'ApproximateSolution =[',x(0),u(0,0)
do j =1, JI
write(10,*) x(j),u(n,j)
 end do
write(10,*)x(JI),u(NI,JI),']'
 close(10)
 
contains

subroutine thomas_algorithm (a,b,c,d,JI)
implicit none
real, intent(inout) :: a(*),b(*),c(*),d(*)
   integer, intent(in) :: JI 
   integer j
  
  do j = 2,JI                 !combined decomposition and forward substitution
    a(j) = a(j)/b(j-1)
      b(j) = b(j)-a(j)*c(j-1)
      d(j) = d(j)-a(j)*d(j-1)
   end do 
 
!back substitution
   d(j) = d(j)/b(j)
   do j = JI-1,1,-1
      d(j) = (d(j)-c(j)*d(j+1))/b(j)
   end do
end subroutine thomas_algorithm
end Program crank_nicolson





 
Upvote 0 Downvote
There were no complier errors. So after i complie and run the program, i got run-time error.
Run-time Error:
Error 112, Reference to undefined variable,array element or function result[/UNDEF]

CRANK_NICOLSON ~THOMAS_ALGORITHM- in file crank_nicolsonf.95 at line 59[+0081] which is a(j) = a(j)/b(j-1)

main - in file crank_nicolson.f95 at line 38 [+0f99] which is
call thomas_algorithm(a,b,c,d,JI+1)


I am trying to solve the 1d heat equation using crank-nicolson scheme. And for that i have used the thomas algorithm in the subroutine. Can you please check my subroutine too, did i missed some codes?? Im trying to connect the subroutine into main program and link it together to generate the value of u(n+1,j) and open the output and graphics into the matlab files.
 
Upvote 0 Downvote
I think the problem is that

[ul][li] in the main program, you are setting JI to 20[/li]
[li]then, you allocate your arrays to have indices from 0 to JI, or 0 to 20[/li]
[li]then, you call your subroutine and pass JI+1 or 21 as the last argument, and, so[/li]
[li]the variable JI local to the subroutine (don't let the name fool you, this is a different variable that JI in the main program))has the value of 21
[li] which on the last iteration of the loop gives the value of 21 to j,
[li] which is then used in a(j)...but a(.) is only defined up to 20![/li][/ul]
 
Upvote 0 Downvote
Salgerman Can you please clarify me by changing the code that i needs to, so can i understand it. Thank you
 
Upvote 0 Downvote
in the subroutine thomas_algorithm (a,b,c,d,JI), i still have JI which is 20.
 
Upvote 0 Downvote
Nope, I bet you don't have JI=20 inside the subroutine...add a write statement and print the value of JI from within your subroutine, you will see.

Other than that, I don't know what you are modeling; I am just paying attention to the Fortran and helping you debug...so, I don't know exactly how to change your code to make it correct...all I know is why you are getting an error message.
 
Upvote 0 Downvote
subroutine thomas_algorithm (a,b,c,d,JI)
! a - sub-diagonal (means it is the diagonal below the main diagonal)
! b - the main diagonal
! c - sup-diagonal (means it is the diagonal above the main diagonal)
! d - right part
implicit none
real, intent(inout),dimension(JI) :: a(0:JI),b(0:JI),c(0:JI),d(0:JI)
integer, intent(in) :: JI
integer j

do j = 2,JI !combined decomposition and forward substitution
a(j) = a(j)/b(j-1)
b(j) = b(j)-a(j)*c(j-1)
d(j) = d(j)-a(j)*d(j-1)
end do

!back substitution
d(JI) = d(JI)/b(JI)
do j = JI-1,1,-1
d(j) = (d(j)-c(j)*d(j+1))/b(j)
end do
end subroutine thomas_algorithm
end Program crank_nicolson_one

I CHANGED THE SUBROUTINE BUT STILL GETTING THE ERRORS.
SINCE I CAN NOT ADD JI+1 IN THE subroutine thomas_algorithm (a,b,c,d,JI+1).
 
Upvote 0 Downvote
Try the following program. I have corrected many things :
- arrays have not to start from 0.
- I kept u from 0 to ji+1 to manage the boundary conditions
- I have deleted the storage of u versus time because the results are printed on a file
- I have inserted the routine THOMAS within the program to check the arguments

I finally understood your program in interpreting that (lambda*S)/(rho*cp) is equal to 1.

Code: 
program crank_nicolson

implicit none
real, allocatable :: x(:),u(:),a(:),b(:),c(:),d(:)
real:: m,dx,dt,tmax
integer:: n,j,ni,ji

print*, 'enter the total number of time steps'
read*, ni
print*, 'enter the final time'
read*, tmax
dt=tmax/ni !the size of timestep
print*, 'this gives stepsize of time dt=',dt 
dx= 0.1 !delta x =0.1
ji= 20
m = dt/(2*dx**2)

allocate (x(ji),u(0:ji+1),a(ji),b(ji),c(ji),d(ji))
open(10,file='crank_nicolsan.m')

!initial condition

do j=1,ji
  x(j)= -1+j*dx
  u(j)=(1+x(j))*(1-x(j))**2 !x=-1+j*dx
end do


!boundary condition

u(0)=0
u(ji+1)=0 

do n =1 , ni !time loop
  do j = 1, ji ! a,b,c,d are the coefficients of c-n scheme 
    a(j) = -m
    b(j) = 1+2.0*m
    c(j) = -m
    d(j) = m*u(j+1)+(1-2.0*m)*u(j)+m*u(j-1) 
  end do
  call thomas(a,b,c,d)
  do j=1,ji
    u(j)=d(j)
  end do
  write(10,*) 'step ',n
  do j =1, ji
    write(10,*) x(j),u(j)
  enddo
end do

contains

subroutine thomas (a,b,c,d)

  real, intent(inout) :: a(:),b(:),c(:),d(:)
  integer j

  do j = 2,ji !combined decomposition and forward substitution
    a(j) = a(j)/b(j-1)
    b(j) = b(j)-a(j)*c(j-1)
    d(j) = d(j)-a(j)*d(j-1)
  end do 

  !back substitution
  d(ji) = d(ji)/b(ji)
  do j = ji-1,1,-1
    d(j) = (d(j)-c(j)*d(j+1))/b(j)
  end do

end subroutine thomas 

end program crank_nicolson

François Jacq
 
Upvote 0 Downvote
Status
Not open for further replies.

Similar threads

Vahid9
  • Locked
  • Question
Printing in a fortran MPI/OpenMP environment
Replies
2
Views
380
Vahid9
Vahid9
R
  • Locked
  • Question
How do I print a real number using GNU gfortran
Replies
3
Views
464
rdflg
R
turbo_arte
  • Locked
  • Question
Multall Fortran 77 help 1
Replies
11
Views
349
mikrom
mikrom
ksqr
  • Locked
  • Question
Looking for a Fortran 77 or newer compiler fo PC 1
2
Replies
28
Views
1K
xwb
xwb
Deniall
  • Locked
  • Question
Determining a running Fortran program's "residence directory" 2
2
Replies
25
Views
822
mikrom
mikrom

Part and Inventory Search

Sponsor

Back
Top