examples/smpi/NAS/SP/initialize.f

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