c
c =====================================================
subroutine setaux(maxmx,maxmy,maxmz,meqn,mbc,ibx,iby,ibz,
& mx,my,mz,q,aux,maux,xc,yc,zc,dx,dy,dz,t,dt)
c =====================================================
c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h, o-z)
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc, 1-ibz*mbc:maxmz+ibz*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc, 1-ibz*mbc:maxmz+ibz*mbc)
c
c
c
c # advection
c # aux(i,j,k,1) is u velocity on left face of cell
c # aux(i,j,k,2) is v velocity on bottom face of cell
c # aux(i,j,k,3) is w velocity on back face of cell
c
c # on input xc, yc, zc contain cell centers
c
c # define the velocity field (deformation flow as in paper):
u(x,y,z) = 2.d0 * dsin(pi*x)**2 * dsin(pi2*y) * dsin(pi2*z) *
$ dcos((t+dt)*pi)
v(x,y,z) = -dsin(pi2*x) * dsin(pi*y)**2 * dsin(pi2*z) *
$ dcos((t+dt)*pi)
w(x,y,z) = -dsin(pi2*x) * dsin(pi2*y) * dsin(pi*z)**2 *
$ dcos((t+dt)*pi)
c
! u(x,y,z) = 2.d0 * dsin(pi*x)**2 * dsin(pi2*y) * dsin(pi2*z)
! v(x,y,z) = -dsin(pi2*x) * dsin(pi*y)**2 * dsin(pi2*z)
! w(x,y,z) = -dsin(pi2*x) * dsin(pi2*y) * dsin(pi*z)**2
c
pi = 4.d0*datan(1.d0)
pi2 = 2.d0*pi
c
c (xc,yc,zc) lower left corner of cell (1-mbc,1-mbc,1-mbc)
c
c (xc+(i-(1-mbc)+.5)*dx, yc+(j-(1-mbc)+.5)*dy, zc+(k-(1-mbc)+.5)*dz)
c center coords of cell (1-mbc,1-mbc,1-mbc)
c
do 10 k=1-ibz*mbc,mz+ibz*mbc
do 10 j=1-iby*mbc,my+iby*mbc
do 10 i=1-ibx*mbc,mx+ibx*mbc
aux(1,i,j,k) = u(xc+(i+ibx*mbc-1.d0)*dx,
& yc+(j+iby*mbc-.5d0)*dy,
& zc+(k+ibz*mbc-.5d0)*dz)
aux(2,i,j,k) = v(xc+(i+ibx*mbc-.5d0)*dx,
& yc+(j+iby*mbc-1.d0)*dy,
& zc+(k+ibz*mbc-.5d0)*dz)
aux(3,i,j,k) = w(xc+(i+ibx*mbc-.5d0)*dx,
& yc+(j+iby*mbc-.5d0)*dy,
& zc+(k+ibz*mbc-1.d0)*dz)
10 continue
c
return
end