!=======================================================
! 1-dimensional Diffusion Equation with simple
! Dirichlet boundary conditions u(0,t)=u(1,t)=0
! 0<= x <= 1 and 0<= t <= tf
! 
! We set initial condition u(x,t=0) that satisfies 
! the given boundary conditions. 
! Nx is the number of points in spatial lattice:
! x = 0 + (j-1)*dx, j=1,...,Nx and dx = (1-0)/(Nx-1)
! Nt is the number of points in temporal lattice:
! t = 0 + (j-1)*dt, j=1,...,Nt and dt = (tf-0)/(Nt-1)
!
! u(x,0) = sin(pi*x) tested against analytical solution
! u(x,t) = sin(pi*x)*exp(-pi*pi*t)
!
!=======================================================
program diffusion_1d
 implicit none 
 integer , parameter   :: dp = 8
 real(dp), allocatable :: u(:), du(:)
 real(dp)              :: t,x,dx,dt,tf,courant
 integer               :: Nx,Nt,i,j
 real(dp), parameter   :: ZERO=0.0_dp, ONE=1.0_dp, HALF=0.5_dp,TWO=2.0_dp
 real(dp), parameter   :: PI=atan2(ZERO,-ONE)
! --- Input:
 print *, '# Enter: Nx, Nt, tf:'
 read  *,           Nx, Nt, tf
 if(Nx <=  3) stop 'Nx <= 3'
 if(Nt <=  2) stop 'Nt <= 2'

 allocate(u(Nx),du(Nx))

! --- Initialize:
 dx      = ONE/(Nx-1)
 dt      = tf /(Nt-1)
 courant = dt /dx**2
 print * ,'# 1d Diffusion Equation: 0<=x<=1, 0<=t<=tf'
 print * ,'# dx= ',dx,' dt= ',dt,' tf= ', tf
 print * ,'# Nx= ',Nx,' Nt= ',Nt
 print * ,'# Courant Number= ',courant
 if(courant > HALF) print *,'# WARNING: courant > 0.5'
 open(unit=11,file='d.dat') ! data file
! --- Initial condition at t=0 ------------------------------
!u(x,0) = sin( pi x)
 do i  = 1, Nx
  x    = (i-1)*dx
  u(i) = sin(PI*x)
 end do
 u(1)  = ZERO
 u(Nx) = ZERO
 do i= 1,Nx
  x    = (i-1)*dx
  write(11,*) ZERO, x, u(i)
 end do
 write (11,*)' '
! ----------------------------------------------------------
! --- Calculate time evolution:
 do   j  = 2, Nt
  t      = (j-1)*dt
! ----- second derivative:
  do  i  = 2, Nx-1
   du(i) = courant*(u(i+1)-TWO*u(i)+u(i-1))
  end do
! ----- update:
  do i   = 2, Nx-1
   u(i)  = u(i) + du(i)
  end do
  do i   = 1, Nx
   x     = (i-1)*dx
   write(11,*) t, x, u(i)
  end do
  write (11,*)' '
       
 end do ! do j=2,Nt

 close(11)
end program diffusion_1d
!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2022)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------