Télécharger phapro.eso

Retour à la liste

Numérotation des lignes :

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

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