program tmol1 implicit none save c..test the method of lines solution of the 1d diffusion equation c..test problems from c..g.d. smith "numerical solution of partial differential c.. equations, finite difference methods" c.. oxford univ. press 1978 c..w.e. schiesser "the numerical method of lines" c.. academic press 1991 c..communicate grid and boundary conditions integer ngridmax parameter (ngridmax = 200) integer jl,jr,ngrid,icall,ijac double precision xgrid(ngridmax),xnl1,xlpow,xnl2,xnr1,xrpow,xnr2 common /tmolc1/ xgrid,xnl1,xlpow,xnl2,xnr1,xrpow,xnr2, 1 jl,jr,ngrid,icall,ijac c..for the ode integration external derv2,jac2,stifbs3 integer ndim,nok,nbad,kount,iprint parameter (ndim=100) double precision start,stptry,stpmin,stopp,bc(ngridmax),tol,dxsav, 1 ttime(ndim),ans(ngridmax,ndim),odescal c..local variables external aroot integer i,j,n,njj,icase,nb,niter double precision xlength,dx,te(ngridmax),xtime, 1 errmax,percent,maxper,xwhere,sum, 2 lo,hi,aroot,xb1(1000),xb2(1000),alfan,alfan2, 3 zbrent,rtol parameter (rtol = 1.0d-12) c..math constants double precision pi parameter (pi = 3.1415926535897932384d0) c..format statements 01 format(1x,a,i6,t15,a,i6,t30,a,i6) 20 format(1x,a,0p12f8.4,/,(8(4x,0p12f8.4))) 21 format(1x,a,0p12f8.4) 22 format(1x,'time=',1pe11.3) 23 format(1x,a,1pe12.4,a,1pe12.4) c..for each test case do icase = 1,6 icall = 0 ijac = 0 c..schiesser page 16, simple dirchlett if (icase .eq. 1) then jl = 0 jr = 0 xnl1 = 0.0d0 xlpow = 1.0d0 xnl2 = 0.0d0 xnr1 = 0.0d0 xrpow = 1.0d0 xnr2 = 0.0d0 xlength = 1.0d0 njj = 1 ngrid = 10*njj + 1 dx = xlength/float(ngrid-1) do i=1,ngrid xgrid(i) = float(i-1) * dx bc(i) = sin(pi*xgrid(i)/xlength) enddo write(6,*) write(6,*) 'example from schiesser page 16 ngrid=',ngrid c..schiesser page 28, simple neumann else if (icase .eq. 2) then jl = 1 jr = 1 xnl1 = 0.0d0 xlpow = 1.0d0 xnl2 = 0.0d0 xnr1 = 0.0d0 xrpow = 1.0d0 xnr2 = 0.0d0 xlength = 1.0d0 njj = 1 ngrid = 10*njj + 1 dx = xlength/float(ngrid-1) do i=1,ngrid xgrid(i) = float(i-1) * dx bc(i) = cos(pi*xgrid(i)/xlength) enddo write(6,*) write(6,*) 'example from schiesser page 28 ngrid=',ngrid c..smith page 21, initial discontinuity at x=0.5 else if (icase .eq. 3) then jl = 0 jr = 0 xnl1 = 0.0d0 xlpow = 1.0d0 xnl2 = 0.0d0 xnr1 = 0.0d0 xrpow = 1.0d0 xnr2 = 0.0d0 xlength = 1.0d0 njj = 1 ngrid = 10*njj + 1 dx = xlength/float(ngrid-1) do i=1,ngrid xgrid(i) = float(i-1) * dx if (xgrid(i) .le. 0.5) then bc(i) = 2.0d0 * xgrid(i) else bc(i) = 2.0d0 * (1.0d0 - xgrid(i)) end if enddo write(6,*) write(6,*) 'example from smith page 21 ngrid=',ngrid c..smith page 35, exercises xnl1 and xnr1 else if (icase .eq. 4) then jl = 1 jr = 1 xnl1 = 1.0d0 xlpow = 1.0d0 xnl2 = 0.0d0 xnr1 = -1.0d0 xrpow = 1.0d0 xnr2 = 0.0d0 xlength = 1.0d0 njj = 1 ngrid = 10*njj + 1 dx = xlength/float(ngrid-1) do i=1,ngrid xgrid(i) = float(i-1) * dx bc(i) = 1.0d0 enddo write(6,*) write(6,*) 'example from smith page 35 ngrid=',ngrid c..smith page 62, different neumann conditions on each end else if (icase .eq. 5) then jl = 1 jr = 1 xnl1 = 0.0d0 xlpow = 1.0d0 xnl2 = 0.0d0 xnr1 = 0.0d0 xrpow = 1.0d0 xnr2 = 1.0d0 xlength = 0.5d0 njj = 1 ngrid = 5*njj + 1 dx = xlength/float(ngrid-1) do i=1,ngrid xgrid(i) = float(i-1) * dx bc(i) = 0.0d0 enddo write(6,*) write(6,*) 'example from smith page 62 ngrid=',ngrid c..scheisser page 54, nonlinear boundary else if (icase .eq. 6) then jl = 0 jr = 1 xnl1 = 0.0d0 xlpow = 1.0d0 xnl2 = 0.0d0 xnr1 = -1.0d0 xrpow = 4.0d0 xnr2 = 1.0d0 xlength = 1.0d0 njj = 10 ngrid = 5*njj + 1 dx = xlength/float(ngrid-1) do i=1,ngrid xgrid(i) = float(i-1) * dx bc(i) = 0.0d0 enddo bc(1) = 1.0d0 write(6,*) write(6,*) 'example from schiesser page 54 ngrid=',ngrid c..end of test case ifs endif c..initialize some ode control parameters start = 0.0d0 stopp = 0.0d0 stptry = 1.0d-6 stpmin = 1.0d-10 odescal = 1.0d0 iprint = 0 dxsav = stopp tol = 1.0d-6 c..get the time to evolve 100 continue write(6,*) write(6,*) 'give time to evolve (0.0 for next case; -1 to end)=>' read(5,*) xtime if (xtime .eq. 0.0) goto 200 if (xtime .eq. -1.0) stop 'normal termination' c..set the start,step stop times start = stopp stopp = start + xtime if (start .ne. 0.0) stptry = stopp dxsav = stopp c..integrate call odeint(start,stptry,stpmin,stopp,bc, 1 tol,dxsav,ndim, 2 ttime,ans,ndim,ngridmax,ndim,ngrid, 3 nok,nbad,kount,odescal,iprint, 4 derv2,jac2,stifbs3) c..number of good and bad (but redone steps) write(6,01) 'nok=',nok,' nbad=',nbad,' kount=',kount write(6,01) 'icall=',icall,' ijac=',ijac c..loop through the output points do j=1,kount c..compute the exact answer c..scheisser page 16, simple dirchlett if (icase .eq. 1) then do i=1,ngrid,njj te(i) = exp((-pi**2/xlength**2)*ttime(j)) 1 * sin(pi*xgrid(i)/xlength) enddo c..scheisser page 28, simple neumann else if (icase .eq. 2) then do i=1,ngrid,njj te(i) = exp((-pi**2/xlength**2)*ttime(j)) 1 * cos(pi*xgrid(i)/xlength) enddo c..smith page 21, initial discontinuity at x=0.5 else if (icase .eq. 3) then do i=1,ngrid,njj sum = 0.0d0 do n=1,100 sum = sum + sin(0.5d0*pi*n)*sin(n*pi*xgrid(i))/n**2 1 * exp(-(n*pi)**2*ttime(j)) enddo te(i) = 8.0d0/pi**2 * sum enddo c..smith page 35, exercises xnl1 and xnr1 c..get the positive roots of alfa*tan(alfa)=0.5 c..ignore the sign changes at n*pi/2 else if (icase .eq. 4) then lo = 0.0 hi = 1000.0d0 nb = 5000 call zbrak(aroot,lo,hi,5000,xb1,xb2,nb) do i=1,ngrid,njj sum = 0.0d0 do n=1,nb/2,2 alfan = zbrent(aroot,xb1(n),xb2(n),rtol,niter) alfan2 = alfan*alfan sum = sum + 1 exp(-4.0d0*alfan2*ttime(j)) 2 / (cos(alfan)*(3.0d0+4.0d0*alfan2)) 3 * cos(2.0d0*alfan*(xgrid(i)-0.5d0)) enddo te(i) = 4.0d0 * sum enddo c..smith page 62, different neumann conditions on each end else if (icase .eq. 5) then do i=1,ngrid,njj sum = 0.0d0 do n=1,100 sum = sum + 1 (-1.0d0)**n * exp(-4.0d0*(pi*n)**2*ttime(j)) 2 * cos(2.0d0*n*pi*xgrid(i))/n**2 enddo te(i) = 2.0d0*ttime(j) + 0.5d0* ( 1 (12.0d0*xgrid(i)**2 - 1.0d0)/6.0d0 2 - 2.0d0/pi**2 * sum) enddo c..scheisser page 54, nonlinear boundary else if (icase .eq. 6) then do i=1,ngrid,njj te(i) = 0.0d0 enddo c..end of test case ifs end if c..get the max percentage difference errmax = -1.0e30 do i=1,ngrid,njj if (abs(te(i)) .gt. 1.0d-6 .and. ans(i,j) .gt. 1.0d-6) then percent = ( (te(i) - ans(i,j))/te(i) ) * 100.0d0 errmax = max(errmax,abs(percent)) if (errmax .eq. abs(percent)) then xwhere = xgrid(i) maxper = percent end if end if enddo c..write out the grid, numerical solution and exact solution write(6,*) write(6,22) ttime(j) if (ngrid .le. 12) then write(6,21) 'x =',(xgrid(i),i=1,ngrid,njj) write(6,21) 't =',(ans(i,j),i=1,ngrid,njj) write(6,21) 'te=',(te(i),i=1,ngrid,njj) else write(6,20) 'x =',(xgrid(i),i=1,ngrid,njj) write(6,20) 't =',(ans(i,j), i=1,ngrid,njj) write(6,20) 'te=',(te(i),i=1,ngrid,njj) end if write(6,23) 'max % error=',maxper,' occurs at x=',xwhere c..end of kount loop enddo goto 100 c..back for another test case 200 continue enddo stop 'normal termination' end c..the next two routines set up the ode's and their jacobian subroutine derv2(t,y,dydt) implicit none save c..this routine implements the method of lines solution to the c..1d diffusion equation c.. c.. de(T)/dt = 1/rho * d/dx (sigma(T) * dT/dx) + s(T) c.. c..where c..T = temperature c..e(T) = specific internal energy c..rho(T) = density c..sigma(T) = thermal diffusion coefficient c..s(T) = energy generation rate c.. c..alternatively, e(t) could be the enthalpy or just the temperature, c..in which case the meaning of rho(T) may change. c.. c..dirchlett or neumann boundary conditions at either end may be enforced, c..and a non-uniform grid may be used. the descritization of the c..right hand side is second order accurate in space. c.. c..input: c..t = time, not used since the pde or odes don't explicitly depend on time c..y = vector of the function e(T) c.. c..output: c..dydt = vector of the odes c..declare the pass double precision t,y(1),dydt(1) c..communicate grid and boundary conditions integer ngridmax parameter (ngridmax = 200) integer jl,jr,ngrid,icall,ijac double precision xgrid(ngridmax),xnl1,xlpow,xnl2,xnr1,xrpow,xnr2 common /tmolc1/ xgrid,xnl1,xlpow,xnl2,xnr1,xrpow,xnr2, 1 jl,jr,ngrid,icall,ijac c..local variables integer i double precision yx(ngridmax),yxx(ngridmax),sx(ngridmax), 1 ca(3,ngridmax),cb(4,ngridmax), 2 zero double precision sigma(ngridmax),dsigmadt(ngridmax), 1 fac(ngridmax),dfacdt(ngridmax), 2 sdot(ngrid),dsdotdt(ngrid) parameter (zero = 0.0d0) c..count the number of calls to this routine icall = icall + 1 c..get the coefficients call pdecof(ngrid,y,sigma,dsigmadt,fac,dfacdt,sdot,dsdotdt) c..2nd order 2nd derivative applying any boundary conditions, nonuniform grid yx(1) = xnl1 yx(2) = xlpow yx(3) = xnl2 yx(ngrid-2) = xnr1 yx(ngrid-1) = xrpow yx(ngrid) = xnr2 call fss004(xgrid,ngrid,y,yx,yxx,cb,jl,jr) c..2nd order first derivatives, apply any neumann boundary conditions call fss002(xgrid,ngrid,y,yx,ca) if (jl .eq. 1) yx(1) = xnl1*y(1)**xlpow + xnl2 if (jr .eq. 1) yx(ngrid) = xnr1*y(ngrid)**xrpow + xnr2 call fss002(xgrid,ngrid,sigma,sx,ca) c..the odes, apply any dirchlett boundary conditions do i=1,ngrid dydt(i) = fac(i)*(sigma(i) * yxx(i) + sx(i)*yx(i)) + sdot(i) enddo if (jl .eq. 0) dydt(1) = zero if (jr .eq. 0) dydt(ngrid) = zero return end subroutine jac2(t,y,dfdy,nlog,nphys) implicit none save c..this routine sets up the pentadiagonal jacobian of the c..odes contained in the routine above. c..formally the matrix is tridiagonal with two extra entries, alfa and beta: c.. b1*u1 + c1*u2 + alfa*u3 c.. a2*u1 + b2*u2 + c2*u3 c.. a3*u2 + b3*u4 + c3*u5 c.. ... c.. beta*u(n-2) + an*u(n-1) + bn*u(n) c.. c..one could use a tridiagonal supplemented with the sherman-woodbury formula c..to solve this system, but it would take at least two, and probably four, c..trigiagonal solutions. i opt to simply solve one pentadiagonal system c..declare the pass integer nlog,nphys double precision t,y(1),dfdy(5,nphys) c..communicate grid and boundary conditions integer ngridmax parameter (ngridmax = 200) integer jl,jr,ngrid,icall,ijac double precision xgrid(ngridmax),xnl1,xlpow,xnl2,xnr1,xrpow,xnr2 common /tmolc1/ xgrid,xnl1,xlpow,xnl2,xnr1,xrpow,xnr2, 1 jl,jr,ngrid,icall,ijac c..local variables integer i,j double precision yx(ngridmax),yxsav(ngridmax),yxx(ngridmax), 1 sx(ngridmax),ca(3,ngridmax),cb(4,ngridmax), 2 sigma(ngridmax),dsigmadt(ngridmax), 3 fac(ngridmax),dfacdt(ngridmax), 4 sdot(ngrid),dsdotdt(ngrid) c..count the number of calls to this routine ijac = ijac + 1 c..zero the jacobian do j=1,nphys do i=1,5 dfdy(i,j) = 0.0d0 enddo enddo c..get the coefficients call pdecof(ngrid,y,sigma,dsigmadt,fac,dfacdt,sdot,dsdotdt) c..2nd order 2nd derivative, nonuniform grid yx(1) = xnl1 yx(2) = xlpow yx(3) = xnl2 yx(ngrid-2) = xnr1 yx(ngrid-1) = xrpow yx(ngrid) = xnr2 call fss004(xgrid,ngrid,y,yx,yxx,cb,jl,jr) c..2nd order first derivatives, apply any neumann boundary conditions call fss002(xgrid,ngrid,y,yx,ca) if (jl .eq. 1) yx(1) = xnl1*y(1)**xlpow + xnl2 if (jr .eq. 1) yx(ngrid) = xnr1*y(ngrid)**xrpow + xnr2 call fss002(xgrid,ngrid,sigma,sx,ca) c..easiest to rescale the coefficents here do i=1,ngrid do j=1,4 cb(j,i) = cb(j,i) * sigma(i) enddo enddo c..apply neumann boundary condition c..the entry dfdy(5,i) is the alfa commented on above if (jl .eq. 1) then i = 1 dfdy(3,i) = fac(i) * (cb(1,i) + dsigmadt(i)*yxx(i) 1 + ca(1,i)*dsigmadt(i)*yx(1) 2 + sx(i)*xnl1*xlpow*y(1)**(xlpow-1.0d0) ) 3 + dfacdt(i)*(sigma(i) * yxx(i) + sx(i)*yx(i)) 4 + dsdotdt(i) dfdy(4,i) = fac(i) * (cb(2,i) + ca(2,i)*dsigmadt(i+1)*yx(1) ) dfdy(5,i) = fac(i) * (cb(3,i) + ca(3,i)*dsigmadt(i+2)*yx(1) ) endif c..middle of the grid do i=2,ngrid-1 dfdy(2,i) = fac(i) * ( 1 cb(1,i) + ca(1,i)* (dsigmadt(i-1)*yx(i) + sx(i)) ) dfdy(3,i) = fac(i) * ( 1 cb(2,i) + dsigmadt(i)*yxx(i) 2 + ca(2,i)*(dsigmadt(i)*yx(i) + sx(i)) ) 3 + dfacdt(i)*(sigma(i) * yxx(i) + sx(i)*yx(i)) 4 + dsdotdt(i) dfdy(4,i) = fac(i) * ( 1 cb(3,i) + ca(3,i)*(dsigmadt(i+1)*yx(i) + sx(i)) ) enddo c..apply neumann boundary condition c..the entry dfdy(1,ngrid) is the beta commented on above if (jr .eq. 1) then i = ngrid dfdy(1,i) = fac(i) * (cb(2,i) + ca(1,i)*dsigmadt(i-2)*yx(i) ) dfdy(2,i) = fac(i) * (cb(3,i) + ca(2,i)*dsigmadt(i-1)*yx(i) ) dfdy(3,i) = fac(i) * (cb(4,i) + dsigmadt(i)*yxx(i) 1 + ca(3,i)*dsigmadt(i)*yx(i) 2 + sx(i)*xnr1*xrpow*y(i)**(xrpow-1.0d0) ) 3 + dfacdt(i)*(sigma(i) * yxx(i) + sx(i)*yx(i)) 4 + dsdotdt(i) endif return end c..the next routine returns the coefficients, and their derivatives, c..of the diffusion equation subroutine pdecof(ngrid,temp,sig,dsigdt,fac,dfacdt,sdot,dsdotdt) implicit none save c..loads the coefficents of the diffusion equation c..input: c..ngrid = number of grid points c..temp = temperature array on the grid points c..output: c..sig = thermal diffusion coefficient c..dsigdt = d(sigma)/d(temperature) c..fac = multiplier in front of the diffusion operator c..dfacdt = d(rho)/d(temperature) c..sdot = energy generation rate c..dsdotdt = d(sdot)/d(temperature) c..declare the pass integer ngrid double precision temp(ngrid), 1 sig(ngrid),dsigdt(ngrid), 2 fac(ngrid),dfacdt(ngrid), 3 sdot(ngrid),dsdotdt(ngrid) c..local variables integer i c..load the constant coefficients for the test problems do i=1,ngrid sig(i) = 1.0d0 dsigdt(i) = 0.0d0 fac(i) = 1.0d0 dfacdt(i) = 0.0d0 sdot(i) = 0.0d0 dsdotdt(i) = 0.0d0 enddo return end c..the next three routines are only needed for the analytic c..solution of one of the test problems double precision function aroot(a) implicit none save c..function for the problem of smith page 35 double precision a aroot = a*tan(a) - 0.5d0 return end double precision function zbrent(func,x1,x2,tol,niter) implicit none save c..using brent's method this routine finds the root of a function func c..between the limits x1 and x2. the root is when accuracy is less than tol. c.. c..note: eps the the machine floating point precision c.. c..declare external func integer niter,itmax,iter parameter (itmax = 100) double precision func,x1,x2,tol,a,b,c,d,e,fa, 1 fb,fc,xm,tol1,p,q,r,s,eps parameter (eps=3.0d-15) c..initialize niter = 0 a = x1 b = x2 fa = func(a) fb = func(b) if ( (fa .gt. 0.0 .and. fb .gt. 0.0) .or. 1 (fa .lt. 0.0 .and. fb .lt. 0.0) ) then write(6,100) x1,fa,x2,fb 100 format(1x,' x1=',1pe11.3,' f(x1)=',1pe11.3,/, 1 1x,' x2=',1pe11.3,' f(x2)=',1pe11.3) stop 'root not bracketed in routine zbrent' end if c = b fc = fb c..rename a,b,c and adjusting bound interval d do iter =1,itmax niter = niter + 1 if ( (fb .gt. 0.0 .and. fc .gt. 0.0) .or. 1 (fb .lt. 0.0 .and. fc .lt. 0.0) ) then c = a fc = fa d = b-a e = d end if if (abs(fc) .lt. abs(fb)) then a = b b = c c = a fa = fb fb = fc fc = fa end if tol1 = 2.0d0 * eps * abs(b) + 0.5d0 * tol xm = 0.5d0 * (c-b) c..convergence check if (abs(xm) .le. tol1 .or. fb .eq. 0.0) then zbrent = b return end if c..attempt quadratic interpolation if (abs(e) .ge. tol1 .and. abs(fa) .gt. abs(fb)) then s = fb/fa if (a .eq. c) then p = 2.0d0 * xm * s q = 1.0d0 - s else q = fa/fc r = fb/fc p = s * (2.0d0 * xm * q *(q-r) - (b-a)*(r - 1.0d0)) q = (q - 1.0d0) * (r - 1.0d0) * (s - 1.0d0) end if c..check if in bounds if (p .gt. 0.0) q = -q p = abs(p) c..accept interpolation if (2.0d0*p .lt. min(3.0d0*xm*q - abs(tol1*q),abs(e*q))) then e = d d = p/q c..or bisect else d = xm e = d end if c..bounds decreasing to slowly use bisection else d = xm e = d end if c..move best guess to a a = b fa = fb if (abs(d) .gt. tol1) then b = b + d else b = b + sign(tol1,xm) end if fb = func(b) enddo stop 'too many iterations in routine zbrent' end subroutine zbrak(fx,x1,x2,n,xb1,xb2,nb) implicit none save c..given a function fx defined on the interval x1 to x2, subdivide the c..interval into n equally spaced segments, and search for zero crossings c..of the function. nb is input as the maximum number of roots sought and c..is reset to the number of bracketing pairs xb1,xb2 that are found. c.. c..declare external fx integer n,nb,nbb,i double precision fx,x1,x2,xb1(nb),xb2(nb),x,dx,fp,fc c..initialize and determine the appropriate mesh spacing nbb = 0 x = x1 dx = (x2-x1)/n fp = fx(x) c..loop over all intervals do i=1,n x = x + dx fc = fx(x) c..if a sign change occurs then record the values for the bounds if (fc*fp .le. 0.0) then nbb = nbb + 1 xb1(nbb) = x - dx xb2(nbb) = x if (nbb .eq. nb) goto 10 end if fp = fc enddo 10 nb = nbb return end c..the next 3 routines implement finite difference formulas c..on non-uniform meshes subroutine fdcoef(mord,nord,x0,grid,coef) implicit none save c..this routine implements simple recursions for calculating the weights c..of finite difference formulas for any order of derivative and any order c..of accuracy on one-dimensional grids with arbitrary spacing. c..from begnt fornberg, math. computation 51, 699, 1988 c.."generation of finite difference formulas on arbitrarily spaced grids" c..input: c..m = order of derivative m=2=first derivative m=3=second derivative ... c..n = order of accuracy n=2 first order, n=3=second order ... c..x0 = point at which to evaluate the derivative c..grid = spatial grid c..output: c..coef = coefficients of the finite difference formula c..declare integer mord,nord,nu,nn,mm,nmmin,mmax,nmax parameter (mmax=8, nmax=10) double precision x0,grid(nord),coef(nord), 1 weight(mmax,nmax,nmax),c1,c2,c3,c4,pfac c..zero the weight array do nu=1,nord do nn=1,nord do mm=1,mord weight(mm,nn,nu) = 0.0d0 enddo enddo enddo weight(1,1,1) = 1.0d0 c1 = 1.0d0 nmmin = min(nord,mord) do nn = 2,nord c2 = 1.0d0 do nu=1,nn-1 c3 = grid(nn) - grid(nu) c2 = c2 * c3 c4 = 1.0d0/c3 pfac = grid(nn) - x0 weight(1,nn,nu) = c4 * ( pfac * weight(1,nn-1,nu) ) do mm=2,nmmin weight(mm,nn,nu) = c4 * ( pfac * weight(mm,nn-1,nu) 1 - dfloat(mm-1)*weight(mm-1,nn-1,nu) ) enddo enddo pfac = (grid(nn-1) - x0) weight(1,nn,nn) = c1/c2 * (-pfac*weight(1,nn-1,nn-1)) c4 = c1/c2 do mm=2,nmmin weight(mm,nn,nn) = c4 * (dfloat(mm-1)*weight(mm-1,nn-1,nn-1) - 1 pfac*weight(mm,nn-1,nn-1)) enddo c1 = c2 enddo c..load the coefficients do nu = 1,nord coef(nu) = weight(mord,nord,nu) enddo return end subroutine fss002(grid,ngrid,u,du,cofs) implicit none save c..computes a second order accurate first derivative on an arbitrary grid. c..input: c..grid = the grid c..ngrid = number of points in the grid c..u = variable on the grid c..output c..du = first derivative of the variable on the grid c..cofs = coefficients of the derivative, essentially the jacobian c..declare integer m,n parameter (m=2, n=3) integer ngrid,i double precision grid(ngrid),u(ngrid),du(ngrid),cofs(3,ngrid), 1 coef(n) c..first point, one sided call fdcoef(m,n,grid(1),grid(1),coef) i = 1 du(i) = coef(1)*u(i) + coef(2)*u(i+1) + coef(3)*u(i+2) cofs(1,i) = coef(1) cofs(2,i) = coef(2) cofs(3,i) = coef(3) c..middle of the grid, central differences do i=2,ngrid-1 call fdcoef(m,n,grid(i),grid(i-1),coef) du(i) = coef(1)*u(i-1) + coef(2)*u(i) + coef(3)*u(i+1) cofs(1,i) = coef(1) cofs(2,i) = coef(2) cofs(3,i) = coef(3) enddo c..last point, one sided call fdcoef(m,n,grid(ngrid),grid(ngrid-2),coef) i = ngrid du(i) = coef(1)*u(i-2) + coef(2)*u(i-1) + coef(3)*u(i) cofs(1,i) = coef(1) cofs(2,i) = coef(2) cofs(3,i) = coef(3) return end subroutine fss004(grid,ngrid,u,du,duu,cofs,nl,nu) implicit none save c..computes a second order accurate second derivative on an arbitrary grid. c..input: c..grid = the grid c..ngrid = number of points in the grid c..u = variable on the grid c..du = derivative on the grid, only used for neumann boundary conditions c.. supported forms are xnl1 * u(1)**xlow + xnl2 c.. supported forms are xnr1 * u(1)**xrow + xnr2 c.. du(1) contains coefficient xnl1 c.. du(2) contains power xlpow c.. du(3) contains coefficient xnl2 c.. du(ngrid-2) contains coefficieicient xnr1 c.. du(ngrid-1) contains power xnr1 c.. du(ngrid) contains coefficient xnr2 c..nl = type of boundary condition at the left nl=0=dirchlett nl=1=neumann c..nu = type of boundary condition at the right nu=0=dirchlett nu=1=neumann c..output c..duu = second derivative of the variable on the grid c..cofs = coefficients of the derivative, essentially the jacobian c..declare integer m,n parameter (m=3, n=3) integer ngrid,nl,nu,i double precision grid(ngrid),u(ngrid),du(ngrid),duu(ngrid), 1 cofs(4,ngrid),coef(n+1),alf(3),dx,u0 c..first point c..dirchlett, one sided if (nl .eq. 0) then call fdcoef(m,n+1,grid(1),grid(1),coef) duu(1) = coef(1)*u(1) + coef(2)*u(2) + 1 coef(3)*u(3) + coef(4)*u(4) cofs(1,1) = coef(1) cofs(2,1) = coef(2) cofs(3,1) = coef(3) cofs(4,1) = coef(4) c..neumann else if (nl .eq. 1) then dx = grid(2) - grid(1) alf(1) = grid(1) - dx alf(2) = grid(1) alf(3) = grid(2) call fdcoef(m,n,grid(1),alf,coef) u0 = u(2) - 2.0d0*dx*(du(1)*u(1)**du(2) + du(3)) duu(1) = coef(1)*u0 + coef(2)*u(1) + coef(3)*u(2) cofs(1,1) = coef(2) - 2.0d0*dx*du(1)*du(2)*u(1)**(du(2)-1.0d0) cofs(2,1) = coef(3) + coef(1) cofs(3,1) = 0.0d0 cofs(4,1) = 0.0d0 end if c..middle of the grid, central differences do i=2,ngrid-1 call fdcoef(m,n,grid(i),grid(i-1),coef) duu(i) = coef(1)*u(i-1) + coef(2)*u(i) + coef(3)*u(i+1) cofs(1,i) = coef(1) cofs(2,i) = coef(2) cofs(3,i) = coef(3) cofs(4,i) = 0.0d0 enddo c..last point c..dirchlett, one sided if (nu .eq. 0) then call fdcoef(m,n+1,grid(ngrid),grid(ngrid-3),coef) duu(ngrid) = coef(1)*u(ngrid-3) + coef(2)*u(ngrid-2) + 1 coef(3)*u(ngrid-1) + coef(4)*u(ngrid) cofs(1,ngrid) = coef(1) cofs(2,ngrid) = coef(2) cofs(3,ngrid) = coef(3) cofs(4,ngrid) = coef(4) c..neumann else if (nu .eq. 1) then i = ngrid dx = grid(i) - grid(i-1) alf(1) = grid(i-1) alf(2) = grid(i) alf(3) = grid(i) + dx call fdcoef(m,n,grid(i),alf,coef) u0 = u(i-1) + 2.0d0*dx*(du(i-2)*u(i)**du(i-1) + du(i)) duu(i) = coef(1)*u(i-1) + coef(2)*u(i) + coef(3)*u0 cofs(1,i) = 0.0d0 cofs(2,i) = 0.0d0 cofs(3,i) = coef(1) + coef(3) cofs(4,i) = coef(2) + 2.0d0*dx*du(i-2) 1 * du(i-1)*u(i)**(du(i-1)-1.0d0) endif return end c..the next routine drives ode integrations subroutine odeint(start,stptry,stpmin,stopp,bc, 1 eps,dxsav,kmax, 2 xrk,yrk,xphys,yphys,xlogi,ylogi, 3 nok,nbad,kount,odescal,iprint, 4 derivs,jakob,steper) implicit none save c..input are the initial value to start the integration from start, the first c..step size stptry, minimum step size stpmin, the ending value of the c..integration stopp, the initial conditions bc, the desired fractional c..accuracy eps, the incremental value dxsav at which to store the solution, c..the maximum number of values to store kmax, the array xrk (physical c..dimension xphys and logical dimension xlogi) where the solution was c..computed between start and stopp, the array yrk (physical dimension c..yphys,xphys and logical dimension ylogi,xlogi) where the solution to the c..ylogi ode's are stored, the error scaling factor odescal, a flag iprint c..to determine if we print the solution to the screen as it evolves, the c..name of the subroutine derivs that sets up the ode's and c..finally the name of the subroutine steper that will do the integrations c..(valid choices are rkqc,bsstep and stiff). c.. c..output are the number of succesful steps taken nok, the number of bad c..(but retried and then succesful) steps taken, the total number of steps c..taken kount and of course the solution arrays xrk and yrk described above. c.. c..if dxsav is zero, the routine stores the solution at every intermediate c..point; if it is not zero, the solution is done at every value of dxsav. c.. c..nmax is the maximum number of differential equations c.. c..declare external derivs,jakob,steper integer nok,nbad,nmax,stpmax,kmax,kount,xphys, 1 yphys,xlogi,ylogi,iprint,i,j,nstp parameter (nmax = 2000, stpmax=10000) double precision bc(yphys),stptry,stpmin,eps,stopp,start, 1 yscal(nmax),y(nmax),dydx(nmax), 2 dxsav,xrk(xphys),yrk(yphys,xphys),odescal, 3 x,xsav,h,hdid,hnext,zero,one,tiny,ttiny parameter (zero=0.0d0, tiny=1.0d-30, ttiny=1.0d-15) c..popular format statements 100 format(1x,i4,1pe10.2) 101 format(1x,1p12e10.2) c..initialize if (ylogi .gt. yphys) stop 'ylogi > yphys in routine odeint' if (yphys .gt. nmax) stop 'yphys > nmax in routine odeint' x = start h = sign(stptry,stopp-start) nok = 0 nbad = 0 kount = 0 c..store the first step do i=1,ylogi y(i) = bc(i) enddo xsav = x - 2.0d0 * dxsav c..take at most stpmax steps do nstp=1,stpmax call derivs(x,y,dydx) c..scaling vector used to monitor accuracy do i=1,ylogi c..constant fractional accuracy c yscal(i) = abs(y(i)) + tiny c..const. frac. cept near zero c.. yscal(i)=abs(y(i)) + abs(h*dydx(i)) + tiny c..step size dependent accuracy c.. yscal(i) = abs(odescal * h * dydx(i)) + tiny c..for stiffs yscal(i) = max(odescal,y(i)) c..strait scaling (decrease to get more time steps) c.. yscal(i) = odescal enddo c..store the intermediate results if (kmax .gt. 0) then if ( (abs(dxsav) - abs(x-xsav)) .le. ttiny) then if ( kount .lt. (kmax-1) ) then kount = kount+1 xrk(kount) = x do i=1,ylogi yrk(i,kount) = y(i) enddo if (iprint .eq. 1) then write(6,100) kount,xrk(kount) write(6,101) (yrk(j,kount), j=1,ylogi) end if xsav=x end if end if end if c..if the step can overshoot the stop point or the dxsav increment then cut it if ((x+h-stopp)*(x+h-start).gt.zero) h = stopp - x if (dxsav.ne.zero .and. h.gt.(xsav-x+dxsav)) h = xsav + dxsav-x c..do an integration step c call steper(y,dydx,ylogi,x,h,eps,yscal,hdid,hnext,derivs) call steper(y,dydx,ylogi,x,h,eps,yscal,hdid,hnext,derivs,jakob) if (hdid.eq.h) then nok = nok+1 else nbad = nbad+1 end if c..this is the normal exit point, save the final step if (nstp .eq. stpmax .or. (x-stopp)*(stopp-start) .ge. zero) then do i=1,ylogi bc(i) = y(i) enddo if (kmax.ne.0) then kount = kount+1 xrk(kount) = x do i=1,ylogi yrk(i,kount) = y(i) enddo if (iprint .eq. 1) then write(6,100) kount,xrk(kount) write(6,101) (yrk(j,kount), j=1,ylogi) end if end if return end if c..set the step size for the next iteration; stay above stpmin h = hnext if (abs(hnext).lt.stpmin) stop 'hnext < stpmin in odeint' c..back for another iteration or death enddo write(6,*) '> than stpmax steps required in odeint' return end c..the next three routines implment the stiff bader-deuflhard stepper c..for pentadiagonal jacobians subroutine stifbs3(y,dydx,nv,x,htry,eps,yscal,hdid,hnext, 1 derivs,jakob) implicit none save c..for pentadiagonal analytic jacobians c.. c..semi-implicit method for integrating stiff ode's with monitoring c..of local truncation error to adjust stepsize. inputs are the dependent c..variable vector y(1:nv) and its derivative dydx(1:nv) at the start of the c..independent variable x. also input are the stepsize to be attempted htry, c..the required accuracy eps, and the vector yscal against which the error is c..scaled. on output, y and x are replaced by their new values, hdid is the c..stepsize actually accomplished, and hnext is the estimated next stepsize. c..derivs is a user supplied function that computes the right hand side of c..the equations. plug into odeint. c.. c..declare external derivs,jakob logical first,reduct integer nv,nmax,kmaxx,imax parameter (nmax = 500, kmaxx=7, imax=kmaxx+1) integer i,iq,k,kk,km,kmax,kopt,nvold,nseq(imax) double precision y(nv),dydx(nv),x,htry,eps,yscal(nv),hdid,hnext, 1 eps1,epsold,errmax,fact,h,red,scale,work,wrkmin, 2 xest,xnew,a(imax),alf(kmaxx,kmaxx),err(kmaxx), 3 yerr(nmax),ysav(nmax),yseq(nmax),safe1,safe2, 4 redmax,redmin,tiny,scalmx, 5 dfdx(nmax),dfdy(5,nmax) parameter (safe1 = 0.25d0, safe2 = 0.7d0, redmax=1.0d-5, 1 redmin = 0.7d0, tiny = 1.0d-30, scalmx = 0.1d0) data first/.true./, epsold/-1.0d0/, nvold/-1/ data nseq /2, 6, 10, 14, 22, 34, 50, 70/ c..assume that the independent variable is not explicit in the odes data dfdx/nmax*0.0d0/ c..a new tolerance or a new number , so reinitialize if (eps .ne. epsold .or. nv .ne. nvold) then hnext = -1.0e29 xnew = -1.0e29 eps1 = safe1 * eps c..compute the work coefficients a_k a(1) = nseq(1) + 1 do k=1,kmaxx a(k+1) = a(k) + nseq(k+1) enddo c..compute alf(k,q) do iq=2,kmaxx do k=1,iq-1 alf(k,iq) = eps1**((a(k+1) - a(iq+1)) / 1 ((a(iq+1) - a(1) + 1.0d0) * (2*k + 1))) enddo enddo epsold = eps nvold = nv c..add cost of jacobians to work coefficients a(1) = nv + a(1) do k=1,kmaxx a(k+1) = a(k) + nseq(k+1) enddo c..determine optimal row number for convergence do kopt=2,kmaxx-1 if (a(kopt+1) .gt. a(kopt)*alf(kopt-1,kopt)) go to 01 enddo 01 kmax = kopt end if c..save the starting values h = htry do i=1,nv ysav(i) = y(i) enddo c..get the analytic pentadiagonal jacobian call jakob(x,y,dfdy,nv,nmax) c..a new stepsize or a new integration, re-establish the order window if (h .ne. hnext .or. x .ne. xnew) then first = .true. kopt = kmax end if reduct = .false. c..evaluate the sequence of semi implicit midpoint rules 02 do k=1,kmax xnew = x + h if (xnew .eq. x) stop 'stepsize too small in routine stiffbs' call simpr3(ysav,dydx,dfdx,dfdy,nmax,nv,x,h,nseq(k),yseq,derivs) xest = (h/nseq(k))**2 call pzextr(k,xest,yseq,y,yerr,nv) c..compute normalized error estimate if (k .ne. 1) then errmax = tiny do i=1,nv errmax = max(errmax,abs(yerr(i)/yscal(i))) enddo errmax = errmax/eps km = k - 1 err(km) = (errmax/safe1)**(1.0d0/(2*km+1)) end if c..in order window if (k .ne. 1 .and. (k .ge. kopt-1 .or. first)) then c..converged if (errmax .lt. 1.0) go to 04 c..possible step size reductions if (k .eq. kmax .or. k .eq. kopt + 1) then red = safe2/err(km) go to 3 else if (k .eq. kopt) then if (alf(kopt-1,kopt) .lt. err(km)) then red = 1.0d0/err(km) go to 03 end if else if (kopt .eq. kmax) then if (alf(km,kmax-1) .lt. err(km)) then red = alf(km,kmax-1) * safe2/err(km) go to 03 end if else if (alf(km,kopt) .lt. err(km)) then red = alf(km,kopt-1)/err(km) go to 03 end if end if c..back for another sequence of midpoints enddo c..reduce stepsize by atleast redmin and at mosr redmax 03 red = min(red,redmin) red = max(red,redmax) h = h * red reduct = .true. go to 2 c..successful step; get optimal row for convergence and corresponding stepsize 04 x = xnew hdid = h first = .false. wrkmin = 1.0e35 do kk=1,km fact = max(err(kk),scalmx) work = fact * a(kk+1) if (work .lt. wrkmin) then scale = fact wrkmin = work kopt = kk + 1 end if enddo c..check for possible order increase, but not if stepsize was just reduced hnext = h/scale if (kopt .ge. k .and. kopt .ne. kmax .and. .not.reduct) then fact = max(scale/alf(kopt-1,kopt),scalmx) if (a(kopt+1)*fact .le. wrkmin) then hnext = h/fact kopt = kopt + 1 end if end if return end subroutine simpr3(y,dydx,dfdx,dfdy,nmax,n,xs,htot,nstep,yout, 1 derivs) implicit none save c..declare external derivs integer nmax,n,nstep,nmaxx parameter (nmaxx=500) integer i,j,nn double precision y(n),dydx(n),dfdx(n),dfdy(5,nmax),xs,htot, 1 yout(n),h,x,a(5,nmaxx),del(nmaxx), 2 ytemp(nmaxx),aa(nmaxx),bb(nmaxx),cc(nmaxx), 3 dd(nmaxx),ee(nmaxx) c..stepsize this trip, and make the a matrix h = htot/nstep do j=1,n do i=1,5 a(i,j) = -h*dfdy(i,j) enddo a(3,j) = a(3,j) + 1.0d0 enddo c..unpack the jacobian into pentadiagonal form do i=1,n aa(i) = a(1,i) bb(i) = a(2,i) cc(i) = a(3,i) dd(i) = a(4,i) ee(i) = a(5,i) enddo c..use yout as temporary storage; the first step do i=1,n yout(i) = h * (dydx(i) + h*dfdx(i)) enddo call pentdag(aa,bb,cc,dd,ee,yout,yout,n) do i=1,n del(i) = yout(i) ytemp(i) = y(i) + del(i) enddo x = xs + h call derivs(x,ytemp,yout) c..use yout as temporary storage; general step do nn=2,nstep do i=1,n yout(i) = h*yout(i) - del(i) enddo call pentdag(aa,bb,cc,dd,ee,yout,yout,n) do i=1,n del(i) = del(i) + 2.0d0 * yout(i) ytemp(i) = ytemp(i) + del(i) enddo x = x + h call derivs(x,ytemp,yout) enddo c..take the last step do i=1,n yout(i) = h * yout(i) - del(i) enddo call pentdag(aa,bb,cc,dd,ee,yout,yout,n) do i=1,n yout(i) = ytemp(i) + yout(i) enddo return end subroutine pzextr(iest,xest,yest,yz,dy,nv) implicit none save c..use polynomial extrapolation to evaluate nv functions at x=0 by fitting c..a polynomial to a sequence of estimates with progressively smaller values c..x=xest, and corresponding function vectors yest(1:nv). the call is number c..iest in the sequence of calls. extrapolated function values are output as c..yz(1:nv), and their estimated error is output as dy(1:nv) c.. c..declare integer iest,nv,j,k1,nmax,imax parameter (nmax=500, imax=13) double precision xest,dy(nv),yest(nv),yz(nv),delta,f1,f2,q, 1 d(nmax),qcol(nmax,imax),x(imax) c..save current independent variables x(iest) = xest do j=1,nv dy(j) = yest(j) yz(j) = yest(j) enddo c..stote first estimate in first column if (iest .eq. 1) then do j=1,nv qcol(j,1) = yest(j) enddo else do j=1,nv d(j) = yest(j) enddo do k1=1,iest-1 delta = 1.0d0/(x(iest-k1) - xest) f1 = xest * delta f2 = x(iest-k1) * delta c..propagate tableu 1 diagonal more do j=1,nv q = qcol(j,k1) qcol(j,k1) = dy(j) delta = d(j) - q dy(j) = f1*delta d(j) = f2*delta yz(j) = yz(j) + dy(j) enddo enddo do j=1,nv qcol(j,iest) = dy(j) enddo end if return end c..the next routine solves a pentadiagonal system of linear equations subroutine pentdag(a,b,c,d,e,f,u,n) implicit none save c..solves for a vector u of length n in the pentadiagonal linear system c.. a_i u_(i-2) + b_i u_(i-1) + c_i u_i + d_i u_(i+1) + e_i u_(i+2) = f_i c..input are the a, b, c, d, e, and f and they are not modified c..in its clearest incarnation, this algorithm uses three storage arrays c..called p, q and r. here, the solution vector u is used for r, cutting c..the extra storage down to two arrays. c..the vectors f and u can be input as the same array c..declare integer n,nmax,i parameter (nmax=500) double precision a(n),b(n),c(n),d(n),e(n),f(n),u(n),p(nmax), 1 q(nmax),bet,den c..initialize elimination and backsubstitution arrays if (c(1) .eq. 0.0) stop 'eliminate u2 trivially' bet = 1.0d0/c(1) p(1) = -d(1) * bet q(1) = -e(1) * bet u(1) = f(1) * bet bet = c(2) + b(2)*p(1) if (bet .eq. 0.0) stop 'singular 1 in pentdag' bet = -1.0d0/bet p(2) = (d(2) + b(2)*q(1)) * bet q(2) = e(2) * bet u(2) = (b(2)*u(1) - f(2)) * bet c..reduce to upper triangular do i=3,n bet = b(i) + a(i) * p(i-2) den = c(i) + a(i)*q(i-2) + bet*p(i-1) if (den .eq. 0.0) stop 'singular 2 in pentdag' den = -1.0d0/den p(i) = (d(i) + bet*q(i-1)) * den q(i) = e(i) * den u(i) = (a(i)*u(i-2) + bet*u(i-1) - f(i)) * den enddo c..backsubstitution u(n-1) = u(n-1) + p(n-1) * u(n) do i=n-2,1,-1 u(i) = u(i) + p(i) * u(i+1) + q(i) * u(i+2) enddo return end