ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Example 3 c c c c this is an example of using MOVCOL to solve the semilinear c c parabolic PDE with blow-up solution c c c c u_t = u_{xx} + u^2 c c c c u(0,t) = 0, u(1,t) = 0, u(x,0) = 20 sin(pi*x) c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c version: Nov. 22, 1995 c c----- c set values for some basic parameters and declare arrays c----- c implicit double precision (a-h, o-z) parameter (npts = 41, npde = 1, ntouta = 2, ipmax = 4) parameter (lrw = 60 + 6*npde + (22+3*ipmax)*npts + + (55+12*ipmax)*npts*npde + + (24+12*ipmax)*npts*npde**2, + liw = 20 + npts + 2*npts*npde) dimension touta(ntouta), par(20), rwork(lrw), iwork(liw) real etime, timaray(2), tcpu c integer p common /com1/iout common /com2/p iout = 1 p = 5 c c----- c set values for some basic parameters c----- c atol = 1.d-6 rtol = 1.d-6 c c prescribe output times. the physical pdes are integrated c in time from touta(1) to touta(ntouta) c touta(1) = 0.00d0 touta(2) = 10.d0 c c set values for par. c the default values are that c par(1) := job = 2 c par(3) := mmpde = 4 c par(4) := tau = 1.e-4 c par(5) := gamma = 1 c par(6) := ip = 2 c par(7) := stpmax = no limit c do 10 i = 1, 20 par(i) = -1.d0 10 continue par(2) = 6 par(3) = 3 par(5) = 1.d-5 par(6) = 0 c c output the information on the screen c c c define the left and right ends of the initial physical region c par(8) = 0.d0 par(9) = 1.d0 c c----- c call MOVCOL c----- c tcpu = etime(timaray) call MOVCOL(npde, npts, atol, rtol, touta, ntouta, par, + iflag, rwork, lrw, iwork, liw) tcpu = etime(timaray) - tcpu c c----- c output cpu times c----- write(6, 99995) tcpu c stop 99995 format(/' **** cpu time for mesh movement: ', f12.2/) end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c define the physical problem c c subroutines: DEFPDE, DEFBCP, DEFIVS, DEFMSH c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine DEFPDE(npde, index, t, x, u, ux, ut, uxt, fg) c c----- c define the physical PDEs, where c fg := F(t, x, u, ux, ut, uxt) if index < 0 c fg := G(t, x, u, ux, ut, uxt) if index > 0 c----- c implicit double precision (a-h,o-z) dimension u(npde), ux(npde), ut(npde), uxt(npde), fg(npde) c integer p common /com2/p if (index.lt.0) fg(1) = ut(1) - u(1)**p if (index.gt.0) fg(1) = ux(1) return end subroutine DEFBCP(npde, index, t, x, xt, u, ux, uxx, ut, + uxt, res) c c----- c define the physical boundary conditions, where c res(k) := B^L_k(t, x, xt, u, ux, uxx, ut, uxt) if index < 0 c res(k) := B^R_k(t, x, xt, u, ux, uxx, ut, uxt) if index > 0 c for k = 1, ..., npde c----- c implicit double precision (a-h,o-z) dimension u(npde), ux(npde), uxx(npde), ut(npde), uxt(npde), + res(npde) if (index.lt.0) res(1) = u(1) - 0.d0 if (index.gt.0) res(1) = u(1) - 0.d0 return end subroutine DEFIVS(npde, x, u, ux) c c----- c define the initial values c----- c implicit double precision (a-h,o-z) dimension u(npde), ux(npde) pi = dacos(-1.d0) u(1) = 20.d0*dsin(pi*x) ux(1) = 20.d0*pi*dcos(pi*x) return end subroutine DEFMSH(npts, xmesh) c c----- c define the initial mesh. this mesh is used only when job = 1 c----- c implicit double precision (a-h,o-z) dimension xmesh(npts) do 10 i = 1, npts c define a uniform mesh xmesh(i) = dfloat(i-1)/dfloat(npts-1) 10 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c define bcs and monitor function for moving mesh c c subroutines: DEFBCM, DEFMNT c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine DEFBCM(npde, index, t, x, xt, u, ux, uxx, ut, + uxt, res) c c----- c define the boundary conditions for the coordinate transformation, c where c res := b^L_M(t, x, xt, u, ux, uxx, ut, uxt) if index < 0 c res := b^R_M(t, x, xt, u, ux, uxx, ut, uxt) if index > 0 c----- c implicit double precision (a-h,o-z) dimension u(npde), ux(npde), uxx(npde), ut(npde), uxt(npde) c fixed endpoints if (index.lt.0) res = x - 0.d0 if (index.gt.0) res = x - 1.d0 return end subroutine DEFMNT(npde, t, x, u, ux, uxx, ut, uxt, fmntr) c c----- c define the monitor function c----- c implicit double precision (a-h,o-z) dimension u(npde), ux(npde), uxx(npde), ut(npde), uxt(npde) c define the arclength monitor function c integer p common /com2/p fmntr = (dabs(u(1))+1.d-15)**(p-1) return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c output solutions c c subroutine: DEFOUT c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine DEFOUT(npde, npts, t, xmesh, xmesht, u, ux, ut, tstep, + touta, ntouta, istop, index, nts) c c----- c output solutions. the time integration may be stopped if istop is set c less than zero. c the code is solving the physical pdes when index > 0 c the code is generating a mesh when index < 0 c----- c implicit double precision (a-h,o-z) dimension xmesh(npts), xmesht(npts), u(npts,npde), ux(npts,npde), + ut(npts,npde), touta(ntouta) c parameter(nrho = 5) dimension rho(nrho) common /com1/iout data rho/1.d8, 1.d9, 1.d10, 2.d10, 5.d10/ c c----- c do not output solutions when generating a mesh with respect to c the given initial values of the physical solution c----- c if (index.lt.0) return c c----- c compute the umax and position c----- c umax = 0.d0 xmin = xmesh(npts)-xmesh(1) do 10 i = 1, npts if (dabs(u(i,1)).ge.umax) then umax = dabs(u(i,1)) im = i end if if (i.le.npts-1) xmin = dmin1(xmin, xmesh(i+1)-xmesh(i)) 10 continue c h = dmin1(xmin/100.d0, 1.d-2) c i1 = 1.d0/(xmin/100.d0) + 1 c i1 = max(101, i1) c im1 = max(1, im-1) c im2 = min(npts, im+1) c do 15 i = 1, i1 c x = xmesh(im1)+(i-1.d0)*h c if (x.gt.xmesh(im2)) go to 16 c call DEFFNC(x, 1, npde, npts, xmesh, u, ux, uval, uxval) c if (dabs(uval).ge.umax) then c umax = dabs(uval) c st = x c end if c 15 continue 16 if (mod(nts,10).eq.0) then write(6, 99997) nts, t, tstep, xmin, umax write(8, 99997) nts, t, umax, dlog(umax), st do 17 i = 1, npts temp = dlog(dabs(xmesh(i)-0.5d0)+1.d-20) temp1 = dlog(umax) write(20+i, 99999) temp1, temp call flush(i+20) 17 continue endif c c----- c output the mesh c----- c if (iout.gt.nrho) then istop = -2 return endif if (umax.ge.rho(iout)) then write(9+iout, 99998) t, umax do 20 i = 2, npts-1 xi = dfloat(i-1)/(npts-1) temp = u(i,1)/umax write(9+iout, 99999) xi, xmesh(i), temp 20 continue iout = iout + 1 close(9+iout) end if if (umax.ge.rho(nrho)*2.d0) istop = -1 c return 99997 format(1x, i10, 5(1x,e12.5)) 99998 format(/'# t = ', e12.5,' umax = ', e12.5) 99999 format(6(1x,e12.5)) end