examples/advect3d/driver.f

   1       PROGRAM ADVECT3D
   2       IMPLICIT NONE
   3       INTEGER nx, ng, tt
   4
   5       parameter (nx = 256)
   6       parameter (ng = 9)
   7       parameter (tt = 1)
   8         REAL x(1-ng:nx+ng), y(1-ng:nx+ng), z(1-ng:nx+ng)
   9         REAL dx, dis
  10         REAL pi, Ad1, Ad2
  11         REAL  A0(1-ng:nx+1+ng, 1-ng:nx+1+ng)
  12         REAL  An(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  13         REAL  Ao(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  14         REAL  Ai(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  15         REAL   u(1-ng:nx  +ng, 1-ng:nx+1+ng,1-ng:nx+  ng)
  16         REAL   v(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+1+ng)
  17         REAL   w(1-ng:nx+1+ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  18         REAL  u0(1-ng:nx  +ng, 1-ng:nx+1+ng,1-ng:nx+  ng)
  19         REAL  v0(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+1+ng)
  20         REAL  w0(1-ng:nx+1+ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  21         REAL  ut(1-ng:nx  +ng, 1-ng:nx+1+ng,1-ng:nx+  ng)
  22         REAL  vt(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+1+ng)
  23         REAL  wt(1-ng:nx+1+ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  24         REAL  uh(1-ng:nx  +ng, 1-ng:nx+1+ng,1-ng:nx+  ng)
  25         REAL  vh(1-ng:nx  +ng, 1-ng:nx+  ng,1-ng:nx+1+ng)
  26         REAL  wh(1-ng:nx+1+ng, 1-ng:nx+  ng,1-ng:nx+  ng)
  27         REAL X0, Y0, Z0, rad
  28         REAL dt
  29         REAL maxi
  30         INTEGER i,j,k,t,m, imgT
  31         INTEGER iloc, jloc, kloc
  32         real xl, yl, zl, third, half
  33         PARAMETER (dt = 0.005)
  34         PARAMETER (pi = 3.1415926)
  35         PARAMETER (X0 = 0.00)
  36         PARAMETER (Y0 = 0.00)
  37         PARAMETER (Z0 = 0.00)
  38         PARAMETER (rad = 0.15)
  39         imgT = Tt/4
  40         dx = 1./REAL(nx)
  41         maxi = 0.0
  42         third = 1./3.
  43         half = 0.5
  44
  45         do i=1-ng,nx+ng
  46            x(i) = -.5 + dx*REAL(i)
  47            y(i) = -.5 + dx*REAL(i)
  48            z(i) = -.5 + dx*REAL(i)
  49         end do
  50         xl = REAL(nx) * dx * .5
  51         yl = REAL(nx) * dx * .5
  52         zl = REAL(nx) * dx * .5
  53
  54         do j=1,nx
  55           do i=1,nx
  56             do k=1,nx
  57               u0(k,i,j) =  1.0 * (x(k) - 0.5*dx)
  58               v0(k,i,j) = 0.
  59               w0(k,i,j) = -1.0 * (z(i) - 0.5*dx)
  60             end do
  61           end do
  62         end do
  63
  64         Ad1 = 0.
  65          do j=1,nx
  66            do i=1,nx
  67              do k=1,nx
  68              dis = SQRT( (x(i) -.5*dx - X0)**2 +
  69      +                   (y(j) -.5*dx - Y0)**2 +
  70      +                   (z(k) -.5*dx - Z0)**2)
  71                 Ao(k,i,j) = (20. * 0.5 * (1. + COS(dis*pi/rad)))*
  72      +                    (.5+sign(.5,rad-dis))
  73               Ai(k,i,j) = Ao(k,i,j)
  74               An(k,i,j) = 0.0
  75               Ad1 = Ad1 + Ao(k,i,j)
  76             end do
  77           end do
  78         end do
  79         do i=1-ng,nx+ng
  80            do k=1-ng,nx+ng
  81               do m = 0, 1-ng
  82                  Ao(k,i,m) = Ao(k,i,1)
  83               enddo
  84               do m = nx + 1, nx + ng
  85                  Ao(k,i,m) = Ao(k,i,nx)
  86               enddo
  87            enddo
  88         enddo
  89         do j=1-ng,nx+ng
  90            do k=1-ng,nx+ng
  91            do m = 0, 1-ng
  92               Ao(k,m,j) = Ao(k,1 ,j)
  93            enddo
  94            do m = nx+1, nx+ng
  95               Ao(k,m,j) = Ao(k,nx,j)
  96            enddo
  97         enddo
  98
  99            do i=1,nx
 100            do k=1-ng,nx+ng
 101               maxi = MAX(maxi,Ao(k,i,j))
 102               if(Ao(k,i,j) .EQ. maxi) then
 103                  iloc = i
 104                  jloc = j
 105                  kloc = k
 106               endif
 107            end do
 108            end do
 109         end do
 110          write(*,*) maxi, 'max val start ',iloc,jloc,kloc
 111         maxi = 0.0
 112         print *, '( nx, ng, tt ) = ', '(', nx, ',', ng, ',', tt, ')'
 113         write (*,*) Tt
 114         do t=1, 5
 115             call twostages(Ao, An, u0, v0, w0, u0, v0, w0,
 116      &                      dt, dx, dx, dx, nx, nx, nx, ng)
 117
 118             do j = 1, nx
 119                do i = 1, nx
 120                   do k = 1, nx
 121                      Ao(k,i,j) = An(k,i,j)
 122                   enddo
 123                enddo
 124             enddo
 125         end do
 126         write(*,*) An(kloc,jloc,iloc), ' val @ finish for initial max '
 127         Ad2 = 0.
 128         do j=1,nx
 129           do i=1,nx
 130            do k=1-ng,nx+ng
 131               maxi = MAX(maxi,Ao(k,i,j))
 132               if(Ao(k,i,j) .EQ. maxi) then
 133                  iloc = i
 134                  jloc = j
 135                  kloc = k
 136               endif
 137               Ad2 = Ad2 + An(k,i,j)
 138            end do
 139           end do
 140         end do
 141         write(*,*) maxi, 'max val finish ',iloc,jloc
 142         maxi = 0.0
 143         do j = 1, nx
 144            do i = 1, nx
 145               do k = 1, nx
 146                  Ai(k, i, j) = An(k, i, k) - Ai(k, i, j)
 147               enddo
 148            enddo
 149         enddo
 150         do j=1,nx
 151           do i=1,nx
 152             do k=1-ng,nx+ng
 153               maxi = MAX(maxi,Ai(k,i,j))
 154               if(Ai(k,i,j) .EQ. maxi) then
 155                  iloc = i
 156                  jloc = j
 157                  kloc = k
 158               endif
 159            end do
 160           end do
 161         end do
 162          write(*,*) 'max error finish: ',maxi,iloc,jloc
 163          write(*,*) 'Starting amount:  ',Ad1
 164          write(*,*) 'Finishing amount: ',Ad2
 165          write(*,*) 'Percent conserved: ',Ad2/Ad1 * 100.
 166
 167        END
 168