1 subroutine e3source(xx, src) 2c----------------------------------------------------------------------- 3c 4c this routine computes the body force term. 5c 6c currently this computes a swirl body with the axis alligned with 7c the z-coordinate 8c 9c----------------------------------------------------------------------- 10 11 include "common.h" 12 13 real*8 xx(npro,nsd), src(npro,nsd) 14 15 real*8 nu 16 17 real*8 r, Stheta, dpdz, rP5 18 if(datmat(2,5,1).eq. 0) then ! calc swirl 19c 20c This is the body force which will drive a swirl in a pipe flow 21c 22 bigR = 0.5d0 23 dpdz = datmat(1,5,1) 24 do iel = 1, npro 25 26 r = sqrt( xx(iel,1)**2 + xx(iel,2)**2) 27 rP5 = (r/bigR)**5 28 29 Stheta = dpdz * sin(0.5*pi*rP5) 30 31 src(iel,1) = -xx(iel,2)/r * Stheta 32 src(iel,2) = xx(iel,1)/r * Stheta 33 src(iel,3) = dpdz 34 enddo 35 else ! contrived test problem 36 37c$$$ nu=datmat(1,2,1) 38c$$$ omag=datmat(3,5,1) ! frame rotation rate 39c$$$ e1 = one/sqrt(3.0d0) ! axis of rotation 40c$$$ e2 = e1 41c$$$ e3 = e1 42c$$$ om1=omag*e1 43c$$$ om2=omag*e2 44c$$$ om3=omag*e3 45c$$$ w=0 46c$$$ 47c$$$ do iel = 1, npro 48c$$$ 49c$$$ x = xx(iel,1) 50c$$$ y = xx(iel,2) 51c$$$ z = xx(iel,3) 52c$$$ 53c$$$c 54c$$$c The following are MAPLE generated forcing functions for 55c$$$c a contrived flow in a rotating reference frame. 56c$$$c 57c$$$ 58c$$$ t1 = x**2 59c$$$ t2 = t1**2 60c$$$ t5 = dexp(14.D0*x) 61c$$$ t6 = t2*x*t5 62c$$$ t7 = y**2 63c$$$ t8 = t7**2 64c$$$ t9 = t8*t7 65c$$$ t12 = t2*t5 66c$$$ t15 = nu*t1 67c$$$ t17 = dexp(7.D0*x) 68c$$$ t18 = t7*y 69c$$$ t19 = t17*t18 70c$$$ t22 = t1*x 71c$$$ t23 = nu*t22 72c$$$ t24 = t17*y 73c$$$ t29 = t17*t7 74c$$$ t34 = nu*x 75c$$$ t43 = nu*t2 76c$$$ t46 = -40.D0*t6*t9-6.D0*t12*t7+92.D0*t15*t19+132.D0*t23*t24+80.D0* 77c$$$ #t6*t18-252.D0*t23*t29+168.D0*t23*t19+96.D0*t34*t29-64.D0*t34*t19+2 78c$$$ #2.D0*t15*t24-138.D0*t15*t29+294.D0*t43*t29 79c$$$ t50 = t2*t22*t5 80c$$$ t57 = nu*t17 81c$$$ t60 = t8*y 82c$$$ t65 = t2**2 83c$$$ t66 = t65*t5 84c$$$ t73 = t22*t5 85c$$$ t77 = t2*t1*t5 86c$$$ t80 = -196.D0*t43*t19-96.D0*t50*t9-122.D0*t43*t24-32.D0*t34*t24-4. 87c$$$ #D0*t57*y+36.D0*t12*t60+192.D0*t50*t18+98.D0*t66*t8+24.D0*t12*t18+2 88c$$$ #8.D0*t66*t9-24.D0*t73*t60-336.D0*t77*t60 89c$$$ t106 = -336.D0*t50*t8-12.D0*t12*t9+8.D0*t73*t9-140.D0*t6*t8+28.D0* 90c$$$ #t73*t8+12.D0*t15*t17+12.D0*t43*t17+288.D0*t50*t60+112.D0*t77*t9-22 91c$$$ #4.D0*t77*t18-56.D0*t66*t18+120.D0*t6*t60 92c$$$ t131 = -8.D0*t57*t18-20.D0*t6*t7+56.D0*t77*t7-42.D0*t12*t8-24.D0*t 93c$$$ #23*t17-48.D0*t50*t7+4.D0*t73*t7-16.D0*t73*t18+392.D0*t77*t8+14.D0* 94c$$$ #t66*t7-84.D0*t66*t60+12.D0*t57*t7 95c$$$ fx = t46+t80+t106+t131 96c$$$ 97c$$$ 98c$$$ t1 = x**2 99c$$$ t2 = t1**2 100c$$$ t3 = nu*t2 101c$$$ t5 = dexp(7.D0*x) 102c$$$ t6 = t5*y 103c$$$ t9 = y**2 104c$$$ t10 = t9**2 105c$$$ t11 = nu*t10 106c$$$ t18 = t1*x 107c$$$ t22 = nu*x 108c$$$ t25 = nu*t18 109c$$$ t32 = dexp(14.D0*x) 110c$$$ t33 = t2*t32 111c$$$ t39 = t2*t1*t32 112c$$$ t40 = t10*t9 113c$$$ t43 = t1*t32 114c$$$ t46 = t9*y 115c$$$ t47 = t10*t46 116c$$$ t50 = -84.D0*t3*t6+343.D0*t11*t5*t2+66.D0*t11*t5*x-98.D0*t11*t5*t1 117c$$$ #8-24.D0*t22*t6+120.D0*t25*t6-287.D0*t11*t5*t1-140.D0*t33*t10+14.D0 118c$$$ #*t3*t5-56.D0*t39*t40-28.D0*t43*t40+8.D0*t43*t47 119c$$$ t52 = t2*x*t32 120c$$$ t55 = t18*t32 121c$$$ t62 = t10*y 122c$$$ t69 = nu*t5 123c$$$ t80 = 120.D0*t52*t10-32.D0*t55*t47+56.D0*t33*t47+30.D0*t11*t5-216. 124c$$$ #D0*t52*t62-144.D0*t55*t62+72.D0*t39*t62-60.D0*t69*t46+30.D0*t69*t9 125c$$$ #-20.D0*t25*t5-20.D0*t43*t10+36.D0*t43*t62 126c$$$ t86 = nu*t1 127c$$$ t89 = t5*t9 128c$$$ t92 = t5*t46 129c$$$ t109 = 112.D0*t55*t40+80.D0*t55*t10-12.D0*t86*t6-275.D0*t86*t89+57 130c$$$ #4.D0*t86*t92-218.D0*t25*t89+196.D0*t25*t92+90.D0*t22*t89-132.D0*t2 131c$$$ #2*t92+427.D0*t3*t89+8.D0*t39*t46+4.D0*t43*t46 132c$$$ t134 = 16.D0*t39*t47+28.D0*t33*t46-48.D0*t52*t47+252.D0*t33*t62+16 133c$$$ #8.D0*t52*t40-24.D0*t52*t46-686.D0*t3*t92+2.D0*t86*t5-40.D0*t39*t10 134c$$$ #-196.D0*t33*t40-16.D0*t55*t46+4.D0*t22*t5 135c$$$ fy = t50+t80+t109+t134 136c$$$ 137c$$$ 138c$$$ t3 = om3*x 139c$$$ t5 = dexp(7.D0*x) 140c$$$ t7 = x**2 141c$$$ t12 = y**2 142c$$$ t15 = (-1.D0+y)**2 143c$$$ fxa = 2.D0*om2*w+2.D0*t3*t5*(2.D0+x-10.D0*t7+7.D0*t7*x)*t12*t15+om 144c$$$ #2*(om1*y-om2*x)-om3*(t3-om1*z) 145c$$$ 146c$$$ t1 = x**2 147c$$$ t4 = (-1.D0+x)**2 148c$$$ t7 = dexp(7.D0*x) 149c$$$ t10 = y**2 150c$$$ fya = 4.D0*om3*t1*t4*t7*y*(1.D0-3.D0*y+2.D0*t10)-2.D0*om1*w+om3*(o 151c$$$ #m2*z-om3*y)-om1*(om1*y-om2*x) 152c$$$ 153c$$$ 154c$$$ t3 = dexp(7.D0*x) 155c$$$ t5 = x**2 156c$$$ t10 = y**2 157c$$$ t13 = (-1.D0+y)**2 158c$$$ t19 = (-1.D0+x)**2 159c$$$ fza = -2.D0*om1*x*t3*(2.D0+x-10.D0*t5+7.D0*t5*x)*t10*t13-4.D0*om2* 160c$$$ #t5*t19*t3*y*(1.D0-3.D0*y+2.D0*t10)+om1*(om3*x-om1*z)-om2*(om2*z-om 161c$$$ #3*y) 162c$$$ 163c$$$ src(iel,1) = fx + fxa 164c$$$ src(iel,2) = fy + fya 165c$$$ src(iel,3) = fza 166c$$$ enddo 167!Analytic lid case 168 169 do iel = 1, npro 170 x = xx(iel,1) 171 y = xx(iel,2) 172 z = xx(iel,3) 173 Re = 1.0/datmat(1,2,1) 174c 175c The following are MAPLE generated forcing functions for 176c a Lid Driven cavity flow with an analytic solution 177c 178 t2 = x**2 179 t3 = t2**2 180 t4 = t3*x 181 t7 = t2*x 182 t8 = 8*t7 183 t13 = y**2 184 t20 = t13**2 185 t21 = t20-t13 186 t26 = t3**2 187 t30 = t3*t2 188 t37 = t13*y 189 t54 = -8/Re*(24.E0/5.E0*t4-12*t3+t8+2*(4*t7-6*t2+2*x)*(12 190 & *t13-2)+(24*x-12)*t21)-64*(t26/2-2*t3*t7+3*t30-2*t4+t3 191 & /2)*(-24*t20*y+8*t37-4*y)+64*t21*(4*t37-2*y)*(-4*t30+12 192 & *t4-14*t3+t8-2*t2) 193 194 src(iel,1) = 0.0 195 src(iel,2) = t54 196 src(iel,3) = 0.0 197 enddo 198 199 endif 200 return 201 end 202 203 204!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 205!****************************************************************** 206!------------------------------------------------------------------- 207 208 209 210 subroutine e3sourceSclr ( Sclr, Sdot, gradS, 211 & dwl, shape_funct, shg, 212 & yl, dxidx, rmu, 213 & u1, u2, u3, xl, 214 & srcR, srcL, uMod, 215 & srcRat) 216 217 218c----------------------------------------------------------------------- 219c 220c calculate the coefficients for the Spalart-Allmaras 221c turbulence model and its jacobian 222c 223c input: 224c Sclr : The scalar that is being transported in this pde. 225c gradS : spatial gradient of Sclr 226c rmu : kinematic viscosity 227c 228c output: 229c rmu : diffusion for eddy viscosity equation 230c srcR : residual terms for 231c (srcR * turb) 232c srcL : jacobian of srcR 233c 234c----------------------------------------------------------------------- 235 use turbSA 236 use turbKE 237 include "common.h" 238c coming in 239 real*8 Sclr(npro), Sdot(npro), 240 & gradS(npro,nsd), dwl(npro,nenl), 241 & shape_funct(npro,nshl), shg(npro,nshl,nsd), 242 & yl(npro,nshl,ndof), dxidx(npro,nsd,nsd), 243 & rmu(npro), u1(npro), 244 & u2(npro), u3(npro), 245 & xl(npro,nenl,nsd) 246c going out 247 real*8 srcR(npro), srcL(npro), 248 & uMod(npro,nsd) 249c used locally 250 251 real*8 gradV(npro,nsd,nsd), absVort(npro), 252 & dwall(npro) 253 real*8 chi, chiP3, fv1, fv2, st, r, 254 & g, fw, s, viscInv, k2d2Inv, 255 & dP2Inv, sixth, tmp, tmp1, p, dp, 256 & fv1p, fv2p, stp, gp, fwp, rp, 257 & chiP2, mytmp(npro), sign_levelset(npro), 258 & sclr_ls(npro) 259!MR CHANGE 260 real*8 SclrNN,qfac,dx,dy,dz,dmax 261 real*8 rd,fd,dwallsqqfact,ep 262!MRCHANGE 263c Kay-Epsilon 264 real*8 Ry(npro), Rt(npro), RtP2(npro), f2(npro), f1(npro), 265 & kay(npro), epsilon(npro), fmu(npro), fmui(npro), 266 & srcjac(npro,4), srcRat(npro) 267 integer e,n 268c---------------------------------------------------------------------- 269c 270c -- -- -- -- 271c | cb2 | 1 | | 272c phi, + |u - -----phi, | phi, - -----|(nu+phi)phi, | 273c t | i sigma i| i sigma| i| 274c -- -- -- --, 275c i 276c 277c ~ phi^2 278c - cb1*S*phi + cw1*fw*----- 279c d^2 280c 281c---------------------------------------------------------------------- 282c 283 if(iRANS.eq.-1) then ! spalart almaras model 284c 285c.... compute the global gradient of u 286c 287 gradV = zero 288 do n = 1, nshl 289c 290c du_i/dx_j 291c 292c i j indices match array where V is the velocity (u in our notes) 293 gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * yl(:,n,2) 294 gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * yl(:,n,3) 295 gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * yl(:,n,4) 296c 297 gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * yl(:,n,2) 298 gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * yl(:,n,3) 299 gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * yl(:,n,4) 300c 301 gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * yl(:,n,2) 302 gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * yl(:,n,3) 303 gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * yl(:,n,4) 304c a j u a i 305c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in yl vector 306c 307 enddo 308c 309c.... magnitude of vorticity 310c 311 absVort = sqrt( (gradV(:,2,3) - gradV(:,3,2)) ** 2 312 & + (gradV(:,3,1) - gradV(:,1,3)) ** 2 313 & + (gradV(:,1,2) - gradV(:,2,1)) ** 2 ) 314 dwall = zero 315 do n = 1, nenl 316 dwall = dwall + shape_funct(:,n) * dwl(:,n) 317 enddo 318 319 sixth = 1.0/6.0 320c 321c.... compute source and its jacobian 322c 323 do e = 1, npro 324 SclrNN= max(Sclr(e),zero) 325c trip after the plane x-z=-1 326c if((xl(e,1,1)-xl(e,1,3)-0.1) .lt. zero) SclrNN=zero 327c trip after the plane whose normal is (0.8660254;0;0.5) and origin is (0.19551661;0;-0.28607881) 328c namely 0.8660254x+0.5z-0.0262829=0 329c if((0.8660254*xl(e,1,1)+0.5*xl(e,1,3)-0.0262829) .lt. zero) 330c & SclrNN=zero 331! if((0.86230254*xl(e,1,1)-0.0842707*xl(e,1,2) 332! +0.49933235*xl(e,1,3)-0.0260094) .lt. zero) 333 334!2W downstream 335! if((0.86230254*xl(e,1,1)-0.0842707*xl(e,1,2) 336! & +0.49933235*xl(e,1,3)-0.0298618) .lt. zero) 337! & SclrNN=zero 338 339!DDES upstream 340! if((0.86230254*xl(e,1,1)-0.0842707*xl(e,1,2) 341! & +0.49933235*xl(e,1,3)+0.03658399) .lt. zero) 342! & SclrNN=zero 343 344 345 346 if(iles.lt.0) then ! for DES97 or DDES 347 348 dx=maxval(xl(e,:,1))-minval(xl(e,:,1)) 349 dy=maxval(xl(e,:,2))-minval(xl(e,:,2)) 350 dz=maxval(xl(e,:,3))-minval(xl(e,:,3)) 351 dmax=max(dx,max(dy,dz)) 352 dmax=0.325*dmax 353 354 if(iles .eq. -1) then ! Original DES97 355 dwall(e) = min(dmax,dwall(e)) 356 357 elseif(iles .eq. -2) then ! DDES 358 359 qfac = sqrt( 360 & gradV(e,1,1)**2+gradV(e,1,2)**2+gradV(e,1,3)**2 361 & + gradV(e,2,1)**2+gradV(e,2,2)**2+gradV(e,2,3)**2 362 & + gradV(e,3,1)**2+gradV(e,3,2)**2+gradV(e,3,3)**2 ) 363 364 ! The following lines fix an issue when all vertices on 365 ! an element are located on a wall with no-slip wall 366 ! conditions imposed, which lead to a divistion by 0 367 368 !rd=SclrNN*saKappaP2Inv/(dwall(e)**2*qfac) 369 dwallsqqfact = max(dwall(e)**2*qfac,1.0d-12) 370 rd = SclrNN*saKappaP2Inv/dwallsqqfact 371 fd = one-tanh((8.0000000000000000d0*rd)**3) 372 dwall(e) = dwall(e)-fd*max(zero,dwall(e)-dmax) 373 374 endif 375 376 endif ! if (iles .lt.0) 377 378 if (dwall(e) .ne. 0) dP2Inv = 1.0 / dwall(e)**2 379 380 k2d2Inv = dP2Inv * saKappaP2Inv 381 382 viscInv = 1.0/datmat(1,2,1) 383 chi = SclrNN * viscInv !skipping the kiy kie for now 384 chiP2 = chi * chi 385 chiP3 = chi * chiP2 386 387 ep=zero 388 if(Sclr(e).gt.zero) ep=one 389 390 tmp = 1.0 / ( chiP3 + saCv1P3 ) 391 fv1 = chiP3 * tmp 392 fv1p = 3.0 * viscInv * saCv1P3 * chiP2 * tmp**2 393 394 tmp = 1.0 / (1.0 + chi * fv1) 395 fv2 = 1.0 - chi * tmp 396 fv2p = (chiP2 * fv1p - viscInv) * tmp**2 397 398 s = absVort(e) !prd(e,2) 399 tmp = SclrNN * k2d2Inv !eOkdP2 400 st = s + fv2 * tmp 401 stp = ep*k2d2Inv * fv2 + tmp*fv2p 402 403 r = zero 404 rp = zero 405 if (st .gt. epsM ) then 406 r = tmp / st 407 r = min( max(r, -8.0d0), 8.0d0) 408 rp = k2d2Inv / st**2 * (ep*st - SclrNN*stp) 409 endif 410 rP5 = r * (r * r) ** 2 411 tmp = 1.0 + saCw2 * (rP5 - 1.0) 412 g = r * tmp 413 gp = rp * (tmp + 5.0 * saCw2 * rP5) 414 415 gP6 = (g * g * g) ** 2 416 tmp = 1.0 / (gP6 + saCw3P6) 417 tmp1 = ( (1.0 + saCw3P6) * tmp ) ** sixth 418 fw = g * tmp1 419 fwp = gp * tmp1 * saCw3P6 * tmp 420 if(abs(r).gt.7.0d0) fwp=zero 421 422 tmp1 = saCw1 * dP2Inv*SclrNN 423 prat = ep*saCb1*st !prodRat 424 srcRat(e)= prat - ep*fw * tmp1 425 426 tmp = zero 427 if(prat .lt. zero) tmp=prat 428 if(stp .lt. zero) tmp=tmp+SaCb1*stp*SclrNN 429 if(fw .gt. zero) tmp=tmp-two*tmp1*fw 430 if(fwp .gt. zero) tmp=tmp-tmp1*SclrNN*fwp 431 432 srcL(e) = -1.0*tmp 433 srcR(e) = srcRat(e)*SclrNN 434 enddo 435c 436c.... One source term has the form (beta_i)*(phi,_i). Since 437c the convective term has (u_i)*(phi,_i), it is useful to treat 438c beta_i as a "correction" to the velocity. In calculating the 439c stabilization terms, the new "modified" velocity (u_i-beta_i) is 440c then used in place of the pure velocity for stabilization terms, 441c and the source term sneaks into the RHS and LHS. Here, the term 442c is given by beta_i=(cb_2)*(phi,_i)/(sigma) 443c 444 445 tmp = saCb2 * saSigmaInv 446 uMod(:,1) = u1(:) - tmp * gradS(:,1) 447 uMod(:,2) = u2(:) - tmp * gradS(:,2) 448 uMod(:,3) = u3(:) - tmp * gradS(:,3) 449 450 elseif(iRANS.eq.-2) then ! K-epsilon 451c.... compute the global gradient of u 452c 453 gradV = zero 454 do n = 1, nshl 455c 456c du_i/dx_j 457c 458c i j indices match array where V is the velocity (u in our notes) 459 gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * yl(:,n,2) 460 gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * yl(:,n,3) 461 gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * yl(:,n,4) 462c 463 gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * yl(:,n,2) 464 gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * yl(:,n,3) 465 gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * yl(:,n,4) 466c 467 gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * yl(:,n,2) 468 gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * yl(:,n,3) 469 gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * yl(:,n,4) 470c a j u a i 471c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in yl vector 472c 473 enddo 474 475 dwall = zero 476 do n = 1, nenl 477 dwall = dwall + shape_funct(:,n) * dwl(:,n) 478 enddo 479 480 kay(:)=zero 481 epsilon(:)=zero 482 do ii=1,npro 483 do jj = 1, nshl 484 kay(ii) = kay(ii) + shape_funct(ii,jj) * yl(ii,jj,6) 485 epsilon(ii) = epsilon(ii) 486 & + shape_funct(ii,jj) * yl(ii,jj,7) 487 enddo 488 enddo 489c no source term in the LHS yet 490 srcL(:)=zero 491 492 call elm3keps (kay, epsilon, 493 & dwall, gradV, 494 & srcRat, srcR, srcJac) 495 if(isclr.eq.1) srcL = srcJac(:,1) 496 if(isclr.eq.2) srcL = srcJac(:,4) 497 iadvdiff=0 ! scalar advection-diffusion flag 498 if(iadvdiff.eq.1)then 499 srcL(:)=zero 500 srcR(:)=zero 501 srcRat(:)=zero 502 endif 503c 504c.... No source terms with the form (beta_i)*(phi,_i) for K or E 505c 506 uMod(:,1) = u1(:) - zero 507 uMod(:,2) = u2(:) - zero 508 uMod(:,3) = u3(:) - zero 509 510 511 elseif (iLSet.ne.0) then 512 if (isclr.eq.1) then 513 514 srcR = zero 515 srcL = zero 516 517 elseif (isclr.eq.2) then !we are redistancing level-sets 518 519CAD If Sclr(:,1).gt.zero, result of sign_term function 1 520CAD If Sclr(:,1).eq.zero, result of sign_term function 0 521CAD If Sclr(:,1).lt.zero, result of sign_term function -1 522 523 sclr_ls = zero !zero out temp variable 524 525 do ii=1,npro 526 527 do jj = 1, nshl ! first find the value of levelset at point ii 528 529 sclr_ls(ii) = sclr_ls(ii) 530 & + shape_funct(ii,jj) * yl(ii,jj,6) 531 532 enddo 533 534 if (sclr_ls(ii) .gt. epsilon_ls)then 535 536 sign_levelset(ii) = one 537 538 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 539 540 sign_levelset(ii) = sclr_ls(ii)/epsilon_ls 541 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 542 543 elseif (sclr_ls(ii) .lt. epsilon_ls) then 544 545 sign_levelset(ii) = - one 546 547 endif 548 549 srcR(ii) = sign_levelset(ii) 550 551 enddo 552 553 srcL = zero 554 555cad The redistancing equation can be written in the following form 556cad 557cad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 558cad 559cad This is rewritten in the form 560cad 561cad d_{,t} + u * d_{,i} = sign(phi) 562cad 563 564CAD For the redistancing equation the "velocity" term is calculated as 565CAD follows 566 567 568 569 mytmp = sign_levelset / sqrt( gradS(:,1)*gradS(:,1) 570 & + gradS(:,2)*gradS(:,2) 571 & + gradS(:,3)*gradS(:,3)) 572 573 uMod(:,1) = mytmp * gradS(:,1) ! These are for the LHS 574 uMod(:,2) = mytmp * gradS(:,2) 575 uMod(:,3) = mytmp * gradS(:,3) 576 577 u1 = umod(:,1) ! These are for the RHS 578 u2 = umod(:,2) 579 u3 = umod(:,3) 580 581 582 endif ! close of scalar 2 of level set 583 584 else ! NOT turbulence and NOT level set so this is a simple 585 ! scalar. If your scalar equation has a source term 586 ! then add your own like the above but leave an unforced case 587 ! as an option like you see here 588 589 srcR = zero 590 srcL = zero 591 endif 592 593 return 594 end 595 596 597 598