Télécharger phapro.eso

Retour à la liste

Numérotation des lignes :

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

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