Télécharger la plaquette Cast3M


fr en

Télécharger phapro.eso

Retour à la liste

Numérotation des lignes : 

phapro

  1. C PHAPRO    SOURCE    CB215821  20/11/25    13:35:44     10792          
  2.  
  3.       SUBROUTINE PHAPRO
  4.  
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7. *
  8. *  Calcule la nouvelle proportion de phase
  9. * en entree
  10. *           : son chamelem de materiau associé
  11. *           : un champ de temperature proposée avec les mult de
  12. *             lagrange
  13. *  en sortie MCHAML de proportion de phase mis a jour
  14. *
  15.  
  16. -INC PPARAM
  17. -INC CCOPTIO
  18. -INC CCREEL
  19.  
  20. -INC SMMODEL
  21. -INC SMCHAML
  22.       POINTEUR mchelp.mchelm, mchelq.mchelm
  23.       POINTEUR melvap.melval, melvaq.melval
  24. -INC SMCHPOI
  25. -INC SMELEME
  26. -INC SMCOORD
  27.  
  28.       segment icpr(npt)
  29.       segment,mtemp
  30.         integer icon(imax),ipos(ino),ivoi(nvoi)
  31.         real*8 qrest(ino),qtot(ino)
  32.       endsegment
  33.       segment iqmoye
  34.         real*8 qmoy(nbel)
  35.       endsegment
  36.       segment ipo2(ino)
  37.  
  38.       segment,mcorr
  39.         integer imelem(npartp),imchap(npartp),imoyep(npartp),
  40.      &          imchaq(npartq)
  41.       endsegment
  42.  
  43. C*      write(ioimp,*) 'Entree dans Phapro'
  44.       xpre   = xpetit
  45.       umxpre = 1.d0 - xpre
  46.  
  47. * Lecture du mchaml de proportions de phase
  48.       call LIROBJ('MCHAML  ',mchel1,1,iretou)
  49.       call ACTOBJ('MCHAML  ',mchel1,1)
  50.       if (ierr.ne.0) return
  51. * Lecture du mchaml de chaleur latente
  52.       call LIROBJ('MCHAML  ',mchelq,1,iretou)
  53.       call ACTOBJ('MCHAML  ',mchelq,1)
  54.       if (ierr.ne.0) return
  55. * Lecture du chpoint des reactions
  56.       call LIROBJ('CHPOINT ',mchpoi,1,iretou)
  57.       call ACTOBJ('CHPOINT ',mchpoi,1)
  58.       if (ierr.ne.0) return
  59.  
  60. * Verification des supports des MCHAMLS (aux noeuds !)
  61.       ipmodl = 0
  62.       call quesup(ipmodl,mchel1,0,0,isupp,iret)
  63.       if (isupp.gt.1) return
  64.       call quesup(ipmodl,mchelq,0,0,isupq,iret)
  65.       if (isupq.gt.1) return
  66.  
  67. * Creation du mchaml de proportions de phase resultat par recopie du
  68. * champ donne en entree
  69.       call copie8(mchel1,mchelp)
  70.  
  71. * On regarde d'abord le chpoint des reactions :
  72.       segact,mchpoi
  73.       ipachp = mchpoi.ipchp(/1)
  74. * S'il est vide, il n'y a rien a faire.
  75.       if (ipachp.eq.0) then
  76. C*        write(ioimp,*) 'chpoint vide, rien a faire.'
  77.         call ACTOBJ('MCHAML  ',mchelp,1)
  78.         call ECROBJ('MCHAML  ',mchelp)
  79.         goto 9900
  80. * Petite verification : une seule partition pour le chpoint des reactions
  81.       else if (ipachp.gt.1) then
  82.         call erreur(180)
  83.         goto 9900
  84.       endif
  85.       msoupo = mchpoi.ipchp(1)
  86.       segact,msoupo
  87. * Recherche de la composante obligatoire 'Q' des reactions
  88.       melchp = igeoc
  89.       mvpchp = ipoval
  90.       icochp = 0
  91.       do 100 ic = 1, nocomp(/2)
  92.         if (nocomp(ic).eq.'Q   ') then
  93.           if (icochp.ne.0) then
  94.             icochp = -1
  95.           else
  96.             icochp = ic
  97.           endif
  98.         endif
  99.  100  continue
  100.       if (icochp.le.0) then
  101.         if (icochp.lt.0)
  102.      &     write(ioimp,*) 'ERREUR : composante Q en double !'
  103.         call erreur(21)
  104.         goto 9900
  105.       endif
  106. * NB : par la suite, on ne travaillera plus que sur melchp et la
  107. *      icochp-eme composante de mvpchp !
  108.  
  109. C Informations de correspondance entre mchelp et mchelq :
  110.       segact,mchelp,mchelq
  111.       npartp = mchelp.imache(/1)
  112.       npartq = mchelq.imache(/1)
  113.  
  114.       segini,mcorr
  115.       do 101 ima = 1, npartp
  116.         meleme = mchelp.imache(ima)
  117.         do ico = 1, npartq
  118.           if (meleme.eq.mchelq.imache(ico)) then
  119.             imelem(ima) = ico
  120.             goto 101
  121.           endif
  122.         enddo
  123.         call erreur(103)
  124.         goto 9900
  125.  101  continue
  126.  
  127.       do 102 ima = 1, npartp
  128.         mchaml = mchelp.ichaml(ima)
  129.         segact,mchaml
  130.         do ico = 1, nomche(/2)
  131.           if (nomche(ico).eq.'PPHA') then
  132.             if (typche(ico).ne.'REAL*8') then
  133.               moterr(1:16) = typche(ico)
  134.               moterr(17:20) = 'PPHA'
  135.               moterr(21:36) = mchelp.titche
  136.               call erreur(552)
  137.               goto 9990
  138.             endif
  139.             imchap(ima) = ico
  140.             goto 102
  141.           endif
  142.         enddo
  143.         moterr(1:4) = 'PPHA'
  144.         call erreur(236)
  145.         goto 9990
  146.  102  continue
  147.  
  148.       do 103 ima = 1, npartq
  149.         mchaml = mchelq.ichaml(ima)
  150.         segact,mchaml
  151.         do ico = 1, nomche(/2)
  152.           if (nomche(ico).eq.'Q   ') then
  153.             if (typche(ico).ne.'REAL*8') then
  154.               moterr(1:16) = typche(ico)
  155.               moterr(17:20) = 'Q   '
  156.               moterr(21:36) = mchelq.titche
  157.               call erreur(552)
  158.               goto 9990
  159.             endif
  160.             imchaq(ima) = ico
  161.             goto 103
  162.           endif
  163.         enddo
  164.         moterr(1:4) = 'Q   '
  165.         call erreur(236)
  166.         goto 9990
  167.  103  continue
  168.  
  169. * Tableau des connexions a partir du mchaml de chaleur latente Q
  170. * icpr : icpr(i) donne le nombre d'elements du mchaml auxquels le noeud
  171. *        num. i (= 0 si le noeud n'est pas dans le mchaml)
  172.       npt = nbpts      
  173.       segini,icpr
  174.       ngrand = 0
  175.       do 110 ima = 1, npartq
  176.         meleme = mchelq.imache(ima)
  177.         segact,meleme
  178.         nbnn = meleme.num(/1)
  179.         nbel = meleme.num(/2)
  180.         ngrand = max(ngrand,nbel)
  181.         do 111 iel = 1, nbel
  182.           do 111 inn = 1, nbnn
  183.             ia = meleme.num(inn,iel)
  184.             icpr(ia) = icpr(ia)+1
  185.  111    continue
  186.  110  continue
  187.       ngrand = ngrand+1
  188. * ino  : nombre total de noeuds (differents) concernes par le mchaml + 1
  189. * imax : somme des icpr(i) pour i = 1 a npt
  190. * nvoi : nombre max d'elements auquel appartient un noeud + 1
  191.       ino  = 1
  192.       imax = 0
  193.       nvoi = 0
  194.       do 112 inn = 1, npt
  195.         ia = icpr(inn)
  196.         if (ia.ne.0) then
  197.           ino = ino + 1
  198.           imax = imax + ia
  199.           nvoi = max(ia,nvoi)
  200.         endif
  201.  112  continue
  202.       nvoi = nvoi+1
  203. * ipos :
  204. * icpr :
  205.       segini,mtemp
  206.       ia = 1
  207.       do 113 inn = 1, npt
  208.         if (icpr(inn).ne.0) then
  209.           ja = ia
  210.           ia = ia + 1
  211.           mtemp.ipos(ia) = mtemp.ipos(ja) + icpr(inn)
  212.           icpr(inn) = ja
  213.         endif
  214.  113  continue
  215. *
  216.       segini,ipo2
  217.       do 114 ima = 1, npartq
  218.         meleme = mchelq.imache(ima)
  219.         nbnn = meleme.num(/1)
  220.         nbel = meleme.num(/2)
  221.         nplu = ngrand*(ima-1)
  222.         do 114 iel = 1, nbel
  223.           do 114 inn = 1, nbnn
  224.             ia = meleme.num(inn,iel)
  225.             ib = icpr(ia)
  226.             ipo2(ib) = ipo2(ib) + 1
  227.             id = ipos(ib) + ipo2(ib)
  228.             icon(id) = iel + nplu
  229.   114 continue
  230.       segsup,ipo2
  231.  
  232.       do 115 ima = 1, npartp
  233.         inomp  = mcorr.imchap(ima)
  234.         meleme = mchelp.imache(ima)
  235.         mchaml = mchelp.ichaml(ima)
  236. c*        segact,meleme,mchaml <- Deja actives
  237.         melvap = mchaml.ielval(inomp)
  238.         ico    = mcorr.imelem(ima)
  239.         inomq  = mcorr.imchaq(ico)
  240. C**     meleme = mchelq.imache(ico) = mchelp.imache(ima)
  241.         mchaml = mchelq.ichaml(ico)
  242. c*        segact,meleme,mchaml <- Deja actives
  243.         melvaq = mchaml.ielval(inomq)
  244.         nbnn = meleme.num(/1)
  245.         nbel = meleme.num(/2)
  246.         segact,melvap*MOD,melvaq
  247.         n1ptel = melvap.velche(/1)
  248.         n1el   = melvap.velche(/2)
  249.         n2ptel = 0
  250.         n2el   = 0
  251. C* Champ de phase est constant :
  252.         if (n1ptel .lt. melvaq.velche(/1)) then
  253.           n1ptel = melvaq.velche(/1)
  254.           n1el   = nbel
  255.           segadj,melvap
  256. c*          write(ioimp,*) 'Segadj complet'
  257.           r_z = melvap.velche(1,1)
  258.           if (r_z.ne.0.D0) then
  259.             do iel = 1, nbel
  260.               do inn = 1, n1ptel
  261.                 melvap.velche(inn,iel) = r_z
  262.               enddo
  263.             enddo
  264.           endif
  265. C* Champ de phase est constant par element :
  266.         else if (n1el .lt. melvaq.velche(/2)) then
  267. c*        n1el   = melvaq.velche(/2)
  268.           n1el   = nbel
  269.           segadj,melvap
  270. c*          write(ioimp,*) 'Segadj element'
  271.           do inn = 1, n1ptel
  272.             r_z = melvap.velche(inn,1)
  273.             if (r_z.ne.0.D0) then
  274.               do iel = 2, n1el
  275.                 melvap.velche(inn,iel) = r_z
  276.               enddo
  277.             endif
  278.           enddo
  279.         endif
  280.         segini,iqmoye
  281.         mcorr.imoyep(ima) = iqmoye
  282.  115  continue
  283.  
  284.       meleme = melchp
  285.       mpoval = mvpchp
  286.       segact,meleme,mpoval
  287.  
  288. * fabrication de la proportion de phase deja mangee par element
  289.       DO 1000 ifois = 1, 2
  290.  
  291. C*        write(ioimp,*) 'IFOIS =',ifois
  292.         if (ifois.eq.2) then
  293.           do i = 1, ino
  294.             qrest(i) = 0.D0
  295.             qtot(i)  = 0.D0
  296.           enddo
  297.         endif
  298.  
  299.         do 570 ima = 1, npartp
  300.  
  301.           inomp  = mcorr.imchap(ima)
  302.           mchaml = mchelp.ichaml(ima)
  303.           melvap = mchaml.ielval(inomp)
  304.           imelp1 = melvap.velche(/1)
  305.           imelp2 = melvap.velche(/2)
  306.  
  307.           ico    = mcorr.imelem(ima)
  308.           inomq  = mcorr.imchaq(ico)
  309.           mchaml = mchelq.ichaml(ico)
  310.           melvaq = mchaml.ielval(inomq)
  311.           imelq1 = melvaq.velche(/1)
  312.           imelq2 = melvaq.velche(/2)
  313.  
  314.           iqmoye = mcorr.imoyep(ima)
  315.  
  316.           meleme = mchelp.imache(ima)
  317.           nbnn = meleme.num(/1)
  318.           nbel = meleme.num(/2)
  319.  
  320.           do 575 iel = 1, nbel
  321.             jmelp2 = min(iel,imelp2)
  322.             jmelq2 = min(iel,imelq2)
  323.             qs = 0.d0
  324.             do 574 inn = 1, nbnn
  325.               is = ABS(icpr(meleme.num(inn,iel)))
  326.               jmelp1 = min(inn,imelp1)
  327.               jmelq1 = min(inn,imelq1)
  328.               r_z  = melvap.velche(jmelp1,jmelp2)
  329.               r_z1 = melvaq.velche(jmelq1,jmelq2)
  330.               qs = qs + r_z
  331.               qtot(is)  = qtot(is) + r_z1
  332.               qrest(is) = qrest(is)+ r_z1*(1.d0-r_z)
  333.   574       continue
  334.             qmoy(iel) = qs/nbnn
  335.   575     continue
  336.   570   continue
  337.  
  338. * boucle sur les noeuds concernés on y passe deux fois, la premiere
  339. * pour traiter les noeuds qui basculent completement
  340.         meleme = melchp
  341.         mpoval = mvpchp
  342.  
  343.         nbechp = meleme.num(/2)
  344.  
  345.         do 600 iel = 1, nbechp
  346.           ia = meleme.num(1,iel)
  347.           ib = icpr(ia)
  348.           if (ib.LE.0) goto 600
  349. * on commence par comparer qreac avec qlat pour voir si tout est mangé
  350.           qreac = vpocha(iel,icochp)
  351.           ic    = ipos(ib)+1
  352.           id    = ipos(ib+1)
  353.           if (ifois.eq.2) go to 620
  354.             qres = qrest(ib)
  355.             if (qreac.lt.0.D0) then
  356.               if (-qreac.lt.qres*umxpre) goto 600
  357.               xval = 1.d0
  358.             else
  359.               qto = qtot(ib)
  360.               if (qreac.lt.(qto-qres)*umxpre) goto 600
  361.               xval = 0.d0
  362.             endif
  363. * tous les elements sont manges : icpr(ia) => -icpr(ia) (<0)
  364.             icpr(ia) = -ABS(ib)
  365.             do 610 ie = ic, id
  366.               il = icon(ie)
  367.               ilv = mod(il,ngrand)
  368.               ieme = (il-ilv) / ngrand + 1
  369.               inomp  = mcorr.imchap(ieme)
  370.               mchaml = mchelp.ichaml(ieme)
  371.               melvap = mchaml.ielval(inomp)
  372.               ipt1   = mchelp.imache(ieme)
  373.               nbnn   = ipt1.num(/1)
  374.               do inn = 1, nbnn
  375.                 if (ipt1.num(inn,ilv).eq.ia) then
  376.                   jmelp1 = min(inn,melvap.velche(/1))
  377.                   jmelp2 = min(ilv,melvap.velche(/2))
  378.                   melvap.velche(jmelp1,jmelp2) = xval
  379.                   goto 610
  380.                 endif
  381.               enddo
  382.   610       continue
  383.             goto 600
  384.   620     continue
  385. * il va falloir en tartiner sur certains elements
  386. * on prend d'abord l'élement dont la plus gande proportion est dans le
  387. * sens qui nous interesse
  388.           i_z = id-ic+1
  389.           do ivo = 1, i_z
  390.             ivoi(ivo) = 0
  391.           enddo
  392.  630      continue
  393.             iiem = -1
  394.             ieg  = 0
  395.             qma  = -1.5D0
  396.             qmi  = +1.5D0
  397. * ifon=1 quand on fond sinon on solidifie
  398.             if (qreac.lt.0.D0) then
  399.               ifon = 1
  400.               qreac = -qreac
  401.             else
  402.               ifon = 0
  403.             endif
  404.             do 502 ie = ic,id
  405.               if (ivoi(ie-ic+1).eq.1) goto 502
  406.               il  = icon(ie)
  407.               ilv = mod(il,ngrand)
  408.               ieme = (il - ilv)/ngrand + 1
  409.               iqmoye = mcorr.imoyep(ieme)
  410.               if (ifon.eq.1) then
  411.                 if (qmoy(ilv).lt.umxpre. and. qmoy(ilv).GT.qma) then
  412.                   iiem = ieme
  413.                   iilv = ilv
  414.                   ieg  = ie
  415.                   qma  = qmoy(ilv)
  416.                 endif
  417.               else
  418.                 if (qmoy(ilv).ge.xpre .and. qmoy(ilv).lt.qmi) then
  419.                   iiem = ieme
  420.                   iilv = ilv
  421.                   ieg  = ie
  422.                   qmi  = qmoy(ilv)
  423.                 endif
  424.               endif
  425.   502       continue
  426.             if (iiem.eq.-1) goto 600
  427. *
  428.             ivoi(ieg-ic+1) = 1
  429. *    mise a jour du point en question dans l'element trouvé
  430.             inomp  = mcorr.imchap(iiem)
  431.             mchaml = mchelp.ichaml(iiem)
  432.             melvap = mchaml.ielval(inomp)
  433.             imelp1 = melvap.velche(/1)
  434.             jmelp2 = min(iilv,melvap.velche(/2))
  435.             ico    = mcorr.imelem(iiem)
  436.             inomq  = mcorr.imchaq(ico)
  437.             mchaml = mchelq.ichaml(ico)
  438.             melvaq = mchaml.ielval(inomq)
  439.             imelq1 = melvaq.velche(/1)
  440.             jmelq2 = min(iilv,melvaq.velche(/2))
  441.  
  442.             ipt1  = mchelp.imache(iiem)
  443.             nbnn = ipt1.num(/1)
  444. C*            iqmoye = mcorr.imoyep(iiem)
  445.             do 505 inn = 1, nbnn
  446.               ir = ipt1.num(inn,iilv)
  447.               if (ir.ne.ia) goto 505
  448.               jmelp1 = min(inn,imelp1)
  449.               jmelq1 = min(inn,imelq1)
  450.               r_p = melvap.velche(jmelp1,jmelp2)
  451.               r_q = melvaq.velche(jmelq1,jmelq2)
  452.               if (ifon.eq.1) then
  453.                 qmax = r_q * (1.d0 - r_p)
  454.                 if (qmax.le.qreac) then
  455.                   melvap.velche(jmelp1,jmelp2) = 1.d0
  456.                   if ((qreac-qmax)/r_q.le.xpre) goto 600
  457.                   qreac = qmax - qreac
  458.                   goto 630
  459.                 else
  460.                   xpro = qreac / r_q
  461.                   melvap.velche(jmelp1,jmelp2) = r_p + xpro
  462.                 endif
  463.               else
  464.                 qmin = r_q * r_p
  465.                 if (qmin.le.qreac) then
  466.                   melvap.velche(jmelp1,jmelp2) = 0.d0
  467.                   if (ABS((qmin-qreac)/r_q).le.xpre) goto 600
  468.                   qreac = qreac - qmin
  469.                   goto 630
  470.                 else
  471.                   xpro = qreac / r_q
  472.                   melvap.velche(jmelp1,jmelp2) = r_p - xpro
  473.                 endif
  474.               endif
  475.               goto 600
  476.   505       continue
  477.   600   continue
  478.  
  479.  1000 CONTINUE
  480.  
  481. *  travail terminé il faut desactiver
  482.       meleme = melchp
  483.       mpoval = mvpchp
  484.  
  485. * desactivation mchelm,mchel1
  486.       do 961 ima = 1, npartp
  487.         inomp  = mcorr.imchap(ima)
  488.         meleme = mchelp.imache(ima)
  489.         mchaml = mchelp.ichaml(ima)
  490.         melvap = mchaml.ielval(inomp)
  491. * test de validité des valeurs 0.d0 <= val <= 1.d0
  492.         nbptp = melvap.velche(/1)
  493.         nbelp = melvap.velche(/2)
  494.         do 976 iel = 1, nbelp
  495.           do 976 inn = 1, nbptp
  496.             r_z  = melvap.velche(inn,iel)
  497.             if (r_z .LT. 0.D0) then
  498.               melvap.velche(inn,iel) = 0.d0
  499.             else if (r_z .GT. 1.D0) then
  500.               melvap.velche(inn,iel) = 1.d0
  501.             endif
  502.   976   continue
  503.   961 continue
  504.  
  505.       call ACTOBJ('MCHAML  ',mchelp,1)
  506.       call ECROBJ('MCHAML  ',mchelp)
  507.  
  508.  9990 continue
  509.       segsup,icpr,mcorr,mtemp
  510.  9900 continue
  511.  
  512.       end
  513.  
  514.  
  515.  
  516.  
  517.  
  518.  

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