[Fortran] 纯文本查看 复制代码
子程序：
module mod_fdm_tool_C3
use iso_fortran_env
implicit none
real(real64), parameter :: pi = 4.d0*datan(1.d0)
contains
subroutine Grid_2D_Type(nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)
implicit none
integer(int32), intent(in) :: nx, nt
real(real64), intent(in) :: xmin, xmax, tmin, tmax
real(real64), intent(out) :: hx, ht, x(0:nx), t(0:nt)
integer(int32) :: i, j
hx = (xmax - xmin) / dble(nx)
ht = (tmax - tmin) / dble(nt)
do concurrent (i = 0 : nx)
x(i) = xmin + dble(i) * hx
end do
x(nx) = xmax
do concurrent (j = 0 : nt)
t(j) = tmin + dble(j) * ht
end do
t(nt) = tmax
return
end subroutine Grid_2D_Type
!===============================================================================
pure elemental function f_rhs(x,t) result(res)
implicit none
real(real64), intent(in) :: x, t
real(real64) :: res
res = 0
return
end function f_rhs
!===============================================================================
pure elemental function Solu_Exact(x,t) result(res)
implicit none
real(real64), intent(in) :: x, t
real(real64) :: res
res = dexp(x+t)
return
end function Solu_Exact
subroutine Apply_DBC(x, t,u)
implicit none
real(real64), intent(in) :: x(:), t(:)
real(real64), intent(out) :: u(:,:)
integer(int32) :: nt
nt = size(t)
u(:, 1) = dexp(x)
u(1,:) = dexp(t)
u(nt,:) = dexp(1.0+t)
return
end subroutine Apply_DBC
!======================================================
subroutine Gene_Alg_Sys_DBC(ht,hx,nx,r,coef1,coef2)
! Algebraic system, i.e., Au = f
implicit none
integer(int32) ,intent(in):: nx
real(real64), intent(in) :: hx,ht
real(real64), intent(out) :: coef1(nx-2,3),coef2(nx-2,3),r
r = ht/hx/hx
! all the following variable are arraies, not scalars
coef1(:,1) = -r/2.d0
coef1(:,2) = 1+r
coef1(:,3) = -r/2.d0
coef2(:,1) = r/2.d0
coef2(:,2) = 1-r
coef2(:,3) = r/2.d0
coef1(1,1) = 0.d0
coef1(nx-2,3) = 0.d0
coef2(1,1) = 0.d0
coef2(nx-2,3) = 0.d0
return
end subroutine Gene_Alg_Sys_DBC
!======================================================
subroutine Solver_Crank_Nicolson (coef1,coef2,L,u)
implicit none
integer(int32) :: nx,nt,i,j
real(real64), intent(in) :: coef1(:,:),coef2(:,:)
real(real64) :: L(:),r
real(real64) :: u(:,:)
real(real64),allocatable::g(:),w(:)
do j = 1, nt-1
L(1)=coef2(1,2)*u(2,j)+coef2(1,3)*u(3,j)+r*(u(1,j)+u(1,j+1))/2
L(nx-2)=coef2(nx-2,1)*u(nx-2,j)+coef2(nx-2,2)*u(nx-1,j)+r*(u(nx,j)+u(nx,j+1))/2
L(2:nx-3)=coef2(2:nx-3,1)*u(2:nx-3,j)+coef2(2:nx-3,2)*u(3:nx-2,j)+coef2(2:nx-3,3)*u(4:nx-1,j)
allocate(g(nx),w(nx))
g(1)=L(1)/coef1(1,2)
w(1)=coef1(1,3)/coef1(1,2)
do i=2,nx-2
w(i)=coef1(i,2)-coef1(i,1)*w(i-1)
g(i)=(L(i)-coef1(i,1)*g(i-1))/w(i)
w(i)=coef1(i,3)/w(i)
end do
u(nx,:)=g(nx)
do i=nx-3,2,-1
u(i,j+1)=g(i)-w(i)*u(i+1,j+1)
end do
end do
deallocate(g,w)
return
end subroutine Solver_Crank_Nicolson
subroutine Make_Directory(dirname)
implicit none
character(len=:), allocatable, intent(in) :: dirname
call system("mkdir "//trim(AdjustL(dirname)))
end subroutine Make_Directory
end module mod_fdm_tool_C3
主程序：
!=============================================================
! cLi, 2022-11-24
!
! Solving example 3.1.1 in p.67:
!
! -u'(t) - u''(x) = 0, x \in (0, 1),t \in(0,1]
! u(x,0) = exp(x), x \in (0, 1)
! u(0,t) = exp(t) ,u(1,t) = exp(1+t) ,t \in(0,1]
!
!
!
!
! u(x,t) = exp(x+t)
!=============================================================
program Crank_Nicolson_2D_DBC
use iso_fortran_env
use mod_fdm_tool_C3
implicit none
integer(int32):: nx, nt, m, k, k_max, fileunit, i,j
real(real64) :: hx, xmin, xmax
real(real64) :: ht, tmin, tmax
real(real64) :: r,iter_tol, iter_fac, error = 0.d0
real(real64), allocatable :: x(:), t(:),L(:), u(:,:), v(:,:),coef1(:,:),coef2(:,:)
character(len=:), allocatable :: dirname, filename, output_fmt
logical :: istat, debug = .false.
namelist /nml_para/ nx, nt, xmin, xmax, tmin, tmax, k_max
print '(1X,A)', 'Example_311'
filename = 'nml_para_DBC.txt'
open(newunit=fileunit, file=filename)
read(unit=fileunit, nml=nml_para)
close(fileunit)
dirname = 'Out'
inquire(file=dirname, exist=istat)
if (.not. istat) then
call Make_Directory(dirname)
end if
filename = dirname//'/results.txt'
open(newunit=fileunit, file=trim(filename))
write(fileunit, '(a)') 'Inf error at x = 0.5, 0.5, 0.5, 0.5, 0.5,'
write(fileunit, '(a)') ' and t = 0.1, 0.2, 0.3, 0.4, 0.5 are'
do k=1,k_max
m = 2**(k-1)
nx = 10*m
nt = 10*m
print '(1X, A, I3, A,I3, A)', 'The spatial step size of 1/', nx, ',the time step size of 1/',nt,' is running...'
if (allocated(u)) deallocate(x, t, u, v)
allocate(x(0:nx), t(0:nt), u(0:nx,0:nt), v(0:nx,0:nt))
! start the FDM solving procedure
call Grid_2D_Type (nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)
call Apply_DBC(x, t, u)
call Gene_Alg_Sys_DBC(ht,hx,nx,r,coef1,coef2)
do concurrent (i = 0 : nx)
v(i,:) = Solu_Exact(x(i),t)
end do
call Solver_Crank_Nicolson (coef1,coef2,L,u)
output_fmt = '(A, I3, 2X, 6(ES9.3, 2X), F5.2)'
write(fileunit, output_fmt) '1/', nx,'1/', nt, &
abs(u(nx/2*1,nt/10*1)-v(nx/2*1,nt/10*1)), &
abs(u(nx/2*1,nt/10*2)-v(nx/2*1,nt/10*2)), &
abs(u(nx/2*1,nt/10*3)-v(nx/2*1,nt/10*3)), &
abs(u(nx/2*1,nt/10*4)-v(nx/2*1,nt/10*4)), &
abs(u(nx/2*1,nt/10*5)-v(nx/2*1,nt/10*5)), &
maxval(abs(u-v)), error/maxval(abs(u-v))
error = maxval(abs(u-v))
end do
close(fileunit)
print *, ''
print *, 'The program is done.'
return
end program Crank_Nicolson_2D_DBC