xlap1e
C XLAP1E SOURCE CB215821 20/11/25 13:43:16 10792 $ MPVOLU,MPNORM,MPSURF,MELEFL, $ KRFACE,KRCENT, $ LCLIMV,KRVIMP,MPVIMP, $ LCLITO,KRTOIM, $ NOMINC, $ MPMUC, $ IJACO, $ IMPR,IRET) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C NOM : XLAP1E C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien C VF 2D. C Ici, on ne calcule que des contributions 'compliquées' à C la matrice jacobienne faisant intervenir les C coefficients pour le calcul des gradients de vitesse C (ICOGRV) C (contributions à (d Res_{\rho e_t} / d var) C var prenant successivement les valeurs : C \rho, \rho u, \rho v, \rho e_t ) C Les contributions sont plus "compliquées" à calculer que C les simples (cf. xlap1c) car on a à dériver des produits C de fonctions de la vitesse. C xlap1d calcule aussi une partie des contributions C 'compliquées'. C C NOTE PERSO : On pourrait programmer ça plus lisiblement en stockant C les contributions dans un tableau de pointeurs (2 C indices, c'est possible ?) et en effectuant des produits C matriciels (coeff. x matrices de dérivées). C On n'y coupera pas si on veut regrouper 2D et 3D... C On ne va pas le faire. 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 : XLAP1A : Calcul de la matrice jacobienne du C résidu du laplacien VF 2D. C*********************************************************************** C ENTREES : ICOGRV (type MCHELM) : coefficients pour le C calcul du gradient de la vitesse aux C interfaces. C MPGRVF (type MPOVAL) : gradient de la vitesse C aux interfaces. C MPROC (type MPOVAL) : masse volumique par C élément. C MPVITC (type MPOVAL) : vitesse par é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 LCLIMV (type logique) : .TRUE. => CL de Dirichlet C sur la vitesse. C KRVIMP (type MLENTI) : tableau de repérage dans C maillage des CL de Dirichlet sur la C vitesse. C MPVIMP (type MPOVAL) : valeurs des CL de C Dirichlet sur la vitesse. C LCLITO (type logique) : .TRUE. => CL de Dirichlet C sur le tenseur des contraintes. C KRTOIM (type MLENTI) : tableau de repérage dans C maillage des CL de Dirichlet sur le tenseur des C contraintes C NOMINC (type MLMOTS) : noms des inconnues. C MPMUC (type MPOVAL) : viscosité dynamique (SI). C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du C résidu du laplacien VF 2D. C SORTIES : - C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 12/10/2001, version initiale C HISTORIQUE : v1, 12/10/2001, 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 MPGRVF.MPOVAL POINTEUR MPMUC.MPOVAL POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL POINTEUR MPVIMP.MPOVAL -INC SMCHAML POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL -INC SMELEME POINTEUR MELEFL.MELEME POINTEUR MCOGRV.MELEME -INC SMLENTI POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI POINTEUR KRCENT.MLENTI,KRFACE.MLENTI -INC SMLMOTS POINTEUR NOMINC.MLMOTS -INC SMMATRIK POINTEUR IJACO.MATRIK * * 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 compliquées 1 de la part du gradient de * vitesse (CCGRV1) POINTEUR CCGRV1.GMATSI SEGMENT MATSIM CHARACTER*8 NOMPRI,NOMDUA REAL*8 VALMA1(1,NPP1,NEL1) REAL*8 VALMA2(2,NPP2,NEL2) ENDSEGMENT POINTEUR RETRHO.MATSIM POINTEUR RETROU.MATSIM POINTEUR RETROV.MATSIM * REAL*8 MU * INTEGER IMPR,IRET * LOGICAL LCLIMV,LCLITO LOGICAL LMUR LOGICAL LCTR1A,LCTR1B * 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,NLCLTO,NLCLV INTEGER NLCED,NLCEG,NLFVI INTEGER NPTEL INTEGER ICOORX,NLCGAU,NLCDRO * REAL*8 DRD,DRG REAL*8 DUXXF,DUXYF,DUYXF,DUYYF REAL*8 UXD,UXF,UXG,UYD,UYF,UYG REAL*8 XD,XF,XG,XFMXD,XFMXG REAL*8 YD,YF,YG,YFMYD,YFMYG REAL*8 ALPHAX,ALPHAY,CNX,CNY REAL*8 SIGNOR,SURFFA,VOLUEL REAL*8 RHOP,UP,VP REAL*8 FACTOR,FACTDU,FACTDV REAL*8 DUDRHO,DUDROU,DVDRHO,DVDROV * * Executable statements * IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap1e.eso' * On calcule les contributions à (d Res_{\rho e_t} / d var) ; var * prenant successivement les valeurs : * \rho, \rho u, \rho v, \rho e_t * On dérive les termes : (\tens{\tau(\grad{u})} \prod \vect{u}) * \pscal \vect{n} * ce qui donne deux contributions. * Contribution 1 : * ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n} * Contribution 2 : * ( \tens{\tau} \prod (d \vect{u} / d var)) \pscal \vect{n} * Note : * pas de contribution à (d Res_{\rho e_t} / d \rho e_t). * Les noms de matrices élémentaires (type MATSIM) associées sont : * RETRHO, RETROU, RETROV IF (LCLIMV) THEN SEGACT KRVIMP SEGACT MPVIMP ENDIF IF (LCLITO) THEN SEGACT KRTOIM ENDIF SEGACT NOMINC SEGACT KRCENT SEGACT KRFACE SEGACT MELEFL SEGACT MPSURF SEGACT MPNORM SEGACT MPVOLU SEGACT MPMUC SEGACT MPROC SEGACT MPVITC SEGACT MPGRVF SEGACT ICOGRV NSOUCH=ICOGRV.IMACHE(/1) DO 1 ISOUCH=1,NSOUCH MCOGRV=ICOGRV.IMACHE(ISOUCH) JCOGRV=ICOGRV.ICHAML(ISOUCH) SEGACT JCOGRV KDUNDX=JCOGRV.IELVAL(1) KDUNDY=JCOGRV.IELVAL(2) SEGDES JCOGRV SEGACT KDUNDX SEGACT KDUNDY SEGACT MCOGRV NELEM=MCOGRV.NUM(/2) NPTEL=MCOGRV.NUM(/1) NPP1=NPTEL-1 NPP2=NPTEL-1 NEL1=NELEM NEL2=NELEM IEL1=1 IEL2=1 SEGINI RETRHO SEGINI RETROU SEGINI RETROV SEGINI CCGRV1 RETRHO.NOMPRI(5:8)=' ' RETRHO.NOMDUA(5:8)=' ' RETROU.NOMPRI(5:8)=' ' RETROU.NOMDUA(5:8)=' ' RETROV.NOMPRI(5:8)=' ' RETROV.NOMDUA(5:8)=' ' DO 12 IELEM=1,NELEM * Le premier point du support de ICOGRV est un point FACE NGFACE=MCOGRV.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 1 à la matrice jacobienne IJACO de la * face NGFACE * (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 tenseur des contraintes sur la face est imposé par les * conditions aux limites, la contribution 1 de la face à IJACO est * nulle. LCTR1A=.TRUE. IF (LCLITO) THEN NLCLTO=KRTOIM.LECT(NGFACE) IF (NLCLTO.NE.0) THEN LCTR1A=.FALSE. ENDIF ENDIF IF (LCTR1A) THEN 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 ICOORX = ((IDIM + 1) * (NGFACE - 1))+1 XF = MCOORD.XCOOR(ICOORX) YF = MCOORD.XCOOR(ICOORX+1) ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1 XG = MCOORD.XCOOR(ICOORX) YG = MCOORD.XCOOR(ICOORX+1) XFMXG = XF - XG YFMYG = YF - YG DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)) ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1 XD = MCOORD.XCOOR(ICOORX) YD = MCOORD.XCOOR(ICOORX+1) XFMXD = XF - XD YFMYD = YF - YD DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)) ELSE ALPHA=0.0D0 UMALPH=1.0D0 ENDIF MU=UMALPH*MPMUC.VPOCHA(NLCGAU,1) + * On calcule tout d'abord une interpolation de la vitesse sur la * face NGFACE de la même manière que dans xlap12.eso. IF (LCLIMV) THEN NLFVI=KRVIMP.LECT(NGFACE) ELSE NLFVI=0 ENDIF * La vitesse est imposée sur la face, rien à calculer IF (NLFVI.NE.0) THEN UXF=MPVIMP.VPOCHA(NLFVI,1) UYF=MPVIMP.VPOCHA(NLFVI,2) ELSE NLCEG=KRCENT.LECT(NGCGAU) NLCED=KRCENT.LECT(NGCDRO) * Cas non au bord IF (.NOT.LMUR) THEN * Paramètres géométriques ICOORX = ((IDIM + 1) * (NGFACE - 1))+1 XF = MCOORD.XCOOR(ICOORX) YF = MCOORD.XCOOR(ICOORX+1) ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1 XG = MCOORD.XCOOR(ICOORX) YG = MCOORD.XCOOR(ICOORX+1) XFMXG = XF - XG YFMYG = YF - YG DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)) ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1 XD = MCOORD.XCOOR(ICOORX) YD = MCOORD.XCOOR(ICOORX+1) XFMXD = XF - XD YFMYD = YF - YD DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)) DUXXF = MPGRVF.VPOCHA(NLFACE,1) DUXYF = MPGRVF.VPOCHA(NLFACE,2) DUYXF = MPGRVF.VPOCHA(NLFACE,3) DUYYF = MPGRVF.VPOCHA(NLFACE,4) * Calcul de la vitesse UXG = MPVITC.VPOCHA(NLCEG,1) UYG = MPVITC.VPOCHA(NLCEG,2) UXD = MPVITC.VPOCHA(NLCED,1) UYD = MPVITC.VPOCHA(NLCED,2) * Correction de la vitesse lineaire exacte UXF = UXF + UYF = UYF + * Cas au bord ELSE * Parametres geometriques ICOORX = ((IDIM + 1) * (NGFACE - 1))+1 XF = MCOORD.XCOOR(ICOORX) YF = MCOORD.XCOOR(ICOORX+1) ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1 XG = MCOORD.XCOOR(ICOORX) YG = MCOORD.XCOOR(ICOORX+1) XFMXG = XF - XG YFMYG = YF - YG DUXXF = MPGRVF.VPOCHA(NLFACE,1) DUXYF = MPGRVF.VPOCHA(NLFACE,2) DUYXF = MPGRVF.VPOCHA(NLFACE,3) DUYYF = MPGRVF.VPOCHA(NLFACE,4) * Calcul de la vitesse UXF = MPVITC.VPOCHA(NLCEG,1) UYF = MPVITC.VPOCHA(NLCEG,2) * Correction de la vitesse lineaire exacte UXF = UXF + (DUXXF * XFMXG ) + (DUXYF * YFMYG ) UYF = UYF + (DUYXF * XFMXG ) + (DUYYF * YFMYG ) ENDIF ENDIF * 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 ELSE NPD=1 ENDIF 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 CCGRV1.POIDU2(IPD,IEL2)=NPDUAL ELSE CCGRV1.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=MCOGRV.NUM(IPP+1,IELEM) LCTR1B=.TRUE. IF (LCLIMV) THEN NLCLV=KRVIMP.LECT(NPPRIM) IF (NLCLV.NE.0) THEN LCTR1B=.FALSE. ENDIF ENDIF IF (.NOT.LCTR1B) THEN * Lorsque une contribution est nulle, on fixe artificiellement le * point primal égal au point dual. IF (.NOT.LMUR) THEN CCGRV1.POIPR2(IPP,IEL2)=NPDUAL RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0 RETROU.VALMA2(IPD,IPP,IEL2)=0.D0 RETROV.VALMA2(IPD,IPP,IEL2)=0.D0 ELSE CCGRV1.POIPR1(IPP,IEL1)=NPDUAL RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0 RETROU.VALMA1(IPD,IPP,IEL1)=0.D0 RETROV.VALMA1(IPD,IPP,IEL1)=0.D0 ENDIF ELSE * Les contributions 1 valent : * (d Res_{\rho e_t})_d / (d var)_p = * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \mu * * [ [ [ (n_x * (( 4/3 * u_f * \alpha_x) + (v_f * \alpha_y))) * + (n_y * ((u_f * \alpha_y) + (-2/3 * v_f * \alpha_x))) * ] * * ((du)_p / (d var)_p) * ] * + [ [ (n_x * ((-2/3 * u_f * \alpha_y) + (v_f * \alpha_x))) * + (n_y * ((u_f * \alpha_x) + ( 4/3 * v_f * \alpha_y))) * ] * * ((dv)_p / (d var)_p) * ] * ] 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) ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM) ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM) RHOP =MPROC.VPOCHA(NLCENP,1) UP =MPVITC.VPOCHA(NLCENP,1) VP =MPVITC.VPOCHA(NLCENP,2) FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*MU FACTDU=(CNX*((( 4.D0/3.D0)*UXF*ALPHAX) $ +(UYF*ALPHAY))) $ + (CNY*((UXF*ALPHAY) $ +((-2.D0/3.D0)*UYF*ALPHAX))) FACTDV=(CNX*(((-2.D0/3.D0)*UXF*ALPHAY) $ +(UYF*ALPHAX))) $ + (CNY*((UXF*ALPHAX) $ +(( 4.D0/3.D0)*UYF*ALPHAY))) DUDRHO=-UP /RHOP DUDROU=1.D0/RHOP DVDRHO=-VP /RHOP DVDROV=1.D0/RHOP IF (.NOT.LMUR) THEN CCGRV1.POIPR2(IPP,IEL2)=NPPRIM RETRHO.VALMA2(IPD,IPP,IEL2)= $ FACTOR*( (FACTDU*DUDRHO) $ +(FACTDV*DVDRHO)) RETROU.VALMA2(IPD,IPP,IEL2)= $ FACTOR* (FACTDU*DUDROU) RETROV.VALMA2(IPD,IPP,IEL2)= $ FACTOR* (FACTDV*DVDROV) ELSE CCGRV1.POIPR1(IPP,IEL1)=NPPRIM RETRHO.VALMA1(IPD,IPP,IEL1)= $ FACTOR*( (FACTDU*DUDRHO) $ +(FACTDV*DVDRHO)) RETROU.VALMA1(IPD,IPP,IEL1)= $ FACTOR* (FACTDU*DUDROU) RETROV.VALMA1(IPD,IPP,IEL1)= $ FACTOR* (FACTDV*DVDROV) ENDIF ENDIF 124 CONTINUE 122 CONTINUE IF (.NOT.LMUR) THEN IEL2=IEL2+1 ELSE IEL1=IEL1+1 ENDIF ENDIF 12 CONTINUE NPP1=NPTEL-1 NPP2=NPTEL-1 NEL1=IEL1-1 NEL2=IEL2-1 SEGADJ RETRHO SEGADJ RETROU SEGADJ RETROV SEGADJ CCGRV1 CCGRV1.LMATSI(**)=RETRHO CCGRV1.LMATSI(**)=RETROU CCGRV1.LMATSI(**)=RETROV * On accumule les matrices résultantes dans IJACO IF (IRET.NE.0) GOTO 9999 SEGSUP RETRHO SEGSUP RETROU SEGSUP RETROV SEGSUP CCGRV1 * SEGDES MCOGRV SEGDES KDUNDY SEGDES KDUNDX 1 CONTINUE SEGDES ICOGRV SEGDES MPGRVF SEGDES MPVITC SEGDES MPMUC SEGDES MPROC SEGDES MPVOLU SEGDES MPNORM SEGDES MPSURF SEGDES MELEFL SEGDES KRFACE SEGDES KRCENT SEGDES NOMINC IF (LCLITO) THEN SEGDES KRTOIM ENDIF IF (LCLIMV) THEN SEGDES KRVIMP SEGDES MPVIMP ENDIF * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine xlap1e' RETURN * * End of subroutine XLAP1E * END
© Cast3M 2003 - Tous droits réservés.
Mentions légales