! pde22_eq.f90  2D uniform grid boundary value PDE 
!               known solution to test method                         
!               u(x,y) = 2^16 * x^4 *(1-x)^4 * y^4 * (1-y)^4
!                                                                      
!               differential equation to solve                       
!               d^2u/dx^2 + d^2u/dy^2 = c(x,y)
!               Uxx(X,Y) := 2^16
!                           12*X^2 * (1-X)^4 * Y^4 * (1-Y)^4 -
!                           32*X^3 * (1-X)^3 * Y^4 * (1-Y)^4 +
!                           12*X^4 * (1-X)^2 * Y^4 * (1-Y)^4
!               Uyy(X,Y) := 2^16
!                           12*X^4 * (1-X)^4 * Y^2 * (1-Y)^4 -
!                           32*X^4 * (1-X)^4 * Y^3 * (1-Y)^3 +
!                           12*X^4 * (1-X)^4 * Y^4 * (1-Y)^2
!               c(x,y) := Uxx(x,y) + Uyy(x,y)
!               uniform grid on rectangle 0.0,1.0 to 0.0,1.0                
!                                                                   
! set up and solve system of linear equations                       
!                                                                   
! create two dimensional array with                                 
! (nx-2)*(ny-2) rows, one for each equation [nx in fortran is last subscript]
!
! uses nuderiv.f90
! uses inverse.f90
! uses simeqb.f90

module global  
  implicit none
  integer, parameter :: nx = 10
  integer, parameter :: ny = 10 ! includes boundary
  integer, parameter :: neqn = (nx-2)*(ny-2)
  integer, parameter :: cs = neqn+1
  double precision :: xmin = 0.0
  double precision :: xmax = 1.0
  double precision :: ymin = 0.0
  double precision :: ymax = 1.0
  double precision :: hx
  double precision :: hy
  double precision, dimension(neqn) :: us ! solution being iterated        
  double precision, dimension(neqn,cs) :: ut ! temporary for partial results 
  double precision, dimension(nx) :: xg
  double precision, dimension(ny) :: yg
  double precision, dimension(nx) :: cxx    ! nuderiv coefficients 
  double precision, dimension(ny) :: cyy
  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

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

    subroutine simeqb(n, B, X)
      implicit none
      integer, intent(in) :: n
      double precision, dimension(n,n+1), intent(inout) :: B
      double precision, dimension(n), intent(out) :: X
    end subroutine simeqb

    subroutine nuderiv(order, npoint, point, X, C)
      implicit none
      integer, intent(in) :: order
      integer, intent(in) :: npoint
      integer, intent(in) :: point
      double precision, dimension(npoint), intent(in) :: X
      double precision, dimension(npoint), intent(inout) :: C
    end subroutine nuderiv

    function S(i,j) result(k)
      implicit none 
      integer, intent(in) :: i
      integer, intent(in) :: j
      integer :: k
    end function S

    subroutine print_B()  ! only for checking equations
    end subroutine print_B

    subroutine check_rows()  ! only for checking
    end subroutine check_rows

    subroutine check_soln()
    end subroutine check_soln

    subroutine write_soln2(n,Ux)
      implicit none
      integer, intent(in) :: n
      double precision, dimension(n), intent(in) :: Ux
    end subroutine write_soln2

  end interface

end module global

program pde22_eq
  use global
  implicit none
  integer :: i, j
  double precision :: start, now
  integer :: clock_count, clock_rate, clock_max
  
  print *, "pde22_eq.f90 running"
  print *, "differential equation to solve"                       
  print *, "d^2u/dx^2 + d^2u/dy^2  = c(x,y)"
  print *, "Uxx(X,Y) := 2^16 *"
  print *, "            12*X^2 * (1-X)^4 * Y^4 * (1-Y)^4 -"
  print *, "            32*X^3 * (1-X)^3 * Y^4 * (1-Y)^4 +"
  print *, "            12*X^4 * (1-X)^2 * Y^4 * (1-Y)^4"
  print *, "Uyy(X,Y) := 2^16 *"
  print *, "            12*X^4 * (1-X)^4 * Y^2 * (1-Y)^4 -"
  print *, "            32*X^4 * (1-X)^4 * Y^3 * (1-Y)^3 +"
  print *, "            12*X^4 * (1-X)^4 * Y^4 * (1-Y)^2"
  print *, "c(x,y) := Uxx(x,y) + Uyy(x,y)"
  print *, " " 
  print *, "uniform grid on rectangle 0.0,1.0 to 0.0,1.0" 
  print *, "known solution to test method"
  print *, "u(x,y) = 2^16 * x^4 *(1-x)^4 * y^4 * (1-y)^4"
  print *, "xmin=",xmin, "xmax=",xmax, "nx=",nx
  print *, "ymin=",ymin, "ymax=",ymax, "ny=",ny

  call system_clock(count=clock_count, count_rate=clock_rate, &
                    count_max=clock_max)
  ! print *, "clock_count=", clock_count, ", clock_rate=", &
  ! clock_rate, ', clock_max=', clock_max
  call system_clock(count=clock_count)
  start = (1.0D0*clock_count)/clock_rate
  print *, "start at ", start, " seconds" ! should be zero

  hx = (xmax-xmin)/(nx-1)
  hy = (ymax-ymin)/(ny-1)
  do i=1,nx
    xg(i) = xmin+(i-1)*hx
    print *, "xg(",i,")=",xg(i)
  end do
  do j=1,ny
    yg(j) = ymin+(j-1)*hy
    print *, "yg(",j,")=",yg(j)
  end do
  print *, "u(0.5,0.5)=",u(0.5d0,0.5d0)
  print *, "c(0.5,0.5)=",c(0.5d0,0.5d0)
  
  call init_matrix()
  print *, "matrix initialized"
  ! call print_B() ! debug print
  call check_rows()
  call simeqb(neqn, ut, us)

  call system_clock(count=clock_count)
  now = (1.0D0*clock_count)/clock_rate
  print "('solving took ',f7.3,' seconds')", now-start
  call check_soln()
  call print_soln()
  call check()
  call write_soln2(neqn, us)
  call system_clock(count=clock_count)
  now = (1.0D0*clock_count)/clock_rate
  print "('pde22_eq.f90 finished in ',f7.3,' seconds')", now-start
end program pde22_eq 

function u(x, y) result(value) ! solution for boundary and check 
  implicit none
  double precision, intent(in) :: x, y
  double precision :: value
  value = 2**16 * x**4 * (x-1.0)**4 * y**4 * (y-1.0)**4
end function u

function c(x, y) result(term) ! from solve
  implicit none
  double precision, intent(in) :: x, y
  double precision :: term
  term =  65536.0*( &
          12.0*X*X* &
          (1.0-X)*(1.0-X)*(1.0-X)*(1.0-X)*Y*Y*Y*Y*(1.0-Y)*(1.0-Y)*(1.0-Y)*(1.0-Y)- &
          32.0*X*X*X* &
          (1.0-X)*(1.0-X)*(1.0-X)*        Y*Y*Y*Y*(1.0-Y)*(1.0-Y)*(1.0-Y)*(1.0-Y)+ &
          12.0*X*X*X*X* &
          (1.0-X)*(1.0-X)*                Y*Y*Y*Y*(1.0-Y)*(1.0-Y)*(1.0-Y)*(1.0-Y)+ &
          12.0*X*X*X*X* &
          (1.0-X)*(1.0-X)*(1.0-X)*(1.0-X)*Y*Y*    (1.0-Y)*(1.0-Y)*(1.0-Y)*(1.0-Y)- &
          32.0*X*X*X*X* &
          (1.0-X)*(1.0-X)*(1.0-X)*(1.0-X)*Y*Y*Y*  (1.0-Y)*(1.0-Y)*(1.0-Y)+ &
          12.0*X*X*X*X* &
          (1.0-X)*(1.0-X)*(1.0-X)*(1.0-X)*Y*Y*Y*Y*(1.0-Y)*(1.0-Y))
end function c

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=2,nx-1
    do j=2,ny-1
      uexact = u(xg(i),yg(j))
      err = abs(us(S(i,j))-uexact)
      error = error + err
      if(err>maxerror) maxerror=err
    end do
  end do
  avgerror = error/neqn
  print "('avgerr=',e15.8,' maxerr=',e15.8)", &
           avgerror, maxerror
end subroutine check 

subroutine print_soln()
  use global
  implicit none
  integer :: i, j
  
  print "(/a)", "exact solution u, computed us, error"
  do i=2,nx-1
    do j=2,ny-1
      print "('xg=',f6.3,', yg=',f6.3,', u=',f6.3,', us=',f9.3,', err=',e15.7)", &
      xg(i), yg(j), u(xg(i),yg(j)), us(S(i,j)), &
      u(xg(i),yg(j))-us(S(i,j))
    end do
  end do
end subroutine print_soln 

subroutine init_matrix()
  use global
  implicit none
  integer :: i, j, k, ii, jj, kk
  integer :: ik, kj, kx, ky, kxky
  
  ! zero internal cells 
  do i=2,nx-1
    do j=2,ny-1
      do ii=2,nx-1
        do jj=2,ny-1
          ut(S(i,j),S(ii,jj)) = 0.0
        end do ! jj
      end do ! ii
      ut(S(i,j),cs) = 0.0
    end do ! j
  end do ! i

  print *, "internal cells zeroed "          
  do i=2,nx-1  ! inside boundary 
    do j=2,ny-1
      ii = S(i,j) ! row of matrix 
      ! for each term 

      ! d^2 u(x,y)/dx^2 
      call nuderiv(2,nx,i,xg,cxx)
      if(i.eq.1 .and. j.eq.1) then
        print *, "cxx"
        do k=0,nx
          print *, cxx(k)
        end do
      end if
      do k=1,nx
	if(k.eq.1 .or. k.eq.nx) then
          ut(ii,cs) = ut(ii,cs) - cxx(k)*u(xg(k),yg(j))
	else
          kj = S(k,j)
          ut(ii,kj) = ut(ii,kj) + cxx(k)
        end if
      end do

      ! d^2 u(x,y)/dy^2 
      call nuderiv(2,ny,j,yg,cyy)
      do k=1,ny
	if(k.eq.1 .or. k.eq.ny) then
          ut(ii,cs) = ut(ii,cs) - cyy(k)*u(xg(i),yg(k))
	else
          ik = S(i,k)
          ut(ii,ik) = ut(ii,ik) + cyy(k)
        end if
      end do

      ! RHS constant 
      ut(ii,cs) = ut(ii,cs) + c(xg(i),yg(j))
      ! end terms 
    end do
  end do
end subroutine init_matrix
  
function S(i,j) result(k)
  use global
  implicit none 
  integer, intent(in) :: i
  integer, intent(in) :: j
  integer :: k
  k = (i-2)*(ny-2)+(j-2)+1  ! i,j in us solution
end function S

subroutine print_B()  ! only for checking equations
  use global
  implicit none
  integer :: i, ii, j, jj 

  print *, "matrix B includes RHS"
  do i=2,nx-1
    do j=2,ny-1
      do ii=2,nx-1
        do jj=2,ny-1
          print "('i=',i2,', j=',i2,', ii=',i2,', jj=',i2,' B=',e15.8)", &
                   i, j, ii, jj, ut(S(i,j),S(ii,jj))
          end do ! jj
        end do ! ii
        print "('i=',i2,', j=',i2,'             RHS=',e15.8)", &
                 i, j, ut(S(i,j),cs)
      end do ! j
    end do ! i
  print *, " "
end subroutine print_B

subroutine check_rows()  ! only for checking
  use global
  implicit none
  ! row is in ut solution coordinates
  integer :: i, j, ii, jj, row, col
  double precision :: sum

  do i=2,nx-1
    do j=2,ny-1
      row = S(i,j)
      sum = 0.0d0
      do ii=2,nx-1
        do jj=2,ny-1
          col = S(ii,jj)
          sum = sum + ut(row,col)*u(xg(ii),yg(jj)) ! u is check solution
        end do ! ii
      end do ! jj
      sum = sum - ut(row,cs)
      if(abs(sum).gt.0.001d0) then
        print *, "i=",i,",j=",j," is bad err=",sum
      end if
    end do ! j
  end do !i
end subroutine check_rows

subroutine check_soln() ! check when solution unknown
  use global
  implicit none
  double precision :: smaxerr, err, x, y, uxx, uyy
  integer :: i, j, ii, jj

  !  ux(x,y) + uy(x,y) = c(x,y)
  print *, "check_soln against PDE"
  smaxerr = 0.0d0
  do i=2,nx-1      ! through dof equations
    call nuderiv(2,nx,i,xg,cxx)
    x = xg(i)
    do ii=2,ny-1
      call nuderiv(2,ny,ii,yg,cyy)
      y = yg(ii)
      err = -c(x,y) ! add other side of equation terms
      uxx = 0.0d0
      do j=1,nx  ! compute numerical derivative at (i,ii)
        if(j.eq.1 .or. j.eq.nx) then
          uxx = uxx + cxx(j)*u(xg(j),yg(ii)) ! boundary
        else
          uxx = uxx + cxx(j)*Us(S(j,ii))
        end if
        end do
        uyy = 0.0
        do jj=1,ny
          if(jj.eq.1 .or. jj.eq.ny) then
            uyy = uyy + cyy(jj)*U(xg(i),yg(jj)) ! boundary
          else
            uyy = uyy + cyy(jj)*Us(S(i,jj))
          end if
        end do
        err = err + uxx + uyy
        if(abs(err) .gt. smaxerr) then
          smaxerr = abs(err)
        end if
      end do ! ii Y
    end do ! i X
  print *, "check soln against PDE max error=", smaxerr
end subroutine check_soln

! write_soln2  for gnuplot
subroutine write_soln2(n, Ux)
  use global
  implicit none
  integer, intent(in) :: n
  double precision, dimension(n), intent(in) :: Ux
  integer :: i, ii

  print *, "writing pde22_eq_f90.dat"
  open(unit=11,file="pde22_eq_f90.dat", action="write")
  do i=1,nx
    do ii=1,ny
      if(i.eq.1 .or. i.eq.nx .or. ii.eq.1 .or. ii.eq. ny) then
        write(unit=11,fmt="(f10.7,' ',f10.7,' ',f12.7)") &
              xg(i), yg(ii), U(xg(i),yg(ii)) ! boundary
      else
        write(unit=11,fmt="(f10.7,' ',f10.7,' ',f12.7)") &
              xg(i), yg(ii), Ux(S(i,ii))     ! computed solution
      end if
    end do
    write(unit=11,fmt="(' ')")
  end do
  write(unit=11,fmt="(' ')")
  close(unit=11)
  print *,  "finished writing pde22_eq_f90.dat"
end subroutine write_soln2