Workshop : Numerical Methods, Fortran, Debugging and Parallelisation

Hi,

Been attending a few lectures by Dr. Swathi. Here the example problems were to solve the diffusion equation, using FORTRAN. I will deposit my code here, if someone wants to use it.

This code solves the code for explicit form.

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

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

End Subroutine set_grid
End Module Grids
!
Module netcdf_setup
use Grids


Contains
subroutine nc_setup (grid_2d,ncid,nc_filename)
integer,intent(inout) ::ncid
type (grid) :: grid_2d
integer :: ier
character(len=20),intent(inout)::nc_filename
character(len=20) ::fname
integer :: xdimid,ydimid,timedimid
integer :: xid,yid
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_2d%ni,xdimid)
ier=NF_DEF_DIM(ncid,'y',grid_2d%nj,ydimid)
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,'y',NF_REAL,1,(/ydimid/),yid)
ier=NF_DEF_VAR(ncid,'time',NF_REAL,1,(/timedimid/),timeid)
! define dependent variables
ier=NF_DEF_VAR(ncid,'Temp',NF_REAL,3,(/xdimid,ydimid,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,yid,'long_name',1,'y')
ier=NF_PUT_ATT_TEXT(ncid,yid,'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_2d%ni/),grid_2d%xt0)
ier=NF_PUT_VARA_REAL(ncid,yid,(/1/),(/grid_2d%nj/),grid_2d%yt0)
ier=nf_close(ncid)
!
End subroutine nc_setup
End Module netcdf_setup

Module update_temp
use Grids
Contains
subroutine update_u(grid_2d,dt)
type(grid) :: grid_2d
real,intent(in) ::dt
real,allocatable,dimension(:,:) :: u,unew
!
allocate(u(0:grid_2d%ni+1,0:grid_2d%nj+1))
allocate(unew(0:grid_2d%ni+1,0:grid_2d%nj+1))

u(:,:) = grid_2d%u(:,:)
dx = grid_2d%dx
dy = grid_2d%dy
do j = 1,grid_2d%nj
do i = 1,grid_2d%ni
unew(i,j) = u(i,j)*(1-2*dt*(dx**2 + dy**2)/(dx**2*dy**2)) + &
dt*((u(i-1,j) + u(i+1,j))/dx**2 + (u(i,j-1) + u(i,j+1))/dy**2)
enddo
enddo
grid_2d%u(:,:) = unew(:,:)
deallocate(u,unew)
End Subroutine update_u
End Module update_temp


!
! Swathi..Test program to solve u<sub>t = delsq(u) on a unit square
! with zero bc's and sin(pi*x)*sin(pi*y) initial condition
PROGRAM heat_eqn
use Grids
use netcdf_setup
use update_temp
implicit none
integer :: ni, nj, nt
real, allocatable,dimension(:,:) :: usave
real :: dt,dx,dy,dt_cfl
include 'netcdf.inc'
integer :: ncid,ier
character (len=20) ::nc_filename

integer timeid,tempid,ltime
!
type(grid) :: my_grid
!
integer :: i,j,n
!
ni=16
nj=16
nt=500
dt=8e-4
nc_filename = 'heat_eqn_2D.nc'
!
write(6,'(a,3i5,es12.4)' )'nj,nj,nt,dt',nj,nj,nt,dt
! Make some checks
if ( nj <= 0 .OR. nj <= 0 .OR. dt <= 0 .OR. nt <= 0) then write(6,*)'=>Error one of ni,nj,nt,dt is <= 0' endif ! Call gird and initialisation call set_grid(my_grid,ni,nj) ! Do a cfl check dx = my_grid%dx dy = my_grid%dy dt_cfl = 0.5*dx**2*dy**2/(dx**2+dy**2) write(6,*)'dt_cfl=',dt_cfl if ( dt > dt_cfl) then
write(6,*)'CFL condition violated: stopping'
stop
endif
! netcdf setup
call nc_setup (my_grid,ncid,nc_filename)
!Loop over time and write output

ltime=0
allocate (usave(ni,nj))
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,1:my_grid%nj) = my_grid%u(1:my_grid%ni,1:my_grid%nj)
ier=NF_PUT_VARA_REAL(ncid,tempid,(/1,1,ltime/),(/my_grid%ni,my_grid%nj,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</sub>