Télécharger murh.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : murh.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. * Thermique
  5. * Conduction dans un mur avec source et échange thermique.
  6. * Comparaison à la solution analytique en régime permanent.
  7. * Teste les opérateurs LAPN ECHI FIMP DFDT
  8. * La source volumique est imposée avec FIMP
  9. *
  10. * |y phi=0
  11. * --------------------
  12. * | |
  13. * | | phi=h(T-T0)
  14. *Tw | |
  15. * | |
  16. * -------------------->x
  17. * O phi=0 L
  18. *
  19. * Solution analytique:
  20. *
  21. * T(x)=Tw - q/2(X**2+LX)/lambda + (Tl-Tw)X/L
  22. *
  23. * avec Tl = (lambdaTw + LhT0)/(lambda+LH) + q/2 L**2/(lambda+LH)
  24. *
  25. *
  26. *$$$$ EXEC
  27. * EXEC PROCEDUR MAGN 03/03/31 21:15:04 4631
  28. *X EXEC (Procedure)
  29. * Procedure EXEC
  30. *
  31. * Objet : Execute un algorithme décrit dans une table RV
  32. * de type EQEX.
  33. * Cette table est créée par l'opérateur EQEX
  34. *
  35. * Syntaxe : EXEC RV ;
  36. *
  37. ***********************************************************************
  38. * VERSION : ????
  39. * HISTORIQUE : 20/12/99: gounand
  40. * Rajout de la gestion de la matrice servant à l'assemblage
  41. * (rv . 'METHINV' . 'MATASS')
  42. *
  43. * HISTORIQUE : 10/11/00: magnaud
  44. * Reecriture de exec pour tirer parti du // (operateur ASSI et RESO)
  45. * Pour l'instant cette procedure est testee experimentalement dans
  46. * la directive EQUA de EXEC
  47. *
  48. * HISTORIQUE : 04/04/01: magnaud
  49. * Calcul et impression tout les FIDT cycles de la norme Linf du residu
  50. *
  51. * HISTORIQUE :
  52. * HISTORIQUE :
  53. ************************************************************************
  54. * DISCPRES = (d'apres KMIC)
  55. * KPRE=3 pression P0 KPRE=4 pression P1 KPRE=2 cas macro 1ère génération
  56. * KPRE=5 pression MSOMMET
  57. ************************************************************************
  58. 'DEBPROC' EXEC ;
  59. 'ARGUMENT' rv*'TABLE ' ;
  60.  
  61. * Logique pilotant le recalcul de la matrice de pression
  62. CALPRE=FAUX ;
  63. Si ('EXIST' rv 'CALPRE') ;
  64. vertytab rv 'CALPRE' 'LOGIQUE' ;
  65. CALPRE=rv.'CALPRE' ;
  66. Finsi ;
  67.  
  68. Si(NON('EXIST' rv 'XEQUA'));
  69. rv.'XEQUA'= FAUX ;
  70. Finsi ;
  71.  
  72. ********************************************************************************
  73. ********************************************************************************
  74. ****** Directive EQUA **********************************************************
  75. ********************************************************************************
  76. ********************************************************************************
  77.  
  78. Si (rv.'XEQUA');
  79.  
  80. IMPKRES=0 ;
  81. Si ('EXIST' rv 'NUMASSI'); numassi=rv.'NUMASSI';
  82. Sinon; numassi=1 ;
  83. Finsi;
  84.  
  85.  
  86. Xdfdt = FAUX ;
  87.  
  88. ITMA=rv . 'ITMA';
  89. 'SI' ('<EG' ITMA 1) ;
  90. ITMA=1 ;
  91. 'FINSI' ;
  92.  
  93. NITER=rv . 'NITER';
  94. 'SI' ('<EG' NITER 1) ;
  95. NITER=1 ;
  96. 'FINSI' ;
  97.  
  98. OMEGA=rv . 'OMEGA';
  99. 'SI' ('<EG' NITER 1) ;
  100. OMEGA=1.;
  101. 'FINSI' ;
  102.  
  103. NUPDT=rv.'PASDETPS'.'NUPASDT';
  104.  
  105. 'SI'((NON (rv.'TRAN')) et ('>' ITMA 1) );
  106. Mess ' Incoherence : Algorithme Permanent et ITMA >1';
  107. QUITTER EXEC ;
  108. 'FINSI' ;
  109.  
  110. 'SI'((rv.'TRAN') et (NON (rv.'PROJ')) et ('EGA' NUPDT 1));
  111. Mess ' Algorithme transitoire standard implicite ou explicite ';
  112. mess '========================================================';
  113. 'FINSI' ;
  114.  
  115. 'SI'((rv.'PROJ') et ('EGA' NUPDT 1));
  116. mess 'Algorithme de Projection';
  117. mess '========================';
  118. 'FINSI' ;
  119.  
  120. 'SI' ( non (rv.'TRAN'));
  121. mess 'Algorithme de recherche d un permanent ';
  122. mess '=======================================';
  123. 'FINSI' ;
  124.  
  125.  
  126. 'REPETER' BCLTPS ITMA ;
  127. si (rv.'PROJ');
  128. testp1 = non ( exist rv 'MATC') ;
  129. sinon ;
  130. testp1 = FAUX ;
  131. finsi ;
  132. 'REPETER' BCLINT NITER ;
  133. Si (non (rv.'TRAN'));
  134. mess 'ITERATION ' &BCLINT;
  135. Finsi ;
  136.  
  137. Si (rv.'XRIG');
  138. MRIG='RGDT ';
  139. sf mau = assi numassi 'KOPS' 'RIGIDITE' ;
  140. sp map = assi numassi 'KOPS' 'RIGIDITE' ;
  141. sr mar = assi numassi 'KOPS' 'RIGIDITE' ;
  142. Sinon ;
  143. MRIG=TEXT ' ';
  144. sf mau = assi numassi 'KOPS' 'MATRIK' ;
  145. sp map = assi numassi 'KOPS' 'MATRIK' ;
  146. sr mar = assi numassi 'KOPS' 'MATRIK' ;
  147. Finsi ;
  148.  
  149. NRES=0 ;
  150.  
  151. *****************************************************
  152. * ._._. \|/
  153. * ._._. CALCUL MATRICES ELEMENTAIRES
  154. * ._._. /|\
  155. *****************************************************
  156.  
  157.  
  158. 'REPETER' bloc2 ('DIME' (rv . 'LISTOPER')) ;
  159. nomper = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ;
  160. * mess ' Début BLOC2 nomper=' nomper ;
  161. notable= 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ;
  162. Si('>EG' (rv.'IMPR') 2); mess 'Operateur ' nomper ; Finsi ;
  163. td=rv.notable ;
  164. tdk=td.'KOPT' ;
  165. iarg= rv.notable.'IARG' ;
  166. arg1 = TEXT ' ' ; arg2 = TEXT ' ' ; arg3 = TEXT ' ' ;
  167. arg4 = TEXT ' ' ; arg5 = TEXT ' ' ; arg6 = TEXT ' ' ;
  168. arg7 = TEXT ' ' ; arg8 = TEXT ' ' ; arg9 = TEXT ' ' ;
  169.  
  170. Si ('>EG' iarg 1 );
  171. Si ( EGA ('TYPE' (td.'ARG1')) 'MOT ' ) ;
  172. arg1=rv.'INCO'.(td.'ARG1') ;
  173. Sinon ;
  174. arg1=(td.'ARG1') ;
  175. Finsi ;
  176. Finsi ;
  177.  
  178. Si ('>EG' iarg 2 );
  179. Si ( EGA ('TYPE' (td.'ARG2')) 'MOT ' ) ;
  180. arg2=rv.'INCO'.(td.'ARG2') ;
  181. Sinon ;
  182. arg2=(td.'ARG2') ;
  183. Finsi ;
  184. Finsi ;
  185.  
  186. Si ('>EG' iarg 3 );
  187. Si ( EGA ('TYPE' (td.'ARG3')) 'MOT ' ) ;
  188. arg3=rv.'INCO'.(td.'ARG3') ;
  189. Sinon ;
  190. arg3=(td.'ARG3') ;
  191. Finsi ;
  192. Finsi ;
  193.  
  194. Si ('>EG' iarg 4 );
  195. Si ( EGA ('TYPE' (td.'ARG4')) 'MOT ' ) ;
  196. arg4=rv.'INCO'.(td.'ARG4') ;
  197. Sinon ;
  198. arg4=(td.'ARG4') ;
  199. Finsi ;
  200. Finsi ;
  201.  
  202. * Traitement du cas transitoire
  203. XTRA =CHAI 'PERM' ; DT = TEXT ' ' ; SCHE=TEXT ' ' ;
  204. inc1 = TEXT ' ' ; inc2 = TEXT ' ' ;
  205. imc1 = TEXT ' ' ; imc2 = TEXT ' ' ;
  206. BETAT= TEXT ' ' ;
  207. XBDF2 = FAUX ;
  208.  
  209. XOT=(('EGA' nomper 'NS ') ou ('EGA' nomper 'TSCA ')
  210. ou ('EGA' nomper 'DFDT ') ou ('EGA' nomper 'TSCAL '));
  211.  
  212. Si (rv.'PROJ' et ('EGA' nomper 'KBBT ')) ;
  213. XTRA=CHAI 'PROJ' ;
  214. Finsi ;
  215.  
  216. Si (rv.'TRAN') ;
  217. Xdfdt = VRAI ;
  218. XTRA=CHAI 'TRANS' ;
  219. Si (rv.'PROJ') ; XTRA=CHAI 'PROJ' ; Finsi ;
  220. DT=(rv.'DELTAT') ;
  221. Si ( EGA ('TYPE' DT) 'MOT ' ) ;
  222. DT=rv.'INCO'.DT ;
  223. Finsi ;
  224. SCHE = rv.'SCHEMAT';
  225. Si (EGA SCHE 'SEMI');
  226. BETAT=(rv.'Betat') ;
  227. Finsi ;
  228. Si (EGA SCHE 'BDF2');
  229. XBDF2 = VRAI ;
  230. Finsi ;
  231. Finsi ;
  232.  
  233. nbic=dime td.'LISTINCO';
  234.  
  235. Si ('>EG' nbic 1 );
  236. ninc1=extr td.'NUMEINCO' 1 ;
  237. a= chai 'FORMAT' '(A3,I1)' 'INC' ninc1 ;
  238. inc1=rv.'INCO'.(rv.a);
  239. Si (XBDF2) ;
  240. a= chai 'FORMAT' '(A3,I1)' 'IMC' ninc1 ;
  241. imc1=rv.'INCO'.(rv.a);
  242. Finsi;
  243. Finsi ;
  244.  
  245. Si ('>EG' nbic 2 );
  246. ninc2=extr td.'NUMEINCO' 2 ;
  247. a= chai 'FORMAT' '(A3,I1)' 'INC' ninc2 ;
  248. inc2=rv.'INCO'.(rv.a);
  249. Si (XBDF2) ;
  250. a= chai 'FORMAT' '(A3,I1)' 'IMC' ninc2 ;
  251. imc2=rv.'INCO'.(rv.a);
  252. Finsi;
  253. Finsi ;
  254.  
  255.  
  256. MM=TEXT ' ';
  257. DCEN='????' ;
  258. Si(EGA (tdk.'MTRMASS') 2 );MM='MMDIAGO';Finsi;
  259. Si(EGA (tdk.'IDCEN') 1 );DCEN='CENTREE';Finsi;
  260. Si(EGA (tdk.'IDCEN') 2 );DCEN='SUPGDC ';Finsi;
  261. Si(EGA (tdk.'IDCEN') 3 );DCEN='SUPG ';Finsi;
  262. Si(EGA (tdk.'IKOMP') 0 );MC=TEXT ' ';Finsi;
  263. Si(EGA (tdk.'IKOMP') 1 );MC='CONS';Finsi;
  264.  
  265. INEFMD= tdk.'INEFMD' ;
  266. tdz=TEXT ' ' ;
  267. Si (EGA nomper 'VNIMP ');
  268. tdz=CHAI 'V' ;
  269. nomper= text 'NS';
  270. Finsi ;
  271.  
  272. * temp place ;
  273. * mess ' Avant ' nomper ' tdz ' tdz ' DT=' DT ' XTRA=' ; list xtra ;
  274. * mess SCHE ; mess BETAT;
  275. msi mai =assi numassi nomper tdz
  276. (td.'DOMZ') (td.'LISTINCO') (td.'TYPEINCO') (rv.'TYPPRESS')
  277. INEFMD MRIG MM 'CONS' 'STABP' (tdk.'STABP') 'CMD' (tdk.'CMD') DCEN
  278. XTRA DT SCHE BETAT
  279. 'INCO' inc1 imc1 inc2 imc2
  280. (td.'IARG') arg1 arg2 arg3 arg4 ;
  281. * temp place ;
  282. * mess ' Apres ' nomper;
  283.  
  284. Si ((rv.'PROJ') et ('EGA' nomper 'KBBT '));
  285. msp =assi numassi msi ;
  286. mac =assi numassi mai ;
  287. 'ITER' bloc2 ;
  288. Finsi ;
  289.  
  290. Si ((rv.'PROJ') et ('EGA' nomper 'DFDT '));
  291. mam =assi numassi mai ;
  292. Finsi ;
  293.  
  294. * Si ((rv.'PROJ') et ('EGA' nomper 'TOIMP '));
  295. * list msi ;
  296. * Finsi ;
  297.  
  298. Si ((rv.'PROJ') et ('EGA' nomper 'NS '));
  299. XTRA=CHAI 'TRANS' ;
  300. DT=(rv.'DELTAT') ;
  301. Si ( EGA ('TYPE' DT) 'MOT ' ) ;
  302. DT=rv.'INCO'.DT ;
  303. Finsi ;
  304. SCHE = rv.'SCHEMAT';
  305. Si (EGA SCHE 'SEMI');
  306. BETAT=(rv.'Betat') ;
  307. Finsi ;
  308. Si (EGA SCHE 'BDF2');
  309. XBDF2 = VRAI ;
  310. Finsi ;
  311. * mess ' Appel DFDT DT=' DT ' XTRA=' ; list xtra ;
  312. msm mam =assi numassi DFDT
  313. (td.'DOMZ') (td.'LISTINCO') (td.'TYPEINCO') (rv.'TYPPRESS')
  314. INEFMD MRIG MM 'CONS' 'CMD' (tdk.'CMD') DCEN
  315. XTRA DT SCHE BETAT
  316. 'INCO' inc1 imc1 inc2 imc2
  317. 1 arg1 ;
  318. oubli msm;
  319. Finsi ;
  320.  
  321. mau = assi numassi mau 'ET' mai ;
  322. sf = assi numassi sf 'ET' msi ;
  323.  
  324.  
  325. rv.'TVNP'=FAUX;
  326. rv.'TVNPC'=FAUX;
  327.  
  328. Si ('EGA' nomper 'VNIMP '); *
  329. mess ' Cas oper VNIMP' ;
  330. Si (ega (rv.'TYPPRESS') 'MSOMMET');
  331. rv.'TVNPC'=VRAI;
  332. mvnpc=mai ;
  333. svnpc=msi ;
  334. Sinon ;
  335. rv.'TVNP'=VRAI ;
  336. mar = assi numassi mar 'ET' mai ;
  337. sr = assi numassi sr 'ET' msi ;
  338. Finsi ;
  339. Sinon ;
  340.  
  341. Finsi ;
  342.  
  343. 'FIN' bloc2 ;
  344. *mess 'Fin Bloc2 ' ;
  345.  
  346. 'SI'((rv.'TRAN') et (NON Xdfdt));
  347. Mess ' Incoherence : Algorithme Transitoire '
  348. ' et pas d increment en temps';
  349. Mess ' Il manque un operateur DFDT NS ou TSCA ';
  350. QUITTER EXEC ;
  351. 'FINSI' ;
  352.  
  353. NRES=NRES+1;
  354. *mess 'NRES=' NRES;
  355. TABRES=chai 'FORMAT' '(A6,I1)' 'TABRES' NRES;
  356. Si (NON('EXIST' rv TABRES)); rv.TABRES=TABLE ; Finsi ;
  357. Si (NON ('EGA' ('TYPE' rv.TABRES) 'TABLE '));
  358. rv.TABRES=TABLE ;
  359. Finsi ;
  360.  
  361. rv.TABRES.'mau' =mau ;
  362. rv.TABRES.'sf' =sf ;
  363. rv.TABRES.'clim'=rv.'CLIM' ;
  364. rv.TABRES.'METHINV'=rv.'METHINV' ;
  365.  
  366. *****************************************************
  367. * __ |\
  368. * \ /_/o | | \|/
  369. * \=> \_\/ \| --O-- PROJECTION
  370. * /|\ /\ /|\
  371. *************/*|*\***********************************
  372.  
  373. 'SI' (rv.'PROJ');
  374. SCHE = rv.'SCHEMAT';
  375. cnvi1=chai 'FORMAT' '(I1,A3)' 1 (extr rv.'LISTINCO' 1) ;
  376. cnvi2=chai 'FORMAT' '(I1,A3)' 2 (extr rv.'LISTINCO' 1) ;
  377. cnvi3=chai 'FORMAT' '(I1,A3)' 3 (extr rv.'LISTINCO' 1) ;
  378. cnppp=chai 'FORMAT' '(I1,A3)' 3 (extr rv.'LISTINCO' 2) ;
  379. 'SI' (testp1) ;
  380. mess ' CALCUL de CMCT ' ;
  381. Diag = doma (rv.'MODELE') 'XXDIAGSI' ;
  382. Digv= (nomc Diag cnvi1) et (nomc Diag cnvi2) ;
  383. Si ('EGA' (VALE DIME) 3);
  384. Digv= Digv et (nomc Diag cnvi3);
  385. Finsi ;
  386.  
  387. * Diag : Matrice diagonale de base ('SCAL')
  388. * Digv : Matrice diagonale de base ('1U' '2U' '3U')
  389. * Digvl: Matrice diagonale contenant les conditions aux limites
  390.  
  391. IDiag = 'INVERSE' Digv ;
  392.  
  393. Digvl= 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1;
  394. IDigvl= 'INVERSE' Digvl ;
  395.  
  396. *
  397. * Calcul matrices MD-1
  398. *
  399. MDM1='KMF' mam IDiag ;
  400. Si ('EGA' SCHE 'BDF2'); MDM1=MDM1 * (((rv.'DELTAT')*2.)/3.) ;
  401. Sinon ; MDM1=MDM1 * (rv.'DELTAT') ;
  402. Finsi ;
  403. *
  404. * Calcul matrices C
  405. *
  406. * msp mac =assi numassi 'KBBT'
  407. * (rv.'MODELE') (rv.'LISTINCO') (rv.'TYPEINCO') (rv.'TYPPRESS')
  408. * MRIG (1) (1.) ;
  409.  
  410. matpc mac = kops 'CMCTSPLT' mac ;
  411.  
  412. Si ('EGA' (rv.'TYPPRESS') 'MSOMMET');
  413. mess ' Cas des pressions continues' ;
  414. incp= mots (extr rv.'LISTINCO' 2) ;
  415. msp matpr =assi numassi 'LAPN'
  416. (rv.'MODELE') incp (mots 'TEMPERAT')
  417. MRIG 'LINE' (1) (1.) ;
  418.  
  419. Sinon ;
  420. mess ' Cas des pressions discontinues' ;
  421.  
  422. Si (rv.'TVNPC'); mac=mac et mvnpc; Finsi;
  423. rv . 'MATC' = mac ;
  424. * list mac 4 ;
  425. scr mcr = 'KOPS' 'MATRIK' ;
  426. Si (rv.'TVNP') ;
  427. rv . 'MBTR' = mar ;
  428. Dunit=Idigv;
  429. crt= kops 'CMCT' mac mar (Dunit) ;
  430. ctr= kops 'CMCT' mar mac (Dunit) ;
  431. rrt= kops 'CMCT' mar mar (Dunit) ;
  432. mcr= mcr et crt et ctr et rrt ;
  433. Finsi ;
  434.  
  435. * list mac 3 ;list IDigvl;
  436. matpr = kops 'CMCT' mac mac IDigvl;
  437. * mess 'APRES CMCT' ; list matpr 3 ;
  438. matpr = matpc et matpr et mcr ;
  439.  
  440. Finsi ;
  441.  
  442. rv.'IDigvl'= IDigvl ;
  443. rv.'MATP' = matpr ;
  444. rv.'MATC' = mac ;
  445. rv.'MDM1' = MDM1 ;
  446.  
  447. 'FINSI' ;
  448.  
  449. IDigvl= rv.'IDigvl';
  450. matpr = rv.'MATP' ;
  451. mac = rv.'MATC' ;
  452. MDM1 = rv.'MDM1' ;
  453.  
  454.  
  455. 'SI' ('EXIST' (rv.'INCO') 'PINM1') ;
  456. **
  457. ** mess ' t n-1' ;
  458. ** mess 'Calcul de C P' ;
  459. **
  460.  
  461. pinm1= rv.'INCO'.'PINM1';
  462. * cpi2=(extr rv.'LISTINCO' 2) ;
  463. * pinm1= nomc (rv.'INCO'.(rv.'INC2')) cpi2;
  464. cpre = 'KMF' (rv . 'MATC') pinm1 'TRAN' ;
  465. Si (rv.'TVNP') ;
  466. cxre = 'KMF' (rv . 'MBTR') pinm1 'TRAN' ;
  467. cpre = cpre et cxre ;
  468. Finsi ;
  469.  
  470. Si ('EGA' (rv.'TYPPRESS') 'MSOMMET');
  471. gradpres = 'KOPS' MDM1 '*' cpre ;
  472. Sinon ;
  473. gradpres = 'KOPS' MDM1 '*' cpre ;
  474. Finsi ;
  475. gradpres = nomc (mots 'UX' 'UY' 'UZ')
  476. (mots cnvi1 cnvi2 cnvi3) gradpres;
  477.  
  478. oublier cpre ;
  479. rv.'INCO'.'GRADPRES' = gradpres ;
  480. sf = sf - (rv.'INCO'.'GRADPRES') ;
  481. rv.'TABRES1'.'sf'=sf ;
  482.  
  483. 'FINSI' ;
  484. clim = exco rv.'CLIM' (mots cnvi1 cnvi2 cnvi3) 'NOID' ;
  485. rv.'TABRES1'.'clim'=clim ;
  486.  
  487. ********** FIN Traitement PRESSION ************************
  488. * mess ' FIN Traitement PRESSION' ;
  489.  
  490. 'FINSI' ;
  491. *******************************************************
  492.  
  493.  
  494. *****************************************************
  495. * \___\ /
  496. * \/ RESOLUTION
  497. * \ \|/
  498. *************************--O--***********************
  499. res mau = assi numassi 'KOPS' 'MATRIK' ;
  500. Si (rv.XRIG);
  501. res = chan res 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
  502. Finsi ;
  503. Repeter Bclres NRES ;
  504. * mess ' Resolution ' NRES ;
  505. TABRES=rv.(chai 'FORMAT' '(A6,I1)' 'TABRES' &Bclres);
  506. MTIV = TABRES.'METHINV';
  507. mau = TABRES.'mau' ;
  508. sf = TABRES.'sf' ;
  509. clim = TABRES.'clim';
  510. *........................................................
  511. Si (rv.'XRIG') ;
  512. gclim=extr clim 'MAILLAGE' ;
  513. rclim=bloq 'T' gclim ;
  514. dclim=depi rclim (nomc 'Q' clim);
  515. mau = mau et rclim ;
  516. sf = sf et dclim ;
  517. * list mau; list sf ;
  518. res1 = 'RESO' mau sf ;
  519. res = ASSI numassi res1 et res ;
  520. Sinon ;
  521. MATASS= mau ;
  522. TYRENU=MTIV.'TYRENU';PCMLAG=MTIV.'PCMLAG';KTYPI =MTIV.'TYPINV';
  523. XINIT=TEXT ' ' ;DI2 = TEXT ' ' ;DI3 = TEXT ' ' ;
  524. DI4 = TEXT ' ' ;DI5 = TEXT ' ' ;
  525. MAPREC=TEXT ' ';KPREC =TEXT ' ';PR3 = TEXT ' ' ;
  526. PR4 = TEXT ' ' ;PR5 = TEXT ' ' ;
  527. Si (('>EG' KTYPI 2) ET ('<EG' KTYPI 5));
  528. XINIT mar = assi numassi 'KOPS' 'MATRIK' ;
  529. DI2=MTIV.'NITMAX';DI3=MTIV.'RESID';DI4=MTIV.'BCGSBTOL';
  530. Si ('EGA' KTYPI 5);DI5=MTIV.'GMRESTRT';Finsi ;
  531. MAPREC= mau ;
  532. KPREC=MTIV.'PRECOND';
  533. Si ('EGA' KPREC 4);PR3=MTIV.'MILURELX';Finsi ;
  534. Si (('EGA' KPREC 5) OU ('EGA' KPREC 6));
  535. PR4=MTIV.'ILUTLFIL';PR5=MTIV.'ILUTDTOL';
  536. Finsi ;
  537. Finsi ;
  538. * list mau 3 ;
  539. * list sf ;
  540. res1 = ASSI numassi 'KRES' 'LL' mau
  541. 'TYPI' MTIV KTYPI 'CLIM' clim 'SMBR' sf
  542. 'IMPR' IMPKRES MATASS TYRENU PCMLAG
  543. XINIT DI2 DI3 DI4 DI5
  544. MAPREC KPREC PR3 PR4 PR5;
  545. res = ASSI numassi res1 et res ;
  546. Finsi;
  547. *........................................................
  548. Fin Bclres ;
  549.  
  550. rv.inco.'Res'=res ;
  551.  
  552. 'SI' (rv.'PROJ');
  553. **
  554. * mess 'RESOLUTION : ETAPE DE PROJECTION ' ;
  555. ** _n
  556. ** C U
  557. cun = ('KMF' (rv . 'MATC') res)*(1./(rv.'DELTAT')) ;
  558. Si (rv.'TVNP');
  559. cxn = ('KMF' (rv . 'MBTR') res)*(1./(rv.'DELTAT')) ;
  560. cun = cun et cxn ;
  561. Finsi ;
  562. clim=exco (rv.'CLIM') (mots (extr rv.'LISTINCO' 2)) 'NOID';
  563. MTIV=rv.'METHINV2' ;
  564.  
  565. *........................................................
  566. Si (rv.'XRIG') ;
  567. gclim=extr clim 'MAILLAGE' ;
  568. rclim=bloq 'T' gclim ;
  569. dclim=depi rclim (nomc 'Q' clim);
  570. matpr = matpr et rclim ;
  571. cun = cun et dclim ;
  572. sl = 'RESO' matpr cun ;
  573. Sinon ;
  574. MATASS= matpr ;
  575. TYRENU=MTIV.'TYRENU';PCMLAG=MTIV.'PCMLAG';KTYPI =MTIV.'TYPINV';
  576. XINIT=TEXT ' ' ;DI2 = TEXT ' ' ;DI3 = TEXT ' ' ;
  577. DI4 = TEXT ' ' ;DI5 = TEXT ' ' ;
  578. MAPREC=TEXT ' ';KPREC =TEXT ' ';PR3 = TEXT ' ' ;
  579. PR4 = TEXT ' ' ;PR5 = TEXT ' ' ;
  580. Si (('>EG' KTYPI 2) ET ('<EG' KTYPI 5));
  581. XINIT mar = assi numassi 'KOPS' 'MATRIK' ;
  582. DI2=MTIV.'NITMAX';DI3=MTIV.'RESID';DI4=MTIV.'BCGSBTOL';
  583. Si ('EGA' KTYPI 5);DI5=MTIV.'GMRESTRT';Finsi ;
  584. MAPREC= matpr;
  585. KPREC=MTIV.'PRECOND';
  586. Si ('EGA' KPREC 4);PR3=MTIV.'MILURELX';Finsi ;
  587. Si (('EGA' KPREC 5) OU ('EGA' KPREC 6));
  588. PR4=MTIV.'ILUTLFIL';PR5=MTIV.'ILUTDTOL';
  589. Finsi ;
  590. Finsi ;
  591. sl = ASSI numassi 'KRES' 'LL' matpr
  592. 'TYPI' MTIV KTYPI 'CLIM' clim 'SMBR' cun
  593. 'IMPR' IMPKRES MATASS TYRENU PCMLAG
  594. XINIT DI2 DI3 DI4 DI5
  595. MAPREC KPREC PR3 PR4 PR5;
  596. Finsi;
  597. *........................................................
  598. ctl = 'KMF' (rv . 'MATC') sl 'TRAN' ;
  599. Si(rv.'TVNP') ;
  600. cxl = 'KMF' (rv . 'MBTR') sl 'TRAN' ;
  601. ctl = ctl + cxl ;
  602. Finsi ;
  603. a=KOPS IDigvl '*' ctl ;
  604. au = nomc (mots 'UX' 'UY' 'UZ')(mots cnvi1 cnvi2 cnvi3) a ;
  605. oublier ctl ;
  606. res = res - (au*(rv.'DELTAT')) ;
  607. 'FINSI' ;
  608.  
  609.  
  610. *************************--O--***********************
  611. * \ /|\
  612. * /\____ RELAXATION
  613. * / \ \
  614. *****************************************************
  615. nbic=dime rv.'LISTINCO';
  616. omega=rv.'OMEGA';
  617. NUPDT=rv.'PASDETPS'.'NUPASDT';
  618. LRA= CHAI 'Pas de temps ' NUPDT ' ERREUR LINF : ' ;
  619. 'SI' ('MULT' NUPDT (rv.'FIDT'));
  620. Si('EXIST' rv.'PASDETPS' 'NUPDT');
  621. rv.'PASDETPS'.'NUPDT'=rv.'PASDETPS'.'NUPDT' et (PROG NUPDT) ;
  622. Sinon ;
  623. rv.'PASDETPS'.'NUPDT'= (PROG NUPDT) ;
  624. Finsi ;
  625. 'FINSI' ;
  626. Repeter Relax nbic;
  627.  
  628. a= chai 'FORMAT' '(A3,I1)' 'INC' &relax ;
  629.  
  630. Si (XBDF2) ;
  631. b= chai 'FORMAT' '(A3,I1)' 'IMC' &relax ;
  632. rv.'INCO'.(rv.b)=rv.'INCO'.(rv.a);
  633. * mess ' On fait ' (rv.a) ' -> ' (rv.b) ;
  634. Finsi ;
  635.  
  636. nc =extr rv.'LISTINCO' &relax;
  637. nc1=mots nc ;
  638. nc2=mots 'SCAL' ;
  639. ty=extr rv.'TYPEINCO' &relax;
  640. * mess ' Composante ' nc ' Type inco ' ty ;
  641. Si (('EGA' ty 'VITESSE ') et ('EGA' ('VALEUR' 'DIME') 2));
  642. nc2=mots 'UX' 'UY' ;
  643. a1= chai 'FORMAT' '(I1,A3)' 1 nc ;
  644. a2= chai 'FORMAT' '(I1,A3)' 2 nc ;
  645. nc1=mots a1 a2 ;
  646. Finsi ;
  647. Si (('EGA' ty 'VITESSE ') et ('EGA' ('VALEUR' 'DIME') 3));
  648. nc2=mots 'UX' 'UY' 'UZ' ;
  649. a1= chai 'FORMAT' '(I1,A3)' 1 nc ;
  650. a2= chai 'FORMAT' '(I1,A3)' 2 nc ;
  651. a3= chai 'FORMAT' '(I1,A3)' 3 nc ;
  652. nc1=mots a1 a2 a3 ;
  653. Finsi ;
  654. ri= exco nc1 res nc2 'NOID';
  655. nbc=dime (extr ri 'COMP') ;
  656. Si('EGA' nbc 0);
  657. Si('EGA' NUPDT 1);
  658. Mess ' ATTENTION La (ou les) composante(s) suivantes ne sont pas'
  659. ' calculees' ;
  660. list nc1 ;
  661. Finsi ;
  662. Sinon ;
  663. 'SI' ('MULT' NUPDT (rv.'FIDT'));
  664. eri=ri - rv.'INCO'.(rv.a);
  665. ELI=MAXI ERI 'ABS';
  666. ELI=(LOG (ELI + 1.0E-10))/(LOG 10.);
  667. LRA= CHAI LRA ' ' (rv.a) ' ' ELI ;
  668. ER = PROG ELI;
  669. ERA=CHAI 'ERLI' (rv.a);
  670. Si('EXIST' rv.'PASDETPS' ERA);
  671. rv.'PASDETPS'.ERA=rv.'PASDETPS'.ERA et ER ;
  672. Sinon ;
  673. rv.'PASDETPS'.ERA= ER ;
  674. Finsi ;
  675. 'FINSI' ;
  676.  
  677. rv.'INCO'.(rv.a)=(omega*ri) + ((1. - omega)*rv.'INCO'.(rv.a)) ;
  678. Finsi ;
  679.  
  680. Fin Relax ;
  681.  
  682. 'SI' ('MULT' NUPDT (rv.'FIDT'));
  683. MESS LRA ;
  684. 'FINSI' ;
  685.  
  686. Fin BCLINT ;
  687.  
  688. ***** Avancement en temps
  689. rv.inco.'Resnm'=res ;
  690. rv.'PASDETPS'.'NUPASDT'=rv.'PASDETPS'.'NUPASDT' + 1 ;
  691. * Algorithme de projection P(n+1) = P(n) + (sl / dt)
  692. 'SI' (rv.'PROJ');
  693.  
  694. 'SI' ('EXIST' (rv.'INCO') 'PINM1') ;
  695.  
  696. * PNM1=(-1.)*rv.'INCO'.'PRESSION' ;
  697. PNM1=(-1.)*rv.'INCO'.(rv.'INC2') ;
  698. si (XBDF2) ;
  699. PN=PNM1 + (1.5*sl) ;
  700. sinon ;
  701. PN=PNM1 + (sl) ;
  702. finsi ;
  703.  
  704. * rv.'INCO'.'PRESSION'=(-1.)*PN ;
  705. rv.'INCO'.(rv.'INC2')=(-1.)*PN ;
  706. rv.'INCO'.'PINM1'=PN ;
  707. 'SINON' ;
  708. * rv.'INCO'.'PRESSION' = (-1.)*(sl) ;
  709. rv.'INCO'.(rv.'INC2')= (-1.)*(sl) ;
  710. rv.'INCO'.'PINM1'=(2.*sl);
  711. 'FINSI' ;
  712. 'OUBLIER' sl ;
  713. 'FINSI' ;
  714. ***** Avancement en temps FIN
  715. Fin BCLTPS ;
  716.  
  717. QUITTER EXEC ;
  718. Finsi ;
  719. ********************************************************************************
  720. ********************************************************************************
  721. ****** Fin directive EQUA ******************************************************
  722. ********************************************************************************
  723. ********************************************************************************
  724.  
  725. 'SI' ('EGA' (rv . 'NAVISTOK') 0) ;
  726. EXAC rv ;
  727. 'QUITTER' EXEC ;
  728. 'FINSI' ;
  729.  
  730. *******PROJ*************************************************************
  731. *******PROJ Trois possibilités *****************************************
  732. *******PROJ TMDM1 (vrai)_> Correction Gresho NON
  733. *******PROJ TGRAD (vrai)-> Formulation Gradient NON
  734. *******PROJ TPNM2 (vrai)-> elimination end of step velocity OUI
  735. *******PROJ TMDM1 = VRAI TGRAD = FAUX TPNM2 = FAUX ancien algo
  736. TYPROJ='VPI1' ;
  737. 'SI' ('EXIST' RV 'TYPROJ') ;
  738. TYPROJ=RV.'TYPROJ' ;
  739. 'SI' (non (exist (MOTS 'PSCT' 'VPI1' 'VPI2') TYPROJ));
  740. Mess '*********************************************************' ;
  741. Mess ' ERREUR ERREUR ERREUR ERREUR ERREUR ERREUR ' ;
  742. Mess ' ' ;
  743. Mess 'Le mot ' TYPROJ ' n existe pas dans la liste.';
  744. Mess 'Les méthodes de projection autorisées sont :';
  745. Mess 'VPI1 et VPI2 ' ;
  746. Mess '*********************************************************' ;
  747. erreur 21 ;
  748. quitter EXEC ;
  749. 'FINSI' ;
  750. 'FINSI' ;
  751.  
  752. 'SI' ('EGA' TYPROJ 'PSCT') ;
  753. TGRAD = FAUX ;
  754. TMDM1 = FAUX ;
  755. TPNM2 = FAUX ;
  756. 'FINSI' ;
  757. 'SI' ('EGA' TYPROJ 'VPI1') ;
  758. TGRAD = FAUX ;
  759. TMDM1 = VRAI ;
  760. TPNM2 = FAUX ;
  761. 'FINSI' ;
  762. 'SI' ('EGA' TYPROJ 'VPI2') ;
  763. TGRAD = FAUX ;
  764. TMDM1 = FAUX ;
  765. TPNM2 = VRAI ;
  766. 'FINSI' ;
  767.  
  768. *******PROJ Trois possibilités *****************************************
  769. *******PROJ*************************************************************
  770.  
  771. nomvi=rv . 'NOMVI' ;
  772.  
  773. 'SI' ('EGA' ('VALEUR' 'DIME') 2) ;
  774. vnul=0.D0 0.D0 ;
  775. vuni=1.D0 1.D0 ;
  776. cnvi1= 'MOT' ('TEXTE' ('CHAINE' 1 nomvi)) ;
  777. cnvi2= 'MOT' ('TEXTE' ('CHAINE' 2 nomvi)) ;
  778. lc= mots cnvi1 cnvi2 ;
  779. lcu=mots 'UX' 'UY' ;
  780. 'SINON' ;
  781. vnul=0.D0 0.D0 0.D0 ;
  782. vuni=1.D0 1.D0 1.D0 ;
  783. cnvi1= 'MOT' ('TEXTE' ('CHAINE' 1 nomvi)) ;
  784. cnvi2= 'MOT' ('TEXTE' ('CHAINE' 2 nomvi)) ;
  785. cnvi3= 'MOT' ('TEXTE' ('CHAINE' 3 nomvi)) ;
  786. lc= mots cnvi1 cnvi2 cnvi3 ;
  787. lcu=mots 'UX' 'UY' 'UZ' ;
  788. 'FINSI' ;
  789. 'SI' ('NON' ('EXISTE' rv 'OMEGA')) ;
  790. omeg=1.D0 ;
  791. 'SINON' ;
  792. omeg=rv . 'OMEGA' ;
  793. 'FINSI' ;
  794.  
  795. testpr ='EXISTE' rv 'PRESSION' ;
  796. testprj ='EXISTE' rv 'PROJ' ;
  797. testran=testpr 'ET' ('EXISTE' rv 'CO') ;
  798.  
  799. si testpr ; rvp = rv.'PRESSION' ;
  800. achp matpr= 'KOPS' 'MATRIK' ;
  801. rvp.'MATP'= matpr;
  802. finsi ;
  803. si testprj; rvp = rv.'PROJ' ; Finsi ;
  804.  
  805. 'SI' (rv.'IMPR' >EG 1) ;
  806. si testpr ;
  807. mess 'Algorithme semi explicite (ANCIEN)';
  808. mess '==================================';
  809. Finsi ;
  810. si testprj;
  811. mess 'Algorithme de Projection';
  812. mess '========================';
  813. finsi ;
  814. si ( non (testpr ou testprj)) ;
  815. mess 'Algorithme standard implicite ou explicite ';
  816. mess '===========================================';
  817. finsi ;
  818. 'FINSI' ;
  819.  
  820. 'SI' ('NON' ('EXISTE' rv 'HIST')) ;
  821. rv . 'HIST' = 'TABLE' ;
  822. 'FINSI' ;
  823.  
  824.  
  825. ITMA=(rv . 'ITMA') ;
  826. IMPTCRR=0 ;
  827. IMPKRES=0 ;
  828. 'SI' ('<EG' ITMA 1) ;
  829. ITMA=1 ;
  830. IMPTCRR=1 ;
  831. 'FINSI' ;
  832. Si (testprj ) ; IMPTCRR=1 ; finsi ;
  833. Si ( non (testpr ou testprj)) ; IMPTCRR=1 ; finsi ;
  834.  
  835.  
  836. * Gestion de la matrice de préconditionnement
  837. *
  838. * calass : doit-on recalculer l'assemblage
  839. * calprec : doit-on recalculer le préconditionneur
  840. * calassp : idem pour la matrice de pression
  841. * calprecp: idem pour la matrice de pression
  842. * fcprect : fréquence de recalcul du préc. en fn du pas de temps
  843. * fcpreci : fréquence de recalcul du préc. dans la boucle
  844. * d'itérations internes pour les non-linéarités
  845. * resmn : le résidu au pas de temps ou à l'itération interne précédente
  846. *
  847. calprec = VRAI ;
  848. fcprect = rv. 'METHINV' . 'FCPRECT' ;
  849. fcpreci = rv. 'METHINV' . 'FCPRECI' ;
  850. resmn maprec = 'KOPS' 'MATRIK' ;
  851. resmn1 maprec1 = 'KOPS' 'MATRIK' ;
  852. resmn2 maprec2 = 'KOPS' 'MATRIK' ;
  853. 'SI'(exist rv 'resmn');resmn=rv.'resmn' ; finsi ;
  854. 'SI'(non (exist rv 'calprec'));rv.'calprec'=VRAI;finsi;
  855. 'SI'(non (exist rv 'calass' ));rv.'calass' =VRAI;calass =VRAI;finsi;
  856. 'REPETER' bloc1 ITMA ;
  857. * mess ' DEBUT BLOC1 ';
  858.  
  859. testp1 = FAUX ;
  860. si (testpr ou testprj) ;
  861. Si CALPRE ; testp1=VRAI ; Finsi ;
  862. Si (non (exist rvp 'MATC')) ; testp1=VRAI ; Finsi ;
  863. finsi ;
  864.  
  865. 'SI' ('MULT' (&bloc1 '+' 1) fcprect) ;
  866. rv.'calprec' = VRAI ;
  867. 'FINSI' ;
  868. 'REPETER' bloci (rv . 'NITER') ;
  869. 'SI' ('MULT' (&bloci '+' 1) fcpreci) ;
  870. rv.'calprec' = VRAI ;
  871. 'FINSI' ;
  872.  
  873. 'SI' (testprj) ;
  874. calprecp=VRAI ;
  875. calassp =VRAI ;
  876. si (exist rvp 'calprecp' ) ; calprecp=rvp.'calprecp';
  877. maprec1=rvp.'maprec1' ;finsi ;
  878. si (exist rvp 'calassp' ) ; calassp=rvp.'calassp';
  879. matass1=rvp.'matass1' ;finsi ;
  880. fcprectp = rvp . 'METHINV' . 'FCPRECT' ;
  881. 'SI' ('MULT' (&bloc1 '+' 1) fcprectp) ;
  882. calprecp= VRAI ;
  883. 'FINSI' ;
  884. 'FINSI' ;
  885.  
  886.  
  887.  
  888. st mat = 'KOPS' 'MATRIK' ;
  889. sf mau = 'KOPS' 'MATRIK' ;
  890.  
  891.  
  892. mdfdt = 0 ;
  893. 'REPETER' bloc2 ('DIME' (rv . 'LISTOPER')) ;
  894. nomper = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ;
  895. notable= 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ;
  896. * mess 'Operateur Bloc2 mdfdt ? ' nomper ;
  897. mdfdt = mdfdt + rv. notable . 'KOPT' . 'KFORM' ;
  898.  
  899. si (ega nomper 'DFDT ');
  900. ISCHT= (rv . notable.kopt.'ISCHT') ;
  901. msi mai=('TEXTE' nomper) (rv . notable) ;
  902. mat = mat 'ET' mai ;
  903. st = st 'ET' msi ;
  904.  
  905. sinon ;
  906.  
  907. * mess ' 1OPER =' nomper;
  908. msi mai=('TEXTE' nomper) (rv . notable) ;
  909. mau = mau 'ET' mai ;
  910. sf = sf 'ET' msi ;
  911.  
  912. finsi ;
  913.  
  914. 'FIN' bloc2 ;
  915.  
  916. s2 = sf et st ;
  917. ma1 = mau 'ET' mat ;
  918.  
  919. ********** Traitement PRESSION ****************************
  920. *mess 'Traitement PRESSION';
  921.  
  922. 'SI' (exist rv 'rvpd') ;
  923. rvpd = rv.'rvpd' ;
  924. IDigv=rv.'IDigv';
  925. Digv=rv.'Digv';
  926. matpr= rvp.'MATP' ;
  927. 'FINSI' ;
  928.  
  929. */1 CALCUL de CMCT pour PRESSION
  930. 'SI' (testpr et testp1) ;
  931.  
  932. rvpd = (rvp.'DOMAINE') ;
  933. Diago = doma rvpd 'XXDIAGSI' ;
  934. si ( ega ('VALEUR' 'DIME') 2) ;
  935. Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) ;
  936. sinon ;
  937. Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) et
  938. ( exco Diago 'SCAL' cnvi3 ) ;
  939. finsi ;
  940. Digv = kcht rvpd vect sommet comp lc Digv ;
  941.  
  942. 'SI' ('EXISTE' rv 'CLIM') ;
  943. rvp . 'CLIM' = rv . 'CLIM' ;
  944. Digv = 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1;
  945. 'FINSI' ;
  946.  
  947. IDigv= 'INVERSE' digv ;
  948.  
  949. 'SI' ('NON' ('EXISTE' rvp 'DIAGV')) ;
  950. rvp . 'DIAGV' = digv ;
  951. 'FINSI' ;
  952.  
  953. rvp . 'MATC' = 'KMAB' rvp ;
  954. rvp . 'PRESSION' = 'KCHT' rvpd 'SCAL'
  955. 'CENTRE' 0.D0 ;
  956. rvp . 'GRADP' = 'KCHT' rvpd 'VECT'
  957. 'SOMMET' vnul ;
  958.  
  959. achp matpr= 'KOPS' 'MATRIK' ;
  960. rvp . 'MATP' = matpr ;
  961. rv.'rvpd'=rvpd;
  962. rv.'IDigv'=IDigv;
  963. rv.'Digv'=Digv;
  964.  
  965. 'FINSI' ;
  966.  
  967. */3 CALCUL de CMCT pour PROJ
  968. 'SI' (testprj et testp1) ;
  969.  
  970. Diago mma = 'KOPS' 'MATRIK' ;
  971.  
  972. idfdt = 0 ;
  973. 'REPETER' blocj ('DIME' (rv . 'LISTOPER')) ;
  974. nomper = 'EXTRAIRE' &blocj (rv . 'LISTOPER') ;
  975. notable= 'MOT' ('TEXTE' ('CHAINE' &blocj nomper)) ;
  976. *? mess 'Operateur ' nomper ;
  977. si (ega nomper 'DFDT ');
  978. idfdt=idfdt + 1 ;
  979. Diago = Diago et (doma (rv. notable . 'DOMZ') 'XXDIAGSI') ;
  980. rvpd=rv. notable . 'DOMZ';
  981. finsi ;
  982. 'FIN' blocj ;
  983. si (ega idfdt 0) ; mess ' Pas de DFDT ?? ' ; erreur 21 ; finsi ;
  984.  
  985. si ( ega ('VALEUR' 'DIME') 2) ;
  986. Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) ;
  987. sinon ;
  988. Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) et
  989. ( exco Diago 'SCAL' cnvi3 ) ;
  990. finsi ;
  991. Digv = kcht rvpd vect sommet comp lc Digv ;
  992. * Matrice diagonale ne contenant pas les conditions aux limites : Diag
  993. Diag = copier Digv;
  994. IDiag = 'INVERSE' Diag ;
  995.  
  996. Digv = 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1;
  997.  
  998. IDigv= 'INVERSE' digv ;
  999.  
  1000. 'SI' ('NON' ('EXISTE' rvp 'DIAGV')) ;
  1001. rvp . 'DIAGV' = digv ;
  1002. 'FINSI' ;
  1003.  
  1004. rvp.'INCO'=rv.'INCO' ;
  1005.  
  1006. sp map = 'KOPS' 'MATRIK' ;
  1007. sr mar = 'KOPS' 'MATRIK' ;
  1008. svnpc mvnpc = 'KOPS' 'MATRIK' ;
  1009. 'REPETER' blocpj ('DIME' (rvp . 'LISTOPER')) ;
  1010. nomper = 'EXTRAIRE' &blocpj (rvp . 'LISTOPER') ;
  1011. notable= 'MOT' ('TEXTE' ('CHAINE' &blocpj nomper)) ;
  1012. mess 'Procedure PROJ Operateur : ' nomper ;
  1013.  
  1014. si (EGA nomper 'KBBT');
  1015. rvp . notable . 'KOPT' . 'IKOMP' = 1 ;
  1016. finsi ;
  1017.  
  1018. msi mai=('TEXTE' nomper) (rvp . notable) ;
  1019.  
  1020. *******PROJ*************************************************************
  1021. *******PROJ TGRAD calcul maig (Formulation en Gradient) DEBUT*** *******
  1022. Si TGRAD ;
  1023. mess 'Procedure PROJCT (Gradient) Operateur : ' nomper ;
  1024. rvp . notable . 'KOPT' . 'IKOMP' = 0 ;
  1025. msig maig=('TEXTE' nomper) (rvp . notable) ;
  1026. Finsi ;
  1027. *******PROJ -> MATG FIN*****
  1028. *******PROJ/////////////////////////////////////////////////////////////
  1029.  
  1030. TVNP=FAUX;
  1031. TVNPC=FAUX;
  1032. Si (('EGA' nomper 'VNIMP ') et (non (ega (rvp.'DISCPRES') 5)));
  1033. TVNP=VRAI ;
  1034. mar = mar 'ET' mai ;
  1035. sr = sr 'ET' msi ;
  1036. Finsi ;
  1037.  
  1038. * mess 'nomper=' nomper (rvp.'DISCPRES');
  1039. Si (('EGA' nomper 'VNIMP ') et (ega (rvp.'DISCPRES') 5));
  1040. TVNPC=VRAI;
  1041. mvnpc=mvnpc et mai ;
  1042. svnpc=svnpc et msi ;
  1043. Finsi ;
  1044.  
  1045. Si (non ('EGA' nomper 'VNIMP '));
  1046. map = map 'ET' mai ;
  1047. sp = sp 'ET' msi ;
  1048. Finsi ;
  1049.  
  1050. 'FIN' blocpj ;
  1051.  
  1052. matpc mac = kops 'CMCTSPLT' map ;
  1053. *******PROJ*************************************************************
  1054. *******PROJ TGRAD (vrai Formulation en Gradient) CMCTSPLT DEBUT*** *******
  1055. Si TGRAD ;
  1056. matpcg macg = kops 'CMCTSPLT' maig ;
  1057. Finsi ;
  1058. *******PROJ -> MATG FIN*****
  1059. *******PROJ/////////////////////////////////////////////////////////////
  1060.  
  1061. Si TVNPC; mac=mac et mvnpc; Finsi;
  1062. rvp . 'MATC' = mac ;
  1063.  
  1064. *******PROJ*************************************************************
  1065. *******PROJ TGRAD (vrai Formulation en Gradient) -> MATG DEBUT*** *******
  1066. SI TGRAD ; rvp . 'MATG' = macg; FINSI ;
  1067. *******PROJ -> MATG FIN*****
  1068. *******PROJ/////////////////////////////////////////////////////////////
  1069.  
  1070. scr mcr = 'KOPS' 'MATRIK' ;
  1071. Si TVNP ;
  1072. rvp . 'MBTR' = mar ;
  1073. Dunit=Idigv;
  1074. crt= kops 'CMCT' mac mar (Dunit) ;
  1075. ctr= kops 'CMCT' mar mac (Dunit) ;
  1076. rrt= kops 'CMCT' mar mar (Dunit) ;
  1077. mcr= mcr et crt et ctr et rrt ;
  1078. Finsi ;
  1079. rvp . 'TVNP' = TVNP ;
  1080.  
  1081. Si( ega (rvp.'DISCPRES') 5 ) ;
  1082. mess ' Cas des pressions continues' ;
  1083. matpr = matpc et mcr ;
  1084. Sinon ;
  1085. mess ' Cas des pressions discontinues' ;
  1086.  
  1087. matpr = kops 'CMCT' mac mac IDigv ;
  1088. matpr = matpc et matpr et mcr ;
  1089.  
  1090. Finsi ;
  1091.  
  1092.  
  1093. rvp . 'MATP' = matpr ;
  1094. rv.'rvpd'=rvpd;
  1095. rv.'IDigv'=IDigv;
  1096. rv.'Digv'=Digv;
  1097. 'FINSI' ;
  1098. *************
  1099. ** Fin calcul CMCT pour testp1 et testprj
  1100. ******************************************
  1101. * t
  1102. *************** Calcul de C p Pression
  1103. 'SI' testpr ;
  1104.  
  1105. dt = (rv . 'PASDETPS' . 'DELTAT') '*' (rv . 'ALFA') ;
  1106. rvp . 'DELTAT' = dt ;
  1107. f = 'COPIER' s2 ;
  1108. u = rv . 'INCO' . nomvi ;
  1109.  
  1110. lc = 'EXTRAIRE' digv 'COMP' ;
  1111. fu = 'KCHT' rvpd 'VECT' 'SOMMET'
  1112. 'COMP' lc ('EXCO' f lc) ;
  1113.  
  1114. 'SI' ('EXISTE' rv 'CLIM') ;
  1115. dti = -1.D0 '/' dt ;
  1116. dm1f= 'KOPS' (dt '*' ('KOPS' fu '*' IDigv))
  1117. 'CLIM' (rv . 'CLIM') 3 ;
  1118. dm1f= dti * dm1f ;
  1119. 'SINON' ;
  1120. dm1f= (-1.D0) '*' ('KOPS' fu '*' IDigv) ;
  1121. 'FINSI' ;
  1122.  
  1123. rvp . 'PRESSION' = 'KMF' (rvp . 'MATC') dm1f ;
  1124.  
  1125. 'KRES' rvp (rvp . 'PRESSION')
  1126. 'BETA' (rvp . 'KBETA') (rvp . 'BETA')
  1127. 'PIMP' (rvp . 'KPIMP') (rvp . 'PIMP') ;
  1128. rvp . 'GRADP' = 'KMTP' 1 (rvp . 'MATC')
  1129. (rvp . 'PRESSION') lc ;
  1130. s2 = s2 + (rvp . 'GRADP') ;
  1131. rv.'INCO'.'PRESSION'=rvp . 'PRESSION';
  1132.  
  1133. 'FINSI' ;
  1134.  
  1135. *******PROJ*************************************************************
  1136. *******PROJ MDM1 Ici on calcule M D-1 pour ensuite calculer Ctp DEBUT*** *******
  1137. *******PROJ on ne calcule M D-1 que si (rv 'MDM1') n'existe pas DEBUT*** *******
  1138. *******PROJ ou CALPRE VRAI DEBUT*** *******
  1139. 'SI' testprj ;
  1140.  
  1141. dt = (rv . 'PASDETPS' . 'DELTAT') '*' (rv . 'ALFA') ;
  1142. rvp . 'DELTAT' = dt ;
  1143.  
  1144. ** Produit M D-1
  1145. 'SI' TMDM1 ;
  1146. 'SI' (('EXIST' rv 'MDM1') et (NON CALPRE)) ;
  1147. MDM1=rv.'MDM1' ;
  1148. 'SINON' ;
  1149. mess ' On calcule MD-1 ' ;
  1150. stn matn= 'KOPS' 'MATRIK' ;
  1151.  
  1152. 'REPETER' blocpj1 ('DIME' (rvp . 'LISTOPER')) ;
  1153. nomper1 = 'EXTRAIRE' &blocpj1 (rvp . 'LISTOPER') ;
  1154. notable1= 'MOT' ('TEXTE' ('CHAINE' &blocpj1 nomper1)) ;
  1155. * mess 'Operateur ' nomper1 ;
  1156.  
  1157. si (ega nomper1 'KBBT ');
  1158. nomiv1= extr (rvp . notable1 . 'LISTINCO') 1 ;
  1159.  
  1160. 'REPETER' blocpj2 ('DIME' (rv . 'LISTOPER')) ;
  1161. nomper2 = 'EXTRAIRE' &blocpj2 (rv . 'LISTOPER') ;
  1162. notable2= 'MOT' ('TEXTE' ('CHAINE' &blocpj2 nomper2)) ;
  1163. * mess 'Operateur ' nomper2 ;
  1164.  
  1165. si (ega nomper2 'DFDT ');
  1166. nomiv2= extr (rv . notable2 . 'LISTINCO') 1 ;
  1167. si ('EGA' nomiv2 nomiv1) ;
  1168.  
  1169. domzp=rv . notable2 . 'DOMZ' ;
  1170. msi mai='DFDT' (rv . notable2) ;
  1171. matn = matn et mai ;
  1172. sti = kcht domzp vect sommet comp lc vuni;
  1173. finsi ;
  1174. finsi ;
  1175. 'FIN' blocpj2 ;
  1176. finsi ;
  1177. 'FIN' blocpj1 ;
  1178.  
  1179. MDM1='KMF' matn IDiag ;
  1180.  
  1181. Si(EGA ISCHT 1);
  1182. MDM1=MDM1 * ((dt*2.)/3.) ;
  1183. Sinon ;
  1184. MDM1=MDM1 * dt ;
  1185. Finsi ;
  1186.  
  1187. rv.'MDM1' = MDM1 ;
  1188.  
  1189. *******PROJ MDM1 -> rv.'MDM1' FIN*****
  1190. *******PROJ/////////////////////////////////////////////////////////////
  1191. 'FINSI';
  1192. 'FINSI';
  1193. *******PROJ*************************************************************
  1194. *******PROJ t n-1 DEBUT*** *******
  1195. *******PROJ Calcul de C P DEBUT*** *******
  1196.  
  1197. TVNP=rvp.'TVNP' ;
  1198.  
  1199. 'SI' ('EXIST' (rv.'INCO') 'PRESSION') ;
  1200. PPI = rv.'INCO'.'PRESSION' ;
  1201. 'SI' ('EXIST' (rv.'INCO') 'PNM2') ;
  1202. * TPNM2 elimination of end of step velocity (2 Pn - Pn-1)
  1203. PPI = 2*(PPI) - rv.'INCO'.'PNM2' ;
  1204. 'FINSI' ;
  1205.  
  1206. *Formulation en u grad p
  1207. Si (Exist rvp 'MATG');
  1208. cpre = 'KMF' (rvp . 'MATG') PPI 'TRAN' ;
  1209. *Formulation en p div u
  1210. Sinon;
  1211. cpre = 'KMF' (rvp . 'MATC') PPI 'TRAN' ;
  1212. Finsi ;
  1213.  
  1214. Si TVNP ;
  1215. cxre = 'KMF' (rvp . 'MBTR') PPI 'TRAN' ;
  1216. cpre = cpre et cxre ;
  1217. Finsi ;
  1218.  
  1219. *******PROJ MDM1
  1220. Si TMDM1 ;
  1221. * consistence selon (Gresho)
  1222. gradpres = 'KOPS' MDM1 '*' cpre ;
  1223. Sinon ;
  1224. * consistence selon (Guermond)
  1225. gradpres = cpre ;
  1226. Finsi ;
  1227.  
  1228. oublier cpre ;
  1229. rv.'INCO'.'GRADPRES' = 'NOMC' lcu lc gradpres;
  1230.  
  1231. *Formulation en u grad p
  1232. Si (Exist rvp 'MATG');
  1233. s2 = sf - (rv.'INCO'.'GRADPRES') + st ;
  1234. Sinon ;
  1235. *Formulation en p div u
  1236. s2 = sf + (rv.'INCO'.'GRADPRES') + st ;
  1237. Finsi ;
  1238.  
  1239. 'FINSI' ;
  1240.  
  1241. *******PROJ t n-1 FIN*****
  1242. *******PROJ MDM1 '*' C P -> rv.'INCO'.'GRADPRES' FIN*****
  1243. *******PROJ t n-1 FIN*****
  1244. *******PROJ si (p div u) -> S2 = F + C P + st FIN*****
  1245. *******PROJ t n-1 FIN*****
  1246. *******PROJ si (u grad p) -> S2 = F - C P + st FIN*****
  1247. *******PROJ/////////////////////////////////////////////////////////////
  1248. 'FINSI' ;
  1249.  
  1250. ********** FIN Traitement PRESSION ************************
  1251. * mess ' FIN Traitement PRESSION' ;
  1252.  
  1253. 'SI' ('EXISTE' rv 'CLIM') ;
  1254. s1 = rv . 'CLIM' ;
  1255. 'SINON' ;
  1256. s1=' ' ;
  1257. 'FINSI' ;
  1258. rv . 'S2' = s2 ;
  1259.  
  1260. *******************************************************************
  1261. * Résolution hors QDM méthode de projection
  1262. 'SI'((NON testprj) ou (EGA mdfdt 0));
  1263. *mess '* Résolution hors QDM méthode de projection ';
  1264.  
  1265. * mess ' ====================> calprec '; list calprec ;
  1266. 'SI' rv.'calass' ;
  1267. rv . 'METHINV' . 'MATASS' =ma1 ;
  1268. rv.'calass'=FAUX ;
  1269. 'FINSI' ;
  1270.  
  1271. 'SI' rv.'calprec' ;
  1272. rv . 'METHINV' . 'MAPREC' = ma1 ;
  1273. rv.'calprec'=FAUX ;
  1274. 'FINSI' ;
  1275.  
  1276. rv . 'METHINV' . 'XINIT' = resmn ;
  1277.  
  1278. res = 'KRES' ma1 'TYPI' (rv . 'METHINV')
  1279. 'CLIM' s1
  1280. 'SMBR' s2
  1281. 'IMPR' IMPKRES ;
  1282.  
  1283. 'FINSI' ;
  1284.  
  1285. 'SI' (testprj et (non(EGA mdfdt 0)));
  1286. *******PROJ*************************************************************
  1287. *******PROJ Résolution QDM DEBUT*** *******
  1288. *******PROJ (On éclate les résolutions) DEBUT*** *******
  1289.  
  1290. * Résolution QDM méthode de projection (On éclate les résolutions)
  1291. *mess ' Résolution QDM méthode de projection';
  1292. res mm1= kops 'MATRIK' ;
  1293. liscom=extr s1 'COMP';
  1294. nbcom=dime liscom ;
  1295. * MEss ' NBCOM=' nbcom ;
  1296.  
  1297. Si(NON(EXIST RV 'TABRES')); RV.'TABRES'=TABLE ; Finsi ;
  1298. Si(NON(EGA ('TYPE' RV.'TABRES') 'TABLE')); RV.'TABRES'=TABLE;Finsi;
  1299.  
  1300. repeter Bclcom nbcom ;
  1301. nmc=extr liscom &Bclcom;
  1302. * Mess ' Composante : ' nmc ;
  1303. ma1i=kops 'EXTRCOMP' nmc 'SCAL' ma1;
  1304. s1i = exco s1 nmc 'NOID' ;
  1305. s2i = exco s2 nmc 'NOID' ;
  1306. resmni = exco resmn nmc 'NOID' ;
  1307.  
  1308. MATASSI =chai 'FORMAT' '(A6,I1)' 'MATASS' &Bclcom;
  1309. MAPRECI =chai 'FORMAT' '(A6,I1)' 'MAPREC' &Bclcom;
  1310.  
  1311. 'SI' (rv.'calass' ou (NON ('EXIST' (RV.'TABRES') MATASSI))) ;
  1312. rv . 'TABRES' . MATASSI = ma1i;
  1313. rv.'calass'=FAUX ;
  1314. 'FINSI' ;
  1315.  
  1316. 'SI' (rv.'calprec' ou (NON ('EXIST' (RV.'TABRES') MAPRECI)));
  1317. mess 'On recalcule le preconditionneur ';
  1318. rv . 'TABRES' . MAPRECI = ma1i;
  1319. rv.'calprec'=FAUX ;
  1320. 'FINSI' ;
  1321.  
  1322. rv . 'METHINV' . 'MATASS' = rv . 'TABRES' . MATASSI ;
  1323. rv . 'METHINV' . 'MAPREC' = rv . 'TABRES' . MAPRECI ;
  1324. rv . 'METHINV' . 'XINIT' = resmni ;
  1325.  
  1326.  
  1327. resi = 'KRES' ma1i 'TYPI' (rv . 'METHINV')
  1328. 'CLIM' s1i
  1329. 'SMBR' s2i
  1330. 'IMPR' IMPKRES ;
  1331. resi = nomc resi nmc ;
  1332. res=res et resi ;
  1333. Fin Bclcom ;
  1334. oubli ma1 ;
  1335. *mess ' Fin Résolution QDM méthode de projection';
  1336. * Fin Résolution QDM méthode de projection
  1337. *******PROJ Résolution QDM -> res FIN*****
  1338. *******PROJ (On éclate les résolutions) FIN*****
  1339. *******PROJ/////////////////////////////////////////////////////////////
  1340. 'FINSI' ;
  1341.  
  1342. 'SI' testprj ;
  1343. *******PROJ*************************************************************
  1344. *******PROJ ETAPE DE PROJECTION DEBUT*** *******
  1345. *******PROJ on calcul cun (alias c U tilde) DEBUT*** *******
  1346. * mess 'ETAPE DE PROJECTION ' ;
  1347. * _n
  1348. * C U
  1349. cun = ('KMF' (rvp . 'MATC') res) ;
  1350. Si TVNP ;
  1351. cxn = ('KMF' (rvp . 'MBTR') res) ;
  1352. cun = cun et cxn ;
  1353. Finsi ;
  1354.  
  1355. * mess ' On calcule les seconds membres de l équation de pression ' ;
  1356. * mess ' s ils existent (opérateurs FIMP) ' ;
  1357.  
  1358. 'REPETER' blocpj ('DIME' (rvp . 'LISTOPER')) ;
  1359. nomper = 'EXTRAIRE' &blocpj (rvp . 'LISTOPER') ;
  1360. notable= 'MOT' ('TEXTE' ('CHAINE' &blocpj nomper)) ;
  1361. si (ega nomper 'FIMP ');
  1362. * mess 'Seconds membres de l équation de pression Opérateur ' nomper ;
  1363. msi mai=('TEXTE' nomper) (rvp . notable) ;
  1364. cun = cun et msi ;
  1365. finsi ;
  1366. 'FIN' blocpj ;
  1367.  
  1368. cun = cun * (-1./dt) ;
  1369.  
  1370. * mess ' ====================> calprecp '; list calprecp ;
  1371. 'SI' calassp;
  1372. matass1=matpr ;
  1373. rvp.'matass1'=matass1 ;
  1374. 'FINSI' ;
  1375.  
  1376. 'SI' calprecp ;
  1377. 'SI' (RV.IMPR >EG 1) ;
  1378. 'MESS' ' On preconditionne CMCT ' ;
  1379. 'FINSI' ;
  1380. maprec1 = matpr ;
  1381. rvp.'maprec1'=maprec1 ;
  1382. 'FINSI' ;
  1383. rvp . 'METHINV' . 'MATASS' = matass1 ;
  1384. rvp . 'METHINV' . 'MAPREC' = maprec1 ;
  1385. rvp . 'METHINV' . 'XINIT' = resmn2 ;
  1386.  
  1387. sl= 'KRES' matpr 'TYPI' (rvp . 'METHINV')
  1388. 'CLIM' (rvp . CLIM )
  1389. 'SMBR' cun
  1390. 'IMPR' IMPKRES ; oublier cun ;
  1391.  
  1392. oublier resmn2 ; resmn2 = sl ;
  1393. ctl = 'KMF' (rvp . 'MATC') sl 'TRAN' ;
  1394. Si TVNP ;
  1395. cxl = 'KMF' (rvp . 'MBTR') sl 'TRAN' ;
  1396. ctl = ctl + cxl ;
  1397. Finsi ;
  1398.  
  1399. * elimination of end of step velocity (Si exist PNM2) (on ne corrige pas)
  1400. 'SI'(NON(TPNM2));
  1401. 'SI' (NON ('EXIST' (rv.'INCO') 'PNM2')) ;
  1402. a=nomc lcu lc ('KOPS' IDigv '*' ctl) ;
  1403. res = res + (a*dt) ;
  1404. 'FINSI' ;
  1405. 'FINSI' ;
  1406.  
  1407. oublier ctl ;
  1408.  
  1409. *******PROJ -> (alias c U tilde) FIN*****
  1410. *******PROJ -> sl (alias lambda) FIN*****
  1411. *******PROJ -> on corrige u FIN*****
  1412. *******PROJ/////////////////////////////////////////////////////////////
  1413. 'FINSI' ;
  1414.  
  1415. ***** Avancement en temps
  1416. *? Si (testprj et ('EGA' mdfdt 0)) ; IMPTCRR=RV.'IMPR' ; finsi ;
  1417. Si ('EGA' mdfdt 0) ; IMPTCRR=RV.'IMPR' ; finsi ;
  1418. IMPKRES=0 ;
  1419. * mess 'On passe dans TCRR ' ;
  1420. eps = 'TCRR' res omeg (rv . 'INCO') 'IMPR' IMPTCRR ;
  1421.  
  1422. 'OUBLIER' resmn ;
  1423. resmn = res ;
  1424. 'MENAGE' ;
  1425. 'FIN' bloci ;
  1426. irt=0 ;
  1427.  
  1428. * mess ' ITMA= ' (rv . 'ITMA') ' mdfdt= ' mdfdt ;
  1429. 'SI' ('EGA' (rv . 'ITMA') 0) ;
  1430. * mess ' ON PASSE PLUS DS TCNM>>>>>>>>>>>>>>' ;
  1431. irt = 'TCNM' rv 'NOUP';
  1432. 'SINON' ;
  1433. * mess ' ON PASSE DS TCNM>>>>' (rv . 'ITMA') mdfdt;
  1434. irt = 'TCNM' rv ;
  1435. 'FINSI' ;
  1436.  
  1437. ***** Avancement en temps algorithme de projection
  1438.  
  1439. 'SI' (testprj) ;
  1440. 'SI' calprecp;
  1441. calprecp= FAUX ;
  1442. 'SINON' ;
  1443. 'SI' ('NON' calassp) ;
  1444. 'OUBLIER' matpr ;
  1445. 'FINSI' ;
  1446. 'FINSI' ;
  1447. calassp=FAUX ;
  1448. rvp.'calassp'=calassp ;
  1449. rvp.'calprecp'=calprecp ;
  1450. 'FINSI' ;
  1451.  
  1452. *******PROJ*************************************************************
  1453. *******PROJ P(n+1) = P(n) + sl DEBUT*** *******
  1454. * Algorithme de projection P(n+1) = P(n) + sl
  1455. 'SI' testprj ;
  1456.  
  1457. 'SI' (NON ('EXIST' (rv.'INCO') 'PRESSION'));
  1458. rv.'INCO'.'PRESSION' = sl;
  1459. Si TPNM2 ; rv.'INCO'.'PNM2' = sl ; Finsi ;
  1460. 'SINON' ;
  1461. PNM1=rv.'INCO'.'PRESSION';
  1462. si ( EGA ISCHT 1) ;
  1463. PN=PNM1 + (1.5*sl) ;
  1464. sinon ;
  1465. PN=PNM1 + sl ;
  1466. finsi ;
  1467.  
  1468. Si TPNM2 ; rv.'INCO'.'PNM2'=PNM1 ; Finsi ;
  1469. rv.'INCO'.'PRESSION'=PN ;
  1470. 'FINSI' ;
  1471.  
  1472. Si (EGA TYPROJ 'PSCT') ;
  1473. rv.'INCO'.'PRESSION' = sl;
  1474. Finsi ;
  1475.  
  1476. 'OUBLIER' sl ;
  1477. 'FINSI' ;
  1478. *******PROJ P(n+1) = P(n) + sl FIN*****
  1479. *******PROJ/////////////////////////////////////////////////////////////
  1480. ***** Avancement en temps FIN
  1481.  
  1482.  
  1483. 'SI' testran ;
  1484. rv . 'CO' . 'VITESSE' ='KOPS' (rv . 'INCO' . nomvi)
  1485. '-' (rv . 'SEDIM') ;
  1486. k = 'ABS' (rv . 'INCO' . 'KN') ;
  1487. e = 'ABS' (rv . 'INCO' . 'EN') ;
  1488. k = 'KOPS' ('KOPS' k '*' k) '/' e ;
  1489. dif = 'KOPS' ('KOPS' k '*' 0.09) '+' (rv . 'COEF') ;
  1490. rv . 'CO' . 'DIFFU' = 'NOEL' rvpd dif ;
  1491. rv . 'CO' . 'TEMPERA' = rv . 'INCO' . 'CN' ;
  1492. 'FINSI' ;
  1493. 'MENAGE' ;
  1494.  
  1495. 'SI' ('EGA' irt 1) ;
  1496. 'MESSAGE' ' Temps final atteint : '
  1497. (rv . 'PASDETPS' . 'TPS') ;
  1498. 'QUITTER' bloc1 ;
  1499. 'FINSI' ;
  1500. 'FIN' bloc1 ;
  1501.  
  1502. 'SI' testpr ;
  1503. rvp . 'PRESSION' = 'KCHT' rvpd 'SCAL' 'CENTRE'
  1504. (rvp . 'PRESSION') ;
  1505. rvp . 'PN' = 'ELNO' rvpd
  1506. (rvp . 'PRESSION') ;
  1507. 'FINSI' ;
  1508. ************************ E X E C ************************************
  1509. 'FINPROC' ;
  1510. opti dime 2 elem QUA4 ;
  1511.  
  1512. GRAPH=VRAI ;
  1513. GRAPH=FAUX ;
  1514. DISCR = LINE ;
  1515. NBIT = 3 ;
  1516.  
  1517. hz= 2. ;
  1518. r0=1.16 ;
  1519. p0=0 0 ;P1= r0 0 ; p2= r0 0.1 ; P3= 0 0.1 ;
  1520. n1 = 16 ;
  1521. n2 = 2 ;
  1522. nz = 2 ;
  1523. ab = p0 d n1 p1 ;
  1524. axe = p1 d n2 p2 ;
  1525. bc = p2 d n1 p3 ;
  1526. wall = P3 d n2 p0 ;
  1527.  
  1528. mtor= daller ab wall bc axe ;
  1529. Mmt= chan QUAF mtor ;
  1530. Maxe = chan QUAF axe ;
  1531. Mwall = chan QUAF wall ;
  1532. elim (Mwall et Mmt et Maxe) 1.e-5 ;
  1533. $wall = mode Mwall 'NAVIER_STOKES' DISCR ;
  1534. wall = doma $wall maillage ;
  1535. $axe = mode Maxe 'NAVIER_STOKES' DISCR ;
  1536. axe = doma $axe maillage ;
  1537. $mt= mode Mmt 'NAVIER_STOKES' DISCR ;
  1538. mt= doma $mt maillage ;
  1539.  
  1540. Si(EGA GRAPH VRAI) ; trace mt ; Finsi ;
  1541.  
  1542. * Beton
  1543. Lambda = 0.92 ;
  1544. SRC= 75 ;
  1545. tw= 14.6 ;
  1546. h=50.;
  1547. T0=8.;
  1548. DT=20. ;
  1549.  
  1550. rt=eqex EQUA 'T' $mt 'TEMPERAT'
  1551. 'TRAN' DT 'EUL_IMPL' 'TN'
  1552. 'OMEGA' 1. 'NITER' 1 ITMA 1 'FIDT' 1
  1553. ZONE $mt 'OPER' 'DFDT' 1. 'INCO' 'T'
  1554. ZONE $mt 'OPER' 'LAPN' lambda 'INCO' 'T'
  1555. ZONE $mt 'OPER' 'FIMP' src 'INCO' 'T'
  1556. ZONE $axe 'OPER' 'ECHI' 'H' 'T0' 'INCO' 'T'
  1557. ;
  1558.  
  1559. RT= eqex RT CLIM
  1560. 'T' TIMP wall tw ;
  1561. ;
  1562.  
  1563. rt.'METHINV'.TYPINV=1 ;
  1564. rt.'METHINV'.IMPINV=0 ;
  1565. rt.'METHINV'.NITMAX=400;
  1566. rt.'METHINV'.PRECOND=3 ;
  1567. rt.'METHINV'.RESID =1.e-8 ;
  1568. rt. 'METHINV' . 'FCPRECT'=1 ;
  1569. rt. 'METHINV' . 'FCPRECI'=1 ;
  1570.  
  1571. rt.'INCO'=table 'INCO';
  1572. rt.inco.'TN'= kcht $mt scal sommet 1. ;
  1573. rt.inco.'TNM'= kcht $mt scal sommet 1. ;
  1574. rt.inco.'H'= kcht $axe scal sommet h ;
  1575. rt.inco.'T0'= kcht $axe scal sommet T0 ;
  1576.  
  1577.  
  1578. xr= coor 1 mt ;
  1579. ilb=1./lambda/2. ;
  1580. Tlt=((lambda*Tw)+(r0*h*T0))/(lambda + (r0*h));
  1581. Tlt=Tlt + ((0.5 * src *r0 *r0)/(lambda + (r0*h)));
  1582. a=(Tlt - Tw + (0.5*src*r0*r0/lambda)/r0 ;
  1583. t =tw - (src * ilb * xr * xr) + (a * xr);
  1584. dg=doma $mt 'XXDIAGSI' ;
  1585.  
  1586. repeter bloc1 nbit ;
  1587. exec rt ;
  1588. tn= rt.inco.'TN' ;
  1589. ec= t - tn ;
  1590. el2= somt (dg * (ec*ec)) ;
  1591. mess 'Erreur L2 ' el2 ;
  1592. tl=extr tn 'SCAL' p1 ;
  1593. mess 'Tl = ' Tl 'Tl theo = ' Tlt;
  1594. fin bloc1 ;
  1595.  
  1596. evb= evol 'CHPO' tn ab;
  1597. evh= evol 'CHPO' tn (inve bc);
  1598. eva= evol 'CHPO' t ab;
  1599.  
  1600. tab1=table ;
  1601. tab1.1 = 'MARQ CROI REGU' ;
  1602. tab1.2 = 'MARQ TRIA REGU' ;
  1603. tab1.3 = 'MARQ LOSA REGU' ;
  1604. tab1.'TITRE' = table ;
  1605. tab1.'TITRE'. 1 = 'MOT' 'bas ' ;
  1606. tab1.'TITRE'. 2 = 'MOT' 'haut' ;
  1607. tab1.'TITRE'. 3 = 'MOT' 'Sol analytique' ;
  1608.  
  1609. Si GRAPH ;
  1610. dess (evh et evb et eva)
  1611. 'TITX' 'X (m)' 'TITY' 'Temperature (oC)'
  1612. LEGE TAB1 DATE MIMA TITR
  1613. ' Temperature ' ;
  1614. Finsi ;
  1615.  
  1616. si ( El2 > 1.E-8 ) ; erreur 5 ; finsi ;
  1617.  
  1618. opti dime 3 elem cub8 ;
  1619.  
  1620. mt0 = mtor volu nz 'TRAN' (0 0 hz) ;
  1621. mt1= syme mt0 'PLAN' p0 p1 (1 0 hZ) ;
  1622. mt = mt1 et mt0 ;
  1623. ax = axe et (syme axe 'POINT' p1) ;
  1624. axe = ax trans nz (0 0 hz) ;
  1625. wal = wall et (syme wall 'POINT' p0) ;
  1626. wall = wal trans nz (0 0 hz) ;
  1627. elim (mt et wall et axe) 1.e-5 ;
  1628.  
  1629. *** Cas 3D
  1630.  
  1631. Mmt= chan QUAF mt ;
  1632. Maxe = chan QUAF axe ;
  1633. Mwall = chan QUAF wall ;
  1634. elim (Maxe et Mmt et Mwall) 1.e-5 ;
  1635. $axe = mode Maxe 'NAVIER_STOKES' DISCR ;
  1636. axe = doma $axe maillage ;
  1637. $wall = mode Mwall 'NAVIER_STOKES' DISCR ;
  1638. wall = doma $wall maillage ;
  1639. $mt= mode Mmt 'NAVIER_STOKES' DISCR ;
  1640. mt= doma $mt maillage ;
  1641.  
  1642. * Granit
  1643. Lambda = 2.5;
  1644. SRC= -250.;
  1645. tw= 20 ;
  1646. t0= 37;
  1647. h=50. ;
  1648. DT=10. ;
  1649.  
  1650. rt=eqex EQUA 'RIGIDITE' 'T' $mt 'TEMPERAT'
  1651. 'TRAN' DT 'EUL_IMPL' 'TN'
  1652. *'PERM' 'TN'
  1653. 'OMEGA' 1. 'NITER' 1 ITMA 1 'FIDT' 1
  1654. * 'OPTI' 'MMDIAGO'
  1655. ZONE $mt 'OPER' 'DFDT' 1. 'INCO' 'T'
  1656. ZONE $mt 'OPER' 'LAPN' lambda 'INCO' 'T'
  1657. ZONE $mt 'OPER' 'FIMP' src 'INCO' 'T'
  1658. 'OPTI' 'MMDIAGO'
  1659. ZONE $axe 'OPER' 'ECHI' 'H' 'T0' 'INCO' 'T'
  1660. ;
  1661.  
  1662. RT= eqex RT CLIM
  1663. 'T' TIMP wall tw ;
  1664.  
  1665.  
  1666. rt.'METHINV'.TYPINV=3 ;
  1667. rt.'METHINV'.IMPINV=0 ;
  1668. rt.'METHINV'.NITMAX=400;
  1669. rt.'METHINV'.PRECOND=3 ;
  1670. rt.'METHINV'.RESID =1.e-8 ;
  1671. rt. 'METHINV' . 'FCPRECT'=1 ;
  1672. rt. 'METHINV' . 'FCPRECI'=1 ;
  1673.  
  1674. rt.'INCO'=table 'INCO';
  1675. rt.inco.'TN'= kcht $mt scal sommet 1. ;
  1676. rt.inco.'TNM'= kcht $mt scal sommet 1. ;
  1677. rt.inco.'H'= kcht $axe scal sommet h ;
  1678. rt.inco.'T0'= kcht $axe scal sommet T0 ;
  1679.  
  1680.  
  1681. xr= coor 1 mt ;
  1682. ilb=1./lambda/2. ;
  1683. Tlt=((lambda*Tw)+(r0*h*T0))/(lambda + (r0*h));
  1684. Tlt=Tlt + ((0.5 * src *r0 *r0)/(lambda + (r0*h)));
  1685. a=(Tlt - Tw + (0.5*src*r0*r0/lambda))/r0 ;
  1686. t =tw - (src * ilb * xr * xr) + (a * xr);
  1687. dg=doma $mt 'XXDIAGSI' ;
  1688.  
  1689. repeter bloc1 nbit ;
  1690. exec rt ;
  1691. tn=rt.inco.'TN' ;
  1692. ec= t - tn ;
  1693. el2= somt (dg * (ec*ec)) ;
  1694. mess 'Erreur L2 ' el2 ;
  1695. fin bloc1 ;
  1696.  
  1697. evb= evol 'CHPO' tn ab;
  1698. evh= evol 'CHPO' tn (inve bc);
  1699. eva= evol 'CHPO' t ab;
  1700.  
  1701. tab1=table ;
  1702. tab1.1 = 'MARQ CROI REGU' ;
  1703. tab1.2 = 'MARQ TRIA REGU' ;
  1704. tab1.3 = 'MARQ LOSA REGU' ;
  1705. tab1.'TITRE' = table ;
  1706. tab1.'TITRE'. 1 = 'MOT' 'bas ' ;
  1707. tab1.'TITRE'. 2 = 'MOT' 'haut' ;
  1708. tab1.'TITRE'. 3 = 'MOT' 'Sol analytique' ;
  1709.  
  1710. Si GRAPH ;
  1711. dess (evh et evb et eva)
  1712. 'TITX' 'X (m)' 'TITY' 'Temperature (oC)'
  1713. LEGE TAB1 DATE MIMA TITR
  1714. ' Temperature ' ;
  1715. Finsi ;
  1716.  
  1717. si ( El2 > 1.E-9 ) ; erreur 5 ; finsi ;
  1718.  
  1719. FIN ;
  1720.  
  1721.  
  1722.  
  1723.  
  1724.  
  1725.  
  1726.  
  1727.  
  1728.  
  1729.  
  1730.  
  1731.  

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