子程序有问题请大家帮忙看下

小小小 · 发表于 2022-11-29 16:02:45

在运行时，显示Solver_Crank-Nicolson (coef1,coef2,f,u)这个子程序每一行都显示有问题，我把代码都附在下面了，麻烦大家帮忙运行看下，非常感谢！

用到的文本信息nml-pata-DBC.txt

&nml_para
nx = 100,
nt = 100,
xmin  = 0.d0,
xmax  = 1.d0,
tmin  = 0.d0,
tmax  = 1.d0,
k_max =5 ,

/

[Fortran] 纯文本查看 复制代码

001module mod_fdm_tool_C3
002    use iso_fortran_env
003    implicit none
004    real(real64), parameter :: pi = 4.d0*datan(1.d0)
005contains
006 
007subroutine Grid_2D_Type(nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)
008    implicit none
009    integer(int32), intent(in)  :: nx, nt
010    real(real64),   intent(in)  :: xmin, xmax, tmin, tmax
011    real(real64),   intent(out) :: hx, ht, x(0:nx), t(0:nt)
012    integer(int32) :: i, j
013 
014    hx = (xmax - xmin) / dble(nx)
015    ht = (tmax - tmin) / dble(nt)
016 
017    do concurrent (i = 0 : nx)
018        x(i) = xmin + dble(i) * hx
019    end do
020    x(nx) = xmax
021 
022    do concurrent (j = 0 : nt)
023        t(j) = tmin + dble(j) * ht
024    end do
025    t(nt) = tmax
026 
027    return
028end subroutine Grid_2D_Type
029 
030 
031 
032!===============================================================================
033pure elemental function f_rhs(x,t) result(res)
034    implicit none
035    real(real64), intent(in) :: x, t
036    real(real64) :: res
037 
038    res = 0
039 
040    return
041end function f_rhs
042 
043!===============================================================================
044pure elemental function Solu_Exact(x,t) result(res)
045    implicit none
046    real(real64), intent(in) :: x, t
047    real(real64) :: res
048 
049    res = dexp(x+t)
050 
051    return
052end function Solu_Exact
053 
054 
055subroutine Apply_DBC(x, t,u)
056    implicit none
057    real(real64), intent(in)  :: x(:), t(:)
058    real(real64), intent(out) :: u(:,:)
059    integer(int32) :: nt
060 
061    nt = size(t)
062 
063    u(:, 1) = dexp(x)
064    u(1,:) = dexp(t)
065    u(nt,:) = dexp(1.0+t)
066 
067 
068 
069    return
070end subroutine Apply_DBC
071 
072!======================================================
073subroutine Gene_Alg_Sys_DBC(ht,hx,nx,r,coef1,coef2)
074    ! Algebraic system, i.e., Au = f
075    implicit none
076    integer(int32) ,intent(in):: nx
077    real(real64), intent(in)  :: hx,ht
078    real(real64), intent(out) :: coef1(nx-1,3),coef2(nx-1,3),r
079 
080    r = ht/hx/hx
081 
082 
083    ! all the following variable are arraies, not scalars
084    coef1(:,1) = -r/2.d0
085    coef1(:,2) = 1+r
086    coef1(:,3) = -r/2.d0
087 
088    coef2(:,1) = r/2.d0
089    coef2(:,2) = 1-r
090    coef2(:,3) = r/2.d0
091 
092 
093 
094    coef1(1,1) = 0.d0
095    coef1(nx-1,3) = 0.d0
096 
097    coef2(1,1) = 0.d0
098    coef2(nx-1,3) = 0.d0
099 
100 
101 
102    return
103end subroutine Gene_Alg_Sys_DBC
104 
105 
106!这个子程序有问题======================================================
107subroutine Solver_Crank-Nicolson (coef1,coef2,f,u)
108 
109    implicit none
110    integer(int32) :: nx,nt,i,j
111    real(real64), intent(in)  :: coef1(:,:),coef2(:,:)
112    real(real64) :: f(:)
113    real(real64) :: u(:,:)
114    real(real64) ::r
115    real(real64),allocatable::g(:),w(:)
116 
117 
118    do j = 1, nt-1
119      f(1)=coef2(1,2)*u(1,j)
120      f(nx-1)=coef2(nx-1,1)*u(nx-2,j)+coef2(nx-1,2)*u(nx-1,j)
121      f(2:nx-1)=coef2(2:nx-2,1)*u(1:nx-3,j)+coef2(2:nx-2,2)*u（2:nx-2,j)+coef2(2,nx-2,3)*u(3:nx-1,j)
122 
123      allocate(g(nx-1),w(nx-1))
124 
125      g(1)=f(1)/coef1(1,2)
126      w(1)=coef1(1,3)/coef1(1,2)
127 
128      do i=2:nx-1
129        w(i)=coef1(i,2)-coef1(i,1)*w(i-1)
130        g(i)=(f(i)-coef1(i,1)*g(i-1))/w(i)
131        w(i)=coef1(i,3)/w(i)
132      end do
133      u(nx-1,:)=g(nx-1)
134 
135      do i=nx-2,1,-1
136        u(i,j+1)=g(i)-w(i)*u(i+1,j+1)
137      end do
138    end do
139 
140      deallocate(g,w)
141 
142    return
143end subroutine Solver_Crank-Nicolson
144 
145end module mod_fdm_tool_C3
146 
147!=============================================================
148! cLi, 2022-11-24
149!
150! Solving example 3.1.1 in p.67:
151!
152!   -u'(t) - u''(x) = 0,  x \in (0, 1),t \in(0,1]
153!   u(x,0) =  exp(x),   x \in (0, 1)
154!   u(0,t) = exp(t) ,u(1,t) = exp(1+t) ,t \in(0,1]
155!
156!
157!   
158!   
159!   u(x,t) = exp(x+t)
160!=============================================================
161program Crank-Nicolson_2D_DBC
162 
163    use iso_fortran_env
164    use mod_fdm_tool_C3
165 
166implicit none
167 
168integer(int32) :: nx, nt, m, k, k_max, fileunit,  i,j
169real(real64)   :: hx, xmin, xmax
170real(real64)   :: ht, tmin, tmax
171 
172real(real64)   :: iter_tol, iter_fac, error = 0.d0
173real(real64), allocatable :: x(:), t(:), f(:,:), u(:,:), v(:,:)
174character(len=:), allocatable :: dirname, filename, output_fmt
175logical :: istat, debug = .false.
176 
177namelist /nml_para/ nx, nt, xmin, xmax, tmin, tmax, k_max
178 
179print '(1X,A)', 'Example_311'
180 
181filename = 'nml_para_DBC.txt'
182open(newunit=fileunit, file=filename)
183read(unit=fileunit, nml=nml_para)
184close(fileunit)
185 
186dirname = 'Out'
187inquire(file=dirname, exist=istat)
188if (.not. istat) then
189    call Make_Directory(dirname)
190end if
191 
192 
193filename = dirname//'/results.txt'
194open(newunit=fileunit, file=trim(filename))
195write(fileunit, '(a)') 'Inf error at x = 0.5, 0.5,  0.5, 0.5,  0.5,'
196write(fileunit, '(a)') '         and t = 0.1, 0.2,  0.3, 0.4,  0.5 are'
197 
198do k=1,k_max
199   m  = 2**(k-1)
200    nx = 10*m
201    nt = 10*m
202    print '(1X, A, I3, A)', 'The spatial step size of 1/', nx, ',the time step size of 1/',nt,' is running...'
203 
204    if (allocated(u)) deallocate(x, t, f, u, v)
205    allocate(x(0:nx), t(0:nt), f(0:nx,0:nt), u(0:nx,0:nt), v(0:nx,0:nt))
206 
207     ! start the FDM solving procedure
208    call Grid_2D_Type (nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)
209 
210    call Apply_DBC(x, t, u)
211 
212    call Gene_Alg_Sys_DBC(ht,hx,nx,coef1,coef2)
213 
214    do concurrent (i = 0 : nx)
215 
216        v(i,:) = Solu_Exact(x(i),t)
217    end do
218 
219 
220    call Solver_Crank-Nicolson (coef1,coef2,f,u)
221 
222     output_fmt = '(A, I3, 2X, 6(ES9.3, 2X), F5.2)'
223    write(fileunit, output_fmt) '1/', nx,'1/', nt,       &
224 
225            abs(u(nx/2*1,nt/10*1)-v(nx/2*1,nt/10*1), &
226 
227            abs(u(nx/2*1,nt/10*2)-v(u(nx/2*1,nt/10*2), &
228            abs(u(nx/2*1,nt/10*3)-v(u(nx/2*1,nt/10*3), &
229            abs(u(nx/2*1,nt/10*4)-v(nx/2*1,nt/10*4), &
230            abs(u(nx/2*1,nt/10*5)-v(nx/2*1,nt/10*5), &
231            maxval(abs(u-v)), error/maxval(abs(u-v))
232 
233    error = maxval(abs(u-v))
234end do
235 
236close(fileunit)
237print *, ''
238print *, 'The program is done.'
239 
240return
241end program Crank-Nicolson_2D_DBC