1 subroutine e3DC (g1yi, g2yi, g3yi, A0, raLS, 2 & rtLS, giju, DC, ri, 3 & rmi, stiff, A0DC) 4c 5c---------------------------------------------------------------------- 6c 7c This routine calculates the contribution of the Discontinuity- 8c Capturing operator to RHS and preconditioner. 9c 10c g1yi (nflow,npro) : grad-y in direction 1 11c g2yi (nflow,npro) : grad-y in direction 2 12c g3yi (nflow,npro) : grad-y in direction 3 13c A0 (nsymdf,npro) : A0 matrix (Symm. storage) 14c raLS (npro) : square of LS residual (A0inv norm) 15c rtLS (npro) : square of LS residual (Tau norm) 16c giju (6,npro) : metric matrix 17c DC (ngauss,npro) : discontinuity-capturing factor 18c intp : integration point number 19c 20c output: 21c ri (nflow*(nsd+1),npro) : partial residual 22c rmi (nflow*(nsd+1),npro) : partial modified residual 23c stiff (nsymdf,6,npro) : diffusivity matrix 24c DC (npro) : discontinuity-capturing factor 25c 26c 27c Zdenek Johan, Summer 1990. (Modified from e2dc.f) 28c Zdenek Johan, Winter 1991. (Recoded) 29c Zdenek Johan, Winter 1991. (Fortran 90) 30c---------------------------------------------------------------------- 31c 32 include "common.h" 33c 34 dimension g1yi(npro,nflow), g2yi(npro,nflow), 35 & g3yi(npro,nflow), A0(npro,5,5), 36 & raLS(npro), rtLS(npro), 37 & giju(npro,6), DC(npro,ngauss), 38 & ri(npro,nflow*(nsd+1)), rmi(npro,nflow*(nsd+1)), 39 & stiff(npro,3*nflow,3*nflow),dtmp(npro) 40c 41 42 dimension ggyi(npro,nflow), gAgyi(npro,15), 43 & gnorm(npro), A0gyi(npro,15), 44 & yiA0DCyj(npro,6), A0DC(npro,4) 45c 46c ... -----------------------> initialize <---------------------------- 47c 48 A0gyi = zero 49 gAgyi = zero 50 yiA0DCyj = zero 51 DC = zero 52c.... -----------------------> global gradient <---------------------- 53c 54c.... calculate (A0 y_,j) --> A0gyi 55c 56c A0 Y_{,1} 57c 58 A0gyi( :,1) = A0(:,1,1)*g1yi(:,1) 59 & + A0(:,1,2)*g1yi(:,2) 60 & + A0(:,1,3)*g1yi(:,3) 61 & + A0(:,1,4)*g1yi(:,4) 62 & + A0(:,1,5)*g1yi(:,5) 63 A0gyi( :,2) = A0(:,2,1)*g1yi(:,1) 64 & + A0(:,2,2)*g1yi(:,2) 65 & + A0(:,2,3)*g1yi(:,3) 66 & + A0(:,2,4)*g1yi(:,4) 67 & + A0(:,2,5)*g1yi(:,5) 68 A0gyi( :,3) = A0(:,3,1)*g1yi(:,1) 69 & + A0(:,3,2)*g1yi(:,2) 70 & + A0(:,3,3)*g1yi(:,3) 71 & + A0(:,3,4)*g1yi(:,4) 72 & + A0(:,3,5)*g1yi(:,5) 73 A0gyi( :,4) = A0(:,4,1)*g1yi(:,1) 74 & + A0(:,4,2)*g1yi(:,2) 75 & + A0(:,4,3)*g1yi(:,3) 76 & + A0(:,4,4)*g1yi(:,4) 77 & + A0(:,4,5)*g1yi(:,5) 78 A0gyi( :,5) = A0(:,5,1)*g1yi(:,1) 79 & + A0(:,5,2)*g1yi(:,2) 80 & + A0(:,5,3)*g1yi(:,3) 81 & + A0(:,5,4)*g1yi(:,4) 82 & + A0(:,5,5)*g1yi(:,5) 83c 84c A0 Y_{,2} 85c 86 A0gyi( :,6) = A0(:,1,1)*g2yi(:,1) 87 & + A0(:,1,2)*g2yi(:,2) 88 & + A0(:,1,3)*g2yi(:,3) 89 & + A0(:,1,4)*g2yi(:,4) 90 & + A0(:,1,5)*g2yi(:,5) 91 A0gyi( :,7) = A0(:,2,1)*g2yi(:,1) 92 & + A0(:,2,2)*g2yi(:,2) 93 & + A0(:,2,3)*g2yi(:,3) 94 & + A0(:,2,4)*g2yi(:,4) 95 & + A0(:,2,5)*g2yi(:,5) 96 A0gyi( :,8) = A0(:,3,1)*g2yi(:,1) 97 & + A0(:,3,2)*g2yi(:,2) 98 & + A0(:,3,3)*g2yi(:,3) 99 & + A0(:,3,4)*g2yi(:,4) 100 & + A0(:,3,5)*g2yi(:,5) 101 A0gyi( :,9) = A0(:,4,1)*g2yi(:,1) 102 & + A0(:,4,2)*g2yi(:,2) 103 & + A0(:,4,3)*g2yi(:,3) 104 & + A0(:,4,4)*g2yi(:,4) 105 & + A0(:,4,5)*g2yi(:,5) 106 A0gyi(:,10) = A0(:,5,1)*g2yi(:,1) 107 & + A0(:,5,2)*g2yi(:,2) 108 & + A0(:,5,3)*g2yi(:,3) 109 & + A0(:,5,4)*g2yi(:,4) 110 & + A0(:,5,5)*g2yi(:,5) 111c 112c A0 Y_{,3} 113c 114 A0gyi(:,11) = A0(:,1,1)*g3yi(:,1) 115 & + A0(:,1,2)*g3yi(:,2) 116 & + A0(:,1,3)*g3yi(:,3) 117 & + A0(:,1,4)*g3yi(:,4) 118 & + A0(:,1,5)*g3yi(:,5) 119 A0gyi(:,12) = A0(:,2,1)*g3yi(:,1) 120 & + A0(:,2,2)*g3yi(:,2) 121 & + A0(:,2,3)*g3yi(:,3) 122 & + A0(:,2,4)*g3yi(:,4) 123 & + A0(:,2,5)*g3yi(:,5) 124 A0gyi(:,13) = A0(:,3,1)*g3yi(:,1) 125 & + A0(:,3,2)*g3yi(:,2) 126 & + A0(:,3,3)*g3yi(:,3) 127 & + A0(:,3,4)*g3yi(:,4) 128 & + A0(:,3,5)*g3yi(:,5) 129 A0gyi(:,14) = A0(:,4,1)*g3yi(:,1) 130 & + A0(:,4,2)*g3yi(:,2) 131 & + A0(:,4,3)*g3yi(:,3) 132 & + A0(:,4,4)*g3yi(:,4) 133 & + A0(:,4,5)*g3yi(:,5) 134 A0gyi(:,15) = A0(:,5,1)*g3yi(:,1) 135 & + A0(:,5,2)*g3yi(:,2) 136 & + A0(:,5,3)*g3yi(:,3) 137 & + A0(:,5,4)*g3yi(:,4) 138 & + A0(:,5,5)*g3yi(:,5) 139c 140c.... calculate (giju A0 y_,j) --> gAgyi 141c 142 143 gAgyi( :,1) = giju(:,1)*A0gyi( :,1) 144 & + giju(:,4)*A0gyi( :,6) 145 & + giju(:,5)*A0gyi(:,11) 146 147 gAgyi( :,2) = giju(:,1)*A0gyi( :,2) 148 & + giju(:,4)*A0gyi( :,7) 149 & + giju(:,5)*A0gyi(:,12) 150 151 gAgyi( :,3) = giju(:,1)*A0gyi( :,3) 152 & + giju(:,4)*A0gyi( :,8) 153 & + giju(:,5)*A0gyi(:,13) 154 155 gAgyi( :,4) = giju(:,1)*A0gyi( :,4) 156 & + giju(:,4)*A0gyi( :,9) 157 & + giju(:,5)*A0gyi(:,14) 158 159 gAgyi( :,5) = giju(:,1)*A0gyi( :,5) 160 & + giju(:,4)*A0gyi(:,10) 161 & + giju(:,5)*A0gyi(:,15) 162 163 gAgyi( :,6) = giju(:,4)*A0gyi( :,1) 164 & + giju(:,2)*A0gyi( :,6) 165 & + giju(:,6)*A0gyi(:,11) 166 167 gAgyi( :,7) = giju(:,4)*A0gyi( :,2) 168 & + giju(:,2)*A0gyi( :,7) 169 & + giju(:,6)*A0gyi(:,12) 170 171 gAgyi( :,8) = giju(:,4)*A0gyi( :,3) 172 & + giju(:,2)*A0gyi( :,8) 173 & + giju(:,6)*A0gyi(:,13) 174 175 gAgyi( :,9) = giju(:,4)*A0gyi( :,4) 176 & + giju(:,2)*A0gyi( :,9) 177 & + giju(:,6)*A0gyi(:,14) 178 179 gAgyi(:,10) = giju(:,4)*A0gyi( :,5) 180 & + giju(:,2)*A0gyi(:,10) 181 & + giju(:,6)*A0gyi(:,15) 182 183 gAgyi(:,11) = giju(:,5)*A0gyi( :,1) 184 & + giju(:,6)*A0gyi( :,6) 185 & + giju(:,3)*A0gyi(:,11) 186 187 gAgyi(:,12) = giju(:,5)*A0gyi( :,2) 188 & + giju(:,6)*A0gyi( :,7) 189 & + giju(:,3)*A0gyi(:,12) 190 191 gAgyi(:,13) = giju(:,5)*A0gyi( :,3) 192 & + giju(:,6)*A0gyi( :,8) 193 & + giju(:,3)*A0gyi(:,13) 194 195 gAgyi(:,14) = giju(:,5)*A0gyi( :,4) 196 & + giju(:,6)*A0gyi( :,9) 197 & + giju(:,3)*A0gyi(:,14) 198 199 gAgyi(:,15) = giju(:,5)*A0gyi( :,5) 200 & + giju(:,6)*A0gyi(:,10) 201 & + giju(:,3)*A0gyi(:,15) 202c 203c... the denominator term of the DC factor 204c... evaluation of the term Y,i.A0DC Y,j 205c 206 yiA0DCyj(:,1) = A0DC(:,1)*g1yi(:,1)**2 207 & + two*g1yi(:,1)*A0DC(:,2)*g1yi(:,5) 208 & + A0DC(:,3)*g1yi(:,2)**2 209 & + A0DC(:,3)*g1yi(:,3)**2 210 & + A0DC(:,3)*g1yi(:,4)**2 211 & + A0DC(:,4)*g1yi(:,5)**2 212 213 yiA0DCyj(:,2) = A0DC(:,1)*g2yi(:,1)**2 214 & + two*g2yi(:,1)*A0DC(:,2)*g2yi(:,5) 215 & + A0DC(:,3)*g2yi(:,2)**2 216 & + A0DC(:,3)*g2yi(:,3)**2 217 & + A0DC(:,3)*g2yi(:,4)**2 218 & + A0DC(:,4)*g2yi(:,5)**2 219 220 yiA0DCyj(:,3) = A0DC(:,1)*g3yi(:,1)**2 221 & + two*g3yi(:,1)*A0DC(:,2)*g3yi(:,5) 222 & + A0DC(:,3)*g3yi(:,2)**2 223 & + A0DC(:,3)*g3yi(:,3)**2 224 & + A0DC(:,3)*g3yi(:,4)**2 225 & + A0DC(:,4)*g3yi(:,5)**2 226 227 yiA0DCyj(:,4) = g1yi(:,1)*A0DC(:,1)*g2yi(:,1) 228 & + g1yi(:,1)*A0DC(:,2)*g2yi(:,5) 229 & + g1yi(:,2)*A0DC(:,3)*g2yi(:,2) 230 & + g1yi(:,3)*A0DC(:,3)*g2yi(:,3) 231 & + g1yi(:,4)*A0DC(:,3)*g2yi(:,4) 232 & + g1yi(:,5)*A0DC(:,2)*g2yi(:,1) 233 & + g1yi(:,5)*A0DC(:,4)*g2yi(:,5) 234 235 yiA0DCyj(:,5) = g1yi(:,1)*A0DC(:,1)*g3yi(:,1) 236 & + g1yi(:,1)*A0DC(:,2)*g3yi(:,5) 237 & + g1yi(:,2)*A0DC(:,3)*g3yi(:,2) 238 & + g1yi(:,3)*A0DC(:,3)*g3yi(:,3) 239 & + g1yi(:,4)*A0DC(:,3)*g3yi(:,4) 240 & + g1yi(:,5)*A0DC(:,2)*g3yi(:,1) 241 & + g1yi(:,5)*A0DC(:,4)*g3yi(:,5) 242 243 yiA0DCyj(:,6) = g2yi(:,1)*A0DC(:,1)*g3yi(:,1) 244 & + g2yi(:,1)*A0DC(:,2)*g3yi(:,5) 245 & + g2yi(:,2)*A0DC(:,3)*g3yi(:,2) 246 & + g2yi(:,3)*A0DC(:,3)*g3yi(:,3) 247 & + g2yi(:,4)*A0DC(:,3)*g3yi(:,4) 248 & + g2yi(:,5)*A0DC(:,2)*g3yi(:,1) 249 & + g2yi(:,5)*A0DC(:,4)*g3yi(:,5) 250c 251c.... -------------------------> DC factor <-------------------------- 252c 253c if ((ires .ne. 2) .or. (Jactyp .eq. 1)) then 254c 255c.... calculate 2-norm of Grad-local-V with respect to A0 256c 257c.... DC-mallet 258c 259 if (iDC .eq. 1) then 260c 261 fact = one 262 if (ipord .eq. 2) fact = 0.9 263 if (ipord .eq. 3) fact = 0.75 264 265c 266 gnorm = one / ( 267 & giju(:,1)*yiA0DCyj(:,1) 268 & + two*giju(:,4)*yiA0DCyj(:,4) 269 & + two*giju(:,5)*yiA0DCyj(:,5) 270 & + giju(:,2)*yiA0DCyj(:,2) 271 & + two*giju(:,6)*yiA0DCyj(:,6) 272 & + giju(:,3)*yiA0DCyj(:,3) 273 & + epsM ) 274c 275c DC(:,intp)=dim((fact*sqrt(raLS*gnorm)),(rtLS*gnorm)) 276 DC(:,intp)=max(zero,(fact*sqrt(raLS*gnorm))-(rtLS*gnorm)) 277c 278c.... flop count 279c 280! flops = flops + 46*npro 281c 282 endif 283c 284c.... DC-quadratic 285c 286 if (iDC .eq. 2) then 287c 288 gnorm = one / ( 289 & giju(:,1)*yiA0DCyj(:,1) 290 & + two*giju(:,4)*yiA0DCyj(:,4) 291 & + two*giju(:,5)*yiA0DCyj(:,5) 292 & + giju(:,2)*yiA0DCyj(:,2) 293 & + two*giju(:,6)*yiA0DCyj(:,6) 294 & + giju(:,3)*yiA0DCyj(:,3) 295 & + epsM ) 296 297c 298 DC(:,intp) = two * rtLS * gnorm 299c 300c.... flop count 301c 302! flops = flops + 36*npro 303c 304 endif 305c 306c.... DC-min 307c 308 if (iDC .eq. 3) then 309c 310 fact = one 311 if (ipord .eq. 2) fact = pt5 312c 313 gnorm = one / ( 314 & giju(:,1)*yiA0DCyj(:,1) 315 & + two*giju(:,4)*yiA0DCyj(:,4) 316 & + two*giju(:,5)*yiA0DCyj(:,5) 317 & + giju(:,2)*yiA0DCyj(:,2) 318 & + two*giju(:,6)*yiA0DCyj(:,6) 319 & + giju(:,3)*yiA0DCyj(:,3) 320 & + epsM ) 321 322c 323c DC(:,intp) = min( dim(fact * sqrt(raLS * gnorm), 324 DC(:,intp) = min( max(zero,fact * sqrt(raLS * gnorm)- 325 & rtLS * gnorm), two * rtLS * gnorm ) 326c 327c.... flop count 328c 329! flops = flops + 48*npro 330c 331 endif 332c 333c endif 334c 335c.... ----------------------------> RHS <---------------------------- 336c 337c.... add the contribution of DC to ri and/or rmi 338c 339c.... ires = 1 or 3 340c 341 if ((ires .eq. 1) .or. (ires .eq. 3)) then 342c 343 ri ( :,1) = ri ( :,1) + DC(:,intp) * gAgyi( :,1) 344 rmi( :,1) = rmi( :,1) + DC(:,intp) * gAgyi( :,1) 345 ri ( :,2) = ri ( :,2) + DC(:,intp) * gAgyi( :,2) 346 rmi( :,2) = rmi( :,2) + DC(:,intp) * gAgyi( :,2) 347 ri ( :,3) = ri ( :,3) + DC(:,intp) * gAgyi( :,3) 348 rmi( :,3) = rmi( :,3) + DC(:,intp) * gAgyi( :,3) 349 ri ( :,4) = ri ( :,4) + DC(:,intp) * gAgyi( :,4) 350 rmi( :,4) = rmi( :,4) + DC(:,intp) * gAgyi( :,4) 351 ri ( :,5) = ri ( :,5) + DC(:,intp) * gAgyi( :,5) 352 rmi( :,5) = rmi( :,5) + DC(:,intp) * gAgyi( :,5) 353c 354 ri ( :,6) = ri ( :,6) + DC(:,intp) * gAgyi( :,6) 355 rmi( :,6) = rmi( :,6) + DC(:,intp) * gAgyi( :,6) 356 ri ( :,7) = ri ( :,7) + DC(:,intp) * gAgyi( :,7) 357 rmi( :,7) = rmi( :,7) + DC(:,intp) * gAgyi( :,7) 358 ri ( :,8) = ri ( :,8) + DC(:,intp) * gAgyi( :,8) 359 rmi( :,8) = rmi( :,8) + DC(:,intp) * gAgyi( :,8) 360 ri ( :,9) = ri ( :,9) + DC(:,intp) * gAgyi( :,9) 361 rmi( :,9) = rmi( :,9) + DC(:,intp) * gAgyi( :,9) 362 ri (:,10) = ri (:,10) + DC(:,intp) * gAgyi(:,10) 363 rmi(:,10) = rmi(:,10) + DC(:,intp) * gAgyi(:,10) 364c 365 ri (:,11) = ri (:,11) + DC(:,intp) * gAgyi(:,11) 366 rmi(:,11) = rmi(:,12) + DC(:,intp) * gAgyi(:,12) 367 ri (:,12) = ri (:,12) + DC(:,intp) * gAgyi(:,12) 368 rmi(:,12) = rmi(:,12) + DC(:,intp) * gAgyi(:,12) 369 ri (:,13) = ri (:,13) + DC(:,intp) * gAgyi(:,13) 370 rmi(:,13) = rmi(:,13) + DC(:,intp) * gAgyi(:,13) 371 ri (:,14) = ri (:,14) + DC(:,intp) * gAgyi(:,14) 372 rmi(:,14) = rmi(:,14) + DC(:,intp) * gAgyi(:,14) 373 ri (:,15) = ri (:,15) + DC(:,intp) * gAgyi(:,15) 374 rmi(:,15) = rmi(:,15) + DC(:,intp) * gAgyi(:,15) 375c 376! flops = flops + 45*npro 377c 378 endif 379c 380c.... ires = 2 381c 382 if (ires .eq. 2) then 383c 384 rmi( :,1) = rmi( :,1) + DC(:,intp) * gAgyi( :,1) 385 rmi( :,2) = rmi( :,2) + DC(:,intp) * gAgyi( :,2) 386 rmi( :,3) = rmi( :,3) + DC(:,intp) * gAgyi( :,3) 387 rmi( :,4) = rmi( :,4) + DC(:,intp) * gAgyi( :,4) 388 rmi( :,5) = rmi( :,5) + DC(:,intp) * gAgyi( :,5) 389c 390 rmi( :,6) = rmi( :,6) + DC(:,intp) * gAgyi( :,6) 391 rmi( :,7) = rmi( :,7) + DC(:,intp) * gAgyi( :,7) 392 rmi( :,8) = rmi( :,8) + DC(:,intp) * gAgyi( :,8) 393 rmi( :,9) = rmi( :,9) + DC(:,intp) * gAgyi( :,9) 394 rmi(:,10) = rmi(:,10) + DC(:,intp) * gAgyi(:,10) 395c 396 rmi(:,11) = rmi(:,11) + DC(:,intp) * gAgyi(:,11) 397 rmi(:,12) = rmi(:,12) + DC(:,intp) * gAgyi(:,12) 398 rmi(:,13) = rmi(:,13) + DC(:,intp) * gAgyi(:,13) 399 rmi(:,14) = rmi(:,14) + DC(:,intp) * gAgyi(:,14) 400 rmi(:,15) = rmi(:,15) + DC(:,intp) * gAgyi(:,15) 401c 402! flops = flops + 30*npro 403c 404 endif 405c 406c.... -------------------------> Stiffness <-------------------------- 407c 408c.... add the contribution of DC to stiff 409c 410 if (iprec .eq. 1) then ! leave out of LHS, when called from itrres 411 nflow2=two*nflow 412 do j = 1, nflow 413 do i = 1, nflow 414 dtmp(:)=A0(:,i,j)*DC(:,intp) 415c 416c.... add (DC g^1 A0) to stiff [1,1] 417c 418 stiff(:,i,j) = stiff(:,i,j) 419 & + dtmp(:)*giju(:,1) 420c 421c.... add (DC g^1 A0) to stiff [1,2] 422c 423 424 stiff(:,i,j+nflow) = stiff(:,i,j+nflow) 425 & + dtmp(:)*giju(:,4) 426c 427c.... add (DC g^1 A0) to stiff [1,3] 428c 429 430 stiff(:,i,j+nflow2) = stiff(:,i,j+nflow2) 431 & + dtmp(:)*giju(:,5) 432 433c.... add (DC g^1 A0) to stiff [2,1] (similarly below) 434c 435 436 stiff(:,i+nflow,j) = stiff(:,i+nflow,j) 437 & + dtmp(:)*giju(:,4) 438 439 stiff(:,i+nflow,j+nflow) = stiff(:,i+nflow,j+nflow) 440 & + dtmp(:)*giju(:,2) 441 442 stiff(:,i+nflow,j+nflow2) = stiff(:,i+nflow,j+nflow2) 443 & + dtmp(:)*giju(:,6) 444 445 stiff(:,i+nflow2,j) = stiff(:,i+nflow2,j) 446 & + dtmp(:)*giju(:,5) 447 448 stiff(:,i+nflow2,j+nflow) = stiff(:,i+nflow2,j+nflow) 449 & + dtmp(:)*giju(:,6) 450 451 stiff(:,i+nflow2,j+nflow2) = stiff(:,i+nflow2,j+nflow2) 452 & + dtmp(:)*giju(:,3) 453 enddo 454 enddo 455c 456c.... flop count 457c 458! flops = flops + 210*npro 459c 460c.... end of stiffness 461c 462 endif 463c 464c.... return 465c 466 return 467 end 468c 469 subroutine e3dcSclr ( g1yti, g2yti, g3yti, 470 & A0t, raLSt, rTLSt, 471 & DCt, giju, 472 & rti, rmti, stifft) 473c 474c 475c---------------------------------------------------------------------- 476c 477c This routine calculates the contribution of the Discontinuity- 478c Capturing operator to RHS and preconditioner for the scalar solve. 479c 480c g1yti (nflow,npro) : grad-y in direction 1 481c g2yti (nflow,npro) : grad-y in direction 2 482c g3yti (nflow,npro) : grad-y in direction 3 483c A0 (nsymdf,npro) : A0 matrix (Symm. storage) 484c raLS (npro) : square of LS residual (A0inv norm) 485c rtLS (npro) : square of LS residual (Tau norm) 486c giju (6,npro) : metric matrix 487c DC (ngauss,npro) : discontinuity-capturing factor 488c intp : integration point number 489c 490c output: 491c ri (nflow*(nsd+1),npro) : partial residual 492c rmi (nflow*(nsd+1),npro) : partial modified residual 493c stiff (nsymdf,6,npro) : diffusivity matrix 494c DC (npro) : discontinuity-capturing factor 495c 496c 497c Zdenek Johan, Summer 1990. (Modified from e2dc.f) 498c Zdenek Johan, Winter 1991. (Recoded) 499c Zdenek Johan, Winter 1991. (Fortran 90) 500c---------------------------------------------------------------------- 501c 502 include "common.h" 503c 504 dimension g1yti(npro), g2yti(npro), 505 & g3yti(npro), A0t(npro), 506 & raLSt(npro), rtLSt(npro), 507 & giju(npro,6), DCt(npro,ngauss), 508 & rti(npro,nsd+1), rmti(npro,nsd+1), 509 & stifft(npro,nsd,nsd), dtmp(npro) 510c 511 512 dimension ggyit(npro,nflow), gAgyit(npro,3), 513 & gnormt(npro), A0gyit(npro,3), 514 & yiA0DCyjt(npro,6), A0DCt(npro) 515c 516c ... -----------------------> initialize <---------------------------- 517c 518 A0gyit = zero 519 gAgyit = zero 520 yiA0DCyjt = zero 521 DCt = zero 522 A0DCt = A0t 523c.... -----------------------> global gradient <---------------------- 524c 525c.... calculate (A0 y_,j) --> A0gyit 526c 527c A0 Y_{,1} 528c 529 A0gyit( :,1) = A0t(:)*g1yti(:) 530c A0 Y_{,2} 531 A0gyit( :,2) = A0t(:)*g2yti(:) 532c A0 Y_{,3} 533 A0gyit( :,3) = A0t(:)*g3yti(:) 534c 535c.... calculate (giju A0 y_,j) --> gAgyit 536c 537 538 gAgyit( :,1) = giju(:,1)*A0gyit( :,1) 539 & + giju(:,4)*A0gyit( :,2) 540 & + giju(:,5)*A0gyit( :,3) 541 542 gAgyit( :,2) = giju(:,4)*A0gyit( :,1) 543 & + giju(:,2)*A0gyit( :,2) 544 & + giju(:,6)*A0gyit( :,3) 545 546 gAgyit( :,3) = giju(:,5)*A0gyit( :,1) 547 & + giju(:,6)*A0gyit( :,2) 548 & + giju(:,3)*A0gyit( :,3) 549c 550c... the denominator term of the DC factor 551c... evaluation of the term Y,i.A0DC Y,j 552c 553 yiA0DCyjt(:,1) = A0DCt(:)*g1yti(:)**2 554c 555 yiA0DCyjt(:,2) = A0DCt(:)*g2yti(:)**2 556c 557 yiA0DCyjt(:,3) = A0DCt(:)*g3yti(:)**2 558c 559 yiA0DCyjt(:,4) = A0DCt(:)*g1yti(:)*g2yti(:) 560c 561 yiA0DCyjt(:,5) = A0DCt(:)*g1yti(:)*g3yti(:) 562c 563 yiA0DCyjt(:,6) = A0DCt(:)*g2yti(:)*g3yti(:) 564c 565c 566c.... -------------------------> DC factor <-------------------------- 567c 568c if ((ires .ne. 2) .or. (Jactyp .eq. 1)) then 569c 570c.... calculate 2-norm of Grad-local-V with respect to A0 571c 572c.... DC-mallet 573c 574 if (iDCsclr(1) .eq. 1) then 575c 576 fact = one 577 if (ipord .eq. 2) fact = 0.9 578 if (ipord .eq. 3) fact = 0.75 579 580c 581 gnormt = one / ( 582 & giju(:,1)*yiA0DCyjt(:,1) 583 & + two*giju(:,4)*yiA0DCyjt(:,4) 584 & + two*giju(:,5)*yiA0DCyjt(:,5) 585 & + giju(:,2)*yiA0DCyjt(:,2) 586 & + two*giju(:,6)*yiA0DCyjt(:,6) 587 & + giju(:,3)*yiA0DCyjt(:,3) 588 & + epsM ) 589c 590c DCt(:,intp)=dim((fact*sqrt(raLSt*gnormt)),(rtLSt*gnormt)) 591 DCt(:,intp)=max(zero,(fact*sqrt(raLSt*gnormt))-(rtLSt*gnormt)) 592c 593c.... flop count 594c 595! flops = flops + 46*npro 596c 597 endif 598c 599c.... DC-quadratic 600c 601 if (iDCSclr(1) .eq. 2) then 602c 603 gnormt = one / ( 604 & giju(:,1)*yiA0DCyjt(:,1) 605 & + two*giju(:,4)*yiA0DCyjt(:,4) 606 & + two*giju(:,5)*yiA0DCyjt(:,5) 607 & + giju(:,2)*yiA0DCyjt(:,2) 608 & + two*giju(:,6)*yiA0DCyjt(:,6) 609 & + giju(:,3)*yiA0DCyjt(:,3) 610 & + epsM ) 611 612c 613 DCt(:,intp) = two * rtLSt * gnormt 614c 615c.... flop count 616c 617! flops = flops + 36*npro 618c 619 endif 620c 621c.... DC-min 622c 623 if (iDCSclr(1) .eq. 3) then 624c 625 fact = one 626 if (ipord .eq. 2) fact = pt5 627c 628 gnormt = one / ( 629 & giju(:,1)*yiA0DCyjt(:,1) 630 & + two*giju(:,4)*yiA0DCyjt(:,4) 631 & + two*giju(:,5)*yiA0DCyjt(:,5) 632 & + giju(:,2)*yiA0DCyjt(:,2) 633 & + two*giju(:,6)*yiA0DCyjt(:,6) 634 & + giju(:,3)*yiA0DCyjt(:,3) 635 & + epsM ) 636 637c 638c DCt(:,intp) = min( dim(fact * sqrt(raLSt * gnormt), 639 DCt(:,intp) = min( max(zero,fact * sqrt(raLSt * gnormt)- 640 & rtLSt * gnormt), two * rtLSt * gnormt ) 641c 642c.... flop count 643c 644! flops = flops + 48*npro 645c 646 endif 647c 648c endif 649c DCt=DCt*two 650c 651c.... ----------------------------> RHS <---------------------------- 652c 653c.... add the contribution of DC to ri and/or rmi 654c 655c.... ires = 1 or 3 656c 657 if ((ires .eq. 1) .or. (ires .eq. 3)) then 658c 659 rti ( :,1) = rti ( :,1) + DCt(:,intp) * gAgyit( :,1) 660 rmti( :,1) = rmti( :,1) + DCt(:,intp) * gAgyit( :,1) 661 rti ( :,2) = rti ( :,2) + DCt(:,intp) * gAgyit( :,2) 662 rmti( :,2) = rmti( :,2) + DCt(:,intp) * gAgyit( :,2) 663 rti ( :,3) = rti ( :,3) + DCt(:,intp) * gAgyit( :,3) 664 rmti( :,3) = rmti( :,3) + DCt(:,intp) * gAgyit( :,3) 665 666c 667! flops = flops + 45*npro 668c 669 endif 670c 671c.... ires = 2 672c 673 if (ires .eq. 2) then 674c 675 rmti( :,1) = rmti( :,1) + DCt(:,intp) * gAgyit( :,1) 676 rmti( :,2) = rmti( :,2) + DCt(:,intp) * gAgyit( :,2) 677 rmti( :,3) = rmti( :,3) + DCt(:,intp) * gAgyit( :,3) 678 679c 680! flops = flops + 30*npro 681c 682 endif 683c 684c.... -------------------------> Stiffness <-------------------------- 685c 686c.... add the contribution of DC to stiff 687c$$$c 688c if (iprec .eq. 1) then !leave out of LHS, if called from itrres 689 !anyway matrix free not implemented for scalar 690 dtmp(:)=A0t(:)*DCt(:,intp) 691c 692c.... add (DC g^1 A0) to stifft [1,1] 693c 694 stifft(:,1,1) = stifft(:,1,1) 695 & + dtmp(:)*giju(:,1) 696c 697c.... add (DC g^1 A0) to stifft [1,2] 698c 699 stifft(:,1,2) = stifft(:,1,2) 700 & + dtmp(:)*giju(:,4) 701c 702c.... add (DC g^1 A0) to stifft [1,3] 703c 704 stifft(:,1,3) = stifft(:,1,3) 705 & + dtmp(:)*giju(:,5) 706 707c.... add (DC g^1 A0) to stifft [2,1] (similarly below) 708c 709 stifft(:,2,1) = stifft(:,2,1) 710 & + dtmp(:)*giju(:,4) 711 712 stifft(:,2,2) = stifft(:,2,2) 713 & + dtmp(:)*giju(:,2) 714 715 stifft(:,2,3) = stifft(:,2,3) 716 & + dtmp(:)*giju(:,6) 717 718 stifft(:,3,1) = stifft(:,3,1) 719 & + dtmp(:)*giju(:,5) 720 721 stifft(:,3,2) = stifft(:,3,2) 722 & + dtmp(:)*giju(:,6) 723 724 stifft(:,3,3) = stifft(:,3,3) 725 & + dtmp(:)*giju(:,3) 726c 727c.... flop count 728c 729! flops = flops + 210*npro 730c 731c.... end of stiffness 732c 733c endif 734c 735c.... return 736c 737 return 738 end 739c 740