1 subroutine e3LHS ( u1, u2, u3, 2 & uBar, WdetJ, rho, 3 & rLui, rmu, 4 & tauC, tauM, tauBar, 5 & shpfun, shg, xKebe, 6 & xGoC ) 7c------------------------------------------------------------------------ 8c 9c This routine computes the left hand side tangent matrix at an 10c integration point. 11c 12c input: 13c u1(npro) : x1-velocity 14c u2(npro) : x2-velocity 15c u3(npro) : x3-velocity 16c uBar(npro,3) : u - tauM * Li 17c WdetJ(npro) : weighted jacobian determinant 18c rLui(npro,3) : total residual of NS equations 19c rmu(npro) : fluid viscosity 20c rho(npro) : fluid density 21c tauC(npro) : continuity tau 22c tauM(npro) : momentum tau 23c tauBar(npro) : additional tau 24c shpfun(npro,nshl) : element shpfun functions 25c shg(npro,nshl,3) : global grad of element shape functions 26c 27c output: 28c xKebe(npro,9,nshl,nshl) : left hand side 29c xGoC(npro,4,nshl,nshl) : left hand side 30c 31c 32c------------------------------------------------------------------------ 33 include "common.h" 34 35 dimension u1(npro), u2(npro), u3(npro), 36 & uBar(npro,3), WdetJ(npro), rho(npro), 37 & rLui(npro,3), rmu(npro), 38 & tauC(npro), tauM(npro), tauBar(npro), 39 & shpfun(npro,nshl),shg(npro,nshl,3) 40 41 dimension xKebe(npro,9,nshl,nshl), xGoC(npro,4,nshl,nshl) 42c 43c.... local declarations 44c 45 dimension t1(npro,3), t2(npro,3), t3(npro,3), 46 & tmp1(npro), tmp2(npro), 47 & tmp(npro), tlW(npro) 48 49 integer aa, b 50 51 real*8 lhmFct, lhsFct, tsFct(npro) 52 53 lhsFct = alfi * gami * Delt(itseq) 54 lhmFct = almi * (one - flmpl) 55c 56c.... scale variables for efficiency 57c 58 tlW = lhsFct * WdetJ 59 tmp1 = tlW * rho 60 tauM = tlW * tauM 61 tauC = tlW * tauC 62 rmu = tlW * rmu 63 tsFct = lhmFct * WdetJ * rho 64 if(iconvflow.eq.2) then ! 2 is ubar form 3 is cons form but ubar tang. 65 tauBar = lhsFct * WdetJ * tauBar 66 uBar(:,1) = tmp1 * uBar(:,1) 67 uBar(:,2) = tmp1 * uBar(:,2) 68 uBar(:,3) = tmp1 * uBar(:,3) 69 else 70 tauBar = zero !lazy tangent...if effective code it 71 uBar(:,1) = tmp1 * u1(:) 72 uBar(:,2) = tmp1 * u2(:) 73 uBar(:,3) = tmp1 * u3(:) 74 endif 75 76c 77c.... compute mass and convection terms 78c 79 do b = 1, nshl 80 t1(:,1) = uBar(:,1) * shg(:,b,1) 81 & + uBar(:,2) * shg(:,b,2) 82 & + uBar(:,3) * shg(:,b,3) 83c 84c t1=ubar_k N^b,k*rho*alpha_f*gamma*deltat*WdetJ 85c 86 do aa = 1, nshl 87 tmp1 = tsFct * shpfun(:,aa) * shpfun(:,b) 88 tmp2 = tmp1 + t1(:,1) * shpfun(:,aa) 89c 90c tmp1=alpha_m*(1-lmp)*WdetJ*N^aN^b*rho the time term CORRECT 91c tmp2=tmp1+N^a*ubar_k N^b,k*rho*alpha_f*gamma*deltat*WdetJ the 92c second term is convective term CORRECT 93c 94 xKebe(:,1,aa,b) = xKebe(:,1,aa,b) + tmp2 95 xKebe(:,5,aa,b) = xKebe(:,5,aa,b) + tmp2 96 xKebe(:,9,aa,b) = xKebe(:,9,aa,b) + tmp2 97 enddo 98 enddo 99c 100c.... compute the rest of K (symmetric terms) 101c 102 do b = 1, nshl 103 104 t1(:,1) = tauC * shg(:,b,1) 105 t1(:,2) = tauC * shg(:,b,2) 106 t1(:,3) = tauC * shg(:,b,3) 107 108c t1 is tauC*N^b_i,j*alpha_f*gamma*deltat*WdetJ 109 110 t2(:,1) = rmu * shg(:,b,1) 111 t2(:,2) = rmu * shg(:,b,2) 112 t2(:,3) = rmu * shg(:,b,3) 113c t2 is mu*N^b_j,k*alpha_f*gamma*deltat*WdetJ 114 115 tmp1 = tauM * ( u1 * shg(:,b,1) 116 & + u2 * shg(:,b,2) 117 & + u3 * shg(:,b,3) )*rho 118c tmp1 is tauM*(rho u_m N^b_j,m)*alpha_f*gamma*deltat*WdetJ 119 120 tmp2 = tauBar * ( rLui(:,1) * shg(:,b,1) 121 & + rLui(:,2) * shg(:,b,2) 122 & + rLui(:,3) * shg(:,b,3) ) 123c tmp2 is taubar*(L_m N^b_j,m)*alpha_f*gamma*deltat*WdetJ 124 t3(:,1) = t2(:,1) + tmp1 * u1 + tmp2 * rLui(:,1) 125 t3(:,2) = t2(:,2) + tmp1 * u2 + tmp2 * rLui(:,2) 126 t3(:,3) = t2(:,3) + tmp1 * u3 + tmp2 * rLui(:,3) 127 128c t3 is (mu*N^b_j,k + u_k tauM*(rho u_m N^b_j,m)+ L_k*taubar*(L_mN^b_j,m ) 129c *alpha_f*gamma*deltat*WdetJ which isline 2 page 40 of whiting 130c ALMOST (waiting to get hit with N^a_{i,k} 131c mu correct NOW (wrong before) and rho weight on tauM term 132c 133c.... first do the (nodal) diagonal blocks 134c 135 aa = b 136 137 tmp = t3(:,1) * shg(:,aa,1) 138 & + t3(:,2) * shg(:,aa,2) 139 & + t3(:,3) * shg(:,aa,3) 140c previous command is the N^a_{i,k} dot product with t3 defined above 141 142 xKebe(:,1,aa,b) = xKebe(:,1,aa,b) + tmp 143 & + t1(:,1) * shg(:,aa,1) 144 & + t2(:,1) * shg(:,aa,1) 145 xKebe(:,5,aa,b) = xKebe(:,5,aa,b) + tmp 146 & + t1(:,2) * shg(:,aa,2) 147 & + t2(:,2) * shg(:,aa,2) 148 xKebe(:,9,aa,b) = xKebe(:,9,aa,b) + tmp 149 & + t1(:,3) * shg(:,aa,3) 150 & + t2(:,3) * shg(:,aa,3) 151c 152 tmp1 = t1(:,1) * shg(:,aa,2) 153 & + t2(:,2) * shg(:,aa,1) 154 xKebe(:,2,aa,b) = xKebe(:,2,aa,b) + tmp1 155 xKebe(:,4,b,aa) = xKebe(:,4,b,aa) + tmp1 156c 157 tmp1 = t1(:,1) * shg(:,aa,3) 158 & + t2(:,3) * shg(:,aa,1) 159 xKebe(:,3,aa,b) = xKebe(:,3,aa,b) + tmp1 160 xKebe(:,7,b,aa) = xKebe(:,7,b,aa) + tmp1 161c 162 tmp1 = t1(:,2) * shg(:,aa,3) 163 & + t2(:,3) * shg(:,aa,2) 164 xKebe(:,6,aa,b) = xKebe(:,6,aa,b) + tmp1 165 xKebe(:,8,b,aa) = xKebe(:,8,b,aa) + tmp1 166c 167c.... now the off-diagonal (nodal) blocks 168c 169 do aa = b+1, nshl 170 tmp = t3(:,1) * shg(:,aa,1) 171 & + t3(:,2) * shg(:,aa,2) 172 & + t3(:,3) * shg(:,aa,3) 173c 174 tmp1 = tmp 175 & + t1(:,1) * shg(:,aa,1) 176 & + t2(:,1) * shg(:,aa,1) 177 xKebe(:,1,aa,b) = xKebe(:,1,aa,b) + tmp1 178 xKebe(:,1,b,aa) = xKebe(:,1,b,aa) + tmp1 179c 180 tmp1 = tmp 181 & + t1(:,2) * shg(:,aa,2) 182 & + t2(:,2) * shg(:,aa,2) 183 xKebe(:,5,aa,b) = xKebe(:,5,aa,b) + tmp1 184 xKebe(:,5,b,aa) = xKebe(:,5,b,aa) + tmp1 185c 186 tmp1 = tmp 187 & + t1(:,3) * shg(:,aa,3) 188 & + t2(:,3) * shg(:,aa,3) 189 xKebe(:,9,aa,b) = xKebe(:,9,aa,b) + tmp1 190 xKebe(:,9,b,aa) = xKebe(:,9,b,aa) + tmp1 191c 192c.... ( i != j ) 193c 194 tmp1 = t1(:,1) * shg(:,aa,2) 195 & + t2(:,2) * shg(:,aa,1) 196 xKebe(:,2,aa,b) = xKebe(:,2,aa,b) + tmp1 197 xKebe(:,4,b,aa) = xKebe(:,4,b,aa) + tmp1 198c 199 tmp1 = t1(:,1) * shg(:,aa,3) 200 & + t2(:,3) * shg(:,aa,1) 201 xKebe(:,3,aa,b) = xKebe(:,3,aa,b) + tmp1 202 xKebe(:,7,b,aa) = xKebe(:,7,b,aa) + tmp1 203c 204 tmp1 = t1(:,2) * shg(:,aa,1) 205 & + t2(:,1) * shg(:,aa,2) 206 xKebe(:,4,aa,b) = xKebe(:,4,aa,b) + tmp1 207 xKebe(:,2,b,aa) = xKebe(:,2,b,aa) + tmp1 208c 209 tmp1 = t1(:,2) * shg(:,aa,3) 210 & + t2(:,3) * shg(:,aa,2) 211 xKebe(:,6,aa,b) = xKebe(:,6,aa,b) + tmp1 212 xKebe(:,8,b,aa) = xKebe(:,8,b,aa) + tmp1 213c 214 tmp1 = t1(:,3) * shg(:,aa,1) 215 & + t2(:,1) * shg(:,aa,3) 216 xKebe(:,7,aa,b) = xKebe(:,7,aa,b) + tmp1 217 xKebe(:,3,b,aa) = xKebe(:,3,b,aa) + tmp1 218c 219 tmp1 = t1(:,3) * shg(:,aa,2) 220 & + t2(:,2) * shg(:,aa,3) 221 xKebe(:,8,aa,b) = xKebe(:,8,aa,b) + tmp1 222 xKebe(:,6,b,aa) = xKebe(:,6,b,aa) + tmp1 223c 224 enddo 225 enddo 226c 227c.... compute G Nai Nbp,j 228c 229 230 do b = 1, nshl 231 t1(:,1) = tlW * shg(:,b,1) 232 t1(:,2) = tlW * shg(:,b,2) 233 t1(:,3) = tlW * shg(:,b,3) 234 do aa = 1, nshl 235 xGoC(:,1,aa,b) = xGoC(:,1,aa,b) + t1(:,1) * shpfun(:,aa) 236 xGoC(:,2,aa,b) = xGoC(:,2,aa,b) + t1(:,2) * shpfun(:,aa) 237 xGoC(:,3,aa,b) = xGoC(:,3,aa,b) + t1(:,3) * shpfun(:,aa) 238 enddo 239 enddo 240c 241c.... compute C 242c we divide by rho because the L on the weight space is density divided 243c form 244c 245 tauM=tauM/rho 246 do b = 1, nshl 247 t1(:,1) = tauM * shg(:,b,1) 248 t1(:,2) = tauM * shg(:,b,2) 249 t1(:,3) = tauM * shg(:,b,3) 250 do aa = b, nshl 251 xGoC(:,4,aa,b) = xGoC(:,4,aa,b) 252 & + t1(:,1) * shg(:,aa,1) 253 & + t1(:,2) * shg(:,aa,2) 254 & + t1(:,3) * shg(:,aa,3) 255 enddo 256 enddo 257 258c 259c.... return 260c 261 return 262 end 263 264 265c------------------------------------------------------------------------ 266c 267c calculate the tangent matrix for the advection-diffusion equation 268c 269c------------------------------------------------------------------------ 270 subroutine e3LHSSclr ( uMod, giju, dcFct, 271 & Sclr, Sdot, gradS, 272 & WdetJ, rLS, tauS, 273 & shpfun, shg, srcL, 274 & diffus, 275 & xSebe ) 276 277c 278 include "common.h" 279 280 real*8 uMod(npro,nsd), 281 & Sclr(npro), Sdot(npro), gradS(npro,nsd), 282 & WdetJ(npro), rLS(npro), rho(npro), 283 & tauS(npro), shpfun(npro,nshl), 284 & srcL(npro), shg(npro,nshl,3), 285 & xSebe(npro,nshl,nshl) 286 287 real*8 diffus(npro), cp, kptmp(npro),tauSo(npro) 288 289c 290c.... local declarations 291c 292 dimension t1(npro,3), tmp1(npro), tmp2(npro), 293 & tmp(npro), dcFct(npro), giju(npro,6) 294 295 integer aa, b 296 297 real*8 lhsFct, tsFct(npro) 298 299 lhsFct = alfi * gami * Delt(itseq) 300c 301c.... scale variables for efficiency 302c 303 tauSo = tauS 304 tauS = lhsFct * WdetJ * tauS 305 kptmp = lhsFct * WdetJ * diffus 306 tsFct = almi * WdetJ * (one - flmpl) 307 srcL = srcL * WdetJ * lhsFct 308c 309c.... compute mass and convection terms 310c 311 do b = 1, nshl 312 t1(:,1) = WdetJ * ( uMod(:,1) * shg(:,b,1) 313 & + uMod(:,2) * shg(:,b,2) 314 & + uMod(:,3) * shg(:,b,3) ) 315 t1(:,2) = t1(:,1) * tauSo 316 do aa = 1, nshl 317 tmp1 = shpfun(:,aa) * shpfun(:,b) 318 tmp2 = shpfun(:,aa) * lhsFct 319 xSebe(:,aa,b) = xSebe(:,aa,b) + tmp1 * (tsFct + srcL) 320 & + tmp2 * t1(:,1) 321c 322c.... compute mass term for stab u_j N_{a,j} tau N_b (note that a and b 323c flipped on both sides below) 324c 325 xSebe(:,b,aa) = xSebe(:,b,aa) + t1(:,2)*shpfun(:,aa) 326 enddo 327 enddo 328c 329c.... compute the rest of S (symmetric terms) 330c 331 do b = 1, nshl 332 tmp = tauS(:) 333 & * ( uMod(:,1) * shg(:,b,1) 334 & + uMod(:,2) * shg(:,b,2) 335 & + uMod(:,3) * shg(:,b,3) ) 336 337 t1(:,1) = kptmp * shg(:,b,1) + uMod(:,1) * tmp 338 t1(:,2) = kptmp * shg(:,b,2) + uMod(:,2) * tmp 339 t1(:,3) = kptmp * shg(:,b,3) + uMod(:,3) * tmp 340 if (idcsclr(1) .ne. 0) then 341 if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or. 342 & (idcsclr(2).eq.2 .and. isclr.eq.2)) then ! scalar with dc 343c 344 tmp = WdetJ * dcFct * lhsFct 345c 346 giju(:,1) = tmp * giju(:,1) 347 giju(:,2) = tmp * giju(:,2) 348 giju(:,3) = tmp * giju(:,3) 349 giju(:,4) = tmp * giju(:,4) 350 giju(:,5) = tmp * giju(:,5) 351 giju(:,6) = tmp * giju(:,6) 352c 353 t1(:,1) = t1(:,1) + giju(:,1) * shg(:,b,1) 354 2 + giju(:,4) * shg(:,b,2) 355 3 + giju(:,6) * shg(:,b,3) 356 t1(:,2) = t1(:,2) + giju(:,4) * shg(:,b,1) 357 2 + giju(:,2) * shg(:,b,2) 358 3 + giju(:,5) * shg(:,b,3) 359 t1(:,3) = t1(:,3) + giju(:,6) * shg(:,b,1) 360 2 + giju(:,5) * shg(:,b,2) 361 3 + giju(:,3) * shg(:,b,3) 362 endif 363 endif !end of idcsclr 364c 365c.... first do the (nodal) diagonal blocks 366c 367 aa = b 368 369 xSebe(:,aa,b) = xSebe(:,aa,b) + t1(:,1) * shg(:,aa,1) 370 & + t1(:,2) * shg(:,aa,2) 371 & + t1(:,3) * shg(:,aa,3) 372 373c 374c.... now the off-diagonal (nodal) blocks 375c 376 do aa = b+1, nshl 377 tmp = t1(:,1) * shg(:,aa,1) 378 & + t1(:,2) * shg(:,aa,2) 379 & + t1(:,3) * shg(:,aa,3) 380c 381 xSebe(:,aa,b) = xSebe(:,aa,b) + tmp 382 xSebe(:,b,aa) = xSebe(:,b,aa) + tmp 383c 384 enddo 385 enddo 386 387c 388c.... return 389c 390 return 391 end 392 393