Télécharger la plaquette Cast3M


fr en

Télécharger yla1bt.eso

Retour à la liste

Numérotation des lignes : 

yla1bt

  1. C YLA1BT    SOURCE    CB215821  20/11/25    13:43:59     10792          
  2. C YLAP1B    SOURCE    LEPOTIER  02/12/17    21:19:20     4510
  3.       SUBROUTINE YLA1BT(ICOGRT,MPTEMC,
  4.      $     MPVOLU,MPNORM,MPSURF,MELEFL,
  5.      $     KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP,
  6.      $     LMIXT,KRMIXT,IJACO,
  7.      $     IMPR,IRET)
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10. C***********************************************************************
  11. C NOM         : YLA1BT
  12. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  13. C               VF 2D.
  14. C               Ici, on ne calcule que les contributions à la matrice
  15. C               jacobienne faisant intervenir les coefficients pour le
  16. C               calcul des gradients de température (ICOGRT).
  17. C               (contributions à d Res_{\rho e_t} / d var
  18. C                var prenant successivement les valeurs :
  19. C                \rho,   \rho u,    \rho v,     \rho e_t)
  20. C
  21. C
  22. C
  23. C
  24. C LANGAGE     : ESOPE
  25. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  26. C               mél : gounand@semt2.smts.cea.fr
  27. C***********************************************************************
  28. C APPELES (UTIL)   : AJMTK : ajoute un objet de type MATSIM (non
  29. C                            standard) à un objet de type MATRIK.
  30. C APPELE PAR       : YLAP1A : Calcul de la matrice jacobienne du
  31. C                    résidu du laplacien VF 2D.
  32. C***********************************************************************
  33. C ENTREES            : ICOGRT (type MCHELM) : coefficients pour le
  34. C                        calcul du gradient de la température aux
  35. C                        interfaces.
  36. C                      MPROC (type MPOVAL) : masse volumique par
  37. C                                            élément.
  38. C                      MPVITC (type MPOVAL) : vitesse par élément.
  39. C                      MPTEMC (type MPOVAL) : température par élément.
  40. C                      MPVOLU (type MPOVAL) : volume des éléments.
  41. C                      MPNORM (type MPOVAL) : normale aux faces.
  42. C                      MPSURF (type MPOVAL) : surface des faces.
  43. C                      MELEFL (type MELEME) : connectivités face-(centre
  44. C                                             gauche, centre droit).
  45. C                      KRFACE (type MLENTI) : tableau de repérage dans
  46. C                        le maillage des faces des éléments.
  47. C                      KRCENT (type MLENTI) : tableau de repérage dans
  48. C                        le maillage des centres des éléments.
  49. C                      LCLIMT (type logique) : .TRUE. => CL de Dirichlet
  50. C                        sur la température.
  51. C                      KRTIMP (type MLENTI) : tableau de repérage dans
  52. C                        maillage des CL de Dirichlet sur la
  53. C                        température.
  54. C                      LCLIMQ (type logique) : .TRUE. => CL de Dirichlet
  55. C                        sur le flux de chaleur.
  56. C                      KRTIMP (type MLENTI) : tableau de repérage dans
  57. C                        maillage des CL de Dirichlet sur la
  58. C                        température.
  59. C                      KRQIMP (type MLENTI) : tableau de repérage dans
  60. C                        maillage des CL de Dirichlet sur le flux de
  61. C                        chaleur.
  62. C                      LMIXT (type logique) : .TRUE. => CL mixtes
  63. C                      KRMIXT (type MLENTI) : tableau de repérage dans
  64. C                        maillage des CL mixtes.
  65. C                      NOMINC (type MLMOTS) : noms des inconnues.
  66. C                      ICLAU               : option pour ne calculer
  67. C                                            que la thermique
  68. C                      KAPPA (type réel) : conductivité thermique (SI)
  69. C                      CV (type réel) : chaleur spécifique à volume
  70. C                                       constant (SI).
  71. C ENTREES/SORTIES    : IJACO (type MATRIK) : matrice jacobienne du
  72. C                        résidu du laplacien VF 2D.
  73. C SORTIES            : -
  74. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  75. C***********************************************************************
  76. C VERSION    : v1, 08/08/2001, version initiale
  77. C HISTORIQUE : v1, 08/08/2001, création
  78. C HISTORIQUE : 11/02/2003, ajout de l'option 'MIXT'
  79. C HISTORIQUE :
  80. C***********************************************************************
  81. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  82. C en cas de modification de ce sous-programme afin de faciliter
  83. C la maintenance !
  84. C***********************************************************************
  85.  
  86. -INC PPARAM
  87. -INC CCOPTIO
  88. -INC SMCHPOI
  89.       POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL
  90.       POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
  91. -INC SMCHAML
  92.       POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML
  93.       POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL
  94. -INC SMELEME
  95.       POINTEUR MELEFL.MELEME
  96.       POINTEUR MCOGRT.MELEME
  97. -INC SMLENTI
  98.       POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI,KRMIXT.MLENTI
  99.       POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
  100. -INC SMLMOTS
  101.       POINTEUR NOMINC.MLMOTS
  103.       POINTEUR IJACO.MATRIK
  104. *
  105. * Objet matrice élémentaire simplifié
  106. *
  107.       SEGMENT GMATSI
  108.       INTEGER POIPR1(NPP1,NEL1)
  109.       INTEGER POIDU1(1,NEL1)
  110.       INTEGER POIPR2(NPP2,NEL2)
  111.       INTEGER POIDU2(2,NEL2)
  112.       POINTEUR LMATSI(0).MATSIM
  113.       ENDSEGMENT
  114. *  Contributions de la part du gradient de température (CTGRT)
  115.       POINTEUR CTGRT.GMATSI
  116.       SEGMENT MATSIM
  117.       CHARACTER*8 NOMPRI,NOMDUA
  118.       REAL*8 VALMA1(1,NPP1,NEL1)
  119.       REAL*8 VALMA2(2,NPP2,NEL2)
  120.       ENDSEGMENT
  121.       POINTEUR RETRET.MATSIM
  122. *
  123.       REAL*8 KAPPA,CV
  124. *
  125.       INTEGER IMPR,IRET,ICLAU
  126. *
  127.       LOGICAL LCLIMT,LCLIMQ,LMIXT
  128.       LOGICAL LMUR
  129.       LOGICAL LCTRB1,LCTRB2
  130. *
  131.       INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
  132.       INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
  133.       INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
  134.       INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT
  135.       INTEGER NPTEL
  136. *
  137.       REAL*8 BETAX,BETAY,CNX,CNY
  138.       REAL*8 SIGNOR,SURFFA,VOLUEL
  139.       REAL*8 RHOP,UP,VP,TP
  140.       REAL*8 FACTOR
  141.       REAL*8 DTDRHO,DTDROU,DTDROV,DTDRET
  142. C
  143. *
  144. * Executable statements
  145. *
  146.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap1b.eso'
  147. *   On calcule la partie de d Res_{\rho e_t} / d var
  148. *   (var prend successivement les valeurs :
  149. *      \rho,   \rho u,    \rho v,     \rho e_t)
  150. *   faisant intervenir les coefficients pour le calcul des gradients de
  151. *   température (ICOGRT).
  152. *   C'est la partie contenant le terme : - \vect{q} \pscal \vect{n}
  153. *   Les noms de matrices élémentaires (type MATSIM) associées sont :
  154. *   RETRHO, RETROU, RETROV, RETRET
  155.       IF (LCLIMT) THEN
  156.          SEGACT KRTIMP
  157.       ENDIF
  158.       IF (LCLIMQ) THEN
  159.          SEGACT KRQIMP
  160.       ENDIF
  161.       IF (LMIXT) THEN
  162.          SEGACT KRMIXT
  163.       ENDIF
  164.       SEGACT KRCENT
  165.       SEGACT KRFACE
  166.       SEGACT MELEFL
  167.       SEGACT MPSURF
  168.       SEGACT MPNORM
  169.       SEGACT MPVOLU
  170.       SEGACT MPTEMC
  171.       SEGACT ICOGRT
  172.       NSOUCH=ICOGRT.IMACHE(/1)
  173.       DO 1 ISOUCH=1,NSOUCH
  174.          MCOGRT=ICOGRT.IMACHE(ISOUCH)
  175.          JCOGRT=ICOGRT.ICHAML(ISOUCH)
  176.          SEGACT JCOGRT
  177.          KCDTDX=JCOGRT.IELVAL(1)
  178.          SEGDES JCOGRT
  179.          SEGACT KCDTDX
  180.          SEGACT MCOGRT
  181.          NELEM=MCOGRT.NUM(/2)
  182.          NPTEL=MCOGRT.NUM(/1)
  183.          NPP1=NPTEL-1
  184.          NPP2=NPTEL-1
  185.          NEL1=NELEM
  186.          NEL2=NELEM
  187.          IEL1=1
  188.          IEL2=1
  189.  
  190.  
  191. c BOUCLE POUR EVALUER LE NOMBRE DE FACES FRONTIERES
  192. c POUR CHAQUE SOUS DOMAINE
  193.          NFRON = 1
  194.          DO 1200 IELEM=1,NELEM
  195.             NGFACE=MCOGRT.NUM(1,IELEM)
  196.             NLFACE=KRFACE.LECT(NGFACE)
  197.             IF (NLFACE.EQ.0) THEN
  198.                WRITE(IOIMP,*) 'Erreur de programmation n1'
  199.                GOTO 9999
  200.             ENDIF
  201.                NGCGAU=MELEFL.NUM(1,NLFACE)
  202.                NGCDRO=MELEFL.NUM(3,NLFACE)
  203.             IF ( NGCGAU.EQ.NGCDRO) THEN
  204.             NFRON = NFRON + 1
  205.             ENDIF
  206.  1200     CONTINUE
  207. c         WRITE(6,*) 'NFRON= ',NFRON
  208.          NEL1 = NFRON
  209. c         WRITE(6,*) 'NPP1= ',NPP1,'NPP2= ',NPP2
  210. c         WRITE(6,*) 'NEL1= ',NEL1,'NEL2= ',NEL2
  211.  
  212.          SEGINI RETRET
  213.          SEGINI CTGRT
  214.          RETRET.NOMPRI(1:4)='RETN'
  215.          RETRET.NOMPRI(5:8)='    '
  216.          RETRET.NOMDUA(1:4)='RETN'
  217.          RETRET.NOMDUA(5:8)='    '
  218.          DO 12 IELEM=1,NELEM
  219. *     Le premier point du support de ICOGRT est un point FACE
  220.             NGFACE=MCOGRT.NUM(1,IELEM)
  221.             NLFACE=KRFACE.LECT(NGFACE)
  222.             IF (NLFACE.EQ.0) THEN
  223.                WRITE(IOIMP,*) 'Erreur de programmation n1'
  224.                GOTO 9999
  225.             ENDIF
  226. *     On calcule la contribution à la matrice jacobienne IJACO de la face
  227. *     NGFAC (points duaux     : centres à gauche et à droite de la face)
  228. *           (points primaux   : une partie (bicoz conditions aux limites)
  229. *            de ceux du stencil pour le calcul du gradient
  230. *            à la face, ils doivent être des points centres)
  231. *       Si le flux de chaleur sur la face est imposé par les conditions
  232. *       aux limites, la contribution de la face à IJACO est nulle.
  233.             LCTRB1=.TRUE.
  234. c            IF (LCLIMQ) THEN
  235. c               NLCLQ=KRQIMP.LECT(NGFACE)
  236. c               IF (NLCLQ.NE.0) THEN
  237. c                  LCTRB1=.FALSE.
  238. c               ENDIF
  239. c            ENDIF
  240. *            WRITE(6,*) 'NGFACE= ',NGFACE,'LCTRB1= ', LCTRB1
  241.             IF (LCTRB1) THEN
  242.                NGCGAU=MELEFL.NUM(1,NLFACE)
  243.                NGCDRO=MELEFL.NUM(3,NLFACE)
  244.                LMUR=(NGCGAU.EQ.NGCDRO)
  245. *     On distingue le cas où la face est un bord du maillage (mur)
  246. *                  du cas où la face est interne au maillage
  247.                IF (.NOT.LMUR) THEN
  248.                   NPD=2
  249.                ELSE
  250.                   NPD=1
  251.                ENDIF
  252.                NPP=NPTEL-1
  253. *     IPD=1 : point à gauche du point NGFACE
  254. *     IPD=2 : point à droite du point NGFACE
  255.                DO 122 IPD=1,NPD
  256.                   NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
  257.                   IF (.NOT.LMUR) THEN
  258.                      CTGRT.POIDU2(IPD,IEL2)=NPDUAL
  259.                   ELSE
  260.                      CTGRT.POIDU1(IPD,IEL1)=NPDUAL
  261.                   ENDIF
  262.                   NLCEND=KRCENT.LECT(NPDUAL)
  263.                   IF (NLCEND.EQ.0) THEN
  264.                      WRITE(IOIMP,*) 'Erreur grave n1'
  265.                      GOTO 9999
  266.                   ENDIF
  267.                   DO 124 IPP=1,NPP
  268.                      NPPRIM=MCOGRT.NUM(IPP+1,IELEM)
  269.                      LCTRB2=.TRUE.
  270.                      IF (LCLIMT) THEN
  271.                         NLCLT=KRTIMP.LECT(NPPRIM)
  272.                         IF (NLCLT.NE.0) THEN
  273.                            LCTRB2=.FALSE.
  274.                         ENDIF
  275.                      ENDIF
  276.                      IF (LCLIMQ) THEN
  277.                         NLCLQ=KRQIMP.LECT(NPPRIM)
  278.                         IF (NLCLQ.NE.0) THEN
  279.                            LCTRB2=.FALSE.
  280.                         ENDIF
  281.                      ENDIF
  282.                      IF (LMIXT) THEN
  283.                         NLCLM=KRMIXT.LECT(NPPRIM)
  284.                         IF (NLCLM.NE.0) THEN
  285.                            LCTRB2=.FALSE.
  286.                         ENDIF
  287.                      ENDIF
  288.                      IF (.NOT.LCTRB2) THEN
  289. *     Lorsque une contribution est nulle, on fixe artificiellement le
  290. *     point primal égal au point dual.
  291.                         IF (.NOT.LMUR) THEN
  292.                            CTGRT.POIPR2(IPP,IEL2)=NPDUAL
  293.                            RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
  294.                         ELSE
  295.                            CTGRT.POIPR1(IPP,IEL1)=NPDUAL
  296.                            RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
  297.                         ENDIF
  298.                      ELSE
  299. *     Les contributions valent :
  300. *     (d Res_{\rho e_t})_d / (d var)_p =
  301. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \kappa
  302. *     * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
  303. *         ((dT)_p / (d var)_p)]
  304. *     avec :
  305. *          (dT)_p / (d \rho)_p = (1 / \rho_p) *
  306. *                ( (((u_p)^2 + (v_p)^2) / (2 * c_v)) - T )
  307. *          (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
  308. *          (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
  309. *          (dT)_p / (d \rho e_t)_p =   1 / (\rho_p * c_v)
  310. *     \beta_x : coefficients pour le calcul de dT/dx
  311. *     \beta_y : coefficients pour le calcul de dT/dy
  312. *
  313.                         NLCENP=KRCENT.LECT(NPPRIM)
  314.                         IF (NLCENP.EQ.0) THEN
  315.                            WRITE(IOIMP,*) 'Erreur grave n2'
  316.                            GOTO 9999
  317.                         ENDIF
  318. *       normale sortante pour IPD=1, rentrante pour IPD=2
  319.                         SIGNOR=(-1.D0)**(IPD+1)
  320.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  321.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  322.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  323.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  324.                         BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
  325.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA
  326.      $                       *(BETAX)
  327.                         DTDRET=1.D0
  328.                         IF (.NOT.LMUR) THEN
  329.                            CTGRT.POIPR2(IPP,IEL2)=NPPRIM
  330.                            RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
  331. c                        WRITE(6,*) 'IPD=',IPD,'IPP=',IPP,'IEL2=',IEL2
  332. c                        WRITE(6,*) 'FACT=',FACTOR*DTDRET
  333. c                        WRITE(6,*) 'SIGNOR=',SIGNOR,'VOLUEL=',VOLUEL,
  334. c     &                  'SURFFA=',SURFFA,'BETAX= ',BETAX
  335.                         ELSE
  336.                            CTGRT.POIPR1(IPP,IEL1)=NPPRIM
  337.                            RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
  338. c                        WRITE(6,*) 'IPD=',IPD,'IPP=',IPP,'IEL1=',IEL1
  339. c                        WRITE(6,*) 'FACT=',FACTOR*DTDRET
  340. c                        WRITE(6,*) 'SIGNOR=',SIGNOR,'VOLUEL=',VOLUEL,
  341. c     &                  'SURFFA=',SURFFA,'BETAX= ',BETAX
  342.                         ENDIF
  343.  
  344.                      ENDIF
  345.  124              CONTINUE
  346.  122           CONTINUE
  347.                IF (.NOT.LMUR) THEN
  348.                   IEL2=IEL2+1
  349.                ELSE
  350.                   IEL1=IEL1+1
  351.                ENDIF
  352.             ENDIF
  353.  12      CONTINUE
  354.          NPP1=NPTEL-1
  355.          NPP2=NPTEL-1
  356.          NEL1=IEL1-1
  357.          NEL2=IEL2-1
  358.  
  359. c         WRITE(6,*) 'APRES SEGADJ'
  360. c         WRITE(6,*) 'NPP1= ',NPP1,'NPP2= ',NPP2
  361. c         WRITE(6,*) 'NEL1= ',NEL1,'NEL2= ',NEL2
  362.          SEGDES MCOGRT
  363.          SEGDES KCDTDX
  364.  
  365.          SEGADJ RETRET
  366.          SEGADJ CTGRT
  367.          CTGRT.LMATSI(**)=RETRET
  368. * On accumule les matrices résultantes dans IJACO
  369.          CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
  370.          IF (IRET.NE.0) GOTO 9999
  371.          SEGSUP RETRET
  372.          SEGSUP CTGRT
  373. *
  374.  1    CONTINUE
  375.       SEGDES ICOGRT
  376.       SEGDES MPTEMC
  377.       SEGDES MPVOLU
  378.       SEGDES MPNORM
  379.       SEGDES MPSURF
  380.       SEGDES MELEFL
  381.       SEGDES KRFACE
  382.       SEGDES KRCENT
  383.       IF (LCLIMQ) THEN
  384.          SEGDES KRQIMP
  385.       ENDIF
  386.       IF (LCLIMT) THEN
  387.          SEGDES KRTIMP
  388.       ENDIF
  389.       IF (LMIXT) THEN
  390.          SEGDES KRMIXT
  391.       ENDIF
  392. *
  393. * Normal termination
  394. *
  395.       IRET=0
  396.       RETURN
  397. *
  398. * Format handling
  399. *
  400. *
  401. * Error handling
  402. *
  403.  9999 CONTINUE
  404.       IRET=1
  405.       WRITE(IOIMP,*) 'An error was detected in subroutine ylap1b'
  406.       RETURN
  407. *
  408. * End of subroutine YLAP1B
  409. *
  410.       END
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales