1 : C KONJR1 SOURCE MAGN 09/07/30 21:18:23 6458 2 : SUBROUTINE KONJR1(ILINC,IRN,IUN,IPN,IGAMN,INORM,ICHPVO 3 : $ ,ICHPSU,IUINF,IUPRI,MELEMC,MELEFE,MELLIM,IMAT) 4 : C 5 : C************************************************************************ 6 : C 7 : C PROJET : CASTEM 2000 8 : C 9 : C NOM : KONJR1 10 : C 11 : C DESCRIPTION : Voir KON12 12 : C Calcul du jacobien du résidu pour la méthode 13 : C RUSANOLM (CENTERED + diff. num) 14 : C 14 : C Cas deux dimensions, gaz "calorically perfect" 15 : C 16 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 17 : C 18 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 19 : C 20 : C************************************************************************ 21 : C 22 : C 23 : C APPELES (Outils 24 : C CASTEM) : KRIPAD, LICHT, ERREUR 25 : C 26 : C APPELES (Calcul) : CENTJ0, CENTJ2, FRBM1 27 : C 28 : C************************************************************************ 29 : C 30 : C ENTREES 31 : C 32 : C ILINC : liste des inconnues (pointeur d'un objet de type LISTMOTS) 33 : C 34 : C 1) Pointeurs des CHPOINT 35 : C 36 : C IRN : CHPOINT CENTRE contenant la masse volumique ; 37 : C 38 : C IUN : CHPOINT CENTRE contenant la vitesse ; 39 : C 40 : C IPN : CHPOINT CENTRE contenant la pression ; 41 : C 42 : C IGAMN : CHPOINT CENTRE contenant le gamma ; 43 : C 44 : C INORM : CHPOINT FACE contenant les normales aux faces ; 45 : C 46 : C ICHPVO : CHPOINT VOLUME contenant le volume 47 : C 48 : C ICHPSU : CHPOINT FACE contenant la surface des faces 49 : C 50 : C 51 : C 2) Pointeurs de MELEME de la table DOMAINE 52 : C 53 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES 54 : C 55 : C MELEFE : MELEME 'FACEL' du connectivité Faces -> Elts 56 : C 57 : C MELLIM : MELEME sur lequel on ne veut pas calculer 58 : C la contribution à la matrice jacobienne 59 : C 60 : C 3) Cas Bas Mach 61 : C 62 : C IUINF : CHPOINT, one component, "cut-off velocity" 63 : C 0 if there is no Bas MACH 64 : C 65 : C IUPRI : CHPOINT, one component, second "cut-off velocity" 66 : C 0 if there is no Bas MACH 67 : C 68 : C 69 : C SORTIES 70 : C 71 : C IMAT : pointeur de la MATRIK du jacobien du residu 72 : C 73 : C************************************************************************ 74 : C 75 : C HISTORIQUE (Anomalies et modifications éventuelles) 76 : C 77 : C HISTORIQUE : 78 : C 79 : C************************************************************************ 80 : C 81 : C 82 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0 83 : C GAMMA \in (1,3) 84 : C Si non il faut le faire!!! 85 : C 86 : C************************************************************************ 87 : C 88 : C 89 : C**** Variables de COOPTIO 90 : C 91 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 92 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 93 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 94 : C & ,IECHO, IIMPI, IOSPI 95 : C & ,IDIM 96 : C & ,MCOORD 97 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 98 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 99 : C & ,NORINC,NORVAL,NORIND,NORVAD 100 : C & ,NUCROU, IPSAUV 101 : C & ,IFICLE,IPREFI,IREFOR,ISAFOR 102 : C 103 : IMPLICIT INTEGER(I-N) 104 : INTEGER ILINC, IRN,IUN,IPN,IGAMN,INORM,ICHPVO,ICHPSU, IUINF, IUPRI 105 : & , IMAT, IGEOMC, IGEOMF 106 : & , NFAC, NBSOUS, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID 107 : & , NKMT, NBME, NBEL, MP, NP 108 : & , IFAC, NGCF, NLCF, NGCG, NGCD, NLCG, NLCD, NLFL, I1 109 : REAL*8 ROG, PG, UXG, UYG, GAMG, VOLG 110 : & , ROD, PD, UXD, UYD, GAMD, VOLD 111 : & , SURF, CNX, CNY, CTX, CTY, FUNCEL 112 : & , DFRO(4), DFRET(4), DFRUN(4), DFRUT(4), VINF 113 : & , UNG, UTG, UND, UTD, ROF, PF, UNF, UTF, GAMF, AMAT(4,4) 114 : & , RLAMMA, BMAT(4,4) 115 : CHARACTER*8 TYPE 116 : C 117 : C**** LES INCLUDES 118 : C 119 : -INC CCOPTIO 120 : -INC SMCHPOI 121 : -INC SMMATRIK 122 : -INC SMELEME 123 : -INC SMLMOTS 124 : -INC SMLENTI 125 : POINTEUR MPRN.MPOVAL, MPUN.MPOVAL, MPPN.MPOVAL, MPGAMN.MPOVAL, 126 : & MPNORM.MPOVAL, MPVOLU.MPOVAL, MPOVSU.MPOVAL, 127 : & MPUPRI.MPOVAL, MPUINF.MPOVAL 128 : POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME, 129 : & MELEDU.MELEME, MELLIM.MELEME 130 : POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLELIM.MLENTI 131 : POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RRET.IZAFM, 132 : & UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXRET.IZAFM, 133 : & UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYRET.IZAFM, 134 : & RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETRET.IZAFM 135 : POINTEUR MLMINC.MLMOTS 136 : C 137 : C**** KRIPAD pour la correspondance global/local des conditions limits 138 : C 138 : CALL KRIPAD(MELLIM,MLELIM) 139 : c SEGACT MELLIM 140 : C 141 : C**** KRIPAD pour la correspondance global/local des centres 142 : C 142 : CALL KRIPAD(MELEMC,MLENTC) 143 : C 144 : C SEGACT MLENTC 145 : SEGACT MELEMC 146 : C 147 : SEGACT MELEFE 148 : C 148 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF) 148 : CALL LICHT(INORM,MPNORM,TYPE,IGEOMF) 148 : CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC) 149 : C 150 : C**** LICHT active les MPOVALs en *MOD 151 : C 152 : C i.e. 153 : C 154 : C SEGACT MPOVSU*MOD 155 : C SEGACT MPOVNO*MOD 156 : C SEGACT MPVOLU*MOD 157 : C 158 : MELEMF = IGEOMF 158 : CALL KRIPAD(MELEMF,MLENTF) 159 : C 160 : C SEGACT MLENTF 161 : SEGACT MELEMF 162 : C 162 : CALL LICHT(IRN,MPRN,TYPE,IGEOMC) 162 : CALL LICHT(IPN,MPPN,TYPE,IGEOMC) 162 : CALL LICHT(IUN,MPUN,TYPE,IGEOMC) 162 : CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC) 163 : C 164 : C SEGACT MPRN*MOD 165 : C SEGACT MPPN*MOD 166 : C SEGACT MPUN*MOD 167 : C SEGACT MPGAMN*MOD 168 : C 169 : C************************************************************** 170 : C Bas Mach 171 : C************************************************************** 171 : CALL LICHT(IUPRI,MPUPRI,TYPE,IGEOMC) 171 : CALL LICHT(IUINF,MPUINF,TYPE,IGEOMC) 172 : C 173 : NFAC = MELEFE.NUM(/2) 174 : C 175 : C**** Maillage des inconnues primales 176 : C 177 : NBSOUS = 0 178 : NBREF = 0 179 : NBELEM = NFAC 180 : NBNN = 2 181 : C 182 : SEGINI MELEDU 183 : C MELEPR = MELEDU 184 : C 185 : C**** MELEDU = 'SEG2' 186 : C 187 : MELEDU.ITYPEL = 2 188 : C 189 : NRIGE = 7 190 : NMATRI = 1 191 : NKID = 9 192 : NKMT = 7 193 : C 194 : SEGINI MATRIK 195 : IMAT = MATRIK 196 : MATRIK.IRIGEL(1,1) = MELEDU 197 : MATRIK.IRIGEL(2,1) = MELEDU 198 : C 199 : C**** Matrice non symetrique 200 : C 201 : MATRIK.IRIGEL(7,1) = 2 202 : C 203 : NBME = 16 204 : NBSOUS = 1 205 : SEGINI IMATRI 206 : MLMINC = ILINC 207 : SEGACT MLMINC 208 : MATRIK.IRIGEL(4,1) = IMATRI 209 : C 210 : IMATRI.LISPRI(1) = MLMINC.MOTS(1) 211 : IMATRI.LISPRI(2) = MLMINC.MOTS(2) 212 : IMATRI.LISPRI(3) = MLMINC.MOTS(3) 213 : IMATRI.LISPRI(4) = MLMINC.MOTS(4) 214 : IMATRI.LISPRI(5) = MLMINC.MOTS(1) 215 : IMATRI.LISPRI(6) = MLMINC.MOTS(2) 216 : IMATRI.LISPRI(7) = MLMINC.MOTS(3) 217 : IMATRI.LISPRI(8) = MLMINC.MOTS(4) 218 : IMATRI.LISPRI(9) = MLMINC.MOTS(1) 219 : IMATRI.LISPRI(10) = MLMINC.MOTS(2) 220 : IMATRI.LISPRI(11) = MLMINC.MOTS(3) 221 : IMATRI.LISPRI(12) = MLMINC.MOTS(4) 222 : IMATRI.LISPRI(13) = MLMINC.MOTS(1) 223 : IMATRI.LISPRI(14) = MLMINC.MOTS(2) 224 : IMATRI.LISPRI(15) = MLMINC.MOTS(3) 225 : IMATRI.LISPRI(16) = MLMINC.MOTS(4) 226 : C 227 : IMATRI.LISDUA(1) = MLMINC.MOTS(1) 228 : IMATRI.LISDUA(2) = MLMINC.MOTS(1) 229 : IMATRI.LISDUA(3) = MLMINC.MOTS(1) 230 : IMATRI.LISDUA(4) = MLMINC.MOTS(1) 231 : IMATRI.LISDUA(5) = MLMINC.MOTS(2) 232 : IMATRI.LISDUA(6) = MLMINC.MOTS(2) 233 : IMATRI.LISDUA(7) = MLMINC.MOTS(2) 234 : IMATRI.LISDUA(8) = MLMINC.MOTS(2) 235 : IMATRI.LISDUA(9) = MLMINC.MOTS(3) 236 : IMATRI.LISDUA(10) = MLMINC.MOTS(3) 237 : IMATRI.LISDUA(11) = MLMINC.MOTS(3) 238 : IMATRI.LISDUA(12) = MLMINC.MOTS(3) 239 : IMATRI.LISDUA(13) = MLMINC.MOTS(4) 240 : IMATRI.LISDUA(14) = MLMINC.MOTS(4) 241 : IMATRI.LISDUA(15) = MLMINC.MOTS(4) 242 : IMATRI.LISDUA(16) = MLMINC.MOTS(4) 243 : C 244 : NBEL = NBELEM 245 : NBSOUS = 1 246 : NP = 2 247 : MP = 2 248 : SEGINI RR , RUX , RUY , RRET , 249 : & UXR , UXUX , UXUY , UXRET , 250 : & UYR , UYUX , UYUY , UYRET , 251 : & RETR , RETUX , RETUY , RETRET 252 : C 253 : C**** Duale = IMATRI.LISDUA(1) = 'RN' 254 : C Primale = IMATRI.LISPRI(1) = 'RN' 255 : C -> IMATRI.LIZAFM(1,1) = RR 256 : C 257 : C Duale = IMATRI.LISDUA(2) = 'RN' 258 : C Primale = IMATRI.LISPRI(1) = 'RUXN' 259 : C -> IMATRI.LIZAFM(1,2) = RUX 260 : C ... 261 : C 262 : IMATRI.LIZAFM(1,1) = RR 263 : IMATRI.LIZAFM(1,2) = RUX 264 : IMATRI.LIZAFM(1,3) = RUY 265 : IMATRI.LIZAFM(1,4) = RRET 266 : IMATRI.LIZAFM(1,5) = UXR 267 : IMATRI.LIZAFM(1,6) = UXUX 268 : IMATRI.LIZAFM(1,7) = UXUY 269 : IMATRI.LIZAFM(1,8) = UXRET 270 : IMATRI.LIZAFM(1,9) = UYR 271 : IMATRI.LIZAFM(1,10) = UYUX 272 : IMATRI.LIZAFM(1,11) = UYUY 273 : IMATRI.LIZAFM(1,12) = UYRET 274 : IMATRI.LIZAFM(1,13) = RETR 275 : IMATRI.LIZAFM(1,14) = RETUX 276 : IMATRI.LIZAFM(1,15) = RETUY 277 : IMATRI.LIZAFM(1,16) = RETRET 278 : C 279 : DO IFAC = 1, NFAC, 1 280 : NGCF = MELEFE.NUM(2,IFAC) 281 : NLCF = MLENTF.LECT(NGCF) 282 : IF(NLCF .NE. IFAC)THEN 283 : WRITE(IOIMP,*) 'Il ne faut pas jouer avec la table domaine' 283 : CALL ERREUR(5) 284 : GOTO 9999 285 : ENDIF 286 : NGCG = MELEFE.NUM(1,IFAC) 287 : NGCD = MELEFE.NUM(3,IFAC) 288 : NLFL = MLELIM.LECT(NGCF) 289 : IF(NLFL .NE. 0)THEN 290 : C 291 : C********** The point belongs on BC -> No contribution to jacobian! 292 : C 293 : MELEDU.NUM(1,IFAC) = NGCG 294 : MELEDU.NUM(2,IFAC) = NGCD 295 : 296 : ELSEIF(NGCG .NE. NGCD)THEN 297 : C 298 : C********** Les MELEMEs 299 : C 300 : MELEDU.NUM(1,IFAC) = NGCG 301 : MELEDU.NUM(2,IFAC) = NGCD 302 : C 303 : C********** Les etats G et D 304 : C 305 : NLCG = MLENTC.LECT(NGCG) 306 : NLCD = MLENTC.LECT(NGCD) 307 : C 308 : ROG = MPRN.VPOCHA(NLCG,1) 309 : PG = MPPN.VPOCHA(NLCG,1) 310 : UXG = MPUN.VPOCHA(NLCG,1) 311 : UYG = MPUN.VPOCHA(NLCG,2) 312 : GAMG = MPGAMN.VPOCHA(NLCG,1) 313 : VOLG = MPVOLU.VPOCHA(NLCG,1) 314 : C 315 : ROD = MPRN.VPOCHA(NLCD,1) 316 : PD = MPPN.VPOCHA(NLCD,1) 317 : UXD = MPUN.VPOCHA(NLCD,1) 318 : UYD = MPUN.VPOCHA(NLCD,2) 319 : GAMD = MPGAMN.VPOCHA(NLCD,1) 320 : VOLD = MPVOLU.VPOCHA(NLCD,1) 321 : C 322 : C********** La normale G->D 323 : C La tangente 324 : C 325 : SURF = MPOVSU.VPOCHA(NLCF,1) 326 : CNX = MPNORM.VPOCHA(NLCF,1) 327 : CNY = MPNORM.VPOCHA(NLCF,2) 328 : CTX = -1.0D0 * CNY 329 : CTY = CNX 330 : C 331 : C********** Cut off 332 : C 333 : VINF=MPUPRI.VPOCHA(NLCG,1) 334 : VINF=MAX(VINF,MPUPRI.VPOCHA(NLCD,1)) 335 : VINF=MAX(VINF,MPUINF.VPOCHA(NLCG,1)) 336 : VINF=MAX(VINF,MPUINF.VPOCHA(NLCD,1)) 337 : C 338 : C********** L'etat moyen au centre 339 : C 340 : UNG=UXG*CNX+UYG*CNY 341 : UTG=UXG*CTX+UYG*CTY 342 : UND=UXD*CNX+UYD*CNY 343 : UTD=UXD*CTX+UYD*CTY 344 : UNF=0.5D0*(UNG+UND) 345 : UTF=0.5D0*(UTG+UTD) 346 : ROF=0.5D0*(ROG+ROD) 347 : PF=0.5D0*(PG+PD) 348 : GAMF=0.5D0*(GAMG+GAMD) 349 : C 350 : C********** We include in the cut-off UNG,UND,UTG,UTD in order to 351 : C avoid low diffusivity in stagnation regions 352 : C 353 : VINF=MAX(VINF,((UNG*UNG+UTG*UTG)**0.5D0)) 354 : VINF=MAX(VINF,((UND*UND+UTD*UTD)**0.5D0)) 355 : C 356 : C********** Calcule de la matrice de diffusivité pour le variables 357 : C conservatives. 358 : C Dans le repaire (n,t) local 359 : C F(UL,Ur) = 0.5 (F(UL)+F(UR)+AMAT*RLAMMA*(UL-UR)) 360 : C 360 : CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA) 361 : C 362 : C Ici U=(\rho,\rho un, \rho ut, \rho et) 363 : C Il faut calculer la contribution sur F dans le repaire 364 : C (n,t) par rapport aux variables 365 : C U=(\rho,\rho ux, \rho uy, \rho et) 366 : C 367 : DO I1 = 1,4,1 368 : BMAT(I1,1) = AMAT(I1,1) 369 : BMAT(I1,2) = AMAT(I1,2) * CNX + AMAT(I1,3) * CTX 370 : BMAT(I1,3) = AMAT(I1,2) * CNY + AMAT(I1,3) * CTY 371 : BMAT(I1,4) = AMAT(I1,4) 372 : ENDDO 373 : C 374 : C********** La contribution de Gauche 375 : C 376 : C Partie centrée 377 : C 377 : CALL CENTJ0(ROG,UXG,UYG,PG,GAMG,CNX,CNY,CTX,CTY, 378 : & DFRO,DFRUN,DFRUT,DFRET) 379 : C 380 : C********** Centrée + Decentrée 381 : C 382 : C DFRO = derivative of the numerical approx. of (\rho un) 383 : C in the frame (x,y) (centered part) 384 : C BMAT(1,*) = derivative of the numerical approx. of (\rho un) 385 : C in the frame (x,y) (diff. part) 386 : C 387 : DO I1=1,4,1 388 : DFRO(I1) = DFRO(I1) + 0.5D0 * RLAMMA * BMAT(1,I1) 389 : DFRUN(I1) = DFRUN(I1) + 0.5D0 * RLAMMA * BMAT(2,I1) 390 : DFRUT(I1) = DFRUT(I1) + 0.5D0 * RLAMMA * BMAT(3,I1) 391 : DFRET(I1) = DFRET(I1) + 0.5D0 * RLAMMA * BMAT(4,I1) 392 : ENDDO 393 : C 394 : C 395 : C********** AB.AM(IFAC,IPRIM,IDUAL) 396 : C A = nom de l'inconnu duale (Ro,rUX,rUY,RET) 397 : C B = nom de l'inconnu primale (Ro,rUX,rUY,RET) 398 : C IPRIM = 1, 2 -> G, D 399 : C IDUAL = 1, 2 -> G, D 400 : C i.e. 401 : C A_IDUAL = AB.AM(IFAC,IPRIM,IDUAL) * B_IPRIM + ... 402 : C 403 : C 404 : C********** Dual RN 405 : C 406 : FUNCEL = SURF * DFRO(1) 407 : RR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 408 : RR.AM(IFAC,1,2) = FUNCEL / VOLD 409 : C 410 : FUNCEL = SURF * DFRO(2) 411 : RUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 412 : RUX.AM(IFAC,1,2) = FUNCEL / VOLD 413 : C 414 : FUNCEL = SURF * DFRO(3) 415 : RUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 416 : RUY.AM(IFAC,1,2) = FUNCEL / VOLD 417 : C 418 : FUNCEL = SURF * DFRO(4) 419 : RRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 420 : RRET.AM(IFAC,1,2) = FUNCEL / VOLD 421 : C 422 : C********** Dual RUXN 423 : C 424 : FUNCEL = SURF * (DFRUN(1) * CNX + DFRUT(1) * CTX) 425 : UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 426 : UXR.AM(IFAC,1,2) = FUNCEL / VOLD 427 : C 428 : FUNCEL = SURF * (DFRUN(2) * CNX + DFRUT(2) * CTX) 429 : UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 430 : UXUX.AM(IFAC,1,2) = FUNCEL / VOLD 431 : C 432 : FUNCEL = SURF * (DFRUN(3) * CNX + DFRUT(3) * CTX) 433 : UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 434 : UXUY.AM(IFAC,1,2) = FUNCEL / VOLD 435 : C 436 : FUNCEL = SURF * (DFRUN(4) * CNX + DFRUT(4) * CTX) 437 : UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 438 : UXRET.AM(IFAC,1,2) = FUNCEL / VOLD 439 : C 440 : C********** Dual RUYN 441 : C 442 : FUNCEL = SURF * (DFRUN(1) * CNY + DFRUT(1) * CTY) 443 : UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 444 : UYR.AM(IFAC,1,2) = FUNCEL / VOLD 445 : C 446 : FUNCEL = SURF * (DFRUN(2) * CNY + DFRUT(2) * CTY) 447 : UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 448 : UYUX.AM(IFAC,1,2) = FUNCEL / VOLD 449 : C 450 : FUNCEL = SURF * (DFRUN(3) * CNY + DFRUT(3) * CTY) 451 : UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 452 : UYUY.AM(IFAC,1,2) = FUNCEL / VOLD 453 : C 454 : FUNCEL = SURF * (DFRUN(4) * CNY + DFRUT(4) * CTY) 455 : UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 456 : UYRET.AM(IFAC,1,2) = FUNCEL / VOLD 457 : C 458 : C********** Dual RETN 459 : C 460 : FUNCEL = SURF * DFRET(1) 461 : RETR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 462 : RETR.AM(IFAC,1,2) = FUNCEL / VOLD 463 : C 464 : FUNCEL = SURF * DFRET(2) 465 : RETUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 466 : RETUX.AM(IFAC,1,2) = FUNCEL / VOLD 467 : C 468 : FUNCEL = SURF * DFRET(3) 469 : RETUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 470 : RETUY.AM(IFAC,1,2) = FUNCEL / VOLD 471 : C 472 : FUNCEL = SURF * DFRET(4) 473 : RETRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 474 : RETRET.AM(IFAC,1,2) = FUNCEL / VOLD 475 : C 476 : C 477 : C********** La contribution de D 478 : C Je change l'orientation de normales 479 : C Dans ce cas F(UD,UG) = 0.5 * (F(UD)+F(UD) 480 : C + RLAMMA * AMAT * (UD - UG)) 481 : C 482 : CNX = -1.0D0 * CNX 483 : CNY = -1.0D0 * CNY 484 : CTX = -1.0D0 * CTX 485 : CTY = -1.0D0 * CTY 486 : C 487 : C********** L'etat moyen au centre 488 : C Meme que avant; mais UNF, UTF change de signe 489 : C 490 : UNF=-1.0D0*UNF 491 : UTF=-1.0D0*UTF 492 : C 493 : C********** On recalcule de la matrice de diffusivité pour le variables 494 : C conservatives (pas d'optimisation illisible!!!) 495 : C 495 : CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA) 496 : C 497 : C Ici U=(\rho,\rho un, \rho ut, \rho et) 498 : C Il faut calculer la contribution sur F dans le repaire 499 : C (n,t) par rapport aux variables 500 : C U=(\rho,\rho ux, \rho uy, \rho et) 501 : C 502 : DO I1 = 1,4,1 503 : BMAT(I1,1) = AMAT(I1,1) 504 : BMAT(I1,2) = AMAT(I1,2) * CNX + AMAT(I1,3) * CTX 505 : BMAT(I1,3) = AMAT(I1,2) * CNY + AMAT(I1,3) * CTY 506 : BMAT(I1,4) = AMAT(I1,4) 507 : ENDDO 508 : C 508 : CALL CENTJ0(ROD,UXD,UYD,PD,GAMD,CNX,CNY,CTX,CTY, 509 : & DFRO,DFRUN,DFRUT,DFRET) 510 : C 511 : C********** Centrée + Decentrée 512 : C 513 : DO I1=1,4,1 514 : DFRO(I1) = DFRO(I1) + 0.5D0 * RLAMMA * BMAT(1,I1) 515 : DFRUN(I1) = DFRUN(I1) + 0.5D0 * RLAMMA * BMAT(2,I1) 516 : DFRUT(I1) = DFRUT(I1) + 0.5D0 * RLAMMA * BMAT(3,I1) 517 : DFRET(I1) = DFRET(I1) + 0.5D0 * RLAMMA * BMAT(4,I1) 518 : ENDDO 519 : C 520 : C 521 : C********** Dual RN 522 : C 523 : FUNCEL = SURF * DFRO(1) 524 : RR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 525 : RR.AM(IFAC,2,1) = FUNCEL / VOLG 526 : C 527 : FUNCEL = SURF * DFRO(2) 528 : RUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 529 : RUX.AM(IFAC,2,1) = FUNCEL / VOLG 530 : C 531 : FUNCEL = SURF * DFRO(3) 532 : RUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 533 : RUY.AM(IFAC,2,1) = FUNCEL / VOLG 534 : C 535 : FUNCEL = SURF * DFRO(4) 536 : RRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 537 : RRET.AM(IFAC,2,1) = FUNCEL / VOLG 538 : C 539 : C********** Dual RUXN 540 : C 541 : FUNCEL = SURF * (DFRUN(1) * CNX + DFRUT(1) * CTX) 542 : UXR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 543 : UXR.AM(IFAC,2,1) = FUNCEL / VOLG 544 : C 545 : FUNCEL = SURF * (DFRUN(2) * CNX + DFRUT(2) * CTX) 546 : UXUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 547 : UXUX.AM(IFAC,2,1) = FUNCEL / VOLG 548 : C 549 : FUNCEL = SURF * (DFRUN(3) * CNX + DFRUT(3) * CTX) 550 : UXUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 551 : UXUY.AM(IFAC,2,1) = FUNCEL / VOLG 552 : C 553 : FUNCEL = SURF * (DFRUN(4) * CNX + DFRUT(4) * CTX) 554 : UXRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 555 : UXRET.AM(IFAC,2,1) = FUNCEL / VOLG 556 : C 557 : C********** Dual RUYN 558 : C 559 : FUNCEL = SURF * (DFRUN(1) * CNY + DFRUT(1) * CTY) 560 : UYR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 561 : UYR.AM(IFAC,2,1) = FUNCEL / VOLG 562 : C 563 : FUNCEL = SURF * (DFRUN(2) * CNY + DFRUT(2) * CTY) 564 : UYUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 565 : UYUX.AM(IFAC,2,1) = FUNCEL / VOLG 566 : C 567 : FUNCEL = SURF * (DFRUN(3) * CNY + DFRUT(3) * CTY) 568 : UYUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 569 : UYUY.AM(IFAC,2,1) = FUNCEL / VOLG 570 : C 571 : FUNCEL = SURF * (DFRUN(4) * CNY + DFRUT(4) * CTY) 572 : UYRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 573 : UYRET.AM(IFAC,2,1) = FUNCEL / VOLG 574 : C 575 : C********** Dual RETN 576 : C 577 : FUNCEL = SURF * DFRET(1) 578 : RETR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 579 : RETR.AM(IFAC,2,1) = FUNCEL / VOLG 580 : C 581 : FUNCEL = SURF * DFRET(2) 582 : RETUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 583 : RETUX.AM(IFAC,2,1) = FUNCEL / VOLG 584 : C 585 : FUNCEL = SURF * DFRET(3) 586 : RETUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 587 : RETUY.AM(IFAC,2,1) = FUNCEL / VOLG 588 : C 589 : FUNCEL = SURF * DFRET(4) 590 : RETRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD 591 : RETRET.AM(IFAC,2,1) = FUNCEL / VOLG 592 : C 593 : ELSE 594 : C 595 : C********** Murs (NGCG = NGCD) 596 : C 597 : C 598 : C********** Les MELEMEs 599 : C 600 : MELEDU.NUM(1,IFAC) = NGCG 601 : MELEDU.NUM(2,IFAC) = NGCD 602 : NLCG = MLENTC.LECT(NGCG) 603 : C 604 : ROG = MPRN.VPOCHA(NLCG,1) 605 : PG = MPPN.VPOCHA(NLCG,1) 606 : UXG = MPUN.VPOCHA(NLCG,1) 607 : UYG = MPUN.VPOCHA(NLCG,2) 608 : GAMG = MPGAMN.VPOCHA(NLCG,1) 609 : VOLG = MPVOLU.VPOCHA(NLCG,1) 610 : C 611 : C********** Cut off 612 : C 613 : VINF=MPUPRI.VPOCHA(NLCG,1) 614 : VINF=MAX(VINF,MPUINF.VPOCHA(NLCG,1)) 615 : C 616 : C********** La normale sortante 617 : C 618 : SURF = MPOVSU.VPOCHA(NLCF,1) 619 : CNX = MPNORM.VPOCHA(NLCF,1) 620 : CNY = MPNORM.VPOCHA(NLCF,2) 621 : C 622 : C********** L'etat moyen au centre 623 : C 624 : UNF=0.0D0 625 : UTF=0.5D0*(UXG*CTX+UYG*CTY) 626 : ROF=ROG 627 : PF=PG 628 : GAMF=GAMG 629 : C 630 : C********** We include in the cut-off UG 631 : C 632 : VINF=MAX(VINF,((UXG*UXG+UYG*UYG)**0.5D0)) 633 : C 634 : C********** Calcule de la matrice de diffusivité pour le variables 635 : C conservatives. 636 : C Dans le repaire (n,t) local 637 : C F(UL,Ur) = 0.5 (F(UL)+F(UR)+AMAT*RLAMMA*(UL-UR)) 638 : C 638 : CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA) 639 : C 640 : C Ici U=(\rho,\rho un, \rho ut, \rho et) 641 : C Il faut calculer la contribution sur F dans le repaire 642 : C (n,t) par rapport aux variables 643 : C U=(\rho,\rho ux, \rho uy, \rho et) 644 : C Dans le cas murs, le seul flux normal donne une contribution 645 : C non nulle 646 : C En plus AMAT(1,2)=AMAT(3,2)=AMAT(4,2)=0 647 : C AMAT(2,2)=1.0D0 648 : C 649 : C write(*,*) 'AMAT(*,2)=',(AMAT(I1,2),I1=1,4) 650 : C 651 : C********** Contribution centrée 652 : C 652 : CALL CENTJ2(ROG,UXG,UYG,GAMG,CNX,CNY, 653 : & DFRUN) 654 : C 655 : C********** Centrée + Decentrée 656 : C 657 : DFRUN(2) = DFRUN(2) + (RLAMMA * CNX) 658 : DFRUN(3) = DFRUN(3) + (RLAMMA * CNY) 659 : C 660 : C********** Dual RN 661 : C 662 : RR.AM(IFAC,1,1) = 0.0D0 663 : RR.AM(IFAC,1,2) = 0.0D0 664 : C 665 : RUX.AM(IFAC,1,1) = 0.0D0 666 : RUX.AM(IFAC,1,2) = 0.0D0 667 : C 668 : RUY.AM(IFAC,1,1) = 0.0D0 669 : RUY.AM(IFAC,1,2) = 0.0D0 670 : C 671 : RRET.AM(IFAC,1,1) = 0.0D0 672 : RRET.AM(IFAC,1,2) = 0.0D0 673 : C 674 : C********** Dual RUXN 675 : C 676 : FUNCEL = SURF * DFRUN(1) * CNX 677 : UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 678 : UXR.AM(IFAC,1,2) = 0.0D0 679 : C 680 : FUNCEL = SURF * DFRUN(2) * CNX 681 : UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 682 : UXUX.AM(IFAC,1,2) = 0.0D0 683 : C 684 : FUNCEL = SURF * DFRUN(3) * CNX 685 : UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 686 : UXUY.AM(IFAC,1,2) = 0.0D0 687 : C 688 : FUNCEL = SURF * DFRUN(4) * CNX 689 : UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 690 : UXRET.AM(IFAC,1,2) = 0.0D0 691 : C 692 : C********** Dual RUYN 693 : C 694 : FUNCEL = SURF * DFRUN(1) * CNY 695 : UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 696 : UYR.AM(IFAC,1,2) = 0.0D0 697 : C 698 : FUNCEL = SURF * DFRUN(2) * CNY 699 : UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 700 : UYUX.AM(IFAC,1,2) = 0.0D0 701 : C 702 : FUNCEL = SURF * DFRUN(3) * CNY 703 : UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 704 : UYUY.AM(IFAC,1,2) = 0.0D0 705 : C 706 : FUNCEL = SURF * DFRUN(4) * CNY 707 : UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG 708 : UYRET.AM(IFAC,1,2) = 0.0D0 709 : C 710 : C********** Dual RETN 711 : C 712 : RETR.AM(IFAC,1,1) = 0.0D0 713 : RETR.AM(IFAC,1,2) = 0.0D0 714 : C 715 : RETUX.AM(IFAC,1,1) = 0.0D0 716 : RETUX.AM(IFAC,1,2) = 0.0D0 717 : C 718 : RETUY.AM(IFAC,1,1) = 0.0D0 719 : RETUY.AM(IFAC,1,2) = 0.0D0 720 : C 721 : RETRET.AM(IFAC,1,1) = 0.0D0 722 : RETRET.AM(IFAC,1,2) = 0.0D0 723 : C 724 : C********** Dual RN 725 : C 726 : RR.AM(IFAC,2,2) = 0.0D0 727 : RR.AM(IFAC,2,1) = 0.0D0 728 : C 729 : RUX.AM(IFAC,2,2) = 0.0D0 730 : RUX.AM(IFAC,2,1) = 0.0D0 731 : C 732 : RUY.AM(IFAC,2,2) = 0.0D0 733 : RUY.AM(IFAC,2,1) = 0.0D0 734 : C 735 : RRET.AM(IFAC,2,2) = 0.0D0 736 : RRET.AM(IFAC,2,1) = 0.0D0 737 : C 738 : C********** Dual RUXN 739 : C 740 : UXR.AM(IFAC,2,2) = 0.0D0 741 : UXR.AM(IFAC,2,1) = 0.0D0 742 : C 743 : UXUX.AM(IFAC,2,2) = 0.0D0 744 : UXUX.AM(IFAC,2,1) = 0.0D0 745 : C 746 : UXUY.AM(IFAC,2,2) = 0.0D0 747 : UXUY.AM(IFAC,2,1) = 0.0D0 748 : C 749 : UXRET.AM(IFAC,2,2) = 0.0D0 750 : UXRET.AM(IFAC,2,1) = 0.0D0 751 : C 752 : C********** Dual RUYN 753 : C 754 : UYR.AM(IFAC,2,2) = 0.0D0 755 : UYR.AM(IFAC,2,1) = 0.0D0 756 : C 757 : UYUX.AM(IFAC,2,2) = 0.0D0 758 : UYUX.AM(IFAC,2,1) = 0.0D0 759 : C 760 : UYUY.AM(IFAC,2,2) = 0.0D0 761 : UYUY.AM(IFAC,2,1) = 0.0D0 762 : C 763 : UYRET.AM(IFAC,2,2) = 0.0D0 764 : UYRET.AM(IFAC,2,1) = 0.0D0 765 : C 766 : C********** Dual RETN 767 : C 768 : RETR.AM(IFAC,2,2) = 0.0D0 769 : RETR.AM(IFAC,2,1) = 0.0D0 770 : C 771 : RETUX.AM(IFAC,2,2) = 0.0D0 772 : RETUX.AM(IFAC,2,1) = 0.0D0 773 : C 774 : RETUY.AM(IFAC,2,2) = 0.0D0 775 : RETUY.AM(IFAC,2,1) = 0.0D0 776 : C 777 : RETRET.AM(IFAC,2,2) = 0.0D0 778 : RETRET.AM(IFAC,2,1) = 0.0D0 779 : C 780 : ENDIF 781 : ENDDO 782 : C 783 : SEGDES MELEMC 784 : SEGDES MELEFE 785 : SEGDES MELEMF 786 : C 787 : SEGDES MPOVSU 788 : SEGDES MPVOLU 789 : SEGDES MPNORM 790 : C 791 : SEGDES MPRN 792 : SEGDES MPPN 793 : SEGDES MPUN 794 : SEGDES MPGAMN 795 : C 796 : SEGDES MELEDU 797 : SEGDES MATRIK 798 : SEGDES IMATRI 799 : C 800 : SEGDES RR , RUX , RUY , RRET , 801 : & UXR , UXUX , UXUY , UXRET , 802 : & UYR , UYUX , UYUY , UYRET , 803 : & RETR , RETUX , RETUY , RETRET 804 : 805 : SEGSUP MLENTC 806 : SEGSUP MLENTF 807 : SEGDES MLMINC 808 : SEGSUP MLELIM 809 : C 810 : IF(MELLIM .NE.0) SEGDES MELLIM 811 : C 812 : SEGDES MPUPRI 813 : SEGDES MPUINF 814 : C 815 : 9999 CONTINUE 816 : RETURN 817 : END 818 : 819 : 820 : 821 : 822 :
© Cast3M 2003 - All rights reserved.
Disclaimer