! pde3.f90  just checking second order approximation         
!      u(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  
!       du/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  
!   d^2u/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  
!       du/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  
!   d^2u/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  
!                                                          
! solve d^2u/dx^2 + d^2u/dy^2 = 12x + 10y = c(x,y)         
! by second order method on square -1,-1 to 1,1            

! second order numerical difference
! d^2u/dx^2 + d^2u/dy^2 =
! 1/h^2(u(x-h,y)-2u(x,y)+u(x+h,y))+
! 1/h^2(u(x,y-h)-2u(x,y)+u(x,y+h)) + O(h^2)
!  = c(x,y)
! solve for new u(x,y) =
! (u(x-h,y)+u(x+h,y)+u(x,y-h)+u(x,y+h)-c(x,y)*h^2)/4

module global  
  implicit none
  integer :: nx = 10 ! one less than 'C' when Fortran subscript starts at zero
  integer :: ny = 10
  double precision, dimension(0:100,0:100) :: us ! solution being iterated        
  double precision, dimension(0:100,0:100) :: ut ! temporary for partial results 
                 ! nx-1 by ny-1 internal, non boundary cells 
  double precision :: xmin = -1.0
  double precision :: ymin = -1.0
  double precision :: xmax =  1.0
  double precision :: ymax =  1.0
  double precision :: h   ! x  and  y step = (xmax-xmin)/(n-1) = (ymax-ymin) 
  double precision :: diff=1.0d100  ! sum abs difference form last iteration 
  double precision :: error         ! sum of absolute exact-computed         
  double precision :: avgerror      ! average error                          
  double precision :: maxerror      ! maximum cell error                     
  interface
    function u(x, y) result(value) ! solution for boundary and check 
      implicit none
      double precision, intent(in) :: x, y
      double precision :: value
    end function u
  end interface
  interface
    function c(x, y) result(term) ! from solve
      implicit none
      double precision, intent(in) :: x, y
      double precision :: term
    end function c
  end interface
end module global

program pde3
  use global
  implicit none
  integer :: iter=0
  
  print *, "pde3.f90 running"
  h = (xmax-xmin)/nx ! step size in x 
  ! = (ymax-ymin)/ny ! step size in y 
  print *, "xmin=",xmin, "xmax=",xmax, "nx=",nx, "h=",h
  print *, "ymin=",ymin, "ymax=",ymax, "ny=",ny, "h=",h
  
  call printexact()
  call init_boundary()
  print "(/a)", "initial condition"
  call printus()
  
  do while(diff>0.0001 .and. iter<200)
    call iterate()
    call check()
    avgerror = error/((nx-1)*(ny-1))
    print "('iter=',i4,' diff=',f8.4,' error=',f9.4,' avg=',f8.4,' max=',f8.4)", &
          iter, diff, error, avgerror, maxerror
    iter = iter+1
    if(mod(iter,50).eq.0) then
      call printus()
      call printdiff()
    end if
  end do
end program pde3 

function u(x, y) result(value) ! solution for boundary and check 
  implicit none
  double precision, intent(in) :: x, y
  double precision :: value
  value = x*x*x + y*y*y + 2.0*x*x*y + 3.0*x*y*y + 4.0*x + 5.0*y +6.0
end function u

function c(x, y) result(term) ! from solve
  implicit none
  double precision, intent(in) :: x, y
  double precision :: term
  term = 12.0*x + 10.0*y
end function c

function u00(i, j) result(value) ! new iterated value 
  use global
  implicit none
  integer, intent(in) :: i, j
  double precision :: value
  value = ( us(i-1,j) + us(i+1,j) + us(i,j-1) + us(i,j+1) - &
           c(xmin+i*h,ymin+j*h)*h*h )/4.0
end function u00

subroutine init_boundary()
  use global
  implicit none
  integer i, j
  
  ! set boundary on four sides 
  do i=0,nx
    us(i,0)  = u(xmin+i*h, ymin)
    us(i,ny) = u(xmin+i*h, ymax)
  end do
  do j=0,ny
    us(0,j)  = u(xmin, ymin+j*h)
    us(nx,j) = u(xmax, ymin+j*h)
  end do
  ! zero internal cells 
  do i=1,nx-1
    do j=1,ny-1
      us(i,j) = 0.0
    end do
  end do
end subroutine init_boundary 

subroutine iterate()
  use global
  implicit none
  integer :: i, j

  interface
    function u00(i, j) result(value) ! new iterated value 
      implicit none
      integer, intent(in) :: i, j
      double precision :: value
    end function u00
  end interface
  
  diff = 0.0
  do i=1,nx-1
    do j=1,ny-1
      ut(i,j) = u00(i,j) 
      diff = diff + abs(ut(i,j)-us(i,j))
    end do
  end do
  ! copy temp to solution 
  do i=1,nx-1
    do j=1,ny-1
      us(i,j) = ut(i,j)
    end do
  end do
end subroutine iterate 

subroutine check() ! typically not known, instrumentation 
  use global
  implicit none
  integer :: i, j
  double precision :: uexact, err
  
  error = 0.0
  maxerror = 0.0
  do i=1,nx-1
    do j=1,ny-1
      uexact = u(xmin+i*h,ymin+j*h)
      err = abs(us(i,j)-uexact)
      error = error + err
      if(err>maxerror) maxerror=err
    end do
  end do
end subroutine check 

subroutine printus()
  use global
  implicit none
  integer :: i, j
  
  print "(/a)", "computed solution"
  do i=0,nx
    print "(11f7.3)", (us(i,j), j=0,min(ny,10))
  end do
end subroutine printus 

subroutine printexact()
  use global
  implicit none
  integer :: i, j
  
  print "(/a)", "exact solution"
  do i=0,nx
    print "(11f7.3)", (u(xmin+i*h,ymin+j*h),j=0,min(ny,10))
  end do
end subroutine printexact 

subroutine printdiff()
  use global
  implicit none
  integer :: i, j
  
  print "(/a)", "difference exact-computed solution" 
  do i=1,nx-1
      print "(7x, 10f7.4)", (u(xmin+i*h,ymin+j*h)-us(i,j), j=1,min(ny-1,10))
  end do
end subroutine printdiff