C$TEST DTTGR2
c main program
common /cstak/ ds
double precision ds(350000)
external handle, bc, af
integer ndx, ndy, istkgt, is(1000), iu, ix
integer iy, nu, kx, nx, ky, ny
integer idumb
real errpar(2), rs(1000)
logical ls(1000)
complex cs(500)
double precision tstart, dt, lx, ly, rx, ry
double precision ws(500), tstop
equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c to solve two coupled, nonlinear heat equations.
c u1 sub t = div . ( u1x, u1y ) - u1*u2 + g1
c u2 sub t = div . ( u2x, u2y ) - u1*u2 + g2
c the port library stack and its aliases.
c initialize the port library stack length.
call istkin(350000, 4)
call enter(1)
nu = 2
lx = 0
rx = 1
ly = 0
ry = 1
kx = 4
ky = 4
ndx = 3
ndy = 3
tstart = 0
tstop = 1
dt = 1e-2
errpar(1) = 1e-2
errpar(2) = 1e-4
c uniform grid.
ix = idumb(lx, rx, ndx, kx, nx)
c uniform grid.
iy = idumb(ly, ry, ndy, ky, ny)
c space for the solution.
iu = istkgt(nu*(nx-kx)*(ny-ky), 4)
call setd(nu*(nx-kx)*(ny-ky), 1d0, ws(iu))
call dttgr(ws(iu), nu, kx, ws(ix), nx, ky, ws(iy), ny, tstart,
1 tstop, dt, af, bc, errpar, handle)
call leave
call wrapup
stop
end
subroutine af(t, x, nx, y, ny, nu, u, ut, ux, uy, uxt, uyt
1 , a, au, aut, aux, auy, auxt, auyt, f, fu, fut, fux, fuy, fuxt,
2 fuyt)
integer nu, nx, ny
double precision t, x(nx), y(ny), u(nx, ny, nu), ut(nx, ny, nu),
1 ux(nx, ny, nu)
double precision uy(nx, ny, nu), uxt(nx, ny, nu), uyt(nx, ny, nu),
1 a(nx, ny, nu, 2), au(nx, ny, nu, nu, 2), aut(nx, ny, nu, nu, 2)
double precision aux(nx, ny, nu, nu, 2), auy(nx, ny, nu, nu, 2),
1 auxt(nx, ny, nu, nu, 2), auyt(nx, ny, nu, nu, 2), f(nx, ny, nu)
2 , fu(nx, ny, nu, nu)
double precision fut(nx, ny, nu, nu), fux(nx, ny, nu, nu), fuy(nx,
1 ny, nu, nu), fuxt(nx, ny, nu, nu), fuyt(nx, ny, nu, nu)
integer p, q
double precision dexp
do 2 q = 1, ny
do 1 p = 1, nx
a(p, q, 1, 1) = ux(p, q, 1)
aux(p, q, 1, 1, 1) = 1
a(p, q, 1, 2) = uy(p, q, 1)
auy(p, q, 1, 1, 2) = 1
f(p, q, 1) = ut(p, q, 1)+u(p, q, 1)*u(p, q, 2)
fu(p, q, 1, 1) = u(p, q, 2)
fu(p, q, 1, 2) = u(p, q, 1)
fut(p, q, 1, 1) = 1
a(p, q, 2, 1) = ux(p, q, 2)
aux(p, q, 2, 2, 1) = 1
a(p, q, 2, 2) = uy(p, q, 2)
auy(p, q, 2, 2, 2) = 1
f(p, q, 2) = ut(p, q, 2)+u(p, q, 1)*u(p, q, 2)
fu(p, q, 2, 1) = u(p, q, 2)
fu(p, q, 2, 2) = u(p, q, 1)
fut(p, q, 2, 2) = 1
f(p, q, 1) = f(p, q, 1)-(dexp(t*(x(p)-y(q)))*(x(p)-y(q)-2d0*
1 t*t)+1d0)
f(p, q, 2) = f(p, q, 2)-(dexp(t*(y(q)-x(p)))*(y(q)-x(p)-2d0*
1 t*t)+1d0)
1 continue
2 continue
return
end
subroutine bc(t, x, nx, y, ny, lx, rx, ly, ry, u, ut, ux,
1 uy, uxt, uyt, nu, b, bu, but, bux, buy, buxt, buyt)
integer nu, nx, ny
double precision t, x(nx), y(ny), lx, rx, ly
double precision ry, u(nx, ny, nu), ut(nx, ny, nu), ux(nx, ny, nu)
1 , uy(nx, ny, nu), uxt(nx, ny, nu)
double precision uyt(nx, ny, nu), b(nx, ny, nu), bu(nx, ny, nu,
1 nu), but(nx, ny, nu, nu), bux(nx, ny, nu, nu), buy(nx, ny, nu
2 , nu)
double precision buxt(nx, ny, nu, nu), buyt(nx, ny, nu, nu)
integer i, j
double precision dexp
do 2 j = 1, ny
do 1 i = 1, nx
bu(i, j, 1, 1) = 1
b(i, j, 1) = u(i, j, 1)-dexp(t*(x(i)-y(j)))
bu(i, j, 2, 2) = 1
b(i, j, 2) = u(i, j, 2)-dexp(t*(y(j)-x(i)))
1 continue
2 continue
return
end
subroutine handle(t0, u0, t, u, nv, dt, tstop)
integer nv
double precision t0, u0(nv), t, u(nv), dt, tstop
common /cstak/ ds
double precision ds(500)
common /d7tgrp/ errpar, nu, mxq, myq
integer nu, mxq, myq
real errpar(2)
common /d7tgrm/ kx, ix, nx, ky, iy, ny
integer kx, ix, nx, ky, iy, ny
integer ifa, ita(2), ixa(2), nta(2), nxa(2), idlumd
integer ixs, iys, nxs, nys, istkgt, i
integer j, iewe, ka(2), ma(2), is(1000), i1mach
real rs(1000)
logical ls(1000)
complex cs(500)
double precision dabs, erru, dmax1, ws(500)
integer temp, temp1, temp2
equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c the port library stack and its aliases.
if (t0 .ne. t) goto 2
write (6, 1) t
1 format (16h restart for t =, 1pe10.2)
return
2 call enter(1)
c find the error in the solution at 2*kx * 2*ky points / mesh rectangle.
c x search grid.
ixs = idlumd(ws(ix), nx, 2*kx, nxs)
c y search grid.
iys = idlumd(ws(iy), ny, 2*ky, nys)
c u search grid values.
iewe = istkgt(nu*nxs*nys, 4)
c the exact solution.
call ewe(t, ws(ixs), nxs, ws(iys), nys, ws(iewe), nu)
ka(1) = kx
ka(2) = ky
ita(1) = ix
ita(2) = iy
nta(1) = nx
nta(2) = ny
ixa(1) = ixs
ixa(2) = iys
nxa(1) = nxs
nxa(2) = nys
ma(1) = 0
c get solution.
ma(2) = 0
c approximate solution values.
ifa = istkgt(nxs*nys, 4)
do 5 j = 1, nu
c evaluate them.
temp = (j-1)*(nx-kx)*(ny-ky)
call dtsd1(2, ka, ws, ita, nta, u(temp+1), ws, ixa, nxa, ma,
1 ws(ifa))
c error in solution values.
erru = 0
temp = nxs*nys
do 3 i = 1, temp
temp2 = iewe+i-1+(j-1)*nxs*nys
temp1 = ifa+i
erru = dmax1(erru, dabs(ws(temp2)-ws(temp1-1)))
3 continue
temp = i1mach(2)
write (temp, 4) t, j, erru
4 format (14h error in u(.,, 1pe10.2, 1h,, i2, 3h) =, 1pe10.2)
5 continue
call leave
return
end
subroutine ewe(t, x, nx, y, ny, u, nu)
integer nu, nx, ny
double precision t, x(nx), y(ny), u(nx, ny, nu)
integer i, j, p
real float
double precision dble, dexp
c the exact solution.
do 3 p = 1, nu
do 2 i = 1, nx
do 1 j = 1, ny
u(i, j, p) = dexp(dble(float((-1)**(p+1)))*t*(x(i)-y(j)))
1 continue
2 continue
3 continue
return
end
.