C$TEST TTGR5
c main program
common /cstak/ ds
double precision ds(350000)
external handle, bc, af
integer ndx, ndy, istkgt, i, is(1000), iu
integer ix, iy, nu, kx, nx, ky
integer ny
real errpar(2), tstart, dt, lx, ly, rx
real ry, ws(1000), rs(1000), float, tstop
logical ls(1000)
complex cs(500)
integer temp, temp1
equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c to solve laplaces equation with real ( z*log(z) ) as solution.
c the port library stack and its aliases.
c initialize the port library stack length.
call istkin(350000, 4)
call enter(1)
nu = 1
lx = 0
rx = 1
ly = 0
ry = 1
kx = 4
ky = 4
ndx = 2
ndy = 2
tstart = 0
tstop = 1
dt = 1
errpar(1) = 1e-2
errpar(2) = 1e-4
nx = ndx+2*(kx-1)
c space for x mesh.
ix = istkgt(nx, 3)
do 1 i = 1, kx
temp = ix+i
ws(temp-1) = 0
temp = ix+nx-i
ws(temp) = rx
1 continue
c 0 and rx mult = kx.
temp = ndx-1
do 2 i = 1, temp
temp1 = ix+kx-2+i
ws(temp1) = rx*(float(i-1)/(float(ndx)-1e0))**kx
2 continue
ny = ndy+2*(ky-1)
c space for y mesh.
iy = istkgt(ny, 3)
do 3 i = 1, ky
temp = iy+i
ws(temp-1) = 0
temp = iy+ny-i
ws(temp) = ry
3 continue
c 0 and ry mult = ky.
temp = ndy-1
do 4 i = 1, temp
temp1 = iy+ky-2+i
ws(temp1) = ry*(float(i-1)/(float(ndy)-1e0))**ky
4 continue
c space for the solution.
iu = istkgt(nu*(nx-kx)*(ny-ky), 3)
call setr(nu*(nx-kx)*(ny-ky), 0e0, ws(iu))
call ttgr(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, xi, nx, yi, ny, nu, u, ut, ux, uy, uxt,
1 uyt, a, au, aut, aux, auy, auxt, auyt, f, fu, fut, fux, fuy,
2 fuxt, fuyt)
integer nu, nx, ny
real t, xi(nx), yi(ny), u(nx, ny, nu), ut(nx, ny, nu), ux(nx, ny
1 , nu)
real uy(nx, ny, nu), uxt(nx, ny, nu), uyt(nx, ny, nu), a(nx, ny,
1 nu, 2), au(nx, ny, nu, nu, 2), aut(nx, ny, nu, nu, 2)
real aux(nx, ny, nu, nu, 2), auy(nx, ny, nu, nu, 2), auxt(nx, ny
1 , nu, nu, 2), auyt(nx, ny, nu, nu, 2), f(nx, ny, nu), fu(nx,
2 ny, nu, nu)
real fut(nx, ny, nu, nu), fux(nx, ny, nu, nu), fuy(nx, ny, nu, nu)
1 , fuxt(nx, ny, nu, nu), fuyt(nx, ny, nu, nu)
integer i, p, q
do 3 i = 1, nu
do 2 q = 1, ny
do 1 p = 1, nx
a(p, q, i, 1) = ux(p, q, i)
a(p, q, i, 2) = uy(p, q, i)
aux(p, q, i, i, 1) = 1
auy(p, q, i, i, 2) = 1
1 continue
2 continue
3 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
real t, x(nx), y(ny), lx, rx, ly
real ry, u(nx, ny, nu), ut(nx, ny, nu), ux(nx, ny, nu), uy(nx, ny,
1 nu), uxt(nx, ny, nu)
real uyt(nx, ny, nu), b(nx, ny, nu), bu(nx, ny, nu, nu), but(nx,
1 ny, nu, nu), bux(nx, ny, nu, nu), buy(nx, ny, nu, nu)
real buxt(nx, ny, nu, nu), buyt(nx, ny, nu, nu)
integer i, j
real cos, sin, r, alog, atan, sqrt
real theta
do 6 j = 1, ny
do 5 i = 1, nx
if (y(j) .ne. ly) goto 1
b(i, j, 1) = uy(i, j, 1)
c neumann data on bottom.
buy(i, j, 1, 1) = 1
goto 4
1 r = sqrt(x(i)**2+y(j)**2)
c dirichlet data.
if (x(i) .le. 0.) goto 2
theta = atan(y(j)/x(i))
goto 3
2 theta = 2.*atan(1e0)
3 b(i, j, 1) = u(i, j, 1)-r*(cos(theta)*alog(r)-theta*sin(
1 theta))
bu(i, j, 1, 1) = 1
4 continue
5 continue
6 continue
return
end
subroutine handle(t0, u0, t, u, nv, dt, tstop)
integer nv
real t0, u0(nv), t, u(nv), dt, tstop
common /a7tgrp/ errpar, nu, mxq, myq
integer nu, mxq, myq
real errpar(2)
common /a7tgrm/ kx, ix, nx, ky, iy, ny
integer kx, ix, nx, ky, iy, ny
if (t0 .ne. t) goto 2
write (6, 1) t
1 format (16h restart for t =, 1pe10.2)
return
2 call gerr(kx, ix, nx, ky, iy, ny, u, nu, t)
return
end
subroutine gerr(kx, ix, nx, ky, iy, ny, u, nu, t)
integer kx, ix, nx, ky, iy, ny
integer nu
real u(1), t
common /cstak/ ds
double precision ds(500)
integer ifa, ita(2), ixa(2), nta(2), nxa(2), ixs
integer iys, nxs, nys, istkgt, i, iewe
integer ka(2), ma(2), is(1000), ilumd, i1mach
real abs, erru, amax1, rs(1000), ws(1000)
logical ls(1000)
complex cs(500)
integer temp, temp1, temp2
equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c to get and print the error at each time-step.
c u(nx-kx,ny-ky,nu).
c the port library stack and its aliases.
call enter(1)
c find the error in the solution at 2*kx * 2*ky points / mesh rectangle.
c x search grid.
ixs = ilumd(ws(ix), nx, 2*kx, nxs)
c y search grid.
iys = ilumd(ws(iy), ny, 2*ky, nys)
c u search grid values.
iewe = istkgt(nxs*nys, 3)
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, 3)
c evaluate them.
call tsd1(2, ka, ws, ita, nta, u, ws, ixa, nxa, ma, ws(ifa))
c error in solution values.
erru = 0
temp = nxs*nys
do 1 i = 1, temp
temp2 = iewe+i
temp1 = ifa+i
erru = amax1(erru, abs(ws(temp2-1)-ws(temp1-1)))
1 continue
temp = i1mach(2)
write (temp, 2) t, erru
2 format (14h error in u(.,, 1pe10.2, 3h) =, 1pe10.2)
call leave
return
end
subroutine ewe(t, x, nx, y, ny, u, nu)
integer nu, nx, ny
real t, x(nx), y(ny), u(nx, ny, nu)
integer i, j, p
real cos, sin, r, alog, atan, sqrt
real theta
c the exact solution.
do 7 p = 1, nu
do 6 i = 1, nx
do 5 j = 1, ny
r = sqrt(x(i)**2+y(j)**2)
if (x(i) .le. 0.) goto 1
theta = atan(y(j)/x(i))
goto 2
1 theta = 2.*atan(1e0)
2 if (r .le. 0.) goto 3
u(i, j, p) = r*(cos(theta)*alog(r)-theta*sin(theta))
goto 4
3 u(i, j, p) = 0
4 continue
5 continue
6 continue
7 continue
return
end
.