1c$$$ subroutine e3DC (g1yi, g2yi, g3yi, 2c$$$ & giju, tauM, A0, raLS, 3c$$$ & rtLS, giju, DC, ri, 4c$$$ & rmi, stiff, A0DC) 5c$$$c 6c$$$c---------------------------------------------------------------------- 7c$$$c 8c$$$c This routine calculates the contribution of the Discontinuity- 9c$$$c Capturing operator to RHS and preconditioner. 10c$$$c 11c$$$c g1yi (nflow,npro) : grad-y in direction 1 12c$$$c g2yi (nflow,npro) : grad-y in direction 2 13c$$$c g3yi (nflow,npro) : grad-y in direction 3 14c$$$c A0 (nsymdf,npro) : A0 matrix (Symm. storage) 15c$$$c raLS (npro) : square of LS residual (A0inv norm) 16c$$$c rtLS (npro) : square of LS residual (Tau norm) 17c$$$c giju (6,npro) : metric matrix 18c$$$c DC (ngauss,npro) : discontinuity-capturing factor 19c$$$c intp : integration point number 20c$$$c 21c$$$c output: 22c$$$c ri (nflow*(nsd+1),npro) : partial residual 23c$$$c rmi (nflow*(nsd+1),npro) : partial modified residual 24c$$$c stiff (nsymdf,6,npro) : diffusivity matrix 25c$$$c DC (npro) : discontinuity-capturing factor 26c$$$c 27c$$$c 28c$$$c Zdenek Johan, Summer 1990. (Modified from e2dc.f) 29c$$$c Zdenek Johan, Winter 1991. (Recoded) 30c$$$c Zdenek Johan, Winter 1991. (Fortran 90) 31c$$$c---------------------------------------------------------------------- 32c$$$c 33c$$$ include "common.h" 34c$$$c 35c$$$ dimension g1yi(npro,nflow), g2yi(npro,nflow), 36c$$$ & g3yi(npro,nflow), A0(npro,5,5), 37c$$$ & raLS(npro), rtLS(npro), 38c$$$ & giju(npro,6), DC(npro,ngauss), 39c$$$ & ri(npro,nflow*(nsd+1)), rmi(npro,nflow*(nsd+1)), 40c$$$ & stiff(npro,3*nflow,3*nflow),itmp(npro) 41c$$$c 42c$$$ 43c$$$ dimension ggyi(npro,nflow), gAgyi(npro,15), 44c$$$ & gnorm(npro), A0gyi(npro,15), 45c$$$ & yiA0DCyj(npro,6), A0DC(npro,4) 46c$$$c 47c$$$c ... -----------------------> initialize <---------------------------- 48c$$$c 49c$$$ A0gyi = zero 50c$$$ gAgyi = zero 51c$$$ yiA0DCyj = zero 52c$$$ DC = zero 53c$$$c.... -----------------------> global gradient <---------------------- 54c$$$c 55c$$$c.... calculate (A0 y_,j) --> A0gyi 56c$$$c 57c$$$c A0 Y_{,1} 58c$$$c 59c$$$ A0gyi( :,1) = A0(:,1,1)*g1yi(:,1) 60c$$$ & + A0(:,1,2)*g1yi(:,2) 61c$$$ & + A0(:,1,3)*g1yi(:,3) 62c$$$ & + A0(:,1,4)*g1yi(:,4) 63c$$$ & + A0(:,1,5)*g1yi(:,5) 64c$$$ A0gyi( :,2) = A0(:,2,1)*g1yi(:,1) 65c$$$ & + A0(:,2,2)*g1yi(:,2) 66c$$$ & + A0(:,2,3)*g1yi(:,3) 67c$$$ & + A0(:,2,4)*g1yi(:,4) 68c$$$ & + A0(:,2,5)*g1yi(:,5) 69c$$$ A0gyi( :,3) = A0(:,3,1)*g1yi(:,1) 70c$$$ & + A0(:,3,2)*g1yi(:,2) 71c$$$ & + A0(:,3,3)*g1yi(:,3) 72c$$$ & + A0(:,3,4)*g1yi(:,4) 73c$$$ & + A0(:,3,5)*g1yi(:,5) 74c$$$ A0gyi( :,4) = A0(:,4,1)*g1yi(:,1) 75c$$$ & + A0(:,4,2)*g1yi(:,2) 76c$$$ & + A0(:,4,3)*g1yi(:,3) 77c$$$ & + A0(:,4,4)*g1yi(:,4) 78c$$$ & + A0(:,4,5)*g1yi(:,5) 79c$$$ A0gyi( :,5) = A0(:,5,1)*g1yi(:,1) 80c$$$ & + A0(:,5,2)*g1yi(:,2) 81c$$$ & + A0(:,5,3)*g1yi(:,3) 82c$$$ & + A0(:,5,4)*g1yi(:,4) 83c$$$ & + A0(:,5,5)*g1yi(:,5) 84c$$$c 85c$$$c A0 Y_{,2} 86c$$$c 87c$$$ A0gyi( :,6) = A0(:,1,1)*g2yi(:,1) 88c$$$ & + A0(:,1,2)*g2yi(:,2) 89c$$$ & + A0(:,1,3)*g2yi(:,3) 90c$$$ & + A0(:,1,4)*g2yi(:,4) 91c$$$ & + A0(:,1,5)*g2yi(:,5) 92c$$$ A0gyi( :,7) = A0(:,2,1)*g2yi(:,1) 93c$$$ & + A0(:,2,2)*g2yi(:,2) 94c$$$ & + A0(:,2,3)*g2yi(:,3) 95c$$$ & + A0(:,2,4)*g2yi(:,4) 96c$$$ & + A0(:,2,5)*g2yi(:,5) 97c$$$ A0gyi( :,8) = A0(:,3,1)*g2yi(:,1) 98c$$$ & + A0(:,3,2)*g2yi(:,2) 99c$$$ & + A0(:,3,3)*g2yi(:,3) 100c$$$ & + A0(:,3,4)*g2yi(:,4) 101c$$$ & + A0(:,3,5)*g2yi(:,5) 102c$$$ A0gyi( :,9) = A0(:,4,1)*g2yi(:,1) 103c$$$ & + A0(:,4,2)*g2yi(:,2) 104c$$$ & + A0(:,4,3)*g2yi(:,3) 105c$$$ & + A0(:,4,4)*g2yi(:,4) 106c$$$ & + A0(:,4,5)*g2yi(:,5) 107c$$$ A0gyi(:,10) = A0(:,5,1)*g2yi(:,1) 108c$$$ & + A0(:,5,2)*g2yi(:,2) 109c$$$ & + A0(:,5,3)*g2yi(:,3) 110c$$$ & + A0(:,5,4)*g2yi(:,4) 111c$$$ & + A0(:,5,5)*g2yi(:,5) 112c$$$c 113c$$$c A0 Y_{,3} 114c$$$c 115c$$$ A0gyi(:,11) = A0(:,1,1)*g3yi(:,1) 116c$$$ & + A0(:,1,2)*g3yi(:,2) 117c$$$ & + A0(:,1,3)*g3yi(:,3) 118c$$$ & + A0(:,1,4)*g3yi(:,4) 119c$$$ & + A0(:,1,5)*g3yi(:,5) 120c$$$ A0gyi(:,12) = A0(:,2,1)*g3yi(:,1) 121c$$$ & + A0(:,2,2)*g3yi(:,2) 122c$$$ & + A0(:,2,3)*g3yi(:,3) 123c$$$ & + A0(:,2,4)*g3yi(:,4) 124c$$$ & + A0(:,2,5)*g3yi(:,5) 125c$$$ A0gyi(:,13) = A0(:,3,1)*g3yi(:,1) 126c$$$ & + A0(:,3,2)*g3yi(:,2) 127c$$$ & + A0(:,3,3)*g3yi(:,3) 128c$$$ & + A0(:,3,4)*g3yi(:,4) 129c$$$ & + A0(:,3,5)*g3yi(:,5) 130c$$$ A0gyi(:,14) = A0(:,4,1)*g3yi(:,1) 131c$$$ & + A0(:,4,2)*g3yi(:,2) 132c$$$ & + A0(:,4,3)*g3yi(:,3) 133c$$$ & + A0(:,4,4)*g3yi(:,4) 134c$$$ & + A0(:,4,5)*g3yi(:,5) 135c$$$ A0gyi(:,15) = A0(:,5,1)*g3yi(:,1) 136c$$$ & + A0(:,5,2)*g3yi(:,2) 137c$$$ & + A0(:,5,3)*g3yi(:,3) 138c$$$ & + A0(:,5,4)*g3yi(:,4) 139c$$$ & + A0(:,5,5)*g3yi(:,5) 140c$$$c 141c$$$c.... calculate (giju A0 y_,j) --> gAgyi 142c$$$c 143c$$$ 144c$$$ gAgyi( :,1) = giju(:,1)*A0gyi( :,1) 145c$$$ & + giju(:,4)*A0gyi( :,6) 146c$$$ & + giju(:,5)*A0gyi(:,11) 147c$$$ 148c$$$ gAgyi( :,2) = giju(:,1)*A0gyi( :,2) 149c$$$ & + giju(:,4)*A0gyi( :,7) 150c$$$ & + giju(:,5)*A0gyi(:,12) 151c$$$ 152c$$$ gAgyi( :,3) = giju(:,1)*A0gyi( :,3) 153c$$$ & + giju(:,4)*A0gyi( :,8) 154c$$$ & + giju(:,5)*A0gyi(:,13) 155c$$$ 156c$$$ gAgyi( :,4) = giju(:,1)*A0gyi( :,4) 157c$$$ & + giju(:,4)*A0gyi( :,9) 158c$$$ & + giju(:,5)*A0gyi(:,14) 159c$$$ 160c$$$ gAgyi( :,5) = giju(:,1)*A0gyi( :,5) 161c$$$ & + giju(:,4)*A0gyi(:,10) 162c$$$ & + giju(:,5)*A0gyi(:,15) 163c$$$ 164c$$$ gAgyi( :,6) = giju(:,4)*A0gyi( :,1) 165c$$$ & + giju(:,2)*A0gyi( :,6) 166c$$$ & + giju(:,6)*A0gyi(:,11) 167c$$$ 168c$$$ gAgyi( :,7) = giju(:,4)*A0gyi( :,2) 169c$$$ & + giju(:,2)*A0gyi( :,7) 170c$$$ & + giju(:,6)*A0gyi(:,12) 171c$$$ 172c$$$ gAgyi( :,8) = giju(:,4)*A0gyi( :,3) 173c$$$ & + giju(:,2)*A0gyi( :,8) 174c$$$ & + giju(:,6)*A0gyi(:,13) 175c$$$ 176c$$$ gAgyi( :,9) = giju(:,4)*A0gyi( :,4) 177c$$$ & + giju(:,2)*A0gyi( :,9) 178c$$$ & + giju(:,6)*A0gyi(:,14) 179c$$$ 180c$$$ gAgyi(:,10) = giju(:,4)*A0gyi( :,5) 181c$$$ & + giju(:,2)*A0gyi(:,10) 182c$$$ & + giju(:,6)*A0gyi(:,15) 183c$$$ 184c$$$ gAgyi(:,11) = giju(:,5)*A0gyi( :,1) 185c$$$ & + giju(:,6)*A0gyi( :,6) 186c$$$ & + giju(:,3)*A0gyi(:,11) 187c$$$ 188c$$$ gAgyi(:,12) = giju(:,5)*A0gyi( :,2) 189c$$$ & + giju(:,6)*A0gyi( :,7) 190c$$$ & + giju(:,3)*A0gyi(:,12) 191c$$$ 192c$$$ gAgyi(:,13) = giju(:,5)*A0gyi( :,3) 193c$$$ & + giju(:,6)*A0gyi( :,8) 194c$$$ & + giju(:,3)*A0gyi(:,13) 195c$$$ 196c$$$ gAgyi(:,14) = giju(:,5)*A0gyi( :,4) 197c$$$ & + giju(:,6)*A0gyi( :,9) 198c$$$ & + giju(:,3)*A0gyi(:,14) 199c$$$ 200c$$$ gAgyi(:,15) = giju(:,5)*A0gyi( :,5) 201c$$$ & + giju(:,6)*A0gyi(:,10) 202c$$$ & + giju(:,3)*A0gyi(:,15) 203c$$$c 204c$$$c... the denominator term of the DC factor 205c$$$c... evaluation of the term Y,i.A0DC Y,j 206c$$$c 207c$$$ yiA0DCyj(:,1) = A0DC(:,1)*g1yi(:,1)**2 208c$$$ & + two*g1yi(:,1)*A0DC(:,2)*g1yi(:,5) 209c$$$ & + A0DC(:,3)*g1yi(:,2)**2 210c$$$ & + A0DC(:,3)*g1yi(:,3)**2 211c$$$ & + A0DC(:,3)*g1yi(:,4)**2 212c$$$ & + A0DC(:,4)*g1yi(:,5)**2 213c$$$ 214c$$$ yiA0DCyj(:,2) = A0DC(:,1)*g2yi(:,1)**2 215c$$$ & + two*g2yi(:,1)*A0DC(:,2)*g2yi(:,5) 216c$$$ & + A0DC(:,3)*g2yi(:,2)**2 217c$$$ & + A0DC(:,3)*g2yi(:,3)**2 218c$$$ & + A0DC(:,3)*g2yi(:,4)**2 219c$$$ & + A0DC(:,4)*g2yi(:,5)**2 220c$$$ 221c$$$ yiA0DCyj(:,3) = A0DC(:,1)*g3yi(:,1)**2 222c$$$ & + two*g3yi(:,1)*A0DC(:,2)*g3yi(:,5) 223c$$$ & + A0DC(:,3)*g3yi(:,2)**2 224c$$$ & + A0DC(:,3)*g3yi(:,3)**2 225c$$$ & + A0DC(:,3)*g3yi(:,4)**2 226c$$$ & + A0DC(:,4)*g3yi(:,5)**2 227c$$$ 228c$$$ yiA0DCyj(:,4) = g1yi(:,1)*A0DC(:,1)*g2yi(:,1) 229c$$$ & + g1yi(:,1)*A0DC(:,2)*g2yi(:,5) 230c$$$ & + g1yi(:,2)*A0DC(:,3)*g2yi(:,2) 231c$$$ & + g1yi(:,3)*A0DC(:,3)*g2yi(:,3) 232c$$$ & + g1yi(:,4)*A0DC(:,3)*g2yi(:,4) 233c$$$ & + g1yi(:,5)*A0DC(:,2)*g2yi(:,1) 234c$$$ & + g1yi(:,5)*A0DC(:,4)*g2yi(:,5) 235c$$$ 236c$$$ yiA0DCyj(:,5) = g1yi(:,1)*A0DC(:,1)*g3yi(:,1) 237c$$$ & + g1yi(:,1)*A0DC(:,2)*g3yi(:,5) 238c$$$ & + g1yi(:,2)*A0DC(:,3)*g3yi(:,2) 239c$$$ & + g1yi(:,3)*A0DC(:,3)*g3yi(:,3) 240c$$$ & + g1yi(:,4)*A0DC(:,3)*g3yi(:,4) 241c$$$ & + g1yi(:,5)*A0DC(:,2)*g3yi(:,1) 242c$$$ & + g1yi(:,5)*A0DC(:,4)*g3yi(:,5) 243c$$$ 244c$$$ yiA0DCyj(:,6) = g2yi(:,1)*A0DC(:,1)*g3yi(:,1) 245c$$$ & + g2yi(:,1)*A0DC(:,2)*g3yi(:,5) 246c$$$ & + g2yi(:,2)*A0DC(:,3)*g3yi(:,2) 247c$$$ & + g2yi(:,3)*A0DC(:,3)*g3yi(:,3) 248c$$$ & + g2yi(:,4)*A0DC(:,3)*g3yi(:,4) 249c$$$ & + g2yi(:,5)*A0DC(:,2)*g3yi(:,1) 250c$$$ & + g2yi(:,5)*A0DC(:,4)*g3yi(:,5) 251c$$$c 252c$$$c.... -------------------------> DC factor <-------------------------- 253c$$$c 254c$$$ if ((ires .ne. 2) .or. (Jactyp .eq. 1)) then 255c$$$c 256c$$$c.... calculate 2-norm of Grad-local-V with respect to A0 257c$$$c 258c$$$c.... DC-mallet 259c$$$c 260c$$$ if (iDC .eq. 1) then 261c$$$c 262c$$$ fact = one 263c$$$ if (ipord .eq. 2) fact = 0.9 264c$$$ if (ipord .eq. 3) fact = 0.75 265c$$$ 266c$$$c 267c$$$ gnorm = one / ( 268c$$$ & giju(:,1)*yiA0DCyj(:,1) 269c$$$ & + two*giju(:,4)*yiA0DCyj(:,4) 270c$$$ & + two*giju(:,5)*yiA0DCyj(:,5) 271c$$$ & + giju(:,2)*yiA0DCyj(:,2) 272c$$$ & + two*giju(:,6)*yiA0DCyj(:,6) 273c$$$ & + giju(:,3)*yiA0DCyj(:,3) 274c$$$ & + epsM ) 275c$$$c 276c$$$ DC(:,intp)=dim((fact*sqrt(raLS*gnorm)),(rtLS*gnorm)) 277c$$$c 278c$$$c.... flop count 279c$$$c 280c$$$ flops = flops + 46*npro 281c$$$c 282c$$$ endif 283c$$$c 284c$$$c.... DC-quadratic 285c$$$c 286c$$$ if (iDC .eq. 2) then 287c$$$c 288c$$$ gnorm = one / ( 289c$$$ & giju(:,1)*yiA0DCyj(:,1) 290c$$$ & + two*giju(:,4)*yiA0DCyj(:,4) 291c$$$ & + two*giju(:,5)*yiA0DCyj(:,5) 292c$$$ & + giju(:,2)*yiA0DCyj(:,2) 293c$$$ & + two*giju(:,6)*yiA0DCyj(:,6) 294c$$$ & + giju(:,3)*yiA0DCyj(:,3) 295c$$$ & + epsM ) 296c$$$ 297c$$$c 298c$$$ DC(:,intp) = two * rtLS * gnorm 299c$$$c 300c$$$c.... flop count 301c$$$c 302c$$$ flops = flops + 36*npro 303c$$$c 304c$$$ endif 305c$$$c 306c$$$c.... DC-min 307c$$$c 308c$$$ if (iDC .eq. 3) then 309c$$$c 310c$$$ fact = one 311c$$$ if (ipord .eq. 2) fact = pt5 312c$$$c 313c$$$ gnorm = one / ( 314c$$$ & giju(:,1)*yiA0DCyj(:,1) 315c$$$ & + two*giju(:,4)*yiA0DCyj(:,4) 316c$$$ & + two*giju(:,5)*yiA0DCyj(:,5) 317c$$$ & + giju(:,2)*yiA0DCyj(:,2) 318c$$$ & + two*giju(:,6)*yiA0DCyj(:,6) 319c$$$ & + giju(:,3)*yiA0DCyj(:,3) 320c$$$ & + epsM ) 321c$$$ 322c$$$c 323c$$$ DC(:,intp) = min( dim(fact * sqrt(raLS * gnorm), 324c$$$ & rtLS * gnorm), two * rtLS * gnorm ) 325c$$$c 326c$$$c.... flop count 327c$$$c 328c$$$ flops = flops + 48*npro 329c$$$c 330c$$$ endif 331c$$$c 332c$$$ endif 333c$$$c 334c$$$c.... ----------------------------> RHS <---------------------------- 335c$$$c 336c$$$c.... add the contribution of DC to ri and/or rmi 337c$$$c 338c$$$c.... ires = 1 or 3 339c$$$c 340c$$$ if ((ires .eq. 1) .or. (ires .eq. 3)) then 341c$$$c 342c$$$ ri ( :,1) = ri ( :,1) + DC(:,intp) * gAgyi( :,1) 343c$$$ rmi( :,1) = rmi( :,1) + DC(:,intp) * gAgyi( :,1) 344c$$$ ri ( :,2) = ri ( :,2) + DC(:,intp) * gAgyi( :,2) 345c$$$ rmi( :,2) = rmi( :,2) + DC(:,intp) * gAgyi( :,2) 346c$$$ ri ( :,3) = ri ( :,3) + DC(:,intp) * gAgyi( :,3) 347c$$$ rmi( :,3) = rmi( :,3) + DC(:,intp) * gAgyi( :,3) 348c$$$ ri ( :,4) = ri ( :,4) + DC(:,intp) * gAgyi( :,4) 349c$$$ rmi( :,4) = rmi( :,4) + DC(:,intp) * gAgyi( :,4) 350c$$$ ri ( :,5) = ri ( :,5) + DC(:,intp) * gAgyi( :,5) 351c$$$ rmi( :,5) = rmi( :,5) + DC(:,intp) * gAgyi( :,5) 352c$$$c 353c$$$ ri ( :,6) = ri ( :,6) + DC(:,intp) * gAgyi( :,6) 354c$$$ rmi( :,6) = rmi( :,6) + DC(:,intp) * gAgyi( :,6) 355c$$$ ri ( :,7) = ri ( :,7) + DC(:,intp) * gAgyi( :,7) 356c$$$ rmi( :,7) = rmi( :,7) + DC(:,intp) * gAgyi( :,7) 357c$$$ ri ( :,8) = ri ( :,8) + DC(:,intp) * gAgyi( :,8) 358c$$$ rmi( :,8) = rmi( :,8) + DC(:,intp) * gAgyi( :,8) 359c$$$ ri ( :,9) = ri ( :,9) + DC(:,intp) * gAgyi( :,9) 360c$$$ rmi( :,9) = rmi( :,9) + DC(:,intp) * gAgyi( :,9) 361c$$$ ri (:,10) = ri (:,10) + DC(:,intp) * gAgyi(:,10) 362c$$$ rmi(:,10) = rmi(:,10) + DC(:,intp) * gAgyi(:,10) 363c$$$c 364c$$$ ri (:,11) = ri (:,11) + DC(:,intp) * gAgyi(:,11) 365c$$$ rmi(:,11) = rmi(:,12) + DC(:,intp) * gAgyi(:,12) 366c$$$ ri (:,12) = ri (:,12) + DC(:,intp) * gAgyi(:,12) 367c$$$ rmi(:,12) = rmi(:,12) + DC(:,intp) * gAgyi(:,12) 368c$$$ ri (:,13) = ri (:,13) + DC(:,intp) * gAgyi(:,13) 369c$$$ rmi(:,13) = rmi(:,13) + DC(:,intp) * gAgyi(:,13) 370c$$$ ri (:,14) = ri (:,14) + DC(:,intp) * gAgyi(:,14) 371c$$$ rmi(:,14) = rmi(:,14) + DC(:,intp) * gAgyi(:,14) 372c$$$ ri (:,15) = ri (:,15) + DC(:,intp) * gAgyi(:,15) 373c$$$ rmi(:,15) = rmi(:,15) + DC(:,intp) * gAgyi(:,15) 374c$$$c 375c$$$ flops = flops + 45*npro 376c$$$c 377c$$$ endif 378c$$$c 379c$$$c.... ires = 2 380c$$$c 381c$$$ if (ires .eq. 2) then 382c$$$c 383c$$$ rmi( :,1) = rmi( :,1) + DC(:,intp) * gAgyi( :,1) 384c$$$ rmi( :,2) = rmi( :,2) + DC(:,intp) * gAgyi( :,2) 385c$$$ rmi( :,3) = rmi( :,3) + DC(:,intp) * gAgyi( :,3) 386c$$$ rmi( :,4) = rmi( :,4) + DC(:,intp) * gAgyi( :,4) 387c$$$ rmi( :,5) = rmi( :,5) + DC(:,intp) * gAgyi( :,5) 388c$$$c 389c$$$ rmi( :,6) = rmi( :,6) + DC(:,intp) * gAgyi( :,6) 390c$$$ rmi( :,7) = rmi( :,7) + DC(:,intp) * gAgyi( :,7) 391c$$$ rmi( :,8) = rmi( :,8) + DC(:,intp) * gAgyi( :,8) 392c$$$ rmi( :,9) = rmi( :,9) + DC(:,intp) * gAgyi( :,9) 393c$$$ rmi(:,10) = rmi(:,10) + DC(:,intp) * gAgyi(:,10) 394c$$$c 395c$$$ rmi(:,11) = rmi(:,11) + DC(:,intp) * gAgyi(:,11) 396c$$$ rmi(:,12) = rmi(:,12) + DC(:,intp) * gAgyi(:,12) 397c$$$ rmi(:,13) = rmi(:,13) + DC(:,intp) * gAgyi(:,13) 398c$$$ rmi(:,14) = rmi(:,14) + DC(:,intp) * gAgyi(:,14) 399c$$$ rmi(:,15) = rmi(:,15) + DC(:,intp) * gAgyi(:,15) 400c$$$c 401c$$$ flops = flops + 30*npro 402c$$$c 403c$$$ endif 404c$$$c 405c$$$c.... -------------------------> Stiffness <-------------------------- 406c$$$c 407c$$$c.... add the contribution of DC to stiff 408c$$$c 409c$$$ if (iprec .eq. 1) then 410c$$$ nflow2=two*nflow 411c$$$ do j = 1, nflow 412c$$$ do i = 1, nflow 413c$$$ itmp(:)=A0(:,i,j)*DC(:,intp) 414c$$$c 415c$$$c.... add (DC g^1 A0) to stiff [1,1] 416c$$$c 417c$$$ stiff(:,i,j) = stiff(:,i,j) 418c$$$ & + itmp(:)*giju(:,1) 419c$$$c 420c$$$c.... add (DC g^1 A0) to stiff [1,2] 421c$$$c 422c$$$ 423c$$$ stiff(:,i,j+nflow) = stiff(:,i,j+nflow) 424c$$$ & + itmp(:)*giju(:,4) 425c$$$c 426c$$$c.... add (DC g^1 A0) to stiff [1,3] 427c$$$c 428c$$$ 429c$$$ stiff(:,i,j+nflow2) = stiff(:,i,j+nflow2) 430c$$$ & + itmp(:)*giju(:,5) 431c$$$ 432c$$$c.... add (DC g^1 A0) to stiff [2,1] (similarly below) 433c$$$c 434c$$$ 435c$$$ stiff(:,i+nflow,j) = stiff(:,i+nflow,j) 436c$$$ & + itmp(:)*giju(:,4) 437c$$$ 438c$$$ stiff(:,i+nflow,j+nflow) = stiff(:,i+nflow,j+nflow) 439c$$$ & + itmp(:)*giju(:,2) 440c$$$ 441c$$$ stiff(:,i+nflow,j+nflow2) = stiff(:,i+nflow,j+nflow2) 442c$$$ & + itmp(:)*giju(:,6) 443c$$$ 444c$$$ stiff(:,i+nflow2,j) = stiff(:,i+nflow2,j) 445c$$$ & + itmp(:)*giju(:,5) 446c$$$ 447c$$$ stiff(:,i+nflow2,j+nflow) = stiff(:,i+nflow2,j+nflow) 448c$$$ & + itmp(:)*giju(:,6) 449c$$$ 450c$$$ stiff(:,i+nflow2,j+nflow2) = stiff(:,i+nflow2,j+nflow2) 451c$$$ & + itmp(:)*giju(:,3) 452c$$$ enddo 453c$$$ enddo 454c$$$c 455c$$$c.... flop count 456c$$$c 457c$$$ flops = flops + 210*npro 458c$$$c 459c$$$c.... end of stiffness 460c$$$c 461c$$$ endif 462c$$$c 463c$$$c.... return 464c$$$c 465c$$$ return 466c$$$ end 467c$$$c 468c 469 subroutine e3dcSclr ( gradS, giju, gGradS, 470 & rLS, tauS, srcR, 471 & dcFct) 472c 473c 474c---------------------------------------------------------------------- 475c 476c This routine calculates the contribution of the Discontinuity- 477c Capturing operator to RHS and preconditioner for the scalar solve. 478c 479c g1yti (nflow,npro) : grad-y in direction 1 480c g2yti (nflow,npro) : grad-y in direction 2 481c g3yti (nflow,npro) : grad-y in direction 3 482c A0 (nsymdf,npro) : A0 matrix (Symm. storage) 483c raLS (npro) : square of LS residual (A0inv norm) 484c rtLS (npro) : square of LS residual (Tau norm) 485c giju (6,npro) : metric matrix 486c DC (ngauss,npro) : discontinuity-capturing factor 487c intp : integration point number 488c 489c output: 490c ri (nflow*(nsd+1),npro) : partial residual 491c rmi (nflow*(nsd+1),npro) : partial modified residual 492c stiff (nsymdf,6,npro) : diffusivity matrix 493c DC (npro) : discontinuity-capturing factor 494c 495c 496c Zdenek Johan, Summer 1990. (Modified from e2dc.f) 497c Zdenek Johan, Winter 1991. (Recoded) 498c Zdenek Johan, Winter 1991. (Fortran 90) 499c---------------------------------------------------------------------- 500c 501 include "common.h" 502c 503 dimension gradS(npro,nsd), gGradS(npro,nsd), 504 & rLS(npro), tauS(npro), 505 & giju(npro,6), dcFct(npro), 506 & srcR(npro) 507c 508c.... Form GijUp gradS and gradS . GijUp gradS (store in dcFct) 509c 510 511 gGradS(:,1) = GijU(:,1) * gradS(:,1) 512 1 + GijU(:,4) * gradS(:,2) 513 2 + GijU(:,6) * gradS(:,3) 514 gGradS(:,2) = GijU(:,4) * gradS(:,1) 515 1 + GijU(:,2) * gradS(:,2) 516 2 + GijU(:,5) * gradS(:,3) 517 gGradS(:,3) = GijU(:,6) * gradS(:,1) 518 1 + GijU(:,5) * gradS(:,2) 519 2 + GijU(:,3) * gradS(:,3) 520c 521 dcFct(:) = gradS(:,1) * gGradS(:,1) 522 1 + gradS(:,2) * gGradS(:,2) 523 2 + gradS(:,3) * gGradS(:,3) 524 3 + epsM 525 526 dcFct(:) = 1.0/ dcFct(:) 527c 528c.... Form pdeRes 2-norm / gradT 2-norm 529c 530 531 dcFct = dcFct * (rLS - srcR) ** 2 532c 533c.... -------------------------> DC factor <------------------------ 534c 535c.... DC-mallet 536c 537 if (idcsclr(1) .eq. 1) then 538c 539 fact = one 540 if (ipord .eq. 2) fact = 0.9 541 if (ipord .eq. 3) fact = 0.75 542c 543c$$$ dcFct(:)=dim((fact*sqrt(dcFct(:))),(tauS(:)*dcFct(:))) !not work 544 !with all compilers 545 dcFct(:)=max(zero,(fact*sqrt(dcFct(:)))-(tauS(:)*dcFct(:))) 546c 547 endif 548c 549c 550c.... DC-quadratic 551c 552 if (idcsclr(1) .eq. 2) then 553c 554 dcFct(:) = two * tauS(:) * dcFct(:) 555c 556 endif 557c 558c.... DC-min 559c 560 if (idcsclr(1) .eq. 3) then 561c 562 fact = one 563 if (ipord .eq. 2) fact = 0.9 564c 565 dcFct(:) = min( max(zero, (fact * sqrt(dcFct(:)) - 566 & tauS(:)*dcFct(:)) ), two * tauS(:) * dcFct(:)) 567c 568 endif 569c 570c.... Scale the gGradT for residual formation 571c 572 gGradS(:,1) = dcFct(:) * gGradS(:,1) 573 gGradS(:,2) = dcFct(:) * gGradS(:,2) 574 gGradS(:,3) = dcFct(:) * gGradS(:,3) 575 576 577 578 return 579 end 580c 581