! pde_spherical_eq.f90     spherical coordinates
! solve del^2 U(radius,theta,phi) + U(radius,theta,phi) = f(r,t,p)
! given Dirchlet boundary values Ub(r,t,p), and f(r,t,p)
! Much checking in this code, analytic solution is
! U(r,t,p) = r^3 + t^3/216 + p^3/27
!
! x = r cos(theta) sin(phi)    0 <= theta <= 2 Pi
! y = r sin(theta) sin(phi)    0 < phi < Pi restricted for derivatives
! z = r cos(phi)               r < 0     restricted for derivatives
! r = sqrt(x^2+y^2+z^2)
! theta = arctan(y,x) if theta<0 then theta = 2 Pi + arctan(y,x)
! phi = arctan(sqrt(x^2+y^2),z)
!
! use existing partial derivative computation in Cartesian coordinates to
! compute partial derivatives in spherical coordinates
!
! compute del U(radius,theta,phi) r=radius, t=theta angle, p=phi angle
! del U  =  (d/dr U, (1/(r sin(p))) d/dt U, 1/r d/dp U)
!
! compute del^2 U(radius,theta,phi) r=radius, t=theta angle, p=phi angle
! del^2 U  = (1/r^2) d/dr(r^2 d/dr U) +
!            (1/(r^2 sin(p)^2)) d^2/dt^2 U +
!            1/(r^2 sin(p)) d/dp(sin(p) d/dp U)
!          = (2/r) d/dr U + d^2/dr^2 U +
!            (1/(r^2 sin^2(p))) d^2/dt^2 U
!            (cos(p)/(r^2 sin(p)) d/dp U +
!            (1/(r^2)) d^2/dp^2 U +
!
! f(r,t,p) = 12*r + 1/36*t/(r^2*sin(p)^2) +
!             1/9*cos(p)*p^2/(r^2*sin(p)) + 2/9*p/r^2 +
!             r^3 + 1/216*t^3 + 1/27*p^3
!
! uses nuderiv.f90
! uses simeq.f90

module global  
  implicit none
  integer, parameter :: Nr = 5 
  integer, parameter :: Nt = 7
  integer, parameter :: Np = 5  
  integer, parameter :: Nrtp = (Nr-2)*(Nt-2)*(Np-2) ! DOF
  double precision, parameter :: Pi = 3.141592653589793238462643383279502884197
  double precision, dimension(Nr) :: Rg ! r grid values
  data Rg / 0.3, 0.5, 0.7, 0.9, 1.0 /
  double precision, dimension(Nt) :: Tg ! theta grid values
  data Tg / 0.1, 0.7, 1.5, 2.35, 3.25, 4.2, 5.25 /
  double precision, dimension(Np) :: Pg ! p grid values
  data Pg / 0.2, 0.9, 1.7, 2.4, 2.65 /
  double precision, dimension(Nr*Nt*Np) :: Ua  
  double precision, dimension(Nr) :: Cr  ! first deriv wrt r
  double precision, dimension(Nr) :: Crr ! second deriv wrt r
  double precision, dimension(Nt) :: Ctt ! second deriv wrt t
  double precision, dimension(Np) :: Cp  ! first deriv wrt p
  double precision, dimension(Np) :: Cpp ! second deriv wrt p
  double precision, dimension(Nrtp,Nrtp+1) :: A ! no bound, has RHS
  double precision, dimension(Nrtp) :: X! computer solution to U

  interface
    function F(r, t, p) result(term) ! RHS
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: term
    end function F
  end interface
  interface
    subroutine simeq(n, B, X)
      implicit none
      integer, intent(in) :: n
      double precision, dimension(n,n+1), intent(in) :: b
      double precision, dimension(n), intent(out) :: X
     end subroutine simeq
  end interface
  interface
    function Ub(r, t, p) result(value) ! for boundary and solution check 
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Ub
  end interface
  interface
    function Ur(r, t, p) result(value) ! analytic first derivative, for check
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Ur
  end interface
  interface
    function Urr(r, t, p) result(value) ! analytic second derivative, for check
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Urr
  end interface
  interface
    function Ut(r, t, p) result(value) ! analytic first derivative, for check
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Ut
  end interface
  interface
    function Utt(r, t, p) result(value) ! analytic second derivative, for check
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Utt
  end interface
  interface
    function Up(r, t, p) result(value) ! analytic first derivative, for check
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Up
  end interface
  interface
    function Upp(r, t, p) result(value) ! analytic second derivative, for check
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Upp
  end interface
  interface
    function Uuu(r, t, p) result(value)
      implicit none
      double precision, intent(in) :: r, t, p
      double precision :: value
    end function Uuu
  end interface
  interface
    function D2(Urv, Urrv, Uttv, Upv, Uppv, R, P) result(value)
      implicit none
      double precision, intent(in) :: Urv, Urrv, Uttv, Upv, Uppv, R, P
      double precision :: value
    end function D2
  end interface
  interface
    function S(I, J, K) result(index) ! for Ua
      implicit none
      integer, intent(in) :: I, J, K
      integer :: index
    end function S
  end interface
  interface
    function Sd(I, J, K) result(index) ! for A and X
      implicit none
      integer, intent(in) :: I, J, K
      integer :: index
    end function Sd
  end interface


  interface ! external, compile with his file
     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
  end interface
end module global

program pde_spherical_eq
  use global
  implicit none
  integer :: I, J, K, Ii, Jj, Kk
  double precision :: R, T, P, Ana, Cmp, Err
  double precision :: Urv, Urrv, Uttv, Upv, Uppv ! computed derivatives
  double precision :: Maxerrur, Maxerrurr, Maxerrutt, Maxerrup, Maxerrupp, &
                      Maxerr
  
  print *, "pde_spherical_eq.f90 running"
  print *, "differential equation to solve"                       
  print *, "del^2 U(radius,theta,phi) + U(radius,theta,phi) = F(r,t,p)"
  print *, "given Dirchlet boundary values Ub(r,t,p), and F(r,t,p)"
  print *, "Much checking in this code, analytic solution is"
  print *, "U(r,t,p) = r^3 + t^3/216 + p^3/27"

  print *, Nr," r gird values"
  do I=1,Nr
     print *, "Rg(",I,")=",Rg(I)
  end do
  print *, Nt," theta gird values"
  do J=1,Nt
     print *, "Tg(",J,")=",Tg(J)
  end do
  print *, Np," phi gird values"
  do K=1,Np
     print *, "Pg(",K,")=",Pg(K)
  end do

  print *, "analytic solution, partial second derivative, F(R,T,P)"
  do I=1,Nr
    R = Rg(I)
    do J=1,Nt
      T = Tg(J)
      do K=1,Np
	P = Pg(K)
        Ua(S(I,J,K)) = Ub(Rg(I),Tg(J),Pg(K))
        print *, "Ua(",I,",",J,",",K,")=",Ua(S(I,J,K)), &
		 ", D2=",D2(Ur(R,T,P),Urr(R,T,P),Utt(R,T,P), &
                            Up(R,T,P),Upp(R,T,P),R,P)
        print *, "Uuu+U=",Uuu(R,T,P)+Ub(R,T,P),", F=",F(R,T,P)
      end do
    end do
  end do

  do I=1,Nr
    R = Rg(I)
    do J=1,Nt
      T = Tg(J)
      do K=1,Np
        P = Pg(K)
        ! compute Ur, Urr, Utt, Upp and call D2
        call Nuderiv(1, Nr, I, Rg, Cr)
        call Nuderiv(2, Nr, I, Rg, Crr)
        Urv = 0.0
        Urrv = 0.0
        do Ii=1,Nr
           Urv = Urv + Cr(Ii)*Ub(Rg(Ii),T,P)
           Urrv = Urrv + Crr(Ii)*Ub(Rg(Ii),T,P)
        end do
        Uttv = 0.0
        call Nuderiv(2, Nt, J, Tg, Ctt)
        do Jj=1,Nt
           Uttv = Uttv + Ctt(Jj)*Ub(R,Tg(Jj),P)
        end do
        Upv = 0.0
        Uppv = 0.0
        call Nuderiv(1, Np, K, Pg, Cp)
        call Nuderiv(2, Np, K, Pg, Cpp)
        do Kk=1,Np
           Upv = Upv + Cp(Kk)*Ub(R,T,Pg(Kk))
           Uppv = Uppv + Cpp(Kk)*Ub(R,T,Pg(Kk))
        end do

        Cmp = D2(Urv, Urrv, Uttv, Upv, Uppv, R, P)
        Ana = Uuu(R,T,P)
        Err = Ana - Cmp
        Maxerr = Max(Maxerr, Err)
        print *, "i=",I,", j=",J,", k=",K, &
                  ", Ana=",Ana,", Cmp=",Cmp,", err=",Err
        Err = Ur(R,T,P)-Urv
        Maxerrur = Max(Maxerrur, Err)
        print *, "i=",I,", j=",J,", k=",K, &
                  ", Ur=",Ur(R,T,P),", Urv=",Urv,", err=",Err
        Err = Urr(R,T,P)-Urrv
        Maxerrurr = Max(Maxerrurr, Err)
        print *, "i=",I,", j=",J,", k=",K, &
                  ", Urr=",Urr(R,T,P),", Urrv=",Urrv,", err=",Err
        Err = Utt(R,T,P)-Uttv
        Maxerrutt = Max(Maxerrutt, Err)
        print *, "i=",I,", j=",J,", k=",K, &
                  ", Utt=",Utt(R,T,P),", Uttv=",Uttv,", err=",Err
        Err = Up(R,T,P)-Upv
        Maxerrup = Max(Maxerrup, Err)
        print *, "i=",I,", j=",J,", k=",K, &
                  ", Up=",Up(R,T,P),", Upv=",Upv,", err=",Err
        Err = Upp(R,T,P)-Uppv
        Maxerrupp = Max(Maxerrupp, Err)
        print *, "i=",I,", j=",J,", k=",K, &
                  ", Upp=",Upp(R,T,P),", Uppv=",Uppv,", err=",Err
        print *, " "
      end do
    end do
  end do
  print *, "Maxerr Ur=",Maxerrur
  print *, "Maxerr Urr=",Maxerrurr
  print *, "Maxerr Utt=",Maxerrutt
  print *, "Maxerr Up=",Maxerrup
  print *, "Maxerr Upp=",Maxerrupp
  print *, "Maxerr deriv=",Maxerr

  call build_A()
  call simeq(Nrtp, A, X) ! X is solution
  print *, "check against analytic solution"
  Maxerr = 0.0
  do I=2,Nr-1
    do J=2,Nt-1
      do K=2,Np-1
        Err = X(Sd(I, J, K))-Ua(S(I,J,K))
        Maxerr = Max(Maxerr, Abs(Err))
        print *, "(",I,",",J,",",K,") comp=",  &
                 X(Sd(I, J, K)),", ana=",Ua(S(I,J,K)), &
                 ", err=",Err
      end do ! K
    end do ! J
  end do ! I
  print *, "Maxerr soln=",Maxerr

  print *, "pde_spherical_eq.f90 finishing"
end program pde_spherical_eq 
  
subroutine simeq(n, B, X)
  implicit none ! tailored version for this application
  integer, intent(in) :: n
  double precision, dimension(n,n+1), intent(inout) :: B
  double precision, dimension(n), intent(out) :: X
  integer, dimension(n) :: ROW           ! ROW INTERCHANGE INDICES 
  integer :: HOLD , I_PIVOT                ! PIVOT INDICES 
  double precision :: PIVOT                ! PIVOT ELEMENT VALUE 
  double precision :: ABS_PIVOT
  integer :: i, j, k, m

  m = n+1

  ! SET UP ROW  INTERCHANGE VECTORS 
  do k=1,n
    ROW(k) = k
  end do ! k

  ! BEGIN MAIN REDUCTION LOOP 
  do k=1,n

    ! FIND LARGEST ELEMENT FOR PIVOT 
    PIVOT = B(ROW(k),k)
    ABS_PIVOT = abs(PIVOT)
    I_PIVOT = k
    do i=k,n
      if( abs(B(ROW(i),k)) > ABS_PIVOT) then
        I_PIVOT = i
        PIVOT = B(ROW(i),k)
        ABS_PIVOT = abs ( PIVOT )
      end if
    end do ! i

    ! HAVE PIVOT, INTERCHANGE ROW POINTERS 
    HOLD = ROW(k)
    ROW(k) = ROW(I_PIVOT)
    ROW(I_PIVOT) = HOLD

    ! CHECK FOR NEAR SINGULAR 
    if( ABS_PIVOT < 1.0E-10 ) then
      do j=k+1,n+1
        B(ROW(k),j) = 0.0
      end do ! j
      print *, 'redundant row (singular) ', ROW(k)
    else

      ! REDUCE ABOUT PIVOT 
      do j=k+1,n+1
        B(ROW(k),j) = B(ROW(k),j) / B(ROW(k),k)
      end do ! j

      ! INNER REDUCTION LOOP 
      do i=1,n
        if( i .ne. k) then
          do j=k+1,n+1
            B(ROW(i),j) = B(ROW(i),j) - B(ROW(i),k) * B(ROW(k),j)
          end do ! j
        end if
      end do ! i
    end if 
    ! FINISHED INNER REDUCTION 
  end do ! k

  ! END OF MAIN REDUCTION LOOP 
  ! BUILD  X  FOR RETURN, UNSCRAMBLING ROWS 
  do i=1,n
    X(i) = B(ROW(i),n+1)
  end do ! i
end subroutine simeq 

function Ub(r, t, p) result(value) ! for boundary and solution check 
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = r**3 + t**3/216.0 + p**3/27.0
end function Ub

function Ur(r, t, p) result(value) ! analytic first derivative, for check
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = 3.0*r**2
end function Ur

function Urr(r, t, p) result(value) ! analytic second derivative, for check
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = 6.0*r
end function Urr

function Ut(r, t, p) result(value) ! analytic first derivative, for check
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = t**2/27.0
end function Ut

function Utt(r, t, p) result(value) ! analytic second derivative, for check
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = t/36.0
end function Utt

function Up(r, t, p) result(value) ! analytic first derivative, for check
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = p**2/9.0
end function Up

function Upp(r, t, p) result(value) ! analytic second derivative, for check
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value = 2.0*p/9.0
end function Upp

function Uuu(r, t, p) result(value)
  use global
  implicit none
  double precision, intent(in) :: r, t, p
  double precision :: value
  value =  2.0*Ur(R,T,P)/R + Urr(R,T,P) + &
           Utt(R,T,P)/(R*R*sin(P)*sin(P)) + &
           cos(P)*Up(R,T,P)/(R*R*sin(P)) + Upp(R,T,P)/(R*R)
end function Uuu

function D2(Urv, Urrv, Uttv, Upv, Uppv, R, P) result(value)
  implicit none
  double precision, intent(in) :: Urv, Urrv, Uttv, Upv, Uppv, R, P
  double precision :: value
  value = 2.0*Urv/R + Urrv + Uttv/(R*R*sin(P)*sin(P)) + &
           cos(P)*Upv/(R*R*sin(P)) + Uppv/(R*R)
end function D2

function S(I, J, K) result(index)
  use global
  implicit none
  integer, intent(in) :: I, J, K
  integer :: index
  index = (I-1)*Nt*Np + (J-1)*Np + K
end function S

function Sd(I, J, K) result(index)
  ! index of each grid point inside boundary, A and X
  use global
  implicit none
  integer, intent(in) :: I, J, K
  integer :: index
  index = (I-2)*(Nt-2)*(Np-2) + (J-2)*(Np-2) + (K-1)
end function Sd

! f(r,t,p)=Urrs(r,t,p)+Utts(r,t,p)+Upps(r,t,p)+Ub(r,t,p)
! f(r,t,p) = 12*r + 1/36*t/(r^2*sin(p)^2) +
!             1/9*cos(p)*p^2/(r^2*sin(p)) + 2/9*p/r^2 +
!             r^3 + 1/216*t^3 + 1/27*p^3
function F(R, T, P) result(term)
  double precision, intent(in) :: R, T, P
  double precision :: term
  term = 12.0*R + (1.0/36.0)*T/(R**2*sin(P)**2) + &
         (1.0/9.0)*cos(P)*P**2/(R*R*sin(P)) + &
         (2.0/9.0)*P/R**2 + &
         R**3 + (1.0/216.0)*T**3 + (1.0/27.0)*P**3
end function F  ! RHS

subroutine build_A() ! A is discrete system of equations, including Y
                        ! solution is X, from A*X=Y
  use global
  implicit none
  integer :: Row, Idof, Cs
  double precision :: Urrs, Utts, Upps, Ud2ck, Err
  double precision :: Err1, Err2, Err3
  integer :: I, J, K, Ii, Jj, Kk
  double precision :: R, T, P

  ! equation Urrs(r,t,p)+Utts(r,t,p)+Upps(r,t,p)+U(r,t,p) = f(r,t,p)
  ! for every (r,t,p) inside boundary
  print *, "building discrete system of ",Nrtp," equations"
  Cs = Nrtp+1 ! RHS
  ! initialize A and X
  do I=2,Nr-1
    do J=2,Nt-1
      do K=2,Np-1
        Row = Sd(I, J, K) ! equation index
        do Ii=2,Nr-1
          do Jj=2,Nt-1
            do Kk=2,Np-1
              A(Row,Sd(Ii,Jj,Kk)) = 0.0
            end do ! Kk
          end do ! Jj
        end do ! Ii
        A(Row,Cs) = F(Rg(I),Tg(J),Pg(K)) ! RHS
        X(Row) = 0.0
      end do ! K
    end do ! J 
  end do ! I
  ! now build A inside boundary (could use sparse)
  do I=2,Nr-1
    R = Rg(I)
    do J=2,Nt-1
      T = Tg(J)
      do K=2,Np-1
        P = Pg(K)
        Row = Sd(I, J, K) ! equation index
        Ud2ck = 0.0 ! only for development checking

        ! equation Urrs(r,t,p)+Utts(r,t,p)+Upps(r,t,p)+U(r,t,p) = f(r,t,p)
        ! Urrs term
        call nuderiv(1, Nr, I, Rg, Cr)
        call nuderiv(2, Nr, I, Rg, Crr)
        do Ii=1,Nr
          Urrs = 2.0*Cr(Ii)/R + Crr(Ii)
          Ud2ck = Ud2ck + Urrs

          if(Ii.eq.1 .or. Ii.eq.Nr) then ! boundary, negative on RHS
             A(Row,Cs) = A(Row,Cs) - Urrs*Ub(Rg(Ii),T,P)
          else
            Idof = Sd(Ii,J,K)
            A(Row,Idof) = A(Row,Idof) + Urrs
          end if
        end do !Ii

        ! Utts term
        call nuderiv(2, Nt, J, Tg, Ctt)
        do Jj=1,Nt
          Utts = Ctt(Jj)/(R*R*Sin(P)*Sin(P))
          Ud2ck = Ud2ck + Utts
          if(Jj.eq.1 .or. Jj.eq.Nt) then ! boundary, negative on RHS
             A(Row,Cs) = A(Row,Cs) - Utts*Ub(R,Tg(Jj),P)
          else
            idof = Sd(I,Jj,K)
            A(Row,Idof) = A(Row,Idof) + Utts
          end if
        end do ! Jj

        ! Upps term
        call nuderiv(1, Np, K, Pg, Cp)
        call nuderiv(2, Np, K, Pg, Cpp)
        do Kk=1,Np
          Upps = cos(P)*Cp(Kk)/(R*R*sin(P)) + Cpp(Kk)/(R*R)
          Ud2ck = Ud2ck + Upps
          if(Kk.eq.1 .or. Kk.eq.Np) then ! boundary, negative on RHS
            A(Row,Cs) = A(Row,Cs) - Upps*Ub(R,T,Pg(Kk))
          else
            Idof = Sd(I,J,Kk)
            A(Row,Idof) = A(Row,Idof) + Upps
          end if
        end do ! Kk

        ! U term
        A(Row,Row) = A(Row,Row) + 1.0

        ! check
        Err = abs(Ud2ck)
        if(Err.gt. 0.001) then
          print *, "(",I,",",J,",",K,") Ud2ck=",Ud2ck,", err=",Err
        end if
      end do ! K
      if(Row.eq.3) then ! debug print
        Err1 = 0.0
        Err2 = 0.0
        Err3 = 0.0
        do Ii=2,Nr-1
          do Jj=2,Nt-1
            do Kk=2,Np-1
              print *, "(",Ii,",",Jj,",",Kk,") A(1)=", &
                       A(1,Sd(Ii,Jj,Kk)),", A(2)=", &
                       A(2,Sd(Ii,Jj,Kk)),", A(3)=", &
                       A(3,Sd(Ii,Jj,Kk))
              Err1 = Err1 + A(1,Sd(Ii,Jj,Kk))*Ub(Rg(Ii),Tg(Jj),Pg(Kk))
              Err2 = Err2 + A(2,Sd(Ii,Jj,Kk))*Ub(Rg(Ii),Tg(Jj),Pg(Kk))
              Err3 = Err3 + A(3,Sd(Ii,Jj,Kk))*Ub(Rg(Ii),Tg(Jj),Pg(Kk))
            end do ! Kk
          end do ! Jj
        end do ! Ii
        print *, "RHS Cs=",Cs," A(1)=",A(1,Cs),", A(2)=", &
                  A(2,Cs),", A(3)=",A(3,Cs)
        Err1 = Err1 - A(1,Cs)
        Err2 = Err2 - A(2,Cs)
        Err3 = Err3 - A(3,Cs)
        print *, "row check"
        print *, "Err1=",Err1,", Err2=",Err2,", Err3=",Err3
      end if ! debug print
    end do ! J
  end do ! I
end subroutine Build_A