Télécharger la plaquette Cast3M


fr en

Télécharger zlap1b.eso

Retour à la liste

Numérotation des lignes : 

zlap1b

  1. C ZLAP1B    SOURCE    CB215821  20/11/25    13:45:05     10792          
  2.       SUBROUTINE ZLAP1B(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         : ZLAP1B
  12. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  13. C               VF 2D (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       : ZLAP1A : Calcul de la matrice jacobienne du
  31. C                    résidu du laplacien VF 2D.
  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 2D.
  55. C SORTIES            : -
  56. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  57. C***********************************************************************
  58. C VERSION    : v1, 22/02/2002, version initiale
  59. C HISTORIQUE : v1, 22/02/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
  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,CNX,CNY
  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,XG,YG,XFMXG,YFMYG,DRG
  150.      &     ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
  151. C
  152. *
  153. * Executable statements
  154. *
  155.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans zlap1b.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.             SEGDES JCOGRY
  191.             SEGACT KDYKDX
  192.             SEGACT KDYKDY
  193.             SEGACT MCOGRY
  194.             NELEM=MCOGRY.NUM(/2)
  195.             NPTEL=MCOGRY.NUM(/1)
  196.             NPP1=NPTEL-1
  197.             NPP2=NPTEL-1
  198.             NEL1=NELEM
  199.             NEL2=NELEM
  200.             IEL1=1
  201.             IEL2=1
  202.             SEGINI RYKRHO
  203.             SEGINI RYKRYK
  204.             SEGINI CTGRY
  205.             RYKRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  206.             RYKRHO.NOMPRI(5:8)='    '
  207.             RYKRHO.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
  208.             RYKRHO.NOMDUA(5:8)='    '
  209.             RYKRYK.NOMPRI(1:4)=NOMINC.MOTS(IDIM+2+IESP)
  210.             RYKRYK.NOMPRI(5:8)='    '
  211.             RYKRYK.NOMDUA(1:4)=NOMINC.MOTS(IDIM+2+IESP)
  212.             RYKRYK.NOMDUA(5:8)='    '
  213.             DO 12 IELEM=1,NELEM
  214. *     Le premier point du support de ICOGRY est un point FACE
  215.                NGFACE=MCOGRY.NUM(1,IELEM)
  216.                NLFACE=KRFACE.LECT(NGFACE)
  217.                IF (NLFACE.EQ.0) THEN
  218.                   WRITE(IOIMP,*) 'Erreur de programmation n°1'
  219.                   GOTO 9999
  220.                ENDIF
  221. *     On calcule la contribution à la matrice jacobienne IJACO de la face
  222. *     NGFAC (points duaux     : centres à gauche et à droite de la face)
  223. *           (points primaux   : une partie (bicoz conditions aux limites)
  224. *            de ceux du stencil pour le calcul du gradient
  225. *            à la face, ils doivent être des points centres)
  226. *       Si le flux de chaleur sur la face est imposé par les conditions
  227. *       aux limites, la contribution de la face à IJACO est nulle.
  228.                NGCGAU=MELEFL.NUM(1,NLFACE)
  229.                NGCDRO=MELEFL.NUM(3,NLFACE)
  230.                NLCGAU=KRCENT.LECT(NGCGAU)
  231.                NLCDRO=KRCENT.LECT(NGCDRO)
  232.                LMUR=(NGCGAU.EQ.NGCDRO)
  233. *     On distingue le cas où la face est un bord du maillage (mur)
  234. *                  du cas où la face est interne au maillage
  235.                IF (.NOT.LMUR) THEN
  236.                   NPD=2
  237.                   ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
  238.                   XF = MCOORD.XCOOR(ICOORX)
  239.                   YF = MCOORD.XCOOR(ICOORX+1)
  240.                   ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
  241.                   XG = MCOORD.XCOOR(ICOORX)
  242.                   YG = MCOORD.XCOOR(ICOORX+1)
  243.                   XFMXG = XF - XG
  244.                   YFMYG = YF - YG
  245.                   DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
  246.                   ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
  247.                   XD = MCOORD.XCOOR(ICOORX)
  248.                   YD = MCOORD.XCOOR(ICOORX+1)
  249.                   XFMXD = XF - XD
  250.                   YFMYD = YF - YD
  251.                   DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
  252.                   ALPHA=DRG/(DRG+DRD)
  253.                   UMALPH= 1.0D0 - ALPHA
  254.                ELSE
  255.                   NPD=1
  256.                   ALPHA=0.0D0
  257.                   UMALPH=1.0D0
  258.                ENDIF
  259.                IF (LCLIMR) THEN
  260.                   NLFRI=KRRIMP.LECT(NGFACE)
  261.                ELSE
  262.                   NLFRI=0
  263.                ENDIF
  264.                IF (NLFRI.GT.0) THEN
  265.                   RHOF=MPRIMP.VPOCHA(NLFRI,1)
  266.                ELSE
  267.                   RHOG = MPROC.VPOCHA(NLCGAU,1)
  268.                   RHOD = MPROC.VPOCHA(NLCDRO,1)
  269.                   RHOF = UMALPH * RHOG + ALPHA * RHOD
  270.                ENDIF
  271.                DKG  = MPDKC.VPOCHA(NLCGAU,1)
  272.                DKD  = MPDKC.VPOCHA(NLCDRO,1)
  273.                DKF  = UMALPH * DKG + ALPHA * DKD
  274.                NPP=NPTEL-1
  275. *     IPD=1 : point à gauche du point NGFACE
  276. *     IPD=2 : point à droite du point NGFACE
  277.                DO 122 IPD=1,NPD
  278.                   NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
  279.                   IF (.NOT.LMUR) THEN
  280.                      CTGRY.POIDU2(IPD,IEL2)=NPDUAL
  281.                   ELSE
  282.                      CTGRY.POIDU1(IPD,IEL1)=NPDUAL
  283.                   ENDIF
  284.                   NLCEND=KRCENT.LECT(NPDUAL)
  285.                   IF (NLCEND.EQ.0) THEN
  286.                      WRITE(IOIMP,*) 'Erreur grave n°1'
  287.                      GOTO 9999
  288.                   ENDIF
  289.                   DO 124 IPP=1,NPP
  290.                      NPPRIM=MCOGRY.NUM(IPP+1,IELEM)
  291.                      LCTRB2=.TRUE.
  292.                      IF (LCLIMY) THEN
  293.                         NLCLY=KRYIMP.LECT(NPPRIM)
  294.                         IF (NLCLY.NE.0) THEN
  295.                            LCTRB2=.FALSE.
  296.                         ENDIF
  297.                      ENDIF
  298.                      IF (.NOT.LCTRB2) THEN
  299. *     Lorsque une contribution est nulle, on fixe artificiellement le
  300. *     point primal égal au point dual.
  301.                         IF (.NOT.LMUR) THEN
  302.                            CTGRY.POIPR2(IPP,IEL2)=NPDUAL
  303.                            RYKRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  304.                            RYKRYK.VALMA2(IPD,IPP,IEL2)=0.D0
  305.                         ELSE
  306.                            CTGRY.POIPR1(IPP,IEL1)=NPDUAL
  307.                            RYKRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  308.                            RYKRYK.VALMA1(IPD,IPP,IEL1)=0.D0
  309.                         ENDIF
  310.                      ELSE
  311. *     Les contributions valent :
  312. *     (d Res_{\rho Yk})_d / (d var)_p =
  313. *     +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \rho_f * Dk_f
  314. *     * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
  315. *         ((dYk)_p / (d var)_p)]
  316. *     avec :
  317. *          (dYk)_p / (d \rho)_p = - Yk_p / \rho_p
  318. *          (dYk)_p / (d \rho Yk)_p =   1 / \rho_p
  319. *     \beta_x : coefficients pour le calcul de dYk/dx
  320. *     \beta_y : coefficients pour le calcul de dYk/dy
  321. *
  322.                         NLCENP=KRCENT.LECT(NPPRIM)
  323.                         IF (NLCENP.EQ.0) THEN
  324.                            WRITE(IOIMP,*) 'Erreur grave n°2'
  325.                            GOTO 9999
  326.                         ENDIF
  327. *       normale sortante pour IPD=1, rentrante pour IPD=2
  328.                         SIGNOR=(-1.D0)**(IPD+1)
  329.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  330.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  331.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  332.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  333.                         BETAX =KDYKDX.VELCHE(IPP+1,IELEM)
  334.                         BETAY =KDYKDY.VELCHE(IPP+1,IELEM)
  335.                         RHOP  =MPROC.VPOCHA(NLCENP,1)
  336.                         YKP   =MPYKC.VPOCHA(NLCENP,1)
  337. *                        YKP   =MPYKC.VPOCHA(NLCEND,1)
  338.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*RHOF*DKF
  339.      $                       *((CNX*BETAX)+(CNY*BETAY))
  340.                         DYKDRO=-YKP/RHOP
  341.                         DYKDRY=1.D0/RHOP
  342.                         IF (.NOT.LMUR) THEN
  343.                            CTGRY.POIPR2(IPP,IEL2)=NPPRIM
  344.                            RYKRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRO
  345.                            RYKRYK.VALMA2(IPD,IPP,IEL2)=FACTOR*DYKDRY
  346.                         ELSE
  347.                            CTGRY.POIPR1(IPP,IEL1)=NPPRIM
  348.                            RYKRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRO
  349.                            RYKRYK.VALMA1(IPD,IPP,IEL1)=FACTOR*DYKDRY
  350.                         ENDIF
  351.                      ENDIF
  352.  124              CONTINUE
  353.  122           CONTINUE
  354.                IF (.NOT.LMUR) THEN
  355.                   IEL2=IEL2+1
  356.                ELSE
  357.                   IEL1=IEL1+1
  358.                ENDIF
  359.  12         CONTINUE
  360.             NPP1=NPTEL-1
  361.             NPP2=NPTEL-1
  362.             NEL1=IEL1-1
  363.             NEL2=IEL2-1
  364.             SEGADJ RYKRHO
  365.             SEGADJ RYKRYK
  366.             SEGADJ CTGRY
  367.             CTGRY.LMATSI(**)=RYKRHO
  368.             CTGRY.LMATSI(**)=RYKRYK
  369. * On accumule les matrices résultantes dans IJACO
  370.             CALL AJMTK(CTGRY,IJACO,IMPR,IRET)
  371.             IF (IRET.NE.0) GOTO 9999
  372.             SEGSUP RYKRHO
  373.             SEGSUP RYKRYK
  374.             SEGSUP CTGRY
  375. *
  376.             SEGDES MCOGRY
  377.             SEGDES KDYKDY
  378.             SEGDES KDYKDX
  379.  1       CONTINUE
  380.          SEGDES ICOGRY
  381.          SEGDES MPYKC
  382.          SEGDES MPROC
  383.          SEGDES MPDKC
  384.          SEGDES MPVOLU
  385.          SEGDES MPNORM
  386.          SEGDES MPSURF
  387.          SEGDES MELEFL
  388.          SEGDES KRFACE
  389.          SEGDES KRCENT
  390.          SEGDES NOMINC
  391.          IF (LCLIMY) THEN
  392.             SEGDES KRYIMP
  393.          ENDIF
  394.          IF (LCLIMR) THEN
  395.             SEGDES KRRIMP
  396.             SEGDES MPRIMP
  397.          ENDIF
  398.       ENDDO
  399.       SEGDES PROPH2
  400.       SEGDES PROPHY
  401. *
  402. * Normal termination
  403. *
  404.       IRET=0
  405.       RETURN
  406. *
  407. * Format handling
  408. *
  409. *
  410. * Error handling
  411. *
  412.  9999 CONTINUE
  413.       IRET=1
  414.       WRITE(IOIMP,*) 'An error was detected in subroutine zlap1b'
  415.       RETURN
  416. *
  417. * End of subroutine ZLAP1B
  418. *
  419.       END
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  

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