Télécharger la plaquette Cast3M


fr en

Télécharger zlap2b.eso

Retour à la liste

Numérotation des lignes : 

zlap2b

  1. C ZLAP2B    SOURCE    CB215821  20/11/25    13:45:09     10792          
  2.       SUBROUTINE ZLAP2B(PROPHY,PROPH2,MPROC,
  3.      $     MPVOLU,MPNORM,MPSURF,MELEFL,
  4.      $     KRFACE,KRCENT,LCLIMR,KRRIMP,MPRIMP,
  5.      $     NOMINC,
  6.      $     IJACO,
  7.      $     IMPR,IRET)
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10. C***********************************************************************
  11. C NOM         : ZLAP2B
  12. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  13. C               VF 3D (termes multi-espèces).
  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 Yk
  17. C               (contributions à d Res_{\rho Yk} / d var
  18. C                var prenant successivement les valeurs :
  19. C                \rho,   \rho Yk)
  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       : ZLAP2A : Calcul de la matrice jacobienne du
  31. C                    résidu du laplacien VF 3D.
  32. C***********************************************************************
  33. C ENTREES            : PROPHY (type PROPHY) : propriétés des espèces
  34. C                      PROPH2 (type PROPH2) : précond. de PROPHY
  35. C                      MPROC (type MPOVAL) : masse volumique par
  36. C                                            élément.
  37. C                      MPVOLU (type MPOVAL) : volume des éléments.
  38. C                      MPNORM (type MPOVAL) : normale aux faces.
  39. C                      MPSURF (type MPOVAL) : surface des faces.
  40. C                      MELEFL (type MELEME) : connectivités face-(centre
  41. C                                             gauche, centre droit).
  42. C                      KRFACE (type MLENTI) : tableau de repérage dans
  43. C                        le maillage des faces des éléments.
  44. C                      KRCENT (type MLENTI) : tableau de repérage dans
  45. C                        le maillage des centres des éléments.
  46. C                      LCLIMR (type logique) : .TRUE. => CL de Dirichlet
  47. C                        sur la densité.
  48. C                      KRRIMP (type MLENTI) : tableau de repérage dans
  49. C                        maillage des CL de Dirichlet sur la densité.
  50. C                      MPRIMP (type MPOVAL) : valeurs des CL de
  51. C                        Dirichlet sur la densité.
  52. C                      NOMINC (type MLMOTS) : noms des inconnues.
  53. C ENTREES/SORTIES    : IJACO (type MATRIK) : matrice jacobienne du
  54. C                        résidu du laplacien VF 3D.
  55. C SORTIES            : -
  56. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  57. C***********************************************************************
  58. C VERSION    : v1, 08/03/2002, version initiale
  59. C HISTORIQUE : v1, 08/03/2002, création
  60. C HISTORIQUE :
  61. C HISTORIQUE :
  62. C***********************************************************************
  63. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  64. C en cas de modification de ce sous-programme afin de faciliter
  65. C la maintenance !
  66. C***********************************************************************
  67.  
  68. -INC PPARAM
  69. -INC CCOPTIO
  70. -INC SMCOORD
  71. -INC SMCHPOI
  72.       POINTEUR MPDKC.MPOVAL
  73.       POINTEUR MPROC.MPOVAL,MPYKC.MPOVAL
  74.       POINTEUR MPRIMP.MPOVAL
  75.       POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
  76. -INC SMCHAML
  77.       POINTEUR ICOGRY.MCHELM,JCOGRY.MCHAML
  78.       POINTEUR KDYKDX.MELVAL,KDYKDY.MELVAL,KDYKDZ.MELVAL
  79. -INC SMELEME
  80.       POINTEUR MELEFL.MELEME
  81.       POINTEUR MCOGRY.MELEME
  82. -INC SMLENTI
  83.       POINTEUR KRYIMP.MLENTI,KRRIMP.MLENTI
  84.       POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
  85. -INC SMLMOTS
  86.       POINTEUR NOMINC.MLMOTS
  88.       POINTEUR IJACO.MATRIK
  89.       INTEGER NESP,IESP
  90.       SEGMENT PROPHY
  91.          CHARACTER*4 NOMESP(NESP+1)
  92.          REAL*8 CV(NESP+1)
  93.          REAL*8 R(NESP+1)
  94.          REAL*8 H0K(NESP+1)
  95.          POINTEUR CDIFF(NESP+1).MCHPOI
  96.          POINTEUR YK(NESP+1).MCHPOI
  97.          POINTEUR GRADYK(NESP+1).MCHPOI
  98.          POINTEUR CGRYK(NESP+1).MCHELM
  99.          POINTEUR CLYK(NESP+1).MCHPOI
  100.       ENDSEGMENT
  101.       SEGMENT PROPH2
  102.          POINTEUR MPDIFF(NESP+1).MPOVAL
  103.          POINTEUR MPVALY(NESP+1).MPOVAL
  104.          POINTEUR MPGRAD(NESP+1).MPOVAL
  105.          LOGICAL LCLIM(NESP+1)
  106.          POINTEUR KRCLIM(NESP+1).MLENTI
  107.       ENDSEGMENT
  108. *
  109. * Objet matrice élémentaire simplifié
  110. *
  111.       SEGMENT GMATSI
  112.       INTEGER POIPR1(NPP1,NEL1)
  113.       INTEGER POIDU1(1,NEL1)
  114.       INTEGER POIPR2(NPP2,NEL2)
  115.       INTEGER POIDU2(2,NEL2)
  116.       POINTEUR LMATSI(0).MATSIM
  117.       ENDSEGMENT
  118. *  Contributions de la part du gradient de YK (CTGRY)
  119.       POINTEUR CTGRY.GMATSI
  120.       SEGMENT MATSIM
  121.       CHARACTER*8 NOMPRI,NOMDUA
  122.       REAL*8 VALMA1(1,NPP1,NEL1)
  123.       REAL*8 VALMA2(2,NPP2,NEL2)
  124.       ENDSEGMENT
  125.       POINTEUR RYKRHO.MATSIM
  126.       POINTEUR RYKRYK.MATSIM
  127. *
  128.       REAL*8 DKG,DKD,DKF,RHOG,RHOD,RHOF
  129. *
  130.       INTEGER IMPR,IRET
  131. *
  132.       LOGICAL LCLIMY,LCLIMR
  133.       LOGICAL LMUR
  134.       LOGICAL LCTRB2
  135. *
  136.       INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
  137.       INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
  138.       INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
  139.       INTEGER NLCENP,NLCEND,NLFACE,NLCLY,NLFRI
  140.       INTEGER NPTEL
  141. *
  142.       REAL*8 BETAX,BETAY,BETAZ,CNX,CNY,CNZ
  143.       REAL*8 SIGNOR,SURFFA,VOLUEL
  144.       REAL*8 RHOP,YKP
  145.       REAL*8 FACTOR
  146.       REAL*8 DYKDRO,DYKDRY
  147. *
  148.       INTEGER ICOORX,NLCGAU,NLCDRO
  149.       REAL*8 XF,YF,ZF,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG
  150.      &     ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH
  151. C
  152. *
  153. * Executable statements
  154. *
  155.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans zlap2b.eso'
  156.       SEGACT PROPHY
  157.       NESP=PROPHY.CV(/1)-1
  158.       SEGACT PROPH2
  159.       DO IESP=1,NESP
  160.          IF (LCLIMR) THEN
  161.             SEGACT KRRIMP
  162.             SEGACT MPRIMP
  163.          ENDIF
  164.          LCLIMY=PROPH2.LCLIM(IESP)
  165.          IF (LCLIMY) THEN
  166.             KRYIMP=PROPH2.KRCLIM(IESP)
  167.             SEGACT KRYIMP
  168.          ENDIF
  169.          SEGACT NOMINC
  170.          SEGACT KRCENT
  171.          SEGACT KRFACE
  172.          SEGACT MELEFL
  173.          SEGACT MPSURF
  174.          SEGACT MPNORM
  175.          SEGACT MPVOLU
  176.          MPDKC=PROPH2.MPDIFF(IESP)
  177.          SEGACT MPDKC
  178.          SEGACT MPROC
  179.          MPYKC=PROPH2.MPVALY(IESP)
  180.          SEGACT MPYKC
  181.          ICOGRY=PROPHY.CGRYK(IESP)
  182.          SEGACT ICOGRY
  183.          NSOUCH=ICOGRY.IMACHE(/1)
  184.          DO 1 ISOUCH=1,NSOUCH
  185.             MCOGRY=ICOGRY.IMACHE(ISOUCH)
  186.             JCOGRY=ICOGRY.ICHAML(ISOUCH)
  187.             SEGACT JCOGRY
  188.             KDYKDX=JCOGRY.IELVAL(1)
  189.             KDYKDY=JCOGRY.IELVAL(2)
  190.             KDYKDZ=JCOGRY.IELVAL(3)
  191.             SEGDES JCOGRY
  192.             SEGACT KDYKDX
  193.             SEGACT KDYKDY
  194.             SEGACT KDYKDZ
  195.             SEGACT MCOGRY
  196.             NELEM=MCOGRY.NUM(/2)
  197.             NPTEL=MCOGRY.NUM(/1)
  198.             NPP1=NPTEL-1
  199.             NPP2=NPTEL-1
  200.             NEL1=NELEM
  201.             NEL2=NELEM
  202.             IEL1=1
  203.             IEL2=1
  204.             SEGINI RYKRHO
  205.             SEGINI RYKRYK
  206.             SEGINI CTGRY
  207.             RYKRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  208.             RYKRHO.NOMPRI(5:8)='    '
  209.             RYKRHO.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
  210.             RYKRHO.NOMDUA(5:8)='    '
  211.             RYKRYK.NOMPRI(1:4)=NOMINC.MOTS(IDIM+2+IESP)
  212.             RYKRYK.NOMPRI(5:8)='    '
  213.             RYKRYK.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
  214.             RYKRYK.NOMDUA(5:8)='    '
  215.             DO 12 IELEM=1,NELEM
  216. *     Le premier point du support de ICOGRY est un point FACE
  217.                NGFACE=MCOGRY.NUM(1,IELEM)
  218.                NLFACE=KRFACE.LECT(NGFACE)
  219.                IF (NLFACE.EQ.0) THEN
  220.                   WRITE(IOIMP,*) 'Erreur de programmation n°1'
  221.                   GOTO 9999
  222.                ENDIF
  223. *     On calcule la contribution à la matrice jacobienne IJACO de la face
  224. *     NGFAC (points duaux     : centres à gauche et à droite de la face)
  225. *           (points primaux   : une partie (bicoz conditions aux limites)
  226. *            de ceux du stencil pour le calcul du gradient
  227. *            à la face, ils doivent être des points centres)
  228. *       Si le flux de chaleur sur la face est imposé par les conditions
  229. *       aux limites, la contribution de la face à IJACO est nulle.
  230.                NGCGAU=MELEFL.NUM(1,NLFACE)
  231.                NGCDRO=MELEFL.NUM(3,NLFACE)
  232.                NLCGAU=KRCENT.LECT(NGCGAU)
  233.                NLCDRO=KRCENT.LECT(NGCDRO)
  234.                LMUR=(NGCGAU.EQ.NGCDRO)
  235. *     On distingue le cas où la face est un bord du maillage (mur)
  236. *                  du cas où la face est interne au maillage
  237.                IF (.NOT.LMUR) THEN
  238.                   NPD=2
  239.                   ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
  240.                   XF = MCOORD.XCOOR(ICOORX)
  241.                   YF = MCOORD.XCOOR(ICOORX+1)
  242.                   ZF = MCOORD.XCOOR(ICOORX+2)
  243.                   ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
  244.                   XG = MCOORD.XCOOR(ICOORX)
  245.                   YG = MCOORD.XCOOR(ICOORX+1)
  246.                   ZG = MCOORD.XCOOR(ICOORX+2)
  247.                   XFMXG = XF - XG
  248.                   YFMYG = YF - YG
  249.                   ZFMZG = ZF - ZG
  250.                   DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
  251.                   ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
  252.                   XD = MCOORD.XCOOR(ICOORX)
  253.                   YD = MCOORD.XCOOR(ICOORX+1)
  254.                   ZD = MCOORD.XCOOR(ICOORX+2)
  255.                   XFMXD = XF - XD
  256.                   YFMYD = YF - YD
  257.                   ZFMZD = ZF - ZD
  258.                   DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
  259.                   ALPHA=DRG/(DRG+DRD)
  260.                   UMALPH= 1.0D0 - ALPHA
  261.                ELSE
  262.                   NPD=1
  263.                   ALPHA=0.0D0
  264.                   UMALPH=1.0D0
  265.                ENDIF
  266.                IF (LCLIMR) THEN
  267.                   NLFRI=KRRIMP.LECT(NGFACE)
  268.                ELSE
  269.                   NLFRI=0
  270.                ENDIF
  271.                IF (NLFRI.GT.0) THEN
  272.                   RHOF=MPRIMP.VPOCHA(NLFRI,1)
  273.                ELSE
  274.                   RHOG = MPROC.VPOCHA(NLCGAU,1)
  275.                   RHOD = MPROC.VPOCHA(NLCDRO,1)
  276.                   RHOF = UMALPH * RHOG + ALPHA * RHOD
  277.                ENDIF
  278.                DKG  = MPDKC.VPOCHA(NLCGAU,1)
  279.                DKD  = MPDKC.VPOCHA(NLCDRO,1)
  280.                DKF  = UMALPH * DKG + ALPHA * DKD
  281.                NPP=NPTEL-1
  282. *     IPD=1 : point à gauche du point NGFACE
  283. *     IPD=2 : point à droite du point NGFACE
  284.                DO 122 IPD=1,NPD
  285.                   NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
  286.                   IF (.NOT.LMUR) THEN
  287.                      CTGRY.POIDU2(IPD,IEL2)=NPDUAL
  288.                   ELSE
  289.                      CTGRY.POIDU1(IPD,IEL1)=NPDUAL
  290.                   ENDIF
  291.                   NLCEND=KRCENT.LECT(NPDUAL)
  292.                   IF (NLCEND.EQ.0) THEN
  293.                      WRITE(IOIMP,*) 'Erreur grave n°1'
  294.                      GOTO 9999
  295.                   ENDIF
  296.                   DO 124 IPP=1,NPP
  297.                      NPPRIM=MCOGRY.NUM(IPP+1,IELEM)
  298.                      LCTRB2=.TRUE.
  299.                      IF (LCLIMY) THEN
  300.                         NLCLY=KRYIMP.LECT(NPPRIM)
  301.                         IF (NLCLY.NE.0) THEN
  302.                            LCTRB2=.FALSE.
  303.                         ENDIF
  304.                      ENDIF
  305.                      IF (.NOT.LCTRB2) THEN
  306. *     Lorsque une contribution est nulle, on fixe artificiellement le
  307. *     point primal égal au point dual.
  308.                         IF (.NOT.LMUR) THEN
  309.                            CTGRY.POIPR2(IPP,IEL2)=NPDUAL
  310.                            RYKRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  311.                            RYKRYK.VALMA2(IPD,IPP,IEL2)=0.D0
  312.                         ELSE
  313.                            CTGRY.POIPR1(IPP,IEL1)=NPDUAL
  314.                            RYKRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  315.                            RYKRYK.VALMA1(IPD,IPP,IEL1)=0.D0
  316.                         ENDIF
  317.                      ELSE
  318. *     Les contributions valent :
  319. *     (d Res_{\rho Yk})_d / (d var)_p =
  320. *     +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \rho_f * Dk_f
  321. *     * [ ((n_x * \beta_x) +(n_y * \beta_y) +(n_z * \beta_z)) *
  322. *         ((dYk)_p / (d var)_p)]
  323. *     avec :
  324. *          (dYk)_p / (d \rho)_p = - Yk_p / \rho_p
  325. *          (dYk)_p / (d \rho Yk)_p =   1 / \rho_p
  326. *     \beta_x : coefficients pour le calcul de dYk/dx
  327. *     \beta_y : coefficients pour le calcul de dYk/dy
  328. *     \beta_z : coefficients pour le calcul de dYk/dz
  329. *
  330.                         NLCENP=KRCENT.LECT(NPPRIM)
  331.                         IF (NLCENP.EQ.0) THEN
  332.                            WRITE(IOIMP,*) 'Erreur grave n°2'
  333.                            GOTO 9999
  334.                         ENDIF
  335. *       normale sortante pour IPD=1, rentrante pour IPD=2
  336.                         SIGNOR=(-1.D0)**(IPD+1)
  337.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  338.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  339.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  340.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  341.                         CNZ   =MPNORM.VPOCHA(NLFACE,3)
  342.                         BETAX =KDYKDX.VELCHE(IPP+1,IELEM)
  343.                         BETAY =KDYKDY.VELCHE(IPP+1,IELEM)
  344.                         BETAZ =KDYKDZ.VELCHE(IPP+1,IELEM)
  345.                         RHOP  =MPROC.VPOCHA(NLCENP,1)
  346.                         YKP   =MPYKC.VPOCHA(NLCENP,1)
  347. *                        YKP   =MPYKC.VPOCHA(NLCEND,1)
  348.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*RHOF*DKF
  349.      $                       *((CNX*BETAX)+(CNY*BETAY)+(CNZ*BETAZ))
  350.                         DYKDRO=-YKP/RHOP
  351.                         DYKDRY=1.D0/RHOP
  352.                         IF (.NOT.LMUR) THEN
  353.                            CTGRY.POIPR2(IPP,IEL2)=NPPRIM
  354.                            RYKRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRO
  355.                            RYKRYK.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRY
  356.                         ELSE
  357.                            CTGRY.POIPR1(IPP,IEL1)=NPPRIM
  358.                            RYKRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRO
  359.                            RYKRYK.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRY
  360.                         ENDIF
  361.                      ENDIF
  362.  124              CONTINUE
  363.  122           CONTINUE
  364.                IF (.NOT.LMUR) THEN
  365.                   IEL2=IEL2+1
  366.                ELSE
  367.                   IEL1=IEL1+1
  368.                ENDIF
  369.  12         CONTINUE
  370.             NPP1=NPTEL-1
  371.             NPP2=NPTEL-1
  372.             NEL1=IEL1-1
  373.             NEL2=IEL2-1
  374.             SEGADJ RYKRHO
  375.             SEGADJ RYKRYK
  376.             SEGADJ CTGRY
  377.             CTGRY.LMATSI(**)=RYKRHO
  378.             CTGRY.LMATSI(**)=RYKRYK
  379. * On accumule les matrices résultantes dans IJACO
  380.             CALL AJMTK(CTGRY,IJACO,IMPR,IRET)
  381.             IF (IRET.NE.0) GOTO 9999
  382.             SEGSUP RYKRHO
  383.             SEGSUP RYKRYK
  384.             SEGSUP CTGRY
  385. *
  386.             SEGDES MCOGRY
  387.             SEGDES KDYKDZ
  388.             SEGDES KDYKDY
  389.             SEGDES KDYKDX
  390.  1       CONTINUE
  391.          SEGDES ICOGRY
  392.          SEGDES MPYKC
  393.          SEGDES MPROC
  394.          SEGDES MPDKC
  395.          SEGDES MPVOLU
  396.          SEGDES MPNORM
  397.          SEGDES MPSURF
  398.          SEGDES MELEFL
  399.          SEGDES KRFACE
  400.          SEGDES KRCENT
  401.          SEGDES NOMINC
  402.          IF (LCLIMY) THEN
  403.             SEGDES KRYIMP
  404.          ENDIF
  405.          IF (LCLIMR) THEN
  406.             SEGDES KRRIMP
  407.             SEGDES MPRIMP
  408.          ENDIF
  409.       ENDDO
  410.       SEGDES PROPH2
  411.       SEGDES PROPHY
  412. *
  413. * Normal termination
  414. *
  415.       IRET=0
  416.       RETURN
  417. *
  418. * Format handling
  419. *
  420. *
  421. * Error handling
  422. *
  423.  9999 CONTINUE
  424.       IRET=1
  425.       WRITE(IOIMP,*) 'An error was detected in subroutine zlap2b'
  426.       RETURN
  427. *
  428. * End of subroutine ZLAP2B
  429. *
  430.       END
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  

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