1 subroutine e3Res ( u1, u2, u3, 2 & uBar, aci, WdetJ, 3 & g1yi, g2yi, g3yi, 4 & rLui, rmu, rho, 5 & tauC, tauM, tauBar, 6 & shpfun, shg, src, 7 & rl, pres, acl, 8 & rlsli) 9c------------------------------------------------------------------------ 10c 11c This routine computes the residual vector at an 12c integration point. 13c 14c input: 15c u1(npro) : x1-velocity 16c u2(npro) : x2-velocity 17c u3(npro) : x3-velocity 18c uBar(npro,3) : u - tauM * Li 19c aci(npro,3) : acceleration 20c rlsli(npro,6) : resolved Leonard stresses 21c WdetJ(npro) : weighted jacobian determinant 22c g1yi(npro,ndof) : x1-gradient of variables 23c g2yi(npro,ndof) : x2-gradient of variables 24c g3yi(npro,ndof) : x3-gradient of variables 25c rLui(npro,3) : total residual of NS equations 26c rmu(npro) : fluid viscosity 27c rho(npro) : density 28c tauC(npro) : continuity tau 29c tauM(npro) : momentum tau 30c tauBar(npro) : additional tau 31c shpfun(npro,nshl) : element shape functions 32c shg(npro,nshl,nsd) : global grad of element shape functions 33c src(npro,nsd) : body force term 34c 35c output: 36c rl(npro,nshl,nflow) 37c 38c------------------------------------------------------------------------ 39 include "common.h" 40 41 dimension u1(npro), u2(npro), u3(npro), 42 & uBar(npro,nsd), aci(npro,nsd), WdetJ(npro), 43 & g1yi(npro,nflow), g2yi(npro,nflow), g3yi(npro,nflow), 44 & rLui(npro,nsd), rmu(npro), rho(npro), 45 & tauC(npro), tauM(npro), tauBar(npro), 46 & shpfun(npro,nshl),shg(npro,nshl,nsd), src(npro,nsd), 47 & pres(npro) 48 49 dimension rl(npro,nshl,nflow), 50 & acl(npro,nshl,ndof), 51 & rlsli(npro,6) 52c 53c.... local declarations 54c 55 real*8 tmp1(npro), tmp2(npro), tmp3(npro), 56 & tmp(npro), rGNa(npro,nsd,nsd), rNa(npro,nsd), 57 & locmass(npro,nshl),omega(3) 58 59 integer aa 60c 61c.... initialize multipliers for Na and Na_{,i} 62c 63 rNa = zero 64 rGNa = zero 65c 66c.... compute the Na multiplier 67c 68 tmps = one-flmpr ! consistant mass factor 69 70c 71c no density yet...it comes later 72c 73 rNa(:,1) = aci(:,1) * tmps 74 & - src(:,1) 75 rNa(:,2) = aci(:,2) * tmps 76 & - src(:,2) 77 rNa(:,3) = aci(:,3) * tmps 78 & - src(:,3) 79 80c 81c.... rotating frame terms if needed 82c 83 84 if(matflg(6,1).eq.1) then ! rotation 85 86 omega(1)=datmat(1,6,1) 87 omega(2)=datmat(2,6,1) 88 omega(3)=datmat(3,6,1) 89c 90c no density yet...it comes later 91c 92 93 rNa(:,1) = rNa(:,1) + (omega(2)-omega(3)) * tauM * rLui(:,1) 94 rNa(:,2) = rNa(:,2) + (omega(3)-omega(1)) * tauM * rLui(:,2) 95 rNa(:,3) = rNa(:,3) + (omega(1)-omega(2)) * tauM * rLui(:,3) 96 endif 97 98 99c 100c.... compute the Na,i multiplier 101c 102 tmp = -pres + tauC * (g1yi(:,2) + g2yi(:,3) + g3yi(:,4)) 103 tmp1 = rmu * ( g2yi(:,2) + g1yi(:,3) ) 104 tmp2 = rmu * ( g3yi(:,3) + g2yi(:,4) ) 105 tmp3 = rmu * ( g1yi(:,4) + g3yi(:,2) ) 106 107 108 if(iconvflow.eq.2) then ! advective form (NO IBP either) 109c 110c no density yet...it comes later 111c 112 rNa(:,1) = rNa(:,1) 113 & + ubar(:,1) * g1yi(:,2) 114 & + ubar(:,2) * g2yi(:,2) 115 & + ubar(:,3) * g3yi(:,2) 116 rNa(:,2) = rNa(:,2) 117 & + ubar(:,1) * g1yi(:,3) 118 & + ubar(:,2) * g2yi(:,3) 119 & + ubar(:,3) * g3yi(:,3) 120 rNa(:,3) = rNa(:,3) 121 & + ubar(:,1) * g1yi(:,4) 122 & + ubar(:,2) * g2yi(:,4) 123 & + ubar(:,3) * g3yi(:,4) 124 125 rGNa(:,1,1) = two * rmu * g1yi(:,2) + tmp 126 rGNa(:,1,2) = tmp1 127 rGNa(:,1,3) = tmp3 128 rGNa(:,2,1) = tmp1 129 rGNa(:,2,2) = two * rmu * g2yi(:,3) + tmp 130 rGNa(:,2,3) = tmp2 131 rGNa(:,3,1) = tmp3 132 rGNa(:,3,2) = tmp2 133 rGNa(:,3,3) = two * rmu * g3yi(:,4) + tmp 134 else ! conservative form (with IBP) 135 136c IBP conservative convection 137c ||||| 138c vvvvv 139 rGNa(:,1,1) = two * rmu * g1yi(:,2) + tmp - u1(:)*u1(:)*rho(:) 140 rGNa(:,1,2) = tmp1 - u1(:)*u2(:)*rho(:) 141 rGNa(:,1,3) = tmp3 - u1(:)*u3(:)*rho(:) 142 rGNa(:,2,1) = tmp1 - u1(:)*u2(:)*rho(:) 143 rGNa(:,2,2) = two * rmu * g2yi(:,3) + tmp - u2(:)*u2(:)*rho(:) 144 rGNa(:,2,3) = tmp2 - u3(:)*u2(:)*rho(:) 145 rGNa(:,3,1) = tmp3 - u1(:)*u3(:)*rho(:) 146 rGNa(:,3,2) = tmp2 - u3(:)*u2(:)*rho(:) 147 rGNa(:,3,3) = two * rmu * g3yi(:,4) + tmp - u3(:)*u3(:)*rho(:) 148 endif 149 if((iLES.gt.10).and.(iLES.lt.20)) then ! bard 150 rGNa(:,1,1) = rGNa(:,1,1) - rlsli(:,1)*rho(:) 151 rGNa(:,1,2) = rGNa(:,1,2) - rlsli(:,4)*rho(:) 152 rGNa(:,1,3) = rGNa(:,1,3) - rlsli(:,5)*rho(:) 153 rGNa(:,2,1) = rGNa(:,2,1) - rlsli(:,4)*rho(:) 154 rGNa(:,2,2) = rGNa(:,2,2) - rlsli(:,2)*rho(:) 155 rGNa(:,2,3) = rGNa(:,2,3) - rlsli(:,6)*rho(:) 156 rGNa(:,3,1) = rGNa(:,3,1) - rlsli(:,5)*rho(:) 157 rGNa(:,3,2) = rGNa(:,3,2) - rlsli(:,6)*rho(:) 158 rGNa(:,3,3) = rGNa(:,3,3) - rlsli(:,3)*rho(:) 159 endif 160 161 tmp1 = tauM * rLui(:,1) 162 tmp2 = tauM * rLui(:,2) 163 tmp3 = tauM * rLui(:,3) 164 165 rGNa(:,1,1) = rGNa(:,1,1) + tmp1 * u1 166 rGNa(:,1,2) = rGNa(:,1,2) + tmp1 * u2 167 rGNa(:,1,3) = rGNa(:,1,3) + tmp1 * u3 168 rGNa(:,2,1) = rGNa(:,2,1) + tmp2 * u1 169 rGNa(:,2,2) = rGNa(:,2,2) + tmp2 * u2 170 rGNa(:,2,3) = rGNa(:,2,3) + tmp2 * u3 171 rGNa(:,3,1) = rGNa(:,3,1) + tmp3 * u1 172 rGNa(:,3,2) = rGNa(:,3,2) + tmp3 * u2 173 rGNa(:,3,3) = rGNa(:,3,3) + tmp3 * u3 174 175 if(iconvflow.eq.1) then 176c 177c... get the u_j w_{i,i} term in there to match A_j^T w_{i,j} tau L_i 178c to match the SUPG of incompressible limit 179c 180 rGNa(:,1,1) = rGNa(:,1,1) + tmp1 * u1 181 rGNa(:,1,2) = rGNa(:,1,2) + tmp2 * u1 182 rGNa(:,1,3) = rGNa(:,1,3) + tmp3 * u1 183 rGNa(:,2,1) = rGNa(:,2,1) + tmp1 * u2 184 rGNa(:,2,2) = rGNa(:,2,2) + tmp2 * u2 185 rGNa(:,2,3) = rGNa(:,2,3) + tmp3 * u2 186 rGNa(:,3,1) = rGNa(:,3,1) + tmp1 * u3 187 rGNa(:,3,2) = rGNa(:,3,2) + tmp2 * u3 188 rGNa(:,3,3) = rGNa(:,3,3) + tmp3 * u3 189 endif 190 191 if(iconvflow.eq.2) then ! advective form has a taubar term to restore con 192 tmp1 = tauBar 193 & * ( rLui(:,1) * g1yi(:,2) 194 & + rLui(:,2) * g2yi(:,2) 195 & + rLui(:,3) * g3yi(:,2) ) 196 tmp2 = tauBar 197 & * ( rLui(:,1) * g1yi(:,3) 198 & + rLui(:,2) * g2yi(:,3) 199 & + rLui(:,3) * g3yi(:,3) ) 200 tmp3 = tauBar 201 & * ( rLui(:,1) * g1yi(:,4) 202 & + rLui(:,2) * g2yi(:,4) 203 & + rLui(:,3) * g3yi(:,4) ) 204 205 rGNa(:,1,1) = rGNa(:,1,1) + tmp1 * rLui(:,1) 206 rGNa(:,1,2) = rGNa(:,1,2) + tmp1 * rLui(:,2) 207 rGNa(:,1,3) = rGNa(:,1,3) + tmp1 * rLui(:,3) 208 rGNa(:,2,1) = rGNa(:,2,1) + tmp2 * rLui(:,1) 209 rGNa(:,2,2) = rGNa(:,2,2) + tmp2 * rLui(:,2) 210 rGNa(:,2,3) = rGNa(:,2,3) + tmp2 * rLui(:,3) 211 rGNa(:,3,1) = rGNa(:,3,1) + tmp3 * rLui(:,1) 212 rGNa(:,3,2) = rGNa(:,3,2) + tmp3 * rLui(:,2) 213 rGNa(:,3,3) = rGNa(:,3,3) + tmp3 * rLui(:,3) 214 endif ! end of advective form 215c 216c.... everything that gets multiplied by rNa was supposed 217c to have density multiplying it. Do it now. 218 219 rNa(:,1) = rNa(:,1) * rho 220 rNa(:,2) = rNa(:,2) * rho 221 rNa(:,3) = rNa(:,3) * rho 222 223c 224c 225c.... multiply the residual pieces by the weight space 226c 227 do aa = 1,nshl 228c 229c.... continuity 230c 231 rl(:,aa,4) = rl(:,aa,4) + WdetJ 232 & * ( shg(:,aa,1) * uBar(:,1) 233 & + shg(:,aa,2) * uBar(:,2) 234 & + shg(:,aa,3) * uBar(:,3) ) 235c 236c.... momentum 237c 238 rl(:,aa,1) = rl(:,aa,1) - WdetJ 239 & * ( shpfun(:,aa) * rNa(:,1) 240 & + shg(:,aa,1) * rGNa(:,1,1) 241 & + shg(:,aa,2) * rGNa(:,1,2) 242 & + shg(:,aa,3) * rGNa(:,1,3) ) 243 rl(:,aa,2) = rl(:,aa,2) - WdetJ 244 & * ( shpfun(:,aa) * rNa(:,2) 245 & + shg(:,aa,1) * rGNa(:,2,1) 246 & + shg(:,aa,2) * rGNa(:,2,2) 247 & + shg(:,aa,3) * rGNa(:,2,3) ) 248 rl(:,aa,3) = rl(:,aa,3) - WdetJ 249 & * ( shpfun(:,aa) * rNa(:,3) 250 & + shg(:,aa,1) * rGNa(:,3,1) 251 & + shg(:,aa,2) * rGNa(:,3,2) 252 & + shg(:,aa,3) * rGNa(:,3,3) ) 253 254 enddo 255c 256c.... return 257c 258 return 259 end 260 261 262 263c------------------------------------------------------------------------ 264c 265c calculate the residual for the advection-diffusion equation 266c 267c------------------------------------------------------------------------ 268 subroutine e3ResSclr ( uMod, gGradS, 269 & Sclr, Sdot, gradS, 270 & WdetJ, rLS, tauS, 271 & shpfun, shg, src, 272 & diffus, 273 & rl ) 274c 275 include "common.h" 276 277 real*8 uMod(npro,nsd), gGradS(npro, nsd), 278 & Sclr(npro), Sdot(npro), gradS(npro,nsd), 279 & WdetJ(npro), rLS(npro), rho(npro), 280 & tauS(npro), shpfun(npro,nshl), src(npro), 281 & shg(npro,nshl,3), rl(npro,nshl) 282 283 real*8 diffus(npro) 284c 285c.... local declarations 286c 287 real*8 rGNa(npro,nsd), rNa(npro), rcp(npro), tmp(npro) 288 289 integer aa 290c 291c.... initialize multipliers for Na and Na_{,i} 292c 293 rNa = zero 294 rGNa = zero 295c 296c.... Na multiplier 297c 298 tmps = one-flmpr ! consistant mass factor 299 rcp = one ! rho * cp 300 301 302 rNa = rcp*(tmps*Sdot + uMod(:,1) * gradS(:,1) 303 & + uMod(:,2) * gradS(:,2) 304 & + uMod(:,3) * gradS(:,3) ) 305 & - src 306 307 308 309 tmp = rcp * tauS * (rLS -src) 310c 311c.... Na,i multiplier 312c 313 rGNa(:,1) = diffus * gradS(:,1) + uMod(:,1) * tmp 314 rGNa(:,2) = diffus * gradS(:,2) + uMod(:,2) * tmp 315 rGNa(:,3) = diffus * gradS(:,3) + uMod(:,3) * tmp 316c 317 if (idcsclr(1) .ne. 0) then 318 if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or. 319 & (idcsclr(2).eq.2 .and. isclr.eq.2)) then ! scalar with dc 320c 321c.... add the contribution of DC to residual 322c 323 rGNa(:,1) = rGNa(:,1) + gGradS(:,1) ! gGradS is 324 rGNa(:,2) = rGNa(:,2) + gGradS(:,2) ! g^{ij}*Y_{j}*dcFct 325 rGNa(:,3) = rGNa(:,3) + gGradS(:,3) ! calculated in e3dc.f 326c 327 endif 328 endif ! end of idcsclr 329c 330c.... multiply the residual pieces by the weight space 331c 332 do aa = 1,nshl 333c 334 rl(:,aa) = rl(:,aa) - WdetJ 335 & * ( shpfun(:,aa) * rNa(:) 336 & + shg(:,aa,1) * rGNa(:,1) 337 & + shg(:,aa,2) * rGNa(:,2) 338 & + shg(:,aa,3) * rGNa(:,3) ) 339 340 enddo 341c 342c.... return 343c 344 return 345 end 346 347 348 349c---------------------------------------------------------------------- 350c Calculate the strong PDE residual. 351c---------------------------------------------------------------------- 352 subroutine e3resStrongPDE( 353 & aci, u1, u2, u3, Temp, rho, xx, 354 & g1yi, g2yi, g3yi, 355 & rLui, src, divqi) 356 include "common.h" 357c INPUTS 358 double precision, intent(in), dimension(npro,nsd) :: 359 & aci, xx 360 double precision, intent(in), dimension(npro,nflow) :: 361 & g1yi, g2yi, g3yi 362 double precision, intent(in), dimension(npro) :: 363 & u1, u2, u3, Temp, rho 364c OUTPUTS 365 double precision, intent(out), dimension(npro,nsd) :: 366 & rLui, src 367c LOCALS 368 double precision, dimension(npro) :: 369 & divu 370 double precision, dimension(npro,nsd) :: 371 & divqi 372 double precision, dimension(nsd) :: 373 & omega 374c.... compute source term 375 src = zero 376 if(matflg(5,1) .ge. 1) then 377c body force contribution to src 378 bfx = datmat(1,5,1) ! Boussinesq, g*alfap 379 bfy = datmat(2,5,1) 380 bfz = datmat(3,5,1) 381 select case ( matflg(5,1) ) 382 case ( 1 ) ! standard linear body force 383 src(:,1) = bfx 384 src(:,2) = bfy 385 src(:,3) = bfz 386 case ( 2 ) ! boussinesq body force 387 Tref = datmat(2,2,1) 388 src(:,1) = bfx * (Temp(:)-Tref) 389 src(:,2) = bfy * (Temp(:)-Tref) 390 src(:,3) = bfz * (Temp(:)-Tref) 391 case ( 3 ) ! user specified f(x,y,z) 392 call e3source(xx, src) 393 end select 394 endif 395c 396 if(matflg(6,1).eq.1) then 397c coriolis force contribution to src 398 omega(1)=datmat(1,6,1) 399 omega(2)=datmat(2,6,1) 400 omega(3)=datmat(3,6,1) 401c note that we calculate f as if it contains the usual source 402c plus the Coriolis and the centrifugal forces taken to the rhs (sign change) 403c as long as we are doing SUPG with no accounting for these terms in the 404c LHS this is the only change (which will find its way to the RHS momentum 405c equation (both Galerkin and SUPG parts)). 406c 407c uncomment later if you want rotation always about z axis 408c orig_src - om x om x r - two om x u 409c 410c$$$ src(:,1)=src(:,1)+omega(3)*omega(3)*xx(:,1)+two*omega(3)*u2 411c$$$ src(:,2)=src(:,2)+omega(3)*omega(3)*xx(:,2)-two*omega(3)*u1 412c 413c more general for testing 414c 415 src(:,1)=src(:,1) 416 & -omega(2)*(omega(1)*xx(:,2)-omega(2)*xx(:,1)) 417 & -omega(3)*(omega(1)*xx(:,3)-omega(3)*xx(:,1)) 418 & -two*(omega(2)*u3-omega(3)*u2) 419 src(:,2)=src(:,2) 420 & -omega(1)*(omega(2)*xx(:,1)-omega(1)*xx(:,2)) 421 & -omega(3)*(omega(2)*xx(:,3)-omega(3)*xx(:,2)) 422 & -two*(omega(3)*u1-omega(1)*u3) 423 src(:,3)=src(:,3) 424 & -omega(1)*(omega(3)*xx(:,1)-omega(1)*xx(:,3)) 425 & -omega(2)*(omega(3)*xx(:,2)-omega(2)*xx(:,3)) 426 & -two*(omega(1)*u2-omega(2)*u1) 427 endif 428c calculate momentum residual 429 rLui(:,1) =(aci(:,1) + u1 * g1yi(:,2) 430 & + u2 * g2yi(:,2) 431 & + u3 * g3yi(:,2) - src(:,1) ) * rho 432 & + g1yi(:,1) 433 & - divqi(:,1) 434 rLui(:,2) =(aci(:,2) + u1 * g1yi(:,3) 435 & + u2 * g2yi(:,3) 436 & + u3 * g3yi(:,3) - src(:,2) ) * rho 437 & + g2yi(:,1) 438 & - divqi(:,2) 439 rLui(:,3) =(aci(:,3) + u1 * g1yi(:,4) 440 & + u2 * g2yi(:,4) 441 & + u3 * g3yi(:,4) - src(:,3) ) * rho 442 & + g3yi(:,1) 443 & - divqi(:,3) 444 if(iconvflow.eq.1) then 445 divu(:) = (g1yi(:,2) + g2yi(:,3) + g3yi(:,4))*rho 446 rLui(:,1)=rlui(:,1)+u1(:)*divu(:) 447 rLui(:,2)=rlui(:,2)+u2(:)*divu(:) 448 rLui(:,3)=rlui(:,3)+u3(:)*divu(:) 449 endif 450c 451 return 452 end subroutine e3resStrongPDE 453