setaux.f.html
Source file:   setaux.f
Directory:   /Users/rjl/clawpack_src/clawpack_master/classic/examples/advection_2d_annulus
Converted:   Mon Feb 19 2024 at 17:58:41   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 
c     ============================================
      subroutine setaux(mbc,mx,my,xlower,ylower,dxc,dyc,
     &                  maux,aux)
c     ============================================
c
c     # set auxiliary arrays for advection on a curvilinear grid
c
c     # on input, (xc(i),yc(j)) gives uniformly spaced computational grid.
c     # on output, 
c     #   aux(i,j,1) is edge velocity at "left" boundary of grid point (i,j)
c     #   aux(i,j,2) is edge velocity at "bottom" boundary of grid point (i,j)
c     #   aux(i,j,3) = kappa  is ratio of cell area to dxc*dyc
c     
      implicit double precision (a-h,o-z)
      double precision, intent(inout), 
     &     dimension(3, 1-mbc:mx+mbc,1-mbc:my+mbc) :: aux
      dimension xccorn(5),yccorn(5),xpcorn(5),ypcorn(5)
c
c
      do 20 j=1-mbc,my+mbc
         do 20 i=1-mbc,mx+mbc
c
c           # computational points (xc,yc) are mapped to physical
c           # coordinates (xp,yp) by mapc2p:
c
c           # lower left corner:
        xccorn(1) = xlower + (i-1)*dxc
        yccorn(1) = ylower + (j-1)*dyc
        call mapc2p(xccorn(1),yccorn(1),xpcorn(1),ypcorn(1))

c           # upper left corner:
        xccorn(2) = xccorn(1)
        yccorn(2) = yccorn(1) + dyc
        call mapc2p(xccorn(2),yccorn(2),xpcorn(2),ypcorn(2))
c
c           # upper right corner:
        xccorn(3) = xccorn(1) + dxc
        yccorn(3) = yccorn(1) + dyc
        call mapc2p(xccorn(3),yccorn(3),xpcorn(3),ypcorn(3))
c
c           # lower right corner:
        xccorn(4) = xccorn(1) + dxc
        yccorn(4) = yccorn(1)
            call mapc2p(xccorn(4),yccorn(4),xpcorn(4),ypcorn(4))
c
c
c           # compute edge velocities by differencing stream function:
c
        aux(1,i,j) = (stream(xpcorn(2),ypcorn(2))
     &            - stream(xpcorn(1),ypcorn(1)))/ dyc
c
        aux(2,i,j) = -(stream(xpcorn(4),ypcorn(4))
     &            - stream(xpcorn(1),ypcorn(1)))/ dxc

c
c
c           # compute area of physical cell from four corners:

        xpcorn(5) = xpcorn(1)
        ypcorn(5) = ypcorn(1)
        area = 0.d0
        do ic=1,4
          area = area + 0.5d0 * (ypcorn(ic)+ypcorn(ic+1)) *
     &               (xpcorn(ic+1)-xpcorn(ic))
        enddo
        
        aux(3,i,j) = area / (dxc*dyc)
c
   20  continue
c
       return

       end