examples/smpi/NAS/LU/jacu.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5       subroutine jacu(k)
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c   compute the upper triangular part of the jacobian matrix
  12 c---------------------------------------------------------------------
  13
  14
  15       implicit none
  16
  17       include 'applu.incl'
  18
  19 c---------------------------------------------------------------------
  20 c  input parameters
  21 c---------------------------------------------------------------------
  22       integer k
  23
  24 c---------------------------------------------------------------------
  25 c  local variables
  26 c---------------------------------------------------------------------
  27       integer i, j
  28       double precision  r43
  29       double precision  c1345
  30       double precision  c34
  31       double precision  tmp1, tmp2, tmp3
  32
  33
  34
  35       r43 = ( 4.0d+00 / 3.0d+00 )
  36       c1345 = c1 * c3 * c4 * c5
  37       c34 = c3 * c4
  38
  39          do j = jst, jend
  40             do i = ist, iend
  41
  42 c---------------------------------------------------------------------
  43 c   form the block daigonal
  44 c---------------------------------------------------------------------
  45                tmp1 = 1.0d+00 / u(1,i,j,k)
  46                tmp2 = tmp1 * tmp1
  47                tmp3 = tmp1 * tmp2
  48
  49                d(1,1,i,j) =  1.0d+00
  50      >                       + dt * 2.0d+00 * (   tx1 * dx1
  51      >                                          + ty1 * dy1
  52      >                                          + tz1 * dz1 )
  53                d(1,2,i,j) =  0.0d+00
  54                d(1,3,i,j) =  0.0d+00
  55                d(1,4,i,j) =  0.0d+00
  56                d(1,5,i,j) =  0.0d+00
  57
  58                d(2,1,i,j) =  dt * 2.0d+00
  59      >          * (  tx1 * ( - r43 * c34 * tmp2 * u(2,i,j,k) )
  60      >             + ty1 * ( -       c34 * tmp2 * u(2,i,j,k) )
  61      >             + tz1 * ( -       c34 * tmp2 * u(2,i,j,k) ) )
  62                d(2,2,i,j) =  1.0d+00
  63      >          + dt * 2.0d+00
  64      >          * (  tx1 * r43 * c34 * tmp1
  65      >             + ty1 *       c34 * tmp1
  66      >             + tz1 *       c34 * tmp1 )
  67      >          + dt * 2.0d+00 * (   tx1 * dx2
  68      >                             + ty1 * dy2
  69      >                             + tz1 * dz2  )
  70                d(2,3,i,j) = 0.0d+00
  71                d(2,4,i,j) = 0.0d+00
  72                d(2,5,i,j) = 0.0d+00
  73
  74                d(3,1,i,j) = dt * 2.0d+00
  75      >      * (  tx1 * ( -       c34 * tmp2 * u(3,i,j,k) )
  76      >         + ty1 * ( - r43 * c34 * tmp2 * u(3,i,j,k) )
  77      >         + tz1 * ( -       c34 * tmp2 * u(3,i,j,k) ) )
  78                d(3,2,i,j) = 0.0d+00
  79                d(3,3,i,j) = 1.0d+00
  80      >         + dt * 2.0d+00
  81      >              * (  tx1 *       c34 * tmp1
  82      >                 + ty1 * r43 * c34 * tmp1
  83      >                 + tz1 *       c34 * tmp1 )
  84      >         + dt * 2.0d+00 * (  tx1 * dx3
  85      >                           + ty1 * dy3
  86      >                           + tz1 * dz3 )
  87                d(3,4,i,j) = 0.0d+00
  88                d(3,5,i,j) = 0.0d+00
  89
  90                d(4,1,i,j) = dt * 2.0d+00
  91      >      * (  tx1 * ( -       c34 * tmp2 * u(4,i,j,k) )
  92      >         + ty1 * ( -       c34 * tmp2 * u(4,i,j,k) )
  93      >         + tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k) ) )
  94                d(4,2,i,j) = 0.0d+00
  95                d(4,3,i,j) = 0.0d+00
  96                d(4,4,i,j) = 1.0d+00
  97      >         + dt * 2.0d+00
  98      >              * (  tx1 *       c34 * tmp1
  99      >                 + ty1 *       c34 * tmp1
 100      >                 + tz1 * r43 * c34 * tmp1 )
 101      >         + dt * 2.0d+00 * (  tx1 * dx4
 102      >                           + ty1 * dy4
 103      >                           + tz1 * dz4 )
 104                d(4,5,i,j) = 0.0d+00
 105
 106                d(5,1,i,j) = dt * 2.0d+00
 107      > * ( tx1 * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
 108      >             - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
 109      >             - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
 110      >             - ( c1345 ) * tmp2 * u(5,i,j,k) )
 111      >   + ty1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
 112      >             - ( r43*c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
 113      >             - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
 114      >             - ( c1345 ) * tmp2 * u(5,i,j,k) )
 115      >   + tz1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
 116      >             - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
 117      >             - ( r43*c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
 118      >             - ( c1345 ) * tmp2 * u(5,i,j,k) ) )
 119                d(5,2,i,j) = dt * 2.0d+00
 120      > * ( tx1 * ( r43*c34 - c1345 ) * tmp2 * u(2,i,j,k)
 121      >   + ty1 * (     c34 - c1345 ) * tmp2 * u(2,i,j,k)
 122      >   + tz1 * (     c34 - c1345 ) * tmp2 * u(2,i,j,k) )
 123                d(5,3,i,j) = dt * 2.0d+00
 124      > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k)
 125      >   + ty1 * ( r43*c34 -c1345 ) * tmp2 * u(3,i,j,k)
 126      >   + tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k) )
 127                d(5,4,i,j) = dt * 2.0d+00
 128      > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
 129      >   + ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
 130      >   + tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k) )
 131                d(5,5,i,j) = 1.0d+00
 132      >   + dt * 2.0d+00 * ( tx1 * c1345 * tmp1
 133      >                    + ty1 * c1345 * tmp1
 134      >                    + tz1 * c1345 * tmp1 )
 135      >   + dt * 2.0d+00 * (  tx1 * dx5
 136      >                    +  ty1 * dy5
 137      >                    +  tz1 * dz5 )
 138
 139 c---------------------------------------------------------------------
 140 c   form the first block sub-diagonal
 141 c---------------------------------------------------------------------
 142                tmp1 = 1.0d+00 / u(1,i+1,j,k)
 143                tmp2 = tmp1 * tmp1
 144                tmp3 = tmp1 * tmp2
 145
 146                a(1,1,i,j) = - dt * tx1 * dx1
 147                a(1,2,i,j) =   dt * tx2
 148                a(1,3,i,j) =   0.0d+00
 149                a(1,4,i,j) =   0.0d+00
 150                a(1,5,i,j) =   0.0d+00
 151
 152                a(2,1,i,j) =  dt * tx2
 153      >          * ( - ( u(2,i+1,j,k) * tmp1 ) ** 2
 154      >     + c2 * 0.50d+00 * (  u(2,i+1,j,k) * u(2,i+1,j,k)
 155      >                        + u(3,i+1,j,k) * u(3,i+1,j,k)
 156      >                        + u(4,i+1,j,k) * u(4,i+1,j,k) ) * tmp2 )
 157      >          - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i+1,j,k) )
 158                a(2,2,i,j) =  dt * tx2
 159      >          * ( ( 2.0d+00 - c2 ) * ( u(2,i+1,j,k) * tmp1 ) )
 160      >          - dt * tx1 * ( r43 * c34 * tmp1 )
 161      >          - dt * tx1 * dx2
 162                a(2,3,i,j) =  dt * tx2
 163      >              * ( - c2 * ( u(3,i+1,j,k) * tmp1 ) )
 164                a(2,4,i,j) =  dt * tx2
 165      >              * ( - c2 * ( u(4,i+1,j,k) * tmp1 ) )
 166                a(2,5,i,j) =  dt * tx2 * c2
 167
 168                a(3,1,i,j) =  dt * tx2
 169      >              * ( - ( u(2,i+1,j,k) * u(3,i+1,j,k) ) * tmp2 )
 170      >         - dt * tx1 * ( - c34 * tmp2 * u(3,i+1,j,k) )
 171                a(3,2,i,j) =  dt * tx2 * ( u(3,i+1,j,k) * tmp1 )
 172                a(3,3,i,j) =  dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
 173      >          - dt * tx1 * ( c34 * tmp1 )
 174      >          - dt * tx1 * dx3
 175                a(3,4,i,j) = 0.0d+00
 176                a(3,5,i,j) = 0.0d+00
 177
 178                a(4,1,i,j) = dt * tx2
 179      >          * ( - ( u(2,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 )
 180      >          - dt * tx1 * ( - c34 * tmp2 * u(4,i+1,j,k) )
 181                a(4,2,i,j) = dt * tx2 * ( u(4,i+1,j,k) * tmp1 )
 182                a(4,3,i,j) = 0.0d+00
 183                a(4,4,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
 184      >          - dt * tx1 * ( c34 * tmp1 )
 185      >          - dt * tx1 * dx4
 186                a(4,5,i,j) = 0.0d+00
 187
 188                a(5,1,i,j) = dt * tx2
 189      >          * ( ( c2 * (  u(2,i+1,j,k) * u(2,i+1,j,k)
 190      >                      + u(3,i+1,j,k) * u(3,i+1,j,k)
 191      >                      + u(4,i+1,j,k) * u(4,i+1,j,k) ) * tmp2
 192      >              - c1 * ( u(5,i+1,j,k) * tmp1 ) )
 193      >          * ( u(2,i+1,j,k) * tmp1 ) )
 194      >          - dt * tx1
 195      >          * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i+1,j,k)**2 )
 196      >              - (     c34 - c1345 ) * tmp3 * ( u(3,i+1,j,k)**2 )
 197      >              - (     c34 - c1345 ) * tmp3 * ( u(4,i+1,j,k)**2 )
 198      >              - c1345 * tmp2 * u(5,i+1,j,k) )
 199                a(5,2,i,j) = dt * tx2
 200      >          * ( c1 * ( u(5,i+1,j,k) * tmp1 )
 201      >             - 0.50d+00 * c2
 202      >             * ( (  3.0d+00*u(2,i+1,j,k)*u(2,i+1,j,k)
 203      >                  + u(3,i+1,j,k)*u(3,i+1,j,k)
 204      >                  + u(4,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 ) )
 205      >           - dt * tx1
 206      >           * ( r43*c34 - c1345 ) * tmp2 * u(2,i+1,j,k)
 207                a(5,3,i,j) = dt * tx2
 208      >           * ( - c2 * ( u(3,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
 209      >           - dt * tx1
 210      >           * (  c34 - c1345 ) * tmp2 * u(3,i+1,j,k)
 211                a(5,4,i,j) = dt * tx2
 212      >           * ( - c2 * ( u(4,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
 213      >           - dt * tx1
 214      >           * (  c34 - c1345 ) * tmp2 * u(4,i+1,j,k)
 215                a(5,5,i,j) = dt * tx2
 216      >           * ( c1 * ( u(2,i+1,j,k) * tmp1 ) )
 217      >           - dt * tx1 * c1345 * tmp1
 218      >           - dt * tx1 * dx5
 219
 220 c---------------------------------------------------------------------
 221 c   form the second block sub-diagonal
 222 c---------------------------------------------------------------------
 223                tmp1 = 1.0d+00 / u(1,i,j+1,k)
 224                tmp2 = tmp1 * tmp1
 225                tmp3 = tmp1 * tmp2
 226
 227                b(1,1,i,j) = - dt * ty1 * dy1
 228                b(1,2,i,j) =   0.0d+00
 229                b(1,3,i,j) =  dt * ty2
 230                b(1,4,i,j) =   0.0d+00
 231                b(1,5,i,j) =   0.0d+00
 232
 233                b(2,1,i,j) =  dt * ty2
 234      >           * ( - ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
 235      >           - dt * ty1 * ( - c34 * tmp2 * u(2,i,j+1,k) )
 236                b(2,2,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
 237      >          - dt * ty1 * ( c34 * tmp1 )
 238      >          - dt * ty1 * dy2
 239                b(2,3,i,j) =  dt * ty2 * ( u(2,i,j+1,k) * tmp1 )
 240                b(2,4,i,j) = 0.0d+00
 241                b(2,5,i,j) = 0.0d+00
 242
 243                b(3,1,i,j) =  dt * ty2
 244      >           * ( - ( u(3,i,j+1,k) * tmp1 ) ** 2
 245      >      + 0.50d+00 * c2 * ( (  u(2,i,j+1,k) * u(2,i,j+1,k)
 246      >                           + u(3,i,j+1,k) * u(3,i,j+1,k)
 247      >                           + u(4,i,j+1,k) * u(4,i,j+1,k) )
 248      >                          * tmp2 ) )
 249      >       - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j+1,k) )
 250                b(3,2,i,j) =  dt * ty2
 251      >                   * ( - c2 * ( u(2,i,j+1,k) * tmp1 ) )
 252                b(3,3,i,j) =  dt * ty2 * ( ( 2.0d+00 - c2 )
 253      >                   * ( u(3,i,j+1,k) * tmp1 ) )
 254      >       - dt * ty1 * ( r43 * c34 * tmp1 )
 255      >       - dt * ty1 * dy3
 256                b(3,4,i,j) =  dt * ty2
 257      >                   * ( - c2 * ( u(4,i,j+1,k) * tmp1 ) )
 258                b(3,5,i,j) =  dt * ty2 * c2
 259
 260                b(4,1,i,j) =  dt * ty2
 261      >              * ( - ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
 262      >       - dt * ty1 * ( - c34 * tmp2 * u(4,i,j+1,k) )
 263                b(4,2,i,j) = 0.0d+00
 264                b(4,3,i,j) =  dt * ty2 * ( u(4,i,j+1,k) * tmp1 )
 265                b(4,4,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
 266      >                        - dt * ty1 * ( c34 * tmp1 )
 267      >                        - dt * ty1 * dy4
 268                b(4,5,i,j) = 0.0d+00
 269
 270                b(5,1,i,j) =  dt * ty2
 271      >          * ( ( c2 * (  u(2,i,j+1,k) * u(2,i,j+1,k)
 272      >                      + u(3,i,j+1,k) * u(3,i,j+1,k)
 273      >                      + u(4,i,j+1,k) * u(4,i,j+1,k) ) * tmp2
 274      >               - c1 * ( u(5,i,j+1,k) * tmp1 ) )
 275      >          * ( u(3,i,j+1,k) * tmp1 ) )
 276      >          - dt * ty1
 277      >          * ( - (     c34 - c1345 )*tmp3*(u(2,i,j+1,k)**2)
 278      >              - ( r43*c34 - c1345 )*tmp3*(u(3,i,j+1,k)**2)
 279      >              - (     c34 - c1345 )*tmp3*(u(4,i,j+1,k)**2)
 280      >              - c1345*tmp2*u(5,i,j+1,k) )
 281                b(5,2,i,j) =  dt * ty2
 282      >          * ( - c2 * ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
 283      >          - dt * ty1
 284      >          * ( c34 - c1345 ) * tmp2 * u(2,i,j+1,k)
 285                b(5,3,i,j) =  dt * ty2
 286      >          * ( c1 * ( u(5,i,j+1,k) * tmp1 )
 287      >          - 0.50d+00 * c2
 288      >          * ( (  u(2,i,j+1,k)*u(2,i,j+1,k)
 289      >               + 3.0d+00 * u(3,i,j+1,k)*u(3,i,j+1,k)
 290      >               + u(4,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 ) )
 291      >          - dt * ty1
 292      >          * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j+1,k)
 293                b(5,4,i,j) =  dt * ty2
 294      >          * ( - c2 * ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
 295      >          - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j+1,k)
 296                b(5,5,i,j) =  dt * ty2
 297      >          * ( c1 * ( u(3,i,j+1,k) * tmp1 ) )
 298      >          - dt * ty1 * c1345 * tmp1
 299      >          - dt * ty1 * dy5
 300
 301 c---------------------------------------------------------------------
 302 c   form the third block sub-diagonal
 303 c---------------------------------------------------------------------
 304                tmp1 = 1.0d+00 / u(1,i,j,k+1)
 305                tmp2 = tmp1 * tmp1
 306                tmp3 = tmp1 * tmp2
 307
 308                c(1,1,i,j) = - dt * tz1 * dz1
 309                c(1,2,i,j) =   0.0d+00
 310                c(1,3,i,j) =   0.0d+00
 311                c(1,4,i,j) = dt * tz2
 312                c(1,5,i,j) =   0.0d+00
 313
 314                c(2,1,i,j) = dt * tz2
 315      >           * ( - ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
 316      >           - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k+1) )
 317                c(2,2,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
 318      >           - dt * tz1 * c34 * tmp1
 319      >           - dt * tz1 * dz2
 320                c(2,3,i,j) = 0.0d+00
 321                c(2,4,i,j) = dt * tz2 * ( u(2,i,j,k+1) * tmp1 )
 322                c(2,5,i,j) = 0.0d+00
 323
 324                c(3,1,i,j) = dt * tz2
 325      >           * ( - ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
 326      >           - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k+1) )
 327                c(3,2,i,j) = 0.0d+00
 328                c(3,3,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
 329      >           - dt * tz1 * ( c34 * tmp1 )
 330      >           - dt * tz1 * dz3
 331                c(3,4,i,j) = dt * tz2 * ( u(3,i,j,k+1) * tmp1 )
 332                c(3,5,i,j) = 0.0d+00
 333
 334                c(4,1,i,j) = dt * tz2
 335      >        * ( - ( u(4,i,j,k+1) * tmp1 ) ** 2
 336      >            + 0.50d+00 * c2
 337      >            * ( ( u(2,i,j,k+1) * u(2,i,j,k+1)
 338      >                + u(3,i,j,k+1) * u(3,i,j,k+1)
 339      >                + u(4,i,j,k+1) * u(4,i,j,k+1) ) * tmp2 ) )
 340      >        - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k+1) )
 341                c(4,2,i,j) = dt * tz2
 342      >             * ( - c2 * ( u(2,i,j,k+1) * tmp1 ) )
 343                c(4,3,i,j) = dt * tz2
 344      >             * ( - c2 * ( u(3,i,j,k+1) * tmp1 ) )
 345                c(4,4,i,j) = dt * tz2 * ( 2.0d+00 - c2 )
 346      >             * ( u(4,i,j,k+1) * tmp1 )
 347      >             - dt * tz1 * ( r43 * c34 * tmp1 )
 348      >             - dt * tz1 * dz4
 349                c(4,5,i,j) = dt * tz2 * c2
 350
 351                c(5,1,i,j) = dt * tz2
 352      >     * ( ( c2 * (  u(2,i,j,k+1) * u(2,i,j,k+1)
 353      >                 + u(3,i,j,k+1) * u(3,i,j,k+1)
 354      >                 + u(4,i,j,k+1) * u(4,i,j,k+1) ) * tmp2
 355      >       - c1 * ( u(5,i,j,k+1) * tmp1 ) )
 356      >            * ( u(4,i,j,k+1) * tmp1 ) )
 357      >       - dt * tz1
 358      >       * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k+1)**2)
 359      >           - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k+1)**2)
 360      >           - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k+1)**2)
 361      >          - c1345 * tmp2 * u(5,i,j,k+1) )
 362                c(5,2,i,j) = dt * tz2
 363      >       * ( - c2 * ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
 364      >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k+1)
 365                c(5,3,i,j) = dt * tz2
 366      >       * ( - c2 * ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
 367      >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k+1)
 368                c(5,4,i,j) = dt * tz2
 369      >       * ( c1 * ( u(5,i,j,k+1) * tmp1 )
 370      >       - 0.50d+00 * c2
 371      >       * ( (  u(2,i,j,k+1)*u(2,i,j,k+1)
 372      >            + u(3,i,j,k+1)*u(3,i,j,k+1)
 373      >            + 3.0d+00*u(4,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 ) )
 374      >       - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k+1)
 375                c(5,5,i,j) = dt * tz2
 376      >       * ( c1 * ( u(4,i,j,k+1) * tmp1 ) )
 377      >       - dt * tz1 * c1345 * tmp1
 378      >       - dt * tz1 * dz5
 379
 380             end do
 381          end do
 382
 383       return
 384       end