Program Two_D_ADI_FPDE
!Algortithm description:
! The 2-D fractional diffusion equation is solved here by an Implicit Euler method with
! Alternating Directions Implicit (ADI) method..
! A shifted Grunwald estimate for the fractional order terms is used for stability.
! The algorithm applies the derivation in the "problem" statement of the paper
! titled "Finite difference methods for two-dimensional fractional dispersion equation"
! to appear in J. of Computational Physics.
!----------
!Start Date: 10/03/04
!Last Edit: 07/04/05
!By: Charles Tadjeran
!Acknowlegement: This project was supported by NSF grant DMS-0139927.
!----------
!The 2-D initial boundary value FADE being solved is :
! del(u)/del(t) = d(x,y,t) * del^alpha(u)/del(x^alpha) +
! e(x,y,t) * del^alpha(u)/del(y^beta) + q(x,y,t)
! over the rectangular region: x_low < x < x_high, y_low < y < y_high,
! and time limits: T_begin < t <= T_end
! where u = u(x,y,t)
!
! The initial condition is specified by: u(x,y, t=T_begin) = f(x,y)
! Dirichlet boundary conditions are assumed on all boundaries:
! u(x=x_low, y, t) = xbl(y,t)
! u(x=x_high, y, t) = xbh(y,t)
! u(x, y=y_low, t) = ybl(x,t)
! u(x, y=y_high, t) = ybh(x,t)
!
! Additionally:
! d = d(x,y,t) >=0 is the x dispersive coefficient function
! e = e(x,y,t) >=0 is the y dispersive coefficient function
! q = q(x,y,t) >=0 is the forcing (source/sink) term
!
! It is assumed that
! 1<= alpha <= 2 , and 1<= beta <= 2 fractional order of the diffusion terms
!User to define/set variables:
! x_low :low limit on x
! x_high :high limit on x ( x_low <= x <= x_high )
! y_low :low limit on y
! y_high :high limit on y ( y_low <= y <= y_high )
!
! NX :number of subintervals in x-direction
! NY :number of subintervals in y-direction
! NT :number of subintervals (or timesteps) in time
! alpha, beta :the fractional orders
IMPLICIT NONE
real :: x_low, x_high, y_low, y_high
real :: alpha, beta !fractional derivative orders for the diffusive terms
integer :: NX, NY !number of subintervals in x and y directions
integer :: NT !number of subintervals in time
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha, beta
real :: T_begin !start integration time for t
real :: T_end !end integration time t : on T_begin <= t <= T_end
real :: d !x direction diffusion coefficient d= d(x,y,t)
real :: e !y direction diffusion coefficient e= e(x,y,t)
real :: q !forcing function in the ADE where q = q(x,y,t)
real :: f !initial condition u(x, y,t=T_begin) = f(x,y).
real :: xbl !x_low (Dirichlet) boundary function u(x=x_low,y,t) = xbl(y,t)
real :: xbh !x_high (Dirichlet) boundary function u(x=x_high,y,t) = xbh(t)
real :: ybl !y_low (Dirichlet) boundary function u(x,y=y_low,t) = ybl(x,t)
real :: ybh !y-high (Dirichlet) boundary function u(x,y=y_high,t) = ybh(x,t)
real :: del_x, del_y, del_t ! set from user inputs NX, NY, NT
integer :: npts = 401 !max number of x or y subintervals.
! ***RUNTIME CONSTRAINT: NX and NY to be <= npts. For larger values change dimensions below
real :: x(0:400), y(0:400), u(0:400, 0:400) !x, y 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), Gbeta_wt(0:400) !array of Grunwald weights- computed once here
real :: G_wt_scale_x, G_wt_scale_y !basically 0-th Grunwald weights
!some local working variables
integer :: i, j, k, nxfix, nyfix, nstep
real :: tn, xterm_wt,yterm_wt, x_diffuse_wt, y_diffuse_wt, ratio, temp, u_exact, sln_err, abs_sln_err
real :: u_int(0:400, 0:400) !temp array of intermediate solution values prior to updates
real :: deriv_term !temp derivative term computation
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)
!*******Begin USER NOTE 1
!Begin SETTING UP FOR USER INPUTS below
!Read prompts for the following user inputs: x_low, x_high, y_low, y_high, T_begin, T_end, alpha, beta
!write (*,*) "alpha, beta?"
!read (*,*) alpha, beta
!write (*,*) "Input values for x_low, x_high, y_low, y_high ?"
!read (*,*) x_low, x_high, y_low, y_high
!write (*,*) "Input values for T_begin, T_end?"
!read (*,*) T_begin, T_end
!FOR GENERAL CASES, un-comment the above user read and write prompt statement, and comment out the following section
! which was designed to easily provide for testing and minimizes typing many inputs
! during each test session to reduce user input errors
!BEGINNING of TEST SESSION INPUTS -->
alpha = 1.8
beta = 1.6
x_low = 0.0
x_high = 1.0
y_low = 0.0
y_high = 1.0
T_begin = 0.0
T_end = 1.0
!END OF TEST SESSION INPUTS <-- which should be commented out if a new test case is constructed
!*******End USER NOTE 1
write (*,*) "Input values for NX, NY, NT ?"
read (*,*) NX, NY, NT
NT = 2 * ((NT+1)/2) !To make sure that we have an even number of timesteps
del_t = (T_end - T_begin)/NT
!run some quick checks (not exhaustive)
if ((T_begin >= T_end) .OR. (NT <2) .OR. (NX <2) .OR. (NY <2) .OR. (x_low >= x_high) .OR. (y_low >= y_high) .OR. (NX>npts) .OR. (NY>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=*) "beta =", beta
Write(UNIT=20,FMT=*) "x_low =", x_low
Write(UNIT=20,FMT=*) "x_high =", x_high
Write(UNIT=20,FMT=*) "NX =", NX
Write(UNIT=20,FMT=*) "y_low =", y_low
Write(UNIT=20,FMT=*) "y_high =", y_high
Write(UNIT=20,FMT=*) "NY =", NY
Write(UNIT=20,FMT=*) "T_begin =", T_begin
Write(UNIT=20,FMT=*) "T_end =", T_end
Write(UNIT=20,FMT=*) "NT =", NT
!Initialize the discretized problem settings
del_x = (x_high - x_low) / real(NX)
del_y = (y_high - y_low) / real(NY)
!Set the gridpoints
do i = 0 , NX
x(i) = x_low + real(i) * del_x
enddo
do i = 0 , NY
y(i) = y_low + real(i) * del_y
enddo
!Set the initial condition on the gridpoints
do i = 0 , NX
do j = 0 , NY
u(i, j) = f(x(i), y(j)) !set initial condition at time t = T_begin
end do
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
!The y-direction
Gbeta_wt(0) = 1.0
do i = 0 , NY-1
Gbeta_wt(i+1) = - Gbeta_wt(i) * (beta - i)/real(i+1)
end do
! !
!For each tn, solve the resulting system: first sweep in x-direction, then in y-direction
!
!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
G_wt_scale_y = 1.0/(del_y**beta)
xterm_wt = del_t * G_wt_scale_x
yterm_wt = del_t * G_wt_scale_y
nstep = 0 !0th time step
do while ( nstep < NT) !main time step forward move: 2 time-steps thru this loop each time
!XXXXXXXXXXXXXXXXXXXXX-DIRECTION
nstep = nstep+1
tn = T_begin + nstep * del_t
do nyfix = 1, NY-1
! Set up matrix equations for the finite differences: Ax U_j^n+1 = U_j^n + delt * Q_j^n+1 = P_n+1 ..
! First the NX x NX coeffs for the super-triangular matrix in x-direction
do i = 1 , NX-1 !compute matrix of coeffs for row i (x-direction implicit)
x_diffuse_wt = d(x(i), y(nyfix),tn) * xterm_wt
! y_diffuse_wt = e(x(i), y(nyfix),tn) * yterm_wt !not sure if best at t_n or t_{n-1}
do k = 0, i+1 !column k
A(i,k) = - x_diffuse_wt * Galpha_wt(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,nyfix) + del_t * q(x(i), y(nyfix),tn)
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(y(nyfix), tn) !this is not a consistent boundary condition, not important as we are expecting a zero boundary condition
!!
!! deriv_term = 0.0
!! do i= 0,nyfix
!! deriv_term = deriv_term + Gbeta_wt(nyfix-i) * xbl(y(i+1), tn)
!! enddo
!! rhs(0) = xbl(y(nyfix), tn) - yterm_wt * e(x_low, y(nyfix),tn) * deriv_term
!!
!!same as below
rhs(0) = 0 !intermediate sln = (1- e(x,y,t) del^\beta U/del y^\beta)U^{n+1}
A(NX,NX-1) = 0.0
A(NX,NX) = 1.0
deriv_term = 0.0
do i= 0,nyfix
deriv_term = deriv_term + Gbeta_wt(nyfix-i) * xbh(y(i+1), tn)
enddo
rhs(NX) = xbh(y(nyfix), tn) - yterm_wt * e(x_high, y(nyfix),tn) * deriv_term
! 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
! Then a forward sweep of the triangularized system to get the "intermediate" solution
u_int(0,nyfix) = rhs(0) !u_int(0,nyfix) = low Dirichlet boundary condition
do i = 1, NX-1
temp = 0.0
do j = 0,i-1
temp = temp + A(i,j) * u_int(j,nyfix)
enddo !done j
u_int(i,nyfix) = (rhs(i) - temp) /A(i,i)
enddo !done i
u_int(NX,nyfix) = rhs(NX) !u_int(NX,nyfix) = Dirichlet boundary condition
!
!-=-=-=-=-
!
! Solution u(x, y ,t) is now available at this fixed y-point (nyfix).
! Now increment nyfix and repeat to get next y-slice
enddo !nyfix
! update the solution, and then alternate the direction:
do i = 0,NX
do j = 1, NY-1
u(i,j) = u_int(i,j) !Note u(0/NX,0/NY) (u at four corners) are never used/solved for/updated in this method
enddo
u(i,0) = ybl(x(i),tn) !At the y boundaries, these values are not solved/are also not needed weep
u(i,NY) = ybh(x(i),tn)
enddo
!YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY
! nstep = nstep + 1 !advance the timestep for the alternate direction
! tn = T_begin + nstep * del_t
do nxfix = 1, NX-1
! Set up matrix equations for the finite differences: Ay U_i^n+1 = U_i^n = rhs_i ..
! First the NY x NY coeffs for the super-triangular matrix in Y-direction
do i = 1 , NY-1 !compute matix coeffs for row i (y-direction implicit)
! x_diffuse_wt = d(x(nxfix), y(i),tn) * xterm_wt !not sure if best at t_n or t_{n-1}
y_diffuse_wt = e(x(nxfix), y(i),tn) * yterm_wt
do k = 0, i+1 !column k
A(i,k) = - y_diffuse_wt * Gbeta_wt(i-k+1)
enddo
A(i,i) = A(i,i) + 1.0
! and the right hand side of this equation
rhs(i) = u(nxfix,i)
enddo
!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) = ybl(x(nxfix),tn) !expecting a zero boundary condition
A(NY,NY-1) = 0.0
A(NY,NY) = 1.0
rhs(NY) = ybh(x(nxfix), tn)
! Solve the resulting matrix system - sum of a lower triangular & a super-diagnoal matrices
! First do a backward sweep to make it a lower triangular system
! Subtract a multiple of the last row with two non-zero entries
ratio = A(NY-1,NY) / A(NY,NY)
A(NY-1,NY-1) = A(NY-1,NY-1)- ratio * A(NY,NY-1)
rhs(NY-1) = rhs(NY-1) - ratio * rhs(NY)
! Remaining rows are not as sparse as the last
do i = NY-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
! Then a forward sweep to get the "intermediate" solution
u_int(nxfix, 0) = rhs(0) !u_int(0,nyfix) = low Dirichlet boundary condition = 0
do i = 1, NY-1
temp = 0.0
do j = 0,i-1
temp = temp + A(i,j) * u_int(nxfix,j)
enddo !done j
u_int(nxfix,i) = (rhs(i) - temp) /A(i,i)
enddo !done i
u_int(nxfix,NY) = rhs(NY) !Note u_int(nxfix,NY) = high Dirichlet boundary condition
! Solution u(x, y ,t) now available at this fixed x-point (nxfix).
enddo !nxfix: increment nxfix and repeat to get next x-slice
! update the solution and then repeat the big "do while in the sky"
do j = 0,NY
do i = 1, NX-1
u(i,j) = u_int(i,j)
enddo
u(0,j) = xbl(y(j),tn) !At the x boundaries, these values are not updated in this sweep
u(NX,j) = xbh(y(j),tn)
enddo
enddo !do while nstep
!Print the final results
!*******Begin USER NOTE 2
! NOTE***_ Since for this test case the exact solution is known, the following
! code also computes the absolute solution error and writes the results to the output file.
! For general problems, where the exact solution is not known, the part of the code below that
! computes and prints the exact solution and the error should be deleted,
! That is delete code for: u_exact, sln_err, abs_sln_err)
100 Format (5(F12.4,3x))
write(UNIT=20,FMT=*) "x(i), y(j), u_exact, u(i,j), sln_err,T_end=", tn
abs_sln_err = 0.0
do i = 1 ,NX-1
do j = 1, NY-1 !do j = 0, NY
u_exact = exp(-tn) * (x(i)**3) * (y(j)**3.6)
sln_err = abs(u_exact - u(i,j))
abs_sln_err = max(sln_err, abs_sln_err)
Write(UNIT=20,FMT = 100) x(i) , y(j) , u_exact, u(i,j), sln_err
write(*,*) "/"
enddo
enddo
Write(UNIT=20,FMT = *) "Absolute error = ", abs_sln_err
!*******End USER NOTE 2
Close(UNIT=20, IOSTAT= checkstatus)
write(*,*) "File Open IOSTAT =", checkstatus
END Program Two_D_ADI_FPDE
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function d(x, y, t)
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: x, y, t, d
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha, beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6, & alpha = 1.8, beta = 1.6
d = 0.18363375 * (x**2.8) * y !gamma(2.2)/6 = 1.10180249088/6 =0.18363375
End Function d
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function e(x, y, t)
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: x, y, t, e
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6, & alpha = 1.8, beta = 1.6
e = 0.1494624672 * x * (y**2.6) !2/gamma(4.6) = 2/13.3812858709 =
End Function e
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function q(x, y, t) !forcing function in the diff eqn, v used/added to set up test cases
real :: q
real :: x_low, y_low !low limit on x and y
real :: x_high, y_high !high limit on x and y
real :: alpha, beta !fractional derivative order for the x and y diffusive terms
integer :: NX, NY !number of subintervals in x and y directions
real :: x, y, t, q
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6; & alpha = 1.8, beta = 1.6
q = -(1.0 + 2*x*y) * exp(-t) * (x**3) * (y**3.6)
End Function q
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function f(x, y) !initial condition for u(x,y,t) at t=0
real :: f
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: x, y, f
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6; & alpha = 1.8, beta = 1.6
f = (x**3) * (y**3.6)
End Function f
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function xbl(y, t) !Dirichlet boundary condition at x_low ; u(x=x_low,y, t)
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: y, t, xbl
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6; & alpha = 1.8, beta = 1.6
xbl = 0.0 ! zero boundary function at x = x_low = 0
End Function xbl
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function xbh(y, t) !Dirichlet boundary condition at x_high: u(x=x_high,y, t)
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: y, t, xbh
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6; & alpha = 1.8, beta = 1.6
xbh = exp(-t)* (x_high**3) * (y**3.6)
End Function xbh
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function ybl(x, t) !lDirichlet boundary condition at y_low: u(x,y=y_low, t)
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: x, t, ybl
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6; & alpha = 1.8, beta = 1.6
ybl = 0.0 ! zero boundary function at y= y_low = 0
End Function ybl
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function ybh(x, t) !Dirichlet boundary condition at y_high: u(x,y=y_high, t)
real :: x_low !low limit on x
real :: x_high !high limit on x
real :: y_low !low limit on y
real :: y_high !high limit on y
real :: alpha !fractional derivative order for the x- diffusive term
real :: beta !fractional derivative order for the y- diffusive term
integer :: NX,NY !number of subintervals in x and y directions
real :: x, t, ybh
COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta
! TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6; & alpha = 1.8, beta = 1.6
ybh = exp(-t)* (x**3) * (y_high**3.6)
End Function ybh