1 : C CLI223 SOURCE MAGN 09/07/30 21:16:08 6458 2 : SUBROUTINE CLI223(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 : CLI223 10 : C 11 : C DESCRIPTION : Subroutine appellée par CLIM22 12 : C Jacobian de residu pres de limit 13 : C OPTION: 'INRI' 2D 14 : C 15 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 16 : C 17 : C AUTEUR : S.Kudriakov, DEN/DM2S/SFME/LTMF 18 : C 19 : C************************************************************************ 20 : C 21 : C APPELES (Calcul) : 22 : C 23 : C************************************************************************ 24 : C 25 : C HISTORIQUE (Anomalies et modifications éventuelles) 26 : C 27 : C HISTORIQUE : 28 : C 29 : C************************************************************************ 30 : C 31 : C---------------------------------------------------- 32 : C**** Variables de COOPTIO 33 : C---------------------------------------------------- 34 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 35 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 36 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 37 : C & ,IECHO, IIMPI, IOSPI 38 : C & ,IDIM, IFICLE, IPREFI 39 : C & ,MCOORD 40 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 41 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 42 : C & ,NORINC,NORVAL,NORIND,NORVAD 43 : C & ,NUCROU, IPSAUV 44 : C 45 : IMPLICIT INTEGER(I-N) 46 : INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC 47 : & ,IGAMC,ICHLIM,ICEL,NFAC,IFAC,MELRES,IJACO 48 : & ,NGF,NGC,NLF,NLC,NLCB 49 : & ,ILIINC,ILIINP,IJAC,II,JJ,K 50 : & ,MP, NBEL, NBME, NBSOUS, NKID, NKMT, NMATRI, NP, NRIGE 51 : & ,NSP,I, IYC,J, LRECP,LRECV,KV 52 : REAL*8 VOLU,SURF,RC,PC,UXC,UYC,GAMC,CNX,CNY,CTX,CTY 53 : & ,RF,PF,UXF,UYF 54 : & ,UNC,UNF,UTF,SF,ASONC,ASONF,ASON 55 : & ,USGM1,G1,G3,ASON2,S,UT,UN,RHO,P,UX,UY 56 : & ,DUNDG1,DASDG1,DRHDG1,DPDG1,DUXDG1,DUYDG1 57 : & ,DFRDG1,DFMXG1,DFMYG1,DFEDG1 58 : & ,DG1DR,DG1DP,DG1DUX,DG1DUY,COEF,DACDG 59 : & ,BR1,BOT,TOP,GAMF 60 : c & ,EPS,CACCA,CACC1,CEL 61 : 62 : CHARACTER*(8) TYPE 63 : -INC CCOPTIO 64 : -INC SMLMOTS 65 : -INC SMELEME 66 : POINTEUR MELEFC.MELEME 67 : -INC SMLENTI 68 : POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI 69 : -INC SMCHPOI 70 : POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL, 71 : & MPVC.MPOVAL, MPPC.MPOVAL, MPLIM.MPOVAL, MPYC.MPOVAL 72 : -INC SMMATRIK 73 : POINTEUR CELL.IZAFM 74 : C------------------------------------------------------- 75 : -INC SMLREEL 76 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL 77 : C------------------------------------------------------- 78 : C********* Les Jacobians ****************************** 79 : C------------------------------------------------------- 80 : SEGMENT JACEL 81 : REAL*8 JAC(3+NSP,3+NSP) 82 : ENDSEGMENT 83 : POINTEUR JTL.JACEL, WL.JACEL, JTT.JACEL 84 : C------------------------------------------------------------- 85 : C******* Les fractionines massiques ************************** 86 : C------------------------------------------------------------- 87 : SEGMENT FRAMAS 88 : REAL*8 YET(NSP) 89 : ENDSEGMENT 90 : POINTEUR YC.FRAMAS, YF.FRAMAS 91 : C------------------------------------------------------- 92 : C********** Les CP's and CV's *********************** 93 : C------------------------------------------------------- 94 : SEGMENT GCONST 95 : REAL*8 GC(NSP) 96 : ENDSEGMENT 97 : POINTEUR CP.GCONST, CV.GCONST 98 : C------------------------------------------------------------- 99 : C********** Segments for the vectors *********************** 100 : C------------------------------------------------------------- 101 : SEGMENT VECEL 102 : REAL*8 VV(NSP) 103 : ENDSEGMENT 104 : POINTEUR DYDG1.VECEL, DFRYG1.VECEL, 105 : & DG1DY.VECEL, DGDYC.VECEL 106 : C---------------------------------------------------- 107 : SEGINI JTL 108 : SEGINI JTT 109 : SEGINI WL 110 : C---------------------------------------------------- 111 : C**** KRIPAD pour la correspondance global/local 112 : C---------------------------------------------------- 112 : CALL KRIPAD(MELEMC,MLEMC) 112 : CALL KRIPAD(MELECB,MLEMCB) 112 : CALL KRIPAD(MELEMF,MLEMF) 113 : C---------------------------------------------------- 114 : C**** CHPOINTs de la table DOMAINE 115 : C---------------------------------------------------- 115 : CALL LICHT(INORM,MPNORM,TYPE,ICEL) 115 : CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL) 115 : CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL) 116 : C---------------------------------------------------- 117 : C**** CHPOINTs des variables 118 : C---------------------------------------------------- 118 : CALL LICHT(IROC,MPRC,TYPE,ICEL) 118 : CALL LICHT(IVITC,MPVC,TYPE,ICEL) 118 : CALL LICHT(IPC,MPPC,TYPE,ICEL) 118 : CALL LICHT(IYC,MPYC,TYPE,ICEL) 118 : CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL) 119 : C-------------------------------------------------------- 120 : C**** Boucle sur le face pour le calcul des invariants de 121 : C Riemann et du flux 122 : C-------------------------------------------------------- 123 : SEGACT MELEFC 124 : NFAC=MELEFC.NUM(/2) 125 : C--------------------------------- 126 : C**** Objet MATRIK 127 : C--------------------------------- 128 : NRIGE = 7 129 : NMATRI = 1 130 : NKID = 9 131 : NKMT = 7 132 : C--------------------------------- 133 : SEGINI MATRIK 134 : IJACO = MATRIK 135 : MATRIK.IRIGEL(1,1) = MELRES 136 : MATRIK.IRIGEL(2,1) = MELRES 137 : C--------------------------------- 138 : C**** Matrice non symetrique 139 : C--------------------------------- 140 : MATRIK.IRIGEL(7,1) = 2 141 : C--------------------------------- 142 : NBME = (3+NSP)*(3+NSP) 143 : NBSOUS = 1 144 : SEGINI IMATRI 145 : IF(IJAC.EQ.1)THEN 146 : MLMOTS=ILIINC 147 : ELSEIF(IJAC.EQ.2)THEN 148 : MLMOTS=ILIINP 149 : ENDIF 150 : SEGACT MLMOTS 151 : MATRIK.IRIGEL(4,1) = IMATRI 152 : C------------------------------------------- 153 : DO 1 J=1,(NSP+3) 154 : KV=(J-1)*(3+NSP) 155 : IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1) 156 : IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2) 157 : IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3) 158 : IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4) 159 : DO 2 I=1,(NSP-1) 160 : IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I) 161 : 2 CONTINUE 162 : 1 CONTINUE 163 : C----------------------------------------------- 164 : SEGDES MLMOTS 165 : MLMOTS=ILIINC 166 : SEGACT MLMOTS 167 : C----------------------------------------------- 168 : DO 3 J=1,(NSP+3) 169 : KV=(J-1)*(3+NSP) 170 : IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j) 171 : IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j) 172 : IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j) 173 : IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j) 174 : DO 4 I=1,(NSP-1) 175 : IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j) 176 : 4 CONTINUE 177 : 3 CONTINUE 178 : C----------------------------------------------- 179 : C----------------------------------------------- 180 : SEGDES MLMOTS 181 : NBEL = NFAC 182 : NBSOUS = 1 183 : NP = 1 184 : MP = 1 185 : C----------------------------------------------------------- 186 : C----------------------------------------------------------- 187 : DO 5 I=1,NBME 188 : SEGINI CELL 189 : IMATRI.LIZAFM(1,I) = CELL 190 : 5 CONTINUE 191 : C--------------------------------- 192 : C--------------------------------- 193 : C--------------------------------- 194 : C**** Fin definition MATRIK 195 : C--------------------------------- 196 : DO IFAC=1,NFAC,1 197 : NGF=MELEFC.NUM(1,IFAC) 198 : NGC=MELEFC.NUM(2,IFAC) 199 : NLF=MLEMF.LECT(NGF) 200 : NLC=MLEMC.LECT(NGC) 201 : NLCB=MLEMCB.LECT(NGF) 202 : VOLU=MPVOL.VPOCHA(NLC,1) 203 : SURF=MPSURF.VPOCHA(NLF,1) 204 : C In CASTEM les normales sont sortantes 205 : CNX=-1*MPNORM.VPOCHA(NLF,1) 206 : CNY=-1*MPNORM.VPOCHA(NLF,2) 207 : CTX=-1.0D0*CNY 208 : CTY=CNX 209 : C---------------------------------------------- 210 : SEGINI CP, CV 211 : MLRECP = LRECP 212 : MLRECV = LRECV 213 : SEGACT MLRECP, MLRECV 214 : DO 10 I=1,(NSP-1) 215 : CP.GC(I)=MLRECP.PROG(I) 216 : CV.GC(I)=MLRECV.PROG(I) 217 : 10 CONTINUE 218 : CP.GC(NSP)=MLRECP.PROG(NSP) 219 : CV.GC(NSP)=MLRECV.PROG(NSP) 220 : C--------------------------------- 221 : C Variables au centre 222 : C--------------------------------- 223 : RC=MPRC.VPOCHA(NLC,1) 224 : PC=MPPC.VPOCHA(NLC,1) 225 : UXC=MPVC.VPOCHA(NLC,1) 226 : UYC=MPVC.VPOCHA(NLC,2) 227 : SEGINI YC 228 : SEGACT MPYC 229 : DO 100 I=1,(NSP-1) 230 : YC.YET(I)=MPYC.VPOCHA(NLC,I) 231 : 100 CONTINUE 232 : C--------------------------------- 233 : C Variables à la face 234 : C--------------------------------- 235 : RF=MPLIM.VPOCHA(NLCB,1) 236 : UXF=MPLIM.VPOCHA(NLCB,2) 237 : UYF=MPLIM.VPOCHA(NLCB,3) 238 : PF=MPLIM.VPOCHA(NLCB,IDIM+2) 239 : SEGINI YF 240 : SEGACT MPLIM 241 : DO 101 I=1,(NSP-1) 242 : YF.YET(I)=MPLIM.VPOCHA(NLCB,IDIM+2+I) 243 : 101 CONTINUE 244 : c------------------------------------------------------------- 245 : c Computing GAMMA at the cell-center 246 : c------------------------------------------------------------- 247 : top=0.0D0 248 : bot=0.0D0 249 : do 102 i=1,(nsp-1) 250 : top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp)) 251 : bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp)) 252 : 102 continue 253 : top=cp.gc(nsp)+top 254 : bot=cv.gc(nsp)+bot 255 : GAMC=top/bot 256 : C------------------------------------------------------------- 257 : SEGINI DGDYC 258 : do 41 i=1,(nsp-1) 259 : dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)- 260 : & GAMC*(cv.gc(i)-cv.gc(nsp)))/bot 261 : 41 continue 262 : c------------------------------------------------------------- 263 : c Computing GAMMA at the face-center 264 : c------------------------------------------------------------- 265 : top=0.0D0 266 : bot=0.0D0 267 : do 103 i=1,(nsp-1) 268 : top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp)) 269 : bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp)) 270 : 103 continue 271 : top=cp.gc(nsp)+top 272 : bot=cv.gc(nsp)+bot 273 : GAMF=top/bot 274 : C-------------------------------------------------------------- 275 : C******* On calcule UN, UT ASON, S 276 : C-------------------------------------------------------------- 277 : UNC=(UXC*CNX)+(UYC*CNY) 278 : UNF=(UXF*CNX)+(UYF*CNY) 279 : UTF=(UXF*CTX)+(UYF*CTY) 280 : C--------------------------------- 281 : ASONC=(GAMC*PC/RC)**0.5D0 282 : ASONF=(GAMF*PF/RF)**0.5D0 283 : C--------------------------------- 284 : SF=PF/(RF**GAMF) 285 : C------------------------------------------------ 286 : C******* Densite, vitesse, pression sur le bord 287 : C------------------------------------------------ 288 : G1=UNC-(2.0D0*ASONC)/(GAMC-1.0D0) 289 : G3=UNF+(2.0D0*ASONF)/(GAMF-1.0D0) 290 : UN=0.5D0*(G1+G3) 291 : ASON=(0.5D0*(G3-G1)) 292 : ASON=ASON*(GAMF-1.0D0)/2.0D0 293 : ASON2=ASON*ASON 294 : S=SF 295 : UT=UTF 296 : RHO=ASON2/(GAMF*S) 297 : RHO=RHO**(1.0D0/(GAMF-1.0D0)) 298 : P=RHO*ASON2/GAMF 299 : UX=(UN*CNX)+(UT*CTX) 300 : UY=(UN*CNY)+(UT*CTY) 301 : C------------------------------ 302 : C******* Derivatives 303 : C------------------------------ 304 : SEGINI DYDG1 305 : DUNDG1=0.5D0 306 : DASDG1=-0.5D0*(GAMF-1.0D0)/2.0D0 307 : DRHDG1=GAMF*S 308 : DRHDG1=1.0D0/DRHDG1 309 : DRHDG1=DRHDG1**(1.0D0/(GAMF-1.0D0)) 310 : DRHDG1=DRHDG1*2.0D0/(GAMF-1.0D0) 311 : DRHDG1=DRHDG1*(ASON**((3.0D0-GAMF)/(GAMF-1.0D0))) 312 : DRHDG1=DRHDG1*DASDG1 313 : DPDG1=((ASON2/GAMF)*DRHDG1)+(((2*ASON*RHO)/GAMF)*DASDG1) 314 : DO 104 I=1,(NSP-1) 315 : DYDG1.VV(I) = 0.0D0 316 : 104 CONTINUE 317 : CC------------------------------------------------------ 318 : CC******* Test 319 : CC------------------------------------------------------ 320 : c USGM1 = 1.0D0/(GAMF-1.0D0) 321 : c DFEDG1=(DUNDG1*GAMF*USGM1*P) + (UN*GAMF*USGM1*DPDG1) + 322 : cc & (0.5D0*DRHDG1*UN*((UN*UN)+(UT*UT))) + 323 : c & (0.5D0*RHO*DUNDG1*((UN*UN)+(UT*UT))) + 324 : c & (RHO*UN*UN*DUNDG1) 325 : c CACCA=UN*GAMF*USGM1*P + 0.5D0*RHO*UN*(UN*UN+UT*UT) 326 : cc--------------- 327 : c DFRYG1.VV(1) = (DRHDG1*UN*YF.YET(1))+ 328 : c & (RHO*DUNDG1*YF.YET(1))+(RHO*UN*DYDG1.VV(1)) 329 : c CACC1 =RHO*UN*YF.YET(1) 330 : cc--------------- 331 : c EPS=1.0D-8 332 : c G1=G1*(1+EPS) 333 : c CEL=UN 334 : c UN=0.5D0*(G1+G3) 335 : c write(*,*) DUNDG1, ((UN - CEL)/(EPS*G1)) 336 : c CEL=ASON 337 : c ASON=(0.5D0*(G3-G1)) 338 : c ASON=ASON*(GAMM-1.0D0)/2.0D0 339 : c write(*,*) DASDG1, ((ASON - CEL)/(EPS*G1)) 340 : c ASON2=ASON*ASON 341 : c S=SF 342 : c UT=UTF 343 : c CEL=RHO 344 : c RHO=ASON2/(GAMF*S) 345 : c RHO=RHO**(1.0D0/(GAMF-1.0D0)) 346 : c write(*,*) DRHDG1, ((RHO - CEL)/(EPS*G1)) 347 : c CEL=P 348 : c P=RHO*ASON2/GAMF 349 : c write(*,*) DPDG1, ((P - CEL)/(EPS*G1)) 350 : c CEL=CACCA 351 : c CACCA=UN*GAMF*USGM1*P + 0.5D0*RHO*UN*(UN*UN+UT*UT) 352 : c write(*,*) DFEDG1, ((CACCA - CEL)/(EPS*G1)) 353 : c CEL=CACC1 354 : c CACCA =RHO*UN*YF.YET(1) 355 : c write(*,*) DFRYG1.VV(1), ((CACCA - CEL)/(EPS*G1)) 356 : CC 357 : CC************************************************************ 358 : CC 359 : DUXDG1=DUNDG1*CNX 360 : DUYDG1=DUNDG1*CNY 361 : C-------------------------------------------------------------- 362 : C ------- Be carefull here for outlet BC !!!!!!!!!!!!!!!!!!!! 363 : C-------------------------------------------------------------- 364 : SEGINI DFRYG1 365 : USGM1 = 1.0D0/(GAMF-1.0D0) 366 : DFRDG1=(DRHDG1*UN)+(RHO*DUNDG1) 367 : DFMXG1=(DRHDG1*UN*UX)+(RHO*DUNDG1*UX)+ 368 : & (RHO*UN*DUXDG1)+(DPDG1*CNX) 369 : DFMYG1=(DRHDG1*UN*UY)+(RHO*DUNDG1*UY)+ 370 : & (RHO*UN*DUYDG1)+(DPDG1*CNY) 371 : DFEDG1=(DUNDG1*GAMF*USGM1*P) + (UN*GAMF*USGM1*DPDG1) + 372 : & (0.5D0*DRHDG1*UN*((UN*UN)+(UT*UT))) + 373 : & (0.5D0*RHO*DUNDG1*((UN*UN)+(UT*UT))) + 374 : & (RHO*UN*UN*DUNDG1) 375 : DO 105 I=1,(NSP-1) 376 : DFRYG1.VV(I) = (DRHDG1*UN*YF.YET(I))+ 377 : & (RHO*DUNDG1*YF.YET(I))+(RHO*UN*DYDG1.VV(I)) 378 : 105 CONTINUE 379 : C----------------------------------------------------- 380 : C******* Jacobian with respect to primitive variables 381 : C----------------------------------------------------- 382 : SEGINI DG1DY 383 : USGM1 = 1.0D0/(GAMC-1.0D0) 384 : DG1DR=USGM1*ASONC/RC 385 : DG1DP=-1.0D0*USGM1*ASONC/PC 386 : DG1DUX=CNX 387 : DG1DUY=CNY 388 : DACDG=0.5D0*ASONC/GAMC 389 : DO 106 I=1,(NSP-1) 390 : DG1DY.VV(I)=2.0D0*USGM1*(ASONC*USGM1-DACDG)* 391 : & dgdyc.vv(i) 392 : 106 CONTINUE 393 : C---------------------------------------- 394 : COEF=SURF/VOLU 395 : C---------------------------------------- 396 : JTL.JAC(1,1) = DFRDG1*DG1DR*COEF 397 : JTL.JAC(1,2) = DFRDG1*DG1DUX*COEF 398 : JTL.JAC(1,3) = DFRDG1*DG1DUY*COEF 399 : JTL.JAC(1,4) = DFRDG1*DG1DP*COEF 400 : DO 107 I=1,(NSP-1) 401 : JTL.JAC(1,4+I) = DFRDG1*DG1DY.VV(I)*COEF 402 : 107 CONTINUE 403 : c CACCA = G1 404 : c top=0.0D0 405 : c bot=0.0D0 406 : c yc.yet(1) = yc.yet(1)+EPS 407 : c do 1023 i=1,(nsp-1) 408 : c top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp)) 409 : c bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp)) 410 : c 1023 continue 411 : c top=cp.gc(nsp)+top 412 : c bot=cv.gc(nsp)+bot 413 : c GAMC=top/bot 414 : c ASONC=(GAMC*PC/RC)**(0.5D0) 415 : c GAMM = 0.5D0*(GAMC+GAMF) 416 : c G1=UNC-(2.0D0*ASONC)/(GAMC-1.0D0) 417 : c G3=UNF+(2.0D0*ASONF)/(GAMF-1.0D0) 418 : c UN=0.5D0*(G1+G3) 419 : c ASON=(0.5D0*(G3-G1)) 420 : c ASON=ASON*(GAMM-1.0D0)/2.0D0 421 : c ASON2=ASON*ASON 422 : c S=SF 423 : c UT=UTF 424 : c RHO=ASON2/(GAMF*S) 425 : c RHO=RHO**(1.0D0/(GAMF-1.0D0)) 426 : c CEL=CACCA 427 : c CACCA=G1 428 : C---------------------------------------- 429 : JTL.JAC(2,1) = DFMXG1*DG1DR*COEF 430 : JTL.JAC(2,2) = DFMXG1*DG1DUX*COEF 431 : JTL.JAC(2,3) = DFMXG1*DG1DUY*COEF 432 : JTL.JAC(2,4) = DFMXG1*DG1DP*COEF 433 : DO 108 I=1,(NSP-1) 434 : JTL.JAC(2,4+I) = DFMXG1*DG1DY.VV(I)*COEF 435 : 108 CONTINUE 436 : C---------------------------------------- 437 : JTL.JAC(3,1) = DFMYG1*DG1DR*COEF 438 : JTL.JAC(3,2) = DFMYG1*DG1DUX*COEF 439 : JTL.JAC(3,3) = DFMYG1*DG1DUY*COEF 440 : JTL.JAC(3,4) = DFMYG1*DG1DP*COEF 441 : DO 109 I=1,(NSP-1) 442 : JTL.JAC(3,4+I) = DFMYG1*DG1DY.VV(I)*COEF 443 : 109 CONTINUE 444 : C---------------------------------------- 445 : JTL.JAC(4,1) = DFEDG1*DG1DR*COEF 446 : JTL.JAC(4,2) = DFEDG1*DG1DUX*COEF 447 : JTL.JAC(4,3) = DFEDG1*DG1DUY*COEF 448 : JTL.JAC(4,4) = DFEDG1*DG1DP*COEF 449 : DO 110 I=1,(NSP-1) 450 : JTL.JAC(4,4+I) = DFEDG1*DG1DY.VV(I)*COEF 451 : 110 CONTINUE 452 : C---------------------------------------- 453 : DO 111 I=1,(NSP-1) 454 : JTL.JAC(4+I,1) = DFRYG1.VV(I)*DG1DR*COEF 455 : JTL.JAC(4+I,2) = DFRYG1.VV(I)*DG1DUX*COEF 456 : JTL.JAC(4+I,3) = DFRYG1.VV(I)*DG1DUY*COEF 457 : JTL.JAC(4+I,4) = DFRYG1.VV(I)*DG1DP*COEF 458 : DO 112 J=1,(NSP-1) 459 : JTL.JAC(4+I,4+J) = DFRYG1.VV(I)*DG1DY.VV(J)*COEF 460 : 112 CONTINUE 461 : 111 CONTINUE 462 : C--------------------------------------------------------- 463 : c matrix wl(i,j) represents the derivative of the i-component 464 : c of the vector of primitive variables of the left state with 465 : c respect to the j-component of the vector of the conservative 466 : c variables of the left state. 467 : c 468 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) - 469 : c vector of primitive variables; 470 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) - 471 : c vector of conservative variables. 472 : c------------------------------------------------------------- 473 : wl.jac(1,1)=1.0d0 474 : wl.jac(1,2)=0.0d0 475 : wl.jac(1,3)=0.0d0 476 : wl.jac(1,4)=0.0d0 477 : do 83 i=1,(nsp-1) 478 : wl.jac(1,4+i)=0.0d0 479 : 83 continue 480 : c------------------------------ 481 : wl.jac(2,1)=-UXC/RC 482 : wl.jac(2,2)=1.0d0/RC 483 : wl.jac(2,3)=0.0d0 484 : wl.jac(2,4)=0.0d0 485 : do 84 i=1,(nsp-1) 486 : wl.jac(2,4+i)=0.0d0 487 : 84 continue 488 : c------------------------------ 489 : wl.jac(3,1)=-UYC/RC 490 : wl.jac(3,2)=0.0d0 491 : wl.jac(3,3)=1.0d0/RC 492 : wl.jac(3,4)=0.0d0 493 : do 85 i=1,(nsp-1) 494 : wl.jac(3,4+i)=0.0d0 495 : 85 continue 496 : c------------------------------ 497 : br1=0.0d0 498 : do 86 i=1,(nsp-1) 499 : br1=br1+dgdyc.vv(i)*yc.yet(i) 500 : 86 continue 501 : br1=br1*PC/(RC*(GAMC-1.0D0)) 502 : wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1 503 : wl.jac(4,2)=-UXC*(GAMC-1.0D0) 504 : wl.jac(4,3)=-UYC*(GAMC-1.0D0) 505 : wl.jac(4,4)=(GAMC-1.0D0) 506 : do 87 i=1,(nsp-1) 507 : wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0)) 508 : 87 continue 509 : c------------------------------ 510 : do 88 i=1,(nsp-1) 511 : do 89 j=1,4 512 : wl.jac(4+i,j)=0.0d0 513 : if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC 514 : 89 continue 515 : c------------ 516 : do 890 j=5,(4+nsp-1) 517 : wl.jac(4+i,j)=0.0d0 518 : if(4+i.eq.j) then 519 : wl.jac(4+i,j)=1.0d0/RC 520 : endif 521 : 890 continue 522 : 88 continue 523 : c------------------------------------------------ 524 : C------------------------------------------------ 525 : do 114 i=1,(3+nsp) 526 : do 115 j=1,(3+nsp) 527 : jtt.jac(i,j)=0.0d0 528 : do 116 k=1,(3+nsp) 529 : jtt.jac(i,j)=jtt.jac(i,j)+jtl.jac(i,k)*wl.jac(k,j) 530 : 116 continue 531 : 115 continue 532 : 114 continue 533 : C---------------------------------------------------------------- 534 : C******* Jacobian with respect to conservative variables 535 : C---------------------------------------------------------------- 536 : IF(IJAC.EQ.1)THEN 537 : DO 9 II = 1,(3+NSP) 538 : DO 15 JJ = 1,(3+NSP) 539 : KV = (II-1)*(3+NSP) 540 : C---------------------------------- 541 : CELL = IMATRI.LIZAFM(1,KV+JJ) 542 : CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ) 543 : 15 CONTINUE 544 : 9 CONTINUE 545 : ELSEIF(IJAC.EQ.2)THEN 546 : DO 20 II = 1,(3+NSP) 547 : DO 25 JJ = 1,(3+NSP) 548 : KV = (II-1)*(3+NSP) 549 : C---------------------------------- 550 : CELL = IMATRI.LIZAFM(1,KV+JJ) 551 : CELL.AM(IFAC,1,1) = JTL.JAC(II,JJ) 552 : 25 CONTINUE 553 : 20 CONTINUE 554 : ENDIF 555 : c-------------------------------------------------- 556 : ENDDO 557 : C 558 : SEGDES MELEFC 559 : C 560 : SEGSUP MLEMC 561 : SEGSUP MLEMCB 562 : SEGSUP MLEMF 563 : C 564 : SEGDES MPNORM 565 : SEGDES MPVOL 566 : SEGDES MPSURF 567 : SEGDES MPRC 568 : SEGDES MPPC 569 : SEGDES MPVC 570 : SEGDES MPYC 571 : SEGDES MPLIM 572 : SEGDES YC 573 : SEGDES YF 574 : SEGDES CP 575 : SEGDES CV 576 : SEGDES JTL 577 : SEGDES JTT 578 : SEGDES WL 579 : SEGDES DYDG1, DFRYG1, 580 : & DG1DY, DGDYC 581 : SEGDES MATRIK 582 : DO 80 II=1,NBME 583 : CELL = IMATRI.LIZAFM(1,II) 584 : SEGDES CELL 585 : 80 CONTINUE 586 : SEGDES IMATRI 587 : C--------------------------------------------- 588 : 9999 CONTINUE 589 : RETURN 590 : END 591 : 592 : 593 : 594 : 595 : 596 : 597 : 598 : 599 : 600 :
© Cast3M 2003 - All rights reserved.
Disclaimer