Program OneD_CrankNicholson

!Algortithm description: 
!***The 1-D fractional difussion equation is solved by a time centered Crank Nicholoson (like)
!		method with extrapolation in space to get a second order accurate solution.
!
!Start Date:	March	10,	2005
!Last Edited:	July	07,	2005 
!By: Charles Tadjeran
!Acknowledgement: This work was supported by NSF grant DMS-0139927.
!
!The 1-D initial boundary value fDE being solved is : 
!	del(u)/del(t) = d(x,t) * del^alpha(u)/del(x^alpha) +  q(x,t)
!	over the region:  x_low < x < x_high
!	and time limits: T_begin < t <= T_end
!	where u = u(x,t)
!	
!	The initial condition is specified by: u(x, t=T_begin) = f(x)
!	Dirichlet boundary conditions are assumed on all boundaries:
!	u(x=x_low,	t)	= xbl(t) 
!	u(x=x_high,	t)	= xbh(t)
!
!	Additionally:
!	d = d(x,t) >=0 is the dispersive coefficient function
!	q = q(x,t)     is the forcing (source/sink) term
!
!	It is assumed that 1<= alpha <= 2 	fractional order of the diffusion term

!Define/set variables:
!		x_low	:low limit on x
!		x_high	:high limit on x	( x_low <=  x <= x_high )
!		NX_input		:number of subintervals in x-direction
!		NT		:number of subintervals (or timesteps) in time
!		alpha	:the fractional order

IMPLICIT NONE

real	:: x_low, x_high
real	:: alpha		!fractional derivative order for the diffusive term
integer	:: NX_input,NX			!number of subintervals in x direction
integer	:: NT			!number of subintervals in time
integer	:: ngrid		!one of the two grids (1=coarse, 2=fine)
COMMON/equation_params/x_low, x_high, NX, alpha		 

real    :: T_begin		!start integration time for t
real	:: T_end		!end integration time t : on T_begin <= t <= T_end

real	:: d			!diffusion coefficient d= d(x,t)
real	:: q			!forcing function in the ADE where q = q(x,t)

real	:: f			!initial condition u(x, t=T_begin) = f(x).

real	:: xbl			!x_low (Dirichlet) boundary function	u(x=x_low,t)	= xbl(t)
real	:: xbh			!x_high (Dirichlet) boundary function	u(x=x_high,t)	= xbh(t) 

real    :: del_x, del_t	! set from user inputs NX, NT

integer :: npts = 401  !max number of x subintervals.
!	***RUNTIME CONSTRAINT: 2*NX_input  to be <= npts. For larger values change dimensions below
real    :: x(0:400),  u(0:400)		!x and solution arrays 
real	:: A(0:400,0:400), rhs(0:400)			!coefficient matrix & high hand side vector of equations		!
real    :: Galpha_wt(0:400)		!array of Grunwald weights- computed once here
real    :: G_wt_scale_x			!basically 0-th Grunwald weights

!some local working variables
integer	:: i, j, k, nstep
real	:: tn, xterm_wt, x_diffuse_wt, ratio, temp, u_exact, sln_err, abs_sln_err, extrap_err
real	:: tnm1, tn_mid
real	:: deriv_term			!temp derivative term computation
real	:: u_coarse(0:400)		!coarse grid solution
real	:: u_ext(0:200)		  !extrapolated solution
integer	:: checkstatus		!file open/close status report
character(40)		:: savefilename !output file containing sim info and output at extraction well

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


!Read prompts for the following user inputs: x_low, x_high
!write (*,*) "Input values for  x_low, x_high ?"
!read (*,*) x_low, x_high
!Read prompts for the following user inputs: T_begin, T_end, alpha	
!write (*,*) "Input values for  T_begin, T_end?"
!read (*,*) T_begin, T_end
write (*,*) "Setting default hard-coded values to save manual input"
	alpha	= 1.8
	x_low	= 0.0
	x_high	= 1.0
	T_begin = 0.0
	T_end	= 1.0


write (*,*) "Input values for  NX, NT ?"
read (*,*) NX_input, NT

del_t = (T_end - T_begin)/NT

!run some quick checks (not exhaustive)
if ((T_begin >= T_end) .OR. (NT <1) .OR. (NX_input <2) .OR. (x_low >= x_high) .OR. (2*NX_input>npts) ) then
	Write (*,*) "*** Invalid user input. Hit <Enter> to terminate!"
	read (*,*)
	stop
endif


write (*,*) "Input a FILENAME to save the solution at 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_low   =", x_low
Write(UNIT=20,FMT=*) "x_high  =", x_high
Write(UNIT=20,FMT=*) "NX_input =", NX_input
Write(UNIT=20,FMT=*) "T_begin  =", T_begin
Write(UNIT=20,FMT=*) "T_end    =", T_end
Write(UNIT=20,FMT=*) "NT       =", NT


do ngrid = 1,2		!do on two gridsizes, the second with twice the number of subintervals
					!ngrid=1 is the coarse grid, and ngrid=2 is the fine grid
	NX = NX_input * ngrid

!Initialize the discretized problem settings
del_x = (x_high - x_low) / real(NX)

!Set the gridpoints
do i = 0 , NX
    x(i)    = x_low + real(i) * del_x
enddo

!Set the initial condition on the gridpoints
do i = 0 , NX
	u(i)	= f(x(i))	!set initial condition at time t = T_begin
end do

if (del_t <= almost_zero) then		!note NT is user input
    write (*,*) "Terminating integration: starting timestep is too small. Hit <Enter> to terminate!"
	read (*,*)
    stop
else
	write (*,*) " Timestep = ",del_t
end if


!We'll use normalized (ie, with 1/del_x**alpha factored out) Grunwald weights:
!The x-direction
Galpha_wt(0) = 1.0
do i = 0 , NX-1
    Galpha_wt(i+1) = - Galpha_wt(i) * (alpha - i)/real(i+1)
end do

!  !
!For each tn, solve the resulting system
!
!First compute some recurring weights to ease computing
!G_wt_scale_x	= 1.0/(del_x**alpha)		!same as G_wt(0),factor used in Grunwald weights
xterm_wt		= del_t / ( 2 * (del_x**alpha))

tn = T_begin					!0th time step
do nstep = 1, NT				!Starts the main timestep loop
	tnm1 = tn
	tn = T_begin + nstep * del_t
	tn_mid = (tn + tnm1)/2
!	Set up matrix equations for the finite differences

		
	do i = 1 , NX-1			!compute matrix of coeffs for row i (x-direction implicit)
		x_diffuse_wt		= d(x(i),tn_mid) * xterm_wt
		rhs(i) = 0.0
		do k = 0, i+1		!column k
			A(i,k)	= - x_diffuse_wt * Galpha_wt(i-k+1)
			rhs(i)  = rhs(i) + Galpha_wt(k) * u(i-k+1)
		enddo		!do k
		A(i,i) = A(i,i) + 1.0
!		the high-hand side of the above equation (row i)
		rhs(i) = u(i) + x_diffuse_wt * rhs(i) + del_t * q(x(i), tn_mid)
	enddo			!do i


!And also the Dirichlet boundary conditions on the low/high boundaries(other entries =0)
	A(0,0)		= 1.0
	A(0,1)		= 0.0
	rhs(0)		= xbl(tn)	!this is not a consistent boundary condition, not important as we are expecting a zero boundary condition

	A(NX,NX-1)	= 0.0
	A(NX,NX)	= 1.0
	rhs(NX)		= xbh(tn)			

!	Solve the resulting matrix system - note A is"sum of a lower triangular and a super-diagnoal matrices
!
!	First do a backward sweep to make it a lower triangular system
!
!	First subtract a multiple of the last row with two non-zero entries	
	ratio			= A(NX-1,NX) / A(NX,NX)	
	A(NX-1,NX-1)	= A(NX-1,NX-1)- ratio * A(NX,NX-1)
	rhs(NX-1) = rhs(NX-1) - ratio * rhs(NX)
!	Remaining rows are not as sparse as the last row
	do i = NX-1, 2, -1
		if (abs(A(i,i)) <= almost_zero) then		!diagonal entry almost zero
			write(*,*) " Singular coeff matrix, at i=", i, " Press <Enter> to terminate! "
			read (*,*) 
			stop
		else
			ratio =	 A(i-1,i) / A(i,i)	
		endif

		do j = 0, i-1
			A(i-1,j) = A(i-1,j) - ratio * A(i,j)
		enddo	!do j
		rhs(i-1) = rhs(i-1) - ratio * rhs(i)
	enddo	 !do i

	u(0) = rhs(0)			
	do i = 1, NX-1
		temp = 0.0
		do j = 0,i-1	
			temp = temp + A(i,j) * u(j)
		enddo	!done j
		u(i) = (rhs(i) - temp) /A(i,i)
	enddo	!done i

	u(NX) = rhs(NX) 		!= Dirichlet boundary condition

enddo		!Ends the main timestep loop


!Print the final coarse/fine girdsize results - compute exact solution and absolute error for the test case
100	Format	 (4(F12.4,3x))
write(*,*)  "// x(i) , u_exact, u(i), sln_err at T_end =", tn
abs_sln_err = 0.0 
do i = 1 ,NX-1
	u_exact = exp(-tn) * (x(i)**3)
	sln_err = abs(u_exact - u(i))
	abs_sln_err = max(sln_err, abs_sln_err)
	Write(UNIT=20,FMT = 100) x(i) , u_exact, u(i), sln_err
	write(*,*) "/"	
enddo

Write(UNIT=20,FMT = *) " On ngrid =", ngrid
Write(UNIT=20,FMT = *) " Absolute  error =", abs_sln_err

!Resutls from one gridsize are available now - save for Tend
if (ngrid == 1) then		!save the results for the coarse grid
	do i = 1, NX-1
		u_coarse(i) = u(i)
	enddo
endif

enddo		!ngrid loop

!Now carry out spatial extrapolation
do i =1, NX_input-1
	u_ext(i) = 2* u(2*i) - u_coarse(i)
enddo



!Check and write-to-file the accuracy of the extrapolated solution vs the exact solution for the test case
write(*,*)  "// x(i) , u_exact, u_ext(i), sln_err at T_end =", tn
extrap_err = 0.0 
do i = 1, NX_input-1
	u_exact = exp(-tn) * ((x(2*i))**3)
	sln_err = abs(u_exact - u_ext(i))
	extrap_err = max(sln_err, extrap_err)
	Write(UNIT=20,FMT = 100) x(i) , u_exact, u_ext(i), sln_err
	write(*,*) "/"	
enddo


Write(UNIT=20,FMT = *) "Absolute extrapolted error = ", extrap_err

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

END Program OneD_CrankNicholson
		
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function d(x, t)		!The diffusivity function in the fractional diffusion equation
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: alpha	!fractional derivative order for the x- diffusive term
	integer	:: NX		!number of subintervals in x  direction
	real	:: x, t, d

COMMON/equation_params/x_low, x_high, NX, alpha		 

! TEST CASE : solution u(x,y,t) = exp(-t)* x^3 & alpha = 1.8
	d = 0.18363375 * (x**2.8)		!gamma(2.2)/6 = 1.10180249088/6 =0.18363375
	
		
End Function d

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function q(x,t)			!The forcing function in the differential eqn
	real	:: q
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: alpha	!fractional derivative order for the x- diffusive term
	integer	:: NX		!number of subintervals in x  direction
	real	:: x, t

COMMON/equation_params/x_low, x_high, NX, alpha		 


!	TEST CASE : solution u(x,t) = exp(-t)* x^3  & alpha = 1.8
	q = -(1.0 + x) * exp(-t) * (x**3) 

End Function q

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function f(x)			!The initial condition u(x,0) = f(x) 
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: alpha	!fractional derivative order for the x- diffusive term
	integer	:: NX		!number of subintervals in x  direction
	real	:: x, t, f

COMMON/equation_params/x_low, x_high, NX, alpha		 

!	TEST CASE : solution u(x,t) = exp(-t)* x^3 ;  & alpha = 1.8
	f = x**3

End Function f

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function xbl(t)			!Dirichlet boundary condition at x_low ; u(x=x_low, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: alpha	!fractional derivative order for the x- diffusive term
	integer	:: NX		!number of subintervals in x direction
	real	:: x, t, xbl

COMMON/equation_params/x_low, x_high, NX, alpha		 
!	TEST CASE : solution u(x,t) = exp(-t)* x^3 & alpha = 1.8
	xbl = 0.0		! zero boundary function at x = x_low = 0

End Function xbl

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function xbh(t)			!Dirichlet boundary condition at x_high: u(x=x_high, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: alpha	!fractional derivative order for the x- diffusive term
	integer	:: NX		!number of subintervals in x  direction
	real	:: x, t, xbh

COMMON/equation_params/x_low, x_high, NX, alpha		 

!	TEST CASE (1): solution u(x,t) = exp(-t)* x^3 ;  & alpha = 1.8
	xbh = exp(-t)* (x_high**3)

End Function xbh