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 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 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 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