Program TwoSided_Test_Imp_Euler

!Algortithm description: 
!		Set up a two-sided fractional PDE (TSFPDE) and solve by implicit Euler method.		
!		The exact solution and the numerical solution are printed to a file that 					
!		to be viewed by Excel to compare outputs
!		
!*NOTE* We use a shifted Grunwald estimate in the implicit Euler formulation for our problem.
!		which is shown to be unconditionally stable! This implicit method results in a coeff
!		matrix that is not sparse.
!

!Start 02/02/04
!By: Charles Tadjeran

!Describe the TSFPDE being solved is : 
!  del(u)/del(t) =  cp(x,t) * del^alpha (u)/del+(x^alpha) + cn(x,t) * del^alpha (u)/del-(x^alpha)+s(x,t), 
! for x_left < x < x_right ,  T_begin < t <= T_end
! and where 
!	u = u(x,t) ,		solution to the TSFPDE
!	cp = cp(x,t)		coeff function for left-handed  (+) fractional derivative
!	cp = cp(x,t)		coeff function for right-handed (-) fractional derivative
!	s = s(x,t)			source/sink forcing function) 
!	1<= alpha <= 2		fractor
!	u(x, 0) = f(x)		Initial Condition function
!	u(x_left,t) = lb(t)	Dirichlet left boundary function (expected to be = 0 for the fractional PDEs)
!	u(x_right,t)= rb(t) Dirichlet boundry function on the right edge of r domain

!Define/set variables:
!		x_left	:left limit on r
!		x_right	:right limit on r	( x_left < x < x_right )
!		N		:number of mesh points in x-direction
!		alpha	:the fractional order (fractor)
!		T_begin	:later may want to add a lower bound on time t(instead of default 0)	 
	real	:: x_left	!left limit on x
	real	:: x_right	!right limit on x
	real	:: alpha	!fractor
	integer	:: N		!number of subintervals in x direction

COMMON/equation_params/x_left, x_right, N, alpha 		 
COMMON/diff_const_test/gamma3_alpha		! gamma3_alpha = gamma(3-alpha)

real    :: T_begin =0.0	!default left limit on time t
real	:: T_end		!final integration time t : on T_begin <= t <= T_end


real	:: f			!initial concentration u(x,t=T_begin) = f(x).
real	:: s			!forcing function s = s(x,t)

real	:: cp, cn		!coeff functions for the left-handed and right-handed derivative respectivley
real	:: lb			!left  (Dirichlet) boundary function	u(x=x_left,t)  = lb(t)
real	:: rb			!right (Dirichlet) boundary function	u(x=x_right,t) = rb(t) 
						!NOT having the above declaration --> caused runtime errors!
real    :: del_x, del_t	!del_x is set by user input N, # of subintervals in x-direction 
						!del_t is  set by  (Tend-Tbegin)/NSTEPS
integer :: npts = 2001  !max number of x-subintervals. Note that N is REQUIRED to be <= npts
real    :: x(0:1000), u(0:1000) , tn
real	:: A(0:1000,0:1000), P(0:1000)	!coefficient matrix and the right hand side vector of equations		!
real    :: G_wt(0:1000) !array of normalized Grunwald weights - computed only once

real    :: almost_zero = 1.0E-6 ! ~machine precision (for 32-bit word)		

integer	:: checkstatus		!file open/close status report
character(40)	:: savefilename		!output file containing sim output

!NOTE:  Provide following user defined functions as Fortran function calls:
!		Function cp(x,t)	:left-handed coefficient function
!		Function cn(x,t)	:left-handed coefficient function
!		Function f(x)		:initial condition - VANISHES when x is "close" to left/right boundaries
!		Function s(x,t)		:forcing function 
!		Function lb(t)		:left (Dirichlet) boundary function u(x=xleft,t)   = lb(t)
!		Function rb(t)		:right(Dirichlet) boundary function u(x=x_right,t) = rb(t) 

!Read prompts for the following user input for x_left, x_right, T_end, alpha, N	
!write (*,*) "Input values for  T_end, NSTEPS, alpha, N ?"
!read (*,*) T_end, NSTEPS, alpha, gamma3_alpha, N

write (*,*) "Input values for  T_end, NSTEPS, N ?"
read (*,*) T_end, NSTEPS, N

!write (*,*) "Input values for  x_left, x_right  ?"
!read (*,*) x_left, x_right
x_left	= 0.0
x_right = 2.0
del_x = (x_right - x_left) / real(N)		

alpha = 1.8
gamma3_alpha= 0.9181687424			!Have gamma(0.2) = 4.590843712
del_t = (T_end - T_begin)/NSTEPS

write (*,*) "Input a FILENAME to save the simulation data at time t=Tend"
read(*,*) savefilename		
open(UNIT = 20, IOSTAT = checkstatus, FILE =savefilename, STATUS = "NEW")
!write(*,*) "File Open IOSTAT =", checkstatus
Write(UNIT=20,FMT=*) "  alpha =", alpha
Write(UNIT=20,FMT=*) "x_left, x_right=",x_left, x_right
Write(UNIT=20,FMT=*) " T_end, del_t, N =", T_end, del_t, N


!Initialize the discretized problem

do i = 0 , N
    x(i)    = x_left + real(i) * del_x 
	u(i)	= f(x(i))			!initial condition at time t = T_begin
end do

!write(*,*)  "/ Debug Print: Initial condition [x(i), u(i)] at t = T_begin "
do i = 0 , N
	write (*, *) x(i), u(i)
end do
!write(*,*) "/"	

if (del_t <= almost_zero) then		!note del_t is user input
    write (*,*) "Terminating integration: starting timestep is too small= ",del_t
    stop
else
	write (*,*) " Timestep = ",del_t
end if

!We'll use normalized (with 1/del_x**alpha factored out) Grunwald weights instead:
G_wt_scale	=1.0/ (del_x ** alpha)
G_wt(0) = 1.0
do i = 0 , N-1
    G_wt(i+1) = - G_wt(i) * (alpha - i)/real(i+1)
end do

!Debug print below
write(*,*)  "Grunwald weights G_wt(i) with "
write(*,*) " del_x=", del_x, ", alpha=",alpha, ", and G_wt_scale=",G_wt_scale
do i = 0 ,N
	write (*, *) i, G_wt(i)
end do
write(*,*) "/"

roh = del_t / (del_x**alpha)
	
!Start the integration process 
!Then apply implicit Euler method at each timestep.

do NT = 1, NSTEPS 
    tn		=  T_begin + del_t * NT       !next integration time       


!	Set up matrix equations for the finite differences: A U_n = U_n-1 + del_t * S_n = P_n ..
!	First the the N x N coeffs matrix (for u(1,.),u(2,.)....u(N,)) for shifted(=1) Grunwald formulation

    do i = 1 , N-1			!compute matix coeffs for row i 
		cp_wt = cp(x(i),tn) * roh
		cn_wt = cn(x(i),tn) * roh

		do j = i-1, 1, -1    !affected by left-handed frac deriv
			A(i,j)	= - G_wt(i-j+1)	* cp_wt
!?			write(*,*)  " A(i,j)  =",i,j, A(i,j)	!	Debug printing
		enddo

		do j = i+1, N-1      !affected by right-handed frac deriv
			A(i,j)	= - G_wt(j-i+1)	* cn_wt
!?			write(*,*)  " A(i,j)  =",i,j, A(i,j)	!	Debug printing
		enddo

! Adjust and/or set the tri-diagonal band (both left/right handed present in these 3 terms)
		A(i,i-1)	=     - G_wt(0)	*  cn_wt + A(i,i-1)	
		A(i,i)		= 1.0 - G_wt(1)	* (cn_wt + cp_wt)	
		A(i,i+1)	=     - G_wt(0)	*  cp_wt + A(i,i+1)

!Set the right hand side for the matrix equation
		P(i) = u(i) + del_t * s(x(i),tn)

	enddo

!And also the Dirichlet boundary conditions on the left and right boundaries
!	A(0,0)	= 1.0
!	A(N,N)	= 1.0
	P(0) = lb(tn)
	P(N) = rb(tn)

!	do i = 1 , N-1 
!		A(0,i) = 0.0
!		A(N,i) = 0.0
!	enddo

! Start Gaussian Elimination - no pivoting
! Forward elimination
	do i = 1, N-1			!note first and last rows are just Dirichlet conditions
		P(i) = P(i) - A(i,0) * P(0)		!for first row operation
	enddo

	do i = 1, N-2			!subtract a multiple of the i-th row from rows i+1,...N-1
		do k = i+1, N-1
			ratio = A(k,i) / A(i,i)
			do j = i+1, N
				A(k,j) = A(k,j) - ratio * A(i,j)
			enddo

			P(k) = P(k) - ratio * P(i)
		enddo

	enddo
!Back substitution
	u(N) = P(N)
	
	do i= N-1, 1, -1
		do j = i+1, N
			P(i) = P(i) - A(i,j) * u(j)
		enddo

		u(i) = P(i) / A(i,i)
	enddo
	u(0) = P(0)

!We now have the solution vector U at time tn

do i = 0 ,N
	write (*, *) x(i), u(i)		!test print
end do


end do !NSTEPS

!Print the final results
100	Format	 (3(F12.4,3x))
	write(*,*)  "// x(i) , u_exact, u(i) at T_end =", tn
	do i = 0 ,N
		u_exact = 4.0 * exp(-tn) * x(i)**2 * (x(i) -2)**2
		Write(UNIT=20,FMT = 100) x(i), u_exact, u(i)
	end do
	write(*,*) "/"	

Close(UNIT=20, IOSTAT= checkstatus)
write(*,*) "File Open IOSTAT =", checkstatus

end Program TwoSided_Test_Imp_Euler


!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function cp(x,t)
	real	:: cp, x, t
	real	:: x_left	!left limit on x
	real	:: x_right	!right limit on x
	real	:: T_end	!final integration time t : on 0<= t <= T_end
	real	:: alpha	!fractional derivative order for the diffusive term
	integer	:: N		!number of subintervals in r direction
COMMON/equation_params/x_left, x_right, N, alpha 		 
COMMON/diff_const_test/gamma3_alpha	! gamma3_alpha = gamma(3-alpha)

!TEST CASE with exact solution u(x,t) = 4 exp(-t)* x^2 * (x-2) ^2 
cp = gamma3_alpha * (x**1.8)

End Function cp			

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function cn(x,t)
	real	:: cn, x ,t
	real	:: xr_left	!left limit on x
	real	:: r_right	!right limit on x
	real	:: T_end	!final integration time t : on 0<= t <= T_end
	real	:: alpha	!fractional derivative order for the diffusive term
	integer	:: N		!number of subintervals in r direction

COMMON/equation_params/x_left, x_right, N, alpha 		 
COMMON/diff_const_test/gamma3_alpha		! gamma3_alpha = gamma(3-alpha)


!TEST CASE with exact solution u(x,t) = 4 exp(-t)* x^2 * (x-2) ^2 
cn = gamma3_alpha * (2-x)**1.8
		
End Function 

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function s(x, t)		!forcing (source/sink) function in the diff eqn

	real	:: s
	real	:: x_left	!left limit on x
	real	:: x_right	!right limit on x
	real	:: T_end	!final integration time t : on 0<= t <= T_end
	real	:: alpha	!fractional derivative order for the diffusive term
	integer	:: N		!number of subintervals in r direction

COMMON/equation_params/x_left, x_right, N, alpha 		 

!TEST CASE with exact solution u(x,t) = 4 exp(-t)* x^2 * (x-2) ^2 
x2=x*x
x3=x2*x
x4=x3*x

tmx=2-x
tmx2=tmx*tmx
tmx3=tmx2*tmx
tmx4=tmx3*tmx

s=-4 * exp(-t)*(x2*tmx2+8*(x2+tmx2)-20*(x3+tmx3)+9.0909090909*(x4+tmx4))

End Function s

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function f(x)			!initial condition for u(x,t) at t=0
	real	:: f, x
	real	:: x_left	!left limit on x
	real	:: x_right	!right limit on x
	real	:: T_end	!final integration time t : on 0<= t <= T_end
	real	:: alpha	!fractional derivative order for the diffusive term
	integer	:: N		!number of subintervals in r direction
COMMON/equation_params/x_left, x_right, N, alpha 		 

!TEST CASE with exact solution u(x,t) = 4 exp(-t)* x^2 * (x-2) ^2 
f = 4.0 * x * x * (x -2) * (x-2)

End Function f

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function lb(t)			!left Dirichlet boundary condition u(x=x_left,t) = lb(t)
	real	:: lb
	real	:: x_left	!left limit on x
	real	:: x_right	!right limit on x
	real	:: T_end	!final integration time t : on 0<= t <= T_end
	real	:: alpha	!fractional derivative order for the diffusive term
	integer	:: N		!number of subintervals in r direction
COMMON/equation_params/x_left, x_right, N,  alpha 		 

lb = 0.0

End Function lb

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function rb(t)			!right Dirichlet boundary condition u(x=x_right,t) = rb(t)

	real	:: rb
	real	:: x_left	!left limit on x
	real	:: x_right	!right limit on r
	real	:: T_end	!final integration time t : on 0<= t <= T_end
	real	:: alpha	!fractional derivative order for the diffusive term
	integer	:: N		!number of subintervals in r direction
COMMON/equation_params/x_left, x_right, N, alpha 		 

rb = 0.0

End Function rb