1 : C CLI262 SOURCE MAGN 09/07/30 21:16:12 6458 2 : SUBROUTINE CLI262(NSP,MELEMF,MELEMC,MELECB,MELEFC,MELRES,INORM, 3 : & ICHPVO,ICHPSU,LRECP,LRECV, 4 : & IROC,IVITC,IPC,IYC,ICHLIM,ILIINC,ILIINP,IJAC,IJACO) 5 : C************************************************************************ 6 : C 7 : C PROJET : CASTEM 2000 8 : C 9 : C NOM : CLI112 10 : C 11 : C DESCRIPTION : Subroutine appellée par CLIM22 12 : C OPTION: 'INSU' -- Jacobian 13 : C 14 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 15 : C 16 : C AUTEUR : S. Kudriakov, DEN/DM2S/SFME/LTMF 17 : C 18 : C************************************************************************ 19 : C 20 : C APPELES (Calcul) : 21 : C 22 : C************************************************************************ 23 : C 24 : C HISTORIQUE (Anomalies et modifications éventuelles) 25 : C 26 : C HISTORIQUE : 27 : C 28 : C************************************************************************ 29 : C 30 : C---------------------------------------------------- 31 : C**** Variables de COOPTIO 32 : C---------------------------------------------------- 33 : c INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 34 : c & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 35 : c & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 36 : c & ,IECHO, IIMPI, IOSPI 37 : c & ,IDIM, IFICLE, IPREFI 38 : c & ,MCOORD 39 : c & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 40 : c & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 41 : c & ,NORINC,NORVAL,NORIND,NORVAD 42 : c & ,NUCROU, IPSAUV 43 : C 44 : IMPLICIT INTEGER(I-N) 45 : INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC 46 : & ,IGAMC,ICHLIM,ICEL,NFAC,IFAC,MELRES,IJACO 47 : & ,NGF,NGC,NLF,NLC,NLCB 48 : & ,ILIINC,ILIINP,IJAC,II,JJ,K 49 : & ,MP, NBEL, NBME, NBSOUS, NKID, NKMT, NMATRI, NP, NRIGE 50 : & ,NSP,I, IYC,J, LRECP,LRECV,KV 51 : REAL*8 VOLU,SURF,RC,PC,UXC,UYC,GAMC,CNX,CNY,CTX,CTY 52 : & ,PSRF,COEF 53 : & ,BR1,BOT,TOP,DUNDUY,GAMF,UNC 54 : & ,COEF5,DPSRUN,DUXDUN,DRDUN,DUYDUN,DUNDUX,DPDUN 55 : & ,GM1,USGM1,RF,UNF,UTF,UXF,UYF,HTF,PF,SF,ECIN 56 : & ,TVECT(2),NVECT(2),WVEC_L(4),WVEC_R(4) 57 : CHARACTER*(8) TYPE 58 : -INC CCOPTIO 59 : -INC SMLMOTS 60 : -INC SMELEME 61 : POINTEUR MELEFC.MELEME 62 : -INC SMLENTI 63 : POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI 64 : -INC SMCHPOI 65 : POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL, 66 : & MPVC.MPOVAL, MPPC.MPOVAL, MPLIM.MPOVAL, MPYC.MPOVAL 67 : -INC SMMATRIK 68 : POINTEUR CELL.IZAFM 69 : C------------------------------------------------------- 70 : -INC SMLREEL 71 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL 72 : C------------------------------------------------------- 73 : C********* Les Jacobians ****************************** 74 : C------------------------------------------------------- 75 : SEGMENT JACEL 76 : REAL*8 JAC(3+NSP,3+NSP) 77 : ENDSEGMENT 78 : POINTEUR JLL.JACEL, WL.JACEL, JRR.JACEL, JTT.JACEL, 79 : & JFL.JACEL, JFR.JACEL, JR.JACEL 80 : C------------------------------------------------------------- 81 : C******* Les fractionines massiques ************************** 82 : C------------------------------------------------------------- 83 : SEGMENT FRAMAS 84 : REAL*8 YET(NSP) 85 : ENDSEGMENT 86 : POINTEUR YC.FRAMAS, YF.FRAMAS 87 : C------------------------------------------------------- 88 : C********** Les CP's and CV's *********************** 89 : C------------------------------------------------------- 90 : SEGMENT GCONST 91 : REAL*8 GC(NSP) 92 : ENDSEGMENT 93 : POINTEUR CP.GCONST, CV.GCONST 94 : C------------------------------------------------------------- 95 : C********** Segments for the vectors *********************** 96 : C------------------------------------------------------------- 97 : SEGMENT VECEL 98 : REAL*8 VV(NSP) 99 : ENDSEGMENT 100 : POINTEUR DGDYC.VECEL 101 : C-------------------------------------------------- 102 : SEGINI JTT 103 : SEGINI JLL 104 : SEGINI JRR 105 : SEGINI JFL 106 : SEGINI JFR 107 : SEGINI WL 108 : C---------------------------------------------------- 109 : C**** KRIPAD pour la correspondance global/local 110 : C---------------------------------------------------- 110 : CALL KRIPAD(MELEMC,MLEMC) 110 : CALL KRIPAD(MELECB,MLEMCB) 110 : CALL KRIPAD(MELEMF,MLEMF) 111 : C---------------------------------------------------- 112 : C**** CHPOINTs de la table DOMAINE 113 : C---------------------------------------------------- 113 : CALL LICHT(INORM,MPNORM,TYPE,ICEL) 113 : CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL) 113 : CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL) 114 : C---------------------------------------------------- 115 : C**** CHPOINTs des variables 116 : C---------------------------------------------------- 116 : CALL LICHT(IROC,MPRC,TYPE,ICEL) 116 : CALL LICHT(IVITC,MPVC,TYPE,ICEL) 116 : CALL LICHT(IPC,MPPC,TYPE,ICEL) 116 : CALL LICHT(IYC,MPYC,TYPE,ICEL) 116 : CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL) 117 : C-------------------------------------------------------- 118 : C**** Boucle sur le face pour le calcul des invariants de 119 : C Riemann et du flux 120 : C-------------------------------------------------------- 121 : SEGACT MELEFC 122 : NFAC=MELEFC.NUM(/2) 123 : C--------------------------------- 124 : C**** Objet MATRIK 125 : C--------------------------------- 126 : NRIGE = 7 127 : NMATRI = 1 128 : NKID = 9 129 : NKMT = 7 130 : C--------------------------------- 131 : SEGINI MATRIK 132 : IJACO = MATRIK 133 : MATRIK.IRIGEL(1,1) = MELRES 134 : MATRIK.IRIGEL(2,1) = MELRES 135 : C--------------------------------- 136 : C**** Matrice non symetrique 137 : C--------------------------------- 138 : MATRIK.IRIGEL(7,1) = 2 139 : C--------------------------------- 140 : NBME = (3+NSP)*(3+NSP) 141 : NBSOUS = 1 142 : SEGINI IMATRI 143 : IF(IJAC.EQ.1)THEN 144 : MLMOTS=ILIINC 145 : ELSEIF(IJAC.EQ.2)THEN 146 : MLMOTS=ILIINP 147 : ENDIF 148 : SEGACT MLMOTS 149 : MATRIK.IRIGEL(4,1) = IMATRI 150 : C------------------------------------------- 151 : DO 1 J=1,(NSP+3) 152 : KV=(J-1)*(3+NSP) 153 : IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1) 154 : IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2) 155 : IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3) 156 : IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4) 157 : DO 2 I=1,(NSP-1) 158 : IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I) 159 : 2 CONTINUE 160 : 1 CONTINUE 161 : C----------------------------------------------- 162 : SEGDES MLMOTS 163 : MLMOTS=ILIINC 164 : SEGACT MLMOTS 165 : C----------------------------------------------- 166 : DO 3 J=1,(NSP+3) 167 : KV=(J-1)*(3+NSP) 168 : IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j) 169 : IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j) 170 : IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j) 171 : IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j) 172 : DO 4 I=1,(NSP-1) 173 : IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j) 174 : 4 CONTINUE 175 : 3 CONTINUE 176 : C----------------------------------------------- 177 : C----------------------------------------------- 178 : SEGDES MLMOTS 179 : NBEL = NFAC 180 : NBSOUS = 1 181 : NP = 1 182 : MP = 1 183 : C----------------------------------------------------------- 184 : C----------------------------------------------------------- 185 : DO 5 I=1,NBME 186 : SEGINI CELL 187 : IMATRI.LIZAFM(1,I) = CELL 188 : 5 CONTINUE 189 : C--------------------------------- 190 : C--------------------------------- 191 : C**** Fin definition MATRIK 192 : C--------------------------------- 193 : DO IFAC=1,NFAC,1 194 : NGF=MELEFC.NUM(1,IFAC) 195 : NGC=MELEFC.NUM(2,IFAC) 196 : NLF=MLEMF.LECT(NGF) 197 : NLC=MLEMC.LECT(NGC) 198 : NLCB=MLEMCB.LECT(NGF) 199 : VOLU=MPVOL.VPOCHA(NLC,1) 200 : SURF=MPSURF.VPOCHA(NLF,1) 201 : C In CASTEM les normales sont sortantes 202 : CNX=-1*MPNORM.VPOCHA(NLF,1) 203 : CNY=-1*MPNORM.VPOCHA(NLF,2) 204 : CTX=-1.0D0*CNY 205 : CTY=CNX 206 : C---------------------------------------------- 207 : SEGINI CP, CV 208 : MLRECP = LRECP 209 : MLRECV = LRECV 210 : SEGACT MLRECP, MLRECV 211 : DO 10 I=1,(NSP-1) 212 : CP.GC(I)=MLRECP.PROG(I) 213 : CV.GC(I)=MLRECV.PROG(I) 214 : 10 CONTINUE 215 : CP.GC(NSP)=MLRECP.PROG(NSP) 216 : CV.GC(NSP)=MLRECV.PROG(NSP) 217 : C--------------------------------- 218 : C Variables au centre 219 : C--------------------------------- 220 : RC=MPRC.VPOCHA(NLC,1) 221 : PC=MPPC.VPOCHA(NLC,1) 222 : UXC=MPVC.VPOCHA(NLC,1) 223 : UYC=MPVC.VPOCHA(NLC,2) 224 : SEGINI YC 225 : SEGACT MPYC 226 : DO 100 I=1,(NSP-1) 227 : YC.YET(I)=MPYC.VPOCHA(NLC,I) 228 : 100 CONTINUE 229 : UNC=(UXC*CNX)+(UYC*CNY) 230 : C--------------------------------- 231 : C Variables à la face 232 : C--------------------------------- 233 : SEGACT MPLIM 234 : HTF=MPLIM.VPOCHA(NLCB,1) 235 : SF=MPLIM.VPOCHA(NLCB,2) 236 : SEGINI YF 237 : DO 101 I=1,(NSP-1) 238 : YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I) 239 : 101 CONTINUE 240 : UNF=UNC 241 : UTF=0.0D0 242 : UXF=UNF*CNX+UTF*CTX 243 : UYF=UNF*CNY+UTF*CTY 244 : c------------------------------------------------------------- 245 : c Computing GAMMA at the face 246 : c------------------------------------------------------------- 247 : top=0.0D0 248 : bot=0.0D0 249 : do 102 i=1,(nsp-1) 250 : top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp)) 251 : bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp)) 252 : 102 continue 253 : top=cp.gc(nsp)+top 254 : bot=cv.gc(nsp)+bot 255 : GAMF=top/bot 256 : GM1=GAMF-1.0D0 257 : USGM1=1.0D0/GM1 258 : c------------------------------------------------------------- 259 : c Computing GAMMA at the cell-centre 260 : c------------------------------------------------------------- 261 : top=0.0D0 262 : bot=0.0D0 263 : do 103 i=1,(nsp-1) 264 : top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp)) 265 : bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp)) 266 : 103 continue 267 : top=cp.gc(nsp)+top 268 : bot=cv.gc(nsp)+bot 269 : GAMC=top/bot 270 : C------------------------------------------------------------- 271 : SEGINI DGDYC 272 : do 41 i=1,(nsp-1) 273 : dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)- 274 : & GAMC*(cv.gc(i)-cv.gc(nsp)))/bot 275 : 41 continue 276 : c------------------------------------------------------------- 277 : C------------------------------------------------ 278 : C******* Densite, vitesse, pression sur le bord 279 : C------------------------------------------------ 280 : ECIN=0.5D0*((UXF*UXF)+(UYF*UYF)) 281 : PSRF=(GM1/GAMF)*(HTF-ECIN) 282 : RF=PSRF/SF 283 : RF=RF**(1.0D0/GM1) 284 : PF=SF*(RF**GAMF) 285 : C------------------------------ 286 : C******* Derivatives w.r.t. PC 287 : C------------------------------ 288 : DPSRUN=-1*(GM1/GAMF)*UNC 289 : DRDUN=USGM1*RF/PSRF*DPSRUN 290 : DPDUN=GAMF*PSRF*DRDUN 291 : DUXDUN=CNX 292 : DUYDUN=CNY 293 : C------------------------------- 294 : DUNDUX=CNX 295 : DUNDUY=CNY 296 : C-------------------------------------------------------------- 297 : C------------------------------ 298 : C******* Derivatives 299 : C------------------------------ 300 : wvec_l(1)=RF 301 : wvec_l(2)=UXF 302 : wvec_l(3)=UYF 303 : wvec_l(4)=PF 304 : C-------------------------- 305 : wvec_r(1)=RC 306 : wvec_r(2)=UXC 307 : wvec_r(3)=UYC 308 : wvec_r(4)=PC 309 : C-------------------------- 310 : nvect(1)=CNX 311 : nvect(2)=CNY 312 : tvect(1)=CTX 313 : tvect(2)=CTY 313 : call copms2(nsp,jfl,jfr,wvec_l,wvec_r,nvect,tvect, 314 : & yf,mpyc,lrecp,lrecv,nlc) 315 : C-------------------------------------------------------------- 316 : COEF=SURF/VOLU 317 : JLL=JFL 318 : JRR=JFR 319 : SEGACT JLL 320 : SEGACT JRR 321 : SEGINI JR 322 : C------------------------------------------------------------- 323 : C Taking in to account the dependance of the variables 324 : C at the face on the variables at the centre (UNC) 325 : C------------------------------------------------------------- 326 : DO 105 I=1,(3+NSP) 327 : DO 1050 J=1,(3+NSP) 328 : IF(J .EQ. 2) THEN 329 : COEF5=(JLL.JAC(I,1)*DRDUN)+(JLL.JAC(I,2)*DUXDUN)+ 330 : & (JLL.JAC(I,3)*DUYDUN)+(JLL.JAC(I,4)*DPDUN) 331 : JR.JAC(I,J)=(COEF5*DUNDUX+JRR.JAC(I,J))*COEF 332 : ELSEIF(J .EQ. 3) THEN 333 : COEF5=(JLL.JAC(I,1)*DRDUN)+(JLL.JAC(I,2)*DUXDUN)+ 334 : & (JLL.JAC(I,3)*DUYDUN)+(JLL.JAC(I,4)*DPDUN) 335 : JR.JAC(I,J)=(COEF5*DUNDUY+JRR.JAC(I,J))*COEF 336 : ELSE 337 : JR.JAC(I,J)=JRR.JAC(I,J)*COEF 338 : ENDIF 339 : 1050 CONTINUE 340 : 105 CONTINUE 341 : C------------------------------------------------------------- 342 : c matrix wl(i,j) represents the derivative of the i-component 343 : c of the vector of primitive variables of the left state with 344 : c respect to the j-component of the vector of the conservative 345 : c variables of the left state. 346 : c 347 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) - 348 : c vector of primitive variables; 349 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) - 350 : c vector of conservative variables. 351 : c------------------------------------------------------------- 352 : wl.jac(1,1)=1.0d0 353 : wl.jac(1,2)=0.0d0 354 : wl.jac(1,3)=0.0d0 355 : wl.jac(1,4)=0.0d0 356 : do 83 i=1,(nsp-1) 357 : wl.jac(1,4+i)=0.0d0 358 : 83 continue 359 : c------------------------------ 360 : wl.jac(2,1)=-UXC/RC 361 : wl.jac(2,2)=1.0d0/RC 362 : wl.jac(2,3)=0.0d0 363 : wl.jac(2,4)=0.0d0 364 : do 84 i=1,(nsp-1) 365 : wl.jac(2,4+i)=0.0d0 366 : 84 continue 367 : c------------------------------ 368 : wl.jac(3,1)=-UYC/RC 369 : wl.jac(3,2)=0.0d0 370 : wl.jac(3,3)=1.0d0/RC 371 : wl.jac(3,4)=0.0d0 372 : do 85 i=1,(nsp-1) 373 : wl.jac(3,4+i)=0.0d0 374 : 85 continue 375 : c------------------------------ 376 : br1=0.0d0 377 : do 86 i=1,(nsp-1) 378 : br1=br1+dgdyc.vv(i)*yc.yet(i) 379 : 86 continue 380 : br1=br1*PC/(RC*(GAMC-1.0D0)) 381 : wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1 382 : wl.jac(4,2)=-UXC*(GAMC-1.0D0) 383 : wl.jac(4,3)=-UYC*(GAMC-1.0D0) 384 : wl.jac(4,4)=(GAMC-1.0D0) 385 : do 87 i=1,(nsp-1) 386 : wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0)) 387 : 87 continue 388 : c------------------------------ 389 : do 88 i=1,(nsp-1) 390 : do 89 j=1,4 391 : wl.jac(4+i,j)=0.0d0 392 : if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC 393 : 89 continue 394 : c------------ 395 : do 890 j=5,(4+nsp-1) 396 : wl.jac(4+i,j)=0.0d0 397 : if(4+i.eq.j) then 398 : wl.jac(4+i,j)=1.0d0/RC 399 : endif 400 : 890 continue 401 : 88 continue 402 : c------------------------------------------------ 403 : C------------------------------------------------ 404 : do 114 i=1,(3+nsp) 405 : do 115 j=1,(3+nsp) 406 : jtt.jac(i,j)=0.0d0 407 : do 116 k=1,(3+nsp) 408 : jtt.jac(i,j)=jtt.jac(i,j)+ 409 : & jr.jac(i,k)*wl.jac(k,j)/COEF 410 : 116 continue 411 : 115 continue 412 : 114 continue 413 : C---------------------------------------------------------------- 414 : C******* Jacobian with respect to conservative variables 415 : C---------------------------------------------------------------- 416 : IF(IJAC.EQ.1)THEN 417 : DO 9 II = 1,(3+NSP) 418 : DO 15 JJ = 1,(3+NSP) 419 : KV = (II-1)*(3+NSP) 420 : C---------------------------------- 421 : CELL = IMATRI.LIZAFM(1,KV+JJ) 422 : CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ)*COEF 423 : 15 CONTINUE 424 : 9 CONTINUE 425 : ELSEIF(IJAC.EQ.2)THEN 426 : DO 20 II = 1,(3+NSP) 427 : DO 25 JJ = 1,(3+NSP) 428 : KV = (II-1)*(3+NSP) 429 : C---------------------------------- 430 : CELL = IMATRI.LIZAFM(1,KV+JJ) 431 : CELL.AM(IFAC,1,1) = JR.JAC(II,JJ) 432 : 25 CONTINUE 433 : 20 CONTINUE 434 : ENDIF 435 : c-------------------------------------------------- 436 : ENDDO 437 : C 438 : SEGDES MELEFC 439 : C 440 : SEGDES MLEMC 441 : SEGDES MLEMCB 442 : SEGDES MLEMF 443 : C 444 : SEGDES MPNORM 445 : SEGDES MPVOL 446 : SEGDES MPSURF 447 : SEGDES MPRC 448 : SEGDES MPPC 449 : SEGDES MPVC 450 : SEGDES MPYC 451 : SEGDES MPLIM 452 : SEGDES YC 453 : SEGDES YF 454 : SEGDES CP 455 : SEGDES CV 456 : SEGDES JLL 457 : SEGDES JRR 458 : SEGDES JTT 459 : SEGDES JFL 460 : SEGDES JFR 461 : SEGDES JR 462 : SEGDES WL 463 : SEGDES DGDYC 464 : SEGDES MATRIK 465 : DO 80 II=1,NBME 466 : CELL = IMATRI.LIZAFM(1,II) 467 : SEGDES CELL 468 : 80 CONTINUE 469 : SEGDES IMATRI 470 : C--------------------------------------------- 471 : 9999 CONTINUE 472 : RETURN 473 : END 474 : 475 : 476 : 477 : 478 : 479 : 480 : 481 : 482 : 483 : 484 : 485 :
© Cast3M 2003 - All rights reserved.
Disclaimer