!--------------------------------------------------------------------------- ! Program solves the 1-d advection equation (u_t - u_x) = 0 ! ! To compile use "f77 -g -n32 -L/usr/localn32/lib advect1d.f \ ! -lp329f -lsvs -lrnpl -lmfhdf -ldf -ljpeg -lz -lsv -lm -o advect1d !--------------------------------------------------------------------------- program advect1d implicit none character *2 itoc !--------------------------------------------------------------------------- ! on is counter for number of time steps ! k is counter for number of space steps ! nc is the number of crossing times !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! u(maxn) is the array for the exact solution ! uh(maxn,3) is the array for the discrete solution ! e2(maxn) is the array for the second-order error ! eh(maxn) is an estimate of the second-order error equal to uh(k)-u(k) !--------------------------------------------------------------------------- integer maxn, level, olevel, nc, on, temp parameter (maxn = 10000) real*8 x(maxn), u(maxn), y(maxn), ex(maxn) real*8 e2(maxn), eh(maxn), uh(maxn,3) real*8 t, lambda, xprim real*8 r8arg integer iargc, i4arg !--------------------------------------------------------------------------- ! Nt is the number of time steps ! Nx is the number of space steps ! outnum is the number of time steps between outputs !--------------------------------------------------------------------------- integer Nx, Nt, outnum, k integer n, nm1, np1, tmp real*8 deltax, deltat, e parameter (e = 2.7182818285) !--------------------------------------------------------------------------- ! Argument parsing !--------------------------------------------------------------------------- if (iargc() .eq. 4) then level = i4arg(1, -1) olevel = i4arg(2, -1) nc = i4arg(3, -1) lambda = r8arg(4, -1.0d0) else print *,"Usage: advect1d " stop endif !--------------------------------------------------------------------------- ! Error checking for the arguments !--------------------------------------------------------------------------- if (level .lt. 0) then print *,"level must be a positive number" stop endif if (olevel .gt. level) then print *,"olevel must be less than level" stop endif if (nc .lt. 0) then print *,"nc must be a positive number" stop endif if (lambda .gt. 1.0d0 .or. lambda .lt. 0) then print *,"lambda must be between 0 and 1" stop endif Nx = (2**(level))+1 outnum = 2**(level-olevel) Nt = (Nx-1)*nc deltax = 1.0d0/Nx deltat = lambda*deltax nm1 = 1 n = 2 np1 = 3 !--------------------------------------------------------------------------- ! Initialize and output (if necessary) the t=0 and t=dt values for u(x,t) ! ! The leapfrog method works with values of u at three different ! time steps. These values are saved in uh(1,nm1), uh(1,n) and uh(1,np1) !--------------------------------------------------------------------------- t=0.0d0 do k=1,(Nx) x(k) = k*deltax uh(k,nm1) = e**(-(((x(k)-0.5d0)/(0.5d-1))**2)) u(k) = e**(-(((x(k)-0.5d0)/(0.5d-1))**2)) e2(k) = 0 eh(k) = uh(k,nm1)-u(k) end do call vsxynt('uh'//itoc(level),t,x,uh(1,nm1),Nx) !--------------------------------------------------------------------------- ! The exact error is e2 = (deltax)^2*t*(1/6)*(1-lambda^2)*u_xxx !--------------------------------------------------------------------------- t=deltat do k=1,Nx xprim = mod((x(k)+deltat),1.0d0) uh(k,n) = e**(-(((xprim-0.5d0)/(0.5d-1))**2)) u(k) = uh(k,n) temp = ((xprim-0.5d0)/(0.5d-1)) e2(k) = t*(deltax**2)*(1.0d0/6.0d0)*(1-(lambda**2)) e2(k) = e2(k)*(3.2d4*(3.0d0*temp + 2.0d0*(temp**3))) e2(k) = e2(k)*u(k) eh(k) = uh(k,n)-u(k) end do if (outnum .eq. 1) then call vsxynt('uh'//itoc(level),t,x,uh(1,n),Nx) end if open(1,file = 'data', form = 'formatted') do on=2,Nt t = on*deltat !--------------------------------------------------------------------------- ! Calculate uh for the end points, x = 0 and x = 1. Use the periodicity ! of the solution !--------------------------------------------------------------------------- uh(1,np1) = uh(1,nm1) + lambda*(uh(2,n)-uh((Nx-1),n)) uh(Nx,np1) = uh(Nx,nm1) + lambda*(uh(2,n)-uh((Nx-1),n)) !--------------------------------------------------------------------------- ! Now calculate the rest of the points !--------------------------------------------------------------------------- do k=2,Nx-1 uh(k,np1) = uh(k,nm1) + lambda*(uh(k+1,n)-uh(k-1,n)) end do !--------------------------------------------------------------------------- ! Calculate the exact solution, the exact error and the estimated error !--------------------------------------------------------------------------- do k=1,Nx xprim = mod((x(k)+t),1.0d0) u(k) = e**(-(((xprim-0.5d0)/(0.5d-1))**2)) temp = ((xprim-0.5d0)/(0.5d-1)) e2(k) = t*(deltax**2)*(1.0d0/6.0d0)*(1-(lambda**2)) e2(k) = e2(k)*(3.2d4*(3.0d0*temp + 2.0d0*(temp**3))) e2(k) = e2(k)*u(k) eh(k) = uh(k,np1)-u(k) end do !--------------------------------------------------------------------------- ! Write the results to .sdf files every outnum time steps !--------------------------------------------------------------------------- if (mod(on, outnum) .eq. 0) then call vsxynt('uh'//itoc(level),t,x,uh(1,np1),Nx) call vsxynt('u'//itoc(level),t,x,u,Nx) call vsxynt('e2'//itoc(level),t,x,e2,Nx) call vsxynt('eh'//itoc(level),t,x,eh,Nx) write (1,*) t,Nx do k=1,Nx write (1,*)x(k),uh(k,nm1),uh(k,n) end do endif tmp = nm1 nm1 = n n = np1 np1 = tmp end do end