examples/smpi/NAS/BT/rhs.f

   1 c---------------------------------------------------------------------
   2 c---------------------------------------------------------------------
   3
   4       subroutine compute_rhs
   5
   6 c---------------------------------------------------------------------
   7 c---------------------------------------------------------------------
   8
   9       include 'header.h'
  10
  11       integer c, i, j, k, m
  12       double precision rho_inv, uijk, up1, um1, vijk, vp1, vm1,
  13      >     wijk, wp1, wm1
  14
  15
  16 c---------------------------------------------------------------------
  17 c     loop over all cells owned by this node
  18 c---------------------------------------------------------------------
  19       do c = 1, ncells
  20
  21 c---------------------------------------------------------------------
  22 c     compute the reciprocal of density, and the kinetic energy,
  23 c     and the speed of sound.
  24 c---------------------------------------------------------------------
  25          do k = -1, cell_size(3,c)
  26             do j = -1, cell_size(2,c)
  27                do i = -1, cell_size(1,c)
  28                   rho_inv = 1.0d0/u(1,i,j,k,c)
  29                   rho_i(i,j,k,c) = rho_inv
  30                   us(i,j,k,c) = u(2,i,j,k,c) * rho_inv
  31                   vs(i,j,k,c) = u(3,i,j,k,c) * rho_inv
  32                   ws(i,j,k,c) = u(4,i,j,k,c) * rho_inv
  33                   square(i,j,k,c)     = 0.5d0* (
  34      >                 u(2,i,j,k,c)*u(2,i,j,k,c) +
  35      >                 u(3,i,j,k,c)*u(3,i,j,k,c) +
  36      >                 u(4,i,j,k,c)*u(4,i,j,k,c) ) * rho_inv
  37                   qs(i,j,k,c) = square(i,j,k,c) * rho_inv
  38                enddo
  39             enddo
  40          enddo
  41
  42 c---------------------------------------------------------------------
  43 c copy the exact forcing term to the right hand side;  because
  44 c this forcing term is known, we can store it on the whole of every
  45 c cell,  including the boundary
  46 c---------------------------------------------------------------------
  47
  48          do k = 0, cell_size(3,c)-1
  49             do j = 0, cell_size(2,c)-1
  50                do i = 0, cell_size(1,c)-1
  51                   do m = 1, 5
  52                      rhs(m,i,j,k,c) = forcing(m,i,j,k,c)
  53                   enddo
  54                enddo
  55             enddo
  56          enddo
  57
  58
  59 c---------------------------------------------------------------------
  60 c     compute xi-direction fluxes
  61 c---------------------------------------------------------------------
  62          do k = start(3,c), cell_size(3,c)-end(3,c)-1
  63             do j = start(2,c), cell_size(2,c)-end(2,c)-1
  64                do i = start(1,c), cell_size(1,c)-end(1,c)-1
  65                   uijk = us(i,j,k,c)
  66                   up1  = us(i+1,j,k,c)
  67                   um1  = us(i-1,j,k,c)
  68
  69                   rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dx1tx1 *
  70      >                 (u(1,i+1,j,k,c) - 2.0d0*u(1,i,j,k,c) +
  71      >                 u(1,i-1,j,k,c)) -
  72      >                 tx2 * (u(2,i+1,j,k,c) - u(2,i-1,j,k,c))
  73
  74                   rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dx2tx1 *
  75      >                 (u(2,i+1,j,k,c) - 2.0d0*u(2,i,j,k,c) +
  76      >                 u(2,i-1,j,k,c)) +
  77      >                 xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
  78      >                 tx2 * (u(2,i+1,j,k,c)*up1 -
  79      >                 u(2,i-1,j,k,c)*um1 +
  80      >                 (u(5,i+1,j,k,c)- square(i+1,j,k,c)-
  81      >                 u(5,i-1,j,k,c)+ square(i-1,j,k,c))*
  82      >                 c2)
  83
  84                   rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dx3tx1 *
  85      >                 (u(3,i+1,j,k,c) - 2.0d0*u(3,i,j,k,c) +
  86      >                 u(3,i-1,j,k,c)) +
  87      >                 xxcon2 * (vs(i+1,j,k,c) - 2.0d0*vs(i,j,k,c) +
  88      >                 vs(i-1,j,k,c)) -
  89      >                 tx2 * (u(3,i+1,j,k,c)*up1 -
  90      >                 u(3,i-1,j,k,c)*um1)
  91
  92                   rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dx4tx1 *
  93      >                 (u(4,i+1,j,k,c) - 2.0d0*u(4,i,j,k,c) +
  94      >                 u(4,i-1,j,k,c)) +
  95      >                 xxcon2 * (ws(i+1,j,k,c) - 2.0d0*ws(i,j,k,c) +
  96      >                 ws(i-1,j,k,c)) -
  97      >                 tx2 * (u(4,i+1,j,k,c)*up1 -
  98      >                 u(4,i-1,j,k,c)*um1)
  99
 100                   rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dx5tx1 *
 101      >                 (u(5,i+1,j,k,c) - 2.0d0*u(5,i,j,k,c) +
 102      >                 u(5,i-1,j,k,c)) +
 103      >                 xxcon3 * (qs(i+1,j,k,c) - 2.0d0*qs(i,j,k,c) +
 104      >                 qs(i-1,j,k,c)) +
 105      >                 xxcon4 * (up1*up1 -       2.0d0*uijk*uijk +
 106      >                 um1*um1) +
 107      >                 xxcon5 * (u(5,i+1,j,k,c)*rho_i(i+1,j,k,c) -
 108      >                 2.0d0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
 109      >                 u(5,i-1,j,k,c)*rho_i(i-1,j,k,c)) -
 110      >                 tx2 * ( (c1*u(5,i+1,j,k,c) -
 111      >                 c2*square(i+1,j,k,c))*up1 -
 112      >                 (c1*u(5,i-1,j,k,c) -
 113      >                 c2*square(i-1,j,k,c))*um1 )
 114                enddo
 115             enddo
 116          enddo
 117
 118 c---------------------------------------------------------------------
 119 c     add fourth order xi-direction dissipation
 120 c---------------------------------------------------------------------
 121          if (start(1,c) .gt. 0) then
 122             do k = start(3,c), cell_size(3,c)-end(3,c)-1
 123                do j = start(2,c), cell_size(2,c)-end(2,c)-1
 124                   i = 1
 125                   do m = 1, 5
 126                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp *
 127      >                    ( 5.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i+1,j,k,c) +
 128      >                    u(m,i+2,j,k,c))
 129                   enddo
 130
 131                   i = 2
 132                   do m = 1, 5
 133                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 134      >                    (-4.0d0*u(m,i-1,j,k,c) + 6.0d0*u(m,i,j,k,c) -
 135      >                    4.0d0*u(m,i+1,j,k,c) + u(m,i+2,j,k,c))
 136                   enddo
 137                enddo
 138             enddo
 139          endif
 140
 141          do k = start(3,c), cell_size(3,c)-end(3,c)-1
 142             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 143                do i = 3*start(1,c),cell_size(1,c)-3*end(1,c)-1
 144                   do m = 1, 5
 145                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 146      >                    (  u(m,i-2,j,k,c) - 4.0d0*u(m,i-1,j,k,c) +
 147      >                    6.0*u(m,i,j,k,c) - 4.0d0*u(m,i+1,j,k,c) +
 148      >                    u(m,i+2,j,k,c) )
 149                   enddo
 150                enddo
 151             enddo
 152          enddo
 153
 154
 155          if (end(1,c) .gt. 0) then
 156             do k = start(3,c), cell_size(3,c)-end(3,c)-1
 157                do j = start(2,c), cell_size(2,c)-end(2,c)-1
 158                   i = cell_size(1,c)-3
 159                   do m = 1, 5
 160                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 161      >                    ( u(m,i-2,j,k,c) - 4.0d0*u(m,i-1,j,k,c) +
 162      >                    6.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i+1,j,k,c) )
 163                   enddo
 164
 165                   i = cell_size(1,c)-2
 166                   do m = 1, 5
 167                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 168      >                    ( u(m,i-2,j,k,c) - 4.d0*u(m,i-1,j,k,c) +
 169      >                    5.d0*u(m,i,j,k,c) )
 170                   enddo
 171                enddo
 172             enddo
 173          endif
 174
 175 c---------------------------------------------------------------------
 176 c     compute eta-direction fluxes
 177 c---------------------------------------------------------------------
 178          do k = start(3,c), cell_size(3,c)-end(3,c)-1
 179             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 180                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 181                   vijk = vs(i,j,k,c)
 182                   vp1  = vs(i,j+1,k,c)
 183                   vm1  = vs(i,j-1,k,c)
 184                   rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dy1ty1 *
 185      >                 (u(1,i,j+1,k,c) - 2.0d0*u(1,i,j,k,c) +
 186      >                 u(1,i,j-1,k,c)) -
 187      >                 ty2 * (u(3,i,j+1,k,c) - u(3,i,j-1,k,c))
 188                   rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dy2ty1 *
 189      >                 (u(2,i,j+1,k,c) - 2.0d0*u(2,i,j,k,c) +
 190      >                 u(2,i,j-1,k,c)) +
 191      >                 yycon2 * (us(i,j+1,k,c) - 2.0d0*us(i,j,k,c) +
 192      >                 us(i,j-1,k,c)) -
 193      >                 ty2 * (u(2,i,j+1,k,c)*vp1 -
 194      >                 u(2,i,j-1,k,c)*vm1)
 195                   rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dy3ty1 *
 196      >                 (u(3,i,j+1,k,c) - 2.0d0*u(3,i,j,k,c) +
 197      >                 u(3,i,j-1,k,c)) +
 198      >                 yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
 199      >                 ty2 * (u(3,i,j+1,k,c)*vp1 -
 200      >                 u(3,i,j-1,k,c)*vm1 +
 201      >                 (u(5,i,j+1,k,c) - square(i,j+1,k,c) -
 202      >                 u(5,i,j-1,k,c) + square(i,j-1,k,c))
 203      >                 *c2)
 204                   rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dy4ty1 *
 205      >                 (u(4,i,j+1,k,c) - 2.0d0*u(4,i,j,k,c) +
 206      >                 u(4,i,j-1,k,c)) +
 207      >                 yycon2 * (ws(i,j+1,k,c) - 2.0d0*ws(i,j,k,c) +
 208      >                 ws(i,j-1,k,c)) -
 209      >                 ty2 * (u(4,i,j+1,k,c)*vp1 -
 210      >                 u(4,i,j-1,k,c)*vm1)
 211                   rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dy5ty1 *
 212      >                 (u(5,i,j+1,k,c) - 2.0d0*u(5,i,j,k,c) +
 213      >                 u(5,i,j-1,k,c)) +
 214      >                 yycon3 * (qs(i,j+1,k,c) - 2.0d0*qs(i,j,k,c) +
 215      >                 qs(i,j-1,k,c)) +
 216      >                 yycon4 * (vp1*vp1       - 2.0d0*vijk*vijk +
 217      >                 vm1*vm1) +
 218      >                 yycon5 * (u(5,i,j+1,k,c)*rho_i(i,j+1,k,c) -
 219      >                 2.0d0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
 220      >                 u(5,i,j-1,k,c)*rho_i(i,j-1,k,c)) -
 221      >                 ty2 * ((c1*u(5,i,j+1,k,c) -
 222      >                 c2*square(i,j+1,k,c)) * vp1 -
 223      >                 (c1*u(5,i,j-1,k,c) -
 224      >                 c2*square(i,j-1,k,c)) * vm1)
 225                enddo
 226             enddo
 227          enddo
 228
 229 c---------------------------------------------------------------------
 230 c     add fourth order eta-direction dissipation
 231 c---------------------------------------------------------------------
 232          if (start(2,c) .gt. 0) then
 233             do k = start(3,c), cell_size(3,c)-end(3,c)-1
 234                j = 1
 235                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 236                   do m = 1, 5
 237                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp *
 238      >                    ( 5.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j+1,k,c) +
 239      >                    u(m,i,j+2,k,c))
 240                   enddo
 241                enddo
 242
 243                j = 2
 244                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 245                   do m = 1, 5
 246                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 247      >                    (-4.0d0*u(m,i,j-1,k,c) + 6.0d0*u(m,i,j,k,c) -
 248      >                    4.0d0*u(m,i,j+1,k,c) + u(m,i,j+2,k,c))
 249                   enddo
 250                enddo
 251             enddo
 252          endif
 253
 254          do k = start(3,c), cell_size(3,c)-end(3,c)-1
 255             do j = 3*start(2,c), cell_size(2,c)-3*end(2,c)-1
 256                do i = start(1,c),cell_size(1,c)-end(1,c)-1
 257                   do m = 1, 5
 258                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 259      >                    (  u(m,i,j-2,k,c) - 4.0d0*u(m,i,j-1,k,c) +
 260      >                    6.0*u(m,i,j,k,c) - 4.0d0*u(m,i,j+1,k,c) +
 261      >                    u(m,i,j+2,k,c) )
 262                   enddo
 263                enddo
 264             enddo
 265          enddo
 266
 267          if (end(2,c) .gt. 0) then
 268             do k = start(3,c), cell_size(3,c)-end(3,c)-1
 269                j = cell_size(2,c)-3
 270                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 271                   do m = 1, 5
 272                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 273      >                    ( u(m,i,j-2,k,c) - 4.0d0*u(m,i,j-1,k,c) +
 274      >                    6.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j+1,k,c) )
 275                   enddo
 276                enddo
 277
 278                j = cell_size(2,c)-2
 279                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 280                   do m = 1, 5
 281                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 282      >                    ( u(m,i,j-2,k,c) - 4.d0*u(m,i,j-1,k,c) +
 283      >                    5.d0*u(m,i,j,k,c) )
 284                   enddo
 285                enddo
 286             enddo
 287          endif
 288
 289 c---------------------------------------------------------------------
 290 c     compute zeta-direction fluxes
 291 c---------------------------------------------------------------------
 292          do k = start(3,c), cell_size(3,c)-end(3,c)-1
 293             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 294                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 295                   wijk = ws(i,j,k,c)
 296                   wp1  = ws(i,j,k+1,c)
 297                   wm1  = ws(i,j,k-1,c)
 298
 299                   rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dz1tz1 *
 300      >                 (u(1,i,j,k+1,c) - 2.0d0*u(1,i,j,k,c) +
 301      >                 u(1,i,j,k-1,c)) -
 302      >                 tz2 * (u(4,i,j,k+1,c) - u(4,i,j,k-1,c))
 303                   rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dz2tz1 *
 304      >                 (u(2,i,j,k+1,c) - 2.0d0*u(2,i,j,k,c) +
 305      >                 u(2,i,j,k-1,c)) +
 306      >                 zzcon2 * (us(i,j,k+1,c) - 2.0d0*us(i,j,k,c) +
 307      >                 us(i,j,k-1,c)) -
 308      >                 tz2 * (u(2,i,j,k+1,c)*wp1 -
 309      >                 u(2,i,j,k-1,c)*wm1)
 310                   rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dz3tz1 *
 311      >                 (u(3,i,j,k+1,c) - 2.0d0*u(3,i,j,k,c) +
 312      >                 u(3,i,j,k-1,c)) +
 313      >                 zzcon2 * (vs(i,j,k+1,c) - 2.0d0*vs(i,j,k,c) +
 314      >                 vs(i,j,k-1,c)) -
 315      >                 tz2 * (u(3,i,j,k+1,c)*wp1 -
 316      >                 u(3,i,j,k-1,c)*wm1)
 317                   rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dz4tz1 *
 318      >                 (u(4,i,j,k+1,c) - 2.0d0*u(4,i,j,k,c) +
 319      >                 u(4,i,j,k-1,c)) +
 320      >                 zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
 321      >                 tz2 * (u(4,i,j,k+1,c)*wp1 -
 322      >                 u(4,i,j,k-1,c)*wm1 +
 323      >                 (u(5,i,j,k+1,c) - square(i,j,k+1,c) -
 324      >                 u(5,i,j,k-1,c) + square(i,j,k-1,c))
 325      >                 *c2)
 326                   rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dz5tz1 *
 327      >                 (u(5,i,j,k+1,c) - 2.0d0*u(5,i,j,k,c) +
 328      >                 u(5,i,j,k-1,c)) +
 329      >                 zzcon3 * (qs(i,j,k+1,c) - 2.0d0*qs(i,j,k,c) +
 330      >                 qs(i,j,k-1,c)) +
 331      >                 zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk +
 332      >                 wm1*wm1) +
 333      >                 zzcon5 * (u(5,i,j,k+1,c)*rho_i(i,j,k+1,c) -
 334      >                 2.0d0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
 335      >                 u(5,i,j,k-1,c)*rho_i(i,j,k-1,c)) -
 336      >                 tz2 * ( (c1*u(5,i,j,k+1,c) -
 337      >                 c2*square(i,j,k+1,c))*wp1 -
 338      >                 (c1*u(5,i,j,k-1,c) -
 339      >                 c2*square(i,j,k-1,c))*wm1)
 340                enddo
 341             enddo
 342          enddo
 343
 344 c---------------------------------------------------------------------
 345 c     add fourth order zeta-direction dissipation
 346 c---------------------------------------------------------------------
 347          if (start(3,c) .gt. 0) then
 348             k = 1
 349             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 350                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 351                   do m = 1, 5
 352                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp *
 353      >                    ( 5.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j,k+1,c) +
 354      >                    u(m,i,j,k+2,c))
 355                   enddo
 356                enddo
 357             enddo
 358
 359             k = 2
 360             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 361                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 362                   do m = 1, 5
 363                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 364      >                    (-4.0d0*u(m,i,j,k-1,c) + 6.0d0*u(m,i,j,k,c) -
 365      >                    4.0d0*u(m,i,j,k+1,c) + u(m,i,j,k+2,c))
 366                   enddo
 367                enddo
 368             enddo
 369          endif
 370
 371          do k = 3*start(3,c), cell_size(3,c)-3*end(3,c)-1
 372             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 373                do i = start(1,c),cell_size(1,c)-end(1,c)-1
 374                   do m = 1, 5
 375                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 376      >                    (  u(m,i,j,k-2,c) - 4.0d0*u(m,i,j,k-1,c) +
 377      >                    6.0*u(m,i,j,k,c) - 4.0d0*u(m,i,j,k+1,c) +
 378      >                    u(m,i,j,k+2,c) )
 379                   enddo
 380                enddo
 381             enddo
 382          enddo
 383
 384          if (end(3,c) .gt. 0) then
 385             k = cell_size(3,c)-3
 386             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 387                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 388                   do m = 1, 5
 389                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 390      >                    ( u(m,i,j,k-2,c) - 4.0d0*u(m,i,j,k-1,c) +
 391      >                    6.0d0*u(m,i,j,k,c) - 4.0d0*u(m,i,j,k+1,c) )
 392                   enddo
 393                enddo
 394             enddo
 395
 396             k = cell_size(3,c)-2
 397             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 398                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 399                   do m = 1, 5
 400                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
 401      >                    ( u(m,i,j,k-2,c) - 4.d0*u(m,i,j,k-1,c) +
 402      >                    5.d0*u(m,i,j,k,c) )
 403                   enddo
 404                enddo
 405             enddo
 406          endif
 407
 408          do k = start(3,c), cell_size(3,c)-end(3,c)-1
 409             do j = start(2,c), cell_size(2,c)-end(2,c)-1
 410                do i = start(1,c), cell_size(1,c)-end(1,c)-1
 411                   do m = 1, 5
 412                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c) * dt
 413                   enddo
 414                enddo
 415             enddo
 416          enddo
 417
 418       enddo
 419
 420       return
 421       end
 422
 423
 424
 425