examples/smpi/NAS/SP/exact_rhs.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5        subroutine exact_rhs
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c compute the right hand side based on exact solution
  12 c---------------------------------------------------------------------
  13
  14        include 'header.h'
  15
  16        double precision dtemp(5), xi, eta, zeta, dtpp
  17        integer          c, m, i, j, k, ip1, im1, jp1,
  18      >                  jm1, km1, kp1
  19
  20 c---------------------------------------------------------------------
  21 c loop over all cells owned by this node
  22 c---------------------------------------------------------------------
  23        do   c = 1, ncells
  24
  25 c---------------------------------------------------------------------
  26 c         initialize
  27 c---------------------------------------------------------------------
  28           do   m = 1, 5
  29              do   k= 0, cell_size(3,c)-1
  30                 do   j = 0, cell_size(2,c)-1
  31                    do   i = 0, cell_size(1,c)-1
  32                       forcing(i,j,k,m,c) = 0.0d0
  33                    end do
  34                 end do
  35              end do
  36           end do
  37
  38 c---------------------------------------------------------------------
  39 c xi-direction flux differences
  40 c---------------------------------------------------------------------
  41           do   k = start(3,c), cell_size(3,c)-end(3,c)-1
  42              zeta = dble(k+cell_low(3,c)) * dnzm1
  43              do   j = start(2,c), cell_size(2,c)-end(2,c)-1
  44                 eta = dble(j+cell_low(2,c)) * dnym1
  45
  46                 do  i=-2*(1-start(1,c)), cell_size(1,c)+1-2*end(1,c)
  47                    xi = dble(i+cell_low(1,c)) * dnxm1
  48
  49                    call exact_solution(xi, eta, zeta, dtemp)
  50                    do  m = 1, 5
  51                       ue(i,m) = dtemp(m)
  52                    end do
  53
  54                    dtpp = 1.0d0 / dtemp(1)
  55
  56                    do  m = 2, 5
  57                       buf(i,m) = dtpp * dtemp(m)
  58                    end do
  59
  60                    cuf(i)   = buf(i,2) * buf(i,2)
  61                    buf(i,1) = cuf(i) + buf(i,3) * buf(i,3) +
  62      >                        buf(i,4) * buf(i,4)
  63                    q(i) = 0.5d0*(buf(i,2)*ue(i,2) + buf(i,3)*ue(i,3) +
  64      >                           buf(i,4)*ue(i,4))
  65
  66                 end do
  67
  68                 do  i = start(1,c), cell_size(1,c)-end(1,c)-1
  69                    im1 = i-1
  70                    ip1 = i+1
  71
  72                    forcing(i,j,k,1,c) = forcing(i,j,k,1,c) -
  73      >                 tx2*( ue(ip1,2)-ue(im1,2) )+
  74      >                 dx1tx1*(ue(ip1,1)-2.0d0*ue(i,1)+ue(im1,1))
  75
  76                    forcing(i,j,k,2,c) = forcing(i,j,k,2,c) - tx2 * (
  77      >                (ue(ip1,2)*buf(ip1,2)+c2*(ue(ip1,5)-q(ip1)))-
  78      >                (ue(im1,2)*buf(im1,2)+c2*(ue(im1,5)-q(im1))))+
  79      >                 xxcon1*(buf(ip1,2)-2.0d0*buf(i,2)+buf(im1,2))+
  80      >                 dx2tx1*( ue(ip1,2)-2.0d0* ue(i,2)+ue(im1,2))
  81
  82                    forcing(i,j,k,3,c) = forcing(i,j,k,3,c) - tx2 * (
  83      >                 ue(ip1,3)*buf(ip1,2)-ue(im1,3)*buf(im1,2))+
  84      >                 xxcon2*(buf(ip1,3)-2.0d0*buf(i,3)+buf(im1,3))+
  85      >                 dx3tx1*( ue(ip1,3)-2.0d0*ue(i,3) +ue(im1,3))
  86
  87                    forcing(i,j,k,4,c) = forcing(i,j,k,4,c) - tx2*(
  88      >                 ue(ip1,4)*buf(ip1,2)-ue(im1,4)*buf(im1,2))+
  89      >                 xxcon2*(buf(ip1,4)-2.0d0*buf(i,4)+buf(im1,4))+
  90      >                 dx4tx1*( ue(ip1,4)-2.0d0* ue(i,4)+ ue(im1,4))
  91
  92                    forcing(i,j,k,5,c) = forcing(i,j,k,5,c) - tx2*(
  93      >                 buf(ip1,2)*(c1*ue(ip1,5)-c2*q(ip1))-
  94      >                 buf(im1,2)*(c1*ue(im1,5)-c2*q(im1)))+
  95      >                 0.5d0*xxcon3*(buf(ip1,1)-2.0d0*buf(i,1)+
  96      >                               buf(im1,1))+
  97      >                 xxcon4*(cuf(ip1)-2.0d0*cuf(i)+cuf(im1))+
  98      >                 xxcon5*(buf(ip1,5)-2.0d0*buf(i,5)+buf(im1,5))+
  99      >                 dx5tx1*( ue(ip1,5)-2.0d0* ue(i,5)+ ue(im1,5))
 100                 end do
 101
 102 c---------------------------------------------------------------------
 103 c Fourth-order dissipation
 104 c---------------------------------------------------------------------
 105                 if (start(1,c) .gt. 0) then
 106                    do   m = 1, 5
 107                       i = 1
 108                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 109      >                    (5.0d0*ue(i,m) - 4.0d0*ue(i+1,m) +ue(i+2,m))
 110                       i = 2
 111                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 112      >                   (-4.0d0*ue(i-1,m) + 6.0d0*ue(i,m) -
 113      >                     4.0d0*ue(i+1,m) +       ue(i+2,m))
 114                    end do
 115                 endif
 116
 117                 do   m = 1, 5
 118                    do  i = start(1,c)*3, cell_size(1,c)-3*end(1,c)-1
 119                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp*
 120      >                   (ue(i-2,m) - 4.0d0*ue(i-1,m) +
 121      >                    6.0d0*ue(i,m) - 4.0d0*ue(i+1,m) + ue(i+2,m))
 122                    end do
 123                 end do
 124
 125                 if (end(1,c) .gt. 0) then
 126                    do   m = 1, 5
 127                       i = cell_size(1,c)-3
 128                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 129      >                   (ue(i-2,m) - 4.0d0*ue(i-1,m) +
 130      >                    6.0d0*ue(i,m) - 4.0d0*ue(i+1,m))
 131                       i = cell_size(1,c)-2
 132                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 133      >                   (ue(i-2,m) - 4.0d0*ue(i-1,m) + 5.0d0*ue(i,m))
 134                    end do
 135                 endif
 136
 137              end do
 138           end do
 139 c---------------------------------------------------------------------
 140 c  eta-direction flux differences
 141 c---------------------------------------------------------------------
 142           do   k = start(3,c), cell_size(3,c)-end(3,c)-1
 143              zeta = dble(k+cell_low(3,c)) * dnzm1
 144              do   i=start(1,c), cell_size(1,c)-end(1,c)-1
 145                 xi = dble(i+cell_low(1,c)) * dnxm1
 146
 147                 do  j=-2*(1-start(2,c)), cell_size(2,c)+1-2*end(2,c)
 148                    eta = dble(j+cell_low(2,c)) * dnym1
 149
 150                    call exact_solution(xi, eta, zeta, dtemp)
 151                    do   m = 1, 5
 152                       ue(j,m) = dtemp(m)
 153                    end do
 154                    dtpp = 1.0d0/dtemp(1)
 155
 156                    do  m = 2, 5
 157                       buf(j,m) = dtpp * dtemp(m)
 158                    end do
 159
 160                    cuf(j)   = buf(j,3) * buf(j,3)
 161                    buf(j,1) = cuf(j) + buf(j,2) * buf(j,2) +
 162      >                        buf(j,4) * buf(j,4)
 163                    q(j) = 0.5d0*(buf(j,2)*ue(j,2) + buf(j,3)*ue(j,3) +
 164      >                           buf(j,4)*ue(j,4))
 165                 end do
 166
 167                 do  j = start(2,c), cell_size(2,c)-end(2,c)-1
 168                    jm1 = j-1
 169                    jp1 = j+1
 170
 171                    forcing(i,j,k,1,c) = forcing(i,j,k,1,c) -
 172      >                ty2*( ue(jp1,3)-ue(jm1,3) )+
 173      >                dy1ty1*(ue(jp1,1)-2.0d0*ue(j,1)+ue(jm1,1))
 174
 175                    forcing(i,j,k,2,c) = forcing(i,j,k,2,c) - ty2*(
 176      >                ue(jp1,2)*buf(jp1,3)-ue(jm1,2)*buf(jm1,3))+
 177      >                yycon2*(buf(jp1,2)-2.0d0*buf(j,2)+buf(jm1,2))+
 178      >                dy2ty1*( ue(jp1,2)-2.0* ue(j,2)+ ue(jm1,2))
 179
 180                    forcing(i,j,k,3,c) = forcing(i,j,k,3,c) - ty2*(
 181      >                (ue(jp1,3)*buf(jp1,3)+c2*(ue(jp1,5)-q(jp1)))-
 182      >                (ue(jm1,3)*buf(jm1,3)+c2*(ue(jm1,5)-q(jm1))))+
 183      >                yycon1*(buf(jp1,3)-2.0d0*buf(j,3)+buf(jm1,3))+
 184      >                dy3ty1*( ue(jp1,3)-2.0d0*ue(j,3) +ue(jm1,3))
 185
 186                    forcing(i,j,k,4,c) = forcing(i,j,k,4,c) - ty2*(
 187      >                ue(jp1,4)*buf(jp1,3)-ue(jm1,4)*buf(jm1,3))+
 188      >                yycon2*(buf(jp1,4)-2.0d0*buf(j,4)+buf(jm1,4))+
 189      >                dy4ty1*( ue(jp1,4)-2.0d0*ue(j,4)+ ue(jm1,4))
 190
 191                    forcing(i,j,k,5,c) = forcing(i,j,k,5,c) - ty2*(
 192      >                buf(jp1,3)*(c1*ue(jp1,5)-c2*q(jp1))-
 193      >                buf(jm1,3)*(c1*ue(jm1,5)-c2*q(jm1)))+
 194      >                0.5d0*yycon3*(buf(jp1,1)-2.0d0*buf(j,1)+
 195      >                              buf(jm1,1))+
 196      >                yycon4*(cuf(jp1)-2.0d0*cuf(j)+cuf(jm1))+
 197      >                yycon5*(buf(jp1,5)-2.0d0*buf(j,5)+buf(jm1,5))+
 198      >                dy5ty1*(ue(jp1,5)-2.0d0*ue(j,5)+ue(jm1,5))
 199                 end do
 200
 201 c---------------------------------------------------------------------
 202 c Fourth-order dissipation
 203 c---------------------------------------------------------------------
 204                 if (start(2,c) .gt. 0) then
 205                    do   m = 1, 5
 206                       j = 1
 207                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 208      >                    (5.0d0*ue(j,m) - 4.0d0*ue(j+1,m) +ue(j+2,m))
 209                       j = 2
 210                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 211      >                   (-4.0d0*ue(j-1,m) + 6.0d0*ue(j,m) -
 212      >                     4.0d0*ue(j+1,m) +       ue(j+2,m))
 213                    end do
 214                 endif
 215
 216                 do   m = 1, 5
 217                    do  j = start(2,c)*3, cell_size(2,c)-3*end(2,c)-1
 218                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp*
 219      >                   (ue(j-2,m) - 4.0d0*ue(j-1,m) +
 220      >                    6.0d0*ue(j,m) - 4.0d0*ue(j+1,m) + ue(j+2,m))
 221                    end do
 222                 end do
 223                 if (end(2,c) .gt. 0) then
 224                    do   m = 1, 5
 225                       j = cell_size(2,c)-3
 226                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 227      >                   (ue(j-2,m) - 4.0d0*ue(j-1,m) +
 228      >                    6.0d0*ue(j,m) - 4.0d0*ue(j+1,m))
 229                       j = cell_size(2,c)-2
 230                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 231      >                   (ue(j-2,m) - 4.0d0*ue(j-1,m) + 5.0d0*ue(j,m))
 232
 233                    end do
 234                 endif
 235
 236              end do
 237           end do
 238
 239 c---------------------------------------------------------------------
 240 c zeta-direction flux differences
 241 c---------------------------------------------------------------------
 242           do  j=start(2,c), cell_size(2,c)-end(2,c)-1
 243              eta = dble(j+cell_low(2,c)) * dnym1
 244              do   i = start(1,c), cell_size(1,c)-end(1,c)-1
 245                 xi = dble(i+cell_low(1,c)) * dnxm1
 246
 247                 do k=-2*(1-start(3,c)), cell_size(3,c)+1-2*end(3,c)
 248                    zeta = dble(k+cell_low(3,c)) * dnzm1
 249
 250                    call exact_solution(xi, eta, zeta, dtemp)
 251                    do   m = 1, 5
 252                       ue(k,m) = dtemp(m)
 253                    end do
 254
 255                    dtpp = 1.0d0/dtemp(1)
 256
 257                    do   m = 2, 5
 258                       buf(k,m) = dtpp * dtemp(m)
 259                    end do
 260
 261                    cuf(k)   = buf(k,4) * buf(k,4)
 262                    buf(k,1) = cuf(k) + buf(k,2) * buf(k,2) +
 263      >                        buf(k,3) * buf(k,3)
 264                    q(k) = 0.5d0*(buf(k,2)*ue(k,2) + buf(k,3)*ue(k,3) +
 265      >                           buf(k,4)*ue(k,4))
 266                 end do
 267
 268                 do    k=start(3,c), cell_size(3,c)-end(3,c)-1
 269                    km1 = k-1
 270                    kp1 = k+1
 271
 272                    forcing(i,j,k,1,c) = forcing(i,j,k,1,c) -
 273      >                 tz2*( ue(kp1,4)-ue(km1,4) )+
 274      >                 dz1tz1*(ue(kp1,1)-2.0d0*ue(k,1)+ue(km1,1))
 275
 276                    forcing(i,j,k,2,c) = forcing(i,j,k,2,c) - tz2 * (
 277      >                 ue(kp1,2)*buf(kp1,4)-ue(km1,2)*buf(km1,4))+
 278      >                 zzcon2*(buf(kp1,2)-2.0d0*buf(k,2)+buf(km1,2))+
 279      >                 dz2tz1*( ue(kp1,2)-2.0d0* ue(k,2)+ ue(km1,2))
 280
 281                    forcing(i,j,k,3,c) = forcing(i,j,k,3,c) - tz2 * (
 282      >                 ue(kp1,3)*buf(kp1,4)-ue(km1,3)*buf(km1,4))+
 283      >                 zzcon2*(buf(kp1,3)-2.0d0*buf(k,3)+buf(km1,3))+
 284      >                 dz3tz1*(ue(kp1,3)-2.0d0*ue(k,3)+ue(km1,3))
 285
 286                    forcing(i,j,k,4,c) = forcing(i,j,k,4,c) - tz2 * (
 287      >                (ue(kp1,4)*buf(kp1,4)+c2*(ue(kp1,5)-q(kp1)))-
 288      >                (ue(km1,4)*buf(km1,4)+c2*(ue(km1,5)-q(km1))))+
 289      >                zzcon1*(buf(kp1,4)-2.0d0*buf(k,4)+buf(km1,4))+
 290      >                dz4tz1*( ue(kp1,4)-2.0d0*ue(k,4) +ue(km1,4))
 291
 292                    forcing(i,j,k,5,c) = forcing(i,j,k,5,c) - tz2 * (
 293      >                 buf(kp1,4)*(c1*ue(kp1,5)-c2*q(kp1))-
 294      >                 buf(km1,4)*(c1*ue(km1,5)-c2*q(km1)))+
 295      >                 0.5d0*zzcon3*(buf(kp1,1)-2.0d0*buf(k,1)
 296      >                              +buf(km1,1))+
 297      >                 zzcon4*(cuf(kp1)-2.0d0*cuf(k)+cuf(km1))+
 298      >                 zzcon5*(buf(kp1,5)-2.0d0*buf(k,5)+buf(km1,5))+
 299      >                 dz5tz1*( ue(kp1,5)-2.0d0*ue(k,5)+ ue(km1,5))
 300                 end do
 301
 302 c---------------------------------------------------------------------
 303 c Fourth-order dissipation
 304 c---------------------------------------------------------------------
 305                 if (start(3,c) .gt. 0) then
 306                    do   m = 1, 5
 307                       k = 1
 308                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 309      >                    (5.0d0*ue(k,m) - 4.0d0*ue(k+1,m) +ue(k+2,m))
 310                       k = 2
 311                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 312      >                   (-4.0d0*ue(k-1,m) + 6.0d0*ue(k,m) -
 313      >                     4.0d0*ue(k+1,m) +       ue(k+2,m))
 314                    end do
 315                 endif
 316
 317                 do   m = 1, 5
 318                    do  k = start(3,c)*3, cell_size(3,c)-3*end(3,c)-1
 319                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp*
 320      >                   (ue(k-2,m) - 4.0d0*ue(k-1,m) +
 321      >                    6.0d0*ue(k,m) - 4.0d0*ue(k+1,m) + ue(k+2,m))
 322                    end do
 323                 end do
 324
 325                 if (end(3,c) .gt. 0) then
 326                    do    m = 1, 5
 327                       k = cell_size(3,c)-3
 328                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 329      >                   (ue(k-2,m) - 4.0d0*ue(k-1,m) +
 330      >                    6.0d0*ue(k,m) - 4.0d0*ue(k+1,m))
 331                       k = cell_size(3,c)-2
 332                       forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
 333      >                   (ue(k-2,m) - 4.0d0*ue(k-1,m) + 5.0d0*ue(k,m))
 334                    end do
 335                 endif
 336
 337              end do
 338           end do
 339 c---------------------------------------------------------------------
 340 c now change the sign of the forcing function,
 341 c---------------------------------------------------------------------
 342           do   m = 1, 5
 343              do   k = start(3,c), cell_size(3,c)-end(3,c)-1
 344                 do   j = start(2,c), cell_size(2,c)-end(2,c)-1
 345                    do   i = start(1,c), cell_size(1,c)-end(1,c)-1
 346                       forcing(i,j,k,m,c) = -1.d0 * forcing(i,j,k,m,c)
 347                    end do
 348                 end do
 349              end do
 350           end do
 351
 352 c---------------------------------------------------------------------
 353 c      cell loop
 354 c---------------------------------------------------------------------
 355        end do
 356
 357        return
 358        end
 359
 360
 361
 362
 363