Program HeatedPlateForwardEuler

  ! 2D Heat Equation
  ! Finite difference
  ! Forward Euler

  implicit none

  integer                 :: ii, jj, kk           ! counters
  real, parameter         :: C  = 1.0             ! thermal diffusivity (m^2/s) (gold is 1.27E-4)
  real, parameter         :: x  = 1.0             ! side of square (m)
  real, parameter         :: t  = 1               ! simulation time (s)
  real, parameter         :: dx = 0.02            ! x grid spacing (0.01=1cm)
  real, parameter         :: dy = dx              ! square grid cells
  real, parameter         :: dt = (dx**2 * dx**2)/(2*C*(dx**2 + dx**2))  ! observe the CFL condition for 2d Heat
  real, parameter         :: mu = C * dt / dx**2  ! weighting
  real(kind=8), parameter :: epsilon = 1.0E-7     ! convergence threshold
  logical                 :: converged = .false.  ! has the simulation converged?
  integer                 :: lastiter = 0         ! last iteration if converged 
  integer, parameter      :: N = (x/dx)*(x/dx)    ! no. of cells in domain
  integer, parameter      :: M = (x/dx)+2         ! including boundary conditions
  integer                 :: start_count, end_count, count_rate, count_max ! for timing
  real(kind=8), target, allocatable, dimension(:,:) :: U1, U2          ! array storage
  real(kind=8), pointer, dimension(:,:)             :: UK, UKP1, temp  ! for fast storage re-use

  !! allocate memory
  allocate(U1(M,M));
  allocate(U2(M,M));

  !! initialise..
  ! pointers
  UK => U1
  UKP1 => U2
  ! values
  UK = 0.0d+0
  UKP1 = 0.0d+0
  ! and boundary conditions
  ! top row
  UK(M,:) = 100.0d+0
  UKP1(M,:) = 100.0d+0
  ! first & last columns
  UK(:,1) = 100.0d+0
  UK(:,M) = 100.0d+0
  UKP1(:,1) = 100.0d+0
  UKP1(:,M) = 100.0d+0

  call system_clock(start_count,count_rate,count_max) ! start timing
  do kk=1,int(t/dt)

     ! calcs
     do jj=2,M-1
        do ii=2,M-1
           UKP1(ii,jj) = UK(ii,jj) & 
                + mu*(UK(ii-1,jj) - 2*UK(ii,jj) + UK(ii+1,jj)) &
                + mu*(UK(ii,jj-1) - 2*UK(ii,jj) + UK(ii,jj+1))
        end do
     end do

#ifdef DEBUG
     ! test for convergence
     if (maxval(ABS(UKP1-UK))/dt <= epsilon) then
        converged = .true.
        lastiter = kk
        exit
     end if
#endif
  
     ! swap pointers
     temp => UK
     UK => UKP1
     UKP1 => temp

  end do
  call system_clock(end_count) ! end timing

#ifdef DEBUG
  !! write out result--row-by-row
  do ii=M,1,-1
     write(*,*) UKP1(ii,:)
  end do
#endif

  !! report
  if (count_rate .ne. 1000) then
     write(*,*), 'ERROR timing is not in milliseconds'
     stop
  endif
  write(*,*) 'solving for', N, 'cells'
  write(*,*) 'timestep (s):', dt, '(', int(1.0/dt), 'iters/simulation second)'
  write(*,*) 'maxiter:', int(t/dt)
  if (converged) then
     write(*,*) 'Solutution converged to', epsilon, 'at iteration', lastiter
  end if
  write(*,*) 'solution took: ', real(end_count-start_count), '(ms)'

end Program HeatedPlateForwardEuler