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
Tuesday, October 6, 2009
Leap Frog
This is the same code but the module for updating the u in time is done by leap frog scheme.
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 for leap frog
subroutine update_u_leap(grid_2d,grid_2d_start,dt)
type(grid) :: grid_2d,grid_2d_start
real,intent(in) :: dt
real,allocatable,dimension(:,:) :: uold,u,unew
allocate(uold(0:grid_2d%ni+1,0:grid_2d%nj+1))
allocate(u(0:grid_2d%ni+1,0:grid_2d%nj+1))
allocate(unew(0:grid_2d%ni+1,0:grid_2d%nj+1))
uold(:,:) = grid_2d_start%u(:,:)
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)=uold(i,j)+u(i,j)*(-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_start%u(:,:)=u(:,:)
grid_2d%u(:,:)=unew(:,:)
deallocate(uold,u,unew)
End Subroutine update_u_leap
!
!
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 ut = 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, my_grid_start
!
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) call set_grid(my_grid_start,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)
! initial run
call update_u(my_grid,dt)
! subsequent runs
do n = 2,nt
call update_u_leap(my_grid,my_grid_start,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
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 for leap frog
subroutine update_u_leap(grid_2d,grid_2d_start,dt)
type(grid) :: grid_2d,grid_2d_start
real,intent(in) :: dt
real,allocatable,dimension(:,:) :: uold,u,unew
allocate(uold(0:grid_2d%ni+1,0:grid_2d%nj+1))
allocate(u(0:grid_2d%ni+1,0:grid_2d%nj+1))
allocate(unew(0:grid_2d%ni+1,0:grid_2d%nj+1))
uold(:,:) = grid_2d_start%u(:,:)
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)=uold(i,j)+u(i,j)*(-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_start%u(:,:)=u(:,:)
grid_2d%u(:,:)=unew(:,:)
deallocate(uold,u,unew)
End Subroutine update_u_leap
!
!
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 ut = 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, my_grid_start
!
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) call set_grid(my_grid_start,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)
! initial run
call update_u(my_grid,dt)
! subsequent runs
do n = 2,nt
call update_u_leap(my_grid,my_grid_start,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
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 ut = 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
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 ut = 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
Tuesday, July 14, 2009
Solving Shallow Water Wave Equations
//Dhruv Balwada(MATS 379)
//Program for solving Shallow water Equations
//Using the Method of Characterstics
//For constant ocean depth(still working on variable depth)
#include
#include
int main()
{
float dt, dx;
float g=9.81;
float w,a,m,r,l;
float u[1000][2],c[1000][2];
float ita, n[1000], h[1000];
int i,j,k;
FILE *fp;
fp=fopen("wave_eqn_data.dat","w");
// Characterstics
dt=1;
dx=100;
a=0.5;
w=4.398e-3;
m=0;
l=6.28/100000;
//Check CFL condtion
r=dt/dx;
if(r>0.5)
exit(0);
//Initialisation of the array
for(i=0; i<2; i++)
{
for(j=0; j<1000; j++)
{
u[j][i]=0;
c[j][i]=0;
}
}
//The values of h
for(i=0; i<1000; i++)
{
h[i]=m*(1000-i)*dx+200;
}
//The values of c
for(i=0; i<1000; i++)
{
c[i][0]=sqrt(g*h[i]);
}
// Big Time Loop
for(k=0; k<4000; k++)
{
//Boundary Condition Left
ita=a*sin(w*dt*k);
u[0][1]=ita*sqrt(g/h[0]);
c[0][1]=c[0][0]+0.5*(u[0][1]-u[0][0])-m*dt/2+(u[0][0]-c[0][0])*(dt/dx)*(0.5*(u[1][0]-u[0][0])-(c[1][0]-c[0][0]));
//Boundary Condition Right
u[999][1]=u[999][0]-(dt/(2*dx))*(3*u[999][0]-m*dt*k)*(u[999][0]-u[998][0])+m*dt;
c[999][1]=c[999][0]-(dt/dx)*(3*c[999][0]+m*k*dt)*(c[999][0]-c[998][0]);
//Chunk Calculation
for(i=1; i<999; i++)
{
u[i][1]=u[i][0]+(dt/dx)*((u[i][0]+c[i][0])*(0.5*u[i-1][0]-0.5*u[i][0]+c[i-1][0]-c[i][0])-(u[i][0]-c[i][0])*(0.5*u[i+1][0]-0.5*u[i][0]-c[i+1][0]+c[i][0])+m*dx);
c[i][1]=c[i][0]+(dt/(2*dx))*((u[i][0]+c[i][0])*(0.5*u[i-1][0]-0.5*u[i][0]+c[i-1][0]-c[i][0])+(u[i][0]-c[i][0])*(0.5*u[i+1][0]-0.5*u[i][0]-c[i+1][0]+c[i][0]));
}
//Find ita at all the points on the grid
for(j=0; j<1000; j++)
{
n[j]=u[j][1];
}
// Write to data file
if(k%100==0)
{
for(i=0; i<1000; i=i+20)
{
fprintf(fp,"%f",n[i]);
fprintf(fp,"\n");
}
fprintf(fp,"\n");
}
//Moving the Data in time
for(j=0; j<1000; j++)
{
u[j][0]=u[j][1];
c[j][0]=c[j][1];
}
}
fclose(fp);
}
Tuesday, June 9, 2009
My First Fortran Code : Temprature Conversion
!THIS CONVERTS celsius to farenheit
PROGRAM temp_conversion
IMPLICIT NONE
real :: far, cel
PRINT *, " Type The Value in Celsius "
READ *, cel
PRINT *, cel
far=1.8*cel+32.0
PRINT *, "Fareneit= ", far
END PROGRAM temp_conversion
PROGRAM temp_conversion
IMPLICIT NONE
real :: far, cel
PRINT *, " Type The Value in Celsius "
READ *, cel
PRINT *, cel
far=1.8*cel+32.0
PRINT *, "Fareneit= ", far
END PROGRAM temp_conversion
A Vibrating String Plucke and .8 distance from one end
This Image Shows the Vibration of a String Over time. It Goes up and comes back and so on..//Dhruv Balwada
#include"stdio.h"
#include"math.h"
int main()
{
float rho=0.01;
float T=40.0;
float x[101][3];
float c;
int i, k, l;
c=T/rho;
float c1=10*c;
float ratio=(c*c/(c1*c1));
FILE *fp;
fp=fopen("datanormal.dat", "w");
//Initialisation of the Array
for(i=0; i<3; i++)
{
for(k=0; k<101; k++)
x[k][i]=0;
}
//Initial Condition
for(i=0; i<81; i++){x[i][0]=0.00125*i;}
for(i=81;i<101;i++){x[i][0]=0.1-0.005*(i-80);}
//for(i=0; i<101; i++)
//{ l= sin(2*3.141*0.01*i);
//x[i][0]=0.001*l;
//}
//Loop to find the first Value
for(i=1; i<100; i++)
{
x[i][1]=x[i][0]+0.5*ratio*(x[i+1][0]+x[i-1][0]-2*x[i][0]);
}
//The Big time Loop Starts Here
for(k=1;k<1500; k++)
{
for(i=1; i<100; i++)
{
x[i][2]=2*x[i][1]-x[i][0]+ratio*(x[i+1][1]-2*x[i][1]+x[i-1][1]);
}
for(i=0; i<101; i++)
{
x[i][0]=x[i][1];
x[i][1]=x[i][2];
}
if((k%5)==0)
{
for(i=0; i<101; i++)
{
fprintf(fp,"%f",x[i][2]);
fprintf(fp,"\n");
}
fprintf(fp,"\n");
}
}//Big Time Loop Ends
}//Main Ends
Subscribe to:
Posts (Atom)
