C ZLAP2B SOURCE CB215821 20/11/25 13:45:09 10792 SUBROUTINE ZLAP2B(PROPHY,PROPH2,MPROC, $ MPVOLU,MPNORM,MPSURF,MELEFL, $ KRFACE,KRCENT,LCLIMR,KRRIMP,MPRIMP, $ NOMINC, $ IJACO, $ IMPR,IRET) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C NOM : ZLAP2B C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien C VF 3D (termes multi-espèces). C Ici, on ne calcule que les contributions à la matrice C jacobienne faisant intervenir les coefficients pour le C calcul des gradients de Yk C (contributions à d Res_{\rho Yk} / d var C var prenant successivement les valeurs : C \rho, \rho Yk) C C C C C LANGAGE : ESOPE C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non C standard) à un objet de type MATRIK. C APPELE PAR : ZLAP2A : Calcul de la matrice jacobienne du C résidu du laplacien VF 3D. C*********************************************************************** C ENTREES : PROPHY (type PROPHY) : propriétés des espèces C PROPH2 (type PROPH2) : précond. de PROPHY C MPROC (type MPOVAL) : masse volumique par C élément. C MPVOLU (type MPOVAL) : volume des éléments. C MPNORM (type MPOVAL) : normale aux faces. C MPSURF (type MPOVAL) : surface des faces. C MELEFL (type MELEME) : connectivités face-(centre C gauche, centre droit). C KRFACE (type MLENTI) : tableau de repérage dans C le maillage des faces des éléments. C KRCENT (type MLENTI) : tableau de repérage dans C le maillage des centres des éléments. C LCLIMR (type logique) : .TRUE. => CL de Dirichlet C sur la densité. C KRRIMP (type MLENTI) : tableau de repérage dans C maillage des CL de Dirichlet sur la densité. C MPRIMP (type MPOVAL) : valeurs des CL de C Dirichlet sur la densité. C NOMINC (type MLMOTS) : noms des inconnues. C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du C résidu du laplacien VF 3D. C SORTIES : - C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 08/03/2002, version initiale C HISTORIQUE : v1, 08/03/2002, création C HISTORIQUE : C HISTORIQUE : C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMCHPOI POINTEUR MPDKC.MPOVAL POINTEUR MPROC.MPOVAL,MPYKC.MPOVAL POINTEUR MPRIMP.MPOVAL POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL -INC SMCHAML POINTEUR ICOGRY.MCHELM,JCOGRY.MCHAML POINTEUR KDYKDX.MELVAL,KDYKDY.MELVAL,KDYKDZ.MELVAL -INC SMELEME POINTEUR MELEFL.MELEME POINTEUR MCOGRY.MELEME -INC SMLENTI POINTEUR KRYIMP.MLENTI,KRRIMP.MLENTI POINTEUR KRCENT.MLENTI,KRFACE.MLENTI -INC SMLMOTS POINTEUR NOMINC.MLMOTS -INC SMMATRIK POINTEUR IJACO.MATRIK INTEGER NESP,IESP SEGMENT PROPHY CHARACTER*4 NOMESP(NESP+1) REAL*8 CV(NESP+1) REAL*8 R(NESP+1) REAL*8 H0K(NESP+1) POINTEUR CDIFF(NESP+1).MCHPOI POINTEUR YK(NESP+1).MCHPOI POINTEUR GRADYK(NESP+1).MCHPOI POINTEUR CGRYK(NESP+1).MCHELM POINTEUR CLYK(NESP+1).MCHPOI ENDSEGMENT SEGMENT PROPH2 POINTEUR MPDIFF(NESP+1).MPOVAL POINTEUR MPVALY(NESP+1).MPOVAL POINTEUR MPGRAD(NESP+1).MPOVAL LOGICAL LCLIM(NESP+1) POINTEUR KRCLIM(NESP+1).MLENTI ENDSEGMENT * * Objet matrice élémentaire simplifié * SEGMENT GMATSI INTEGER POIPR1(NPP1,NEL1) INTEGER POIDU1(1,NEL1) INTEGER POIPR2(NPP2,NEL2) INTEGER POIDU2(2,NEL2) POINTEUR LMATSI(0).MATSIM ENDSEGMENT * Contributions de la part du gradient de YK (CTGRY) POINTEUR CTGRY.GMATSI SEGMENT MATSIM CHARACTER*8 NOMPRI,NOMDUA REAL*8 VALMA1(1,NPP1,NEL1) REAL*8 VALMA2(2,NPP2,NEL2) ENDSEGMENT POINTEUR RYKRHO.MATSIM POINTEUR RYKRYK.MATSIM * REAL*8 DKG,DKD,DKF,RHOG,RHOD,RHOF * INTEGER IMPR,IRET * LOGICAL LCLIMY,LCLIMR LOGICAL LMUR LOGICAL LCTRB2 * INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2 INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2 INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL INTEGER NLCENP,NLCEND,NLFACE,NLCLY,NLFRI INTEGER NPTEL * REAL*8 BETAX,BETAY,BETAZ,CNX,CNY,CNZ REAL*8 SIGNOR,SURFFA,VOLUEL REAL*8 RHOP,YKP REAL*8 FACTOR REAL*8 DYKDRO,DYKDRY * INTEGER ICOORX,NLCGAU,NLCDRO REAL*8 XF,YF,ZF,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG & ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH C * * Executable statements * IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans zlap2b.eso' SEGACT PROPHY NESP=PROPHY.CV(/1)-1 SEGACT PROPH2 DO IESP=1,NESP IF (LCLIMR) THEN SEGACT KRRIMP SEGACT MPRIMP ENDIF LCLIMY=PROPH2.LCLIM(IESP) IF (LCLIMY) THEN KRYIMP=PROPH2.KRCLIM(IESP) SEGACT KRYIMP ENDIF SEGACT NOMINC SEGACT KRCENT SEGACT KRFACE SEGACT MELEFL SEGACT MPSURF SEGACT MPNORM SEGACT MPVOLU MPDKC=PROPH2.MPDIFF(IESP) SEGACT MPDKC SEGACT MPROC MPYKC=PROPH2.MPVALY(IESP) SEGACT MPYKC ICOGRY=PROPHY.CGRYK(IESP) SEGACT ICOGRY NSOUCH=ICOGRY.IMACHE(/1) DO 1 ISOUCH=1,NSOUCH MCOGRY=ICOGRY.IMACHE(ISOUCH) JCOGRY=ICOGRY.ICHAML(ISOUCH) SEGACT JCOGRY KDYKDX=JCOGRY.IELVAL(1) KDYKDY=JCOGRY.IELVAL(2) KDYKDZ=JCOGRY.IELVAL(3) SEGDES JCOGRY SEGACT KDYKDX SEGACT KDYKDY SEGACT KDYKDZ SEGACT MCOGRY NELEM=MCOGRY.NUM(/2) NPTEL=MCOGRY.NUM(/1) NPP1=NPTEL-1 NPP2=NPTEL-1 NEL1=NELEM NEL2=NELEM IEL1=1 IEL2=1 SEGINI RYKRHO SEGINI RYKRYK SEGINI CTGRY RYKRHO.NOMPRI(1:4)=NOMINC.MOTS(1) RYKRHO.NOMPRI(5:8)=' ' RYKRHO.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP) RYKRHO.NOMDUA(5:8)=' ' RYKRYK.NOMPRI(1:4)=NOMINC.MOTS(IDIM+2+IESP) RYKRYK.NOMPRI(5:8)=' ' RYKRYK.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP) RYKRYK.NOMDUA(5:8)=' ' DO 12 IELEM=1,NELEM * Le premier point du support de ICOGRY est un point FACE NGFACE=MCOGRY.NUM(1,IELEM) NLFACE=KRFACE.LECT(NGFACE) IF (NLFACE.EQ.0) THEN WRITE(IOIMP,*) 'Erreur de programmation n°1' GOTO 9999 ENDIF * On calcule la contribution à la matrice jacobienne IJACO de la face * NGFAC (points duaux : centres à gauche et à droite de la face) * (points primaux : une partie (bicoz conditions aux limites) * de ceux du stencil pour le calcul du gradient * à la face, ils doivent être des points centres) * Si le flux de chaleur sur la face est imposé par les conditions * aux limites, la contribution de la face à IJACO est nulle. NGCGAU=MELEFL.NUM(1,NLFACE) NGCDRO=MELEFL.NUM(3,NLFACE) NLCGAU=KRCENT.LECT(NGCGAU) NLCDRO=KRCENT.LECT(NGCDRO) LMUR=(NGCGAU.EQ.NGCDRO) * On distingue le cas où la face est un bord du maillage (mur) * du cas où la face est interne au maillage IF (.NOT.LMUR) THEN NPD=2 ICOORX = ((IDIM + 1) * (NGFACE - 1))+1 XF = MCOORD.XCOOR(ICOORX) YF = MCOORD.XCOOR(ICOORX+1) ZF = MCOORD.XCOOR(ICOORX+2) ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1 XG = MCOORD.XCOOR(ICOORX) YG = MCOORD.XCOOR(ICOORX+1) ZG = MCOORD.XCOOR(ICOORX+2) XFMXG = XF - XG YFMYG = YF - YG ZFMZG = ZF - ZG DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG)) ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1 XD = MCOORD.XCOOR(ICOORX) YD = MCOORD.XCOOR(ICOORX+1) ZD = MCOORD.XCOOR(ICOORX+2) XFMXD = XF - XD YFMYD = YF - YD ZFMZD = ZF - ZD DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD)) ALPHA=DRG/(DRG+DRD) UMALPH= 1.0D0 - ALPHA ELSE NPD=1 ALPHA=0.0D0 UMALPH=1.0D0 ENDIF IF (LCLIMR) THEN NLFRI=KRRIMP.LECT(NGFACE) ELSE NLFRI=0 ENDIF IF (NLFRI.GT.0) THEN RHOF=MPRIMP.VPOCHA(NLFRI,1) ELSE RHOG = MPROC.VPOCHA(NLCGAU,1) RHOD = MPROC.VPOCHA(NLCDRO,1) RHOF = UMALPH * RHOG + ALPHA * RHOD ENDIF DKG = MPDKC.VPOCHA(NLCGAU,1) DKD = MPDKC.VPOCHA(NLCDRO,1) DKF = UMALPH * DKG + ALPHA * DKD NPP=NPTEL-1 * IPD=1 : point à gauche du point NGFACE * IPD=2 : point à droite du point NGFACE DO 122 IPD=1,NPD NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE) IF (.NOT.LMUR) THEN CTGRY.POIDU2(IPD,IEL2)=NPDUAL ELSE CTGRY.POIDU1(IPD,IEL1)=NPDUAL ENDIF NLCEND=KRCENT.LECT(NPDUAL) IF (NLCEND.EQ.0) THEN WRITE(IOIMP,*) 'Erreur grave n°1' GOTO 9999 ENDIF DO 124 IPP=1,NPP NPPRIM=MCOGRY.NUM(IPP+1,IELEM) LCTRB2=.TRUE. IF (LCLIMY) THEN NLCLY=KRYIMP.LECT(NPPRIM) IF (NLCLY.NE.0) THEN LCTRB2=.FALSE. ENDIF ENDIF IF (.NOT.LCTRB2) THEN * Lorsque une contribution est nulle, on fixe artificiellement le * point primal égal au point dual. IF (.NOT.LMUR) THEN CTGRY.POIPR2(IPP,IEL2)=NPDUAL RYKRHO.VALMA2(IPD,IPP,IEL2)=0.D0 RYKRYK.VALMA2(IPD,IPP,IEL2)=0.D0 ELSE CTGRY.POIPR1(IPP,IEL1)=NPDUAL RYKRHO.VALMA1(IPD,IPP,IEL1)=0.D0 RYKRYK.VALMA1(IPD,IPP,IEL1)=0.D0 ENDIF ELSE * Les contributions valent : * (d Res_{\rho Yk})_d / (d var)_p = * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \rho_f * Dk_f * * [ ((n_x * \beta_x) +(n_y * \beta_y) +(n_z * \beta_z)) * * ((dYk)_p / (d var)_p)] * avec : * (dYk)_p / (d \rho)_p = - Yk_p / \rho_p * (dYk)_p / (d \rho Yk)_p = 1 / \rho_p * \beta_x : coefficients pour le calcul de dYk/dx * \beta_y : coefficients pour le calcul de dYk/dy * \beta_z : coefficients pour le calcul de dYk/dz * NLCENP=KRCENT.LECT(NPPRIM) IF (NLCENP.EQ.0) THEN WRITE(IOIMP,*) 'Erreur grave n°2' GOTO 9999 ENDIF * normale sortante pour IPD=1, rentrante pour IPD=2 SIGNOR=(-1.D0)**(IPD+1) VOLUEL=MPVOLU.VPOCHA(NLCEND,1) SURFFA=MPSURF.VPOCHA(NLFACE,1) CNX =MPNORM.VPOCHA(NLFACE,1) CNY =MPNORM.VPOCHA(NLFACE,2) CNZ =MPNORM.VPOCHA(NLFACE,3) BETAX =KDYKDX.VELCHE(IPP+1,IELEM) BETAY =KDYKDY.VELCHE(IPP+1,IELEM) BETAZ =KDYKDZ.VELCHE(IPP+1,IELEM) RHOP =MPROC.VPOCHA(NLCENP,1) YKP =MPYKC.VPOCHA(NLCENP,1) * YKP =MPYKC.VPOCHA(NLCEND,1) FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*RHOF*DKF $ *((CNX*BETAX)+(CNY*BETAY)+(CNZ*BETAZ)) DYKDRO=-YKP/RHOP DYKDRY=1.D0/RHOP IF (.NOT.LMUR) THEN CTGRY.POIPR2(IPP,IEL2)=NPPRIM RYKRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRO RYKRYK.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRY ELSE CTGRY.POIPR1(IPP,IEL1)=NPPRIM RYKRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRO RYKRYK.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRY ENDIF ENDIF 124 CONTINUE 122 CONTINUE IF (.NOT.LMUR) THEN IEL2=IEL2+1 ELSE IEL1=IEL1+1 ENDIF 12 CONTINUE NPP1=NPTEL-1 NPP2=NPTEL-1 NEL1=IEL1-1 NEL2=IEL2-1 SEGADJ RYKRHO SEGADJ RYKRYK SEGADJ CTGRY CTGRY.LMATSI(**)=RYKRHO CTGRY.LMATSI(**)=RYKRYK * On accumule les matrices résultantes dans IJACO CALL AJMTK(CTGRY,IJACO,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGSUP RYKRHO SEGSUP RYKRYK SEGSUP CTGRY * SEGDES MCOGRY SEGDES KDYKDZ SEGDES KDYKDY SEGDES KDYKDX 1 CONTINUE SEGDES ICOGRY SEGDES MPYKC SEGDES MPROC SEGDES MPDKC SEGDES MPVOLU SEGDES MPNORM SEGDES MPSURF SEGDES MELEFL SEGDES KRFACE SEGDES KRCENT SEGDES NOMINC IF (LCLIMY) THEN SEGDES KRYIMP ENDIF IF (LCLIMR) THEN SEGDES KRRIMP SEGDES MPRIMP ENDIF ENDDO SEGDES PROPH2 SEGDES PROPHY * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine zlap2b' RETURN * * End of subroutine ZLAP2B * END