C     KhSa2012-Ex2

! U(1) = u_1(x,t) [Solution on [0, 1/2], reverse mapped
! U(2) = u_2(x,t) [Solution on [1/2, 1]
! U(3) = w(x,t) [Forward Integral from 0 to x] (Since the domain [0, 1/2]
! is reverse mapped onto [0,1], 0 in the computational domain equals 1/2
! in the physical domain. That is, w(0,t) = int_0^1/2 u ds.

C-----------------------------------------------------------------------
      SUBROUTINE F(T, X, U, UX, UXX, FVAL, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C       THIS SUBROUTINE DEFINES THE RIGHT HAND SIDE VECTOR OF THE
C       NPDE DIMENSIONAL PARABOLIC PARTIAL DIFFERENTIAL EQUATION
C                        UT = F(T, X, U, UX, UXX).
C
C-----------------------------------------------------------------------
        IMPLICIT NONE
C SUBROUTINE PARAMETERS:
C INPUT:
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
        DOUBLE PRECISION        T
C                               THE CURRENT TIME COORDINATE.
C
        DOUBLE PRECISION        X
C                               THE CURRENT SPATIAL COORDINATE.
C
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE APPROXIMATION OF THE
C                               SOLUTION AT THE POINT (T,X).
C
        DOUBLE PRECISION        UX(NPDE)
C                               UX(1:NPDE) IS THE APPROXIMATION OF THE
C                               SPATIAL DERIVATIVE OF THE SOLUTION AT
C                               THE POINT (T,X).
C
        DOUBLE PRECISION        UXX(NPDE)
C                               UXX(1:NPDE) IS THE APPROXIMATION OF THE
C                               SECOND SPATIAL DERIVATIVE OF THE
C                               SOLUTION AT THE POINT (T,X).
C
C OUTPUT:
        DOUBLE PRECISION        FVAL(NPDE)
C                               FVAL(1:NPDE) IS THE RIGHT HAND SIDE
C                               VECTOR F(T, X, U, UX, UXX) OF THE PDE.
C-----------------------------------------------------------------------
C
C     ASSIGN FVAL(1:NPDE) ACCORDING TO THE RIGHT HAND SIDE OF THE PDE
C     IN TERMS OF U(1:NPDE), UX(1:NPDE), UXX(1:NPDE).
C
C     Identify solution components: 
C     u1 = U(1), u2 = U(2), w = U(3)

    ! RHS of PDE: Equation 1: u1_t = 4 * u1_xx
      FVAL(1) = 4.0D0 * UXX(1)
    
    ! RHS of PDE: Equation 2: u2_t = 4 * u2_xx
      FVAL(2) = 4.0D0 * UXX(2)
    
    ! ODE for w: Equation 3: w_xi = -(1/2) * u1, in residual form
    ! (The 1/2 factor comes from the change of variables ds = 1/2 d_xi)
      FVAL(3) = UX(3) + 0.5D0 * U(1)

      RETURN
      END

      SUBROUTINE BNDXA(T, U, UX, BVAL, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C       THE SUBROUTINE IS USED TO DEFINE THE BOUNDARY CONDITIONS AT THE
C       LEFT SPATIAL END POINT X = XA.
C                           B(T, U, UX) = 0
C
C-----------------------------------------------------------------------
        IMPLICIT NONE

C SUBROUTINE PARAMETERS:
C INPUT:
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
        DOUBLE PRECISION        T
C                               THE CURRENT TIME COORDINATE.
C
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE APPROXIMATION OF THE
C                               SOLUTION AT THE POINT (T,XA).
C
        DOUBLE PRECISION        UX(NPDE)
C                               UX(1:NPDE) IS THE APPROXIMATION OF THE
C                               SPATIAL DERIVATIVE OF THE SOLUTION AT
C                               THE POINT (T,XA).
C
C OUTPUT:
        DOUBLE PRECISION        BVAL(NPDE)
C                               BVAL(1:NPDE) IS THE BOUNDARY CONTIDITION
C                               AT THE LEFT BOUNDARY POINT.
C-----------------------------------------------------------------------
   
        REAL(8), PARAMETER :: PI = 3.141592653589793D0
        REAL(8) :: TARGET_INT

        ! Nonlocal constraint value (integral from 0 to 1/2)
        TARGET_INT = (1.0D0 / (2.0D0 * PI)) * 
     #            (2.0D0 * SQRT(2.0D0) * EXP(-PI**2 * T / 4.0D0) + 
     #             3.0D0 * EXP(-PI**2 * T / 9.0D0))


C --- Computational Left Boundary (Physical Midpoint x=0.5) ---

        BVAL(1) = U(3) - TARGET_INT        ! Integral value at x=0.5

        BVAL(2) = U(1) - U(2)              ! Continuity: u1(1) = u2(0) [Interface]

        BVAL(3) = UX(1) + UX(2)            ! Smoothness: u1_x(1) = u2_x(0) [Interface]

      RETURN
      END

      SUBROUTINE BNDXB(T, U, UX, BVAL, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C       THE SUBROUTINE IS USED TO DEFINE THE BOUNDARY CONDITIONS AT THE
C       RIGHT SPATIAL END POINT X = XB.
C                           B(T, U, UX) = 0
C
C-----------------------------------------------------------------------
        IMPLICIT NONE
C SUBROUTINE PARAMETERS:
C INPUT:
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
        DOUBLE PRECISION        T
C                               THE CURRENT TIME COORDINATE.
C
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE APPROXIMATION OF THE
C                               SOLUTION AT THE POINT (T,XB).
C
        DOUBLE PRECISION        UX(NPDE)
C                               UX(1:NPDE) IS THE APPROXIMATION OF THE
C                               SPATIAL DERIVATIVE OF THE SOLUTION AT
C                               THE POINT (T,XB).
C
C OUTPUT:
        DOUBLE PRECISION        BVAL(NPDE)
C                               BVAL(1:NPDE) IS THE BOUNDARY CONTIDITION
C                               AT THE RIGHT BOUNDARY POINT.
C-----------------------------------------------------------------------
    
        REAL(8), PARAMETER :: PI = 3.141592653589793D0
        REAL(8) ::  RHS_U2

        ! Dirichlet value for u(1,t)
        RHS_U2 = 0.5D0 * EXP(-PI**2 * T / 9.0D0)

C --- Computational Right Boundary (Physical x=0 for [0,1/2] and x=1 for [1/2,1])

        BVAL(1) = U(3)                     ! Boundary condition for w(0) = 0

        BVAL(2) = U(2) - RHS_U2            ! Boundary condition for u(1,t) = RHS_U1

        BVAL(3) = UX(3) + 0.5D0 * U(1)     ! w-ODE residual 

      RETURN
      END


      SUBROUTINE UINIT(X, U, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C       THIS SUBROUTINE IS USED TO RETURN THE NPDE-VECTOR OF INITIAL
C       CONDITIONS OF THE UNKNOWN FUNCTION AT THE INITIAL TIME T = T0
C       AT THE SPATIAL COORDINATE X.
C
C-----------------------------------------------------------------------
        IMPLICIT NONE
C SUBROUTINE PARAMETERS:
C INPUT:
        DOUBLE PRECISION        X
C                               THE SPATIAL COORDINATE.
C
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
C OUTPUT:
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS VECTOR OF INITIAL VALUES OF
C                               THE UNKNOWN FUNCTION AT T = T0 AND THE
C                               GIVEN VALUE OF X.

        REAL(8), PARAMETER :: PI = 3.141592653589793D0
        REAL(8) :: PX1, PX2
   
        PX1 = 0.5D0 - 0.5D0 * X  ! Segment 1 (xi=0 -> x=0.5, xi=1 -> x=0)
        PX2 = 0.5D0 + 0.5D0 * X  ! Segment 2 (xi=0 -> x=0.5, xi=1 -> x=1)

c-----------------------------------------------------------------------

C
C     ASSIGN U(1:NPDE) THE INITIAL VALUES OF U(T0,X).
C

C     Initial Condition for u_1
      U(1) = COS(PI * PX1 / 2.0D0) + COS(PI * PX1 / 3.0D0)

C     Initial Condition for u_2
      U(2) = COS(PI * PX2 / 2.0D0) + COS(PI * PX2 / 3.0D0)

C     Initial Condition for w
      U(3) = (2.0D0/PI)*SIN(PI * PX1 / 2.0D0) + (3.0D0/PI)*SIN(PI 
     #       * PX1 / 3.0D0)

      RETURN
      END

      SUBROUTINE TRUU(T, X, U, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C     THIS FUNCTION PROVIDES THE EXACT SOLUTIONS
C-----------------------------------------------------------------------
C SUBROUTINE PARAMETERS:
C INPUT:
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
        DOUBLE PRECISION        T
C                               THE CURRENT TIME COORDINATE.
C
        DOUBLE PRECISION        X
C                               THE CURRENT SPATIAL COORDINATE.
C
C OUTPUT:
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE EXACT SOLUTION AT THE
C                               POINT (T,X).
C-----------------------------------------------------------------------
  
        REAL(8), PARAMETER :: PI = 3.141592653589793D0
        REAL(8) :: PX1, PX2, E1, E2

        E1 = EXP(-(PI**2 * T) / 4.0D0)
        E2 = EXP(-(PI**2 * T) / 9.0D0)
        PX1 = 0.5D0 - 0.5D0 * X
        PX2 = 0.5D0 + 0.5D0 * X

C       Exact solution for u_1; this is reversed mapped so it will have
C       to be plotted vs 1 - x
        U(1) = E1 * COS(PI * PX1 / 2.0D0) + E2 * COS(PI * PX1 / 3.0D0)

C       Exact solution for u_2; standard plotting
        U(2) = E1 * COS(PI * PX2 / 2.0D0) + E2 * COS(PI * PX2 / 3.0D0)

C       Exact solution for w 
        U(3) = (2.0D0/PI) * E1 * SIN(PI * PX1 / 2.0D0) + 
     #             (3.0D0/PI) * E2 * SIN(PI * PX1 / 3.0D0)

      RETURN
      END
