c Example B c c This is a simplification of Example 1 of B.P. Sommeijer, L.F. Shampine, c and J.G. Verwer, RKC: an Explicit Solver for Parabolic PDEs that shows c the use of RKC on a substantial problem. Semi-discretization of the c heat equation in three space variables results in 19**3 = 6859 equations. c The inhomogeneous term and the boundary conditions have been specified so c that there is an analytical solution evaluated in sol(x,y,z,t). The c maximum error of the numerical solution at TEND is measured by comparison c to a reference solution computed to high accuracy, so the error reported c is the error of the time integration, not the difference between the c solutions of the ODEs and the PDE. Some statistics about the integration c are also displayed. c c WARNING: This program expects the file exb.ref containing the reference c solution to be present in the same directory. c integer ndim parameter (ndim=19*19*19) integer info(4),idid double precision t,tend,rtol,atol double precision y(ndim), work(8+4*ndim) integer neqn,i double precision yref(ndim),error integer nfe,nsteps,naccpt,nrejct,nfesig,maxm common /rkcdid/ nfe,nsteps,naccpt,nrejct,nfesig,maxm integer nx,ny,nz common /grid/ nx,ny,nz external f c t = 0d0 tend = 0.7d0 c----------------------------------------- c Define the mesh and the number of ODEs. c Define the initial values. c----------------------------------------- nx = 19 ny = 19 nz = 19 neqn = nx*ny*nz call exact(neqn,t,y) c-------------------------------------- c Load the reference solution at TEND. c-------------------------------------- open(10,file='exb.ref') read(10,*) yref c--------------------------------------------------------------- c info(1) = 1 -- compute a solution at TEND only. c info(2) = 1 -- SPCRAD returns a bound on the spectral radius. c info(3) = 1 -- the Jacobian is constant. c info(4) = 0 -- ATOL is a scalar. c--------------------------------------------------------------- info(1) = 1 info(2) = 1 info(3) = 1 info(4) = 0 c rtol = 1d-2 atol = rtol c idid = 0 call rkc(neqn,f,y,t,tend,rtol,atol,info,work,idid) c--------------------------------- c Was the integration successful? c--------------------------------- if(idid .ne. 1) then write(*,*) ' Failed at t = ',t,' with idid = ',idid stop endif c error = 0d0 do 10 i = 1,neqn error = max(error,abs(y(i) - yref(i))) 10 continue write(*,'(/a,d8.1,a,f6.1,a,d8.2)') ' With rtol = atol =',rtol, & ', the maximum error at tend =',tend,' was',error write(*,'(a,i5,a)') ' The integration cost',nfe, & ' function evaluations.' write(*,'(a,i4,a,i3,a)') ' There were',nsteps,' steps (', & nrejct,' rejected).' write(*,'(a,i4/)') ' The maximum number of stages used was', & maxm end subroutine exact(neqn,t,y) integer neqn double precision t,y(neqn) integer i,j,k,l double precision dx,dy,dz,sol integer nx,ny,nz common /grid/ nx,ny,nz c dx = 1d0/(nx+1) dy = 1d0/(ny+1) dz = 1d0/(nz+1) do 30 i = 1,nx do 20 j = 1,ny do 10 k = 1,nz l = i + (j-1)*nx + (k-1)*nx*ny y(l) = sol(i*dx,j*dy,k*dz,t) 10 continue 20 continue 30 continue return end double precision function sol(x,y,z,t) double precision x,y,z,t double precision arg arg = 5d0*(x + 2d0*y + 1.5d0*z - 0.5d0 - t) sol = tanh(arg) return end double precision function spcrad(neqn,t,y) integer neqn double precision t,y(neqn) integer nx,ny,nz common /grid/ nx,ny,nz spcrad = 4d0*((nx+1)**2 + (ny+1)**2 + (nz+1)**2) return end subroutine f(neqn,t,y,dydt) integer neqn double precision t,y(neqn),dydt(neqn) integer i,j,k,l double precision u(0:20,0:20,0:20),dx,dy,dz,dxsq,dysq,dzsq, & arg,sh,ch,sol integer nx,ny,nz common /grid/ nx,ny,nz c dx = 1d0/(nx+1) dy = 1d0/(ny+1) dz = 1d0/(nz+1) dxsq = dx*dx dysq = dy*dy dzsq = dz*dz do 30 i = 1,nx do 20 j = 1,ny do 10 k = 1,nz u(i,j,k) = y(i + (j-1)*nx + (k-1)*nx*ny) 10 continue 20 continue 30 continue c do 50 i = 1,nx do 40 j = 1,ny u(i,j,0) = sol(i*dx,j*dy,0d0,t) u(i,j,nz+1) = sol(i*dx,j*dy,1d0,t) 40 continue 50 continue c do 70 i = 1,nx do 60 k = 1,nz u(i,0,k) = sol(i*dx,0d0,k*dz,t) u(i,ny+1,k) = sol(i*dx,1d0,k*dz,t) 60 continue 70 continue c do 90 j = 1,ny do 80 k = 1,nz u(0,j,k) = sol(0d0,j*dy,k*dz,t) u(nx+1,j,k) = sol(1d0,j*dy,k*dz,t) 80 continue 90 continue c do 120 i = 1,nx do 110 j = 1,ny do 100 k = 1,nz arg = 5d0*(i*dx + 2d0*j*dy + 1.5d0*k*dz - 0.5d0 - t) sh = sinh(arg) ch = cosh(arg) l = i + (j-1)*nx + (k-1)*nx*ny dydt(l) = (u(i-1,j,k) - 2d0*u(i,j,k) + u(i+1,j,k))/dxsq + & (u(i,j-1,k) - 2d0*u(i,j,k) + u(i,j+1,k))/dysq + & (u(i,j,k-1) - 2d0*u(i,j,k) + u(i,j,k+1))/dzsq + & (-5d0*ch + 362.5d0*sh)/(ch**3) 100 continue 110 continue 120 continue c return end