yfpu
C YFPU SOURCE CB215821 23/01/25 21:15:40 11573 SUBROUTINE YFPU IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C SYNTAXE : C C FPU NU UET YP <,VPAROI> C C 1------2 C (R1,AL1) LEF FLUIDE NOEUDS 1 2 C C C ANU VISCOSITE CINEMATIQUE C UET U* C YP DISTANCE A LA PAROI C VPAROI VITESSE DE LA PAROI (PAR DEFAUT 0.) C C CAS TRIDIMENSIONNEL C 4 ________ 3 C / FLUIDE / C 1 /________/2 C C C*********************************************************************** -INC CCVQUA4 -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMCOORD -INC SIZFFB POINTEUR IZF1.IZFFM,IZH2.IZHR -INC SMLREEL POINTEUR MLREE4.MLREEL -INC SMLENTI POINTEUR IPADS.MLENTI,IPADT.MLENTI -INC SMMATRIK -INC SMELEME POINTEUR MELEMS.MELEME,MELEKE.MELEME,MELEMX.MELEME,MLEMST.MELEME -INC SMCHPOI POINTEUR IZD.MPOVAL,IZDS.MPOVAL,IZDK.MPOVAL,IZDE.MPOVAL POINTEUR IZG1.MCHPOI, IZGG1.MPOVAL POINTEUR IZG2.MCHPOI, IZGG2.MPOVAL POINTEUR IZG3.MCHPOI, IZGG3.MPOVAL POINTEUR IZGV.MCHPOI, IZGGV.MPOVAL POINTEUR IZTU1.MPOVAL POINTEUR MUE.MCHPOI,MZUE.MPOVAL POINTEUR MYP.MCHPOI,MZYP.MPOVAL POINTEUR MZNU.MPOVAL POINTEUR MZRO.MPOVAL,MZMU.MPOVAL POINTEUR IZVOL.MPOVAL POINTEUR MNT.MCHPOI,MZNT.MPOVAL -INC SMLMOTS POINTEUR LINCO.MLMOTS CHARACTER*8 NOMZ,NOMI,TYPE,TYPC,NOM0 CHARACTER*(LOCOMP) NOM4(5) LOGICAL LOGI PARAMETER (NTB=1) CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB),IXV(3),AJ(270),XYZ(90) DATA LTAB/'KIZX '/ C***************************************************************************** CFPU segact mcoord *pv ca le fait mieux en initialisant ikel *pv si on ne le fait pas et que par hasard il vaut 1 *pv les resultats ne sont pas ceux attendus ikel=0 c write(6,*)' debut FPU ' izdk=0 izde=0 IF (IERR.NE.0) RETURN MTABX=KTAB(1) C C- Récupération de la table EQEX (pointeur MTAB1) C IF(MTAB1.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' EQEX ' MOTERR( 9:16) = ' EQEX ' MOTERR(17:24) = ' KIZX ' RETURN ENDIF IF(NASTOK.EQ.0)THEN RETURN ENDIF C C- Récupération de la table INCO (pointeur KINC) C IF(KINC.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' INCO ' MOTERR( 9:16) = ' INCO ' MOTERR(17:24) = ' EQEX ' RETURN ENDIF C***************************************************************************** C OPTIONS C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG IKOMP=0 IAXI=0 IF(IFOMOD.EQ.0)IAXI=2 C C- Récupération de la table des options KOPT (pointeur KOPTI) C IF (KOPTI.EQ.0) THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' KOPT ' MOTERR( 9:16) = ' KOPT ' MOTERR(17:24) = ' KIZX ' RETURN ENDIF C write(6,*)' Avant les options ' TYPE=' ' TYPE=' ' IF(KFORM.NE.0.AND.KFORM.NE.1)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = 'EF/EFM1 ' RETURN ENDIF C***************************************************************************** C C- Récupération de la table DOMAINE associée au domaine local C TYPE=' ' IF(TYPE.NE.'MMODEL')THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' DOMZ ' MOTERR( 9:16) = ' DOMZ ' MOTERR(17:24) = ' KIZX ' RETURN ENDIF C E/ MMODEL : Pointeur de la table contenant l'information cherchée C /S IPOINT : Pointeur sur la table DOMAINE C /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE C INEFMD=4 LINB IF (IERR.NE.0) RETURN SEGACT MELEME IF (IERR.NE.0) RETURN C*** TYPE='LISTMOTS' SEGACT LINCO IKOMP=0 IF (KFORM.EQ.0.AND.IARG.NE.3)THEN write(6,*)' FPU : nb d arguments incorrect :',IARG GO TO 90 ENDIF IF (KFORM.EQ.1.AND.IARG.NE.5.AND.IARG.NE.6)THEN write(6,*)' FPU : nb d arguments incorrect :',IARG GO TO 90 ENDIF IKRS=0 IKR=1 N=1 NC=1 SEGINI MZRO MZRO.VPOCHA(1,1)=1.D0 C--Cas incompressible IF(IARG.EQ.3)THEN IF(IKOMP.EQ.0)THEN C 1er coefficient : nu IXV(1)=0 IXV(2)=1 IXV(3)=0 IRET=0 & MTABX,KINC,1,IXV,MNU,MZNU,NPT1,NC1,IKM,IRET) IF(IRET.EQ.0)RETURN C 2ème coefficient : uet IXV(1)=MELEMC IF(KFORM.EQ.1)IXV(1)=MELEMS IXV(2)=0 IXV(3)=0 IRET=1 & MTABX,KINC,2,IXV,MUE,MZUE,NPT2,NC2,IK2,IRET) IF(IRET.EQ.0)RETURN SEGACT MZUE*MOD C 3ème coefficient : yp IXV(1)=0 IXV(2)=1 IXV(3)=0 IRET=0 & MTABX,KINC,3,IXV,MYP,MZYP,NPT3,NC3,IK3,IRET) IF(IRET.EQ.0)RETURN ENDIF C--Cas compressible ELSEIF(IKOMP.EQ.1)THEN C 1er coefficient : mu IXV(1)=MELEMC IXV(2)=1 IXV(3)=0 IRET=0 & MTABX,KINC,1,IXV,MNU,MZNU,NPT1,NC1,IKM,IRET) IF(IRET.EQ.0)RETURN C 2ème coefficient : uet IXV(1)=MELEMC IXV(2)=0 IXV(3)=0 IRET=0 & MTABX,KINC,2,IXV,MUE,MZUE,NPT2,NC2,IK2,IRET) IF(IRET.EQ.0)RETURN SEGACT MZUE*MOD ENDIF NPTS=MELEMS.NUM(/2) C***************************************************************************** C VERIFICATIONS SUR LES INCONNUES WRITE(6,*)' On attend 1 ou 3' C Indice %m1:8 : nombre d'arguments incorrect MOTERR(1:8) = 'INCO ' RETURN ENDIF C --> 1 ere Inconnue DO 15 I=1,IDIM WRITE(NOM4(I),FMT='(I1,A3)')I,NOMI(1:3) 15 CONTINUE TYPE=' ' IF(TYPE.NE.'CHPOINT ')THEN WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table' RETURN ELSE ENDIF C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ C--Cas Complet: ro un mu uet yp IF((IARG.EQ.5.OR.IARG.EQ.6).AND.KFORM.EQ.1)THEN C 1er coefficient : ro IKRS=1 IXV(1)=IGEOM0 IXV(2)=1 IXV(3)=0 IRET=0 & MTABX,KINC,1,IXV,MRO,MZRO,NPT1,NC1,IKR,IRET) IF(IRET.EQ.0)RETURN C 2ème coefficient : un IXV(1)=-IGEOM0 IXV(2)=0 IXV(3)=0 IRET=1 & MTABX,KINC,2,IXV,MUE,MZUE,NPT2,NC2,IKU,IRET) IF(IRET.EQ.0)RETURN SEGACT MZUE*MOD C 3ème coefficient : mu IXV(1)=IGEOM0 IXV(2)=1 IXV(3)=0 IRET=0 & MTABX,KINC,3,IXV,MMU,MZMU,NPT3,NC3,IKM,IRET) IF(IRET.EQ.0)RETURN C 4ème coefficient : uet IXV(1)=MELEMC IF(KFORM.EQ.1)IXV(1)=MELEMS IXV(2)=0 IXV(3)=0 IRET=1 & MTABX,KINC,4,IXV,MUE,MZUE,NPT2,NC2,IK2,IRET) IF(IRET.EQ.0)RETURN SEGACT MZUE*MOD TYPE=' ' c####################################################################### c################# IKEL = 1 => Syntaxe II ############################## IF(TYPE.EQ.'MMODEL')THEN IKEL=1 TYPE=' ' IF(TYPE.NE.'TABLE')THEN KPREFPU=0 write(6,*)'Operateur FPU : On Preconditionne ' ELSE KPREFPU=1 ENDIF C 6ème coefficient : NUT IXV(1)=IGEOM0 IXV(2)=0 IXV(3)=0 IRET=0 & MTABX,KINC,6,IXV,MNT,MZNT,NPT3,NC3,IK3,IRET) IF(IRET.EQ.0)RETURN C On calcule le Gradient de U * IGEOM0 Support géomtrique domaine complet * MDTOT modèle pour le maillage complet * MTABT table domaine en résultant * MTABZ table domaine de la zone frontière C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c on extrait les élements s'appuyant largement sur la frontière c -> MELEMX IF(KPREFPU.EQ.0)THEN SEGACT MUE,MZUE NSOUPO=MUE.IPCHP(/1) IF(NSOUPO.GT.1)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = 'NSOUPO>1' RETURN ENDIF NAT=MUE.JATTRI(/1) N =MZUE.VPOCHA(/1) NC=MZUE.VPOCHA(/2) SEGINI MYP,MSOUPO,MZYP DO 651 I=1,NAT MYP.JATTRI(I)=MUE.JATTRI(I) 651 CONTINUE SEGINI MSOUPO MYP.IPCHP(1)=MSOUPO MSOUP1=MUE.IPCHP(1) SEGACT MSOUP1 DO 652 I=1,NC NOCOMP(I)=MSOUP1.NOCOMP(I) NOHARM(I)=MSOUP1.NOHARM(I) 652 CONTINUE IGEOC =MSOUP1.IGEOC IPOVAL=MZYP CALL PRCHAN CALL PREXTR CALL NBNO ELSE type=' ' SEGACT MYP MSOUPO=MYP.IPCHP(1) SEGACT MSOUPO MZYP=MSOUPO.IPOVAL SEGACT MZYP type=' ' ENDIF c -> MELEMX c on extrait les elements s'appuyant largement sur la frontiere C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. c write(6,*)' FPU : NBSOUS=',NBSOUS c CALL ECROBJ('MAILLAGE',MELEMT) c CALL PRTRAC c CALL LIROBJ('MAILLAGE',MELEMT,1,IRET) c CALL ECROBJ('MAILLAGE',MELEMX) c CALL ECROBJ('MAILLAGE',MELEMS) c CALL PRFUSE c CALL PRTRAC c CALL LIROBJ('MAILLAGE',MELEMX,1,IRET) SEGACT MLEMST NPTT=MLEMST.NUM(/2) SEGACT MELEMS c nls=MELEMS.NUM(/2) c write(6,*)'MELEMS=' c write(6,1001)(MELEMS.NUM(1,ii),ii=1,nls) c write(6,*)' I P A D S ' c npds=ipads.lect(/1) c write(6,1001)(ipads.lect(ii),ii=1,npds) NC=IZTU1.VPOCHA(/2) NPTI=IZTU1.VPOCHA(/1) SEGACT MELEMX NBSOUS=MELEMX.LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 NUTOEL=0 IES=IDIM SEGACT,MCOORD DO 19 L=1,NBSOUS IPT1=MELEMX IF(NBSOUS.NE.1)IPT1=MELEMX.LISOUS(L) SEGACT IPT1 NP =IPT1.NUM(/1) IF(KPREFPU.EQ.0)THEN SEGINI MLREE1,MLREE3 SEGINI MLREE2 ELSE TYPE=' ' TYPE=' ' TYPE=' ' SEGACT MLREE1,MLREE2,MLREE3 ENDIF N=IZD.VPOCHA(/1) NC=IZD.VPOCHA(/2) SEGINI IZDS IZGG2=IZDS c write(6,*)' On met les conditions limites sur K et epsilon' TYPE='SOMMET' NC=2 c CALL KRCHPT(TYPE,MELEMX,NC,IZG2,NOM4(4)) ENDIF & IZTU1.VPOCHA,NPTI,NC,IPADS.LECT,NPTS, & MZMU.VPOCHA,IKM,MZUE.VPOCHA,MZYP.VPOCHA,IPT1.ITYPEL, 19 CONTINUE SEGDES,MCOORD c CALL ECROBJ('MAILLAGE',MELEMX) c CALL PRTRAC c CALL LIROBJ('MAILLAGE',MELEMX,1,IRET) NC=IZTU1.VPOCHA(/2) NPTI=IZTU1.VPOCHA(/1) TYPE='SOMMET' C --> Vitesse seulement on l'impose a 0. TYPE='SOMMET' NC=IDIM TYPE=' ' c write(6,*)' Y F P U :: KCLIM = ',KCLIM IF(KCLIM.EQ.0)THEN ELSE c On remet à 0 les conditions limites précédentes MCHPO4=KCLIM SEGACT MCHPO4 NSP1=MCHPO4.IPCHP(/1) DO 882 L=1,NSP1 MSOUPO=MCHPO4.IPCHP(L) SEGACT MSOUPO NC=NOCOMP(/2) MELEKE=IGEOC MPOVA4=IPOVAL SEGACT MELEKE,MPOVA4*MOD NBPKE=MELEKE.NUM(/2) DO 883 N=1,NC IF(NOCOMP(N).NE.NOM4(1).AND.NOCOMP(N).NE.NOM4(2) & .AND.NOCOMP(N).NE.NOM4(3))GO TO 883 DO 884 I=1,NBPKE NKE=IPADS.LECT(MELEKE.NUM(1,I)) IF(NKE.NE.0)MPOVA4.VPOCHA(I,N)=0.D0 884 CONTINUE 883 CONTINUE SEGDES MSOUPO,MPOVA4,MELEKE 882 CONTINUE CALL PRFUSE ENDIF C --> FIN Vitesse seulement on l'impose a 0. IZGG2=IZGG1 C --> On impose K et Epsilon c On remet à 0 la condition limite précédente pour NOM4(4) et NOM4(5) TYPE=' ' MCHPO4=KCLIM SEGACT MCHPO4 NSP1=MCHPO4.IPCHP(/1) DO 982 L=1,NSP1 MSOUPO=MCHPO4.IPCHP(L) SEGACT MSOUPO NC=NOCOMP(/2) MELEKE=IGEOC MPOVA4=IPOVAL SEGACT MELEKE,MPOVA4*MOD NBPKE=MELEKE.NUM(/2) DO 983 N=1,NC IF(NOCOMP(N).NE.NOM4(4).AND.NOCOMP(N).NE.NOM4(5))GO TO 983 DO 984 I=1,NBPKE NKE=IPADS.LECT(MELEKE.NUM(1,I)) IF(NKE.NE.0)MPOVA4.VPOCHA(I,N)=0.D0 984 CONTINUE 983 CONTINUE SEGDES MSOUPO,MPOVA4,MELEKE 982 CONTINUE CALL PRFUSE c write(6,*)' On remet les CL sur KN et EN CLIM=',MCHPOI c CALL ECROBJ('CHPOINT',MCHPOI) c CALL PRLIST ENDIF SEGSUP IZDS IF(IKRS.EQ.0)SEGSUP MZRO SEGDES IZTU1 SEGDES IZG1,IZGG1 SEGDES IZVOL SEGDES LINCO SEGSUP MLENTI,IPADS c On crée un MATRIK vide NRIGE=7 NKID =9 NKMT =7 NMATRI=0 SEGINI MATRIK c write(6,*)' Fin FPU Syntaxe II' RETURN c################# IKEL = 1 => Syntaxe II ## FIN ####################### c####################################################################### ELSE IKEL=0 C 5ème coefficient : yp IXV(1)=0 IXV(2)=1 IXV(3)=0 IRET=0 & MTABX,KINC,5,IXV,MYP,MZYP,NPT3,NC3,IK3,IRET) IF(IRET.EQ.0)RETURN ENDIF ENDIF C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ NC=IZTU1.VPOCHA(/2) NPTI=IZTU1.VPOCHA(/1) TYPE='SOMMET' IZGG2=IZGG1 C --> 2 eme Inconnue TYPE='SOMMET' NC=2 TYPE=' ' IF(KCLIM.NE.0)THEN MCHPO4=KCLIM SEGACT MCHPO4 NSP1=MCHPO4.IPCHP(/1) DO 782 L=1,NSP1 MSOUPO=MCHPO4.IPCHP(L) SEGACT MSOUPO NC=NOCOMP(/2) MELEKE=IGEOC MPOVA4=IPOVAL SEGACT MELEKE,MPOVA4*MOD NBPKE=MELEKE.NUM(/2) DO 783 N=1,NC IF(NOCOMP(N).NE.NOM4(4).AND.NOCOMP(N).NE.NOM4(5))GO TO 783 DO 784 I=1,NBPKE NKE=IPADS.LECT(MELEKE.NUM(1,I)) IF(NKE.NE.0)MPOVA4.VPOCHA(I,N)=0.D0 784 CONTINUE 783 CONTINUE SEGDES MSOUPO,MPOVA4,MELEKE 782 CONTINUE ENDIF SEGDES MCHPO4 ENDIF N=IZD.VPOCHA(/1) NC=IZD.VPOCHA(/2) SEGINI IZDS SEGINI IZDE,IZDK ENDIF SEGACT MELEME NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 NUTOEL=0 IES=IDIM DO 1 L=1,NBSOUS IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(L) SEGACT IPT1 NOM0=NOMS(IPT1.ITYPEL)//' ' SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD IZH2=KZHR(2) IZF1=KTP(1) NES=GR(/1) NPG=GR(/3) NP =IPT1.NUM(/1) C SUBROUTINE XCVFPU(NEL,K0,NP,IES,IAXI,IPADL, C & LEF,XYZ, ----> IPT1,COOR C & VOLF, ----> IZVOL.T, C & UN,TK,TE, ----> IZTU1.T,IZTU2.T,IZTU3.T, C & F, ----> IZG1, C & DK,DE, ----> IZD2,IZD3 C & ANU,IKC,UET,YP, ----> IZTG1.T,IKM,IZTG2.T,IZTG3.T, C & VPAROI,IKV, IZTG4.T,IK4, C & PORO,NPR,IPOR) ----> IZPORO,NPOR,IOP7 if(izde.eq.0) izde=izds if(izdk.eq.0) izdk=izds SEGACT,MCOORD IF(KFORM.EQ.0)THEN & NPTS,IPADS.LECT, & IPT1.NUM,XCOOR, & IZVOL.VPOCHA, & IZTU1.VPOCHA,IZGG2.VPOCHA(1,1),IZGG2.VPOCHA(1,2), & IZGG1.VPOCHA, & MZNU.VPOCHA,IKM,MZUE.VPOCHA,MZYP.VPOCHA) SEGACT,MCOORD ELSEIF(KFORM.EQ.1.AND.IKEL.EQ.0)THEN & NPTI,LECT,IZTU1.VPOCHA,NPTS,IPADS.LECT,MZUE.VPOCHA, & IZGG1.VPOCHA,MZMU.VPOCHA,IKM,MZRO.VPOCHA,IKR, ENDIF SEGDES,MCOORD SEGSUP IZFFM,IZF1,IZHR,IZH2 1 CONTINUE IF(KFORM.EQ.1)THEN IF(IKEL.EQ.0)THEN DO 124 I=1,NPTS MZUE.VPOCHA(I,1)=IZDS.VPOCHA(I,1)/IZD.VPOCHA(I,1) 124 CONTINUE ENDIF DO 125 I=1,NPTS IZGG2.VPOCHA(I,1)=IZDK.VPOCHA(I,1)/IZD.VPOCHA(I,1) IZGG2.VPOCHA(I,2)=IZDE.VPOCHA(I,1)/IZD.VPOCHA(I,1) 125 CONTINUE ENDIF ENDIF SEGSUP IZDS SEGSUP IZDE,IZDK ENDIF IF(IKRS.EQ.0)SEGSUP MZRO SEGDES IZTU1 SEGDES IZG1,IZGG1 SEGDES IZVOL SEGDES LINCO SEGSUP MLENTI,IPADS NRIGE=7 NKID =9 NKMT =7 NMATRI=0 SEGINI MATRIK IF(KCLIM.EQ.0)THEN ELSE CALL PRFUSE ENDIF ENDIF c write(6,*)' Fin FPU' RETURN 90 CONTINUE WRITE(6,*)' Interuption anormale de FPU ' RETURN 1001 FORMAT(20(1X,I5)) 1002 FORMAT(10(1X,1PE11.4)) END
© Cast3M 2003 - Tous droits réservés.
Mentions légales