ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Example 2 c c c c this is an example of using MOVCOL to solve Dwyer-Sanders flame c c propagation model (two components) c c c c u_t + f(u,v) = u_{xx} c c v_t - f(u,v) = v_{xx}, 0= 2e-4) c c c c u(x,0) = 1, v(x,0) = 0.2 c c c c f(u,v) = 3.52e6 u e^{-4/v} 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 = 2, ntouta = 4, 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 c----- c open file 'soln.dat' to store the solutions c----- c common /com1/nprnt nprnt = 8 open (nprnt, file = 'soln.dat') c c----- c set values for some basic parameters c----- c atol = 1.d-7 rtol = 1.d-7 c c prescribe output times. the physical pdes are integrated c in time from touta(1) to touta(ntouta) c touta(1) = 0.d0 touta(2) = 2d-3 touta(3) = 4d-3 touta(4) = 6d-3 c c set values for par. in this example, default values are c used for par(j), j = 3, ..., 7. 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(1) = 1 c c output the information on the screen c par(2) = 6 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 close(nprnt) 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) if (index.lt.0) then temp = u(1)*3.52d6*dexp(-4.d0/u(2)) fg(1) = ut(1) + temp fg(2) = ut(2) - temp end if if (index.gt.0) then fg(1) = ux(1) fg(2) = ux(2) end if 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) then res(1) = ux(1) res(2) = ux(2) end if if (index.gt.0) then res(1) = ux(1) temp = 1.2d0 if (t.le.2d-4) temp = 0.2d0+t/2d-4 res(2) = u(2) - temp end if 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) u(1) = 1.d0 ux(1) = 0.d0 u(2) = 0.2d0 ux(2) = 0.d0 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) temp = 1.d0/dfloat(npts-1) do 10 i = 1, npts c define a uniform mesh xmesh(i) = dfloat(i-1)*temp 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 fmntr = dsqrt(1.d0+ux(1)*ux(1)+ux(2)*ux(2)) 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) common /com1/nprnt 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 output solutions and errors at t = touta(1), ..., touta(ntouta) c----- c if (mod(nts,30).eq.0) then write(nprnt, 99997) t do 10 i = 1, npts x = xmesh(i) u1 = u(i,1) u2 = u(i,2) if (dabs(u1).le.1d-10) u1 = 1d-10 if (dabs(u2).le.1d-10) u2 = 1d-10 write(nprnt, 99999) x, u1, u2 10 continue end if return 99997 format(/'# t = ', f10.4/) 99999 format(6(1x,e12.5)) end