ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Example 1 c c c c this is an example of using MOVCOL to solve Burgers' equation with c c an exact solution : c c c c u_t = [ 10^(-4) u_x - u^2/2 ]_x, 0 < x < 1, 0 < t <= 1 c c u(0, t) = exact(0, t), u(1, t) = exact(1,t) c c u(x, 0) = exact(x, 0) 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 = 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-3 rtol = 1.d-5 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) = 0.25d0 touta(3) = 0.55d0 touta(4) = 1.00d0 c c set values for par. in this example, default values are c used for par(j), j = 1, ..., 7, except for par(2). 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 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) eps = 1.d-4 if (index.lt.0) fg(1) = ut(1) if (index.gt.0) fg(1) = eps*ux(1)-u(1)**2/2.d0 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) - exact(0.d0,t) if (index.gt.0) res(1) = u(1) - exact(1.d0,t) 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) = exact(x, 0.d0) ux(1) = dexact(x, 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) 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 fmntr = dsqrt(1.d0+ux(1)**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 do 20 n = 1, ntouta if (dabs(t-touta(n)).le.1.d-15) then write(nprnt, 99997) t do 10 i = 1, npts x = xmesh(i) c call DEFFNC(x, 1, npde, npts, xmesh, u, ux, uval, uxval) uval = u(i,1) uexac = exact(x,t) write(nprnt, 99999) x, uval, uexac, dabs(uval-uexac) 10 continue end if 20 continue return 99997 format(/'# t = ', f10.4/) 99999 format(6(1x,e12.5)) end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c define the exact solution c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function exact(x,t) c c----- c define the exact solution for the current problem c----- c implicit double precision (a-h,o-z) eps = 1.d-4 a = - 0.05d0/eps*(x-0.5d0+4.95d0*t) b = - 0.25d0/eps*(x-0.5d0+0.75d0*t) c = - 0.50d0/eps*(x-0.375d0) temp = dmax1(a, b, c) ea = 0.d0 eb = 0.d0 ec = 0.d0 if (a-temp.ge.-35.d0) ea = dexp(a-temp) if (b-temp.ge.-35.d0) eb = dexp(b-temp) if (c-temp.ge.-35.d0) ec = dexp(c-temp) exact = (0.1d0*ea+0.5d0*eb+ec)/(ea+eb+ec) return end double precision function dexact(x,t) c c----- c define the derivative of the exact solution for the current c problem c----- c implicit double precision (a-h,o-z) eps = 1.d-4 a = - 0.05d0/eps*(x-0.5d0+4.95d0*t) b = - 0.25d0/eps*(x-0.5d0+0.75d0*t) c = - 0.50d0/eps*(x-0.375d0) temp = dmax1(a, b, c) r1 = 0.d0 r2 = 0.d0 r3 = 0.d0 if (a-temp.ge.-35.d0) r1 = dexp(a-temp) if (b-temp.ge.-35.d0) r2 = dexp(b-temp) if (c-temp.ge.-35.d0) r3 = dexp(c-temp) r1x = - 0.05d0*r1/eps r2x = - 0.25d0*r2/eps r3x = - 0.50d0*r3/eps dexact = ((-0.4d0*r2-0.9d0*r3)*r1x+(0.4d0*r1-0.5d0*r3)*r2x ++(0.9d0*r1+0.5d0*r2)*r3x)/(r1+r2+r3)**2 return end