(15) | ||

(16) |

The spatial tiling we have chosen gives us two independent locations for
storing density ,one atop the horizontal velocity *u*,
the other atop the vertical velocity *v*.
It may seem non physical to have one density
for vertical accelerations and another for horizontal ones,
but it is worth maintaining the two densities
as separate entities.
Waves arise from many physical paradigms besides the
acoustical one and some of these
may require that we maintain this possibility to introduce anisotropy.
Another way to introduce anisotropy in the present codes
is to have two compressibilities *K* in the pressure subroutine.
These extra parameters might also be helpful if
we were trying to improve the numerical representation
of the differential operators.
Subroutine `velocity()` uses
`rhou(,)` for at *u* and
`rhov(,)` for at *v*.
Equation (15) in discrete form

(17) |

# Operator of momentum change from gradient of pressure # subroutine velocity( adj, add, rhou,rhov, q,nx,ny, r ) integer adj, add, nx,ny, x, y, p, u, v real rhou(nx,ny), rhov(nx,ny), q(nx,ny,3), r(nx,ny,3) call adjnull( adj, add, q,nx*ny*3, r,nx*ny*3) p=1; u=2; v=3 do x= 2, nx-1 { do y= 2, ny-1 { if( adj == 0) { r(x,y,u)= r(x,y,u) + q(x,y,u) - ( q(x+1,y ,p) - q(x,y,p)) * rhou(x,y) r(x,y,v)= r(x,y,v) + q(x,y,v) - ( q(x ,y+1,p) - q(x,y,p)) * rhov(x,y) r(x,y,p)= r(x,y,p) + q(x,y,p) } else { q(x ,y ,u)= q(x ,y ,u) + r(x,y,u) q(x ,y ,p)= q(x ,y ,p) + r(x,y,u) * rhou(x,y) + r(x,y,p) q(x+1,y ,p)= q(x+1,y ,p) - r(x,y,u) * rhou(x,y) # q(x ,y ,v)= q(x ,y ,v) + r(x,y,v) q(x ,y ,p)= q(x ,y ,p) + r(x,y,v) * rhov(x,y) q(x ,y+1,p)= q(x ,y+1,p) - r(x,y,v) * rhov(x,y) } }} return; end

This time we pass through the pressure unchanged so we also see an added pressure term in the adjoint.

11/12/1997