1 : C KONJA5 SOURCE MAGN 09/07/30 21:17:53 6458 2 : SUBROUTINE KONJA5(IIMPL,ILINC,IROF,IPF,IVITF,IGAMF, 3 : & IRN,IPN,IVN,IGAMN, 4 : $ ICHPVO,ICHPSU, ICHPNO, MELEMC, MELEFE, MELEMF, MELTFA, 5 : & IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN, 6 : & ICRN, ICVN, ICPN, IVLIM,MELLIM,IMAT) 7 : C 8 : C************************************************************************ 9 : C 10 : C PROJET : CASTEM 2000 11 : C 12 : C NOM : KONJA5 13 : C 14 : C DESCRIPTION : Voir KONV11 15 : C Calcul du jacobien du résidu "2nd order accurate" 16 : C 16 : C Cas deux dimensions, gaz "calorically perfect" 17 : C 18 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 19 : C 20 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 21 : C 22 : C************************************************************************ 23 : C 24 : C APPELES (Outils 25 : C CASTEM) : KRIPAD, LICHT, ERREUR 26 : C 27 : C APPELES (Calcul) : VLHJ3, AUSMW, VLHJ1, AUSMF 28 : C 29 : C************************************************************************ 30 : C 31 : C ENTREES 32 : C 33 : C IIMPL : INTEGER 34 : C = 4 -> methode VLH 35 : C = 5 -> methode AUSM+ 36 : C 37 : C ILINC : liste des inconnues (pointeur d'un objet de type 38 : C LISTMOTS) 39 : C 40 : C 1) Pointeurs des CHPOINT/CHAMELEM 41 : C 42 : C IROF : MCHAML 'FACEL' densité 43 : C 44 : C IPF : MCHAML 'FACEL' pression 45 : C 46 : C IVITF : MCHAML 'FACEL' vitesse/normales 47 : C 48 : C IGAMF : MCHAML 'FACEL' gamma 49 : C 50 : C IRN : CHPOINT 'CENTRE' densité 51 : C 52 : C IPN : CHPOINT 'CENTRE' pression 53 : C 54 : C IVN : CHPOINT 'CENTRE' vitesse 55 : C 56 : C IGAMN : CHPOINT 'CENTRE' gamma 57 : C 58 : C ICHPVO : CHPOINT VOLUME contenant le volume 59 : C 60 : C ICHPSU : CHPOINT FACE contenant la surface des faces 61 : C 62 : C ICHPNO : CHPOINT 'FACE' contenant les normales aux faces 63 : C 64 : C IGRN : CHPOINT 'CENTRE' gradient de densité 65 : C 66 : C IGVN : CHPOINT 'CENTRE' gradient de vitesse 67 : C 68 : C IGPN : CHPOINT 'CENTRE' gradient de pression 69 : C 70 : C IGLRN : CHPOINT 'CENTRE' limiteur du gradient de densité 71 : C 72 : C IGLVN : CHPOINT 'CENTRE' limiteur du gradient de vitesse 73 : C 74 : C IGLPN : CHPOINT 'CENTRE' limiteur du gradient de pression 75 : C 76 : C ICRN : MCHAML contenant les coeff. pour calculer le gradient 77 : C de densité 78 : C 79 : C ICVN : MCHAML contenant les coeff. pour calculer le gradient 80 : C de vitesse 81 : C 82 : C ICPN : MCHAML contenant les coeff. pour calculer le gradient 83 : C de pression 84 : C 85 : C IVLIM : CHPOINT des conditions aux limites sur la vitesse 86 : C 87 : C 2) Pointeurs de MELEME de la table DOMAINE 88 : C 89 : C MELEMC : MELEME 'CENTRE', SPG des CENTRES des elts 90 : C 91 : C MELEFE : MELEME 'FACEL' (centre elt gauche, centre de face 92 : C centre elt droite) 93 : C 94 : C MELEMF : MELEME 'FACE', SPG des CENTRES des faces 95 : C 96 : C MELTFA : MELEME 'ELTFA' (centre d'elt, centres de faces) 97 : C 98 : C MELLIM : MELEME SPG des conditions aux bords 99 : C 100 : C SORTIES 101 : C 102 : C IMAT : pointeur de la MATRIK du jacobien du residu 103 : C 104 : C************************************************************************ 105 : C 106 : C HISTORIQUE (Anomalies et modifications éventuelles) 107 : C 108 : C HISTORIQUE : 109 : C 110 : C************************************************************************ 111 : C 112 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0 113 : C GAMMA \in (1,3) 114 : C Si non il faut le faire!!! 115 : C 116 : C************************************************************************ 117 : C 118 : C 119 : C**** Variables de COOPTIO 120 : C 121 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 122 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 123 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 124 : C & ,IECHO, IIMPI, IOSPI 125 : C & ,IDIM, IFICLE, IPREFI 126 : CC & ,MCOORD 127 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 128 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 129 : C & ,NORINC,NORVAL,NORIND,NORVAD 130 : C & ,NUCROU, IPSAUV 131 : C 132 : IMPLICIT INTEGER(I-N) 133 : INTEGER IIMPL,ILINC, ICHPVO, ICHPNO, ICHPSU 134 : & ,IROF,IVITF,IPF,IGAMF, IRN,IPN,IVN,IGAMN 135 : & , IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN 136 : & , ICRN, ICVN, ICPN, IVLIM , IMAT 137 : & , NBSOUP, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID 138 : & , NKMT, NBME, NBEL, MP, NP, NP1 139 : & , IFAC, NGCF, NLCF, NGCF1, NGCG, NGCD 140 : & , NGC, NLC, NLCD, NGNO, NGN1,NGN2, NLNO, IELEM, IELPRI 141 : & , JG, NBSOUS, NBSOU1, ISOUS, NBEL1 142 : & , ISPG1, ISPG2 143 : & , I1, IGEOMC, IGEOMF, IGEOML 144 : & , IFACGR, NGGRA, NLGRA, NGFGR, NLFGR, NG, ND, ICOOR, IFAC2 145 : & , NLVIMP 146 : & , NGPRI,NLPRI,IELETR,IDUATR,IPRITR, NLFL 147 : REAL*8 ROG, PG, UNG, UTG, UXG, UYG, RETG, GAMG, VOLG, VOLD 148 : & , ROD, PD, UND, UTD, UXD, UYD 149 : & , SURF, CNX, CNY, CTX, CTY, ECIN, FUNCEL 150 : & , DFRO(4), DFRET(4), DFRUN(4), DFRUT(4),DFRUX(4),DFRUY(4) 151 : & , ALR, ALP, ALVX, ALVY, XCG, YCG, DXCF, DYCF, CELLR, CELLP 152 : & , CVXVX, CVXVY, CVYVX, CVYVY, CNX1, CNY1, CTX1, CTY1 153 : & , DDP, DDRO, DDUX, DDUY 154 : & , GAMGM1, CELCO1, CELCO2, CELCO3, CELCO4, CELCO5 155 : 156 : CHARACTER*8 TYPE 157 : C 158 : C**** LES INCLUDES 159 : C 160 : -INC SMCOORD 161 : -INC CCOPTIO 162 : -INC SMCHAML 163 : POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, 164 : & MELT1X.MELVAL, MELT1Y.MELVAL 165 : POINTEUR MELVUN.MELVAL, MELVUT.MELVAL 166 : POINTEUR MELRO.MELVAL, MELP.MELVAL, 167 : & MELGAM.MELVAL 168 : POINTEUR MELCR1.MELVAL, MELCR2.MELVAL, MELCP1.MELVAL, 169 : & MELCP2.MELVAL, MELCV1.MELVAL, MELCV2.MELVAL 170 : C 171 : -INC SMCHPOI 172 : -INC SMELEME 173 : -INC SMLMOTS 174 : -INC SMLENTI 175 : -INC SMMATRIK 176 : POINTEUR MPVOLU.MPOVAL, MPNORM.MPOVAL, MPOVSU.MPOVAL, 177 : $ MPVGRN.MPOVAL,MPGLRN.MPOVAL, MPVGPN.MPOVAL, MPGLPN.MPOVAL, 178 : $ MPVGVN.MPOVAL,MPGLVN.MPOVAL, MPVLIM.MPOVAL,MPOVRN.MPOVAL,MPOVVN 179 : $ .MPOVAL,MPOVPN.MPOVAL,MPOGAM.MPOVAL 180 : POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME, 181 : & MELTFA.MELEME,MELPRI.MELEME,MELGRR.MELEME,MELGRP.MELEME 182 : $ ,MELGRV.MELEME, MELLIM.MELEME 183 : POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLESPG.MLENTI,MLEPRI.MLENTI 184 : & ,MLEVL.MLENTI,MLEGRR.MLENTI,MLEGRP.MLENTI,MLEGRV.MLENTI 185 : & ,MLELIM.MLENTI 186 : POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RRET.IZAFM, 187 : & UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXRET.IZAFM, 188 : & UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYRET.IZAFM, 189 : & RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETRET.IZAFM 190 : POINTEUR MLMINC.MLMOTS 191 : C 192 : C**** KRIPAD pour la correspondance global/local des conditions limits 193 : C 193 : CALL KRIPAD(MELLIM,MLELIM) 194 : C SEGACT MELLIM 195 : C 196 : C**** Masse volumique 197 : C 198 : MCHEL1 = IROF 199 : SEGACT MCHEL1 200 : MCHAM1 = MCHEL1.ICHAML(1) 201 : SEGACT MCHAM1 202 : MELRO = MCHAM1.IELVAL(1) 203 : SEGDES MCHEL1 204 : SEGDES MCHAM1 205 : C 206 : C**** Pression 207 : C 208 : MCHEL1 = IPF 209 : SEGACT MCHEL1 210 : MCHAM1 = MCHEL1.ICHAML(1) 211 : SEGACT MCHAM1 212 : MELP = MCHAM1.IELVAL(1) 213 : SEGDES MCHEL1 214 : SEGDES MCHAM1 215 : C 216 : C**** Gamma 217 : C 218 : MCHEL1 = IGAMF 219 : SEGACT MCHEL1 220 : MCHAM1 = MCHEL1.ICHAML(1) 221 : SEGACT MCHAM1 222 : MELGAM = MCHAM1.IELVAL(1) 223 : SEGDES MCHEL1 224 : SEGDES MCHAM1 225 : C 226 : C**** Vitesse et cosinus directeurs du repere (n,t) 227 : C 228 : MCHEL1 = IVITF 229 : SEGACT MCHEL1 230 : C 231 : C**** La vitesse a comme SPG MELEFE 232 : C Le cosinus directeurs ont comme SPG MELEMF 233 : C 234 : C MCHAM1 -> Cosinus directeurs 235 : C MCHAM2 -> Vitesse 236 : C 237 : ISPG1 = MCHEL1.IMACHE(1) 238 : ISPG2 = MCHEL1.IMACHE(2) 239 : IF((ISPG1 .EQ. MELEMF) .AND. (ISPG2 .EQ. MELEFE))THEN 240 : MCHAM1 = MCHEL1.ICHAML(1) 241 : MCHAM2 = MCHEL1.ICHAML(2) 242 : ELSEIF((ISPG1 .EQ. MELEFE) .AND. (ISPG2 .EQ. MELEMF))THEN 243 : MCHAM1 = MCHEL1.ICHAML(2) 244 : MCHAM2 = MCHEL1.ICHAML(1) 245 : ELSE 246 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 246 : CALL ERREUR(5) 247 : GOTO 9999 248 : ENDIF 249 : SEGACT MCHAM1 250 : MELVNX = MCHAM1.IELVAL(1) 251 : MELVNY = MCHAM1.IELVAL(2) 252 : MELT1X = MCHAM1.IELVAL(3) 253 : MELT1Y = MCHAM1.IELVAL(4) 254 : SEGDES MCHAM1 255 : SEGACT MCHAM2 256 : MELVUN = MCHAM2.IELVAL(1) 257 : MELVUT = MCHAM2.IELVAL(2) 258 : SEGDES MCHAM2 259 : SEGDES MCHEL1 260 : C 261 : C**** Activation des MCHAMLs 262 : C 263 : SEGACT MELRO 264 : SEGACT MELP 265 : SEGACT MELGAM 266 : SEGACT MELVUN 267 : SEGACT MELVUT 268 : SEGACT MELVNX 269 : SEGACT MELVNY 270 : SEGACT MELT1X 271 : SEGACT MELT1Y 272 : C 273 : C**************** CHPOINT gradients, limiteurs, vitesse aux bords 274 : C 274 : CALL LICHT(IRN,MPOVRN,TYPE,IGEOMC) 275 : IF(IERR.NE.0)GOTO 9999 276 : C SEGACT MPOVRN 276 : CALL LICHT(IVN,MPOVVN,TYPE,IGEOMC) 277 : IF(IERR.NE.0)GOTO 9999 278 : C SEGACT MPOVVN 278 : CALL LICHT(IPN,MPOVPN,TYPE,IGEOMC) 279 : IF(IERR.NE.0)GOTO 9999 280 : C SEGACT MPOVPN 280 : CALL LICHT(IGAMN,MPOGAM,TYPE,IGEOMC) 281 : IF(IERR.NE.0)GOTO 9999 282 : C SEGACT MPOGAM 282 : CALL LICHT(IGRN,MPVGRN,TYPE,IGEOMC) 283 : IF(IERR.NE.0)GOTO 9999 284 : C SEGACT MPVGRN 284 : CALL LICHT(IGLRN,MPGLRN,TYPE,IGEOMC) 285 : IF(IERR.NE.0)GOTO 9999 286 : C SEGACT MPGLRN 286 : CALL LICHT(IGPN,MPVGPN,TYPE,IGEOMC) 287 : IF(IERR.NE.0)GOTO 9999 288 : C SEGACT MPVGPN 288 : CALL LICHT(IGLPN,MPGLPN,TYPE,IGEOMC) 289 : IF(IERR.NE.0)GOTO 9999 290 : C SEGACT MPGLPN 290 : CALL LICHT(IGVN,MPVGVN,TYPE,IGEOMC) 291 : IF(IERR.NE.0)GOTO 9999 292 : C SEGACT MPVGPN 292 : CALL LICHT(IGLVN,MPGLVN,TYPE,IGEOMC) 293 : IF(IERR.NE.0)GOTO 9999 294 : C SEGACT MPGLVN 294 : CALL LICHT(IVLIM,MPVLIM,TYPE,IGEOML) 295 : IF(IERR.NE.0)GOTO 9999 296 : C SEGACT MPVLIM 297 : C 298 : C**************** KRIPAD pour la correspondance global/local des centres d'elts 299 : C 299 : CALL KRIPAD(MELEMC,MLENTC) 300 : IF(IERR.NE.0)GOTO 9999 301 : C 302 : C SEGACT MLENTC 303 : SEGACT MELEMC 304 : C 305 : C**** KRIPAD pour la correspondance global/local des centres des faces 306 : C 306 : CALL KRIPAD(MELEMF,MLENTF) 307 : IF(IERR.NE.0)GOTO 9999 308 : C 309 : C SEGACT MLENTF 310 : C 311 : C**** KRIPAD pour la correspondance global/local des B.C. sur le vitesse 312 : C 312 : CALL KRIPAD(IGEOML,MLEVL) 313 : IF(IERR.NE.0)GOTO 9999 314 : C 315 : C SEGACT MLEVL 316 : C 317 : SEGACT MELEFE 318 : C 318 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF) 319 : IF(IERR.NE.0)GOTO 9999 319 : CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC) 320 : IF(IERR.NE.0)GOTO 9999 320 : CALL LICHT(ICHPNO,MPNORM,TYPE,IGEOMC) 321 : IF(IERR.NE.0)GOTO 9999 322 : C 323 : C**** LICHT active les MPOVALs en *MOD 324 : C 325 : C i.e. 326 : C 327 : C SEGACT MPOVSU*MOD 328 : C SEGACT MPVOLU*MOD 329 : C SEGACT MPNORM*MOD 330 : C 331 : C 332 : C**** Les gradients ont tous le meme SPG (mais le numero 333 : C du pointeur peut-etre different). 334 : C Ils ont le meme nombre des SPGs que MELTFA 335 : C En plus le SPGs ont le meme ordre: 336 : C le iem elt de MELTFA correspond a l'elt qui a 337 : C comme centre le iem noeud de CENTRE 338 : C 339 : SEGACT MELTFA 340 : NBSOUP=MELTFA.LISOUS(/1) 341 : IF(NBSOUP.EQ.0) NBSOUP=1 342 : JG=NBSOUP 343 : SEGINI MLESPG 344 : IF(NBSOUP .EQ. 1)THEN 345 : MLESPG.LECT(1)=MELTFA 346 : ELSE 347 : DO I1 = 1, NBSOUP, 1 348 : MLESPG.LECT(I1)=MELTFA.LISOUS(I1) 349 : ENDDO 350 : ENDIF 351 : C 352 : MCHEL1=ICRN 353 : SEGACT MCHEL1 354 : NBSOU1=MCHEL1.IMACHE(/1) 355 : IF(NBSOU1 .NE. NBSOUP)THEN 356 : WRITE(IOIMP,*) 'SPG coeffs grads = ???' 357 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 357 : CALL ERREUR(5) 358 : GOTO 9999 359 : ENDIF 360 : C 361 : MCHEL2=ICPN 362 : SEGACT MCHEL2 363 : NBSOU1=MCHEL2.IMACHE(/1) 364 : IF(NBSOU1 .NE. NBSOUP)THEN 365 : WRITE(IOIMP,*) 'SPG coeffs grads = ???' 366 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 366 : CALL ERREUR(5) 367 : GOTO 9999 368 : ENDIF 369 : C 370 : MCHEL2=ICVN 371 : SEGACT MCHEL2 372 : NBSOU1=MCHEL2.IMACHE(/1) 373 : IF(NBSOU1 .NE. NBSOUP)THEN 374 : WRITE(IOIMP,*) 'SPG coeffs grads = ???' 375 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 375 : CALL ERREUR(5) 376 : GOTO 9999 377 : ENDIF 378 : C 379 : C**** N.B.: ICRN, ICPN, ICVN sont activés!!! 380 : C 381 : C 382 : C****************** L'objet matrik 383 : C 384 : C 385 : NRIGE = 7 386 : NMATRI = 1 387 : NKID = 9 388 : NKMT = 7 389 : C 390 : SEGINI MATRIK 391 : IMAT = MATRIK 392 : C 393 : C**** Maillage des inconnues duales: MELEMC 394 : C Maillage des inconnues primales: 395 : C même structure que MCHEL1.IMACHE(ISOUS) 396 : C Donc on partitionne les deux dans la meme 397 : C maniere. 398 : C 399 : JG=NBSOUP 400 : SEGINI MLEPRI 401 : SEGINI MLEGRR 402 : SEGINI MLEGRP 403 : SEGINI MLEGRV 404 : IF(NBSOUP .EQ. 1)THEN 405 : MCHEL1=ICRN 406 : IPT1=MCHEL1.IMACHE(1) 407 : SEGACT IPT1 408 : SEGINI, MELPRI = IPT1 409 : MLEPRI.LECT(1) = MELPRI 410 : MLEGRR.LECT(1) = IPT1 411 : C 412 : MCHEL1=ICPN 413 : IPT1=MCHEL1.IMACHE(1) 414 : SEGACT IPT1 415 : MLEGRP.LECT(1) = IPT1 416 : C 417 : MCHEL1=ICVN 418 : IPT1=MCHEL1.IMACHE(1) 419 : SEGACT IPT1 420 : MLEGRV.LECT(1) = IPT1 421 : ELSE 422 : NBSOUS=NBSOUP 423 : NBREF=0 424 : NBNN=0 425 : NBELEM=0 426 : SEGINI MELPRI 427 : MCHEL1=ICRN 428 : MCHEL2=ICPN 429 : MCHEL3=ICVN 430 : DO ISOUS = 1, NBSOUP, 1 431 : IPT1=MCHEL1.IMACHE(ISOUS) 432 : IPT2=MCHEL2.IMACHE(ISOUS) 433 : IPT3=MCHEL3.IMACHE(ISOUS) 434 : SEGACT IPT1 435 : SEGACT IPT2 436 : SEGACT IPT3 437 : SEGINI, MELEME=IPT1 438 : MELPRI.LISOUS(ISOUS) = MELEME 439 : MLEPRI.LECT(ISOUS) = MELEME 440 : MLEGRR.LECT(ISOUS) = IPT1 441 : MLEGRP.LECT(ISOUS) = IPT2 442 : MLEGRV.LECT(ISOUS) = IPT3 443 : ENDDO 444 : ENDIF 445 : MATRIK.IRIGEL(2,1) = MELPRI 446 : MATRIK.IRIGEL(1,1) = MELPRI 447 : C 448 : C**** Matrice non symetrique 449 : C 450 : MATRIK.IRIGEL(7,1) = 2 451 : NBME = 16 452 : NBSOUS=NBSOUP 453 : SEGINI IMATRI 454 : MLMINC = ILINC 455 : SEGACT MLMINC 456 : MATRIK.IRIGEL(4,1) = IMATRI 457 : C 458 : IMATRI.LISPRI(1) = MLMINC.MOTS(1) 459 : IMATRI.LISPRI(2) = MLMINC.MOTS(2) 460 : IMATRI.LISPRI(3) = MLMINC.MOTS(3) 461 : IMATRI.LISPRI(4) = MLMINC.MOTS(4) 462 : IMATRI.LISPRI(5) = MLMINC.MOTS(1) 463 : IMATRI.LISPRI(6) = MLMINC.MOTS(2) 464 : IMATRI.LISPRI(7) = MLMINC.MOTS(3) 465 : IMATRI.LISPRI(8) = MLMINC.MOTS(4) 466 : IMATRI.LISPRI(9) = MLMINC.MOTS(1) 467 : IMATRI.LISPRI(10) = MLMINC.MOTS(2) 468 : IMATRI.LISPRI(11) = MLMINC.MOTS(3) 469 : IMATRI.LISPRI(12) = MLMINC.MOTS(4) 470 : IMATRI.LISPRI(13) = MLMINC.MOTS(1) 471 : IMATRI.LISPRI(14) = MLMINC.MOTS(2) 472 : IMATRI.LISPRI(15) = MLMINC.MOTS(3) 473 : IMATRI.LISPRI(16) = MLMINC.MOTS(4) 474 : C 475 : IMATRI.LISDUA(1) = MLMINC.MOTS(1) 476 : IMATRI.LISDUA(2) = MLMINC.MOTS(1) 477 : IMATRI.LISDUA(3) = MLMINC.MOTS(1) 478 : IMATRI.LISDUA(4) = MLMINC.MOTS(1) 479 : IMATRI.LISDUA(5) = MLMINC.MOTS(2) 480 : IMATRI.LISDUA(6) = MLMINC.MOTS(2) 481 : IMATRI.LISDUA(7) = MLMINC.MOTS(2) 482 : IMATRI.LISDUA(8) = MLMINC.MOTS(2) 483 : IMATRI.LISDUA(9) = MLMINC.MOTS(3) 484 : IMATRI.LISDUA(10) = MLMINC.MOTS(3) 485 : IMATRI.LISDUA(11) = MLMINC.MOTS(3) 486 : IMATRI.LISDUA(12) = MLMINC.MOTS(3) 487 : IMATRI.LISDUA(13) = MLMINC.MOTS(4) 488 : IMATRI.LISDUA(14) = MLMINC.MOTS(4) 489 : IMATRI.LISDUA(15) = MLMINC.MOTS(4) 490 : IMATRI.LISDUA(16) = MLMINC.MOTS(4) 491 : C 492 : C**** Boucle sur les LISOUS de MELPRI 493 : C 494 : IELEM=0 495 : C 496 : DO ISOUS=1,NBSOUP,1 497 : C 498 : C******* Les coeffs des grad. 499 : C 500 : MCHEL1=ICRN 501 : MCHAM1=MCHEL1.ICHAML(ISOUS) 502 : SEGACT MCHAM1 503 : MELCR1=MCHAM1.IELVAL(1) 504 : MELCR2=MCHAM1.IELVAL(2) 505 : SEGACT MELCR1, MELCR2 506 : SEGDES MCHAM1 507 : C 508 : MCHEL1=ICPN 509 : MCHAM1=MCHEL1.ICHAML(ISOUS) 510 : SEGACT MCHAM1 511 : MELCP1=MCHAM1.IELVAL(1) 512 : MELCP2=MCHAM1.IELVAL(2) 513 : SEGACT MELCP1, MELCP2 514 : SEGDES MCHAM1 515 : C 516 : C 517 : MCHEL1=ICVN 518 : MCHAM1=MCHEL1.ICHAML(ISOUS) 519 : SEGACT MCHAM1 520 : MELCV1=MCHAM1.IELVAL(1) 521 : MELCV2=MCHAM1.IELVAL(2) 522 : SEGACT MELCV1, MELCV2 523 : SEGDES MCHAM1 524 : C 525 : C******* MELEME : maillage elementaire de ELTFA 526 : C 527 : MELEME=MLESPG.LECT(ISOUS) 528 : SEGACT MELEME 529 : NBEL1=MELEME.NUM(/2) 530 : NP1=MELEME.NUM(/1) 531 : C 532 : C******* IPT1 : SPG des inconnues primales et duales 533 : C 534 : IPT1=MLEPRI.LECT(ISOUS) 535 : NBEL = IPT1.NUM(/2) 536 : NP = IPT1.NUM(/1) 537 : C 538 : C******* MELGRR: SPG elementaire des coeffs de GRAD 539 : C Il resemble a IPT1, mais IPT1 peut changer! 540 : C 541 : MELGRR=MLEGRR.LECT(ISOUS) 542 : MELGRP=MLEGRP.LECT(ISOUS) 543 : MELGRV=MLEGRV.LECT(ISOUS) 544 : C 545 : C******* MLESPG et MLEPRI doivent avoir le meme nombre 546 : C d'elts et de noeuds! 547 : C 548 : IF((NBEL1 .NE. NBEL) .OR. (NP1 .NE. (NP-1)))THEN 549 : WRITE(IOIMP,*) 'ELTFA, SPGs gradients' 550 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 550 : CALL ERREUR(5) 551 : GOTO 9999 552 : ENDIF 553 : C 554 : MP = NP 555 : SEGINI RR , RUX , RUY , RRET , 556 : & UXR , UXUX , UXUY , UXRET , 557 : & UYR , UYUX , UYUY , UYRET , 558 : & RETR , RETUX , RETUY , RETRET 559 : C 560 : C**** Duale = IMATRI.LISDUA(1) = 'RN' 561 : C Primale = IMATRI.LISPRI(1) = 'RN' 562 : C -> IMATRI.LIZAFM(1,1) = RR 563 : C 564 : C Duale = IMATRI.LISDUA(2) = 'RN' 565 : C Primale = IMATRI.LISPRI(2) = 'RUXN' 566 : C -> IMATRI.LIZAFM(1,2) = RUX 567 : C ... 568 : C 569 : C NB Pour le moment on y stoke le jacobien du flux 570 : C par rapport aux variables primitives (R,UX,UY,P) 571 : C 572 : IMATRI.LIZAFM(ISOUS,1) = RR 573 : IMATRI.LIZAFM(ISOUS,2) = RUX 574 : IMATRI.LIZAFM(ISOUS,3) = RUY 575 : IMATRI.LIZAFM(ISOUS,4) = RRET 576 : IMATRI.LIZAFM(ISOUS,5) = UXR 577 : IMATRI.LIZAFM(ISOUS,6) = UXUX 578 : IMATRI.LIZAFM(ISOUS,7) = UXUY 579 : IMATRI.LIZAFM(ISOUS,8) = UXRET 580 : IMATRI.LIZAFM(ISOUS,9) = UYR 581 : IMATRI.LIZAFM(ISOUS,10) = UYUX 582 : IMATRI.LIZAFM(ISOUS,11) = UYUY 583 : IMATRI.LIZAFM(ISOUS,12) = UYRET 584 : IMATRI.LIZAFM(ISOUS,13) = RETR 585 : IMATRI.LIZAFM(ISOUS,14) = RETUX 586 : IMATRI.LIZAFM(ISOUS,15) = RETUY 587 : IMATRI.LIZAFM(ISOUS,16) = RETRET 588 : C 589 : C******* Boucle sur les elts de MLEPRI.LECT 590 : C 591 : 592 : DO IELPRI=1,NBEL,1 593 : IELEM=IELEM+1 594 : NGC=IPT1.NUM(NP,IELPRI) 595 : NLC=MLENTC.LECT(NGC) 596 : IF(NLC.NE.IELEM)THEN 597 : WRITE(IOIMP,*) 'NLC != IELEM' 598 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 598 : CALL ERREUR(5) 599 : GOTO 9999 600 : ENDIF 601 : VOLG = MPVOLU.VPOCHA(NLC,1) 602 : C 603 : C************* Pour la contribution des gradients 604 : C 1) les limiteurs 605 : C 606 : ALR = MPGLRN.VPOCHA(NLC,1) 607 : ALP = MPGLPN.VPOCHA(NLC,1) 608 : ALVX = MPGLVN.VPOCHA(NLC,1) 609 : ALVY = MPGLVN.VPOCHA(NLC,2) 610 : ICOOR= (NGC - 1) * 3 611 : XCG = MCOORD.XCOOR(ICOOR + 1) 612 : YCG = MCOORD.XCOOR(ICOOR + 2) 613 : C 614 : C******** Boucle sur le noeuds de MLESPG.LECT 615 : C 616 : DO IFAC=1,(NP - 1),1 617 : NGCF = MELEME.NUM(IFAC,IELPRI) 618 : ICOOR= (NGCF - 1) * 3 619 : DXCF = MCOORD.XCOOR(ICOOR + 1) - XCG 620 : DYCF = MCOORD.XCOOR(ICOOR + 2) - YCG 621 : NLCF = MLENTF.LECT(NGCF) 622 : NGCF1 = MELEFE.NUM(2,NLCF) 623 : IF(NGCF .NE. NGCF1)THEN 624 : WRITE(IOIMP,*) 625 : $ 'FACE ET FACEL = ???' 626 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 626 : CALL ERREUR(5) 627 : GOTO 9999 628 : ENDIF 629 : NGCG = MELEFE.NUM(1,NLCF) 630 : NGCD = MELEFE.NUM(3,NLCF) 631 : SURF = MPOVSU.VPOCHA(NLCF,1) 632 : NLFL = MLELIM.LECT(NGCF) 633 : IF(NLFL .NE. 0)THEN 634 : C 635 : C************* The point belongs on BC -> No contribution to jacobian! 636 : C 637 : DFRO(1)=0.0D0 638 : DFRO(2)=0.0D0 639 : DFRO(3)=0.0D0 640 : DFRO(4)=0.0D0 641 : DFRUN(1)=0.0D0 642 : DFRUN(2)=0.0D0 643 : DFRUN(3)=0.0D0 644 : DFRUN(4)=0.0D0 645 : DFRUT(1)=0.0D0 646 : DFRUT(2)=0.0D0 647 : DFRUT(3)=0.0D0 648 : DFRUT(4)=0.0D0 649 : DFRET(1)=0.0D0 650 : DFRET(2)=0.0D0 651 : DFRET(3)=0.0D0 652 : DFRET(4)=0.0D0 653 : C 654 : C************* On est sur G, sur D, ou sur un mur? 655 : C 656 : ELSEIF(NGCG .EQ. NGCD)THEN 657 : C 658 : C***************** Murs 659 : C Aucune contribution à l'elt droite! 660 : C 661 : CNX = MELVNX.VELCHE(1,NLCF) 662 : CNY = MELVNY.VELCHE(1,NLCF) 663 : CTX = MELT1X.VELCHE(1,NLCF) 664 : CTY = MELT1Y.VELCHE(1,NLCF) 665 : ROG = MELRO.VELCHE(1,NLCF) 666 : UNG = MELVUN.VELCHE(1,NLCF) 667 : UTG = MELVUT.VELCHE(1,NLCF) 668 : UXG = UNG * CNX + UTG * CTX 669 : UYG = UNG * CNY + UTG * CTY 670 : PG = MELP.VELCHE(1,NLCF) 671 : GAMG = MELGAM.VELCHE(1,NLCF) 672 : C 673 : C**************** Jacobien par rapport aux variables primitives sur 674 : C la face 675 : C 676 : IF(IIMPL .EQ. 4)THEN 677 : C 678 : C******************* VLH 679 : C 679 : CALL VLHJ3(ROG,UXG,UYG,PG,GAMG,CNX,CNY,DFRUN) 680 : C 681 : ELSEIF(IIMPL .EQ. 5)THEN 682 : C 683 : C******************* AUSM+ 684 : C 685 : ROD=ROG 686 : PD=PG 687 : UND=-1.0D0*UNG 688 : UTD=UTG 689 : UXD = UND * CNX + UTD * CTX 690 : UYD = UND * CNY + UTD * CTY 690 : CALL AUSMW(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, 691 : & GAMG,CNX,CNY,DFRUN) 692 : ELSE 692 : CALL ERREUR(5) 693 : ENDIF 694 : DFRO(1)=0.0D0 695 : DFRO(2)=0.0D0 696 : DFRO(3)=0.0D0 697 : DFRO(4)=0.0D0 698 : DFRUT(1)=0.0D0 699 : DFRUT(2)=0.0D0 700 : DFRUT(3)=0.0D0 701 : DFRUT(4)=0.0D0 702 : DFRET(1)=0.0D0 703 : DFRET(2)=0.0D0 704 : DFRET(3)=0.0D0 705 : DFRET(4)=0.0D0 706 : C 707 : ELSEIF(NGC .EQ. NGCG)THEN 708 : C 709 : C************* On est sur G 710 : C 711 : CNX = MELVNX.VELCHE(1,NLCF) 712 : CNY = MELVNY.VELCHE(1,NLCF) 713 : CTX = MELT1X.VELCHE(1,NLCF) 714 : CTY = MELT1Y.VELCHE(1,NLCF) 715 : C 716 : ROG = MELRO.VELCHE(1,NLCF) 717 : UNG = MELVUN.VELCHE(1,NLCF) 718 : UTG = MELVUT.VELCHE(1,NLCF) 719 : UXG = UNG * CNX + UTG * CTX 720 : UYG = UNG * CNY + UTG * CTY 721 : PG = MELP.VELCHE(1,NLCF) 722 : GAMG = MELGAM.VELCHE(1,NLCF) 723 : ECIN=0.5D0*ROG*(UNG*UNG+UTG*UTG) 724 : RETG=PG/(GAMG-1.0D0)+ECIN 725 : C 726 : ROD = MELRO.VELCHE(3,NLCF) 727 : UND = MELVUN.VELCHE(3,NLCF) 728 : UTD = MELVUT.VELCHE(3,NLCF) 729 : UXD = UND * CNX + UTD * CTX 730 : UYD = UND * CNY + UTD * CTY 731 : PD = MELP.VELCHE(3,NLCF) 732 : C 733 : NLCD=MLENTC.LECT(NGCD) 734 : VOLD=MPVOLU.VPOCHA(NLCD,1) 735 : IF(IIMPL .EQ. 4)THEN 736 : C 737 : C******************* VLH 738 : C 738 : CALL VLHJ1(ROG,UXG,UYG,PG,RETG,GAMG,CNX,CNY,CTX,CTY 739 : $ ,DFRO,DFRUN,DFRUT,DFRET) 740 : ELSEIF(IIMPL .EQ. 5)THEN 741 : C 742 : C******************* AUSM+ 743 : C 743 : CALL AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, 744 : & GAMG,CNX,CNY,CTX,CTY, 745 : & DFRO,DFRUN,DFRUT,DFRET) 746 : ELSE 746 : CALL ERREUR(5) 747 : ENDIF 748 : ELSEIF(NGC .EQ. NGCD)THEN 749 : CNX = -1.0D0 * MELVNX.VELCHE(1,NLCF) 750 : CNY = -1.0D0 * MELVNY.VELCHE(1,NLCF) 751 : CTX = -1.0D0 * MELT1X.VELCHE(1,NLCF) 752 : CTY = -1.0D0 * MELT1Y.VELCHE(1,NLCF) 753 : C 754 : ROG = MELRO.VELCHE(3,NLCF) 755 : UNG = MELVUN.VELCHE(3,NLCF) 756 : UTG = MELVUT.VELCHE(3,NLCF) 757 : UXG = -1.0D0 * (UNG * CNX + UTG * CTX) 758 : UYG = -1.0D0 * (UNG * CNY + UTG * CTY) 759 : PG = MELP.VELCHE(3,NLCF) 760 : GAMG = MELGAM.VELCHE(3,NLCF) 761 : ECIN=0.5D0*ROG*(UNG*UNG+UTG*UTG) 762 : RETG=PG/(GAMG-1.0D0)+ECIN 763 : C 764 : ROD = MELRO.VELCHE(1,NLCF) 765 : UND = MELVUN.VELCHE(1,NLCF) 766 : UTD = MELVUT.VELCHE(1,NLCF) 767 : UXD = -1.0D0 * (UND * CNX + UTD * CTX) 768 : UYD = -1.0D0 * (UND * CNY + UTD * CTY) 769 : PD = MELP.VELCHE(1,NLCF) 770 : C 771 : NLCD=MLENTC.LECT(NGCG) 772 : VOLD=MPVOLU.VPOCHA(NLCD,1) 773 : IF(IIMPL .EQ. 4)THEN 774 : C 775 : C******************* VLH 776 : C 776 : CALL VLHJ1(ROG,UXG,UYG,PG,RETG,GAMG,CNX,CNY,CTX,CTY 777 : $ ,DFRO,DFRUN,DFRUT,DFRET) 778 : ELSEIF(IIMPL .EQ. 5)THEN 779 : C 780 : C******************* AUSM+ 781 : C 781 : CALL AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, 782 : & GAMG,CNX,CNY,CTX,CTY, 783 : & DFRO,DFRUN,DFRUT,DFRET) 784 : ELSE 784 : CALL ERREUR(5) 785 : ENDIF 786 : ELSE 787 : WRITE(IOIMP,*) 'FACEL = ???' 788 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 788 : CALL ERREUR(5) 789 : GOTO 9999 790 : ENDIF 791 : C 792 : DO I1=1,4 793 : DFRUX(I1)=DFRUN(I1)*CNX+DFRUT(I1)*CTX 794 : DFRUY(I1)=DFRUN(I1)*CNY+DFRUT(I1)*CTY 795 : ENDDO 796 : C 797 : NGNO=MELGRR.NUM(IFAC,IELPRI) 798 : NGN1=MELGRP.NUM(IFAC,IELPRI) 799 : NGN2=MELGRV.NUM(IFAC,IELPRI) 800 : C 801 : IF((NGN1 .NE. NGNO) .OR. (NGN2 .NE. NGNO))THEN 802 : WRITE(IOIMP,*) 'SPG COEFF GRAD =???' 803 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 803 : CALL ERREUR(5) 804 : GOTO 9999 805 : ENDIF 806 : C 807 : NLNO=MLENTC.LECT(NGNO) 808 : IF(NLNO .EQ. 0)THEN 809 : C 810 : C**************** MURS 811 : C 812 : IPT1.NUM(IFAC,IELPRI)=NGC 813 : 814 : C Primal variable in IPT1.NUM(IFAC,IELPRI)=IPT1.NUM(NP,IELPRI) 815 : C Dual variable in IPT1.NUM(NP,IELPRI) 816 : C 817 : FUNCEL = -1.0D0 * SURF / VOLG 818 : UXR.AM(IELPRI,IFAC,NP)=UXR.AM(IELPRI,IFAC,NP) 819 : & +DFRUX(1)*FUNCEL 820 : UYR.AM(IELPRI,IFAC,NP)=UYR.AM(IELPRI,IFAC,NP) 821 : & +DFRUY(1)*FUNCEL 822 : UXUX.AM(IELPRI,IFAC,NP)=UXUX.AM(IELPRI,IFAC,NP) 823 : & +DFRUX(2)*FUNCEL 824 : UYUX.AM(IELPRI,IFAC,NP)=UYUX.AM(IELPRI,IFAC,NP) 825 : & +DFRUY(2)*FUNCEL 826 : UXUY.AM(IELPRI,IFAC,NP)=UXUY.AM(IELPRI,IFAC,NP) 827 : & +DFRUX(3)*FUNCEL 828 : UYUY.AM(IELPRI,IFAC,NP)=UYUY.AM(IELPRI,IFAC,NP) 829 : & +DFRUY(3)*FUNCEL 830 : UXRET.AM(IELPRI,IFAC,NP)=UXRET.AM(IELPRI,IFAC,NP) 831 : & +DFRUX(4)*FUNCEL 832 : UYRET.AM(IELPRI,IFAC,NP)=UYRET.AM(IELPRI,IFAC,NP) 833 : & +DFRUY(4)*FUNCEL 834 : 835 : ELSE 836 : C 837 : C Primal variable in IPT1.NUM(NP,IELPRI) 838 : C Dual variable in IPT1.NUM(NP,IELPRI) 839 : C 840 : FUNCEL=DFRO(1)*SURF 841 : RR.AM(IELPRI,NP,NP)=RR.AM(IELPRI,NP,NP)- 842 : & FUNCEL/VOLG 843 : C 844 : C Primal variable in NP 845 : C Dual variable in IFAC 846 : C 847 : RR.AM(IELPRI,NP,IFAC)=RR.AM(IELPRI,NP,IFAC)+FUNCEL 848 : $ /VOLD 849 : C 850 : FUNCEL=DFRO(2)*SURF 851 : RUX.AM(IELPRI,NP,NP)=RUX.AM(IELPRI,NP,NP) - 852 : & FUNCEL/VOLG 853 : RUX.AM(IELPRI,NP,IFAC)=RUX.AM(IELPRI,NP,IFAC)+ FUNCEL 854 : $ /VOLD 855 : C 856 : FUNCEL=DFRO(3)*SURF 857 : RUY.AM(IELPRI,NP,NP)=RUY.AM(IELPRI,NP,NP) - 858 : & FUNCEL/VOLG 859 : RUY.AM(IELPRI,NP,IFAC)=RUY.AM(IELPRI,NP,IFAC)+ FUNCEL 860 : $ /VOLD 861 : C 862 : FUNCEL=DFRO(4)*SURF 863 : RRET.AM(IELPRI,NP,NP)=RRET.AM(IELPRI,NP,NP) - 864 : & FUNCEL/VOLG 865 : RRET.AM(IELPRI,NP,IFAC)= RRET.AM(IELPRI,NP,IFAC) 866 : $ +FUNCEL/VOLD 867 : C 868 : C 869 : FUNCEL=DFRUX(1)*SURF 870 : UXR.AM(IELPRI,NP,NP)=UXR.AM(IELPRI,NP,NP) - 871 : & FUNCEL/VOLG 872 : UXR.AM(IELPRI,NP,IFAC)=UXR.AM(IELPRI,NP,IFAC)+FUNCEL 873 : $ /VOLD 874 : C 875 : FUNCEL=DFRUX(2)*SURF 876 : UXUX.AM(IELPRI,NP,NP)=UXUX.AM(IELPRI,NP,NP) - 877 : & FUNCEL/VOLG 878 : UXUX.AM(IELPRI,NP,IFAC)=UXUX.AM(IELPRI,NP,IFAC)+ 879 : $ FUNCEL/VOLD 880 : C 881 : FUNCEL=DFRUX(3)*SURF 882 : UXUY.AM(IELPRI,NP,NP)=UXUY.AM(IELPRI,NP,NP) - 883 : & FUNCEL/VOLG 884 : UXUY.AM(IELPRI,NP,IFAC)=UXUY.AM(IELPRI,NP,IFAC)+ 885 : $ FUNCEL/VOLD 886 : C 887 : FUNCEL=DFRUX(4)*SURF 888 : UXRET.AM(IELPRI,NP,NP)=UXRET.AM(IELPRI,NP,NP) - 889 : & FUNCEL/VOLG 890 : UXRET.AM(IELPRI,NP,IFAC)=UXRET.AM(IELPRI,NP,IFAC)+ 891 : $ FUNCEL/VOLD 892 : C 893 : C 894 : FUNCEL=DFRUY(1)*SURF 895 : UYR.AM(IELPRI,NP,NP)=UYR.AM(IELPRI,NP,NP) - 896 : & FUNCEL/VOLG 897 : UYR.AM(IELPRI,NP,IFAC)=UYR.AM(IELPRI,NP,IFAC)+FUNCEL 898 : $ /VOLD 899 : C 900 : FUNCEL=DFRUY(2)*SURF 901 : UYUX.AM(IELPRI,NP,NP)=UYUX.AM(IELPRI,NP,NP) - 902 : & FUNCEL/VOLG 903 : UYUX.AM(IELPRI,NP,IFAC)=UYUX.AM(IELPRI,NP,IFAC)+ 904 : $ FUNCEL/VOLD 905 : C 906 : FUNCEL=DFRUY(3)*SURF 907 : UYUY.AM(IELPRI,NP,NP)=UYUY.AM(IELPRI,NP,NP) - 908 : & FUNCEL/VOLG 909 : UYUY.AM(IELPRI,NP,IFAC)=UYUY.AM(IELPRI,NP,IFAC)+ 910 : $ FUNCEL/VOLD 911 : C 912 : FUNCEL=DFRUY(4)*SURF 913 : UYRET.AM(IELPRI,NP,NP)=UYRET.AM(IELPRI,NP,NP) - 914 : & FUNCEL/VOLG 915 : UYRET.AM(IELPRI,NP,IFAC)=UYRET.AM(IELPRI,NP,IFAC)+ 916 : $ FUNCEL/VOLD 917 : C 918 : C 919 : FUNCEL=DFRET(1)*SURF 920 : RETR.AM(IELPRI,NP,NP)=RETR.AM(IELPRI,NP,NP) - 921 : & FUNCEL/VOLG 922 : RETR.AM(IELPRI,NP,IFAC)=RETR.AM(IELPRI,NP,IFAC)+FUNCEL 923 : $ /VOLD 924 : C 925 : FUNCEL=DFRET(2)*SURF 926 : RETUX.AM(IELPRI,NP,NP)=RETUX.AM(IELPRI,NP,NP) - 927 : & FUNCEL/VOLG 928 : RETUX.AM(IELPRI,NP,IFAC)=RETUX.AM(IELPRI,NP,IFAC)+ 929 : $ FUNCEL/VOLD 930 : C 931 : FUNCEL=DFRET(3)*SURF 932 : RETUY.AM(IELPRI,NP,NP)=RETUY.AM(IELPRI,NP,NP) - 933 : & FUNCEL/VOLG 934 : RETUY.AM(IELPRI,NP,IFAC)=RETUY.AM(IELPRI,NP,IFAC)+ 935 : $ FUNCEL/VOLD 936 : C 937 : FUNCEL=DFRET(4)*SURF 938 : RETRET.AM(IELPRI,NP,NP)=RETRET.AM(IELPRI,NP,NP) - 939 : & FUNCEL/VOLG 940 : RETRET.AM(IELPRI,NP,IFAC)=RETRET.AM(IELPRI,NP,IFAC)+ 941 : $ FUNCEL/VOLD 942 : ENDIF 943 : C 944 : C************* Gradients computed with EULESCAL or EULEVECT 945 : C Computation of gradients contributions 946 : 947 : DO IFAC2 = 1, (NP-1), 1 948 : C 949 : C**************** On controlle que les melemes de ELTFA et 950 : C le SPG du grad ont le meme ordre 951 : C 952 : NGFGR=MELEME.NUM(IFAC2,IELPRI) 953 : NLFGR=MLENTF.LECT(NGFGR) 954 : NG=MELEFE.NUM(1,NLFGR) 955 : ND=MELEFE.NUM(3,NLFGR) 956 : C 957 : IF(NG .EQ. ND)THEN 958 : C 959 : C******************* On est sur un mur: EULESCAL 960 : C 961 : DO IFACGR = 1, (NP -1) , 1 962 : C 963 : C********************** Boucle sur les elts de MELGRR pour trouver la position 964 : C correspondante 965 : C 966 : NGGRA=MELGRR.NUM(IFACGR,IELPRI) 967 : IF(NGFGR .EQ. NGGRA) GOTO 10 968 : ENDDO 969 : WRITE(IOIMP,*) 'ELTFA et SPG grad= ???' 970 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 970 : CALL ERREUR(5) 971 : GOTO 9999 972 : 10 CONTINUE 973 : C 974 : C******************* Les coeff qui nous interesse sont donc dans la 975 : C position IFACGR du gradient 976 : C 977 : CELLR = ALR * (DXCF * MELCR1.VELCHE(IFACGR,IELPRI) 978 : $ +DYCF * MELCR2.VELCHE(IFACGR,IELPRI)) 979 : CELLP = ALP * (DXCF * MELCP1.VELCHE(IFACGR,IELPRI) 980 : $ + DYCF * MELCP2.VELCHE(IFACGR,IELPRI)) 981 : C 982 : C******************* La vitesse est un peux plus difficile 983 : C (HP EULEVECT ou vitesse imposé) 984 : C 985 : C CViVj contribution de vj sur le grad de vi 986 : C 987 : NLVIMP = MLEVL.LECT(NGGRA) 988 : IF(NLVIMP .NE.0)THEN 989 : C 990 : C*********************** NGGRA spg de conditions limites -> no 991 : C contribution 992 : C 993 : CVXVX=0.0D0 994 : CVXVY=0.0D0 995 : CVYVX=0.0D0 996 : CVYVY=0.0D0 997 : ELSE 998 : CNX1=MPNORM.VPOCHA(NLFGR,1) 999 : CNY1=MPNORM.VPOCHA(NLFGR,2) 1000 : CTX1=-1.0D0*CNY1 1001 : CTY1=CNX1 1002 : FUNCEL=(DXCF * MELCV1.VELCHE(IFACGR,IELPRI)+DYCF 1003 : $ * MELCV2.VELCHE(IFACGR,IELPRI)) 1004 : CVXVX=FUNCEL*(CTX1*CTX1-CNX1*CNX1)*ALVX 1005 : CVXVY=FUNCEL*(CTX1*CTY1-CNX1*CNY1)*ALVX 1006 : CVYVX=FUNCEL*(CTX1*CTY1-CNX1*CNY1)*ALVY 1007 : CVYVY=FUNCEL*(CTY1*CTY1-CNY1*CNY1)*ALVY 1008 : ENDIF 1009 : C 1010 : ELSE 1011 : DO IFACGR = 1, (NP -1) , 1 1012 : C 1013 : C********************** Boucle sur les elts de MELGRR pour trouver la position 1014 : C correspondante 1015 : C 1016 : NGGRA=MELGRR.NUM(IFACGR,IELPRI) 1017 : IF((NGGRA .EQ. NG) .OR. (NGGRA .EQ. ND)) GOTO 20 1018 : ENDDO 1019 : WRITE(IOIMP,*) 'ELTFA et SPG grad= ???' 1020 : WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' 1020 : CALL ERREUR(5) 1021 : GOTO 9999 1022 : 20 CONTINUE 1023 : CELLR = ALR * (DXCF * MELCR1.VELCHE(IFACGR,IELPRI) 1024 : $ +DYCF * MELCR2.VELCHE(IFACGR,IELPRI)) 1025 : CELLP = ALP * (DXCF * MELCP1.VELCHE(IFACGR,IELPRI) 1026 : $ + DYCF * MELCP2.VELCHE(IFACGR,IELPRI)) 1027 : CVXVX = ALVX * (DXCF * MELCV1.VELCHE(IFACGR,IELPRI) 1028 : $ + DYCF * MELCV2.VELCHE(IFACGR,IELPRI)) 1029 : CVXVY=0.0D0 1030 : CVYVY = ALVY * (DXCF * MELCV1.VELCHE(IFACGR,IELPRI) 1031 : $ + DYCF * MELCV2.VELCHE(IFACGR,IELPRI)) 1032 : CVYVX = 0.0D0 1033 : 1034 : ENDIF 1035 : C 1036 : C**************** Dans le cas d'un mur (NLN0=0) , IFACGR = primale qui donne 1037 : C contribution en IELPRI,IFACGR,IFAC OU (IELPRI,IFACGR,NP) 1038 : C 1039 : C Dans le cas non mur (NLN0=0) , IFACGR = primale qui donne 1040 : C contribution en IELPRI,IFACGR,NP et en IELPRI,IFAC,NP 1041 : C 1042 : IF(NLNO .EQ. 0)THEN 1043 : C 1044 : C**************** MURS 1045 : C 1046 : FUNCEL = -1.0D0 * SURF / VOLG * CELLR 1047 : UXR.AM(IELPRI,IFACGR,NP)=UXR.AM(IELPRI,IFACGR,NP) 1048 : & +DFRUX(1)*FUNCEL 1049 : UYR.AM(IELPRI,IFACGR,NP)=UYR.AM(IELPRI,IFACGR,NP) 1050 : & +DFRUY(1)*FUNCEL 1051 : C 1052 : FUNCEL = -1.0D0 * SURF / VOLG 1053 : UXUX.AM(IELPRI,IFACGR,NP)=UXUX.AM(IELPRI,IFACGR,NP) 1054 : & +FUNCEL * (DFRUX(2)*CVXVX + DFRUX(3) *CVYVX) 1055 : UYUX.AM(IELPRI,IFACGR,NP)=UYUX.AM(IELPRI,IFACGR,NP) 1056 : & +FUNCEL * (DFRUY(2)*CVXVX + DFRUY(3) *CVYVX) 1057 : C 1058 : FUNCEL = -1.0D0 * SURF / VOLG 1059 : UXUY.AM(IELPRI,IFACGR,NP)=UXUY.AM(IELPRI,IFACGR,NP) 1060 : & +FUNCEL * (DFRUX(3)*CVYVY + DFRUX(2) *CVXVY) 1061 : UYUY.AM(IELPRI,IFACGR,NP)=UYUY.AM(IELPRI,IFACGR,NP) 1062 : & +FUNCEL * (DFRUY(3)*CVYVY + DFRUY(2) *CVXVY) 1063 : 1064 : C 1065 : FUNCEL = -1.0D0 * SURF / VOLG * CELLP 1066 : UXRET.AM(IELPRI,IFACGR,NP)=UXRET.AM(IELPRI,IFACGR 1067 : $ ,NP)+DFRUX(4)*FUNCEL 1068 : UYRET.AM(IELPRI,IFACGR,NP)=UYRET.AM(IELPRI,IFACGR 1069 : $ ,NP)+DFRUY(4)*FUNCEL 1070 : ELSE 1071 : C 1072 : C******************* dR/ 1073 : C 1074 : FUNCEL=DFRO(1)*SURF*CELLR 1075 : RR.AM(IELPRI,IFACGR,NP)=RR.AM(IELPRI,IFACGR,NP)- 1076 : & FUNCEL/VOLG 1077 : RR.AM(IELPRI,IFACGR,IFAC)=RR.AM(IELPRI,IFACGR,IFAC) 1078 : $ +FUNCEL/VOLD 1079 : C 1080 : FUNCEL=DFRO(4)*SURF*CELLP 1081 : RRET.AM(IELPRI,IFACGR,NP)=RRET.AM(IELPRI,IFACGR,NP) 1082 : $ -FUNCEL/VOLG 1083 : RRET.AM(IELPRI,IFACGR,IFAC)=RRET.AM(IELPRI,IFACGR 1084 : $ ,IFAC)+FUNCEL/VOLD 1085 : C 1086 : FUNCEL=SURF * (DFRO(2)*CVXVX + DFRO(3) *CVYVX) 1087 : RUX.AM(IELPRI,IFACGR,NP)=RUX.AM(IELPRI,IFACGR,NP) 1088 : & - FUNCEL / VOLG 1089 : RUX.AM(IELPRI,IFACGR,IFAC)=RUX.AM(IELPRI,IFACGR 1090 : $ ,IFAC) + FUNCEL / VOLD 1091 : C 1092 : FUNCEL=SURF * (DFRO(3)*CVYVY + DFRO(2) *CVXVY) 1093 : RUY.AM(IELPRI,IFACGR,NP)=RUY.AM(IELPRI,IFACGR,NP) 1094 : & - FUNCEL / VOLG 1095 : RUY.AM(IELPRI,IFACGR,IFAC)=RUY.AM(IELPRI,IFACGR 1096 : $ ,IFAC) + FUNCEL / VOLD 1097 : C 1098 : C******************* dUX/ 1099 : C 1100 : FUNCEL=DFRUX(1)*SURF*CELLR 1101 : UXR.AM(IELPRI,IFACGR,NP)=UXR.AM(IELPRI,IFACGR,NP) - 1102 : & FUNCEL/VOLG 1103 : UXR.AM(IELPRI,IFACGR,IFAC)= UXR.AM(IELPRI,IFACGR 1104 : $ ,IFAC)+FUNCEL/VOLD 1105 : C 1106 : FUNCEL=DFRUX(4)*SURF*CELLP 1107 : UXRET.AM(IELPRI,IFACGR,NP)=UXRET.AM(IELPRI,IFACGR 1108 : $ ,NP) -FUNCEL/VOLG 1109 : UXRET.AM(IELPRI,IFACGR,IFAC)= UXRET.AM(IELPRI 1110 : $ ,IFACGR,IFAC)+FUNCEL/VOLD 1111 : C 1112 : FUNCEL=SURF * (DFRUX(2)*CVXVX + DFRUX(3) *CVYVX) 1113 : UXUX.AM(IELPRI,IFACGR,NP)=UXUX.AM(IELPRI,IFACGR,NP) 1114 : & - FUNCEL / VOLG 1115 : UXUX.AM(IELPRI,IFACGR,IFAC)=UXUX.AM(IELPRI,IFACGR 1116 : $ ,IFAC) + FUNCEL / VOLD 1117 : C 1118 : FUNCEL=SURF * (DFRUX(3)*CVYVY + DFRUX(2) *CVXVY) 1119 : UXUY.AM(IELPRI,IFACGR,NP)=UXUY.AM(IELPRI,IFACGR,NP) 1120 : & - FUNCEL / VOLG 1121 : UXUY.AM(IELPRI,IFACGR,IFAC)=UXUY.AM(IELPRI,IFACGR 1122 : $ ,IFAC) + FUNCEL / VOLD 1123 : C 1124 : C******************* dUY/ 1125 : C 1126 : FUNCEL=DFRUY(1)*SURF*CELLR 1127 : UYR.AM(IELPRI,IFACGR,NP)=UYR.AM(IELPRI,IFACGR,NP) - 1128 : & FUNCEL/VOLG 1129 : UYR.AM(IELPRI,IFACGR,IFAC)=UYR.AM(IELPRI,IFACGR 1130 : $ ,IFAC) + FUNCEL/VOLD 1131 : C 1132 : FUNCEL=DFRUY(4)*SURF*CELLP 1133 : UYRET.AM(IELPRI,IFACGR,NP)=UYRET.AM(IELPRI,IFACGR 1134 : $ ,NP) -FUNCEL/VOLG 1135 : UYRET.AM(IELPRI,IFACGR,IFAC)=UYRET.AM(IELPRI,IFACGR 1136 : $ ,IFAC) + FUNCEL/VOLD 1137 : C 1138 : FUNCEL=SURF * (DFRUY(2)*CVXVX + DFRUY(3) *CVYVX) 1139 : UYUX.AM(IELPRI,IFACGR,NP)=UYUX.AM(IELPRI,IFACGR,NP) 1140 : & - FUNCEL / VOLG 1141 : UYUX.AM(IELPRI,IFACGR,IFAC)=UYUX.AM(IELPRI,IFACGR 1142 : $ ,IFAC) + FUNCEL / VOLD 1143 : C 1144 : FUNCEL=SURF * (DFRUY(3)*CVYVY + DFRUY(2) *CVXVY) 1145 : UYUY.AM(IELPRI,IFACGR,NP)=UYUY.AM(IELPRI,IFACGR,NP) 1146 : & - FUNCEL / VOLG 1147 : UYUY.AM(IELPRI,IFACGR,IFAC)=UYUY.AM(IELPRI,IFACGR 1148 : $ ,IFAC) + FUNCEL / VOLD 1149 : C 1150 : C******************* dRET/ 1151 : C 1152 : FUNCEL=DFRET(1)*SURF*CELLR 1153 : RETR.AM(IELPRI,IFACGR,NP)=RETR.AM(IELPRI,IFACGR,NP) 1154 : $ -FUNCEL/VOLG 1155 : RETR.AM(IELPRI,IFACGR,IFAC)=RETR.AM(IELPRI,IFACGR 1156 : $ ,IFAC)+FUNCEL/VOLD 1157 : C 1158 : FUNCEL=DFRET(4)*SURF*CELLP 1159 : RETRET.AM(IELPRI,IFACGR,NP)=RETRET.AM(IELPRI,IFACGR 1160 : $ ,NP)-FUNCEL/VOLG 1161 : RETRET.AM(IELPRI,IFACGR,IFAC)=RETRET.AM(IELPRI 1162 : $ ,IFACGR,IFAC)+FUNCEL/VOLD 1163 : C 1164 : FUNCEL=SURF * (DFRET(2)*CVXVX + DFRET(3) *CVYVX) 1165 : RETUX.AM(IELPRI,IFACGR,NP)=RETUX.AM(IELPRI,IFACGR 1166 : $ ,NP)- FUNCEL / VOLG 1167 : RETUX.AM(IELPRI,IFACGR,IFAC)=RETUX.AM(IELPRI,IFACGR 1168 : $ ,IFAC) + FUNCEL / VOLD 1169 : C 1170 : FUNCEL=SURF * (DFRET(3)*CVYVY + DFRET(2) *CVXVY) 1171 : RETUY.AM(IELPRI,IFACGR,NP)=RETUY.AM(IELPRI,IFACGR 1172 : $ ,NP)- FUNCEL / VOLG 1173 : RETUY.AM(IELPRI,IFACGR,IFAC)=RETUY.AM(IELPRI,IFACGR 1174 : $ ,IFAC) + FUNCEL / VOLD 1175 : C 1176 : ENDIF 1177 : ENDDO 1178 : ENDDO 1179 : ENDDO 1180 : C 1181 : C 1182 : C******* Transformation: 1183 : C jacobien par rapport aux variables primitives -> 1184 : C jacobien par rapport aux variables conservatives 1185 : C Optimisation du programme: 1186 : C on pourrais faire ça avant 1187 : C 1188 : NBEL=RR.AM(/1) 1189 : NP=RR.AM(/2) 1190 : MP=RR.AM(/3) 1191 : C 1192 : C******* Boucle sur les variables primales 1193 : C 1194 : DO IELETR=1,NBEL,1 1195 : DO IPRITR=1,NP,1 1196 : NGPRI=IPT1.NUM(IPRITR,IELETR) 1197 : NLPRI=MLENTC.LECT(NGPRI) 1198 : ROG=MPOVRN.VPOCHA(NLPRI,1) 1199 : PG=MPOVPN.VPOCHA(NLPRI,1) 1200 : GAMG=MPOGAM.VPOCHA(NLPRI,1) 1201 : UXG=MPOVVN.VPOCHA(NLPRI,1) 1202 : UYG=MPOVVN.VPOCHA(NLPRI,2) 1203 : GAMGM1=GAMG-1.0D0 1204 : CELCO1=UXG/ROG 1205 : CELCO2=UYG/ROG 1206 : CELCO3=0.5D0*GAMGM1*(UXG*UXG+UYG*UYG) 1207 : CELCO4=GAMGM1*UXG 1208 : CELCO5=GAMGM1*UYG 1209 : DO IDUATR=1,MP,1 1210 : DDRO=RR.AM(IELETR,IPRITR,IDUATR) 1211 : DDUX=RUX.AM(IELETR,IPRITR,IDUATR) 1212 : DDUY=RUY.AM(IELETR,IPRITR,IDUATR) 1213 : DDP=RRET.AM(IELETR,IPRITR,IDUATR) 1214 : RR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) 1215 : & -(CELCO2*DDUY)+(CELCO3*DDP) 1216 : RUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - 1217 : & CELCO4*DDP 1218 : RUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - 1219 : & CELCO5*DDP 1220 : RRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP 1221 : C 1222 : DDRO=UXR.AM(IELETR,IPRITR,IDUATR) 1223 : DDUX=UXUX.AM(IELETR,IPRITR,IDUATR) 1224 : DDUY=UXUY.AM(IELETR,IPRITR,IDUATR) 1225 : DDP=UXRET.AM(IELETR,IPRITR,IDUATR) 1226 : UXR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) 1227 : & -(CELCO2*DDUY)+(CELCO3*DDP) 1228 : UXUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - 1229 : & CELCO4*DDP 1230 : UXUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - 1231 : & CELCO5*DDP 1232 : UXRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP 1233 : C 1234 : DDRO=UYR.AM(IELETR,IPRITR,IDUATR) 1235 : DDUX=UYUX.AM(IELETR,IPRITR,IDUATR) 1236 : DDUY=UYUY.AM(IELETR,IPRITR,IDUATR) 1237 : DDP=UYRET.AM(IELETR,IPRITR,IDUATR) 1238 : UYR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) 1239 : & -(CELCO2*DDUY)+(CELCO3*DDP) 1240 : UYUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - 1241 : & CELCO4*DDP 1242 : UYUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - 1243 : & CELCO5*DDP 1244 : UYRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP 1245 : C 1246 : DDRO=RETR.AM(IELETR,IPRITR,IDUATR) 1247 : DDUX=RETUX.AM(IELETR,IPRITR,IDUATR) 1248 : DDUY=RETUY.AM(IELETR,IPRITR,IDUATR) 1249 : DDP=RETRET.AM(IELETR,IPRITR,IDUATR) 1250 : RETR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) 1251 : & -(CELCO2*DDUY)+(CELCO3*DDP) 1252 : RETUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - 1253 : & CELCO4*DDP 1254 : RETUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - 1255 : & CELCO5*DDP 1256 : RETRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP 1257 : ENDDO 1258 : ENDDO 1259 : ENDDO 1260 : C 1261 : MELEME = MLESPG.LECT(ISOUS) 1262 : SEGDES MELEME 1263 : MELEME = MLEPRI.LECT(ISOUS) 1264 : SEGDES MELEME 1265 : MELEME = MLEGRR.LECT(ISOUS) 1266 : SEGDES MELEME 1267 : MELEME = MLEGRP.LECT(ISOUS) 1268 : SEGDES MELEME 1269 : MELEME = MLEGRV.LECT(ISOUS) 1270 : SEGDES MELEME 1271 : C 1272 : SEGDES RR , RUX , RUY , RRET , 1273 : & UXR , UXUX , UXUY , UXRET , 1274 : & UYR , UYUX , UYUY , UYRET , 1275 : & RETR , RETUX , RETUY , RETRET 1276 : SEGDES MELCR1, MELCR2 1277 : SEGDES MELCP1, MELCP2 1278 : SEGDES MELCV1, MELCV2 1279 : ENDDO 1280 : C 1281 : C 1282 : C**** Les MELVALs 1283 : C 1284 : SEGDES MELRO 1285 : SEGDES MELP 1286 : SEGDES MELGAM 1287 : SEGDES MELVUN 1288 : SEGDES MELVUT 1289 : SEGDES MELVNX 1290 : SEGDES MELVNY 1291 : SEGDES MELT1X 1292 : SEGDES MELT1Y 1293 : C 1294 : C**** Les MPOVALs 1295 : C 1296 : C 1297 : SEGDES MPOVRN 1298 : SEGDES MPOVVN 1299 : SEGDES MPOVPN 1300 : SEGDES MPOGAM 1301 : SEGDES MPVGRN 1302 : SEGDES MPGLRN 1303 : SEGDES MPVGPN 1304 : SEGDES MPGLPN 1305 : SEGDES MPVGPN 1306 : SEGDES MPGLPN 1307 : IF(MPVLIM .NE. 0) SEGDES MPVLIM 1308 : C 1309 : C**** Les MELEMEs, les MLENTIs, les volumes, les surfaces 1310 : C 1311 : SEGDES MELEMC 1312 : SEGSUP MLENTC 1313 : SEGSUP MLENTF 1314 : SEGSUP MLEVL 1315 : SEGDES MELEFE 1316 : SEGDES MPVOLU 1317 : SEGDES MPNORM 1318 : SEGDES MPOVSU 1319 : SEGDES MELPRI 1320 : SEGSUP MLEPRI 1321 : SEGSUP MLEGRR 1322 : SEGSUP MLEGRP 1323 : SEGSUP MLEGRV 1324 : 1325 : C 1326 : C**** Coeff GRAD 1327 : C 1328 : MCHEL1=ICRN 1329 : SEGDES MCHEL1 1330 : MCHEL1=ICPN 1331 : SEGDES MCHEL1 1332 : MCHEL1=ICVN 1333 : SEGDES MCHEL1 1334 : C 1335 : SEGSUP MLESPG 1336 : SEGDES MELTFA 1337 : C 1338 : C**** MATRIK, liste des inconnues 1339 : C 1340 : SEGDES MATRIK 1341 : SEGDES IMATRI 1342 : SEGDES MLMINC 1343 : 1344 : C 1345 : C 1346 : SEGDES RR , RUX , RUY , RRET , 1347 : & UXR , UXUX , UXUY , UXRET , 1348 : & UYR , UYUX , UYUY , UYRET , 1349 : & RETR , RETUX , RETUY , RETRET 1350 : 1351 : SEGSUP MLENTC 1352 : SEGSUP MLENTF 1353 : SEGDES MLMINC 1354 : SEGDES MLELIM 1355 : 1356 : 9999 CONTINUE 1357 : RETURN 1358 : END 1359 : 1360 : 1361 : 1362 : 1363 : 1364 : 1365 : 1366 :
© Cast3M 2003 - All rights reserved.
Disclaimer