C KONJA5 SOURCE OF166741 24/12/13 21:16:32 12097 SUBROUTINE KONJA5(IIMPL,ILINC,IROF,IPF,IVITF,IGAMF, & IRN,IPN,IVN,IGAMN, $ ICHPVO,ICHPSU, ICHPNO, MELEMC, MELEFE, MELEMF, MELTFA, & IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN, & ICRN, ICVN, ICPN, IVLIM,MELLIM,IMAT) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : KONJA5 C C DESCRIPTION : Voir KONV11 C Calcul du jacobien du résidu "2nd order accurate" C C Cas deux dimensions, gaz "calorically perfect" C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF C C************************************************************************ C C APPELES (Outils C CASTEM) : KRIPAD, LICHT, ERREUR C C APPELES (Calcul) : VLHJ3, AUSMW, VLHJ1, AUSMF C C************************************************************************ C C ENTREES C C IIMPL : INTEGER C = 4 -> methode VLH C = 5 -> methode AUSM+ C C ILINC : liste des inconnues (pointeur d'un objet de type C LISTMOTS) C C 1) Pointeurs des CHPOINT/CHAMELEM C C IROF : MCHAML 'FACEL' densité C C IPF : MCHAML 'FACEL' pression C C IVITF : MCHAML 'FACEL' vitesse/normales C C IGAMF : MCHAML 'FACEL' gamma C C IRN : CHPOINT 'CENTRE' densité C C IPN : CHPOINT 'CENTRE' pression C C IVN : CHPOINT 'CENTRE' vitesse C C IGAMN : CHPOINT 'CENTRE' gamma C C ICHPVO : CHPOINT VOLUME contenant le volume C C ICHPSU : CHPOINT FACE contenant la surface des faces C C ICHPNO : CHPOINT 'FACE' contenant les normales aux faces C C IGRN : CHPOINT 'CENTRE' gradient de densité C C IGVN : CHPOINT 'CENTRE' gradient de vitesse C C IGPN : CHPOINT 'CENTRE' gradient de pression C C IGLRN : CHPOINT 'CENTRE' limiteur du gradient de densité C C IGLVN : CHPOINT 'CENTRE' limiteur du gradient de vitesse C C IGLPN : CHPOINT 'CENTRE' limiteur du gradient de pression C C ICRN : MCHAML contenant les coeff. pour calculer le gradient C de densité C C ICVN : MCHAML contenant les coeff. pour calculer le gradient C de vitesse C C ICPN : MCHAML contenant les coeff. pour calculer le gradient C de pression C C IVLIM : CHPOINT des conditions aux limites sur la vitesse C C 2) Pointeurs de MELEME de la table DOMAINE C C MELEMC : MELEME 'CENTRE', SPG des CENTRES des elts C C MELEFE : MELEME 'FACEL' (centre elt gauche, centre de face C centre elt droite) C C MELEMF : MELEME 'FACE', SPG des CENTRES des faces C C MELTFA : MELEME 'ELTFA' (centre d'elt, centres de faces) C C MELLIM : MELEME SPG des conditions aux bords C C SORTIES C C IMAT : pointeur de la MATRIK du jacobien du residu C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : C C************************************************************************ C C N.B.: On suppose qu'on a déjà controllé RO, P > 0 C GAMMA \in (1,3) C Si non il faut le faire!!! C C************************************************************************ C IMPLICIT INTEGER(I-N) INTEGER IIMPL,ILINC, ICHPVO, ICHPNO, ICHPSU & ,IROF,IVITF,IPF,IGAMF, IRN,IPN,IVN,IGAMN & , IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN & , ICRN, ICVN, ICPN, IVLIM , IMAT & , NBSOUP, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID & , NKMT, NBME, NBEL, MP, NP, NP1 & , IFAC, NGCF, NLCF, NGCF1, NGCG, NGCD & , NGC, NLC, NLCD, NGNO, NGN1,NGN2, NLNO, IELEM, IELPRI & , JG, NBSOUS, NBSOU1, ISOUS, NBEL1 & , ISPG1, ISPG2 & , I1, IGEOMC, IGEOMF, IGEOML & , IFACGR, NGGRA, NLGRA, NGFGR, NLFGR, NG, ND, ICOOR, IFAC2 & , NLVIMP & , NGPRI,NLPRI,IELETR,IDUATR,IPRITR, NLFL REAL*8 ROG, PG, UNG, UTG, UXG, UYG, RETG, GAMG, VOLG, VOLD & , ROD, PD, UND, UTD, UXD, UYD & , SURF, CNX, CNY, CTX, CTY, ECIN, FUNCEL & , DFRO(4), DFRET(4), DFRUN(4), DFRUT(4),DFRUX(4),DFRUY(4) & , ALR, ALP, ALVX, ALVY, XCG, YCG, DXCF, DYCF, CELLR, CELLP & , CVXVX, CVXVY, CVYVX, CVYVY, CNX1, CNY1, CTX1, CTY1 & , DDP, DDRO, DDUX, DDUY & , GAMGM1, CELCO1, CELCO2, CELCO3, CELCO4, CELCO5 CHARACTER*8 TYPE C C**** LES INCLUDES C -INC SMCOORD -INC PPARAM -INC CCOPTIO -INC SMCHAML POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, & MELT1X.MELVAL, MELT1Y.MELVAL POINTEUR MELVUN.MELVAL, MELVUT.MELVAL POINTEUR MELRO.MELVAL, MELP.MELVAL, & MELGAM.MELVAL POINTEUR MELCR1.MELVAL, MELCR2.MELVAL, MELCP1.MELVAL, & MELCP2.MELVAL, MELCV1.MELVAL, MELCV2.MELVAL C -INC SMCHPOI -INC SMELEME -INC SMLMOTS -INC SMLENTI -INC SMMATRIK POINTEUR MPVOLU.MPOVAL, MPNORM.MPOVAL, MPOVSU.MPOVAL, $ MPVGRN.MPOVAL,MPGLRN.MPOVAL, MPVGPN.MPOVAL, MPGLPN.MPOVAL, $ MPVGVN.MPOVAL,MPGLVN.MPOVAL, MPVLIM.MPOVAL,MPOVRN.MPOVAL,MPOVVN $ .MPOVAL,MPOVPN.MPOVAL,MPOGAM.MPOVAL POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME, & MELTFA.MELEME,MELPRI.MELEME,MELGRR.MELEME,MELGRP.MELEME $ ,MELGRV.MELEME, MELLIM.MELEME POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLESPG.MLENTI,MLEPRI.MLENTI & ,MLEVL.MLENTI,MLEGRR.MLENTI,MLEGRP.MLENTI,MLEGRV.MLENTI & ,MLELIM.MLENTI POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RRET.IZAFM, & UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXRET.IZAFM, & UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYRET.IZAFM, & RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETRET.IZAFM POINTEUR MLMINC.MLMOTS C C**** KRIPAD pour la correspondance global/local des conditions limits C CALL KRIPAD(MELLIM,MLELIM) C SEGACT MELLIM C C**** Masse volumique C MCHEL1 = IROF SEGACT MCHEL1 MCHAM1 = MCHEL1.ICHAML(1) SEGACT MCHAM1 MELRO = MCHAM1.IELVAL(1) SEGDES MCHEL1 SEGDES MCHAM1 C C**** Pression C MCHEL1 = IPF SEGACT MCHEL1 MCHAM1 = MCHEL1.ICHAML(1) SEGACT MCHAM1 MELP = MCHAM1.IELVAL(1) SEGDES MCHEL1 SEGDES MCHAM1 C C**** Gamma C MCHEL1 = IGAMF SEGACT MCHEL1 MCHAM1 = MCHEL1.ICHAML(1) SEGACT MCHAM1 MELGAM = MCHAM1.IELVAL(1) SEGDES MCHEL1 SEGDES MCHAM1 C C**** Vitesse et cosinus directeurs du repere (n,t) C MCHEL1 = IVITF SEGACT MCHEL1 C C**** La vitesse a comme SPG MELEFE C Le cosinus directeurs ont comme SPG MELEMF C C MCHAM1 -> Cosinus directeurs C MCHAM2 -> Vitesse C ISPG1 = MCHEL1.IMACHE(1) ISPG2 = MCHEL1.IMACHE(2) IF((ISPG1 .EQ. MELEMF) .AND. (ISPG2 .EQ. MELEFE))THEN MCHAM1 = MCHEL1.ICHAML(1) MCHAM2 = MCHEL1.ICHAML(2) ELSEIF((ISPG1 .EQ. MELEFE) .AND. (ISPG2 .EQ. MELEMF))THEN MCHAM1 = MCHEL1.ICHAML(2) MCHAM2 = MCHEL1.ICHAML(1) ELSE WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF SEGACT MCHAM1 MELVNX = MCHAM1.IELVAL(1) MELVNY = MCHAM1.IELVAL(2) MELT1X = MCHAM1.IELVAL(3) MELT1Y = MCHAM1.IELVAL(4) SEGDES MCHAM1 SEGACT MCHAM2 MELVUN = MCHAM2.IELVAL(1) MELVUT = MCHAM2.IELVAL(2) SEGDES MCHAM2 SEGDES MCHEL1 C C**** Activation des MCHAMLs C SEGACT MELRO SEGACT MELP SEGACT MELGAM SEGACT MELVUN SEGACT MELVUT SEGACT MELVNX SEGACT MELVNY SEGACT MELT1X SEGACT MELT1Y C C**************** CHPOINT gradients, limiteurs, vitesse aux bords C CALL LICHT(IRN,MPOVRN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOVRN CALL LICHT(IVN,MPOVVN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOVVN CALL LICHT(IPN,MPOVPN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOVPN CALL LICHT(IGAMN,MPOGAM,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOGAM CALL LICHT(IGRN,MPVGRN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPVGRN CALL LICHT(IGLRN,MPGLRN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPGLRN CALL LICHT(IGPN,MPVGPN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPVGPN CALL LICHT(IGLPN,MPGLPN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPGLPN CALL LICHT(IGVN,MPVGVN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPVGPN CALL LICHT(IGLVN,MPGLVN,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C SEGACT MPGLVN CALL LICHT(IVLIM,MPVLIM,TYPE,IGEOML) IF(IERR.NE.0)GOTO 9999 C SEGACT MPVLIM C C**************** KRIPAD pour la correspondance global/local des centres d'elts C CALL KRIPAD(MELEMC,MLENTC) IF(IERR.NE.0)GOTO 9999 C C SEGACT MLENTC SEGACT MELEMC C C**** KRIPAD pour la correspondance global/local des centres des faces C CALL KRIPAD(MELEMF,MLENTF) IF(IERR.NE.0)GOTO 9999 C C SEGACT MLENTF C C**** KRIPAD pour la correspondance global/local des B.C. sur le vitesse C CALL KRIPAD(IGEOML,MLEVL) IF(IERR.NE.0)GOTO 9999 C C SEGACT MLEVL C SEGACT MELEFE C CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF) IF(IERR.NE.0)GOTO 9999 CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 CALL LICHT(ICHPNO,MPNORM,TYPE,IGEOMC) IF(IERR.NE.0)GOTO 9999 C C**** LICHT active les MPOVALs en *MOD C C i.e. C C SEGACT MPOVSU*MOD C SEGACT MPVOLU*MOD C SEGACT MPNORM*MOD C C C**** Les gradients ont tous le meme SPG (mais le numero C du pointeur peut-etre different). C Ils ont le meme nombre des SPGs que MELTFA C En plus le SPGs ont le meme ordre: C le iem elt de MELTFA correspond a l'elt qui a C comme centre le iem noeud de CENTRE C SEGACT MELTFA NBSOUP=MELTFA.LISOUS(/1) IF(NBSOUP.EQ.0) NBSOUP=1 JG=NBSOUP SEGINI MLESPG IF(NBSOUP .EQ. 1)THEN MLESPG.LECT(1)=MELTFA ELSE DO I1 = 1, NBSOUP, 1 MLESPG.LECT(I1)=MELTFA.LISOUS(I1) ENDDO ENDIF C MCHEL1=ICRN SEGACT MCHEL1 NBSOU1=MCHEL1.IMACHE(/1) IF(NBSOU1 .NE. NBSOUP)THEN WRITE(IOIMP,*) 'SPG coeffs grads = ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF C MCHEL2=ICPN SEGACT MCHEL2 NBSOU1=MCHEL2.IMACHE(/1) IF(NBSOU1 .NE. NBSOUP)THEN WRITE(IOIMP,*) 'SPG coeffs grads = ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF C MCHEL2=ICVN SEGACT MCHEL2 NBSOU1=MCHEL2.IMACHE(/1) IF(NBSOU1 .NE. NBSOUP)THEN WRITE(IOIMP,*) 'SPG coeffs grads = ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF C C**** N.B.: ICRN, ICPN, ICVN sont activés!!! C C C****************** L'objet matrik C C NRIGE = 7 NMATRI = 1 NKID = 9 NKMT = 7 C SEGINI MATRIK IMAT = MATRIK C C**** Maillage des inconnues duales: MELEMC C Maillage des inconnues primales: C même structure que MCHEL1.IMACHE(ISOUS) C Donc on partitionne les deux dans la meme C maniere. C JG=NBSOUP SEGINI MLEPRI SEGINI MLEGRR SEGINI MLEGRP SEGINI MLEGRV IF(NBSOUP .EQ. 1)THEN MCHEL1=ICRN IPT1=MCHEL1.IMACHE(1) SEGACT IPT1 SEGINI, MELPRI = IPT1 MLEPRI.LECT(1) = MELPRI MLEGRR.LECT(1) = IPT1 C MCHEL1=ICPN IPT1=MCHEL1.IMACHE(1) SEGACT IPT1 MLEGRP.LECT(1) = IPT1 C MCHEL1=ICVN IPT1=MCHEL1.IMACHE(1) SEGACT IPT1 MLEGRV.LECT(1) = IPT1 ELSE NBSOUS=NBSOUP NBREF=0 NBNN=0 NBELEM=0 SEGINI MELPRI MCHEL1=ICRN MCHEL2=ICPN MCHEL3=ICVN DO ISOUS = 1, NBSOUP, 1 IPT1=MCHEL1.IMACHE(ISOUS) IPT2=MCHEL2.IMACHE(ISOUS) IPT3=MCHEL3.IMACHE(ISOUS) SEGACT IPT1 SEGACT IPT2 SEGACT IPT3 SEGINI, MELEME=IPT1 MELPRI.LISOUS(ISOUS) = MELEME MLEPRI.LECT(ISOUS) = MELEME MLEGRR.LECT(ISOUS) = IPT1 MLEGRP.LECT(ISOUS) = IPT2 MLEGRV.LECT(ISOUS) = IPT3 ENDDO ENDIF MATRIK.IRIGEL(2,1) = MELPRI MATRIK.IRIGEL(1,1) = MELPRI C C**** Matrice non symetrique C MATRIK.IRIGEL(7,1) = 2 NBME = 16 NBSOUS=NBSOUP SEGINI IMATRI MLMINC = ILINC SEGACT MLMINC MATRIK.IRIGEL(4,1) = IMATRI C IMATRI.LISPRI(1) = MLMINC.MOTS(1) IMATRI.LISPRI(2) = MLMINC.MOTS(2) IMATRI.LISPRI(3) = MLMINC.MOTS(3) IMATRI.LISPRI(4) = MLMINC.MOTS(4) IMATRI.LISPRI(5) = MLMINC.MOTS(1) IMATRI.LISPRI(6) = MLMINC.MOTS(2) IMATRI.LISPRI(7) = MLMINC.MOTS(3) IMATRI.LISPRI(8) = MLMINC.MOTS(4) IMATRI.LISPRI(9) = MLMINC.MOTS(1) IMATRI.LISPRI(10) = MLMINC.MOTS(2) IMATRI.LISPRI(11) = MLMINC.MOTS(3) IMATRI.LISPRI(12) = MLMINC.MOTS(4) IMATRI.LISPRI(13) = MLMINC.MOTS(1) IMATRI.LISPRI(14) = MLMINC.MOTS(2) IMATRI.LISPRI(15) = MLMINC.MOTS(3) IMATRI.LISPRI(16) = MLMINC.MOTS(4) C IMATRI.LISDUA(1) = MLMINC.MOTS(1) IMATRI.LISDUA(2) = MLMINC.MOTS(1) IMATRI.LISDUA(3) = MLMINC.MOTS(1) IMATRI.LISDUA(4) = MLMINC.MOTS(1) IMATRI.LISDUA(5) = MLMINC.MOTS(2) IMATRI.LISDUA(6) = MLMINC.MOTS(2) IMATRI.LISDUA(7) = MLMINC.MOTS(2) IMATRI.LISDUA(8) = MLMINC.MOTS(2) IMATRI.LISDUA(9) = MLMINC.MOTS(3) IMATRI.LISDUA(10) = MLMINC.MOTS(3) IMATRI.LISDUA(11) = MLMINC.MOTS(3) IMATRI.LISDUA(12) = MLMINC.MOTS(3) IMATRI.LISDUA(13) = MLMINC.MOTS(4) IMATRI.LISDUA(14) = MLMINC.MOTS(4) IMATRI.LISDUA(15) = MLMINC.MOTS(4) IMATRI.LISDUA(16) = MLMINC.MOTS(4) C C**** Boucle sur les LISOUS de MELPRI C IELEM=0 C DO ISOUS=1,NBSOUP,1 C C******* Les coeffs des grad. C MCHEL1=ICRN MCHAM1=MCHEL1.ICHAML(ISOUS) SEGACT MCHAM1 MELCR1=MCHAM1.IELVAL(1) MELCR2=MCHAM1.IELVAL(2) SEGACT MELCR1, MELCR2 SEGDES MCHAM1 C MCHEL1=ICPN MCHAM1=MCHEL1.ICHAML(ISOUS) SEGACT MCHAM1 MELCP1=MCHAM1.IELVAL(1) MELCP2=MCHAM1.IELVAL(2) SEGACT MELCP1, MELCP2 SEGDES MCHAM1 C C MCHEL1=ICVN MCHAM1=MCHEL1.ICHAML(ISOUS) SEGACT MCHAM1 MELCV1=MCHAM1.IELVAL(1) MELCV2=MCHAM1.IELVAL(2) SEGACT MELCV1, MELCV2 SEGDES MCHAM1 C C******* MELEME : maillage elementaire de ELTFA C MELEME=MLESPG.LECT(ISOUS) SEGACT MELEME NBEL1=MELEME.NUM(/2) NP1=MELEME.NUM(/1) C C******* IPT1 : SPG des inconnues primales et duales C IPT1=MLEPRI.LECT(ISOUS) NBEL = IPT1.NUM(/2) NP = IPT1.NUM(/1) C C******* MELGRR: SPG elementaire des coeffs de GRAD C Il resemble a IPT1, mais IPT1 peut changer! C MELGRR=MLEGRR.LECT(ISOUS) MELGRP=MLEGRP.LECT(ISOUS) MELGRV=MLEGRV.LECT(ISOUS) C C******* MLESPG et MLEPRI doivent avoir le meme nombre C d'elts et de noeuds! C IF((NBEL1 .NE. NBEL) .OR. (NP1 .NE. (NP-1)))THEN WRITE(IOIMP,*) 'ELTFA, SPGs gradients' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF C MP = NP SEGINI RR , RUX , RUY , RRET , & UXR , UXUX , UXUY , UXRET , & UYR , UYUX , UYUY , UYRET , & RETR , RETUX , RETUY , RETRET C C**** Duale = IMATRI.LISDUA(1) = 'RN' C Primale = IMATRI.LISPRI(1) = 'RN' C -> IMATRI.LIZAFM(1,1) = RR C C Duale = IMATRI.LISDUA(2) = 'RN' C Primale = IMATRI.LISPRI(2) = 'RUXN' C -> IMATRI.LIZAFM(1,2) = RUX C ... C C NB Pour le moment on y stoke le jacobien du flux C par rapport aux variables primitives (R,UX,UY,P) C IMATRI.LIZAFM(ISOUS,1) = RR IMATRI.LIZAFM(ISOUS,2) = RUX IMATRI.LIZAFM(ISOUS,3) = RUY IMATRI.LIZAFM(ISOUS,4) = RRET IMATRI.LIZAFM(ISOUS,5) = UXR IMATRI.LIZAFM(ISOUS,6) = UXUX IMATRI.LIZAFM(ISOUS,7) = UXUY IMATRI.LIZAFM(ISOUS,8) = UXRET IMATRI.LIZAFM(ISOUS,9) = UYR IMATRI.LIZAFM(ISOUS,10) = UYUX IMATRI.LIZAFM(ISOUS,11) = UYUY IMATRI.LIZAFM(ISOUS,12) = UYRET IMATRI.LIZAFM(ISOUS,13) = RETR IMATRI.LIZAFM(ISOUS,14) = RETUX IMATRI.LIZAFM(ISOUS,15) = RETUY IMATRI.LIZAFM(ISOUS,16) = RETRET C C******* Boucle sur les elts de MLEPRI.LECT C DO IELPRI=1,NBEL,1 IELEM=IELEM+1 NGC=IPT1.NUM(NP,IELPRI) NLC=MLENTC.LECT(NGC) IF(NLC.NE.IELEM)THEN WRITE(IOIMP,*) 'NLC != IELEM' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF VOLG = MPVOLU.VPOCHA(NLC,1) C C************* Pour la contribution des gradients C 1) les limiteurs C ALR = MPGLRN.VPOCHA(NLC,1) ALP = MPGLPN.VPOCHA(NLC,1) ALVX = MPGLVN.VPOCHA(NLC,1) ALVY = MPGLVN.VPOCHA(NLC,2) ICOOR= (NGC - 1) * 3 XCG = MCOORD.XCOOR(ICOOR + 1) YCG = MCOORD.XCOOR(ICOOR + 2) C C******** Boucle sur le noeuds de MLESPG.LECT C DO IFAC=1,(NP - 1),1 NGCF = MELEME.NUM(IFAC,IELPRI) ICOOR= (NGCF - 1) * 3 DXCF = MCOORD.XCOOR(ICOOR + 1) - XCG DYCF = MCOORD.XCOOR(ICOOR + 2) - YCG NLCF = MLENTF.LECT(NGCF) NGCF1 = MELEFE.NUM(2,NLCF) IF(NGCF .NE. NGCF1)THEN WRITE(IOIMP,*) $ 'FACE ET FACEL = ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF NGCG = MELEFE.NUM(1,NLCF) NGCD = MELEFE.NUM(3,NLCF) SURF = MPOVSU.VPOCHA(NLCF,1) NLFL = MLELIM.LECT(NGCF) IF(NLFL .NE. 0)THEN C C************* The point belongs on BC -> No contribution to jacobian! C DFRO(1)=0.0D0 DFRO(2)=0.0D0 DFRO(3)=0.0D0 DFRO(4)=0.0D0 DFRUN(1)=0.0D0 DFRUN(2)=0.0D0 DFRUN(3)=0.0D0 DFRUN(4)=0.0D0 DFRUT(1)=0.0D0 DFRUT(2)=0.0D0 DFRUT(3)=0.0D0 DFRUT(4)=0.0D0 DFRET(1)=0.0D0 DFRET(2)=0.0D0 DFRET(3)=0.0D0 DFRET(4)=0.0D0 C C************* On est sur G, sur D, ou sur un mur? C ELSEIF(NGCG .EQ. NGCD)THEN C C***************** Murs C Aucune contribution à l'elt droite! C CNX = MELVNX.VELCHE(1,NLCF) CNY = MELVNY.VELCHE(1,NLCF) CTX = MELT1X.VELCHE(1,NLCF) CTY = MELT1Y.VELCHE(1,NLCF) ROG = MELRO.VELCHE(1,NLCF) UNG = MELVUN.VELCHE(1,NLCF) UTG = MELVUT.VELCHE(1,NLCF) UXG = UNG * CNX + UTG * CTX UYG = UNG * CNY + UTG * CTY PG = MELP.VELCHE(1,NLCF) GAMG = MELGAM.VELCHE(1,NLCF) C C**************** Jacobien par rapport aux variables primitives sur C la face C IF(IIMPL .EQ. 4)THEN C C******************* VLH C CALL VLHJ3(ROG,UXG,UYG,PG,GAMG,CNX,CNY,DFRUN) C ELSEIF(IIMPL .EQ. 5)THEN C C******************* AUSM+ C ROD=ROG PD=PG UND=-1.0D0*UNG UTD=UTG UXD = UND * CNX + UTD * CTX UYD = UND * CNY + UTD * CTY CALL AUSMW(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, & GAMG,CNX,CNY,DFRUN) ELSE CALL ERREUR(5) ENDIF DFRO(1)=0.0D0 DFRO(2)=0.0D0 DFRO(3)=0.0D0 DFRO(4)=0.0D0 DFRUT(1)=0.0D0 DFRUT(2)=0.0D0 DFRUT(3)=0.0D0 DFRUT(4)=0.0D0 DFRET(1)=0.0D0 DFRET(2)=0.0D0 DFRET(3)=0.0D0 DFRET(4)=0.0D0 C ELSEIF(NGC .EQ. NGCG)THEN C C************* On est sur G C CNX = MELVNX.VELCHE(1,NLCF) CNY = MELVNY.VELCHE(1,NLCF) CTX = MELT1X.VELCHE(1,NLCF) CTY = MELT1Y.VELCHE(1,NLCF) C ROG = MELRO.VELCHE(1,NLCF) UNG = MELVUN.VELCHE(1,NLCF) UTG = MELVUT.VELCHE(1,NLCF) UXG = UNG * CNX + UTG * CTX UYG = UNG * CNY + UTG * CTY PG = MELP.VELCHE(1,NLCF) GAMG = MELGAM.VELCHE(1,NLCF) ECIN=0.5D0*ROG*(UNG*UNG+UTG*UTG) RETG=PG/(GAMG-1.0D0)+ECIN C ROD = MELRO.VELCHE(3,NLCF) UND = MELVUN.VELCHE(3,NLCF) UTD = MELVUT.VELCHE(3,NLCF) UXD = UND * CNX + UTD * CTX UYD = UND * CNY + UTD * CTY PD = MELP.VELCHE(3,NLCF) C NLCD=MLENTC.LECT(NGCD) VOLD=MPVOLU.VPOCHA(NLCD,1) IF(IIMPL .EQ. 4)THEN C C******************* VLH C CALL VLHJ1(ROG,UXG,UYG,PG,RETG,GAMG,CNX,CNY,CTX,CTY $ ,DFRO,DFRUN,DFRUT,DFRET) ELSEIF(IIMPL .EQ. 5)THEN C C******************* AUSM+ C CALL AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, & GAMG,CNX,CNY,CTX,CTY, & DFRO,DFRUN,DFRUT,DFRET) ELSE CALL ERREUR(5) ENDIF ELSEIF(NGC .EQ. NGCD)THEN CNX = -1.0D0 * MELVNX.VELCHE(1,NLCF) CNY = -1.0D0 * MELVNY.VELCHE(1,NLCF) CTX = -1.0D0 * MELT1X.VELCHE(1,NLCF) CTY = -1.0D0 * MELT1Y.VELCHE(1,NLCF) C ROG = MELRO.VELCHE(3,NLCF) UNG = MELVUN.VELCHE(3,NLCF) UTG = MELVUT.VELCHE(3,NLCF) UXG = -1.0D0 * (UNG * CNX + UTG * CTX) UYG = -1.0D0 * (UNG * CNY + UTG * CTY) PG = MELP.VELCHE(3,NLCF) GAMG = MELGAM.VELCHE(3,NLCF) ECIN=0.5D0*ROG*(UNG*UNG+UTG*UTG) RETG=PG/(GAMG-1.0D0)+ECIN C ROD = MELRO.VELCHE(1,NLCF) UND = MELVUN.VELCHE(1,NLCF) UTD = MELVUT.VELCHE(1,NLCF) UXD = -1.0D0 * (UND * CNX + UTD * CTX) UYD = -1.0D0 * (UND * CNY + UTD * CTY) PD = MELP.VELCHE(1,NLCF) C NLCD=MLENTC.LECT(NGCG) VOLD=MPVOLU.VPOCHA(NLCD,1) IF(IIMPL .EQ. 4)THEN C C******************* VLH C CALL VLHJ1(ROG,UXG,UYG,PG,RETG,GAMG,CNX,CNY,CTX,CTY $ ,DFRO,DFRUN,DFRUT,DFRET) ELSEIF(IIMPL .EQ. 5)THEN C C******************* AUSM+ C CALL AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, & GAMG,CNX,CNY,CTX,CTY, & DFRO,DFRUN,DFRUT,DFRET) ELSE CALL ERREUR(5) ENDIF ELSE WRITE(IOIMP,*) 'FACEL = ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF C DO I1=1,4 DFRUX(I1)=DFRUN(I1)*CNX+DFRUT(I1)*CTX DFRUY(I1)=DFRUN(I1)*CNY+DFRUT(I1)*CTY ENDDO C NGNO=MELGRR.NUM(IFAC,IELPRI) NGN1=MELGRP.NUM(IFAC,IELPRI) NGN2=MELGRV.NUM(IFAC,IELPRI) C IF((NGN1 .NE. NGNO) .OR. (NGN2 .NE. NGNO))THEN WRITE(IOIMP,*) 'SPG COEFF GRAD =???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 ENDIF C NLNO=MLENTC.LECT(NGNO) IF(NLNO .EQ. 0)THEN C C**************** MURS C IPT1.NUM(IFAC,IELPRI)=NGC C Primal variable in IPT1.NUM(IFAC,IELPRI)=IPT1.NUM(NP,IELPRI) C Dual variable in IPT1.NUM(NP,IELPRI) C FUNCEL = -1.0D0 * SURF / VOLG UXR.AM(IELPRI,IFAC,NP)=UXR.AM(IELPRI,IFAC,NP) & +DFRUX(1)*FUNCEL UYR.AM(IELPRI,IFAC,NP)=UYR.AM(IELPRI,IFAC,NP) & +DFRUY(1)*FUNCEL UXUX.AM(IELPRI,IFAC,NP)=UXUX.AM(IELPRI,IFAC,NP) & +DFRUX(2)*FUNCEL UYUX.AM(IELPRI,IFAC,NP)=UYUX.AM(IELPRI,IFAC,NP) & +DFRUY(2)*FUNCEL UXUY.AM(IELPRI,IFAC,NP)=UXUY.AM(IELPRI,IFAC,NP) & +DFRUX(3)*FUNCEL UYUY.AM(IELPRI,IFAC,NP)=UYUY.AM(IELPRI,IFAC,NP) & +DFRUY(3)*FUNCEL UXRET.AM(IELPRI,IFAC,NP)=UXRET.AM(IELPRI,IFAC,NP) & +DFRUX(4)*FUNCEL UYRET.AM(IELPRI,IFAC,NP)=UYRET.AM(IELPRI,IFAC,NP) & +DFRUY(4)*FUNCEL ELSE C C Primal variable in IPT1.NUM(NP,IELPRI) C Dual variable in IPT1.NUM(NP,IELPRI) C FUNCEL=DFRO(1)*SURF RR.AM(IELPRI,NP,NP)=RR.AM(IELPRI,NP,NP)- & FUNCEL/VOLG C C Primal variable in NP C Dual variable in IFAC C RR.AM(IELPRI,NP,IFAC)=RR.AM(IELPRI,NP,IFAC)+FUNCEL $ /VOLD C FUNCEL=DFRO(2)*SURF RUX.AM(IELPRI,NP,NP)=RUX.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RUX.AM(IELPRI,NP,IFAC)=RUX.AM(IELPRI,NP,IFAC)+ FUNCEL $ /VOLD C FUNCEL=DFRO(3)*SURF RUY.AM(IELPRI,NP,NP)=RUY.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RUY.AM(IELPRI,NP,IFAC)=RUY.AM(IELPRI,NP,IFAC)+ FUNCEL $ /VOLD C FUNCEL=DFRO(4)*SURF RRET.AM(IELPRI,NP,NP)=RRET.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RRET.AM(IELPRI,NP,IFAC)= RRET.AM(IELPRI,NP,IFAC) $ +FUNCEL/VOLD C C FUNCEL=DFRUX(1)*SURF UXR.AM(IELPRI,NP,NP)=UXR.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UXR.AM(IELPRI,NP,IFAC)=UXR.AM(IELPRI,NP,IFAC)+FUNCEL $ /VOLD C FUNCEL=DFRUX(2)*SURF UXUX.AM(IELPRI,NP,NP)=UXUX.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UXUX.AM(IELPRI,NP,IFAC)=UXUX.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C FUNCEL=DFRUX(3)*SURF UXUY.AM(IELPRI,NP,NP)=UXUY.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UXUY.AM(IELPRI,NP,IFAC)=UXUY.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C FUNCEL=DFRUX(4)*SURF UXRET.AM(IELPRI,NP,NP)=UXRET.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UXRET.AM(IELPRI,NP,IFAC)=UXRET.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C C FUNCEL=DFRUY(1)*SURF UYR.AM(IELPRI,NP,NP)=UYR.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UYR.AM(IELPRI,NP,IFAC)=UYR.AM(IELPRI,NP,IFAC)+FUNCEL $ /VOLD C FUNCEL=DFRUY(2)*SURF UYUX.AM(IELPRI,NP,NP)=UYUX.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UYUX.AM(IELPRI,NP,IFAC)=UYUX.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C FUNCEL=DFRUY(3)*SURF UYUY.AM(IELPRI,NP,NP)=UYUY.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UYUY.AM(IELPRI,NP,IFAC)=UYUY.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C FUNCEL=DFRUY(4)*SURF UYRET.AM(IELPRI,NP,NP)=UYRET.AM(IELPRI,NP,NP) - & FUNCEL/VOLG UYRET.AM(IELPRI,NP,IFAC)=UYRET.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C C FUNCEL=DFRET(1)*SURF RETR.AM(IELPRI,NP,NP)=RETR.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RETR.AM(IELPRI,NP,IFAC)=RETR.AM(IELPRI,NP,IFAC)+FUNCEL $ /VOLD C FUNCEL=DFRET(2)*SURF RETUX.AM(IELPRI,NP,NP)=RETUX.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RETUX.AM(IELPRI,NP,IFAC)=RETUX.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C FUNCEL=DFRET(3)*SURF RETUY.AM(IELPRI,NP,NP)=RETUY.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RETUY.AM(IELPRI,NP,IFAC)=RETUY.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD C FUNCEL=DFRET(4)*SURF RETRET.AM(IELPRI,NP,NP)=RETRET.AM(IELPRI,NP,NP) - & FUNCEL/VOLG RETRET.AM(IELPRI,NP,IFAC)=RETRET.AM(IELPRI,NP,IFAC)+ $ FUNCEL/VOLD ENDIF C C************* Gradients computed with EULESCAL or EULEVECT C Computation of gradients contributions DO IFAC2 = 1, (NP-1), 1 C C**************** On controlle que les melemes de ELTFA et C le SPG du grad ont le meme ordre C NGFGR=MELEME.NUM(IFAC2,IELPRI) NLFGR=MLENTF.LECT(NGFGR) NG=MELEFE.NUM(1,NLFGR) ND=MELEFE.NUM(3,NLFGR) C IF(NG .EQ. ND)THEN C C******************* On est sur un mur: EULESCAL C DO IFACGR = 1, (NP -1) , 1 C C********************** Boucle sur les elts de MELGRR pour trouver la position C correspondante C NGGRA=MELGRR.NUM(IFACGR,IELPRI) IF(NGFGR .EQ. NGGRA) GOTO 10 ENDDO WRITE(IOIMP,*) 'ELTFA et SPG grad= ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 10 CONTINUE C C******************* Les coeff qui nous interesse sont donc dans la C position IFACGR du gradient C CELLR = ALR * (DXCF * MELCR1.VELCHE(IFACGR,IELPRI) $ +DYCF * MELCR2.VELCHE(IFACGR,IELPRI)) CELLP = ALP * (DXCF * MELCP1.VELCHE(IFACGR,IELPRI) $ + DYCF * MELCP2.VELCHE(IFACGR,IELPRI)) C C******************* La vitesse est un peux plus difficile C (HP EULEVECT ou vitesse imposé) C C CViVj contribution de vj sur le grad de vi C NLVIMP = MLEVL.LECT(NGGRA) IF(NLVIMP .NE.0)THEN C C*********************** NGGRA spg de conditions limites -> no C contribution C CVXVX=0.0D0 CVXVY=0.0D0 CVYVX=0.0D0 CVYVY=0.0D0 ELSE CNX1=MPNORM.VPOCHA(NLFGR,1) CNY1=MPNORM.VPOCHA(NLFGR,2) CTX1=-1.0D0*CNY1 CTY1=CNX1 FUNCEL=(DXCF * MELCV1.VELCHE(IFACGR,IELPRI)+DYCF $ * MELCV2.VELCHE(IFACGR,IELPRI)) CVXVX=FUNCEL*(CTX1*CTX1-CNX1*CNX1)*ALVX CVXVY=FUNCEL*(CTX1*CTY1-CNX1*CNY1)*ALVX CVYVX=FUNCEL*(CTX1*CTY1-CNX1*CNY1)*ALVY CVYVY=FUNCEL*(CTY1*CTY1-CNY1*CNY1)*ALVY ENDIF C ELSE DO IFACGR = 1, (NP -1) , 1 C C********************** Boucle sur les elts de MELGRR pour trouver la position C correspondante C NGGRA=MELGRR.NUM(IFACGR,IELPRI) IF((NGGRA .EQ. NG) .OR. (NGGRA .EQ. ND)) GOTO 20 ENDDO WRITE(IOIMP,*) 'ELTFA et SPG grad= ???' WRITE(IOIMP,*) 'SUBROUTINE konja5.eso' CALL ERREUR(5) GOTO 9999 20 CONTINUE CELLR = ALR * (DXCF * MELCR1.VELCHE(IFACGR,IELPRI) $ +DYCF * MELCR2.VELCHE(IFACGR,IELPRI)) CELLP = ALP * (DXCF * MELCP1.VELCHE(IFACGR,IELPRI) $ + DYCF * MELCP2.VELCHE(IFACGR,IELPRI)) CVXVX = ALVX * (DXCF * MELCV1.VELCHE(IFACGR,IELPRI) $ + DYCF * MELCV2.VELCHE(IFACGR,IELPRI)) CVXVY=0.0D0 CVYVY = ALVY * (DXCF * MELCV1.VELCHE(IFACGR,IELPRI) $ + DYCF * MELCV2.VELCHE(IFACGR,IELPRI)) CVYVX = 0.0D0 ENDIF C C**************** Dans le cas d'un mur (NLN0=0) , IFACGR = primale qui donne C contribution en IELPRI,IFACGR,IFAC OU (IELPRI,IFACGR,NP) C C Dans le cas non mur (NLN0=0) , IFACGR = primale qui donne C contribution en IELPRI,IFACGR,NP et en IELPRI,IFAC,NP C IF(NLNO .EQ. 0)THEN C C**************** MURS C FUNCEL = -1.0D0 * SURF / VOLG * CELLR UXR.AM(IELPRI,IFACGR,NP)=UXR.AM(IELPRI,IFACGR,NP) & +DFRUX(1)*FUNCEL UYR.AM(IELPRI,IFACGR,NP)=UYR.AM(IELPRI,IFACGR,NP) & +DFRUY(1)*FUNCEL C FUNCEL = -1.0D0 * SURF / VOLG UXUX.AM(IELPRI,IFACGR,NP)=UXUX.AM(IELPRI,IFACGR,NP) & +FUNCEL * (DFRUX(2)*CVXVX + DFRUX(3) *CVYVX) UYUX.AM(IELPRI,IFACGR,NP)=UYUX.AM(IELPRI,IFACGR,NP) & +FUNCEL * (DFRUY(2)*CVXVX + DFRUY(3) *CVYVX) C FUNCEL = -1.0D0 * SURF / VOLG UXUY.AM(IELPRI,IFACGR,NP)=UXUY.AM(IELPRI,IFACGR,NP) & +FUNCEL * (DFRUX(3)*CVYVY + DFRUX(2) *CVXVY) UYUY.AM(IELPRI,IFACGR,NP)=UYUY.AM(IELPRI,IFACGR,NP) & +FUNCEL * (DFRUY(3)*CVYVY + DFRUY(2) *CVXVY) C FUNCEL = -1.0D0 * SURF / VOLG * CELLP UXRET.AM(IELPRI,IFACGR,NP)=UXRET.AM(IELPRI,IFACGR $ ,NP)+DFRUX(4)*FUNCEL UYRET.AM(IELPRI,IFACGR,NP)=UYRET.AM(IELPRI,IFACGR $ ,NP)+DFRUY(4)*FUNCEL ELSE C C******************* dR/ C FUNCEL=DFRO(1)*SURF*CELLR RR.AM(IELPRI,IFACGR,NP)=RR.AM(IELPRI,IFACGR,NP)- & FUNCEL/VOLG RR.AM(IELPRI,IFACGR,IFAC)=RR.AM(IELPRI,IFACGR,IFAC) $ +FUNCEL/VOLD C FUNCEL=DFRO(4)*SURF*CELLP RRET.AM(IELPRI,IFACGR,NP)=RRET.AM(IELPRI,IFACGR,NP) $ -FUNCEL/VOLG RRET.AM(IELPRI,IFACGR,IFAC)=RRET.AM(IELPRI,IFACGR $ ,IFAC)+FUNCEL/VOLD C FUNCEL=SURF * (DFRO(2)*CVXVX + DFRO(3) *CVYVX) RUX.AM(IELPRI,IFACGR,NP)=RUX.AM(IELPRI,IFACGR,NP) & - FUNCEL / VOLG RUX.AM(IELPRI,IFACGR,IFAC)=RUX.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C FUNCEL=SURF * (DFRO(3)*CVYVY + DFRO(2) *CVXVY) RUY.AM(IELPRI,IFACGR,NP)=RUY.AM(IELPRI,IFACGR,NP) & - FUNCEL / VOLG RUY.AM(IELPRI,IFACGR,IFAC)=RUY.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C C******************* dUX/ C FUNCEL=DFRUX(1)*SURF*CELLR UXR.AM(IELPRI,IFACGR,NP)=UXR.AM(IELPRI,IFACGR,NP) - & FUNCEL/VOLG UXR.AM(IELPRI,IFACGR,IFAC)= UXR.AM(IELPRI,IFACGR $ ,IFAC)+FUNCEL/VOLD C FUNCEL=DFRUX(4)*SURF*CELLP UXRET.AM(IELPRI,IFACGR,NP)=UXRET.AM(IELPRI,IFACGR $ ,NP) -FUNCEL/VOLG UXRET.AM(IELPRI,IFACGR,IFAC)= UXRET.AM(IELPRI $ ,IFACGR,IFAC)+FUNCEL/VOLD C FUNCEL=SURF * (DFRUX(2)*CVXVX + DFRUX(3) *CVYVX) UXUX.AM(IELPRI,IFACGR,NP)=UXUX.AM(IELPRI,IFACGR,NP) & - FUNCEL / VOLG UXUX.AM(IELPRI,IFACGR,IFAC)=UXUX.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C FUNCEL=SURF * (DFRUX(3)*CVYVY + DFRUX(2) *CVXVY) UXUY.AM(IELPRI,IFACGR,NP)=UXUY.AM(IELPRI,IFACGR,NP) & - FUNCEL / VOLG UXUY.AM(IELPRI,IFACGR,IFAC)=UXUY.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C C******************* dUY/ C FUNCEL=DFRUY(1)*SURF*CELLR UYR.AM(IELPRI,IFACGR,NP)=UYR.AM(IELPRI,IFACGR,NP) - & FUNCEL/VOLG UYR.AM(IELPRI,IFACGR,IFAC)=UYR.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL/VOLD C FUNCEL=DFRUY(4)*SURF*CELLP UYRET.AM(IELPRI,IFACGR,NP)=UYRET.AM(IELPRI,IFACGR $ ,NP) -FUNCEL/VOLG UYRET.AM(IELPRI,IFACGR,IFAC)=UYRET.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL/VOLD C FUNCEL=SURF * (DFRUY(2)*CVXVX + DFRUY(3) *CVYVX) UYUX.AM(IELPRI,IFACGR,NP)=UYUX.AM(IELPRI,IFACGR,NP) & - FUNCEL / VOLG UYUX.AM(IELPRI,IFACGR,IFAC)=UYUX.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C FUNCEL=SURF * (DFRUY(3)*CVYVY + DFRUY(2) *CVXVY) UYUY.AM(IELPRI,IFACGR,NP)=UYUY.AM(IELPRI,IFACGR,NP) & - FUNCEL / VOLG UYUY.AM(IELPRI,IFACGR,IFAC)=UYUY.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C C******************* dRET/ C FUNCEL=DFRET(1)*SURF*CELLR RETR.AM(IELPRI,IFACGR,NP)=RETR.AM(IELPRI,IFACGR,NP) $ -FUNCEL/VOLG RETR.AM(IELPRI,IFACGR,IFAC)=RETR.AM(IELPRI,IFACGR $ ,IFAC)+FUNCEL/VOLD C FUNCEL=DFRET(4)*SURF*CELLP RETRET.AM(IELPRI,IFACGR,NP)=RETRET.AM(IELPRI,IFACGR $ ,NP)-FUNCEL/VOLG RETRET.AM(IELPRI,IFACGR,IFAC)=RETRET.AM(IELPRI $ ,IFACGR,IFAC)+FUNCEL/VOLD C FUNCEL=SURF * (DFRET(2)*CVXVX + DFRET(3) *CVYVX) RETUX.AM(IELPRI,IFACGR,NP)=RETUX.AM(IELPRI,IFACGR $ ,NP)- FUNCEL / VOLG RETUX.AM(IELPRI,IFACGR,IFAC)=RETUX.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C FUNCEL=SURF * (DFRET(3)*CVYVY + DFRET(2) *CVXVY) RETUY.AM(IELPRI,IFACGR,NP)=RETUY.AM(IELPRI,IFACGR $ ,NP)- FUNCEL / VOLG RETUY.AM(IELPRI,IFACGR,IFAC)=RETUY.AM(IELPRI,IFACGR $ ,IFAC) + FUNCEL / VOLD C ENDIF ENDDO ENDDO ENDDO C C C******* Transformation: C jacobien par rapport aux variables primitives -> C jacobien par rapport aux variables conservatives C Optimisation du programme: C on pourrais faire ça avant C NBEL=RR.AM(/1) NP=RR.AM(/2) MP=RR.AM(/3) C C******* Boucle sur les variables primales C DO IELETR=1,NBEL,1 DO IPRITR=1,NP,1 NGPRI=IPT1.NUM(IPRITR,IELETR) NLPRI=MLENTC.LECT(NGPRI) ROG=MPOVRN.VPOCHA(NLPRI,1) PG=MPOVPN.VPOCHA(NLPRI,1) GAMG=MPOGAM.VPOCHA(NLPRI,1) UXG=MPOVVN.VPOCHA(NLPRI,1) UYG=MPOVVN.VPOCHA(NLPRI,2) GAMGM1=GAMG-1.0D0 CELCO1=UXG/ROG CELCO2=UYG/ROG CELCO3=0.5D0*GAMGM1*(UXG*UXG+UYG*UYG) CELCO4=GAMGM1*UXG CELCO5=GAMGM1*UYG DO IDUATR=1,MP,1 DDRO=RR.AM(IELETR,IPRITR,IDUATR) DDUX=RUX.AM(IELETR,IPRITR,IDUATR) DDUY=RUY.AM(IELETR,IPRITR,IDUATR) DDP=RRET.AM(IELETR,IPRITR,IDUATR) RR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) & -(CELCO2*DDUY)+(CELCO3*DDP) RUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - & CELCO4*DDP RUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - & CELCO5*DDP RRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP C DDRO=UXR.AM(IELETR,IPRITR,IDUATR) DDUX=UXUX.AM(IELETR,IPRITR,IDUATR) DDUY=UXUY.AM(IELETR,IPRITR,IDUATR) DDP=UXRET.AM(IELETR,IPRITR,IDUATR) UXR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) & -(CELCO2*DDUY)+(CELCO3*DDP) UXUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - & CELCO4*DDP UXUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - & CELCO5*DDP UXRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP C DDRO=UYR.AM(IELETR,IPRITR,IDUATR) DDUX=UYUX.AM(IELETR,IPRITR,IDUATR) DDUY=UYUY.AM(IELETR,IPRITR,IDUATR) DDP=UYRET.AM(IELETR,IPRITR,IDUATR) UYR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) & -(CELCO2*DDUY)+(CELCO3*DDP) UYUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - & CELCO4*DDP UYUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - & CELCO5*DDP UYRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP C DDRO=RETR.AM(IELETR,IPRITR,IDUATR) DDUX=RETUX.AM(IELETR,IPRITR,IDUATR) DDUY=RETUY.AM(IELETR,IPRITR,IDUATR) DDP=RETRET.AM(IELETR,IPRITR,IDUATR) RETR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX) & -(CELCO2*DDUY)+(CELCO3*DDP) RETUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG - & CELCO4*DDP RETUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG - & CELCO5*DDP RETRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP ENDDO ENDDO ENDDO C MELEME = MLESPG.LECT(ISOUS) SEGDES MELEME MELEME = MLEPRI.LECT(ISOUS) SEGDES MELEME MELEME = MLEGRR.LECT(ISOUS) SEGDES MELEME MELEME = MLEGRP.LECT(ISOUS) SEGDES MELEME MELEME = MLEGRV.LECT(ISOUS) SEGDES MELEME C SEGDES RR , RUX , RUY , RRET , & UXR , UXUX , UXUY , UXRET , & UYR , UYUX , UYUY , UYRET , & RETR , RETUX , RETUY , RETRET SEGDES MELCR1, MELCR2 SEGDES MELCP1, MELCP2 SEGDES MELCV1, MELCV2 ENDDO C C C**** Les MELVALs C SEGDES MELRO SEGDES MELP SEGDES MELGAM SEGDES MELVUN SEGDES MELVUT SEGDES MELVNX SEGDES MELVNY SEGDES MELT1X SEGDES MELT1Y C C**** Les MPOVALs C C SEGDES MPOVRN SEGDES MPOVVN SEGDES MPOVPN SEGDES MPOGAM SEGDES MPVGRN SEGDES MPGLRN SEGDES MPVGPN SEGDES MPGLPN SEGDES MPVGPN SEGDES MPGLPN IF(MPVLIM .NE. 0) SEGDES MPVLIM C C**** Les MELEMEs, les MLENTIs, les volumes, les surfaces C SEGDES MELEMC SEGSUP MLENTC SEGSUP MLENTF SEGSUP MLEVL SEGDES MELEFE SEGDES MPVOLU SEGDES MPNORM SEGDES MPOVSU SEGDES MELPRI SEGSUP MLEPRI SEGSUP MLEGRR SEGSUP MLEGRP SEGSUP MLEGRV C C**** Coeff GRAD C MCHEL1=ICRN SEGDES MCHEL1 MCHEL1=ICPN SEGDES MCHEL1 MCHEL1=ICVN SEGDES MCHEL1 C SEGSUP MLESPG SEGDES MELTFA C C**** MATRIK, liste des inconnues C SEGDES MATRIK SEGDES IMATRI SEGDES MLMINC C C SEGDES RR , RUX , RUY , RRET , & UXR , UXUX , UXUY , UXRET , & UYR , UYUX , UYUY , UYRET , & RETR , RETUX , RETUY , RETRET SEGSUP MLENTC SEGSUP MLENTF SEGDES MLMINC SEGDES MLELIM 9999 CONTINUE RETURN END