1 : C CLI272 SOURCE MAGN 09/07/30 21:16:14 6458 2 : SUBROUTINE CLI272(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: 'INJE' -- 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 52 : & ,PSRF,RHOUF,DECIDP,DUNDPC,DUXDPC,DUYDPC 53 : & ,UN,RHO,P,COEF 54 : & ,BR1,BOT,TOP 55 : 56 : CHARACTER*(8) TYPE 57 : -INC CCOPTIO 58 : -INC SMLMOTS 59 : -INC SMELEME 60 : POINTEUR MELEFC.MELEME 61 : -INC SMLENTI 62 : POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI 63 : -INC SMCHPOI 64 : POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL, 65 : & MPVC.MPOVAL, MPPC.MPOVAL, MPLIM.MPOVAL, MPYC.MPOVAL 66 : -INC SMMATRIK 67 : POINTEUR CELL.IZAFM 68 : C------------------------------------------------------- 69 : -INC SMLREEL 70 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL 71 : C------------------------------------------------------- 72 : C********* Les Jacobians ****************************** 73 : C------------------------------------------------------- 74 : SEGMENT JACEL 75 : REAL*8 JAC(3+NSP,3+NSP) 76 : ENDSEGMENT 77 : POINTEUR JTL.JACEL, WL.JACEL, JTT.JACEL 78 : C------------------------------------------------------------- 79 : C******* Les fractionines massiques ************************** 80 : C------------------------------------------------------------- 81 : SEGMENT FRAMAS 82 : REAL*8 YET(NSP) 83 : ENDSEGMENT 84 : POINTEUR YC.FRAMAS, YF.FRAMAS 85 : C------------------------------------------------------- 86 : C********** Les CP's and CV's *********************** 87 : C------------------------------------------------------- 88 : SEGMENT GCONST 89 : REAL*8 GC(NSP) 90 : ENDSEGMENT 91 : POINTEUR CP.GCONST, CV.GCONST 92 : C------------------------------------------------------------- 93 : C********** Segments for the vectors *********************** 94 : C------------------------------------------------------------- 95 : SEGMENT VECEL 96 : REAL*8 VV(NSP) 97 : ENDSEGMENT 98 : POINTEUR DGDYC.VECEL 99 : C-------------------------------------------------- 100 : SEGINI JTL 101 : SEGINI JTT 102 : SEGINI WL 103 : C---------------------------------------------------- 104 : C**** KRIPAD pour la correspondance global/local 105 : C---------------------------------------------------- 105 : CALL KRIPAD(MELEMC,MLEMC) 105 : CALL KRIPAD(MELECB,MLEMCB) 105 : CALL KRIPAD(MELEMF,MLEMF) 106 : C---------------------------------------------------- 107 : C**** CHPOINTs de la table DOMAINE 108 : C---------------------------------------------------- 108 : CALL LICHT(INORM,MPNORM,TYPE,ICEL) 108 : CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL) 108 : CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL) 109 : C---------------------------------------------------- 110 : C**** CHPOINTs des variables 111 : C---------------------------------------------------- 111 : CALL LICHT(IROC,MPRC,TYPE,ICEL) 111 : CALL LICHT(IVITC,MPVC,TYPE,ICEL) 111 : CALL LICHT(IPC,MPPC,TYPE,ICEL) 111 : CALL LICHT(IYC,MPYC,TYPE,ICEL) 111 : CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL) 112 : C-------------------------------------------------------- 113 : C**** Boucle sur le face pour le calcul des invariants de 114 : C Riemann et du flux 115 : C-------------------------------------------------------- 116 : SEGACT MELEFC 117 : NFAC=MELEFC.NUM(/2) 118 : C--------------------------------- 119 : C**** Objet MATRIK 120 : C--------------------------------- 121 : NRIGE = 7 122 : NMATRI = 1 123 : NKID = 9 124 : NKMT = 7 125 : C--------------------------------- 126 : SEGINI MATRIK 127 : IJACO = MATRIK 128 : MATRIK.IRIGEL(1,1) = MELRES 129 : MATRIK.IRIGEL(2,1) = MELRES 130 : C--------------------------------- 131 : C**** Matrice non symetrique 132 : C--------------------------------- 133 : MATRIK.IRIGEL(7,1) = 2 134 : C--------------------------------- 135 : NBME = (3+NSP)*(3+NSP) 136 : NBSOUS = 1 137 : SEGINI IMATRI 138 : IF(IJAC.EQ.1)THEN 139 : MLMOTS=ILIINC 140 : ELSEIF(IJAC.EQ.2)THEN 141 : MLMOTS=ILIINP 142 : ENDIF 143 : SEGACT MLMOTS 144 : MATRIK.IRIGEL(4,1) = IMATRI 145 : C------------------------------------------- 146 : DO 1 J=1,(NSP+3) 147 : KV=(J-1)*(3+NSP) 148 : IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1) 149 : IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2) 150 : IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3) 151 : IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4) 152 : DO 2 I=1,(NSP-1) 153 : IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I) 154 : 2 CONTINUE 155 : 1 CONTINUE 156 : C----------------------------------------------- 157 : SEGDES MLMOTS 158 : MLMOTS=ILIINC 159 : SEGACT MLMOTS 160 : C----------------------------------------------- 161 : DO 3 J=1,(NSP+3) 162 : KV=(J-1)*(3+NSP) 163 : IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j) 164 : IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j) 165 : IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j) 166 : IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j) 167 : DO 4 I=1,(NSP-1) 168 : IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j) 169 : 4 CONTINUE 170 : 3 CONTINUE 171 : C----------------------------------------------- 172 : C----------------------------------------------- 173 : SEGDES MLMOTS 174 : NBEL = NFAC 175 : NBSOUS = 1 176 : NP = 1 177 : MP = 1 178 : C----------------------------------------------------------- 179 : C----------------------------------------------------------- 180 : DO 5 I=1,NBME 181 : SEGINI CELL 182 : IMATRI.LIZAFM(1,I) = CELL 183 : 5 CONTINUE 184 : C--------------------------------- 185 : C--------------------------------- 186 : C**** Fin definition MATRIK 187 : C--------------------------------- 188 : DO IFAC=1,NFAC,1 189 : NGF=MELEFC.NUM(1,IFAC) 190 : NGC=MELEFC.NUM(2,IFAC) 191 : NLF=MLEMF.LECT(NGF) 192 : NLC=MLEMC.LECT(NGC) 193 : NLCB=MLEMCB.LECT(NGF) 194 : VOLU=MPVOL.VPOCHA(NLC,1) 195 : SURF=MPSURF.VPOCHA(NLF,1) 196 : C In CASTEM les normales sont sortantes 197 : CNX=-1*MPNORM.VPOCHA(NLF,1) 198 : CNY=-1*MPNORM.VPOCHA(NLF,2) 199 : C---------------------------------------------- 200 : SEGINI CP, CV 201 : MLRECP = LRECP 202 : MLRECV = LRECV 203 : SEGACT MLRECP, MLRECV 204 : DO 10 I=1,(NSP-1) 205 : CP.GC(I)=MLRECP.PROG(I) 206 : CV.GC(I)=MLRECV.PROG(I) 207 : 10 CONTINUE 208 : CP.GC(NSP)=MLRECP.PROG(NSP) 209 : CV.GC(NSP)=MLRECV.PROG(NSP) 210 : C--------------------------------- 211 : C Variables au centre 212 : C--------------------------------- 213 : RC=MPRC.VPOCHA(NLC,1) 214 : PC=MPPC.VPOCHA(NLC,1) 215 : UXC=MPVC.VPOCHA(NLC,1) 216 : UYC=MPVC.VPOCHA(NLC,2) 217 : SEGINI YC 218 : SEGACT MPYC 219 : DO 100 I=1,(NSP-1) 220 : YC.YET(I)=MPYC.VPOCHA(NLC,I) 221 : 100 CONTINUE 222 : C--------------------------------- 223 : C Variables à la face 224 : C--------------------------------- 225 : RHOUF=MPLIM.VPOCHA(NLCB,1) 226 : PSRF=MPLIM.VPOCHA(NLCB,2) 227 : SEGINI YF 228 : SEGACT MPLIM 229 : DO 101 I=1,(NSP-1) 230 : YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I) 231 : 101 CONTINUE 232 : c------------------------------------------------------------- 233 : c Computing GAMMA at the cell-center 234 : c------------------------------------------------------------- 235 : top=0.0D0 236 : bot=0.0D0 237 : do 102 i=1,(nsp-1) 238 : top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp)) 239 : bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp)) 240 : 102 continue 241 : top=cp.gc(nsp)+top 242 : bot=cv.gc(nsp)+bot 243 : GAMC=top/bot 244 : C------------------------------------------------------------- 245 : SEGINI DGDYC 246 : do 41 i=1,(nsp-1) 247 : dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)- 248 : & GAMC*(cv.gc(i)-cv.gc(nsp)))/bot 249 : 41 continue 250 : c------------------------------------------------------------- 251 : C------------------------------------------------ 252 : C******* Densite, vitesse, pression sur le bord 253 : C------------------------------------------------ 254 : P=PC 255 : RHO=P/PSRF 256 : UN=RHOUF/RHO 257 : C------------------------------ 258 : C******* Derivatives w.r.t. PC 259 : C------------------------------ 260 : DUNDPC=-1*RHOUF/(RHO*RHO*PSRF) 261 : DUXDPC=CNX*DUNDPC 262 : DUYDPC=CNY*DUNDPC 263 : DECIDP=UN*DUNDPC 264 : C-------------------------------------------------------------- 265 : COEF=SURF/VOLU 266 : C------------------------------------------------------------- 267 : JTL.JAC(1,1) = 0.0D0 268 : JTL.JAC(1,2) = 0.0D0 269 : JTL.JAC(1,3) = 0.0D0 270 : JTL.JAC(1,4) = 0.0D0 271 : DO 107 I=1,(NSP-1) 272 : JTL.JAC(1,4+I) = 0.0D0 273 : 107 CONTINUE 274 : C---------------------------------------- 275 : JTL.JAC(2,1) = 0.0D0 276 : JTL.JAC(2,2) = 0.0D0 277 : JTL.JAC(2,3) = 0.0D0 278 : JTL.JAC(2,4) = (RHOUF*DUXDPC+CNX)*COEF 279 : DO 108 I=1,(NSP-1) 280 : JTL.JAC(2,4+I) = 0.0D0 281 : 108 CONTINUE 282 : C---------------------------------------- 283 : JTL.JAC(3,1) = 0.0D0 284 : JTL.JAC(3,2) = 0.0D0 285 : JTL.JAC(3,3) = 0.0D0 286 : JTL.JAC(3,4) = (RHOUF*DUYDPC+CNY)*COEF 287 : DO 109 I=1,(NSP-1) 288 : JTL.JAC(3,4+I) = 0.0D0 289 : 109 CONTINUE 290 : C---------------------------------------- 291 : JTL.JAC(4,1) = 0.0D0 292 : JTL.JAC(4,2) = 0.0D0 293 : JTL.JAC(4,3) = 0.0D0 294 : JTL.JAC(4,4) = (RHOUF*DECIDP)*COEF 295 : DO 110 I=1,(NSP-1) 296 : JTL.JAC(4,4+I) = 0.0D0 297 : 110 CONTINUE 298 : C---------------------------------------- 299 : DO 111 I=1,(NSP-1) 300 : JTL.JAC(4+I,1) = 0.0D0 301 : JTL.JAC(4+I,2) = 0.0D0 302 : JTL.JAC(4+I,3) = 0.0D0 303 : JTL.JAC(4+I,4) = 0.0D0 304 : DO 112 J=1,(NSP-1) 305 : JTL.JAC(4+I,4+J) = 0.0D0 306 : 112 CONTINUE 307 : 111 CONTINUE 308 : C--------------------------------------------------------- 309 : c matrix wl(i,j) represents the derivative of the i-component 310 : c of the vector of primitive variables of the left state with 311 : c respect to the j-component of the vector of the conservative 312 : c variables of the left state. 313 : c 314 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) - 315 : c vector of primitive variables; 316 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) - 317 : c vector of conservative variables. 318 : c------------------------------------------------------------- 319 : wl.jac(1,1)=1.0d0 320 : wl.jac(1,2)=0.0d0 321 : wl.jac(1,3)=0.0d0 322 : wl.jac(1,4)=0.0d0 323 : do 83 i=1,(nsp-1) 324 : wl.jac(1,4+i)=0.0d0 325 : 83 continue 326 : c------------------------------ 327 : wl.jac(2,1)=-UXC/RC 328 : wl.jac(2,2)=1.0d0/RC 329 : wl.jac(2,3)=0.0d0 330 : wl.jac(2,4)=0.0d0 331 : do 84 i=1,(nsp-1) 332 : wl.jac(2,4+i)=0.0d0 333 : 84 continue 334 : c------------------------------ 335 : wl.jac(3,1)=-UYC/RC 336 : wl.jac(3,2)=0.0d0 337 : wl.jac(3,3)=1.0d0/RC 338 : wl.jac(3,4)=0.0d0 339 : do 85 i=1,(nsp-1) 340 : wl.jac(3,4+i)=0.0d0 341 : 85 continue 342 : c------------------------------ 343 : br1=0.0d0 344 : do 86 i=1,(nsp-1) 345 : br1=br1+dgdyc.vv(i)*yc.yet(i) 346 : 86 continue 347 : br1=br1*PC/(RC*(GAMC-1.0D0)) 348 : wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1 349 : wl.jac(4,2)=-UXC*(GAMC-1.0D0) 350 : wl.jac(4,3)=-UYC*(GAMC-1.0D0) 351 : wl.jac(4,4)=(GAMC-1.0D0) 352 : do 87 i=1,(nsp-1) 353 : wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0)) 354 : 87 continue 355 : c------------------------------ 356 : do 88 i=1,(nsp-1) 357 : do 89 j=1,4 358 : wl.jac(4+i,j)=0.0d0 359 : if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC 360 : 89 continue 361 : c------------ 362 : do 890 j=5,(4+nsp-1) 363 : wl.jac(4+i,j)=0.0d0 364 : if(4+i.eq.j) then 365 : wl.jac(4+i,j)=1.0d0/RC 366 : endif 367 : 890 continue 368 : 88 continue 369 : c------------------------------------------------ 370 : C------------------------------------------------ 371 : do 114 i=1,(3+nsp) 372 : do 115 j=1,(3+nsp) 373 : jtt.jac(i,j)=0.0d0 374 : do 116 k=1,(3+nsp) 375 : jtt.jac(i,j)=jtt.jac(i,j)+jtl.jac(i,k)*wl.jac(k,j) 376 : 116 continue 377 : 115 continue 378 : 114 continue 379 : C---------------------------------------------------------------- 380 : C******* Jacobian with respect to conservative variables 381 : C---------------------------------------------------------------- 382 : IF(IJAC.EQ.1)THEN 383 : DO 9 II = 1,(3+NSP) 384 : DO 15 JJ = 1,(3+NSP) 385 : KV = (II-1)*(3+NSP) 386 : C---------------------------------- 387 : CELL = IMATRI.LIZAFM(1,KV+JJ) 388 : CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ) 389 : 15 CONTINUE 390 : 9 CONTINUE 391 : ELSEIF(IJAC.EQ.2)THEN 392 : DO 20 II = 1,(3+NSP) 393 : DO 25 JJ = 1,(3+NSP) 394 : KV = (II-1)*(3+NSP) 395 : C---------------------------------- 396 : CELL = IMATRI.LIZAFM(1,KV+JJ) 397 : CELL.AM(IFAC,1,1) = JTL.JAC(II,JJ) 398 : 25 CONTINUE 399 : 20 CONTINUE 400 : ENDIF 401 : c-------------------------------------------------- 402 : ENDDO 403 : C 404 : SEGDES MELEFC 405 : C 406 : SEGSUP MLEMC 407 : SEGSUP MLEMCB 408 : SEGSUP MLEMF 409 : C 410 : SEGDES MPNORM 411 : SEGDES MPVOL 412 : SEGDES MPSURF 413 : SEGDES MPRC 414 : SEGDES MPPC 415 : SEGDES MPVC 416 : SEGDES MPYC 417 : SEGDES MPLIM 418 : SEGDES YC 419 : SEGDES YF 420 : SEGDES CP 421 : SEGDES CV 422 : SEGDES JTL 423 : SEGDES JTT 424 : SEGDES WL 425 : SEGDES DGDYC 426 : SEGDES MATRIK 427 : DO 80 II=1,NBME 428 : CELL = IMATRI.LIZAFM(1,II) 429 : SEGDES CELL 430 : 80 CONTINUE 431 : SEGDES IMATRI 432 : C--------------------------------------------- 433 : 9999 CONTINUE 434 : RETURN 435 : END 436 : 437 : 438 : 439 : 440 : 441 : 442 : 443 : 444 : 445 :
© Cast3M 2003 - All rights reserved.
Disclaimer