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