1!--------------------------------------------------------------- 2! The following subroutines are from node6t.f in the original 3! code. 4! Last Modified - D. K. Kaushik (1/24/97) 5!--------------------------------------------------------------- 6! 7!=================================== FILLA =================================== 8! 9! Fill the nonzero term of the A matrix 10! - viscous terms are not included right now 11! 12! New version put in by Eric for David Keyes' stuff 13! Includes complete "by hand" differentiation 14! Modified - D. K. Kaushik (1/9/98) 15! cfl is passed as a parameter 16! 17!============================================================================= 18 subroutine FILLA(nnodes,nedge,evec, & 19 & nsface,isface,fxn,fyn,fzn,sxn,syn,szn, & 20 & nsnode,nvnode,nfnode,isnode,ivnode,ifnode, & 21 & qvec,A,cdt, & 22 & rl,vol,xn,yn,zn,cfl,irank,nvertices) 23! 24#include <petsc/finclude/petscsys.h> 25#include <petsc/finclude/petscvec.h> 26#include <petsc/finclude/petscmat.h> 27#include <petsc/finclude/petscpc.h> 28#include <petsc/finclude/petscsnes.h> 29! 30 common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,ipv, 31 & nstage,ncyct,iramp,nitfo,ncyc 32 common/info/title(20),xmach,alpha,yaw,Re,dt,tot,res0,resc, 33 1 ntt,mseq,ivisc,irest,icyc,ihane,ntturb 34 common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1 35 common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0 36 common/gmcom/gtol,icycle,nsrch,ilu0,ifcn 37! 38 real xn(nedge),yn(nedge),zn(nedge),rl(nedge),vol(1) 39 real sxn(1),syn(1),szn(1) 40 real fxn(1),fyn(1),fzn(1), pio, cio, pbo 41 real cdt(1) 42 real val(50), val1(10,5) 43! real A(nnz,4,4) 44 Mat A 45 MatStructure flag 46 47 48 integer isface(1) 49 integer isnode(1),ivnode(1),ifnode(1) 50! integer ia(1),ja(1),iau(1),fhelp(nedge,2) 51 integer irow(25), icol(25) 52 integer i, j, k, in 53 PetscLogDouble flops 54! 55#if defined(INTERLACING) 56 real qvec(5,nvertices) 57 integer evec(2,nedge) 58#define qnode(i,j) qvec(i,j) 59#define dfp(a,b) val1(b,a) 60#define dfm(a,b) val1(b+5,a) 61#define eptr(j,i) evec(i,j) 62#else 63 real qvec(nvertices,5) 64 integer evec(nedge,2) 65#define qnode(i,j) qvec(j,i) 66#define dfp(a,b) val1(b,a) 67#define dfm(a,b) val1(b+5,a) 68#define eptr(i,j) evec(i,j) 69#endif 70! 71! Loop over the edges and fill the matrix 72! First just zero it out 73! 74 flops = 0.0 75 call MatZeroEntries(A,ierr) 76! do 1000 j = 1,4 77! do 1000 k = 1,4 78! do 1000 i = 1,nnz 79! A(i,j,k) = 0.0 80!1000 continue 81! print *, "I am in FILLA" 82! 83! Take care of the time step term on the diagonal 84! 85! cfl = abs(dt) 86! if (iramp.lt.0) then 87! cfl = cfl1*res0/resc 88! if (cfl.gt.cfl2) cfl=cfl2 89! if (ntt.eq.1.and.irest.eq.0) cfl=10 90! else 91! if (dt.lt.0.0) cfl = cfl1 92! if (dt.lt.0.0.and.iramp.gt.0.and.ntt.gt.iramp) cfl = cfl2 93! if (dt.lt.0.0.and.iramp.gt.0.and.ntt.le.iramp) then 94! cfl = cfl1 + (cfl2 - cfl1)*(ntt - 1)/iramp 95! end if 96! end if 97! 98! write(6,555)res0,resc,cfl,cfl1,cfl2 99! 100#if defined(INTERLACING) 101 do i = 1,nnodes 102 temp = vol(i)/(cfl*cdt(i)) 103 do j = 1,5 104 ir = j - 1 + 5*(i-1) 105 call MatSetValuesLocal(A,1,ir,1,ir,temp,ADD_VALUES,ierr) 106 enddo 107 enddo 108 flops = flops + 2*nnodes 109#else 110 do j = 1,5 111 do i = 1,nnodes 112 temp = vol(i)/(cfl*cdt(i)) 113 ir = i - 1 + nnodes*(j-1) 114 call MatSetValues(A,1,ir,1,ir,temp,ADD_VALUES,ierr) 115 enddo 116 enddo 117 flops = flops + 5*2*nnodes 118#endif 119 120! call MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY,ierr) 121! call MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY,ierr) 122! print *, "Finished doing time stepping part to diagonal term" 123! 124! Now fill A from interior contributions 125! We will fix the boundaries later 126! 127! First convert to primitive variables 128! 129 call ETOP(nvertices,qvec) 130! 131 vconst = 1.0 132 if (ivisc.eq.0) vconst = 0.0 133 xmr = xmach/Re 134 c512 = 5./12.*xmr*vconst 135 flops = flops + 4.0 136! 137 do 1030 n = 1, nedge 138 node1 = eptr(n,1) 139 node2 = eptr(n,2) 140 if ((node1 .le. nnodes).or.(node2 .le. nnodes)) then 141! 142! Get unit normals and area 143! 144 xnorm = xn(n) 145 ynorm = yn(n) 146 znorm = zn(n) 147 area = rl(n) 148! 149! Get variables on "left" side of face 150! 151 rhol = qnode(1,node1) 152 rholrL = 1. 153 rholuL = 0. 154 rholvL = 0. 155 rholwL = 0. 156 rholpL = 0. 157! 158 rholrR = 0. 159 rholuR = 0. 160 rholvR = 0. 161 rholwR = 0. 162 rholpR = 0. 163 ul = qnode(2,node1) 164 ulrL = 0.0 165 uluL = 1.0 166 ulvL = 0.0 167 ulwL = 0.0 168 ulpL = 0.0 169! 170 ulrR = 0.0 171 uluR = 0.0 172 ulvR = 0.0 173 ulwR = 0.0 174 ulpR = 0.0 175 vl = qnode(3,node1) 176 vlrL = 0.0 177 vluL = 0.0 178 vlvL = 1.0 179 vlwL = 0.0 180 vlpL = 0.0 181! 182 vlrR = 0.0 183 vluR = 0.0 184 vlvR = 0.0 185 vlwR = 0.0 186 vlpR = 0.0 187 wl = qnode(4,node1) 188 wlrL = 0.0 189 wluL = 0.0 190 wlvL = 0.0 191 wlwL = 1.0 192 wlpL = 0.0 193! 194 wlrR = 0.0 195 wluR = 0.0 196 wlvR = 0.0 197 wlwR = 0.0 198 wlpR = 0.0 199! 200 q2l = ul*ul + vl*vl + wl*wl 201 q2lrL = 2.*ul*ulrL + 2.*vl*vlrL + 2.*wl*wlrL 202 q2luL = 2.*ul*uluL + 2.*vl*vluL + 2.*wl*wluL 203 q2lvL = 2.*ul*ulvL + 2.*vl*vlvL + 2.*wl*wlvL 204 q2lwL = 2.*ul*ulwL + 2.*vl*vlwL + 2.*wl*wlwL 205 q2lpL = 2.*ul*ulpL + 2.*vl*vlpL + 2.*wl*wlpL 206! 207 q2lrR = 2.*ul*ulrR + 2.*vl*vlrR + 2.*wl*wlrR 208 q2luR = 2.*ul*uluR + 2.*vl*vluR + 2.*wl*wluR 209 q2lvR = 2.*ul*ulvR + 2.*vl*vlvR + 2.*wl*wlvR 210 q2lwR = 2.*ul*ulwR + 2.*vl*vlwR + 2.*wl*wlwR 211 q2lpR = 2.*ul*ulpR + 2.*vl*vlpR + 2.*wl*wlpR 212! 213 pressl = qnode(5,node1) 214 presslrL = 0. 215 pressluL = 0. 216 presslvL = 0. 217 presslwL = 0. 218 presslpL = 1. 219! 220 presslrR = 0. 221 pressluR = 0. 222 presslvR = 0. 223 presslwR = 0. 224 presslpR = 0. 225! 226 enrgyl = pressl/gm1 + .5*rhol*q2l 227! 228 enrgylrL = presslrL/gm1 + .5*(rhol*q2lrL + q2l*rholrL) 229 enrgyluL = pressluL/gm1 + .5*(rhol*q2luL + q2l*rholuL) 230 enrgylvL = presslvL/gm1 + .5*(rhol*q2lvL + q2l*rholvL) 231 enrgylwL = presslwL/gm1 + .5*(rhol*q2lwL + q2l*rholwL) 232 enrgylpL = presslpL/gm1 + .5*(rhol*q2lpL + q2l*rholpL) 233! 234 enrgylrR = presslrR/gm1 + .5*(rhol*q2lrR + q2l*rholrR) 235 enrgyluR = pressluR/gm1 + .5*(rhol*q2luR + q2l*rholuR) 236 enrgylvR = presslvR/gm1 + .5*(rhol*q2lvR + q2l*rholvR) 237 enrgylwR = presslwR/gm1 + .5*(rhol*q2lwR + q2l*rholwR) 238 enrgylpR = presslpR/gm1 + .5*(rhol*q2lpR + q2l*rholpR) 239! 240 Hl = (enrgyl + pressl)/rhol 241! 242 HlrL = (rhol*(enrgylrL+presslrL) - (enrgyl+pressl)*rholrL) 243 & / rhol / rhol 244 HluL = (rhol*(enrgyluL+pressluL) - (enrgyl+pressl)*rholuL) 245 & / rhol / rhol 246 HlvL = (rhol*(enrgylvL+presslvL) - (enrgyl+pressl)*rholvL) 247 & / rhol / rhol 248 HlwL = (rhol*(enrgylwL+presslwL) - (enrgyl+pressl)*rholwL) 249 & / rhol / rhol 250 HlpL = (rhol*(enrgylpL+presslpL) - (enrgyl+pressl)*rholpL) 251 & / rhol / rhol 252! 253 HlrR = (rhol*(enrgylrR+presslrR) - (enrgyl+pressl)*rholrR) 254 & / rhol / rhol 255 HluR = (rhol*(enrgyluR+pressluR) - (enrgyl+pressl)*rholuR) 256 & / rhol / rhol 257 HlvR = (rhol*(enrgylvR+presslvR) - (enrgyl+pressl)*rholvR) 258 & / rhol / rhol 259 HlwR = (rhol*(enrgylwR+presslwR) - (enrgyl+pressl)*rholwR) 260 & / rhol / rhol 261 HlpR = (rhol*(enrgylpR+presslpR) - (enrgyl+pressl)*rholpR) 262 & / rhol / rhol 263! 264! 265 cl = sqrt(gamma*pressl/rhol) 266! 267 clrL = 0.5 / sqrt(gamma*pressl/rhol) * 268 & gamma * (rhol*presslrL - pressl*rholrL) 269 & / rhol / rhol 270 cluL = 0.5 / sqrt(gamma*pressl/rhol) * 271 & gamma * (rhol*pressluL - pressl*rholuL) 272 & / rhol / rhol 273 clvL = 0.5 / sqrt(gamma*pressl/rhol) * 274 & gamma * (rhol*presslvL - pressl*rholvL) 275 & / rhol / rhol 276 clwL = 0.5 / sqrt(gamma*pressl/rhol) * 277 & gamma * (rhol*presslwL - pressl*rholwL) 278 & / rhol / rhol 279 clpL = 0.5 / sqrt(gamma*pressl/rhol) * 280 & gamma * (rhol*presslpL - pressl*rholpL) 281 & / rhol / rhol 282! 283 clrR = 0.5 / sqrt(gamma*pressl/rhol) * 284 & gamma * (rhol*presslrR - pressl*rholrR) 285 & / rhol / rhol 286 cluR = 0.5 / sqrt(gamma*pressl/rhol) * 287 & gamma * (rhol*pressluR - pressl*rholuR) 288 & / rhol / rhol 289 clvR = 0.5 / sqrt(gamma*pressl/rhol) * 290 & gamma * (rhol*presslvR - pressl*rholvR) 291 & / rhol / rhol 292 clwR = 0.5 / sqrt(gamma*pressl/rhol) * 293 & gamma * (rhol*presslwR - pressl*rholwR) 294 & / rhol / rhol 295 clpR = 0.5 / sqrt(gamma*pressl/rhol) * 296 & gamma * (rhol*presslpR - pressl*rholpR) 297 & / rhol / rhol 298! 299! 300 ubarl = xnorm*ul + ynorm*vl + znorm*wl 301! 302 ubarlrL = xnorm*ulrL + ynorm*vlrL + znorm*wlrL 303 ubarluL = xnorm*uluL + ynorm*vluL + znorm*wluL 304 ubarlvL = xnorm*ulvL + ynorm*vlvL + znorm*wlvL 305 ubarlwL = xnorm*ulwL + ynorm*vlwL + znorm*wlwL 306 ubarlpL = xnorm*ulpL + ynorm*vlpL + znorm*wlpL 307! 308 ubarlrR = xnorm*ulrR + ynorm*vlrR + znorm*wlrR 309 ubarluR = xnorm*uluR + ynorm*vluR + znorm*wluR 310 ubarlvR = xnorm*ulvR + ynorm*vlvR + znorm*wlvR 311 ubarlwR = xnorm*ulwR + ynorm*vlwR + znorm*wlwR 312 ubarlpR = xnorm*ulpR + ynorm*vlpR + znorm*wlpR 313! 314! Get variables on "right" side of face 315! 316 rhor = qnode(1,node2) 317 rhorrL = 0. 318 rhoruL = 0. 319 rhorvL = 0. 320 rhorwL = 0. 321 rhorpL = 0. 322! 323 rhorrR = 1. 324 rhoruR = 0. 325 rhorvR = 0. 326 rhorwR = 0. 327 rhorpR = 0. 328! 329 ur = qnode(2,node2) 330 urrL = 0.0 331 uruL = 0.0 332 urvL = 0.0 333 urwL = 0.0 334 urpL = 0.0 335! 336 urrR = 0.0 337 uruR = 1.0 338 urvR = 0.0 339 urwR = 0.0 340 urpR = 0.0 341! 342 vr = qnode(3,node2) 343 vrrL = 0.0 344 vruL = 0.0 345 vrvL = 0.0 346 vrwL = 0.0 347 vrpL = 0.0 348! 349 vrrR = 0.0 350 vruR = 0.0 351 vrvR = 1.0 352 vrwR = 0.0 353 vrpR = 0.0 354! 355 wr = qnode(4,node2) 356 wrrL = 0.0 357 wruL = 0.0 358 wrvL = 0.0 359 wrwL = 0.0 360 wrpL = 0.0 361! 362 wrrR = 0.0 363 wruR = 0.0 364 wrvR = 0.0 365 wrwR = 1.0 366 wrpR = 0.0 367! 368 q2r = ur*ur + vr*vr + wr*wr 369 q2rrL = 2.*ur*urrL + 2.*vr*vrrL + 2.*wr*wrrL 370 q2ruL = 2.*ur*uruL + 2.*vr*vruL + 2.*wr*wruL 371 q2rvL = 2.*ur*urvL + 2.*vr*vrvL + 2.*wr*wrvL 372 q2rwL = 2.*ur*urwL + 2.*vr*vrwL + 2.*wr*wrwL 373 q2rpL = 2.*ur*urpL + 2.*vr*vrpL + 2.*wr*wrpL 374! 375 q2rrR = 2.*ur*urrR + 2.*vr*vrrR + 2.*wr*wrrR 376 q2ruR = 2.*ur*uruR + 2.*vr*vruR + 2.*wr*wruR 377 q2rvR = 2.*ur*urvR + 2.*vr*vrvR + 2.*wr*wrvR 378 q2rwR = 2.*ur*urwR + 2.*vr*vrwR + 2.*wr*wrwR 379 q2rpR = 2.*ur*urpR + 2.*vr*vrpR + 2.*wr*wrpR 380! 381 pressr = qnode(5,node2) 382 pressrrL = 0. 383 pressruL = 0. 384 pressrvL = 0. 385 pressrwL = 0. 386 pressrpL = 0. 387! 388 pressrrR = 0. 389 pressruR = 0. 390 pressrvR = 0. 391 pressrwR = 0. 392 pressrpR = 1. 393! 394 enrgyr = pressr/gm1 + .5*rhor*q2r 395! 396 enrgyrrL = pressrrL/gm1 + .5*(rhor*q2rrL + q2r*rhorrL) 397 enrgyruL = pressruL/gm1 + .5*(rhor*q2ruL + q2r*rhoruL) 398 enrgyrvL = pressrvL/gm1 + .5*(rhor*q2rvL + q2r*rhorvL) 399 enrgyrwL = pressrwL/gm1 + .5*(rhor*q2rwL + q2r*rhorwL) 400 enrgyrpL = pressrpL/gm1 + .5*(rhor*q2rpL + q2r*rhorpL) 401! 402 enrgyrrR = pressrrR/gm1 + .5*(rhor*q2rrR + q2r*rhorrR) 403 enrgyruR = pressruR/gm1 + .5*(rhor*q2ruR + q2r*rhoruR) 404 enrgyrvR = pressrvR/gm1 + .5*(rhor*q2rvR + q2r*rhorvR) 405 enrgyrwR = pressrwR/gm1 + .5*(rhor*q2rwR + q2r*rhorwR) 406 enrgyrpR = pressrpR/gm1 + .5*(rhor*q2rpR + q2r*rhorpR) 407! 408 Hr = (enrgyr + pressr)/rhor 409! 410 HrrL = (rhor*(enrgyrrL+pressrrL) - (enrgyr+pressr)*rhorrL) 411 & / rhor / rhor 412 HruL = (rhor*(enrgyruL+pressruL) - (enrgyr+pressr)*rhoruL) 413 & / rhor / rhor 414 HrvL = (rhor*(enrgyrvL+pressrvL) - (enrgyr+pressr)*rhorvL) 415 & / rhor / rhor 416 HrwL = (rhor*(enrgyrwL+pressrwL) - (enrgyr+pressr)*rhorwL) 417 & / rhor / rhor 418 HrpL = (rhor*(enrgyrpL+pressrpL) - (enrgyr+pressr)*rhorpL) 419 & / rhor / rhor 420! 421 HrrR = (rhor*(enrgyrrR+pressrrR) - (enrgyr+pressr)*rhorrR) 422 & / rhor / rhor 423 HruR = (rhor*(enrgyruR+pressruR) - (enrgyr+pressr)*rhoruR) 424 & / rhor / rhor 425 HrvR = (rhor*(enrgyrvR+pressrvR) - (enrgyr+pressr)*rhorvR) 426 & / rhor / rhor 427 HrwR = (rhor*(enrgyrwR+pressrwR) - (enrgyr+pressr)*rhorwR) 428 & / rhor / rhor 429 HrpR = (rhor*(enrgyrpR+pressrpR) - (enrgyr+pressr)*rhorpR) 430 & / rhor / rhor 431! 432! 433 cr = sqrt(gamma*pressr/rhor) 434! 435 crrL = 0.5 / sqrt(gamma*pressr/rhor) * 436 & gamma * (rhor*pressrrL - pressr*rhorrL) 437 & / rhor / rhor 438 cruL = 0.5 / sqrt(gamma*pressr/rhor) * 439 & gamma * (rhor*pressruL - pressr*rhoruL) 440 & / rhor / rhor 441 crvL = 0.5 / sqrt(gamma*pressr/rhor) * 442 & gamma * (rhor*pressrvL - pressr*rhorvL) 443 & / rhor / rhor 444 crwL = 0.5 / sqrt(gamma*pressr/rhor) * 445 & gamma * (rhor*pressrwL - pressr*rhorwL) 446 & / rhor / rhor 447 crpL = 0.5 / sqrt(gamma*pressr/rhor) * 448 & gamma * (rhor*pressrpL - pressr*rhorpL) 449 & / rhor / rhor 450! 451 crrR = 0.5 / sqrt(gamma*pressr/rhor) * 452 & gamma * (rhor*pressrrR - pressr*rhorrR) 453 & / rhor / rhor 454 cruR = 0.5 / sqrt(gamma*pressr/rhor) * 455 & gamma * (rhor*pressruR - pressr*rhoruR) 456 & / rhor / rhor 457 crvR = 0.5 / sqrt(gamma*pressr/rhor) * 458 & gamma * (rhor*pressrvR - pressr*rhorvR) 459 & / rhor / rhor 460 crwR = 0.5 / sqrt(gamma*pressr/rhor) * 461 & gamma * (rhor*pressrwR - pressr*rhorwR) 462 & / rhor / rhor 463 crpR = 0.5 / sqrt(gamma*pressr/rhor) * 464 & gamma * (rhor*pressrpR - pressr*rhorpR) 465 & / rhor / rhor 466! 467! 468 ubarr = xnorm*ur + ynorm*vr + znorm*wr 469! 470 ubarrrL = xnorm*urrL + ynorm*vrrL + znorm*wrrL 471 ubarruL = xnorm*uruL + ynorm*vruL + znorm*wruL 472 ubarrvL = xnorm*urvL + ynorm*vrvL + znorm*wrvL 473 ubarrwL = xnorm*urwL + ynorm*vrwL + znorm*wrwL 474 ubarrpL = xnorm*urpL + ynorm*vrpL + znorm*wrpL 475! 476 ubarrrR = xnorm*urrR + ynorm*vrrR + znorm*wrrR 477 ubarruR = xnorm*uruR + ynorm*vruR + znorm*wruR 478 ubarrvR = xnorm*urvR + ynorm*vrvR + znorm*wrvR 479 ubarrwR = xnorm*urwR + ynorm*vrwR + znorm*wrwR 480 ubarrpR = xnorm*urpR + ynorm*vrpR + znorm*wrpR 481! 482! Compute rho averages 483! 484 rho = sqrt(rhol*rhor) 485! 486 rhorL = 0.5 / sqrt(rhol*rhor) * (rhol*rhorrL + 487 & rhor*rholrL) 488 rhouL = 0.5 / sqrt(rhol*rhor) * (rhol*rhoruL + 489 & rhor*rholuL) 490 rhovL = 0.5 / sqrt(rhol*rhor) * (rhol*rhorvL + 491 & rhor*rholvL) 492 rhowL = 0.5 / sqrt(rhol*rhor) * (rhol*rhorwL + 493 & rhor*rholwL) 494 rhopL = 0.5 / sqrt(rhol*rhor) * (rhol*rhorpL + 495 & rhor*rholpL) 496! 497 rhorR = 0.5 / sqrt(rhol*rhor) * (rhol*rhorrR + 498 & rhor*rholrR) 499 rhouR = 0.5 / sqrt(rhol*rhor) * (rhol*rhoruR + 500 & rhor*rholuR) 501 rhovR = 0.5 / sqrt(rhol*rhor) * (rhol*rhorvR + 502 & rhor*rholvR) 503 rhowR = 0.5 / sqrt(rhol*rhor) * (rhol*rhorwR + 504 & rhor*rholwR) 505 rhopR = 0.5 / sqrt(rhol*rhor) * (rhol*rhorpR + 506 & rhor*rholpR) 507! 508 wat = rho/(rho + rhor) 509! 510 watrL = ((rho + rhor)*rhorL - rho*(rhorL + rhorrL)) / 511 & (rho + rhor) / (rho + rhor) 512 watuL = ((rho + rhor)*rhouL - rho*(rhouL + rhoruL)) / 513 & (rho + rhor) / (rho + rhor) 514 watvL = ((rho + rhor)*rhovL - rho*(rhovL + rhorvL)) / 515 & (rho + rhor) / (rho + rhor) 516 watwL = ((rho + rhor)*rhowL - rho*(rhowL + rhorwL)) / 517 & (rho + rhor) / (rho + rhor) 518 watpL = ((rho + rhor)*rhopL - rho*(rhopL + rhorpL)) / 519 & (rho + rhor) / (rho + rhor) 520! 521 watrR = ((rho + rhor)*rhorR - rho*(rhorR + rhorrR)) / 522 & (rho + rhor) / (rho + rhor) 523 watuR = ((rho + rhor)*rhouR - rho*(rhouR + rhoruR)) / 524 & (rho + rhor) / (rho + rhor) 525 watvR = ((rho + rhor)*rhovR - rho*(rhovR + rhorvR)) / 526 & (rho + rhor) / (rho + rhor) 527 watwR = ((rho + rhor)*rhowR - rho*(rhowR + rhorwR)) / 528 & (rho + rhor) / (rho + rhor) 529 watpR = ((rho + rhor)*rhopR - rho*(rhopR + rhorpR)) / 530 & (rho + rhor) / (rho + rhor) 531! 532!234567890c234567890c234567890c234567890c234567890c234567890c23456789012 533! 534 u = ul*wat + ur*(1. - wat) 535! 536 urL = ul*watrL + wat*ulrL + ur*(-watrL) + (1. - wat)*urrL 537 uuL = ul*watuL + wat*uluL + ur*(-watuL) + (1. - wat)*uruL 538 uvL = ul*watvL + wat*ulvL + ur*(-watvL) + (1. - wat)*urvL 539 uwL = ul*watwL + wat*ulwL + ur*(-watwL) + (1. - wat)*urwL 540 upL = ul*watpL + wat*ulpL + ur*(-watpL) + (1. - wat)*urpL 541! 542 urR = ul*watrR + wat*ulrR + ur*(-watrR) + (1. - wat)*urrR 543 uuR = ul*watuR + wat*uluR + ur*(-watuR) + (1. - wat)*uruR 544 uvR = ul*watvR + wat*ulvR + ur*(-watvR) + (1. - wat)*urvR 545 uwR = ul*watwR + wat*ulwR + ur*(-watwR) + (1. - wat)*urwR 546 upR = ul*watpR + wat*ulpR + ur*(-watpR) + (1. - wat)*urpR 547! 548 v = vl*wat + vr*(1. - wat) 549! 550 vrL = vl*watrL + wat*vlrL + vr*(-watrL) + (1. - wat)*vrrL 551 vuL = vl*watuL + wat*vluL + vr*(-watuL) + (1. - wat)*vruL 552 vvL = vl*watvL + wat*vlvL + vr*(-watvL) + (1. - wat)*vrvL 553 vwL = vl*watwL + wat*vlwL + vr*(-watwL) + (1. - wat)*vrwL 554 vpL = vl*watpL + wat*vlpL + vr*(-watpL) + (1. - wat)*vrpL 555! 556 vrR = vl*watrR + wat*vlrR + vr*(-watrR) + (1. - wat)*vrrR 557 vuR = vl*watuR + wat*vluR + vr*(-watuR) + (1. - wat)*vruR 558 vvR = vl*watvR + wat*vlvR + vr*(-watvR) + (1. - wat)*vrvR 559 vwR = vl*watwR + wat*vlwR + vr*(-watwR) + (1. - wat)*vrwR 560 vpR = vl*watpR + wat*vlpR + vr*(-watpR) + (1. - wat)*vrpR 561! 562 w = wl*wat + wr*(1. - wat) 563! 564 wrL = wl*watrL + wat*wlrL + wr*(-watrL) + (1. - wat)*wrrL 565 wuL = wl*watuL + wat*wluL + wr*(-watuL) + (1. - wat)*wruL 566 wvL = wl*watvL + wat*wlvL + wr*(-watvL) + (1. - wat)*wrvL 567 wwL = wl*watwL + wat*wlwL + wr*(-watwL) + (1. - wat)*wrwL 568 wpL = wl*watpL + wat*wlpL + wr*(-watpL) + (1. - wat)*wrpL 569! 570 wrR = wl*watrR + wat*wlrR + wr*(-watrR) + (1. - wat)*wrrR 571 wuR = wl*watuR + wat*wluR + wr*(-watuR) + (1. - wat)*wruR 572 wvR = wl*watvR + wat*wlvR + wr*(-watvR) + (1. - wat)*wrvR 573 wwR = wl*watwR + wat*wlwR + wr*(-watwR) + (1. - wat)*wrwR 574 wpR = wl*watpR + wat*wlpR + wr*(-watpR) + (1. - wat)*wrpR 575! 576 H = Hl*wat + Hr*(1. - wat) 577! 578 HrL = Hl*watrL + wat*HlrL + Hr*(-watrL) + (1. - wat)*HrrL 579 HuL = Hl*watuL + wat*HluL + Hr*(-watuL) + (1. - wat)*HruL 580 HvL = Hl*watvL + wat*HlvL + Hr*(-watvL) + (1. - wat)*HrvL 581 HwL = Hl*watwL + wat*HlwL + Hr*(-watwL) + (1. - wat)*HrwL 582 HpL = Hl*watpL + wat*HlpL + Hr*(-watpL) + (1. - wat)*HrpL 583! 584 HrR = Hl*watrR + wat*HlrR + Hr*(-watrR) + (1. - wat)*HrrR 585 HuR = Hl*watuR + wat*HluR + Hr*(-watuR) + (1. - wat)*HruR 586 HvR = Hl*watvR + wat*HlvR + Hr*(-watvR) + (1. - wat)*HrvR 587 HwR = Hl*watwR + wat*HlwR + Hr*(-watwR) + (1. - wat)*HrwR 588 HpR = Hl*watpR + wat*HlpR + Hr*(-watpR) + (1. - wat)*HrpR 589! 590 q2 = u*u + v*v + w*w 591! 592 q2rL = 2.*u*urL + 2.*v*vrL + 2.*w*wrL 593 q2uL = 2.*u*uuL + 2.*v*vuL + 2.*w*wuL 594 q2vL = 2.*u*uvL + 2.*v*vvL + 2.*w*wvL 595 q2wL = 2.*u*uwL + 2.*v*vwL + 2.*w*wwL 596 q2pL = 2.*u*upL + 2.*v*vpL + 2.*w*wpL 597! 598 q2rR = 2.*u*urR + 2.*v*vrR + 2.*w*wrR 599 q2uR = 2.*u*uuR + 2.*v*vuR + 2.*w*wuR 600 q2vR = 2.*u*uvR + 2.*v*vvR + 2.*w*wvR 601 q2wR = 2.*u*uwR + 2.*v*vwR + 2.*w*wwR 602 q2pR = 2.*u*upR + 2.*v*vpR + 2.*w*wpR 603! 604 c = sqrt(gm1*(H - 0.5*q2)) 605! 606 crL = 0.5 / c * gm1*(HrL - 0.5*q2rL) 607 cuL = 0.5 / c * gm1*(HuL - 0.5*q2uL) 608 cvL = 0.5 / c * gm1*(HvL - 0.5*q2vL) 609 cwL = 0.5 / c * gm1*(HwL - 0.5*q2wL) 610 cpL = 0.5 / c * gm1*(HpL - 0.5*q2pL) 611! 612 crR = 0.5 / c * gm1*(HrR - 0.5*q2rR) 613 cuR = 0.5 / c * gm1*(HuR - 0.5*q2uR) 614 cvR = 0.5 / c * gm1*(HvR - 0.5*q2vR) 615 cwR = 0.5 / c * gm1*(HwR - 0.5*q2wR) 616 cpR = 0.5 / c * gm1*(HpR - 0.5*q2pR) 617! 618 ubar = xnorm*u + ynorm*v + znorm*w 619! 620 ubarrL = xnorm*urL + ynorm*vrL + znorm*wrL 621 ubaruL = xnorm*uuL + ynorm*vuL + znorm*wuL 622 ubarvL = xnorm*uvL + ynorm*vvL + znorm*wvL 623 ubarwL = xnorm*uwL + ynorm*vwL + znorm*wwL 624 ubarpL = xnorm*upL + ynorm*vpL + znorm*wpL 625! 626 ubarrR = xnorm*urR + ynorm*vrR + znorm*wrR 627 ubaruR = xnorm*uuR + ynorm*vuR + znorm*wuR 628 ubarvR = xnorm*uvR + ynorm*vvR + znorm*wvR 629 ubarwR = xnorm*uwR + ynorm*vwR + znorm*wwR 630 ubarpR = xnorm*upR + ynorm*vpR + znorm*wpR 631! 632! Now compute eigenvalues, eigenvectors, and strengths 633! 634 eig1 = abs(ubar + c) 635 eig2 = abs(ubar - c) 636 eig3 = abs(ubar) 637 flops = flops + 1521.0 638! 639 if (ubar+c.gt.0) then 640 eig1rL = ubarrL + crL 641 eig1uL = ubaruL + cuL 642 eig1vL = ubarvL + cvL 643 eig1wL = ubarwL + cwL 644 eig1pL = ubarpL + cpL 645 eig1rR = ubarrR + crR 646 eig1uR = ubaruR + cuR 647 eig1vR = ubarvR + cvR 648 eig1wR = ubarwR + cwR 649 eig1pR = ubarpR + cpR 650 flops = flops + 10.0 651 else 652 eig1rL = -(ubarrL + crL) 653 eig1uL = -(ubaruL + cuL) 654 eig1vL = -(ubarvL + cvL) 655 eig1wL = -(ubarwL + cwL) 656 eig1pL = -(ubarpL + cpL) 657 eig1rR = -(ubarrR + crR) 658 eig1uR = -(ubaruR + cuR) 659 eig1vR = -(ubarvR + cvR) 660 eig1wR = -(ubarwR + cwR) 661 eig1pR = -(ubarpR + cpR) 662 flops = flops + 20.0 663 endif 664! 665 if (ubar-c.gt.0) then 666 eig2rL = ubarrL - crL 667 eig2uL = ubaruL - cuL 668 eig2vL = ubarvL - cvL 669 eig2wL = ubarwL - cwL 670 eig2pL = ubarpL - cpL 671 eig2rR = ubarrR - crR 672 eig2uR = ubaruR - cuR 673 eig2vR = ubarvR - cvR 674 eig2wR = ubarwR - cwR 675 eig2pR = ubarpR - cpR 676 flops = flops + 10.0 677 else 678 eig2rL = -(ubarrL - crL) 679 eig2uL = -(ubaruL - cuL) 680 eig2vL = -(ubarvL - cvL) 681 eig2wL = -(ubarwL - cwL) 682 eig2pL = -(ubarpL - cpL) 683 eig2rR = -(ubarrR - crR) 684 eig2uR = -(ubaruR - cuR) 685 eig2vR = -(ubarvR - cvR) 686 eig2wR = -(ubarwR - cwR) 687 eig2pR = -(ubarpR - cpR) 688 flops = flops + 20.0 689 endif 690! 691 if (ubar.gt.0) then 692 eig3rL = ubarrL 693 eig3uL = ubaruL 694 eig3vL = ubarvL 695 eig3wL = ubarwL 696 eig3pL = ubarpL 697 eig3rR = ubarrR 698 eig3uR = ubaruR 699 eig3vR = ubarvR 700 eig3wR = ubarwR 701 eig3pR = ubarpR 702 else 703 eig3rL = -(ubarrL) 704 eig3uL = -(ubaruL) 705 eig3vL = -(ubarvL) 706 eig3wL = -(ubarwL) 707 eig3pL = -(ubarpL) 708 eig3rR = -(ubarrR) 709 eig3uR = -(ubaruR) 710 eig3vR = -(ubarvR) 711 eig3wR = -(ubarwR) 712 eig3pR = -(ubarpR) 713 flops = flops + 10.0 714 endif 715! 716 drho = rhor - rhol 717! 718 drhorL = rhorrL - rholrL 719 drhouL = rhoruL - rholuL 720 drhovL = rhorvL - rholvL 721 drhowL = rhorwL - rholwL 722 drhopL = rhorpL - rholpL 723! 724 drhorR = rhorrR - rholrR 725 drhouR = rhoruR - rholuR 726 drhovR = rhorvR - rholvR 727 drhowR = rhorwR - rholwR 728 drhopR = rhorpR - rholpR 729! 730 dpress = pressr - pressl 731! 732 dpressrL = pressrrL - presslrL 733 dpressuL = pressruL - pressluL 734 dpressvL = pressrvL - presslvL 735 dpresswL = pressrwL - presslwL 736 dpresspL = pressrpL - presslpL 737! 738 dpressrR = pressrrR - presslrR 739 dpressuR = pressruR - pressluR 740 dpressvR = pressrvR - presslvR 741 dpresswR = pressrwR - presslwR 742 dpresspR = pressrpR - presslpR 743! 744 du = ur - ul 745! 746 durL = urrL - ulrL 747 duuL = uruL - uluL 748 duvL = urvL - ulvL 749 duwL = urwL - ulwL 750 dupL = urpL - ulpL 751! 752 durR = urrR - ulrR 753 duuR = uruR - uluR 754 duvR = urvR - ulvR 755 duwR = urwR - ulwR 756 dupR = urpR - ulpR 757! 758 dv = vr - vl 759! 760 dvrL = vrrL - vlrL 761 dvuL = vruL - vluL 762 dvvL = vrvL - vlvL 763 dvwL = vrwL - vlwL 764 dvpL = vrpL - vlpL 765! 766 dvrR = vrrR - vlrR 767 dvuR = vruR - vluR 768 dvvR = vrvR - vlvR 769 dvwR = vrwR - vlwR 770 dvpR = vrpR - vlpR 771! 772 dw = wr - wl 773! 774 dwrL = wrrL - wlrL 775 dwuL = wruL - wluL 776 dwvL = wrvL - wlvL 777 dwwL = wrwL - wlwL 778 dwpL = wrpL - wlpL 779! 780 dwrR = wrrR - wlrR 781 dwuR = wruR - wluR 782 dwvR = wrvR - wlvR 783 dwwR = wrwR - wlwR 784 dwpR = wrpR - wlpR 785! 786 dubar = ubarr - ubarl 787! 788 dubarrL = ubarrrL - ubarlrL 789 dubaruL = ubarruL - ubarluL 790 dubarvL = ubarrvL - ubarlvL 791 dubarwL = ubarrwL - ubarlwL 792 dubarpL = ubarrpL - ubarlpL 793! 794 dubarrR = ubarrrR - ubarlrR 795 dubaruR = ubarruR - ubarluR 796 dubarvR = ubarrvR - ubarlvR 797 dubarwR = ubarrwR - ubarlwR 798 dubarpR = ubarrpR - ubarlpR 799! 800 c2 = c*c 801! 802 c2rL = 2. * c * crL 803 c2uL = 2. * c * cuL 804 c2vL = 2. * c * cvL 805 c2wL = 2. * c * cwL 806 c2pL = 2. * c * cpL 807! 808 c2rR = 2. * c * crR 809 c2uR = 2. * c * cuR 810 c2vR = 2. * c * cvR 811 c2wR = 2. * c * cwR 812 c2pR = 2. * c * cpR 813! 814! jumps have units of density 815! 816 dv1 = 0.5*(dpress + rho*c*dubar)/c2 817! 818 dv1rL = 0.5*(c2*(dpressrL + rho*(c*dubarrL + dubar*crL) 819 & + c*dubar*rhorL) - (dpress + rho*c*dubar)*c2rL) 820 & / c2 / c2 821 dv1uL = 0.5*(c2*(dpressuL + rho*(c*dubaruL + dubar*cuL) 822 & + c*dubar*rhouL) - (dpress + rho*c*dubar)*c2uL) 823 & / c2 / c2 824 dv1vL = 0.5*(c2*(dpressvL + rho*(c*dubarvL + dubar*cvL) 825 & + c*dubar*rhovL) - (dpress + rho*c*dubar)*c2vL) 826 & / c2 / c2 827 dv1wL = 0.5*(c2*(dpresswL + rho*(c*dubarwL + dubar*cwL) 828 & + c*dubar*rhowL) - (dpress + rho*c*dubar)*c2wL) 829 & / c2 / c2 830 dv1pL = 0.5*(c2*(dpresspL + rho*(c*dubarpL + dubar*cpL) 831 & + c*dubar*rhopL) - (dpress + rho*c*dubar)*c2pL) 832 & / c2 / c2 833! 834 dv1rR = 0.5*(c2*(dpressrR + rho*(c*dubarrR + dubar*crR) 835 & + c*dubar*rhorR) - (dpress + rho*c*dubar)*c2rR) 836 & / c2 / c2 837 dv1uR = 0.5*(c2*(dpressuR + rho*(c*dubaruR + dubar*cuR) 838 & + c*dubar*rhouR) - (dpress + rho*c*dubar)*c2uR) 839 & / c2 / c2 840 dv1vR = 0.5*(c2*(dpressvR + rho*(c*dubarvR + dubar*cvR) 841 & + c*dubar*rhovR) - (dpress + rho*c*dubar)*c2vR) 842 & / c2 / c2 843 dv1wR = 0.5*(c2*(dpresswR + rho*(c*dubarwR + dubar*cwR) 844 & + c*dubar*rhowR) - (dpress + rho*c*dubar)*c2wR) 845 & / c2 / c2 846 dv1pR = 0.5*(c2*(dpresspR + rho*(c*dubarpR + dubar*cpR) 847 & + c*dubar*rhopR) - (dpress + rho*c*dubar)*c2pR) 848 & / c2 / c2 849! 850! 851 dv2 = 0.5*(dpress - rho*c*dubar)/c2 852! 853 dv2rL = 0.5*(c2*(dpressrL - rho*(c*dubarrL + dubar*crL) 854 & - c*dubar*rhorL) - (dpress - rho*c*dubar)*c2rL) 855 & / c2 / c2 856 dv2uL = 0.5*(c2*(dpressuL - rho*(c*dubaruL + dubar*cuL) 857 & - c*dubar*rhouL) - (dpress - rho*c*dubar)*c2uL) 858 & / c2 / c2 859 dv2vL = 0.5*(c2*(dpressvL - rho*(c*dubarvL + dubar*cvL) 860 & - c*dubar*rhovL) - (dpress - rho*c*dubar)*c2vL) 861 & / c2 / c2 862 dv2wL = 0.5*(c2*(dpresswL - rho*(c*dubarwL + dubar*cwL) 863 & - c*dubar*rhowL) - (dpress - rho*c*dubar)*c2wL) 864 & / c2 / c2 865 dv2pL = 0.5*(c2*(dpresspL - rho*(c*dubarpL + dubar*cpL) 866 & - c*dubar*rhopL) - (dpress - rho*c*dubar)*c2pL) 867 & / c2 / c2 868! 869 dv2rR = 0.5*(c2*(dpressrR - rho*(c*dubarrR + dubar*crR) 870 & - c*dubar*rhorR) - (dpress - rho*c*dubar)*c2rR) 871 & / c2 / c2 872 dv2uR = 0.5*(c2*(dpressuR - rho*(c*dubaruR + dubar*cuR) 873 & - c*dubar*rhouR) - (dpress - rho*c*dubar)*c2uR) 874 & / c2 / c2 875 dv2vR = 0.5*(c2*(dpressvR - rho*(c*dubarvR + dubar*cvR) 876 & - c*dubar*rhovR) - (dpress - rho*c*dubar)*c2vR) 877 & / c2 / c2 878 dv2wR = 0.5*(c2*(dpresswR - rho*(c*dubarwR + dubar*cwR) 879 & - c*dubar*rhowR) - (dpress - rho*c*dubar)*c2wR) 880 & / c2 / c2 881 dv2pR = 0.5*(c2*(dpresspR - rho*(c*dubarpR + dubar*cpR) 882 & - c*dubar*rhopR) - (dpress - rho*c*dubar)*c2pR) 883 & / c2 / c2 884! 885 dv3 = rho 886! 887 dv3rL = rhorL 888 dv3uL = rhouL 889 dv3vL = rhovL 890 dv3wL = rhowL 891 dv3pL = rhopL 892! 893 dv3rR = rhorR 894 dv3uR = rhouR 895 dv3vR = rhovR 896 dv3wR = rhowR 897 dv3pR = rhopR 898! 899 dv4 = (c*c*drho - dpress)/c2 900! 901 dv4rL = (c2*((c*(c*drhorL+drho*crL)+c*drho*crL) - dpressrL) 902 & - (c*c*drho - dpress)*c2rL) / c2 / c2 903 dv4uL = (c2*((c*(c*drhouL+drho*cuL)+c*drho*cuL) - dpressuL) 904 & - (c*c*drho - dpress)*c2uL) / c2 / c2 905 dv4vL = (c2*((c*(c*drhovL+drho*cvL)+c*drho*cvL) - dpressvL) 906 & - (c*c*drho - dpress)*c2vL) / c2 / c2 907 dv4wL = (c2*((c*(c*drhowL+drho*cwL)+c*drho*cwL) - dpresswL) 908 & - (c*c*drho - dpress)*c2wL) / c2 / c2 909 dv4pL = (c2*((c*(c*drhopL+drho*cpL)+c*drho*cpL) - dpresspL) 910 & - (c*c*drho - dpress)*c2pL) / c2 / c2 911! 912 dv4rR = (c2*((c*(c*drhorR+drho*crR)+c*drho*crR) - dpressrR) 913 & - (c*c*drho - dpress)*c2rR) / c2 / c2 914 dv4uR = (c2*((c*(c*drhouR+drho*cuR)+c*drho*cuR) - dpressuR) 915 & - (c*c*drho - dpress)*c2uR) / c2 / c2 916 dv4vR = (c2*((c*(c*drhovR+drho*cvR)+c*drho*cvR) - dpressvR) 917 & - (c*c*drho - dpress)*c2vR) / c2 / c2 918 dv4wR = (c2*((c*(c*drhowR+drho*cwR)+c*drho*cwR) - dpresswR) 919 & - (c*c*drho - dpress)*c2wR) / c2 / c2 920 dv4pR = (c2*((c*(c*drhopR+drho*cpR)+c*drho*cpR) - dpresspR) 921 & - (c*c*drho - dpress)*c2pR) / c2 / c2 922! 923 r21 = u + c*xnorm 924! 925 r21rL = urL + xnorm*crL 926 r21uL = uuL + xnorm*cuL 927 r21vL = uvL + xnorm*cvL 928 r21wL = uwL + xnorm*cwL 929 r21pL = upL + xnorm*cpL 930! 931 r21rR = urR + xnorm*crR 932 r21uR = uuR + xnorm*cuR 933 r21vR = uvR + xnorm*cvR 934 r21wR = uwR + xnorm*cwR 935 r21pR = upR + xnorm*cpR 936! 937 r31 = v + c*ynorm 938! 939 r31rL = vrL + ynorm*crL 940 r31uL = vuL + ynorm*cuL 941 r31vL = vvL + ynorm*cvL 942 r31wL = vwL + ynorm*cwL 943 r31pL = vpL + ynorm*cpL 944! 945 r31rR = vrR + ynorm*crR 946 r31uR = vuR + ynorm*cuR 947 r31vR = vvR + ynorm*cvR 948 r31wR = vwR + ynorm*cwR 949 r31pR = vpR + ynorm*cpR 950! 951 r41 = w + c*znorm 952! 953 r41rL = wrL + znorm*crL 954 r41uL = wuL + znorm*cuL 955 r41vL = wvL + znorm*cvL 956 r41wL = wwL + znorm*cwL 957 r41pL = wpL + znorm*cpL 958! 959 r41rR = wrR + znorm*crR 960 r41uR = wuR + znorm*cuR 961 r41vR = wvR + znorm*cvR 962 r41wR = wwR + znorm*cwR 963 r41pR = wpR + znorm*cpR 964! 965 r51 = H + c*ubar 966! 967 r51rL = HrL + c*ubarrL + ubar*crL 968 r51uL = HuL + c*ubaruL + ubar*cuL 969 r51vL = HvL + c*ubarvL + ubar*cvL 970 r51wL = HwL + c*ubarwL + ubar*cwL 971 r51pL = HpL + c*ubarpL + ubar*cpL 972! 973 r51rR = HrR + c*ubarrR + ubar*crR 974 r51uR = HuR + c*ubaruR + ubar*cuR 975 r51vR = HvR + c*ubarvR + ubar*cvR 976 r51wR = HwR + c*ubarwR + ubar*cwR 977 r51pR = HpR + c*ubarpR + ubar*cpR 978! 979 r22 = u - c*xnorm 980! 981 r22rL = urL - xnorm*crL 982 r22uL = uuL - xnorm*cuL 983 r22vL = uvL - xnorm*cvL 984 r22wL = uwL - xnorm*cwL 985 r22pL = upL - xnorm*cpL 986! 987 r22rR = urR - xnorm*crR 988 r22uR = uuR - xnorm*cuR 989 r22vR = uvR - xnorm*cvR 990 r22wR = uwR - xnorm*cwR 991 r22pR = upR - xnorm*cpR 992! 993 r32 = v - c*ynorm 994! 995 r32rL = vrL - ynorm*crL 996 r32uL = vuL - ynorm*cuL 997 r32vL = vvL - ynorm*cvL 998 r32wL = vwL - ynorm*cwL 999 r32pL = vpL - ynorm*cpL 1000! 1001 r32rR = vrR - ynorm*crR 1002 r32uR = vuR - ynorm*cuR 1003 r32vR = vvR - ynorm*cvR 1004 r32wR = vwR - ynorm*cwR 1005 r32pR = vpR - ynorm*cpR 1006! 1007 r42 = w - c*znorm 1008! 1009 r42rL = wrL - znorm*crL 1010 r42uL = wuL - znorm*cuL 1011 r42vL = wvL - znorm*cvL 1012 r42wL = wwL - znorm*cwL 1013 r42pL = wpL - znorm*cpL 1014! 1015 r42rR = wrR - znorm*crR 1016 r42uR = wuR - znorm*cuR 1017 r42vR = wvR - znorm*cvR 1018 r42wR = wwR - znorm*cwR 1019 r42pR = wpR - znorm*cpR 1020! 1021 r52 = H - c*ubar 1022! 1023 r52rL = HrL - c*ubarrL - ubar*crL 1024 r52uL = HuL - c*ubaruL - ubar*cuL 1025 r52vL = HvL - c*ubarvL - ubar*cvL 1026 r52wL = HwL - c*ubarwL - ubar*cwL 1027 r52pL = HpL - c*ubarpL - ubar*cpL 1028! 1029 r52rR = HrR - c*ubarrR - ubar*crR 1030 r52uR = HuR - c*ubaruR - ubar*cuR 1031 r52vR = HvR - c*ubarvR - ubar*cvR 1032 r52wR = HwR - c*ubarwR - ubar*cwR 1033 r52pR = HpR - c*ubarpR - ubar*cpR 1034! 1035 r23 = du - dubar*xnorm 1036! 1037 r23rL = durL - xnorm*dubarrL 1038 r23uL = duuL - xnorm*dubaruL 1039 r23vL = duvL - xnorm*dubarvL 1040 r23wL = duwL - xnorm*dubarwL 1041 r23pL = dupL - xnorm*dubarpL 1042! 1043 r23rR = durR - xnorm*dubarrR 1044 r23uR = duuR - xnorm*dubaruR 1045 r23vR = duvR - xnorm*dubarvR 1046 r23wR = duwR - xnorm*dubarwR 1047 r23pR = dupR - xnorm*dubarpR 1048! 1049 r33 = dv - dubar*ynorm 1050! 1051 r33rL = dvrL - ynorm*dubarrL 1052 r33uL = dvuL - ynorm*dubaruL 1053 r33vL = dvvL - ynorm*dubarvL 1054 r33wL = dvwL - ynorm*dubarwL 1055 r33pL = dvpL - ynorm*dubarpL 1056! 1057 r33rR = dvrR - ynorm*dubarrR 1058 r33uR = dvuR - ynorm*dubaruR 1059 r33vR = dvvR - ynorm*dubarvR 1060 r33wR = dvwR - ynorm*dubarwR 1061 r33pR = dvpR - ynorm*dubarpR 1062! 1063 r43 = dw - dubar*znorm 1064! 1065 r43rL = dwrL - znorm*dubarrL 1066 r43uL = dwuL - znorm*dubaruL 1067 r43vL = dwvL - znorm*dubarvL 1068 r43wL = dwwL - znorm*dubarwL 1069 r43pL = dwpL - znorm*dubarpL 1070! 1071 r43rR = dwrR - znorm*dubarrR 1072 r43uR = dwuR - znorm*dubaruR 1073 r43vR = dwvR - znorm*dubarvR 1074 r43wR = dwwR - znorm*dubarwR 1075 r43pR = dwpR - znorm*dubarpR 1076! 1077 r53 = u*du + v*dv + w*dw - ubar*dubar 1078! 1079 r53rL = u*durL+du*urL + v*dvrL+dv*vrL + w*dwrL+dw*wrL 1080 & - ubar*dubarrL - dubar*ubarrL 1081 r53uL = u*duuL+du*uuL + v*dvuL+dv*vuL + w*dwuL+dw*wuL 1082 & - ubar*dubaruL - dubar*ubaruL 1083 r53vL = u*duvL+du*uvL + v*dvvL+dv*vvL + w*dwvL+dw*wvL 1084 & - ubar*dubarvL - dubar*ubarvL 1085 r53wL = u*duwL+du*uwL + v*dvwL+dv*vwL + w*dwwL+dw*wwL 1086 & - ubar*dubarwL - dubar*ubarwL 1087 r53pL = u*dupL+du*upL + v*dvpL+dv*vpL + w*dwpL+dw*wpL 1088 & - ubar*dubarpL - dubar*ubarpL 1089! 1090 r53rR = u*durR+du*urR + v*dvrR+dv*vrR + w*dwrR+dw*wrR 1091 & - ubar*dubarrR - dubar*ubarrR 1092 r53uR = u*duuR+du*uuR + v*dvuR+dv*vuR + w*dwuR+dw*wuR 1093 & - ubar*dubaruR - dubar*ubaruR 1094 r53vR = u*duvR+du*uvR + v*dvvR+dv*vvR + w*dwvR+dw*wvR 1095 & - ubar*dubarvR - dubar*ubarvR 1096 r53wR = u*duwR+du*uwR + v*dvwR+dv*vwR + w*dwwR+dw*wwR 1097 & - ubar*dubarwR - dubar*ubarwR 1098 r53pR = u*dupR+du*upR + v*dvpR+dv*vpR + w*dwpR+dw*wpR 1099 & - ubar*dubarpR - dubar*ubarpR 1100! 1101 r24 = u 1102! 1103 r24rL = urL 1104 r24uL = uuL 1105 r24vL = uvL 1106 r24wL = uwL 1107 r24pL = upL 1108! 1109 r24rR = urR 1110 r24uR = uuR 1111 r24vR = uvR 1112 r24wR = uwR 1113 r24pR = upR 1114! 1115 r34 = v 1116! 1117 r34rL = vrL 1118 r34uL = vuL 1119 r34vL = vvL 1120 r34wL = vwL 1121 r34pL = vpL 1122! 1123 r34rR = vrR 1124 r34uR = vuR 1125 r34vR = vvR 1126 r34wR = vwR 1127 r34pR = vpR 1128 1129 r44 = w 1130! 1131 r44rL = wrL 1132 r44uL = wuL 1133 r44vL = wvL 1134 r44wL = wwL 1135 r44pL = wpL 1136! 1137 r44rR = wrR 1138 r44uR = wuR 1139 r44vR = wvR 1140 r44wR = wwR 1141 r44pR = wpR 1142! 1143 r54 = 0.5*q2 1144! 1145 r54rL = 0.5*q2rL 1146 r54uL = 0.5*q2uL 1147 r54vL = 0.5*q2vL 1148 r54wL = 0.5*q2wL 1149 r54pL = 0.5*q2pL 1150! 1151 r54rR = 0.5*q2rR 1152 r54uR = 0.5*q2uR 1153 r54vR = 0.5*q2vR 1154 r54wR = 0.5*q2wR 1155 r54pR = 0.5*q2pR 1156! 1157 t1 = eig1*dv1 + eig2*dv2 1158 & + eig3*dv4 1159! 1160 t1rL = eig1*dv1rL+dv1*eig1rL + eig2*dv2rL+dv2*eig2rL 1161 & + eig3*dv4rL+dv4*eig3rL 1162 t1uL = eig1*dv1uL+dv1*eig1uL + eig2*dv2uL+dv2*eig2uL 1163 & + eig3*dv4uL+dv4*eig3uL 1164 t1vL = eig1*dv1vL+dv1*eig1vL + eig2*dv2vL+dv2*eig2vL 1165 & + eig3*dv4vL+dv4*eig3vL 1166 t1wL = eig1*dv1wL+dv1*eig1wL + eig2*dv2wL+dv2*eig2wL 1167 & + eig3*dv4wL+dv4*eig3wL 1168 t1pL = eig1*dv1pL+dv1*eig1pL + eig2*dv2pL+dv2*eig2pL 1169 & + eig3*dv4pL+dv4*eig3pL 1170! 1171 t1rR = eig1*dv1rR+dv1*eig1rR + eig2*dv2rR+dv2*eig2rR 1172 & + eig3*dv4rR+dv4*eig3rR 1173 t1uR = eig1*dv1uR+dv1*eig1uR + eig2*dv2uR+dv2*eig2uR 1174 & + eig3*dv4uR+dv4*eig3uR 1175 t1vR = eig1*dv1vR+dv1*eig1vR + eig2*dv2vR+dv2*eig2vR 1176 & + eig3*dv4vR+dv4*eig3vR 1177 t1wR = eig1*dv1wR+dv1*eig1wR + eig2*dv2wR+dv2*eig2wR 1178 & + eig3*dv4wR+dv4*eig3wR 1179 t1pR = eig1*dv1pR+dv1*eig1pR + eig2*dv2pR+dv2*eig2pR 1180 & + eig3*dv4pR+dv4*eig3pR 1181! 1182 t2 = eig1*r21*dv1 + eig2*r22*dv2 1183 & + eig3*r23*dv3 + eig3*r24*dv4 1184! 1185 t2rL = eig1*(r21*dv1rL+dv1*r21rL)+r21*dv1*eig1rL 1186 & + eig2*(r22*dv2rL+dv2*r22rL)+r22*dv2*eig2rL 1187 & + eig3*(r23*dv3rL+dv3*r23rL)+r23*dv3*eig3rL 1188 & + eig3*(r24*dv4rL+dv4*r24rL)+r24*dv4*eig3rL 1189! 1190 t2uL = eig1*(r21*dv1uL+dv1*r21uL)+r21*dv1*eig1uL 1191 & + eig2*(r22*dv2uL+dv2*r22uL)+r22*dv2*eig2uL 1192 & + eig3*(r23*dv3uL+dv3*r23uL)+r23*dv3*eig3uL 1193 & + eig3*(r24*dv4uL+dv4*r24uL)+r24*dv4*eig3uL 1194! 1195 t2vL = eig1*(r21*dv1vL+dv1*r21vL)+r21*dv1*eig1vL 1196 & + eig2*(r22*dv2vL+dv2*r22vL)+r22*dv2*eig2vL 1197 & + eig3*(r23*dv3vL+dv3*r23vL)+r23*dv3*eig3vL 1198 & + eig3*(r24*dv4vL+dv4*r24vL)+r24*dv4*eig3vL 1199! 1200 t2wL = eig1*(r21*dv1wL+dv1*r21wL)+r21*dv1*eig1wL 1201 & + eig2*(r22*dv2wL+dv2*r22wL)+r22*dv2*eig2wL 1202 & + eig3*(r23*dv3wL+dv3*r23wL)+r23*dv3*eig3wL 1203 & + eig3*(r24*dv4wL+dv4*r24wL)+r24*dv4*eig3wL 1204! 1205 t2pL = eig1*(r21*dv1pL+dv1*r21pL)+r21*dv1*eig1pL 1206 & + eig2*(r22*dv2pL+dv2*r22pL)+r22*dv2*eig2pL 1207 & + eig3*(r23*dv3pL+dv3*r23pL)+r23*dv3*eig3pL 1208 & + eig3*(r24*dv4pL+dv4*r24pL)+r24*dv4*eig3pL 1209! 1210! 1211 t2rR = eig1*(r21*dv1rR+dv1*r21rR)+r21*dv1*eig1rR 1212 & + eig2*(r22*dv2rR+dv2*r22rR)+r22*dv2*eig2rR 1213 & + eig3*(r23*dv3rR+dv3*r23rR)+r23*dv3*eig3rR 1214 & + eig3*(r24*dv4rR+dv4*r24rR)+r24*dv4*eig3rR 1215! 1216 t2uR = eig1*(r21*dv1uR+dv1*r21uR)+r21*dv1*eig1uR 1217 & + eig2*(r22*dv2uR+dv2*r22uR)+r22*dv2*eig2uR 1218 & + eig3*(r23*dv3uR+dv3*r23uR)+r23*dv3*eig3uR 1219 & + eig3*(r24*dv4uR+dv4*r24uR)+r24*dv4*eig3uR 1220! 1221 t2vR = eig1*(r21*dv1vR+dv1*r21vR)+r21*dv1*eig1vR 1222 & + eig2*(r22*dv2vR+dv2*r22vR)+r22*dv2*eig2vR 1223 & + eig3*(r23*dv3vR+dv3*r23vR)+r23*dv3*eig3vR 1224 & + eig3*(r24*dv4vR+dv4*r24vR)+r24*dv4*eig3vR 1225! 1226 t2wR = eig1*(r21*dv1wR+dv1*r21wR)+r21*dv1*eig1wR 1227 & + eig2*(r22*dv2wR+dv2*r22wR)+r22*dv2*eig2wR 1228 & + eig3*(r23*dv3wR+dv3*r23wR)+r23*dv3*eig3wR 1229 & + eig3*(r24*dv4wR+dv4*r24wR)+r24*dv4*eig3wR 1230! 1231 t2pR = eig1*(r21*dv1pR+dv1*r21pR)+r21*dv1*eig1pR 1232 & + eig2*(r22*dv2pR+dv2*r22pR)+r22*dv2*eig2pR 1233 & + eig3*(r23*dv3pR+dv3*r23pR)+r23*dv3*eig3pR 1234 & + eig3*(r24*dv4pR+dv4*r24pR)+r24*dv4*eig3pR 1235! 1236! 1237 t3 = eig1*r31*dv1 + eig2*r32*dv2 1238 & + eig3*r33*dv3 + eig3*r34*dv4 1239! 1240 t3rL = eig1*(r31*dv1rL+dv1*r31rL)+r31*dv1*eig1rL 1241 & + eig2*(r32*dv2rL+dv2*r32rL)+r32*dv2*eig2rL 1242 & + eig3*(r33*dv3rL+dv3*r33rL)+r33*dv3*eig3rL 1243 & + eig3*(r34*dv4rL+dv4*r34rL)+r34*dv4*eig3rL 1244! 1245 t3uL = eig1*(r31*dv1uL+dv1*r31uL)+r31*dv1*eig1uL 1246 & + eig2*(r32*dv2uL+dv2*r32uL)+r32*dv2*eig2uL 1247 & + eig3*(r33*dv3uL+dv3*r33uL)+r33*dv3*eig3uL 1248 & + eig3*(r34*dv4uL+dv4*r34uL)+r34*dv4*eig3uL 1249! 1250 t3vL = eig1*(r31*dv1vL+dv1*r31vL)+r31*dv1*eig1vL 1251 & + eig2*(r32*dv2vL+dv2*r32vL)+r32*dv2*eig2vL 1252 & + eig3*(r33*dv3vL+dv3*r33vL)+r33*dv3*eig3vL 1253 & + eig3*(r34*dv4vL+dv4*r34vL)+r34*dv4*eig3vL 1254! 1255 t3wL = eig1*(r31*dv1wL+dv1*r31wL)+r31*dv1*eig1wL 1256 & + eig2*(r32*dv2wL+dv2*r32wL)+r32*dv2*eig2wL 1257 & + eig3*(r33*dv3wL+dv3*r33wL)+r33*dv3*eig3wL 1258 & + eig3*(r34*dv4wL+dv4*r34wL)+r34*dv4*eig3wL 1259! 1260 t3pL = eig1*(r31*dv1pL+dv1*r31pL)+r31*dv1*eig1pL 1261 & + eig2*(r32*dv2pL+dv2*r32pL)+r32*dv2*eig2pL 1262 & + eig3*(r33*dv3pL+dv3*r33pL)+r33*dv3*eig3pL 1263 & + eig3*(r34*dv4pL+dv4*r34pL)+r34*dv4*eig3pL 1264! 1265! 1266 t3rR = eig1*(r31*dv1rR+dv1*r31rR)+r31*dv1*eig1rR 1267 & + eig2*(r32*dv2rR+dv2*r32rR)+r32*dv2*eig2rR 1268 & + eig3*(r33*dv3rR+dv3*r33rR)+r33*dv3*eig3rR 1269 & + eig3*(r34*dv4rR+dv4*r34rR)+r34*dv4*eig3rR 1270! 1271 t3uR = eig1*(r31*dv1uR+dv1*r31uR)+r31*dv1*eig1uR 1272 & + eig2*(r32*dv2uR+dv2*r32uR)+r32*dv2*eig2uR 1273 & + eig3*(r33*dv3uR+dv3*r33uR)+r33*dv3*eig3uR 1274 & + eig3*(r34*dv4uR+dv4*r34uR)+r34*dv4*eig3uR 1275! 1276 t3vR = eig1*(r31*dv1vR+dv1*r31vR)+r31*dv1*eig1vR 1277 & + eig2*(r32*dv2vR+dv2*r32vR)+r32*dv2*eig2vR 1278 & + eig3*(r33*dv3vR+dv3*r33vR)+r33*dv3*eig3vR 1279 & + eig3*(r34*dv4vR+dv4*r34vR)+r34*dv4*eig3vR 1280! 1281 t3wR = eig1*(r31*dv1wR+dv1*r31wR)+r31*dv1*eig1wR 1282 & + eig2*(r32*dv2wR+dv2*r32wR)+r32*dv2*eig2wR 1283 & + eig3*(r33*dv3wR+dv3*r33wR)+r33*dv3*eig3wR 1284 & + eig3*(r34*dv4wR+dv4*r34wR)+r34*dv4*eig3wR 1285! 1286 t3pR = eig1*(r31*dv1pR+dv1*r31pR)+r31*dv1*eig1pR 1287 & + eig2*(r32*dv2pR+dv2*r32pR)+r32*dv2*eig2pR 1288 & + eig3*(r33*dv3pR+dv3*r33pR)+r33*dv3*eig3pR 1289 & + eig3*(r34*dv4pR+dv4*r34pR)+r34*dv4*eig3pR 1290! 1291! 1292 t4 = eig1*r41*dv1 + eig2*r42*dv2 1293 & + eig3*r43*dv3 + eig3*r44*dv4 1294! 1295 t4rL = eig1*(r41*dv1rL+dv1*r41rL)+r41*dv1*eig1rL 1296 & + eig2*(r42*dv2rL+dv2*r42rL)+r42*dv2*eig2rL 1297 & + eig3*(r43*dv3rL+dv3*r43rL)+r43*dv3*eig3rL 1298 & + eig3*(r44*dv4rL+dv4*r44rL)+r44*dv4*eig3rL 1299! 1300 t4uL = eig1*(r41*dv1uL+dv1*r41uL)+r41*dv1*eig1uL 1301 & + eig2*(r42*dv2uL+dv2*r42uL)+r42*dv2*eig2uL 1302 & + eig3*(r43*dv3uL+dv3*r43uL)+r43*dv3*eig3uL 1303 & + eig3*(r44*dv4uL+dv4*r44uL)+r44*dv4*eig3uL 1304! 1305 t4vL = eig1*(r41*dv1vL+dv1*r41vL)+r41*dv1*eig1vL 1306 & + eig2*(r42*dv2vL+dv2*r42vL)+r42*dv2*eig2vL 1307 & + eig3*(r43*dv3vL+dv3*r43vL)+r43*dv3*eig3vL 1308 & + eig3*(r44*dv4vL+dv4*r44vL)+r44*dv4*eig3vL 1309! 1310 t4wL = eig1*(r41*dv1wL+dv1*r41wL)+r41*dv1*eig1wL 1311 & + eig2*(r42*dv2wL+dv2*r42wL)+r42*dv2*eig2wL 1312 & + eig3*(r43*dv3wL+dv3*r43wL)+r43*dv3*eig3wL 1313 & + eig3*(r44*dv4wL+dv4*r44wL)+r44*dv4*eig3wL 1314! 1315 t4pL = eig1*(r41*dv1pL+dv1*r41pL)+r41*dv1*eig1pL 1316 & + eig2*(r42*dv2pL+dv2*r42pL)+r42*dv2*eig2pL 1317 & + eig3*(r43*dv3pL+dv3*r43pL)+r43*dv3*eig3pL 1318 & + eig3*(r44*dv4pL+dv4*r44pL)+r44*dv4*eig3pL 1319! 1320! 1321 t4rR = eig1*(r41*dv1rR+dv1*r41rR)+r41*dv1*eig1rR 1322 & + eig2*(r42*dv2rR+dv2*r42rR)+r42*dv2*eig2rR 1323 & + eig3*(r43*dv3rR+dv3*r43rR)+r43*dv3*eig3rR 1324 & + eig3*(r44*dv4rR+dv4*r44rR)+r44*dv4*eig3rR 1325! 1326 t4uR = eig1*(r41*dv1uR+dv1*r41uR)+r41*dv1*eig1uR 1327 & + eig2*(r42*dv2uR+dv2*r42uR)+r42*dv2*eig2uR 1328 & + eig3*(r43*dv3uR+dv3*r43uR)+r43*dv3*eig3uR 1329 & + eig3*(r44*dv4uR+dv4*r44uR)+r44*dv4*eig3uR 1330! 1331 t4vR = eig1*(r41*dv1vR+dv1*r41vR)+r41*dv1*eig1vR 1332 & + eig2*(r42*dv2vR+dv2*r42vR)+r42*dv2*eig2vR 1333 & + eig3*(r43*dv3vR+dv3*r43vR)+r43*dv3*eig3vR 1334 & + eig3*(r44*dv4vR+dv4*r44vR)+r44*dv4*eig3vR 1335! 1336 t4wR = eig1*(r41*dv1wR+dv1*r41wR)+r41*dv1*eig1wR 1337 & + eig2*(r42*dv2wR+dv2*r42wR)+r42*dv2*eig2wR 1338 & + eig3*(r43*dv3wR+dv3*r43wR)+r43*dv3*eig3wR 1339 & + eig3*(r44*dv4wR+dv4*r44wR)+r44*dv4*eig3wR 1340! 1341 t4pR = eig1*(r41*dv1pR+dv1*r41pR)+r41*dv1*eig1pR 1342 & + eig2*(r42*dv2pR+dv2*r42pR)+r42*dv2*eig2pR 1343 & + eig3*(r43*dv3pR+dv3*r43pR)+r43*dv3*eig3pR 1344 & + eig3*(r44*dv4pR+dv4*r44pR)+r44*dv4*eig3pR 1345! 1346! 1347 t5 = eig1*r51*dv1 + eig2*r52*dv2 1348 & + eig3*r53*dv3 + eig3*r54*dv4 1349! 1350 t5rL = eig1*(r51*dv1rL+dv1*r51rL)+r51*dv1*eig1rL 1351 & + eig2*(r52*dv2rL+dv2*r52rL)+r52*dv2*eig2rL 1352 & + eig3*(r53*dv3rL+dv3*r53rL)+r53*dv3*eig3rL 1353 & + eig3*(r54*dv4rL+dv4*r54rL)+r54*dv4*eig3rL 1354! 1355 t5uL = eig1*(r51*dv1uL+dv1*r51uL)+r51*dv1*eig1uL 1356 & + eig2*(r52*dv2uL+dv2*r52uL)+r52*dv2*eig2uL 1357 & + eig3*(r53*dv3uL+dv3*r53uL)+r53*dv3*eig3uL 1358 & + eig3*(r54*dv4uL+dv4*r54uL)+r54*dv4*eig3uL 1359! 1360 t5vL = eig1*(r51*dv1vL+dv1*r51vL)+r51*dv1*eig1vL 1361 & + eig2*(r52*dv2vL+dv2*r52vL)+r52*dv2*eig2vL 1362 & + eig3*(r53*dv3vL+dv3*r53vL)+r53*dv3*eig3vL 1363 & + eig3*(r54*dv4vL+dv4*r54vL)+r54*dv4*eig3vL 1364! 1365 t5wL = eig1*(r51*dv1wL+dv1*r51wL)+r51*dv1*eig1wL 1366 & + eig2*(r52*dv2wL+dv2*r52wL)+r52*dv2*eig2wL 1367 & + eig3*(r53*dv3wL+dv3*r53wL)+r53*dv3*eig3wL 1368 & + eig3*(r54*dv4wL+dv4*r54wL)+r54*dv4*eig3wL 1369! 1370 t5pL = eig1*(r51*dv1pL+dv1*r51pL)+r51*dv1*eig1pL 1371 & + eig2*(r52*dv2pL+dv2*r52pL)+r52*dv2*eig2pL 1372 & + eig3*(r53*dv3pL+dv3*r53pL)+r53*dv3*eig3pL 1373 & + eig3*(r54*dv4pL+dv4*r54pL)+r54*dv4*eig3pL 1374! 1375! 1376 t5rR = eig1*(r51*dv1rR+dv1*r51rR)+r51*dv1*eig1rR 1377 & + eig2*(r52*dv2rR+dv2*r52rR)+r52*dv2*eig2rR 1378 & + eig3*(r53*dv3rR+dv3*r53rR)+r53*dv3*eig3rR 1379 & + eig3*(r54*dv4rR+dv4*r54rR)+r54*dv4*eig3rR 1380! 1381 t5uR = eig1*(r51*dv1uR+dv1*r51uR)+r51*dv1*eig1uR 1382 & + eig2*(r52*dv2uR+dv2*r52uR)+r52*dv2*eig2uR 1383 & + eig3*(r53*dv3uR+dv3*r53uR)+r53*dv3*eig3uR 1384 & + eig3*(r54*dv4uR+dv4*r54uR)+r54*dv4*eig3uR 1385! 1386 t5vR = eig1*(r51*dv1vR+dv1*r51vR)+r51*dv1*eig1vR 1387 & + eig2*(r52*dv2vR+dv2*r52vR)+r52*dv2*eig2vR 1388 & + eig3*(r53*dv3vR+dv3*r53vR)+r53*dv3*eig3vR 1389 & + eig3*(r54*dv4vR+dv4*r54vR)+r54*dv4*eig3vR 1390! 1391 t5wR = eig1*(r51*dv1wR+dv1*r51wR)+r51*dv1*eig1wR 1392 & + eig2*(r52*dv2wR+dv2*r52wR)+r52*dv2*eig2wR 1393 & + eig3*(r53*dv3wR+dv3*r53wR)+r53*dv3*eig3wR 1394 & + eig3*(r54*dv4wR+dv4*r54wR)+r54*dv4*eig3wR 1395! 1396 t5pR = eig1*(r51*dv1pR+dv1*r51pR)+r51*dv1*eig1pR 1397 & + eig2*(r52*dv2pR+dv2*r52pR)+r52*dv2*eig2pR 1398 & + eig3*(r53*dv3pR+dv3*r53pR)+r53*dv3*eig3pR 1399 & + eig3*(r54*dv4pR+dv4*r54pR)+r54*dv4*eig3pR 1400! 1401! Compute flux using variables from left side of face 1402! 1403 fluxp1 = area*rhol*ubarl 1404! 1405 fluxp1rL = area*(rhol*ubarlrL + ubarl*rholrL) 1406 fluxp1uL = area*(rhol*ubarluL + ubarl*rholuL) 1407 fluxp1vL = area*(rhol*ubarlvL + ubarl*rholvL) 1408 fluxp1wL = area*(rhol*ubarlwL + ubarl*rholwL) 1409 fluxp1pL = area*(rhol*ubarlpL + ubarl*rholpL) 1410! 1411 fluxp1rR = area*(rhol*ubarlrR + ubarl*rholrR) 1412 fluxp1uR = area*(rhol*ubarluR + ubarl*rholuR) 1413 fluxp1vR = area*(rhol*ubarlvR + ubarl*rholvR) 1414 fluxp1wR = area*(rhol*ubarlwR + ubarl*rholwR) 1415 fluxp1pR = area*(rhol*ubarlpR + ubarl*rholpR) 1416! 1417 fluxp2 = area*(rhol*ul*ubarl + xnorm*pressl) 1418! 1419 fluxp2rL = area*(rhol*(ul*ubarlrL+ubarl*ulrL) + 1420 & ul*ubarl*rholrL + xnorm*presslrL) 1421 fluxp2uL = area*(rhol*(ul*ubarluL+ubarl*uluL) + 1422 & ul*ubarl*rholuL + xnorm*pressluL) 1423 fluxp2vL = area*(rhol*(ul*ubarlvL+ubarl*ulvL) + 1424 & ul*ubarl*rholvL + xnorm*presslvL) 1425 fluxp2wL = area*(rhol*(ul*ubarlwL+ubarl*ulwL) + 1426 & ul*ubarl*rholwL + xnorm*presslwL) 1427 fluxp2pL = area*(rhol*(ul*ubarlpL+ubarl*ulpL) + 1428 & ul*ubarl*rholpL + xnorm*presslpL) 1429! 1430 fluxp2rR = area*(rhol*(ul*ubarlrR+ubarl*ulrR) + 1431 & ul*ubarl*rholrR + xnorm*presslrR) 1432 fluxp2uR = area*(rhol*(ul*ubarluR+ubarl*uluR) + 1433 & ul*ubarl*rholuR + xnorm*pressluR) 1434 fluxp2vR = area*(rhol*(ul*ubarlvR+ubarl*ulvR) + 1435 & ul*ubarl*rholvR + xnorm*presslvR) 1436 fluxp2wR = area*(rhol*(ul*ubarlwR+ubarl*ulwR) + 1437 & ul*ubarl*rholwR + xnorm*presslwR) 1438 fluxp2pR = area*(rhol*(ul*ubarlpR+ubarl*ulpR) + 1439 & ul*ubarl*rholpR + xnorm*presslpR) 1440! 1441! 1442 fluxp3 = area*(rhol*vl*ubarl + ynorm*pressl) 1443! 1444 fluxp3rL = area*(rhol*(vl*ubarlrL+ubarl*vlrL) + 1445 & vl*ubarl*rholrL + ynorm*presslrL) 1446 fluxp3uL = area*(rhol*(vl*ubarluL+ubarl*vluL) + 1447 & vl*ubarl*rholuL + ynorm*pressluL) 1448 fluxp3vL = area*(rhol*(vl*ubarlvL+ubarl*vlvL) + 1449 & vl*ubarl*rholvL + ynorm*presslvL) 1450 fluxp3wL = area*(rhol*(vl*ubarlwL+ubarl*vlwL) + 1451 & vl*ubarl*rholwL + ynorm*presslwL) 1452 fluxp3pL = area*(rhol*(vl*ubarlpL+ubarl*vlpL) + 1453 & vl*ubarl*rholpL + ynorm*presslpL) 1454! 1455 fluxp3rR = area*(rhol*(vl*ubarlrR+ubarl*vlrR) + 1456 & vl*ubarl*rholrR + ynorm*presslrR) 1457 fluxp3uR = area*(rhol*(vl*ubarluR+ubarl*vluR) + 1458 & vl*ubarl*rholuR + ynorm*pressluR) 1459 fluxp3vR = area*(rhol*(vl*ubarlvR+ubarl*vlvR) + 1460 & vl*ubarl*rholvR + ynorm*presslvR) 1461 fluxp3wR = area*(rhol*(vl*ubarlwR+ubarl*vlwR) + 1462 & vl*ubarl*rholwR + ynorm*presslwR) 1463 fluxp3pR = area*(rhol*(vl*ubarlpR+ubarl*vlpR) + 1464 & vl*ubarl*rholpR + ynorm*presslpR) 1465! 1466! 1467 fluxp4 = area*(rhol*wl*ubarl + znorm*pressl) 1468! 1469 fluxp4rL = area*(rhol*(wl*ubarlrL+ubarl*wlrL) + 1470 & wl*ubarl*rholrL + znorm*presslrL) 1471 fluxp4uL = area*(rhol*(wl*ubarluL+ubarl*wluL) + 1472 & wl*ubarl*rholuL + znorm*pressluL) 1473 fluxp4vL = area*(rhol*(wl*ubarlvL+ubarl*wlvL) + 1474 & wl*ubarl*rholvL + znorm*presslvL) 1475 fluxp4wL = area*(rhol*(wl*ubarlwL+ubarl*wlwL) + 1476 & wl*ubarl*rholwL + znorm*presslwL) 1477 fluxp4pL = area*(rhol*(wl*ubarlpL+ubarl*wlpL) + 1478 & wl*ubarl*rholpL + znorm*presslpL) 1479! 1480 fluxp4rR = area*(rhol*(wl*ubarlrR+ubarl*wlrR) + 1481 & wl*ubarl*rholrR + znorm*presslrR) 1482 fluxp4uR = area*(rhol*(wl*ubarluR+ubarl*wluR) + 1483 & wl*ubarl*rholuR + znorm*pressluR) 1484 fluxp4vR = area*(rhol*(wl*ubarlvR+ubarl*wlvR) + 1485 & wl*ubarl*rholvR + znorm*presslvR) 1486 fluxp4wR = area*(rhol*(wl*ubarlwR+ubarl*wlwR) + 1487 & wl*ubarl*rholwR + znorm*presslwR) 1488 fluxp4pR = area*(rhol*(wl*ubarlpR+ubarl*wlpR) + 1489 & wl*ubarl*rholpR + znorm*presslpR) 1490! 1491! 1492 fluxp5 = area*(enrgyl + pressl)*ubarl 1493! 1494 fluxp5rL = area*((enrgyl + pressl)*ubarlrL + 1495 & ubarl*(enrgylrL + presslrL)) 1496 fluxp5uL = area*((enrgyl + pressl)*ubarluL + 1497 & ubarl*(enrgyluL + pressluL)) 1498 fluxp5vL = area*((enrgyl + pressl)*ubarlvL + 1499 & ubarl*(enrgylvL + presslvL)) 1500 fluxp5wL = area*((enrgyl + pressl)*ubarlwL + 1501 & ubarl*(enrgylwL + presslwL)) 1502 fluxp5pL = area*((enrgyl + pressl)*ubarlpL + 1503 & ubarl*(enrgylpL + presslpL)) 1504! 1505 fluxp5rR = area*((enrgyl + pressl)*ubarlrR + 1506 & ubarl*(enrgylrR + presslrR)) 1507 fluxp5uR = area*((enrgyl + pressl)*ubarluR + 1508 & ubarl*(enrgyluR + pressluR)) 1509 fluxp5vR = area*((enrgyl + pressl)*ubarlvR + 1510 & ubarl*(enrgylvR + presslvR)) 1511 fluxp5wR = area*((enrgyl + pressl)*ubarlwR + 1512 & ubarl*(enrgylwR + presslwR)) 1513 fluxp5pR = area*((enrgyl + pressl)*ubarlpR + 1514 & ubarl*(enrgylpR + presslpR)) 1515! 1516! Now the right side 1517! 1518 fluxm1 = area*rhor*ubarr 1519! 1520 fluxm1rL = area*(rhor*ubarrrL + ubarr*rhorrL) 1521 fluxm1uL = area*(rhor*ubarruL + ubarr*rhoruL) 1522 fluxm1vL = area*(rhor*ubarrvL + ubarr*rhorvL) 1523 fluxm1wL = area*(rhor*ubarrwL + ubarr*rhorwL) 1524 fluxm1pL = area*(rhor*ubarrpL + ubarr*rhorpL) 1525! 1526 fluxm1rR = area*(rhor*ubarrrR + ubarr*rhorrR) 1527 fluxm1uR = area*(rhor*ubarruR + ubarr*rhoruR) 1528 fluxm1vR = area*(rhor*ubarrvR + ubarr*rhorvR) 1529 fluxm1wR = area*(rhor*ubarrwR + ubarr*rhorwR) 1530 fluxm1pR = area*(rhor*ubarrpR + ubarr*rhorpR) 1531! 1532! 1533 fluxm2 = area*(rhor*ur*ubarr + xnorm*pressr) 1534! 1535 fluxm2rL = area*(rhor*(ur*ubarrrL+ubarr*urrL) + 1536 & ur*ubarr*rhorrL + xnorm*pressrrL) 1537 fluxm2uL = area*(rhor*(ur*ubarruL+ubarr*uruL) + 1538 & ur*ubarr*rhoruL + xnorm*pressruL) 1539 fluxm2vL = area*(rhor*(ur*ubarrvL+ubarr*urvL) + 1540 & ur*ubarr*rhorvL + xnorm*pressrvL) 1541 fluxm2wL = area*(rhor*(ur*ubarrwL+ubarr*urwL) + 1542 & ur*ubarr*rhorwL + xnorm*pressrwL) 1543 fluxm2pL = area*(rhor*(ur*ubarrpL+ubarr*urpL) + 1544 & ur*ubarr*rhorpL + xnorm*pressrpL) 1545! 1546 fluxm2rR = area*(rhor*(ur*ubarrrR+ubarr*urrR) + 1547 & ur*ubarr*rhorrR + xnorm*pressrrR) 1548 fluxm2uR = area*(rhor*(ur*ubarruR+ubarr*uruR) + 1549 & ur*ubarr*rhoruR + xnorm*pressruR) 1550 fluxm2vR = area*(rhor*(ur*ubarrvR+ubarr*urvR) + 1551 & ur*ubarr*rhorvR + xnorm*pressrvR) 1552 fluxm2wR = area*(rhor*(ur*ubarrwR+ubarr*urwR) + 1553 & ur*ubarr*rhorwR + xnorm*pressrwR) 1554 fluxm2pR = area*(rhor*(ur*ubarrpR+ubarr*urpR) + 1555 & ur*ubarr*rhorpR + xnorm*pressrpR) 1556! 1557! 1558 fluxm3 = area*(rhor*vr*ubarr + ynorm*pressr) 1559! 1560 fluxm3rL = area*(rhor*(vr*ubarrrL+ubarr*vrrL) + 1561 & vr*ubarr*rhorrL + ynorm*pressrrL) 1562 fluxm3uL = area*(rhor*(vr*ubarruL+ubarr*vruL) + 1563 & vr*ubarr*rhoruL + ynorm*pressruL) 1564 fluxm3vL = area*(rhor*(vr*ubarrvL+ubarr*vrvL) + 1565 & vr*ubarr*rhorvL + ynorm*pressrvL) 1566 fluxm3wL = area*(rhor*(vr*ubarrwL+ubarr*vrwL) + 1567 & vr*ubarr*rhorwL + ynorm*pressrwL) 1568 fluxm3pL = area*(rhor*(vr*ubarrpL+ubarr*vrpL) + 1569 & vr*ubarr*rhorpL + ynorm*pressrpL) 1570! 1571 fluxm3rR = area*(rhor*(vr*ubarrrR+ubarr*vrrR) + 1572 & vr*ubarr*rhorrR + ynorm*pressrrR) 1573 fluxm3uR = area*(rhor*(vr*ubarruR+ubarr*vruR) + 1574 & vr*ubarr*rhoruR + ynorm*pressruR) 1575 fluxm3vR = area*(rhor*(vr*ubarrvR+ubarr*vrvR) + 1576 & vr*ubarr*rhorvR + ynorm*pressrvR) 1577 fluxm3wR = area*(rhor*(vr*ubarrwR+ubarr*vrwR) + 1578 & vr*ubarr*rhorwR + ynorm*pressrwR) 1579 fluxm3pR = area*(rhor*(vr*ubarrpR+ubarr*vrpR) + 1580 & vr*ubarr*rhorpR + ynorm*pressrpR) 1581! 1582! 1583 fluxm4 = area*(rhor*wr*ubarr + znorm*pressr) 1584! 1585 fluxm4rL = area*(rhor*(wr*ubarrrL+ubarr*wrrL) + 1586 & wr*ubarr*rhorrL + znorm*pressrrL) 1587 fluxm4uL = area*(rhor*(wr*ubarruL+ubarr*wruL) + 1588 & wr*ubarr*rhoruL + znorm*pressruL) 1589 fluxm4vL = area*(rhor*(wr*ubarrvL+ubarr*wrvL) + 1590 & wr*ubarr*rhorvL + znorm*pressrvL) 1591 fluxm4wL = area*(rhor*(wr*ubarrwL+ubarr*wrwL) + 1592 & wr*ubarr*rhorwL + znorm*pressrwL) 1593 fluxm4pL = area*(rhor*(wr*ubarrpL+ubarr*wrpL) + 1594 & wr*ubarr*rhorpL + znorm*pressrpL) 1595! 1596 fluxm4rR = area*(rhor*(wr*ubarrrR+ubarr*wrrR) + 1597 & wr*ubarr*rhorrR + znorm*pressrrR) 1598 fluxm4uR = area*(rhor*(wr*ubarruR+ubarr*wruR) + 1599 & wr*ubarr*rhoruR + znorm*pressruR) 1600 fluxm4vR = area*(rhor*(wr*ubarrvR+ubarr*wrvR) + 1601 & wr*ubarr*rhorvR + znorm*pressrvR) 1602 fluxm4wR = area*(rhor*(wr*ubarrwR+ubarr*wrwR) + 1603 & wr*ubarr*rhorwR + znorm*pressrwR) 1604 fluxm4pR = area*(rhor*(wr*ubarrpR+ubarr*wrpR) + 1605 & wr*ubarr*rhorpR + znorm*pressrpR) 1606! 1607 fluxm5 = area*(enrgyr + pressr)*ubarr 1608! 1609 fluxm5rL = area*((enrgyr + pressr)*ubarrrL + 1610 & ubarr*(enrgyrrL + pressrrL)) 1611 fluxm5uL = area*((enrgyr + pressr)*ubarruL + 1612 & ubarr*(enrgyruL + pressruL)) 1613 fluxm5vL = area*((enrgyr + pressr)*ubarrvL + 1614 & ubarr*(enrgyrvL + pressrvL)) 1615 fluxm5wL = area*((enrgyr + pressr)*ubarrwL + 1616 & ubarr*(enrgyrwL + pressrwL)) 1617 fluxm5pL = area*((enrgyr + pressr)*ubarrpL + 1618 & ubarr*(enrgyrpL + pressrpL)) 1619! 1620 fluxm5rR = area*((enrgyr + pressr)*ubarrrR + 1621 & ubarr*(enrgyrrR + pressrrR)) 1622 fluxm5uR = area*((enrgyr + pressr)*ubarruR + 1623 & ubarr*(enrgyruR + pressruR)) 1624 fluxm5vR = area*((enrgyr + pressr)*ubarrvR + 1625 & ubarr*(enrgyrvR + pressrvR)) 1626 fluxm5wR = area*((enrgyr + pressr)*ubarrwR + 1627 & ubarr*(enrgyrwR + pressrwR)) 1628 fluxm5pR = area*((enrgyr + pressr)*ubarrpR + 1629 & ubarr*(enrgyrpR + pressrpR)) 1630! 1631 flux1 = 0.5*(fluxp1 + fluxm1 - area*t1) 1632 flux2 = 0.5*(fluxp2 + fluxm2 - area*t2) 1633 flux3 = 0.5*(fluxp3 + fluxm3 - area*t3) 1634 flux4 = 0.5*(fluxp4 + fluxm4 - area*t4) 1635 flux5 = 0.5*(fluxp5 + fluxm5 - area*t5) 1636! 1637 flux1rL = 0.5*(fluxp1rL + fluxm1rL - area*t1rL) 1638 flux1uL = 0.5*(fluxp1uL + fluxm1uL - area*t1uL) 1639 flux1vL = 0.5*(fluxp1vL + fluxm1vL - area*t1vL) 1640 flux1wL = 0.5*(fluxp1wL + fluxm1wL - area*t1wL) 1641 flux1pL = 0.5*(fluxp1pL + fluxm1pL - area*t1pL) 1642 flux1rR = 0.5*(fluxp1rR + fluxm1rR - area*t1rR) 1643 flux1uR = 0.5*(fluxp1uR + fluxm1uR - area*t1uR) 1644 flux1vR = 0.5*(fluxp1vR + fluxm1vR - area*t1vR) 1645 flux1wR = 0.5*(fluxp1wR + fluxm1wR - area*t1wR) 1646 flux1pR = 0.5*(fluxp1pR + fluxm1pR - area*t1pR) 1647! 1648 flux2rL = 0.5*(fluxp2rL + fluxm2rL - area*t2rL) 1649 flux2uL = 0.5*(fluxp2uL + fluxm2uL - area*t2uL) 1650 flux2vL = 0.5*(fluxp2vL + fluxm2vL - area*t2vL) 1651 flux2wL = 0.5*(fluxp2wL + fluxm2wL - area*t2wL) 1652 flux2pL = 0.5*(fluxp2pL + fluxm2pL - area*t2pL) 1653 flux2rR = 0.5*(fluxp2rR + fluxm2rR - area*t2rR) 1654 flux2uR = 0.5*(fluxp2uR + fluxm2uR - area*t2uR) 1655 flux2vR = 0.5*(fluxp2vR + fluxm2vR - area*t2vR) 1656 flux2wR = 0.5*(fluxp2wR + fluxm2wR - area*t2wR) 1657 flux2pR = 0.5*(fluxp2pR + fluxm2pR - area*t2pR) 1658! 1659 flux3rL = 0.5*(fluxp3rL + fluxm3rL - area*t3rL) 1660 flux3uL = 0.5*(fluxp3uL + fluxm3uL - area*t3uL) 1661 flux3vL = 0.5*(fluxp3vL + fluxm3vL - area*t3vL) 1662 flux3wL = 0.5*(fluxp3wL + fluxm3wL - area*t3wL) 1663 flux3pL = 0.5*(fluxp3pL + fluxm3pL - area*t3pL) 1664 flux3rR = 0.5*(fluxp3rR + fluxm3rR - area*t3rR) 1665 flux3uR = 0.5*(fluxp3uR + fluxm3uR - area*t3uR) 1666 flux3vR = 0.5*(fluxp3vR + fluxm3vR - area*t3vR) 1667 flux3wR = 0.5*(fluxp3wR + fluxm3wR - area*t3wR) 1668 flux3pR = 0.5*(fluxp3pR + fluxm3pR - area*t3pR) 1669! 1670 flux4rL = 0.5*(fluxp4rL + fluxm4rL - area*t4rL) 1671 flux4uL = 0.5*(fluxp4uL + fluxm4uL - area*t4uL) 1672 flux4vL = 0.5*(fluxp4vL + fluxm4vL - area*t4vL) 1673 flux4wL = 0.5*(fluxp4wL + fluxm4wL - area*t4wL) 1674 flux4pL = 0.5*(fluxp4pL + fluxm4pL - area*t4pL) 1675 flux4rR = 0.5*(fluxp4rR + fluxm4rR - area*t4rR) 1676 flux4uR = 0.5*(fluxp4uR + fluxm4uR - area*t4uR) 1677 flux4vR = 0.5*(fluxp4vR + fluxm4vR - area*t4vR) 1678 flux4wR = 0.5*(fluxp4wR + fluxm4wR - area*t4wR) 1679 flux4pR = 0.5*(fluxp4pR + fluxm4pR - area*t4pR) 1680! 1681 flux5rL = 0.5*(fluxp5rL + fluxm5rL - area*t5rL) 1682 flux5uL = 0.5*(fluxp5uL + fluxm5uL - area*t5uL) 1683 flux5vL = 0.5*(fluxp5vL + fluxm5vL - area*t5vL) 1684 flux5wL = 0.5*(fluxp5wL + fluxm5wL - area*t5wL) 1685 flux5pL = 0.5*(fluxp5pL + fluxm5pL - area*t5pL) 1686 flux5rR = 0.5*(fluxp5rR + fluxm5rR - area*t5rR) 1687 flux5uR = 0.5*(fluxp5uR + fluxm5uR - area*t5uR) 1688 flux5vR = 0.5*(fluxp5vR + fluxm5vR - area*t5vR) 1689 flux5wR = 0.5*(fluxp5wR + fluxm5wR - area*t5wR) 1690 flux5pR = 0.5*(fluxp5pR + fluxm5pR - area*t5pR) 1691!/* 1692! These last 50 equations are the dfp's and dfm's we're 1693! looking for, but they're w.r.t. primitive variables. 1694! 1695! We need to convert these last results into conservative 1696! variables, using one last chain-rule. 1697! 1698! q = primitive 1699! Q = conservative 1700! 1701! First the left transformation 1702! */ 1703 q1Q1L = 1. 1704 q1Q2L = 0. 1705 q1Q3L = 0. 1706 q1Q4L = 0. 1707 q1Q5L = 0. 1708! 1709 q2Q1L = - ul / rhol 1710 q2Q2L = 1. / rhol 1711 q2Q3L = 0. 1712 q2Q4L = 0. 1713 q2Q5L = 0. 1714! 1715 q3Q1L = - vl / rhol 1716 q3Q2L = 0. 1717 q3Q3L = 1. / rhol 1718 q3Q4L = 0. 1719 q3Q5L = 0. 1720! 1721 q4Q1L = - wl / rhol 1722 q4Q2L = 0. 1723 q4Q3L = 0. 1724 q4Q4L = 1. / rhol 1725 q4Q5L = 0. 1726! 1727 q5Q1L = gm1 / 2. * (ul*ul + vl*vl + wl*wl) 1728 q5Q2L = - gm1 * ul 1729 q5Q3L = - gm1 * vl 1730 q5Q4L = - gm1 * wl 1731 q5Q5L = gm1 1732! 1733 dfp(1,1) = flux1rL*q1Q1L + flux1uL*q2Q1L + flux1vL*q3Q1L 1734 & + flux1wL*q4Q1L + flux1pL*q5Q1L 1735! 1736 dfp(1,2) = flux1rL*q1Q2L + flux1uL*q2Q2L + flux1vL*q3Q2L 1737 & + flux1wL*q4Q2L + flux1pL*q5Q2L 1738! 1739 dfp(1,3) = flux1rL*q1Q3L + flux1uL*q2Q3L + flux1vL*q3Q3L 1740 & + flux1wL*q4Q3L + flux1pL*q5Q3L 1741! 1742 dfp(1,4) = flux1rL*q1Q4L + flux1uL*q2Q4L + flux1vL*q3Q4L 1743 & + flux1wL*q4Q4L + flux1pL*q5Q4L 1744! 1745 dfp(1,5) = flux1rL*q1Q5L + flux1uL*q2Q5L + flux1vL*q3Q5L 1746 & + flux1wL*q4Q5L + flux1pL*q5Q5L 1747! 1748! 1749 dfp(2,1) = flux2rL*q1Q1L + flux2uL*q2Q1L + flux2vL*q3Q1L 1750 & + flux2wL*q4Q1L + flux2pL*q5Q1L 1751! 1752 dfp(2,2) = flux2rL*q1Q2L + flux2uL*q2Q2L + flux2vL*q3Q2L 1753 & + flux2wL*q4Q2L + flux2pL*q5Q2L 1754! 1755 dfp(2,3) = flux2rL*q1Q3L + flux2uL*q2Q3L + flux2vL*q3Q3L 1756 & + flux2wL*q4Q3L + flux2pL*q5Q3L 1757! 1758 dfp(2,4) = flux2rL*q1Q4L + flux2uL*q2Q4L + flux2vL*q3Q4L 1759 & + flux2wL*q4Q4L + flux2pL*q5Q4L 1760! 1761 dfp(2,5) = flux2rL*q1Q5L + flux2uL*q2Q5L + flux2vL*q3Q5L 1762 & + flux2wL*q4Q5L + flux2pL*q5Q5L 1763! 1764! 1765 dfp(3,1) = flux3rL*q1Q1L + flux3uL*q2Q1L + flux3vL*q3Q1L 1766 & + flux3wL*q4Q1L + flux3pL*q5Q1L 1767! 1768 dfp(3,2) = flux3rL*q1Q2L + flux3uL*q2Q2L + flux3vL*q3Q2L 1769 & + flux3wL*q4Q2L + flux3pL*q5Q2L 1770! 1771 dfp(3,3) = flux3rL*q1Q3L + flux3uL*q2Q3L + flux3vL*q3Q3L 1772 & + flux3wL*q4Q3L + flux3pL*q5Q3L 1773! 1774 dfp(3,4) = flux3rL*q1Q4L + flux3uL*q2Q4L + flux3vL*q3Q4L 1775 & + flux3wL*q4Q4L + flux3pL*q5Q4L 1776! 1777 dfp(3,5) = flux3rL*q1Q5L + flux3uL*q2Q5L + flux3vL*q3Q5L 1778 & + flux3wL*q4Q5L + flux3pL*q5Q5L 1779! 1780! 1781 dfp(4,1) = flux4rL*q1Q1L + flux4uL*q2Q1L + flux4vL*q3Q1L 1782 & + flux4wL*q4Q1L + flux4pL*q5Q1L 1783! 1784 dfp(4,2) = flux4rL*q1Q2L + flux4uL*q2Q2L + flux4vL*q3Q2L 1785 & + flux4wL*q4Q2L + flux4pL*q5Q2L 1786! 1787 dfp(4,3) = flux4rL*q1Q3L + flux4uL*q2Q3L + flux4vL*q3Q3L 1788 & + flux4wL*q4Q3L + flux4pL*q5Q3L 1789! 1790 dfp(4,4) = flux4rL*q1Q4L + flux4uL*q2Q4L + flux4vL*q3Q4L 1791 & + flux4wL*q4Q4L + flux4pL*q5Q4L 1792! 1793 dfp(4,5) = flux4rL*q1Q5L + flux4uL*q2Q5L + flux4vL*q3Q5L 1794 & + flux4wL*q4Q5L + flux4pL*q5Q5L 1795! 1796! 1797 dfp(5,1) = flux5rL*q1Q1L + flux5uL*q2Q1L + flux5vL*q3Q1L 1798 & + flux5wL*q4Q1L + flux5pL*q5Q1L 1799! 1800 dfp(5,2) = flux5rL*q1Q2L + flux5uL*q2Q2L + flux5vL*q3Q2L 1801 & + flux5wL*q4Q2L + flux5pL*q5Q2L 1802! 1803 dfp(5,3) = flux5rL*q1Q3L + flux5uL*q2Q3L + flux5vL*q3Q3L 1804 & + flux5wL*q4Q3L + flux5pL*q5Q3L 1805! 1806 dfp(5,4) = flux5rL*q1Q4L + flux5uL*q2Q4L + flux5vL*q3Q4L 1807 & + flux5wL*q4Q4L + flux5pL*q5Q4L 1808! 1809 dfp(5,5) = flux5rL*q1Q5L + flux5uL*q2Q5L + flux5vL*q3Q5L 1810 & + flux5wL*q4Q5L + flux5pL*q5Q5L 1811! 1812! Now the right transformation 1813! 1814 q1Q1R = 1. 1815 q1Q2R = 0. 1816 q1Q3R = 0. 1817 q1Q4R = 0. 1818 q1Q5R = 0. 1819! 1820 q2Q1R = - ur / rhor 1821 q2Q2R = 1. / rhor 1822 q2Q3R = 0. 1823 q2Q4R = 0. 1824 q2Q5R = 0. 1825! 1826 q3Q1R = - vr / rhor 1827 q3Q2R = 0. 1828 q3Q3R = 1. / rhor 1829 q3Q4R = 0. 1830 q3Q5R = 0. 1831! 1832 q4Q1R = - wr / rhor 1833 q4Q2R = 0. 1834 q4Q3R = 0. 1835 q4Q4R = 1. / rhor 1836 q4Q5R = 0. 1837! 1838 q5Q1R = gm1 / 2. * (ur*ur + vr*vr + wr*wr) 1839 q5Q2R = - gm1 * ur 1840 q5Q3R = - gm1 * vr 1841 q5Q4R = - gm1 * wr 1842 q5Q5R = gm1 1843! 1844! 1845 dfm(1,1) = flux1rR*q1Q1R + flux1uR*q2Q1R + flux1vR*q3Q1R 1846 & + flux1wR*q4Q1R + flux1pR*q5Q1R 1847! 1848 dfm(1,2) = flux1rR*q1Q2R + flux1uR*q2Q2R + flux1vR*q3Q2R 1849 & + flux1wR*q4Q2R + flux1pR*q5Q2R 1850! 1851 dfm(1,3) = flux1rR*q1Q3R + flux1uR*q2Q3R + flux1vR*q3Q3R 1852 & + flux1wR*q4Q3R + flux1pR*q5Q3R 1853! 1854 dfm(1,4) = flux1rR*q1Q4R + flux1uR*q2Q4R + flux1vR*q3Q4R 1855 & + flux1wR*q4Q4R + flux1pR*q5Q4R 1856! 1857 dfm(1,5) = flux1rR*q1Q5R + flux1uR*q2Q5R + flux1vR*q3Q5R 1858 & + flux1wR*q4Q5R + flux1pR*q5Q5R 1859! 1860! 1861 dfm(2,1) = flux2rR*q1Q1R + flux2uR*q2Q1R + flux2vR*q3Q1R 1862 & + flux2wR*q4Q1R + flux2pR*q5Q1R 1863! 1864 dfm(2,2) = flux2rR*q1Q2R + flux2uR*q2Q2R + flux2vR*q3Q2R 1865 & + flux2wR*q4Q2R + flux2pR*q5Q2R 1866! 1867 dfm(2,3) = flux2rR*q1Q3R + flux2uR*q2Q3R + flux2vR*q3Q3R 1868 & + flux2wR*q4Q3R + flux2pR*q5Q3R 1869! 1870 dfm(2,4) = flux2rR*q1Q4R + flux2uR*q2Q4R + flux2vR*q3Q4R 1871 & + flux2wR*q4Q4R + flux2pR*q5Q4R 1872! 1873 dfm(2,5) = flux2rR*q1Q5R + flux2uR*q2Q5R + flux2vR*q3Q5R 1874 & + flux2wR*q4Q5R + flux2pR*q5Q5R 1875! 1876! 1877 dfm(3,1) = flux3rR*q1Q1R + flux3uR*q2Q1R + flux3vR*q3Q1R 1878 & + flux3wR*q4Q1R + flux3pR*q5Q1R 1879! 1880 dfm(3,2) = flux3rR*q1Q2R + flux3uR*q2Q2R + flux3vR*q3Q2R 1881 & + flux3wR*q4Q2R + flux3pR*q5Q2R 1882! 1883 dfm(3,3) = flux3rR*q1Q3R + flux3uR*q2Q3R + flux3vR*q3Q3R 1884 & + flux3wR*q4Q3R + flux3pR*q5Q3R 1885! 1886 dfm(3,4) = flux3rR*q1Q4R + flux3uR*q2Q4R + flux3vR*q3Q4R 1887 & + flux3wR*q4Q4R + flux3pR*q5Q4R 1888! 1889 dfm(3,5) = flux3rR*q1Q5R + flux3uR*q2Q5R + flux3vR*q3Q5R 1890 & + flux3wR*q4Q5R + flux3pR*q5Q5R 1891! 1892! 1893 dfm(4,1) = flux4rR*q1Q1R + flux4uR*q2Q1R + flux4vR*q3Q1R 1894 & + flux4wR*q4Q1R + flux4pR*q5Q1R 1895! 1896 dfm(4,2) = flux4rR*q1Q2R + flux4uR*q2Q2R + flux4vR*q3Q2R 1897 & + flux4wR*q4Q2R + flux4pR*q5Q2R 1898! 1899 dfm(4,3) = flux4rR*q1Q3R + flux4uR*q2Q3R + flux4vR*q3Q3R 1900 & + flux4wR*q4Q3R + flux4pR*q5Q3R 1901! 1902 dfm(4,4) = flux4rR*q1Q4R + flux4uR*q2Q4R + flux4vR*q3Q4R 1903 & + flux4wR*q4Q4R + flux4pR*q5Q4R 1904! 1905 dfm(4,5) = flux4rR*q1Q5R + flux4uR*q2Q5R + flux4vR*q3Q5R 1906 & + flux4wR*q4Q5R + flux4pR*q5Q5R 1907! 1908! 1909 dfm(5,1) = flux5rR*q1Q1R + flux5uR*q2Q1R + flux5vR*q3Q1R 1910 & + flux5wR*q4Q1R + flux5pR*q5Q1R 1911! 1912 dfm(5,2) = flux5rR*q1Q2R + flux5uR*q2Q2R + flux5vR*q3Q2R 1913 & + flux5wR*q4Q2R + flux5pR*q5Q2R 1914! 1915 dfm(5,3) = flux5rR*q1Q3R + flux5uR*q2Q3R + flux5vR*q3Q3R 1916 & + flux5wR*q4Q3R + flux5pR*q5Q3R 1917! 1918 dfm(5,4) = flux5rR*q1Q4R + flux5uR*q2Q4R + flux5vR*q3Q4R 1919 & + flux5wR*q4Q4R + flux5pR*q5Q4R 1920! 1921 dfm(5,5) = flux5rR*q1Q5R + flux5uR*q2Q5R + flux5vR*q3Q5R 1922 & + flux5wR*q4Q5R + flux5pR*q5Q5R 1923 flops = flops + 4446.0 1924! 1925! Now take care of contribution to node 1 1926! 1927! idiag = iau(node1) 1928! 1929! Diagonal piece 1930! 1931 if (node1 .le. nnodes) then 1932! do j = 1,5 1933! irow(j) = 5*(node1-1)+j-1 1934! icol(j) = irow(j) 1935! do k = 1,5 1936! in = k + 10*(j-1) 1937! val(in) = dfp(j,k) 1938! enddo 1939! enddo 1940! call MatSetValuesLocal(A,5,irow,5,irow,val,ADD_VALUES,ierr) 1941! call CHK_ERR(irank,ierr,irow(1),irow(1)) 1942 1943! 1944! A(idiag,1,1) = A(idiag,1,1) + dfp(1,1) 1945! A(idiag,1,2) = A(idiag,1,2) + dfp(1,2) 1946! A(idiag,1,3) = A(idiag,1,3) + dfp(1,3) 1947! A(idiag,1,4) = A(idiag,1,4) + dfp(1,4) 1948! 1949! A(idiag,2,1) = A(idiag,2,1) + dfp(2,1) 1950! A(idiag,2,2) = A(idiag,2,2) + dfp(2,2) 1951! A(idiag,2,3) = A(idiag,2,3) + dfp(2,3) 1952! A(idiag,2,4) = A(idiag,2,4) + dfp(2,4) 1953! 1954! A(idiag,3,1) = A(idiag,3,1) + dfp(3,1) 1955! A(idiag,3,2) = A(idiag,3,2) + dfp(3,2) 1956! A(idiag,3,3) = A(idiag,3,3) + dfp(3,3) 1957! A(idiag,3,4) = A(idiag,3,4) + dfp(3,4) 1958! 1959! A(idiag,4,1) = A(idiag,4,1) + dfp(4,1) 1960! A(idiag,4,2) = A(idiag,4,2) + dfp(4,2) 1961! A(idiag,4,3) = A(idiag,4,3) + dfp(4,3) 1962! A(idiag,4,4) = A(idiag,4,4) + dfp(4,4) 1963! 1964! A(idiag,5,1) = A(idiag,5,1) + dfp(5,1) 1965! A(idiag,5,2) = A(idiag,5,2) + dfp(5,2) 1966! A(idiag,5,3) = A(idiag,5,3) + dfp(5,3) 1967! A(idiag,5,4) = A(idiag,5,4) + dfp(5,4) 1968! A(idiag,5,5) = A(idiag,5,5) + dfp(5,5) 1969! do j = 1,5 1970! icol(5+j) = 5*(node2-1)+j-1 1971! do k = 1,5 1972! in = k + 5 + 10*(j-1) 1973! val(in) = dfm(j,k) 1974! enddo 1975! enddo 1976#if defined(INTERLACING) 1977#if defined(BLOCKING) 1978 irow(1) = node1 - 1 1979 icol(1) = node1 - 1 1980 icol(2) = node2 - 1 1981 call MatSetValuesBlockedLocal(A,1,irow,2,icol, 1982 > val1,ADD_VALUES,ierr) 1983#else 1984 do j = 1,5 1985 irow(j) = 5*(node1-1)+j-1 1986 icol(j) = irow(j) 1987 icol(5+j) = 5*(node2-1)+j-1 1988 enddo 1989 call MatSetValuesLocal(A,5,irow,10,icol,val1,ADD_VALUES,ierr) 1990#endif 1991#else 1992 do j = 1,5 1993 irow(j) = (node1-1)+(j-1)*nnodes 1994 icol(j) = irow(j) 1995 icol(5+j) = (node2-1)+(j-1)*nnodes 1996 enddo 1997 call MatSetValues(A,5,irow,10,icol,val1,ADD_VALUES,ierr) 1998#endif 1999 endif 2000 2001! 2002! Now grab the offdiagonal 2003! 2004! ioff = fhelp(n,1) 2005! A(ioff,1,1) = A(ioff,1,1) + dfm(1,1) 2006! A(ioff,1,2) = A(ioff,1,2) + dfm(1,2) 2007! A(ioff,1,3) = A(ioff,1,3) + dfm(1,3) 2008! A(ioff,1,4) = A(ioff,1,4) + dfm(1,4) 2009! 2010! A(ioff,2,1) = A(ioff,2,1) + dfm(2,1) 2011! A(ioff,2,2) = A(ioff,2,2) + dfm(2,2) 2012! A(ioff,2,3) = A(ioff,2,3) + dfm(2,3) 2013! A(ioff,2,4) = A(ioff,2,4) + dfm(2,4) 2014! 2015! A(ioff,3,1) = A(ioff,3,1) + dfm(3,1) 2016! A(ioff,3,2) = A(ioff,3,2) + dfm(3,2) 2017! A(ioff,3,3) = A(ioff,3,3) + dfm(3,3) 2018! A(ioff,3,4) = A(ioff,3,4) + dfm(3,4) 2019! 2020! A(ioff,4,1) = A(ioff,4,1) + dfm(4,1) 2021! A(ioff,4,2) = A(ioff,4,2) + dfm(4,2) 2022! A(ioff,4,3) = A(ioff,4,3) + dfm(4,3) 2023! A(ioff,4,4) = A(ioff,4,4) + dfm(4,4) 2024! 2025! A(ioff,5,1) = A(ioff,5,1) + dfm(5,1) 2026! A(ioff,5,2) = A(ioff,5,2) + dfm(5,2) 2027! A(ioff,5,3) = A(ioff,5,3) + dfm(5,3) 2028! A(ioff,5,4) = A(ioff,5,4) + dfm(5,4) 2029! A(ioff,5,5) = A(ioff,5,5) + dfm(5,5) 2030! 2031! Now do the second node 2032! 2033 if (node2 .le. nnodes) then 2034! do j = 1,5 2035! irow(j) = 5*(node2-1)+j-1 2036! icol(j) = irow(j) 2037! do k = 1,5 2038! in = k + 10*(j-1) 2039! val(in) = -dfm(j,k) 2040! enddo 2041! enddo 2042! call MatSetValuesLocal(A,5,irow,5,irow,val,ADD_VALUES,ierr) 2043! call CHK_ERR(irank,ierr,irow(1),irow(1)) 2044 2045! idiag = iau(node2) 2046! A(idiag,1,1) = A(idiag,1,1) - dfm(1,1) 2047! A(idiag,1,2) = A(idiag,1,2) - dfm(1,2) 2048! A(idiag,1,3) = A(idiag,1,3) - dfm(1,3) 2049! A(idiag,1,4) = A(idiag,1,4) - dfm(1,4) 2050! 2051! A(idiag,2,1) = A(idiag,2,1) - dfm(2,1) 2052! A(idiag,2,2) = A(idiag,2,2) - dfm(2,2) 2053! A(idiag,2,3) = A(idiag,2,3) - dfm(2,3) 2054! A(idiag,2,4) = A(idiag,2,4) - dfm(2,4) 2055! 2056! A(idiag,3,1) = A(idiag,3,1) - dfm(3,1) 2057! A(idiag,3,2) = A(idiag,3,2) - dfm(3,2) 2058! A(idiag,3,3) = A(idiag,3,3) - dfm(3,3) 2059! A(idiag,3,4) = A(idiag,3,4) - dfm(3,4) 2060! 2061! A(idiag,4,1) = A(idiag,4,1) - dfm(4,1) 2062! A(idiag,4,2) = A(idiag,4,2) - dfm(4,2) 2063! A(idiag,4,3) = A(idiag,4,3) - dfm(4,3) 2064! A(idiag,4,4) = A(idiag,4,4) - dfm(4,4) 2065! 2066! A(ioff,5,1) = A(ioff,5,1) - dfp(5,1) 2067! A(ioff,5,2) = A(ioff,5,2) - dfp(5,2) 2068! A(ioff,5,3) = A(ioff,5,3) - dfp(5,3) 2069! A(ioff,5,4) = A(ioff,5,4) - dfp(5,4) 2070! A(ioff,5,5) = A(ioff,5,5) - dfp(5,5) 2071! 2072! do j = 1,5 2073! icol(5+j) = 5*(node1-1)+j-1 2074! do k = 1,5 2075! in = k + 5 + 10*(j-1) 2076! val(in) = -dfp(j,k) 2077! enddo 2078! enddo 2079! Exchange elements in place 2080 do j = 1,5 2081 do k = 1,10 2082! temp = -val1(k,j) 2083! val1(k,j) = -val1(k+5,j) 2084! val1(k+5,j) = temp 2085 val1(k,j) = -val1(k,j) 2086 enddo 2087 enddo 2088! 2089! call CHK_ERR(irank,ierr,irow(1),icol(1)) 2090#if defined(INTERLACING) 2091#if defined(BLOCKING) 2092 irow(1) = node2 - 1 2093 icol(1) = node1 - 1 2094 icol(2) = node2 - 1 2095 call MatSetValuesBlockedLocal(A,1,irow,2,icol, 2096 > val1,ADD_VALUES,ierr) 2097#else 2098 do j = 1,5 2099 irow(j) = 5*(node2-1)+j-1 2100 icol(j) = 5*(node1-1)+j-1 2101 icol(5+j) = irow(j) 2102 enddo 2103 call MatSetValuesLocal(A,5,irow,10,icol,val1,ADD_VALUES,ierr) 2104#endif 2105#else 2106 do j = 1,5 2107 irow(j) = (node2-1)+(j-1)*nnodes 2108 icol(j) = (node1-1)+(j-1)*nnodes 2109 icol(5+j) = irow(j) 2110 enddo 2111 call MatSetValues(A,5,irow,10,icol,val1,ADD_VALUES,ierr) 2112#endif 2113 endif 2114 2115! 2116! Now grab the offdiagonal 2117! 2118! ioff = fhelp(n,2) 2119! A(ioff,1,1) = A(ioff,1,1) - dfp(1,1) 2120! A(ioff,1,2) = A(ioff,1,2) - dfp(1,2) 2121! A(ioff,1,3) = A(ioff,1,3) - dfp(1,3) 2122! A(ioff,1,4) = A(ioff,1,4) - dfp(1,4) 2123! 2124! A(ioff,2,1) = A(ioff,2,1) - dfp(2,1) 2125! A(ioff,2,2) = A(ioff,2,2) - dfp(2,2) 2126! A(ioff,2,3) = A(ioff,2,3) - dfp(2,3) 2127! A(ioff,2,4) = A(ioff,2,4) - dfp(2,4) 2128! 2129! A(ioff,3,1) = A(ioff,3,1) - dfp(3,1) 2130! A(ioff,3,2) = A(ioff,3,2) - dfp(3,2) 2131! A(ioff,3,3) = A(ioff,3,3) - dfp(3,3) 2132! A(ioff,3,4) = A(ioff,3,4) - dfp(3,4) 2133! 2134! A(ioff,4,1) = A(ioff,4,1) - dfp(4,1) 2135! A(ioff,4,2) = A(ioff,4,2) - dfp(4,2) 2136! A(ioff,4,3) = A(ioff,4,3) - dfp(4,3) 2137! A(ioff,4,4) = A(ioff,4,4) - dfp(4,4) 2138 2139 endif 2140 1030 continue 2141! 2142! Convert back to conserved variables 2143! 2144 call PTOE(nvertices,qvec) 2145! 2146! Now we have to close the boundaries 2147! Use this if you want the 5/6-1/6 weighting 2148! 2149! fivesix = 5./6. 2150! onesix = 1./6. 2151! do 1040 i = 1,nsface 2152! iface = isface(i) 2153! node1 = eptr(iface,1) 2154! node2 = eptr(iface,2) 2155! 2156! rho1 = qnode(1,node1) 2157! u1 = qnode(2,node1)/rho1 2158! v1 = qnode(3,node1)/rho1 2159! w1 = qnode(4,node1)/rhol 2160! e1 = qnode(5,node1) 2161! q21 = u1*u1 + v1*v1 + w1*w1 2162! 2163! rho2 = qnode(1,node2) 2164! u2 = qnode(2,node2)/rho2 2165! v2 = qnode(3,node2)/rho2 2166! w2 = qnode(4,node2)/rho2 2167! e2 = qnode(5,node2) 2168! q22 = u2*u2 + v2*v2 + w2*w2 2169! 2170! xmean = .5*(x(node1) + x(node2)) 2171! ymean = .5*(y(node1) + y(node2)) 2172! zmean = .5*(z(node1) + z(node2)) 2173! 2174! xnorm = ymean - y(node1) 2175! ynorm = -(xmean - x(node1)) 2176! rlen = sqrt(xnorm*xnorm + ynorm*ynorm) 2177! xnorm = xnorm/rlen 2178! ynorm = ynorm/rlen 2179! 2180! idiag = iau(node1) 2181! ioff = fhelp(iface,1) 2182! 2183! constx = fivesix*rlen*gm1*xnorm 2184! consty = fivesix*rlen*gm1*ynorm 2185! 2186! A(idiag,2,1) = A(idiag,2,1) + constx*q21/2. 2187! A(idiag,2,2) = A(idiag,2,2) - constx*u1 2188! A(idiag,2,3) = A(idiag,2,3) - constx*v1 2189! A(idiag,2,4) = A(idiag,2,4) + constx 2190! 2191! A(idiag,3,1) = A(idiag,3,1) + consty*q21/2. 2192! A(idiag,3,2) = A(idiag,3,2) - consty*u1 2193! A(idiag,3,3) = A(idiag,3,3) - consty*v1 2194! A(idiag,3,4) = A(idiag,3,4) + consty 2195! 2196! Offdiagonal term 2197! 2198! constx = onesix*rlen*gm1*xnorm 2199! consty = onesix*rlen*gm1*ynorm 2200! 2201! A(ioff,2,1) = A(ioff,2,1) + constx*q22/2. 2202! A(ioff,2,2) = A(ioff,2,2) - constx*u2 2203! A(ioff,2,3) = A(ioff,2,3) - constx*v2 2204! A(ioff,2,4) = A(ioff,2,4) + constx 2205! 2206! A(ioff,3,1) = A(ioff,3,1) + consty*q22/2. 2207! A(ioff,3,2) = A(ioff,3,2) - consty*u2 2208! A(ioff,3,3) = A(ioff,3,3) - consty*v2 2209! A(ioff,3,4) = A(ioff,3,4) + consty 2210! 2211! Second node 2212! 2213! idiag = iau(node2) 2214! ioff = fhelp(iface,2) 2215! 2216! constx = fivesix*rlen*gm1*xnorm 2217! consty = fivesix*rlen*gm1*ynorm 2218! 2219! A(idiag,2,1) = A(idiag,2,1) + constx*q22/2. 2220! A(idiag,2,2) = A(idiag,2,2) - constx*u2 2221! A(idiag,2,3) = A(idiag,2,3) - constx*v2 2222! A(idiag,2,4) = A(idiag,2,4) + constx 2223! 2224! A(idiag,3,1) = A(idiag,3,1) + consty*q22/2. 2225! A(idiag,3,2) = A(idiag,3,2) - consty*u2 2226! A(idiag,3,3) = A(idiag,3,3) - consty*v2 2227! A(idiag,3,4) = A(idiag,3,4) + consty 2228! 2229! Offdiagonal term 2230! 2231! constx = onesix*rlen*gm1*xnorm 2232! consty = onesix*rlen*gm1*ynorm 2233! 2234! A(ioff,2,1) = A(ioff,2,1) + constx*q21/2. 2235! A(ioff,2,2) = A(ioff,2,2) - constx*u1 2236! A(ioff,2,3) = A(ioff,2,3) - constx*v1 2237! A(ioff,2,4) = A(ioff,2,4) + constx 2238! 2239! A(ioff,3,1) = A(ioff,3,1) + consty*q21/2. 2240! A(ioff,3,2) = A(ioff,3,2) - consty*u1 2241! A(ioff,3,3) = A(ioff,3,3) - consty*v1 2242! A(ioff,3,4) = A(ioff,3,4) + consty 2243! 2244!1040 continue 2245! 2246! Use this piece if you dont use the 5/6-1/6 weighting 2247! 2248 do 1080 i = 1,nsnode 2249 inode = isnode(i) 2250 if (inode .le. nnodes) then 2251 xnorm = sxn(i) 2252 ynorm = syn(i) 2253 znorm = szn(i) 2254 rlen = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm) 2255 xnorm = xnorm/rlen 2256 ynorm = ynorm/rlen 2257 znorm = znorm/rlen 2258! 2259 rho = qnode(1,inode) 2260 u = qnode(2,inode)/rho 2261 v = qnode(3,inode)/rho 2262 w = qnode(4,inode)/rho 2263 e = qnode(5,inode) 2264 q2 = u*u + v*v + w*w 2265! 2266 constx = rlen*xnorm*gm1 2267 consty = rlen*ynorm*gm1 2268 constz = rlen*znorm*gm1 2269! 2270 val(1) = constx*q2/2.0 2271 val(2) = - constx*u 2272 val(3) = - constx*v 2273 val(4) = - constx*w 2274 val(5) = constx 2275! 2276 val(6) = consty*q2/2.0 2277 val(7) = - consty*u 2278 val(8) = - consty*v 2279 val(9) = - consty*w 2280 val(10) = consty 2281! 2282 val(11) = constz*q2/2.0 2283 val(12) = - constz*u 2284 val(13) = - constz*v 2285 val(14) = - constz*w 2286 val(15) = constz 2287! 2288#if defined(INTERLACING) 2289 irow(1) = 5*(inode-1) + 1 2290 irow(2) = 5*(inode-1) + 2 2291 irow(3) = 5*(inode-1) + 3 2292 icol(1) = 5*(inode - 1) 2293 icol(2) = 5*(inode-1) + 1 2294 icol(3) = 5*(inode-1) + 2 2295 icol(4) = 5*(inode-1) + 3 2296 icol(5) = 5*(inode-1) + 4 2297 call MatSetValuesLocal(A,3,irow,5,icol,val,ADD_VALUES,ierr) 2298#else 2299 irow(1) = inode - 1 + nnodes 2300 irow(2) = inode - 1 + nnodes*2 2301 irow(3) = inode - 1 + nnodes*3 2302 icol(1) = inode - 1 2303 icol(2) = inode - 1 + nnodes 2304 icol(3) = inode - 1 + nnodes*2 2305 icol(4) = inode - 1 + nnodes*3 2306 icol(5) = inode - 1 + nnodes*4 2307 call MatSetValues(A,3,irow,5,icol,val,ADD_VALUES,ierr) 2308#endif 2309 flops = flops + 47.0 2310 endif 2311 2312! 2313! idiag = iau(inode) 2314! A(idiag,2,1) = A(idiag,2,1) + constx*q2/2. 2315! A(idiag,2,2) = A(idiag,2,2) - constx*u 2316! A(idiag,2,3) = A(idiag,2,3) - constx*v 2317! A(idiag,2,4) = A(idiag,2,4) - constx*w 2318! A(idiag,2,5) = A(idiag,2,5) + constx 2319! 2320! A(idiag,3,1) = A(idiag,3,1) + consty*q2/2. 2321! A(idiag,3,2) = A(idiag,3,2) - consty*u 2322! A(idiag,3,3) = A(idiag,3,3) - consty*v 2323! A(idiag,3,4) = A(idiag,3,4) - consty*w 2324! A(idiag,3,5) = A(idiag,3,5) + consty 2325! 2326! A(idiag,4,1) = A(idiag,4,1) + constz*q2/2. 2327! A(idiag,4,2) = A(idiag,4,2) - constz*u 2328! A(idiag,4,3) = A(idiag,4,3) - constz*v 2329! A(idiag,4,4) = A(idiag,4,4) - constz*w 2330! A(idiag,4,5) = A(idiag,4,5) + constz 2331 1080 continue 2332! print *, "Finished doing inviscid nodes" 2333! 2334! Now do viscous faces 2335! 2336 Prandtl = 0.72 2337 Twall = 1.0 2338 Twall = 1.0 + .5*sqrt(Prandtl)*gm1*xmach*xmach 2339 const = twall/ggm1 2340! 2341! do 1050 i = 1,nvnode 2342! inode = ivnode(i) 2343! idiag = iau(inode) 2344! 2345! First zero out all the ones on the row and then fill them in 2346! 2347! jstart = ia(inode) 2348! jend = ia(inode+1) - 1 2349! 2350! do 1060 j = jstart,jend 2351! 2352! If this isn't the diagonal zero it out 2353! (This way we dont disturb the row for the continuity equation 2354! 2355! if (j.ne.idiag) then 2356! A(j,1,1) = 0.0 2357! A(j,1,2) = 0.0 2358! A(j,1,3) = 0.0 2359! A(j,1,4) = 0.0 2360! A(j,1,5) = 0.0 2361! 2362! A(j,2,1) = 0.0 2363! A(j,2,2) = 0.0 2364! A(j,2,3) = 0.0 2365! A(j,2,4) = 0.0 2366! A(j,2,5) = 0.0 2367! 2368! A(j,3,1) = 0.0 2369! A(j,3,2) = 0.0 2370! A(j,3,3) = 0.0 2371! A(j,3,4) = 0.0 2372! A(j,3,5) = 0.0 2373! 2374! A(j,4,1) = 0.0 2375! A(j,4,2) = 0.0 2376! A(j,4,3) = 0.0 2377! A(j,4,4) = 0.0 2378! A(j,4,5) = 0.0 2379! 2380! A(j,5,1) = 0.0 2381! A(j,5,2) = 0.0 2382! A(j,5,3) = 0.0 2383! A(j,5,4) = 0.0 2384! A(j,5,5) = 0.0 2385! 2386! end if 2387!1060 continue 2388! 2389! Now set the diagonal 2390! 2391! A(idiag,2,1) = 0.0 2392! A(idiag,2,2) = 1.0 2393! A(idiag,2,3) = 0.0 2394! A(idiag,2,4) = 0.0 2395! A(idiag,2,5) = 0.0 2396! 2397! A(idiag,3,1) = 0.0 2398! A(idiag,3,2) = 0.0 2399! A(idiag,3,3) = 1.0 2400! A(idiag,3,4) = 0.0 2401! A(idiag,3,5) = 0.0 2402! 2403! A(idiag,4,1) = 0.0 2404! A(idiag,4,2) = 0.0 2405! A(idiag,4,3) = 0.0 2406! A(idiag,4,4) = 1.0 2407! A(idiag,4,5) = 0.0 2408! 2409! A(idiag,5,1) = -const 2410! A(idiag,5,2) = 0.0 2411! A(idiag,5,3) = 0.0 2412! A(idiag,5,4) = 0.0 2413! A(idiag,5,5) = 1.0 2414! 2415!1050 continue 2416! 2417! Now do the farfield 2418! 2419 s0 = c0*c0/(gamma*rho0**gm1) 2420 xgm1 = 1./gm1 2421 flops = flops + 5.0 2422 do 1070 i = 1,nfnode 2423 inode = ifnode(i) 2424 if (inode .le. nnodes) then 2425 xnorm = fxn(i) 2426 ynorm = fyn(i) 2427 znorm = fzn(i) 2428 rlen = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm) 2429 xnorm = xnorm/rlen 2430 ynorm = ynorm/rlen 2431 znorm = znorm/rlen 2432! 2433 rhoi = qnode(1,inode) 2434 ui = qnode(2,inode)/rhoi 2435 vi = qnode(3,inode)/rhoi 2436 wi = qnode(4,inode)/rhoi 2437 ei = qnode(5,inode) 2438 pi = gm1*(ei - .5*rhoi*(ui*ui + vi*vi + wi*wi)) 2439! 2440 unormi = xnorm*ui + ynorm*vi + znorm*wi 2441 unormo = xnorm*u0 + ynorm*v0 + znorm*w0 2442 ci = sqrt(gamma*pi/rhoi) 2443! 2444 rplus = unormi + 2.*ci/gm1 2445 rminus = unormo - 2.*c0/gm1 2446 unormb = .5*(rplus + rminus) 2447 cb = .25*gm1*(rplus - rminus) 2448! 2449! If unormb > 0 this is outflow: take entropy from inside 2450! If unormb < 0 this is inflow: take entropy from outside 2451! 2452 con = .5*gamma*gm1/(rhoi*ci) 2453 cir = con*(ui*ui + vi*vi + wi*wi - ei/rhoi) 2454 cim = -con*ui 2455 cin = -con*vi 2456 cio = -con*wi 2457 cie = con 2458! 2459 pir = .5*gm1*(ui*ui + vi*vi + wi*wi) 2460 pim = -gm1*ui 2461 pin = -gm1*vi 2462 pio = -gm1*wi 2463 pie = gm1 2464! 2465 unormir = -unormi/rhoi 2466 unormim = xnorm/rhoi 2467 unormin = ynorm/rhoi 2468 unormio = znorm/rhoi 2469 unormie = 0.0 2470! 2471 rplusr = unormir + 2.*cir/gm1 2472 rplusm = unormim + 2.*cim/gm1 2473 rplusn = unormin + 2.*cin/gm1 2474 rpluso = unormio + 2.*cio/gm1 2475 rpluse = unormie + 2.*cie/gm1 2476! 2477 rminusr = 0.0 2478 rminusm = 0.0 2479 rminusn = 0.0 2480 rminuso = 0.0 2481 rminuse = 0.0 2482! 2483 unormbr = .5*(rplusr + rminusr) 2484 unormbm = .5*(rplusm + rminusm) 2485 unormbn = .5*(rplusn + rminusn) 2486 unormbo = .5*(rpluso + rminuso) 2487 unormbe = .5*(rpluse + rminuse) 2488! 2489 cbr = .25*gm1*(rplusr - rminusr) 2490 cbm = .25*gm1*(rplusm - rminusm) 2491 cbn = .25*gm1*(rplusn - rminusn) 2492 cbo = .25*gm1*(rpluso - rminuso) 2493 cbe = .25*gm1*(rpluse - rminuse) 2494 flops = flops + 121.0 2495! 2496 if (unormb .gt. 0.0) then 2497 ub = ui + xnorm*(unormb - unormi) 2498 vb = vi + ynorm*(unormb - unormi) 2499 wb = wi + znorm*(unormb - unormi) 2500 sb = ci*ci/(gamma*rhoi**gm1) 2501 rhob = (cb*cb/(gamma*sb))**xgm1 2502 q2b = ub*ub + vb*vb + wb*wb 2503 pb = rhob*cb*cb/gamma 2504 eb = pb/gm1 + .5*rhob*q2b 2505! 2506! Some derivatives 2507! 2508 ubr = -ui/rhoi + xnorm*(unormbr - unormir) 2509 ubm = 1./rhoi + xnorm*(unormbm - unormim) 2510 ubn = xnorm*(unormbn - unormin) 2511 ubo = xnorm*(unormbo - unormio) 2512 ube = xnorm*(unormbe - unormie) 2513! 2514 vbr = -vi/rhoi + ynorm*(unormbr - unormir) 2515 vbm = ynorm*(unormbm - unormim) 2516 vbn = 1./rhoi + ynorm*(unormbn - unormin) 2517 vbo = ynorm*(unormbo - unormio) 2518 vbe = ynorm*(unormbe - unormie) 2519! 2520 wbr = -wi/rhoi + znorm*(unormbr - unormir) 2521 wbm = znorm*(unormbm - unormim) 2522 wbn = znorm*(unormbn - unormin) 2523 wbo = 1./rhoi + znorm*(unormbo - unormio) 2524 wbe = znorm*(unormbe - unormie) 2525! 2526 sbr = (pir - ci*ci)/(rhoi**gamma) 2527 sbm = pim/(rhoi**gamma) 2528 sbn = pin/(rhoi**gamma) 2529 sbo = pio/(rhoi**gamma) 2530 sbe = pie/(rhoi**gamma) 2531! 2532 con = rhob/(cb*gm1) 2533 rhobr = con*(2.*cbr - cb/sb*sbr) 2534 rhobm = con*(2.*cbm - cb/sb*sbm) 2535 rhobn = con*(2.*cbn - cb/sb*sbn) 2536 rhobo = con*(2.*cbo - cb/sb*sbo) 2537 rhobe = con*(2.*cbe - cb/sb*sbe) 2538! 2539 pbr = (2.*rhob*cb*cbr + cb*cb*rhobr)/gamma 2540 pbm = (2.*rhob*cb*cbm + cb*cb*rhobm)/gamma 2541 pbn = (2.*rhob*cb*cbn + cb*cb*rhobn)/gamma 2542 pbo = (2.*rhob*cb*cbo + cb*cb*rhobo)/gamma 2543 pbe = (2.*rhob*cb*cbe + cb*cb*rhobe)/gamma 2544! 2545!234567890c234567890c234567890c234567890c234567890c234567890c23456789012 2546 ebr = pbr/gm1 + .5*(2.*rhob*(ub*ubr + vb*vbr + wb*wbr) 2547 > + q2b*rhobr) 2548 ebm = pbm/gm1 + .5*(2.*rhob*(ub*ubm + vb*vbm + wb*wbm) 2549 > + q2b*rhobm) 2550 ebn = pbn/gm1 + .5*(2.*rhob*(ub*ubn + vb*vbn + wb*wbn) 2551 > + q2b*rhobn) 2552 ebo = pbo/gm1 + .5*(2.*rhob*(ub*ubo + vb*vbo + wb*wbo) 2553 > + q2b*rhobo) 2554 ebe = pbe/gm1 + .5*(2.*rhob*(ub*ube + vb*vbe + wb*wbe) 2555 > + q2b*rhobe) 2556 flops = flops + 208.0 2557! 2558 else 2559 ub = u0 + xnorm*(unormb - unormo) 2560 vb = v0 + ynorm*(unormb - unormo) 2561 wb = w0 + znorm*(unormb - unormo) 2562 sb = s0 2563 rhob = (cb*cb/(gamma*sb))**xgm1 2564 q2b = ub*ub + vb*vb + wb*wb 2565 pb = rhob*cb*cb/gamma 2566 eb = pb/gm1 + .5*rhob*q2b 2567! 2568! Some derivatives 2569! 2570 unormor = 0.0 2571 unormom = 0.0 2572 unormon = 0.0 2573 unormoo = 0.0 2574 unormoe = 0.0 2575! 2576 ubr = xnorm*(unormbr - unormor) 2577 ubm = xnorm*(unormbm - unormom) 2578 ubn = xnorm*(unormbn - unormon) 2579 ubo = xnorm*(unormbo - unormoo) 2580 ube = xnorm*(unormbe - unormoe) 2581! 2582 vbr = ynorm*(unormbr - unormor) 2583 vbm = ynorm*(unormbm - unormom) 2584 vbn = ynorm*(unormbn - unormon) 2585 vbo = ynorm*(unormbo - unormoo) 2586 vbe = ynorm*(unormbe - unormoe) 2587! 2588 wbr = znorm*(unormbr - unormor) 2589 wbm = znorm*(unormbm - unormom) 2590 wbn = znorm*(unormbn - unormon) 2591 wbo = znorm*(unormbo - unormoo) 2592 wbe = znorm*(unormbe - unormoe) 2593! 2594 sbr = 0.0 2595 sbm = 0.0 2596 sbn = 0.0 2597 sbo = 0.0 2598 sbe = 0.0 2599! 2600 con = rhob/(cb*gm1) 2601 rhobr = con*(2.*cbr - cb/sb*sbr) 2602 rhobm = con*(2.*cbm - cb/sb*sbm) 2603 rhobn = con*(2.*cbn - cb/sb*sbn) 2604 rhobo = con*(2.*cbo - cb/sb*sbo) 2605 rhobe = con*(2.*cbe - cb/sb*sbe) 2606! 2607 pbr = (2.*rhob*cb*cbr + cb*cb*rhobr)/gamma 2608 pbm = (2.*rhob*cb*cbm + cb*cb*rhobm)/gamma 2609 pbn = (2.*rhob*cb*cbn + cb*cb*rhobn)/gamma 2610 pbo = (2.*rhob*cb*cbo + cb*cb*rhobo)/gamma 2611 pbe = (2.*rhob*cb*cbe + cb*cb*rhobe)/gamma 2612! 2613 ebr = pbr/gm1 + .5*(2.*rhob*(ub*ubr + vb*vbr + wb*wbr) 2614 > + q2b*rhobr) 2615 ebm = pbm/gm1 + .5*(2.*rhob*(ub*ubm + vb*vbm + wb*wbm) 2616 > + q2b*rhobm) 2617 ebn = pbn/gm1 + .5*(2.*rhob*(ub*ubn + vb*vbn + wb*wbn) 2618 > + q2b*rhobn) 2619 ebo = pbo/gm1 + .5*(2.*rhob*(ub*ubo + vb*vbo + wb*wbo) 2620 > + q2b*rhobo) 2621 ebe = pbe/gm1 + .5*(2.*rhob*(ub*ube + vb*vbe + wb*wbe) 2622 > + q2b*rhobe) 2623 flops = flops + 177.0 2624! 2625 end if 2626! 2627! Now add contribution to lhs 2628! 2629 2630 val(1) = rlen*(rhob*unormbr + unormb*rhobr) 2631 val(2) = rlen*(rhob*unormbm + unormb*rhobm) 2632 val(3) = rlen*(rhob*unormbn + unormb*rhobn) 2633 val(4) = rlen*(rhob*unormbo + unormb*rhobo) 2634 val(5) = rlen*(rhob*unormbe + unormb*rhobe) 2635! 2636 val(6) = rlen*(rhob*ub*unormbr 2637 1 + unormb*(rhob*ubr + ub*rhobr) + xnorm*pbr) 2638 val(7) = rlen*(rhob*ub*unormbm 2639 1 + unormb*(rhob*ubm + ub*rhobm) + xnorm*pbm) 2640 val(8) = rlen*(rhob*ub*unormbn 2641 1 + unormb*(rhob*ubn + ub*rhobn) + xnorm*pbn) 2642 val(9) = rlen*(rhob*ub*unormbo 2643 1 + unormb*(rhob*ubo + ub*rhobo) + xnorm*pbo) 2644 val(10) = rlen*(rhob*ub*unormbe 2645 1 + unormb*(rhob*ube + ub*rhobe) + xnorm*pbe) 2646! 2647 val(11) = rlen*(rhob*vb*unormbr 2648 1 + unormb*(rhob*vbr + vb*rhobr) + ynorm*pbr) 2649 val(12) = rlen*(rhob*vb*unormbm 2650 1 + unormb*(rhob*vbm + vb*rhobm) + ynorm*pbm) 2651 val(13) = rlen*(rhob*vb*unormbn 2652 1 + unormb*(rhob*vbn + vb*rhobn) + ynorm*pbn) 2653 val(14) = rlen*(rhob*vb*unormbo 2654 1 + unormb*(rhob*vbo + vb*rhobo) + ynorm*pbo) 2655 val(15) = rlen*(rhob*vb*unormbe 2656 1 + unormb*(rhob*vbe + vb*rhobe) + ynorm*pbe) 2657! 2658 val(16) = rlen*(rhob*wb*unormbr 2659 1 + unormb*(rhob*wbr + wb*rhobr) + znorm*pbr) 2660 val(17) = rlen*(rhob*wb*unormbm 2661 1 + unormb*(rhob*wbm + wb*rhobm) + znorm*pbm) 2662 val(18) = rlen*(rhob*wb*unormbn 2663 1 + unormb*(rhob*wbn + wb*rhobn) + znorm*pbn) 2664 val(19) = rlen*(rhob*wb*unormbo 2665 1 + unormb*(rhob*wbo + wb*rhobo) + znorm*pbo) 2666 val(20) = rlen*(rhob*wb*unormbe 2667 1 + unormb*(rhob*wbe + wb*rhobe) + znorm*pbe) 2668! 2669 val(21) = rlen*((eb + pb)*unormbr 2670 1 + unormb*(ebr + pbr)) 2671 val(22) = rlen*((eb + pb)*unormbm 2672 1 + unormb*(ebm + pbm)) 2673 val(23) = rlen*((eb + pb)*unormbn 2674 1 + unormb*(ebn + pbn)) 2675 val(24) = rlen*((eb + pb)*unormbo 2676 1 + unormb*(ebo + pbo)) 2677 val(25) = rlen*((eb + pb)*unormbe 2678 1 + unormb*(ebe + pbe)) 2679! 2680#if defined(INTERLACING) 2681#if defined(BLOCKING) 2682 irow(1) = inode - 1 2683 call MatSetValuesBlockedLocal(A,1,irow,1,irow, 2684 > val,ADD_VALUES,ierr) 2685#else 2686 do k = 1,5 2687 irow(k) = 5*(inode-1)+k-1 2688 enddo 2689 call MatSetValuesLocal(A,5,irow,5,irow,val,ADD_VALUES,ierr) 2690#endif 2691#else 2692 do k = 1,5 2693 irow(k) = inode - 1 + nnodes*(k-1) 2694 enddo 2695 call MatSetValues(A,5,irow,5,irow,val,ADD_VALUES,ierr) 2696#endif 2697 2698 flops = flops + 220.0 2699 endif 2700! 2701 1070 continue 2702 2703! print *, "Finished doing far field nodes" 2704 2705! Assemble matrix 2706 call PetscLogFlops(flops,ierr) 2707 2708 call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) 2709 call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) 2710! call MatView(A, VIEWER_STDOUT_SELF,ierr) 2711 flag = SAME_NONZERO_PATTERN 2712! 2713! End of subroutine FILLA 2714! 2715 return 2716 end 2717! 2718! 2719 2720 subroutine CHK_ERR(irank, ierr, irow, icol) 2721 if (ierr .gt. 0) then 2722 write(*,*) 'On processor ',irank, ': Non-zero entry in row ', 2723 1 irow, ' and column ',icol,' is beyond the pre-allocated memory' 2724 endif 2725 return 2726 end 2727 2728 2729 2730