! P-1D-UCDR-HOC-2021+Pro1.f90
!------------------------------------------------------------------------------------
! The progrom for 1-D Unsteady Nonlinear Advection Diffusion Reaction Equation using Thomas method
!--------------------------------------------------------------
! Four-order accurcy in both time and space with full implicit scheme
!--------------------------------------------------------------
! Lin Zhang and Yongbin Ge, Numerical solution of nonlinear advection
! diffusion reaction equation using high-order compact difference method,
! Applied Numerical Mathematics,2021, 166: 127–145.
!--------------------------------------------------------------
! Model equation:
! ut+u*ux=eps*uxx
!-------------------------------------------------------------------------------------------------------------------------------------
!$debug
! use PORTLIB
implicit double precision (a-h,o-z)
parameter(it=100000)
common/vari/u(it),u1(it),ub(it),u2(it),u3(it),u4(it),fu(it),ue(it),ux3(it),uxx3(it),uxx(it),ux(it),Rx(it),Rxx(it),Rx1(it),Rxx1(it),Rx3(it),Rxx3(it),Rx2(it),Rxx2(it),uxx1(it),ux1(it),uxx2(it),ux2(it)
common/const1/N,itfa,eps,Rmin,Rmax
common/const2/ni1
character(len=1) endchar
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----default parameters
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
!------------------default parameters------------------------------------------------
write(*,*)'***************************************************'
write(*,*)'Definition:'
write(*,*)' N: Number of grid subdomain in space direction'
write(*,*)' Nnode: Number of grid points'
write(*,*)' h: =1/N, space step-length'
write(*,*)' T: Total time for computation to be reached'
write(*,*)' dt: Time step-length'
write(*,*)' Tstep: =T/dt,total time steps to reach T'
write(*,*)' v: Diffusion coefficient'
write(*,*)' u: Numerical solution'
write(*,*)' errmax: Maximum error on the final time step'
write(*,*)' errL2: L2-norm error on the final time step'
write(*,*)' errRMS: RMS error on the final time step'
write(*,*)' erraver: Average error on the final time step'
! write(*,*)' CPUtime: CPUtime exclude input and output time'
write(*,*)'***************************************************'
open(5,file='output.dat',status='unknown')
write(*,*) ' '
write(*,*)'Please input parameters as follows:'
!-------------------------------------------------------------------------------------------------------------------------------------
1000 continue
write(*,*) 'N= T= '
read(*,*) N,T
!-------------------------------------------------------------------------------------------------------------------------------------
write(*,*)'Display the computing process? input "y" to be,or "n" not to:'
read(*,*) endchar
IF (endchar.GE.'a'.AND.endchar.LE.'z') then
endchar = CHAR(ICHAR(endchar) -ICHAR('a') +ICHAR('A'));
end if;
if(endchar.eq.'Y') then
write(*,30)
end if
!-------------------------------------------------------------------------------------------------------------------------------------
CPUtime1=TIMEF()
pi=4.d0*datan(1.d0)
pii=pi*pi
Rmin=0.d0
Rmax=1.d0
h=(Rmax-Rmin)/float(N)
tol=1.0e-13
eps=0.001d0 !v
dt=0.01
itfas=idnint(T/dt)
!-------------------------------------------------------------------------------------------------------------------------------------
!----- grid numbers for all levles---------
!
ni1=N+1
!-----initial ----------------
t=0.d0
do 26 i=1,ni1
x=Rmin+(i-1)*h
ub(i)=2.d0*eps*pi*dsin(pi*x)/(100.d0+dcos(pi*x)) !!!
26 continue
do 70 m=1,itfas
!-----assembly coefficients------------------
!*************************To Calculated value of u(i)(u(dt/2,h)) when dt=dt/2************************************
t=m*dt
u(1)=0.d0
u(ni1)=0.d0
num0=0
300 num0=num0+1
errux=0.0
erruxx=0.0
do 200 i=1,ni1
Rx(i)=ux(i)
Rxx(i)=uxx(i)
200 continue
call adi2(ux,ub,uxx,ni1,it)
call adi3(uxx,ub,ux,ni1,it)
do 201 i=1,ni1
errx=dabs(Rx(i)-ux(i))
errxx=dabs(Rxx(i)-uxx(i))
errux=dmax1(errux,errx)
erruxx=dmax1(erruxx,errxx)
201 continue
if(errux.gt.tol.and.erruxx.gt.tol) then
goto 300
end if
do 289 i=1,ni1
u1(i)=ub(i)+(1.d0/2.d0)*dt*(eps*uxx(i)-ub(i)*ux(i))
289 continue
num1=0
400 num1=num1+1
err1ux=0.0
err1uxx=0.0
do 202 i=1,ni1
Rx1(i)=ux1(i)
Rxx1(i)=uxx1(i)
202 continue
call adi2(ux1,u1,uxx1,ni1,it)
call adi3(uxx1,u1,ux1,ni1,it)
do 203 i=1,ni1
err1x=dabs(Rx1(i)-ux1(i))
err1xx=dabs(Rxx1(i)-uxx1(i))
err1ux=dmax1(err1ux,err1x)
err1uxx=dmax1(err1uxx,err1xx)
203 continue
if(err1ux.gt.tol.and.err1uxx.gt.tol) then
goto 400
end if
do 291 i=1,ni1
u2(i)=u1(i)+1.d0/2.d0*dt*(eps*uxx1(i)-u1(i)*ux1(i))
291 continue
num2=0
500 num2=num2+1
err2ux=0.0
err2uxx=0.0
do 204 i=1,ni1
Rx2(i)=ux2(i)
Rxx2(i)=uxx2(i)
204 continue
call adi2(ux2,u2,uxx2,ni1,it)
call adi3(uxx2,u2,ux2,ni1,it)
do 205 i=1,ni1
err2x=dabs(Rx2(i)-ux2(i))
err2xx=dabs(Rxx2(i)-uxx2(i))
err2ux=dmax1(err2ux,err2x)
err2uxx=dmax1(err2uxx,err2xx)
205 continue
if(err2ux.gt.tol.and.err2uxx.gt.tol) then
goto 500
end if
do 293 i=1,ni1
u3(i)=2.d0/3.d0*ub(i)+1.d0/3.d0*u2(i)+1.d0/6.d0*dt*(eps*uxx2(i)-u2(i)*ux2(i))
293 continue
num3=0
502 num3=num3+1
err3ux=0.0
err3uxx=0.0
do 402 i=1,ni1
Rx3(i)=ux3(i)
Rxx3(i)=uxx3(i)
402 continue
call adi2(ux3,u3,uxx3,ni1,it)
call adi3(uxx3,u3,ux3,ni1,it)
do 206 i=1,ni1
err3x=dabs(Rx3(i)-ux3(i))
err3xx=dabs(Rxx3(i)-uxx3(i))
err3ux=dmax1(err3ux,err3x)
err3uxx=dmax1(err3uxx,err3xx)
206 continue
if(err3ux.gt.tol.and.err3uxx.gt.tol) then
goto 502
end if
do 292 i=1,ni1
u4(i)=u3(i)+1.d0/2.d0*dt*(eps*uxx3(i)-u3(i)*ux3(i))
292 continue
do 294 i=2,ni1-1
u(i)=u4(i)
294 continue
do 295 i=1,ni1
ub(i)=u(i)
295 continue
if(endchar.eq.'Y') then
do 61 i=1,ni1
t=m*dt
x=Rmin+h*(i-1)
ue(i)=2.d0*eps*pi*dexp(-pii*eps*t)*dsin(pi*x)/(100.d0+dexp(-pii*eps*t)*dcos(pi*x))
write(*,122) t,x,u(i),ue(i),dabs(u(i)-ue(i))
61 continue
end if
70 continue
CPUtime2=TIMEF()
CPUtime=CPUtime2-CPUtime1
if(endchar.eq.'Y') then
write(*,30)
end if
!
!-----Calculate the value of various errors in the final time step-------------------
!
errmax=0.d0
errL2=0.d0
RMS=0.d0
erraver=0.d0
do 2000 i=1,ni1
x=Rmin+h*(i-1)
t=dt*itfas
ue(i)=2.d0*eps*pi*dexp(-pii*eps*t)*dsin(pi*x)/(100.d0+dexp(-pii*eps*t)*dcos(pi*x))
errf=dabs(u(i)-ue(i))
erraver=erraver+errf
errmax=dmax1(errmax,errf)
errL2=errL2+errf*errf
RMS=RMS+errf*errf
2000 continue
write(*,*) ue(2),u(2),u(2)-ue(2)
errL2=dsqrt(h*errL2)
RMS=dsqrt(RMS/(N+1))
erraver=erraver/(N+1)
!
!-----keep stata for Tecplot to draw pictures------------
!
open(1,file='fig_u_ue_err.dat',status='unknown')
do 2004 i=1,ni1
x=Rmin+(i-1)*h
write(1,222) x,u(i),ue(i),dabs(u(i)-ue(i))
2004 continue
close(1)
222 format(3x,f14.8,3(3x,e14.8))
open(2,file='fig_u.dat',status='unknown')
write(2,*) 'variable=x u'
write(2,*) 'Zone T=grids F=point I=',ni1
do 2005 i=1,ni1
x=Rmin+(i-1)*h
write(2,*) x,u(i)
2005 continue
close(2)
open(3,file='fig_ue.dat',status='unknown')
write(3,*) 'variable=x ue'
write(3,*) 'Zone T=grids F=point I=',ni1
do 3006 i=1,ni1
x=Rmin+(i-1)*h
write(3,*) x,ue(i)
3006 continue
close(3)
open(4,file='fig_error.dat',status='unknown')
write(4,*) 'variable=x error'
write(4,*) 'Zone T=grids F=point I=',ni1
do 2007 i=1,ni1
x=Rmin+(i-1)*h
write(4,*) x,dabs(u(i)-ue(i))
2007 continue
close(4)
!*****************************Print input and default parameters and Print ouput and computed results******************
open(5,file='output.dat',status='old')
write(5,*)'****************************************************'
write(*,*)'*************Print input and default parameters*****'
write(*,*)'N= ',N
write(5,*)'N= ',N
write(5,*)'Nnode= ',N+1
write(*,*)'Nnode= ',N+1
write(*,*)'h= ',h
write(5,*)'h= ',h
write(*,*)'eps= ',eps
write(5,*)'eps= ',eps
write(*,*)'tol= ',tol
write(5,*)'tol= ',tol
write(*,*)'T= ',T
write(5,*)'T= ',T
write(*,*)'itfas= ',itfas
write(5,*)'itfas= ',itfas
write(*,*)'**************Print ouput and computed results*****'
write(*,*)'errmax= ',errmax
write(5,*)'errmax= ',errmax
write(*,*)'errL2= ',errL2
write(5,*)'errL2= ',errL2
write(*,*)'errRMS= ',RMS
write(5,*)'errRMS= ',RMS
write(*,*)'erraver= ',erraver
write(5,*)'erraver= ',erraver
write(*,*)'CPUtime= ',CPUtime
write(5,*)'CPUtime= ',CPUtime
write(*,*)'continued? input y to be,or n to stop:'
read(*,*) endchar
if (endchar.GE.'a'.AND.endchar.LE.'z') then
endchar = CHAR(ICHAR(endchar) -ICHAR('a') +ICHAR('A'));
end if;
if(endchar.eq.'Y') then
goto 1000
else
goto 2001
end if
2001 write(*,*)'Please input <0-9> to return!'
read(*,*) endchar
30 format (7x,'t',11x,'x',10x,'Numerical',12x,'Exact',14x,'Error')
122 format(3x,f8.5,3x,f8.5,3x,e16.8,3x,e16.8,3x,e16.8)
stop
end
!**************************************************************
! subroutine adi1(uxx,ub,ni,it)
! implicit double precision(a-h,o-z)
! dimension ub(it),uxx(it)
! dimension a(100000),b(100000),c(100000),d(100000)
! pi=4.d0*datan(1.d0)
! h=1.d0/float(ni)
! h2=h*h
! a(1)=1.d0
! b(1)=51.d0/52.d0
! c(1)=0.d0
! d(1)=12293.d0*ub(1)/(2340.d0*h2)-18903.d0*ub(2)/(1040.d0*h2)+2891.d0*ub(3)/(104.d0*h2)-23941.d0*ub(4)/(936.d0*h2)+387.d0*ub(5)/(26.d0*h2)-5063.d0*ub(6)/(1040.d0*h2)+494.d0*ub(7)/(720.d0*h2)
! do 2023 i=2,ni-1
! a(i)=5.d0/6.d0
! b(i)=1.d0/12.d0
! c(i)=1.d0/12.d0
! d(i)=(-2.d0*ub(i)+ub(i+1)+ub(i-1))/h2
!2023 continue
! a(ni)=1.d0
! b(ni)=0.d0
! c(ni)=-51.d0/52.d0
! d(ni)=17599.d0*ub(ni)/(4680.d0*h2)-17237.d0*ub(ni-1)/(1040.d0*h2)+795.d0*ub(ni-2)/(26.d0*h2)-28735.d0*ub(ni-3)/(936.d0*h2)+1871.d0*ub(ni-4)/(104.d0*h2)-6117.d0*ub(ni-5)/(1040.d0*h2)+596.d0*ub(ni-6)/(720.d0*h2)
! call trid(1,ni,a,b,c,d)
! do 2024 i=1,ni
! uxx(i)=d(i)
!2024 continue
! return
! end
subroutine adi2(ux,ub,uxx,ni,it)
implicit double precision(a-h,o-z)
dimension ux(it),ub(it),uxx(it)
dimension a(100000),b(100000),c(100000),d(100000),f(100000),e(100000)
common/const1/N,itfas,eps,Rmin,Rmax
h=(Rmax-Rmin)/float(ni-1)
d(1)=(1.d0/h)*(-2218.d0*ub(1)/735.d0+72001.d0*ub(2)/12600.d0-95.d0*ub(3)/28.d0+1193.d0*ub(4)/126.d0-919.d0*ub(5)/126.d0+2981.d0*ub(6)/840.d0-6263.d0*ub(7)/6300.d0+6.d0*ub(8)/49.d0)
d(ni+1)=(1.d0/h)*(134.d0*ub(ni+1)/245.d0-25153.d0*ub(ni)/4200.d0+877.d0*ub(ni-1)/140.d0-50.d0*ub(ni-2)/20.d0-161.d0*ub(ni-3)/126.d0+571.d0*ub(ni-4)/280.d0-2011.d0*ub(ni-5)/2100.d0-727.d0*ub(ni-6)/4410.d0)
d(ni)=(1.d0/h)*(-5803.d0*ub(ni+1)/2520.d0-226.d0*ub(ni)/45.d0+3521.d0*ub(ni-1)/720.d0-31.d0*ub(ni-2)/18.d0-121.d0*ub(ni-3)/72.d0+84.d0*ub(ni-4)/45.d0-607.d0*ub(ni-5)/720.d0-91.d0*ub(ni-6)/630.d0)
d(2)=(1.d0/h)*(-6541.d0*0.d0+211.d0*ub(2)/30.d0-1505.d0*ub(3)/144.d0+104.d0*ub(4)/9.d0-209.d0*ub(5)/24.d0-377.d0*ub(6)/90.d0-839.d0*ub(7)/720.d0+1.d0*ub(8)/7.d0)
a(1)=1.d0
b(1)=187.d0/210.d0
c(1)=0.d0
f(1)=0.d0
e(1)=0.d0
a(2)=0.d0
b(2)=-1.d0/12.d0
c(2)=1.d0
f(2)=0.d0
e(2)=0.d0
do 20001 i=3,ni-1
a(i)=1.d0
b(i)=1.d0/3.d0
c(i)=1.d0/3.d0
f(i)=0.d0
e(i)=0.d0
d(i)=1.d0*ub(i+2)/(36.d0*h)-1.d0*ub(i-2)/(36.d0*h)+14.d0*ub(i+1)/(18.d0*h)-14.d0*ub(i-1)/(18.d0*h)
20001 continue
a(ni+1)=1.d0
c(ni+1)=-187.d0/210.d0
b(ni+1)=0.d0
f(ni+1)=0.d0
e(ni+1)=0.d0
a(ni)=0.d0
c(ni)=1.d0/12.d0
b(ni)=1.d0
f(ni)=0.d0
e(ni)=0.d0
call pentad(1,ni,e,f,a,b,c,d)
do 20004 i=1,ni
ux(i)=d(i)
20004 continue
return
end
!---------------------------------------------------------------------------------------------------------------------
subroutine adi3(uxx,ub,ux,ni,it)
implicit double precision(a-h,o-z)
dimension ux(it),ub(it),uxx(it)
dimension a(100000),b(100000),c(100000),d(100000),f(100000),e(100000)
common/const1/N,itfas,eps,Rmin,Rmax
h=(Rmax-Rmin)/float(ni-1)
h2=h*h
d(1)=(1.d0/h2)*(64447.d0*ub(1)/11088.d0-13401973.d0*ub(2)/498960.d0+479153.d0*ub(3)/7920.d0-1359929.d0*ub(4)/15840.d0+1176541.d0*ub(5)/14256.d0-211403.d0*ub(6)/3960.d0+176807.d0*ub(7)/7920.d0-495793.d0*ub(8)/90720.d0+751.d0*ub(9)/1260.d0)
d(ni+1)=(1.d0/h2)*(2583661.d0*ub(ni+1)/443520.d0-18855749.d0*ub(ni)/686070.d0+328961.d0*ub(ni-1)/5280.d0-1416017.d0*ub(ni-2)/15840.d0+4945481.d0*ub(ni-3)/57024.d0-149477.d0*ub(ni-4)/2640.d0+378583.d0*ub(ni-5)/15840.d0-5894719.d0*ub(ni-6)/997920.d0+96167.d0*ub(ni-7)/147840.d0)
d(2)=(1.d0/h2)*(140531.d0*ub(1)/23688.d0-18529.d0*ub(2)/658.d0+532891.d0*ub(3)/8460.d0-125317.d0*ub(4)/1410.d0+144553.d0*ub(5)/1692.d0-235567.d0*ub(6)/4230.d0+349969.d0*ub(7)/8460.d0-171791.d0*ub(8)/29610.d0+1609.d0*ub(9)/2520.d0)
d(ni)=(1.d0/h2)*(645671.d0*ub(ni+1)/108570.d0-8750519.d0*ub(ni)/325710.d0+5683279.d0*ub(ni-1)/93060.d0-827195.d0*ub(ni-2)/9306.d0+153419.d0*ub(ni-3)/282.d0-273665.d0*ub(ni-4)/9306.d0+2232799.d0*ub(ni-5)/93060.d0-1931399.d0*ub(ni-6)/325710.d0+905.d0*ub(ni-7)/1386.d0)
a(1)=1.d0
b(1)=41.d0/792.d0
c(1)=0.d0
f(1)=0.d0
e(1)=0.d0
a(2)=0.d0
b(2)=-26.d0/47.d0
c(2)=1.d0
f(2)=0.d0
e(2)=0.d0
do 20001 i=3,ni-1
a(i)=1.d0
b(i)=2.d0/11.d0
c(i)=2.d0/11.d0
f(i)=0.d0
e(i)=0.d0
d(i)=3.d0*(ub(i+2)-2.d0*ub(i)+ub(i-2))/(44.d0*h2)+12.d0*(ub(i+1)-2*ub(i)+ub(i-1))/(11.d0*h2)
20001 continue
a(ni+1)=1.d0
c(ni+1)=-41.d0/729.d0
b(ni+1)=0.d0
f(ni+1)=0.d0
e(ni+1)=0.d0
a(ni)=0.d0
c(ni)=26.d0/47.d0
b(ni)=1.d0
f(ni)=0.d0
e(ni)=0.d0
call pentad(1,ni,e,f,a,b,c,d)
do 20004 i=1,ni
uxx(i)=d(i)
20004 continue
return
end
!-----------------------------------------------------------------
subroutine pentad(n0,n,a,b,c,d)
parameter(iu=100000)
implicit real*8(a-h,o-z)
dimension a(iu),b(iu),c(iu),d(iu),f(iu),e(iu)
dimension alpha(iu), beta(iu), gamma(iu), rho(iu), y(iu)
!ccccccccccccccccccccccccccccccccccccccccccccccc
! e = 下二对角 (i-2)
! f = 下主对角 (i-1)
! a = 主对角 (i)
! b = 上主对角 (i+1)
! c = 上二对角 (i+2)
! d = 右端向量
!cccccccccccccccccccccccccccccccccccccccccccccccccc
if(n0 .ge. n) return
! 前向消去
do k = n0, n
if(k .eq. n0) then
! 首行
alpha(k) = a(k)
beta(k) = b(k)
gamma(k) = c(k)
rho(k) = d(k)
else if(k .eq. n0+1) then
! 第二行
r = f(k) / alpha(k-1)
alpha(k) = a(k) - r * beta(k-1)
beta(k) = b(k)
gamma(k) = c(k)
rho(k) = d(k) - r * rho(k-1)
else
! 一般行处理
r1 = f(k) / alpha(k-1)
r2 = e(k) / alpha(k-2)
alpha(k) = a(k) - r1 * beta(k-1) - r2 * gamma(k-2)
beta(k) = b(k) - r1 * gamma(k-1)
gamma(k) = c(k)
rho(k) = d(k) - r1 * rho(k-1) - r2 * rho(k-2)
end if
end do
! 回代
y(n) = rho(n) / alpha(n)
y(n-1) = (rho(n-1) - beta(n-1)*y(n)) / alpha(n-1)
do k = n-2, n0, -1
y(k) = (rho(k) - beta(k)*y(k+1) - gamma(k)*y(k+2)) / alpha(k)
end do
! 将解复制回d数组
do k = n0, n
d(k) = y(k)
end do
return
end
将上面这个程序转换为MATLAB程序
最新发布
扫一扫
3966
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
