examples/smpi/NAS/BT/initialize.f

   1 c---------------------------------------------------------------------
   2 c---------------------------------------------------------------------
   3
   4       subroutine  initialize
   5
   6 c---------------------------------------------------------------------
   7 c---------------------------------------------------------------------
   8
   9 c---------------------------------------------------------------------
  10 c     This subroutine initializes the field variable u using
  11 c     tri-linear transfinite interpolation of the boundary values
  12 c---------------------------------------------------------------------
  13
  14       include 'header.h'
  15
  16       integer c, i, j, k, m, ii, jj, kk, ix, iy, iz
  17       double precision  xi, eta, zeta, Pface(5,3,2), Pxi, Peta,
  18      >     Pzeta, temp(5)
  19
  20 c---------------------------------------------------------------------
  21 c  Later (in compute_rhs) we compute 1/u for every element. A few of
  22 c  the corner elements are not used, but it convenient (and faster)
  23 c  to compute the whole thing with a simple loop. Make sure those
  24 c  values are nonzero by initializing the whole thing here.
  25 c---------------------------------------------------------------------
  26       do c = 1, ncells
  27          do kk = -1, KMAX
  28             do jj = -1, JMAX
  29                do ii = -1, IMAX
  30                   do m = 1, 5
  31                      u(m, ii, jj, kk, c) = 1.0
  32                   end do
  33                end do
  34             end do
  35          end do
  36       end do
  37 c---------------------------------------------------------------------
  38
  39
  40
  41 c---------------------------------------------------------------------
  42 c     first store the "interpolated" values everywhere on the grid
  43 c---------------------------------------------------------------------
  44       do c=1, ncells
  45          kk = 0
  46          do k = cell_low(3,c), cell_high(3,c)
  47             zeta = dble(k) * dnzm1
  48             jj = 0
  49             do j = cell_low(2,c), cell_high(2,c)
  50                eta = dble(j) * dnym1
  51                ii = 0
  52                do i = cell_low(1,c), cell_high(1,c)
  53                   xi = dble(i) * dnxm1
  54
  55                   do ix = 1, 2
  56                      call exact_solution(dble(ix-1), eta, zeta,
  57      >                    Pface(1,1,ix))
  58                   enddo
  59
  60                   do iy = 1, 2
  61                      call exact_solution(xi, dble(iy-1) , zeta,
  62      >                    Pface(1,2,iy))
  63                   enddo
  64
  65                   do iz = 1, 2
  66                      call exact_solution(xi, eta, dble(iz-1),
  67      >                    Pface(1,3,iz))
  68                   enddo
  69
  70                   do m = 1, 5
  71                      Pxi   = xi   * Pface(m,1,2) +
  72      >                    (1.0d0-xi)   * Pface(m,1,1)
  73                      Peta  = eta  * Pface(m,2,2) +
  74      >                    (1.0d0-eta)  * Pface(m,2,1)
  75                      Pzeta = zeta * Pface(m,3,2) +
  76      >                    (1.0d0-zeta) * Pface(m,3,1)
  77
  78                      u(m,ii,jj,kk,c) = Pxi + Peta + Pzeta -
  79      >                    Pxi*Peta - Pxi*Pzeta - Peta*Pzeta +
  80      >                    Pxi*Peta*Pzeta
  81
  82                   enddo
  83                   ii = ii + 1
  84                enddo
  85                jj = jj + 1
  86             enddo
  87             kk = kk+1
  88          enddo
  89       enddo
  90
  91 c---------------------------------------------------------------------
  92 c     now store the exact values on the boundaries
  93 c---------------------------------------------------------------------
  94
  95 c---------------------------------------------------------------------
  96 c     west face
  97 c---------------------------------------------------------------------
  98       c = slice(1,1)
  99       ii = 0
 100       xi = 0.0d0
 101       kk = 0
 102       do k = cell_low(3,c), cell_high(3,c)
 103          zeta = dble(k) * dnzm1
 104          jj = 0
 105          do j = cell_low(2,c), cell_high(2,c)
 106             eta = dble(j) * dnym1
 107             call exact_solution(xi, eta, zeta, temp)
 108             do m = 1, 5
 109                u(m,ii,jj,kk,c) = temp(m)
 110             enddo
 111             jj = jj + 1
 112          enddo
 113          kk = kk + 1
 114       enddo
 115
 116 c---------------------------------------------------------------------
 117 c     east face
 118 c---------------------------------------------------------------------
 119       c  = slice(1,ncells)
 120       ii = cell_size(1,c)-1
 121       xi = 1.0d0
 122       kk = 0
 123       do k = cell_low(3,c), cell_high(3,c)
 124          zeta = dble(k) * dnzm1
 125          jj = 0
 126          do j = cell_low(2,c), cell_high(2,c)
 127             eta = dble(j) * dnym1
 128             call exact_solution(xi, eta, zeta, temp)
 129             do m = 1, 5
 130                u(m,ii,jj,kk,c) = temp(m)
 131             enddo
 132             jj = jj + 1
 133          enddo
 134          kk = kk + 1
 135       enddo
 136
 137 c---------------------------------------------------------------------
 138 c     south face
 139 c---------------------------------------------------------------------
 140       c = slice(2,1)
 141       jj = 0
 142       eta = 0.0d0
 143       kk = 0
 144       do k = cell_low(3,c), cell_high(3,c)
 145          zeta = dble(k) * dnzm1
 146          ii = 0
 147          do i = cell_low(1,c), cell_high(1,c)
 148             xi = dble(i) * dnxm1
 149             call exact_solution(xi, eta, zeta, temp)
 150             do m = 1, 5
 151                u(m,ii,jj,kk,c) = temp(m)
 152             enddo
 153             ii = ii + 1
 154          enddo
 155          kk = kk + 1
 156       enddo
 157
 158
 159 c---------------------------------------------------------------------
 160 c     north face
 161 c---------------------------------------------------------------------
 162       c = slice(2,ncells)
 163       jj = cell_size(2,c)-1
 164       eta = 1.0d0
 165       kk = 0
 166       do k = cell_low(3,c), cell_high(3,c)
 167          zeta = dble(k) * dnzm1
 168          ii = 0
 169          do i = cell_low(1,c), cell_high(1,c)
 170             xi = dble(i) * dnxm1
 171             call exact_solution(xi, eta, zeta, temp)
 172             do m = 1, 5
 173                u(m,ii,jj,kk,c) = temp(m)
 174             enddo
 175             ii = ii + 1
 176          enddo
 177          kk = kk + 1
 178       enddo
 179
 180 c---------------------------------------------------------------------
 181 c     bottom face
 182 c---------------------------------------------------------------------
 183       c = slice(3,1)
 184       kk = 0
 185       zeta = 0.0d0
 186       jj = 0
 187       do j = cell_low(2,c), cell_high(2,c)
 188          eta = dble(j) * dnym1
 189          ii = 0
 190          do i =cell_low(1,c), cell_high(1,c)
 191             xi = dble(i) *dnxm1
 192             call exact_solution(xi, eta, zeta, temp)
 193             do m = 1, 5
 194                u(m,ii,jj,kk,c) = temp(m)
 195             enddo
 196             ii = ii + 1
 197          enddo
 198          jj = jj + 1
 199       enddo
 200
 201 c---------------------------------------------------------------------
 202 c     top face
 203 c---------------------------------------------------------------------
 204       c = slice(3,ncells)
 205       kk = cell_size(3,c)-1
 206       zeta = 1.0d0
 207       jj = 0
 208       do j = cell_low(2,c), cell_high(2,c)
 209          eta = dble(j) * dnym1
 210          ii = 0
 211          do i =cell_low(1,c), cell_high(1,c)
 212             xi = dble(i) * dnxm1
 213             call exact_solution(xi, eta, zeta, temp)
 214             do m = 1, 5
 215                u(m,ii,jj,kk,c) = temp(m)
 216             enddo
 217             ii = ii + 1
 218          enddo
 219          jj = jj + 1
 220       enddo
 221
 222       return
 223       end
 224
 225
 226 c---------------------------------------------------------------------
 227 c---------------------------------------------------------------------
 228
 229       subroutine lhsinit
 230
 231 c---------------------------------------------------------------------
 232 c---------------------------------------------------------------------
 233
 234       include 'header.h'
 235
 236       integer i, j, k, d, c, m, n
 237
 238 c---------------------------------------------------------------------
 239 c     loop over all cells
 240 c---------------------------------------------------------------------
 241       do c = 1, ncells
 242
 243 c---------------------------------------------------------------------
 244 c     first, initialize the start and end arrays
 245 c---------------------------------------------------------------------
 246          do d = 1, 3
 247             if (cell_coord(d,c) .eq. 1) then
 248                start(d,c) = 1
 249             else
 250                start(d,c) = 0
 251             endif
 252             if (cell_coord(d,c) .eq. ncells) then
 253                end(d,c) = 1
 254             else
 255                end(d,c) = 0
 256             endif
 257          enddo
 258
 259 c---------------------------------------------------------------------
 260 c     zero the whole left hand side for starters
 261 c---------------------------------------------------------------------
 262          do k = 0, cell_size(3,c)-1
 263             do j = 0, cell_size(2,c)-1
 264                do i = 0, cell_size(1,c)-1
 265                   do m = 1,5
 266                      do n = 1, 5
 267                         lhsc(m,n,i,j,k,c) = 0.0d0
 268                      enddo
 269                   enddo
 270                enddo
 271             enddo
 272          enddo
 273
 274       enddo
 275
 276       return
 277       end
 278
 279
 280 c---------------------------------------------------------------------
 281 c---------------------------------------------------------------------
 282
 283       subroutine lhsabinit(lhsa, lhsb, size)
 284       implicit none
 285
 286       integer size
 287       double precision lhsa(5, 5, -1:size), lhsb(5, 5, -1:size)
 288
 289       integer i, m, n
 290
 291 c---------------------------------------------------------------------
 292 c     next, set all diagonal values to 1. This is overkill, but convenient
 293 c---------------------------------------------------------------------
 294       do i = 0, size
 295          do m = 1, 5
 296             do n = 1, 5
 297                lhsa(m,n,i) = 0.0d0
 298                lhsb(m,n,i) = 0.0d0
 299             enddo
 300             lhsb(m,m,i) = 1.0d0
 301          enddo
 302       enddo
 303
 304       return
 305       end
 306
 307
 308