Crank Nicholsan

This is the code for doing it using Crank Nicholsan.
This is an implicit scheme and is always stable.
Also earlier we were doing using explicit schemes, this one is implicit and unconditionally stable.

Module Grids
Type grid
integer :: ni
real :: dx
real,dimension(:),allocatable :: xt0
real,dimension(:),allocatable :: u
End type grid

Contains
Subroutine set_grid (grid_1d,ni)
type(grid) :: grid_1d
integer, intent(in) :: ni
integer :: i
real :: pi
!
pi = acos(-1.0)
grid_1d%ni = ni
grid_1d%dx = 1.0/(ni+1)
allocate (grid_1d%xt0(ni))
allocate (grid_1d%u(0:ni+1))
do i = 1,ni
grid_1d%xt0(i) = i*grid_1d%dx
enddo
grid_1d%u(:) = 0.0
! set initial conditions as well
do j = 1,nj
do i = 1,ni
grid_1d%u(i) = sin(pi*grid_1d%xt0(i))
enddo
enddo

End Subroutine set_grid
End Module Grids


Module netcdf_setup
use Grids


Contains
subroutine nc_setup (grid_1d,ncid,nc_filename)
integer,intent(inout) ::ncid
type (grid) :: grid_1d
integer :: ier
character(len=20),intent(inout)::nc_filename
character(len=20) ::fname
integer :: xdimid,timedimid
integer :: xid
integer :: timeid,tempid
include 'netcdf.inc'
!
fname = trim(nc_filename)

! For netcdf output
! define dimensions
ier=nf_create(fname,nf_clobber,ncid)
ier=NF_DEF_DIM(ncid,'x',grid_1d%ni,xdimid)
ier=NF_DEF_DIM(ncid,'time',nf_unlimited,timedimid)
ier=NF_DEF_VAR(ncid,'x',NF_REAL,1,(/xdimid/),xid)
ier=NF_DEF_VAR(ncid,'time',NF_REAL,1,(/timedimid/),timeid)
! define dependent variables
ier=NF_DEF_VAR(ncid,'Temp',NF_REAL,2,(/xdimid,timedimid/),tempid)
! attributes
ier=NF_PUT_ATT_TEXT(ncid,xid,'long_name',1,'x')
ier=NF_PUT_ATT_TEXT(ncid,xid,'units',1,'m')
ier=NF_PUT_ATT_TEXT(ncid,timeid,'long_name',10,'Julian Day')
ier=NF_PUT_ATT_TEXT(ncid,timeid,'units',31,'days since 1994-01-01 00:00:00')
ier=NF_PUT_ATT_TEXT(ncid,tempid,'units',6,'Kelvin')
!
ier=NF_ENDDEF(ncid)
!
ier=NF_PUT_VARA_REAL(ncid,xid,(/1/),(/grid_1d%ni/),grid_1d%xt0)
ier=nf_close(ncid)
!
End subroutine nc_setup
End Module netcdf_setup

! Module for time stepping using the Crank Nicholsan Scheme
Module update_temp
use Grids

Contains
subroutine update_u(grid_1d,dt)
type(grid) :: grid_1d
real,intent(in) :: dt
real :: D
real,allocatable,dimension(:) :: u,unew,C,e,f
real,allocatable,dimension(:,:) :: A,B

allocate(u(0:grid_1d%ni+1))
allocate(unew(0:grid_1d%ni+1))

! assignment of local variables
u(:) = grid_1d%u(:)
dx = grid_1d%dx
D=dt/(dx**2)

! We will solve a system of the form AU^n+1=Bu^n
! Defining A and B
allocate(A(1:grid_1d%ni,1:grid_1d%ni))
allocate(B(1:grid_1d%ni,1:grid_1d%ni))

do i = 2,grid_1d%ni-1
A(i,i) = 1+D
B(i,i) = 1-D
A(i,i-1) = -D/2
A(i,i+1) = -D/2
B(i,i-1) = D/2
B(i,i+1) = D/2
enddo
A(1,1) = 1+D
A(1,2) = -D/2
B(1,1) = 1-D
B(1,2) = D/2
A(grid_1d%ni,grid_1d%ni) = 1+D
A(grid_1d%ni,grid_1d%ni-1) = -D/2
B(grid_1d%ni,grid_1d%ni) = 1-D
B(grid_1d%ni,grid_1d%ni-1) = D/2

! Calculating Bu^n
allocate(C(1:grid_1d%ni))
do i = 2,grid_1d%ni-1
C(i) = B(i,i)*u(i)+B(i,i-1)*u(i-1)+B(i,i+1)*u(i+1)
enddo
C(1) = B(1,1)*u(1)+B(1,2)*u(2)
C(grid_1d%ni) = B(grid_1d%ni,grid_1d%ni)*u(grid_1d%ni)+B(grid_1d%ni,grid_1d%ni-1)*u(grid_1d%ni-1)

! Calculating U^n+1 using the TDMA(thomas algo)
allocate(e(1:grid_1d%ni))
allocate(f(1:grid_1d%ni))

do i = 2,grid_1d%ni
m = A(i,i-1)/A(i-1,i-1)
e(i) = A(i,i)-m*A(i-1,i)
f(i) = C(i)-m*C(i-1)
enddo

unew(grid_1d%ni) = f(grid_1d%ni)/e(grid_1d%ni)

do i = grid_1d%ni-1,1
unew(i)=(f(i)-A(i,i+1)*unew(i+1))/e(i)
enddo

grid_1d%u(:) = unew(:)
deallocate(u,unew,A,B,C,e,f)
End Subroutine update_u
End Module update_temp

! Dhruv ... Program to solve the diffusion equations using the Crank Nicholsan Method
! with zero bc's and sin(pi*x) initial condition
Program heat_eqn
use Grids
use netcdf_setup
use update_temp
implicit none

integer :: ni, nt
real, allocatable,dimension(:) :: usave
real :: dt,dx
include 'netcdf.inc'
integer :: ncid,ier
character (len=20) ::nc_filename

integer timeid,tempid,ltime
!
type(grid) :: my_grid
!
integer :: i,n
!
ni=16
nt=500
dt=8e-4
nc_filename = 'heat_eqn_1D.nc'
!
write(6,'(a,2i5,es12.4)' )'ni,nt,dt',ni,nt,dt
! Make some checks
if ( ni <= 0 .OR. dt <= 0 .OR. nt <= 0) then write(6,*)'=>Error one of ni,nt,dt is <= 0'
endif
! Call gird and initialisation
call set_grid(my_grid,ni)

! netcdf setup
call nc_setup (my_grid,ncid,nc_filename)
!Loop over time and write output

ltime=0
allocate (usave(ni))
ier=nf_open(trim(nc_filename),NF_WRITE,ncid)
ier=nf_inq_varid(ncid,'Temp',tempid)
ier=nf_inq_varid(ncid,'time',timeid)
do n = 1,nt
call update_u(my_grid,dt)
if (mod(n,10) .eq.0) then
ltime=ltime+1
usave(1:my_grid%ni) = my_grid%u(1:my_grid%ni)
ier=NF_PUT_VARA_REAL(ncid,tempid,(/1,ltime/),(/my_grid%ni,1/),usave)
ier=NF_PUT_VARA_REAL(ncid,timeid,(/ltime/),(/1/),n*1.0)

endif
enddo ! for n loop
ier=nf_close(ncid)
END PROGRAM heat_eqn