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