Télécharger exec.procedur

Retour à la liste

Numérotation des lignes :

  1. * EXEC PROCEDUR GOUNAND 22/11/22 21:15:01 11504
  2. ************************************************************************
  3. * EXECUTE UN ALGORITHME DÉCRIT DANS UNE TABL RV DE TYPE EQEX
  4. * CETTE TABL EST CRÉÉE PAR L'OPÉRATEUR EQEX
  5. ************************************************************************
  6. * HISTORIQUE
  7. * 20/12/99 = Rajout de la gestion de la matrice servant à l'assemblage
  8. * (rv.'METHINV'.'MATASS')
  9. * 12/05/06 = ajout STOPITER, NUITER et STOPPDT
  10. * 21/12/07 = ajout projection algébrique incrémentale (cf. GRESP)
  11. * 06/05/14 = ajout test DFDT : qd niter>1, on vérifie que le 2eme
  12. * argument est distinct de linconnue
  13. * 05/06/14 = on teste qd niter>1 et itma>1, car certains cas utilisent
  14. * DFDT pour pénaliser (dvisi et ccar3)
  15. * 05/02/16 = ajout de l'historique RV.'INCO'.'HIST'
  16. * 16/05/18 = corrections mineures + lifting
  17. ************************************************************************
  18. DEBP EXEC RV*'TABLE' NPARTI/'ENTIER' LOPP/'LISTMOTS' ;
  19. *
  20. *
  21. *
  22. *
  23. * +-------------------------------------------------------------------+
  24. * | |
  25. * | I N I T I A L I S A T I O N S E T V É R I F I C A T I O N S |
  26. * | |
  27. * +-------------------------------------------------------------------+
  28. *
  29. *
  30. ************************************************************************
  31. *
  32. * Ajout test DFDT : qd niter>1, on vérifie que le 2eme argument est
  33. * distinct de linconnue (anomalie #8050)
  34. SI (((RV.'NITER') > 1) ET ((RV.'ITMA') > 1)) ;
  35. NOP = DIME (RV.'LISTOPER') ;
  36. REPE IOP NOP ;
  37. NOMPER = EXTR (RV.'LISTOPER') &IOP ;
  38. SI (EGA NOMPER 'DFDT') ;
  39. NOMTAB = CHAI &IOP NOMPER ;
  40. RVN = RV.NOMTAB ;
  41. TYPARG2 = TYPE RVN.'ARG2' ;
  42. SI (EGA TYPARG2 MOT) ;
  43. NOMUN = EXTR (RVN.'LISTINCO') 1 ;
  44. SI (EGA NOMUN RVN.'ARG2') ;
  45. ERRE 1038 ;
  46. FINS ;
  47. FINS ;
  48. FINS ;
  49. FIN IOP ;
  50. FINS ;
  51. *
  52. ************************************************************************
  53. *
  54. * Utilisation de RESOU plutot que KRES ?
  55. TRESOU = FAUX ;
  56. SI (EXIS RV 'TRESOU') ;
  57. TRESOU = RV.'TRESOU' ;
  58. FINS ;
  59. *
  60. ************************************************************************
  61. *
  62. * Parallélisation des résolutions lorsqu'elles sont découplées
  63. TKPR = FAUX ;
  64. SI (EXIS RV 'TKPR') ;
  65. TKPR = RV.'TKPR' ;
  66. FINS ;
  67. *
  68. ************************************************************************
  69. *
  70. * Impressions liées au parallélisme
  71. IMPARA = FAUX ;
  72. SI (EXIS RV 'IMPARA') ;
  73. IMPARA = RV.'IMPARA' ;
  74. FINS ;
  75. *
  76. ************************************************************************
  77. *
  78. * Niveau d'impression global
  79. IIMPR = RV.'IMPR' ;
  80. *
  81. ************************************************************************
  82. *
  83. * Liste des opérateurs potentiellement parallélisables
  84. * LOPERPAR = MOTS 'NS' 'TSCA' 'DFDT' 'KONV' 'LAPN' 'FROT' 'FIMP' 'TOIMP' ;
  85. * LOPERPAR = MOTS 'NS' 'TSCA' 'DFDT' 'LAPN' 'KONV' ;
  86. LOPERPAR = MOTS 'NS' 'TSCA' 'DFDT' ;
  87. * LOPERPAR = MOTS 'NSQT' ;
  88. *
  89. * Parallélisation du calcul des matrices élémentaires si :
  90. * 1) RV.'TCPT' = VRAI
  91. * 2) NPART > 1
  92. * 3) On a bien des objets MMODEL et non pas des tables 'DOMAINE'
  93. *
  94. TCPT = FAUX ;
  95. SI (EXIS RV 'TCPT') ;
  96. TCPT = RV.'TCPT' ;
  97. FINS ;
  98. *
  99. * Préparation pour paralléliser le calcul des matrices élémentaires
  100. SI TCPT ;
  101. *
  102. * On partitionne en fonction du nombre d'assstants (= processeurs)
  103. * par défaut ou bien N (si précisé dans l'instruction EXEC RV N)
  104. SI (EXIS NPARTI) ;
  105. NPART = NPARTI ;
  106. SINON ;
  107. NPART = VALE 'ASSI' ;
  108. FINS ;
  109. *
  110. SI (NON (EXIS RV 'TABPARTM')) ;
  111. TABPARTM = TABL ;
  112. RV.'TABPARTM' = TABPARTM ;
  113. *
  114. REPE IOP (DIME (RV.'LISTOPER')) ;
  115. NOMPER = EXTR &IOP (RV.'LISTOPER') ;
  116. NOMTAB = CHAI &IOP NOMPER ;
  117. RVN = RV.NOMTAB ;
  118. $MD = RVN.'DOMZ' ;
  119. TYM = TYPE $MD ;
  120. *
  121. SI (NON (EGA TYM 'MMODEL')) ;
  122. TCPT = FAUX ;
  123. QUIT IOP ;
  124. FINS ;
  125. *
  126. SI (EXIS LOPERPAR NOMPER) ;
  127. SI IMPARA ;
  128. MESS 'On parallélise ' NOMPER ' NPART =' NPART ;
  129. FINS ;
  130. SI (NON (EXIS TABPARTM $MD)) ;
  131. $MDP = PART 'ARLE' NPART $MD ;
  132. *
  133. $BID = TABL 'ESCLAVE' ;
  134. REPE BIBI NPART ;
  135. * $BID.&BIBI = TABL 'KIZX' ;
  136. $BID.&BIBI = COPI RVN ;
  137. $BID.&BIBI.'DOMZ'= $MDP.&BIBI ;
  138. $TDZ = DOMA ($BID.&BIBI.'DOMZ') 'TABL' ;
  139. $BID.&BIBI.'TDOMZ'= $TDZ ;
  140. FIN BIBI ;
  141. *
  142. TABPARTM.'$bid' = $BID ;
  143. * TABPARTM.'$mdp' = $MDP ;
  144. RVN.'$mdp' = $MDP ;
  145. RVN.'$bid' = $BID ;
  146. FINS ;
  147. SINON ;
  148. SI IMPARA ;
  149. MESS 'On ne parallélise pas ' NOMPER ;
  150. FINS ;
  151. FINS ;
  152. *
  153. FIN IOP ;
  154. *
  155. SINON ;
  156. TABPARTM = RV.'TABPARTM' ;
  157. FINS ;
  158. *
  159. SI (EGA (DIME TABPARTM) 0) ;
  160. TCPT=FAUX ;
  161. SI IMPARA ;
  162. MESS 'Il n y a pas d opérateurs à paralléliser !' ;
  163. FINS ;
  164. FINS ;
  165. *
  166. SI (NPART <EG 1) ;
  167. SI IMPARA ;
  168. MESS 'NPART <EG 1 : On ne partitionne pas ' ;
  169. FINS ;
  170. TCPT = FAUX ;
  171. FINS ;
  172. *
  173. FINS ;
  174. *
  175. ************************************************************************
  176. *
  177. * Initialisation des critères d'arrêt (type LOGIQUE) positionnables
  178. * par l'utilisateur pendant la simulation
  179. *
  180. * (boucle interne)
  181. SI (NON (EXIS RV 'STOPITER')) ;
  182. RV.'STOPITER' = FAUX ;
  183. FINS ;
  184. *
  185. * (boucle externe)
  186. SI (NON (EXIS RV 'STOPPDT')) ;
  187. RV.'STOPPDT' = FAUX ;
  188. FINS ;
  189. *
  190. ************************************************************************
  191. *
  192. * LOGIQUE forçant le recalcul de la matrice de pression a chaque pas
  193. * de temps (boucle externe)
  194. CALPRE = FAUX ;
  195. SI (EXIS RV 'CALPRE') ;
  196. VERTYTAB RV 'CALPRE' 'LOGIQUE' ;
  197. CALPRE = RV.'CALPRE' ;
  198. FINS ;
  199. *
  200. ************************************************************************
  201. *
  202. * Option 'XEQUA' (non documentée...)
  203. SI (NON (EXIS RV 'XEQUA')) ;
  204. RV.'XEQUA' = FAUX ;
  205. FINS ;
  206. *
  207. ************************************************************************
  208. ************************************************************************
  209. ************************************************************************
  210. ************************************************************************
  211. *
  212. * Bifurcation gérée en interne (l'indice 'NAVISTOK' n'a pas vocation
  213. * a être ajouté par l'utilisateur)
  214. SI (EGA (RV.'NAVISTOK') 0) ;
  215. EXAC RV ;
  216. QUIT EXEC ;
  217. FINS ;
  218. *
  219. ************************************************************************
  220. ************************************************************************
  221. ************************************************************************
  222. ************************************************************************
  223. *
  224. * Schéma temporel pour Navier-Stokes
  225. TESTPR = EXIS RV 'PRESSION' ;
  226. TESTPRJ = EXIS RV 'PROJ' ;
  227. *
  228. SI TESTPR ;
  229. RVP = RV.'PRESSION' ;
  230. ACHP MATPR= KOPS 'MATRIK' ;
  231. RVP.'MATP' = MATPR ;
  232. FINS ;
  233. *
  234. SI TESTPRJ ;
  235. RVP = RV.'PROJ' ;
  236. FINS ;
  237. *
  238. SI (IIMPR > 0) ;
  239. SI TESTPR ;
  240. MESS 'Algorithme semi explicite (ANCIEN)' ;
  241. MESS '==================================' ;
  242. FINS ;
  243. SI TESTPRJ ;
  244. MESS 'Algorithme de Projection' ;
  245. MESS '========================' ;
  246. FINS ;
  247. SI (NON (TESTPR OU TESTPRJ)) ;
  248. MESS 'Algorithme standard implicite ou explicite ' ;
  249. MESS '===========================================' ;
  250. FINS ;
  251. FINS ;
  252. *
  253. ************************************************************************
  254. *
  255. * (ancien algo)
  256. * Choix du type de schéma de projection : |
  257. * PSCT VPI1 VPI2 PENA
  258. * -----------------------------------------+-----+-----+-----++------
  259. * TGRAD = Formulation en gradient | NON | NON | NON || NON
  260. * TMDM1 = Correction Gresho | NON | OUI | NON || NON
  261. * TPNM2 = Elimination end of step velocity | NON | NON | OUI || NON
  262. * |
  263. * (nouvel algo)
  264. TYPROJ = 'VPI1' ;
  265. ROW = 1. ;
  266. SI (EXIS RV 'TYPROJ') ;
  267. TYPROJ=RV.'TYPROJ' ;
  268. SI (NON (EXIS (MOTS 'PSCT' 'VPI1' 'VPI2' 'PENA') TYPROJ)) ;
  269. MESS '*****************************************************' ;
  270. MESS ' ERREUR ERREUR ERREUR ERREUR ERREUR ERREUR ' ;
  271. MESS ' ' ;
  272. MESS ' Le mot ' TYPROJ ' n existe pas dans la liste. ' ;
  273. MESS ' Les méthodes de projection autorisées sont : ' ;
  274. MESS ' VPI1, VPI2, PSCT et PENA ' ;
  275. MESS '*****************************************************' ;
  276. ERRE 21 ;
  277. QUIT EXEC ;
  278. FINS ;
  279. FINS ;
  280. *
  281. SI (EGA TYPROJ 'PSCT') ;
  282. TGRAD = FAUX ;
  283. TMDM1 = FAUX ;
  284. TPNM2 = FAUX ;
  285. FINS ;
  286. *
  287. SI (EGA TYPROJ 'VPI1') ;
  288. TGRAD = FAUX ;
  289. TMDM1 = VRAI ;
  290. TPNM2 = FAUX ;
  291. FINS ;
  292. *
  293. SI (EGA TYPROJ 'VPI2') ;
  294. TGRAD = FAUX ;
  295. TMDM1 = FAUX ;
  296. TPNM2 = VRAI ;
  297. FINS ;
  298. *
  299. SI (EGA TYPROJ 'PENA') ;
  300. TGRAD = FAUX ;
  301. TMDM1 = FAUX ;
  302. TPNM2 = FAUX ;
  303. FINS ;
  304. *
  305. ************************************************************************
  306. *
  307. * Affectations différentes selon que l'on est en 2D ou en 3D
  308. IDIM1 = VALE 'DIME' ;
  309. NOMVI = RV.'NOMVI' ;
  310. SI (EGA IDIM1 2) ;
  311. VNUL = 0.D0 0.D0 ;
  312. VUNI = 1.D0 1.D0 ;
  313. CNVI1 = CHAI '1' NOMVI ;
  314. CNVI2 = CHAI '2' NOMVI ;
  315. LC = MOTS CNVI1 CNVI2 ;
  316. LCU = MOTS 'UX' 'UY' ;
  317. SINON ;
  318. VNUL = 0.D0 0.D0 0.D0 ;
  319. VUNI = 1.D0 1.D0 1.D0 ;
  320. CNVI1 = CHAI '1' NOMVI ;
  321. CNVI2 = CHAI '2' NOMVI ;
  322. CNVI3 = CHAI '3' NOMVI ;
  323. LC = MOTS CNVI1 CNVI2 CNVI3 ;
  324. LCU = MOTS 'UX' 'UY' 'UZ' ;
  325. FINS ;
  326. *
  327. ************************************************************************
  328. *
  329. * Facteur de relaxation
  330. SI (NON (EXIS RV 'OMEGA')) ;
  331. OMEG = 1.D0 ;
  332. SINON ;
  333. OMEG = RV.'OMEGA' ;
  334. FINS ;
  335. *
  336. ************************************************************************
  337. *
  338. * Traitement d'un cas particulier = verrue (non documenté...)
  339. TESTRAN = TESTPR ET (EXIS RV 'CO') ;
  340. *
  341. ************************************************************************
  342. *
  343. * Initialisation de la table d'historique
  344. SI (NON (EXIS RV 'HIST')) ;
  345. RV.'HIST' = TABL ;
  346. FINS ;
  347. *
  348. ************************************************************************
  349. *
  350. * Nombre d'itérations externes (pas de temps)
  351. ITMAX = RV.'ITMA' ;
  352. SI (ITMAX < 1) ;
  353. ITMAX = 1 ;
  354. FINS ;
  355. *
  356. ************************************************************************
  357. *
  358. * Fréquence d'impression du solveur KRES
  359. IMPKRES = 0 ;
  360. *
  361. ************************************************************************
  362. *
  363. * Fréquence d'impression des pas de temps (opérateur TCRR)
  364. IMPTCRR = 0 ;
  365. SI (IIMPR >EG 0) ;
  366. SI (ITMAX EGA 1) ;
  367. IMPTCRR = 1 ;
  368. FINS ;
  369. SI TESTPRJ ;
  370. IMPTCRR = 1 ;
  371. FINS ;
  372. SI (NON (TESTPR OU TESTPRJ)) ;
  373. IMPTCRR = 1 ;
  374. FINS ;
  375. FINS ;
  376. *
  377. ************************************************************************
  378. *
  379. * FCPRECT : fréquence de recalcul des préc. dans la boucle externe
  380. * FCPRECI : fréquence de recalcul des préc. dans la boucle interne
  381. * FCPRECTP : idem pour la matrice de pression
  382. * FCPRECIP : idem pour la matrice de pression
  383. FCPRECT = RV. 'METHINV'.'FCPRECT' ;
  384. FCPRECI = RV. 'METHINV'.'FCPRECI' ;
  385. SI TESTPRJ ;
  386. FCPRECTP = RVP.'METHINV'.'FCPRECT' ;
  387. FCPRECIP = RVP.'METHINV'.'FCPRECI' ;
  388. FINS ;
  389. *
  390. * Force le calcul des préconditionneurs pour la première itération
  391. SI (NON (EXIS (RV.'METHINV') 'CALPREC')) ;
  392. RV.'METHINV'.'CALPREC' = VRAI ;
  393. FINS ;
  394. SI TESTPRJ ;
  395. SI (NON (EXIS (RVP.'METHINV') 'CALPREC')) ;
  396. RVP.'METHINV'.'CALPREC' = VRAI ;
  397. FINS ;
  398. FINS ;
  399. *
  400. ************************************************************************
  401. *
  402. * RESMN : résidu au pas de temps ou à l'itération interne précédente
  403. * RESMNP : idem pour la pression
  404. SI (NON (EXIS RV 'resmn')) ;
  405. RESMN MAPREC = KOPS 'MATRIK' ;
  406. RV.'resmn' = RESMN ;
  407. FINS ;
  408. SI (NON (EXIS RV 'resmnp')) ;
  409. RESMNP MAPRECP = KOPS 'MATRIK' ;
  410. RV.'resmnp' = RESMNP ;
  411. FINS ;
  412. *
  413. *
  414. *
  415. *
  416. * +-------------------------------------------------------------------+
  417. * | |
  418. * | A L G O R I T H M E D E R É S O L U T I O N |
  419. * | |
  420. * +-------------------------------------------------------------------+
  421. *
  422. * *******************************************
  423. * *******************************************
  424. * DEBUT DE LA BOUCLE EXTERNE (= PAS DE TEMPS)
  425. * *******************************************
  426. * *******************************************
  427. REPE BLOC1 ITMAX ;
  428. *
  429. * Doit-on recalculer les préconditionneurs ?
  430. SI (MULT &BLOC1 FCPRECT) ;
  431. RV.'METHINV'.'CALPREC' = VRAI ;
  432. FINS ;
  433. SI TESTPRJ ;
  434. SI (MULT &BLOC1 FCPRECTP) ;
  435. RVP.'METHINV'.'CALPREC' = VRAI ;
  436. FINS ;
  437. FINS ;
  438. *
  439. * Doit-on recalculer la matrice de pression ?
  440. TESTP1 = FAUX ;
  441. SI (TESTPR OU TESTPRJ) ;
  442. SI CALPRE ;
  443. TESTP1 = VRAI ;
  444. FINS ;
  445. SI (NON (EXIS RVP 'MATC')) ;
  446. TESTP1 = VRAI ;
  447. FINS ;
  448. FINS ;
  449. *
  450. *
  451. * ********************************************
  452. * ********************************************
  453. * DEBUT DE LA BOUCLE INTERNE (= NON-LINEARITE)
  454. * ********************************************
  455. * ********************************************
  456. REPE BLOCI (RV.'NITER') ;
  457. *
  458. * Numéro de l'itération interne courante
  459. RV.'NUITER' = &BLOCI ;
  460. *
  461. * Doit-on recalculer les préconditionneurs ?
  462. SI (MULT &BLOCI FCPRECI) ;
  463. RV.'METHINV'.'CALPREC' = VRAI ;
  464. FINS ;
  465. SI TESTPRJ ;
  466. SI (MULT &BLOCI FCPRECIP) ;
  467. RVP.'METHINV'.'CALPREC' = VRAI ;
  468. FINS ;
  469. FINS ;
  470. *
  471. *
  472. * ==========================================================
  473. * EXÉCUTION DE TOUS LES OPÉRATEURS/PROCÉDURES DE LA TABLE RV
  474. * => CONSTRUCTION DU SYSTÈME D'ÉQUATIONS DE Q.D.M.
  475. * ==========================================================
  476. *
  477. ST MAT = KOPS 'MATRIK' ;
  478. SF MAU = KOPS 'MATRIK' ;
  479. MDFDT = 0 ;
  480. REPE BLOC2 (DIME (RV.'LISTOPER')) ;
  481. NOMPER = EXTR &BLOC2 (RV.'LISTOPER') ;
  482. NOMTAB = CHAI &BLOC2 NOMPER ;
  483. RVN = RV.NOMTAB ;
  484. *
  485. * KFORM vaut 0 (EFM1) ou 1 (EF) ou 2 (VF) ou 3 (EFMC)
  486. MDFDT = MDFDT + RVN.'KOPT'.'KFORM' ;
  487. *
  488. * Exécution séquentielle
  489. SI (NON TCPT) ;
  490. MSI MAI = (TEXTE NOMPER) RVN ;
  491. *
  492. * Exécution parallèle (mémoire partagée)
  493. SINON ;
  494. SI (EXIS LOPERPAR NOMPER) ;
  495. SI IMPARA ;
  496. MESS ' On parallélise ' NOMPER ' NPART =' NPART ;
  497. FINS ;
  498. $MDP = RVN.'$mdp' ;
  499. $BID = RVN.'$bid' ;
  500. MSI MAI = ASSI 'TOUS' (TEXTE NOMPER) RVN $BID ;
  501. MSI = ETG MSI ;
  502. MAI = ETG MAI ;
  503. SINON ;
  504. SI IMPARA ;
  505. MESS NOMPER ' n est pas parallélisable' ;
  506. FINS ;
  507. MSI MAI = (TEXTE NOMPER) RVN ;
  508. FINS ;
  509. FINS ;
  510. *
  511. SI (EGA NOMPER 'DFDT') ;
  512. * ISCHT vaut 1 (BDF2) ou 2 (BDF4) ou 0 sinon (ordre 1)
  513. ISCHT = RVN.KOPT.'ISCHT' ;
  514. MAT = MAT ET MAI ;
  515. ST = ST ET MSI ;
  516. SINON ;
  517. MAU = MAU ET MAI ;
  518. SF = SF ET MSI ;
  519. FINS ;
  520. FIN BLOC2 ;
  521. *
  522. * ST/MAT contient le second membre et la matrice pour DFDT
  523. * SF/MAU contient le second membre et la matrice pour les autres opérateurs
  524. * S2/MA1 contient le second membre et la matrice pour l'ensemble
  525. S2 = SF ET ST ;
  526. MA1 = MAU ET MAT ;
  527. *
  528. *
  529. * ====================================================
  530. * TRAITEMENTS SPÉCIFIQUES POUR LA PRESSION (SI REQUIS)
  531. * ====================================================
  532. *
  533. SI (EXIS RV 'rvpd') ;
  534. RVPD = RV.'rvpd' ;
  535. IDIGV = RV.'IDigv' ;
  536. DIGV = RV.'Digv' ;
  537. MATPR = RVP.'MATP' ;
  538. FINS ;
  539. *
  540. *
  541. * 1) Calcul de C*D^(-1)*Ct pour l'algorithme semi-explicite (RV.'PRESSION')
  542. * **********************************************************************
  543. *
  544. SI (TESTPR ET TESTP1) ;
  545. RVPD = RVP.'DOMAINE' ;
  546. DIAGO = DOMA RVPD 'XXDIAGSI' ;
  547. SI (IDIM1 EGA 2) ;
  548. DIGV = (NOMC DIAGO CNVI1) ET (NOMC DIAGO CNVI2) ;
  549. SINON ;
  550. DIGV = (NOMC DIAGO CNVI1) ET (NOMC DIAGO CNVI2) ET (NOMC DIAGO CNVI3) ;
  551. FINS ;
  552. DIGV = KCHT RVPD 'VECT' 'SOMMET' 'COMP' LC DIGV ;
  553. *
  554. * Les noeuds de DIGV sur lesquels existe une CLIM en vitesse sont mis a 1.E30
  555. SI (EXIS RV 'CLIM') ;
  556. RVP.'CLIM' = RV.'CLIM' ;
  557. DIGV = KOPS DIGV 'CLIM' (RV.'CLIM') -1 ;
  558. FINS ;
  559. *
  560. IDIGV = INVE DIGV ;
  561. *
  562. SI (NON (EXIS RVP 'DIAGV')) ;
  563. RVP.'DIAGV' = DIGV ;
  564. FINS ;
  565. *
  566. * KMAB => discrétisation de l'opérateur div(U)
  567. RVP.'MATC' = KMAB RVP ;
  568. RVP.'PRESSION' = KCHT RVPD 'SCAL' 'CENTRE' 0.D0 ;
  569. RVP.'GRADP' = KCHT RVPD 'VECT' 'SOMMET' VNUL ;
  570. ACHP MATPR = KOPS 'MATRIK' ;
  571. RVP.'MATP' = MATPR ;
  572. RV.'rvpd' = RVPD ;
  573. RV.'IDigv' = IDIGV ;
  574. RV.'Digv' = DIGV ;
  575. *
  576. FINS ;
  577. *
  578. *
  579. * 2) Calcul de C*D^(-1)*Ct pour l'algorithme de projection (RV.'PROJ')
  580. * *****************************************************************
  581. *
  582. SI (TESTPRJ ET TESTP1) ;
  583. *
  584. * -------------------------------------------------------------------
  585. * Construction de la matrice masse diagonale à partir des tables DFDT
  586. * -------------------------------------------------------------------
  587. *
  588. DIAGO MMA = KOPS 'MATRIK' ;
  589. IDFDT = 0 ;
  590. REPE BLOCJ (DIME (RV.'LISTOPER')) ;
  591. NOMPER = EXTR &BLOCJ (RV.'LISTOPER') ;
  592. NOMTAB = CHAI &BLOCJ NOMPER ;
  593. SI ((EGA NOMPER 'DFDT') ET (EXIS (RV.NOMTAB) 'LISTINCO')) ;
  594. SI (EXIS (RV.NOMTAB.'LISTINCO') NOMVI) ;
  595. IDFDT = IDFDT + 1 ;
  596. RVPD = RV.NOMTAB.'DOMZ' ;
  597. DIAGO = DIAGO ET (DOMA RVPD 'XXDIAGSI') ;
  598. FINS ;
  599. FINS ;
  600. FIN BLOCJ ;
  601. *
  602. SI (EGA IDFDT 0) ;
  603. MESS 'Pas de DFDT ????' ;
  604. ERRE 21 ;
  605. FINS ;
  606. *
  607. SI (IDIM1 EGA 2) ;
  608. DIGV = (NOMC DIAGO CNVI1) ET (NOMC DIAGO CNVI2) ;
  609. SINON ;
  610. DIGV = (NOMC DIAGO CNVI1) ET (NOMC DIAGO CNVI2) ET (NOMC DIAGO CNVI3) ;
  611. FINS ;
  612. *
  613. DIGV = KCHT RVPD 'VECT' 'SOMMET' 'COMP' LC DIGV ;
  614.  
  615. * DIG1 = matrice diagonale indépendante des conditions aux limites
  616. DIG1 = COPI DIGV ;
  617. IDIAG = INVE DIG1 ;
  618. *
  619. * Les noeuds de DIGV sur lesquels existe une CLIM en vitesse sont mis a 1.E30
  620. DIGV = KOPS DIGV 'CLIM' (RV.'CLIM') -1 ;
  621. *
  622. IDIGV = INVE DIGV ;
  623. *
  624. SI (NON (EXIS RVP 'DIAGV')) ;
  625. RVP.'DIAGV' = DIGV ;
  626. FINS ;
  627. RVP.'INCO' = RV.'INCO' ;
  628. *
  629. *
  630. * -----------------------------------------------------------
  631. * Exécution de tous les opérateurs/procédures de la table RVP
  632. * -----------------------------------------------------------
  633. *
  634. * TVNPC et TVNP indiquent que l'opérateur VNIMP est présent
  635. * (respectivement avec MSOMMET ou CENTREP0/CENTREP1)
  636. TVNP = FAUX ;
  637. TVNPC = FAUX ;
  638. *
  639. * SVNPC/MVNPC contient le second membre et la matrice pour VNIMP+MSOMMET
  640. * SR/MAR contient le second membre et la matrice pour les VNIMP+CENTREPx
  641. * SP/MAP contient le second membre et la matrice pour les autres opérateurs
  642. SVNPC MVNPC = KOPS 'MATRIK' ;
  643. SR MAR = KOPS 'MATRIK' ;
  644. SP MAP = KOPS 'MATRIK' ;
  645. *
  646. REPE BLOCPJ (DIME (RVP.'LISTOPER')) ;
  647. NOMPER = EXTR &BLOCPJ (RVP.'LISTOPER') ;
  648. NOMTAB = CHAI &BLOCPJ NOMPER ;
  649. *
  650. SI (EGA NOMPER 'KBBT') ;
  651. * IKOMP vaut 0 (NOCONS) ou 1 (CONS) ou 2 (CONS2)
  652. RVP.NOMTAB.'KOPT'.'IKOMP' = 1 ;
  653. FINS ;
  654. MSI MAI = (TEXTE NOMPER) (RVP.NOMTAB) ;
  655. *
  656. *######################## ATTENTION !!!!!!!!
  657. *######################## MSIG/MAIG SONT UTILISES HORS DE BLOCPJ
  658. *######################## LEUR VALEUR EST INCERTAINE......
  659. SI TGRAD ;
  660. SI (IIMPR > 0) ;
  661. MESS 'Procedure PROJCT (Gradient) Operateur : ' NOMPER ;
  662. FINS ;
  663. * IKOMP vaut 0 (NOCONS) ou 1 (CONS) ou 2 (CONS2)
  664. RVP.NOMTAB.'KOPT'.'IKOMP' = 0 ;
  665. MSIG MAIG = (TEXTE NOMPER) (RVP.NOMTAB) ;
  666. FINS ;
  667. *########################
  668. *########################
  669. *########################
  670. *
  671. SI (EGA NOMPER 'VNIMP') ;
  672. * DISCPRES vaut 3 (CENTREP0) ou 4 (CENTREP1) ou 5 (MSOMMET)
  673. SI (EGA (RVP.'DISCPRES') 5) ;
  674. TVNPC = VRAI ;
  675. MVNPC = MVNPC ET MAI ;
  676. SVNPC = SVNPC ET MSI ;
  677. SINON ;
  678. TVNP = VRAI ;
  679. MAR = MAR ET MAI ;
  680. SR = SR ET MSI ;
  681. FINS ;
  682. SINON ;
  683. MAP = MAP ET MAI ;
  684. SP = SP ET MSI ;
  685. FINS ;
  686. FIN BLOCPJ ;
  687. *
  688. * Séparation des sous-matrices type C*Ct (MAC) des autres sous-matrices (MATPC)
  689. MATPC MAC = KOPS 'CMCTSPLT' MAP ;
  690. *
  691. *######################## ATTENTION !!!!!!!!
  692. *######################## LA VALEUR DE MAIG EST INCERTAINE......
  693. SI TGRAD ;
  694. MATPCG MACG = KOPS 'CMCTSPLT' MAIG ;
  695. RVP.'MATG' = MACG ;
  696. FINS ;
  697. *########################
  698. *########################
  699. *
  700. SI TVNPC ;
  701. MAC = MAC ET MVNPC ;
  702. FINS ;
  703. RVP.'MATC' = MAC ;
  704. *
  705. SCR MCR = KOPS 'MATRIK' ;
  706. SI TVNP ;
  707. RVP.'MBTR' = MAR ;
  708. DUNIT = IDIGV ;
  709. * KOPS 'CMCT' A B C calcule la matrice A*C^(-1)*Bt
  710. CRT = KOPS 'CMCT' MAC MAR DUNIT ;
  711. CTR = KOPS 'CMCT' MAR MAC DUNIT ;
  712. RRT = KOPS 'CMCT' MAR MAR DUNIT ;
  713. * MCR = MCR ET CRT ET CTR ET RRT ;
  714. MCR = CRT ET CTR ET RRT ;
  715. FINS ;
  716. RVP.'TVNP' = TVNP ;
  717. *
  718. SI (EGA (RVP.'DISCPRES') 5) ;
  719. SI (IIMPR > 0) ;
  720. MESS 'Cas des pressions CONTINUES' ;
  721. FINS ;
  722. * MATPR = MATPC ET MCR ;
  723. * MCR est vide !
  724. MATPR = MATPC ;
  725. SINON ;
  726. SI (IIMPR > 0) ;
  727. MESS 'Cas des pressions DISCONTINUES' ;
  728. FINS ;
  729. MATPR = KOPS 'CMCT' MAC MAC IDIGV ;
  730. MATPR = MATPC ET MATPR ET MCR ;
  731. *
  732. SI TRESOU ;
  733. MATPR = KOPS 'RIMA' MATPR 'NSYM' ;
  734.  
  735. INCU = EXTR MAC 'COMP' ;
  736. INCP = EXTR MAC 'COMP' 'DUAL' ;
  737. INCUP = INCU ET INCP ;
  738. OPTI 'INCO' INCUP INCUP ;
  739. *
  740. * RIGIDITE/CHPOINT associés aux C.L. en pression
  741. Si('EXIS' rvp 'rigp') ;
  742. rigp=rvp.'rigp' ;
  743. fpi =rvp.'climreso' ;
  744. Sinon ;
  745. incpres=extr incp 1 ;
  746. mip=extr rvp.clim 'MAILLAGE' ;
  747. rigp=BLOQ mip incpres ;
  748. fpi =DEPI RIGP rvp.'CLIM' ;
  749. rvp.'rigp'=rigp ;
  750. rvp.'climreso'=fpi ;
  751. Finsi ;
  752. matpr = matpr et rigp ;
  753.  
  754.  
  755. SI (IIMPR > 0) ;
  756. MESS '***** ON UTILISE RESOU POUR LA PRESSION *****' ;
  757. FINS ;
  758. SINON ;
  759. SI (IIMPR > 0) ;
  760. MESS '***** ON UTILISE KRES POUR LA PRESSION *****' ;
  761. FINS ;
  762. FINS ;
  763. FINS ;
  764. *
  765. RV.'rvpd' = RVPD ;
  766. RV.'IDigv' = IDIGV ;
  767. RV.'Digv' = DIGV ;
  768. RVP.'MATP' = MATPR ;
  769. *
  770. FINS ;
  771. *
  772. *
  773. * 3) Calcul de Ct*P pour l'algorithme semi-explicite (RV.'PRESSION')
  774. * ***************************************************************
  775. *
  776. SI TESTPR ;
  777. DT = (RV.'PASDETPS'.'DELTAT') * (RV.'ALFA') ;
  778. RVP.'DELTAT' = DT ;
  779. *
  780. * F = COPI S2 ;
  781. * LC = EXTR DIGV 'COMP' ;
  782. F = S2 ;
  783. FU = KCHT RVPD 'VECT' 'SOMMET' 'COMP' LC (EXCO F LC) ;
  784. *
  785. * DM1F vaut (-1 * RV.'CLIM' / DT) aux noeuds supportant une C.L.
  786. * DM1F vaut (-1 * FU * IDIGV) partout ailleurs
  787. SI (EXIS RV 'CLIM') ;
  788. DM1F = KOPS (DT '*' (KOPS FU '*' IDIGV)) 'CLIM' (RV.'CLIM') 3 ;
  789. DTI = -1.D0 / DT ;
  790. DM1F = DTI * DM1F ;
  791. SINON ;
  792. DM1F = (-1.D0) * (KOPS FU '*' IDIGV) ;
  793. FINS ;
  794. *
  795. * Calcul de -C*D^(-1)*F
  796. RVP.'PRESSION' = KMF (RVP.'MATC') DM1F ;
  797. *
  798. * /\/\/\/\/\/\/\/\/\/\/
  799. * Calcul de la pression (vieille syntaxe KRES)
  800. * /\/\/\/\/\/\/\/\/\/\/
  801. KRES RVP (RVP.'PRESSION') 'BETA' (RVP.'KBETA') (RVP.'BETA')
  802. 'PIMP' (RVP.'KPIMP') (RVP.'PIMP') ;
  803. RV.'INCO'.'PRESSION' = RVP.'PRESSION' ;
  804.  
  805. * Calcul de CtP, c'est-à-dire grad(p)
  806. RVP.'GRADP' = KMTP 1 (RVP.'MATC') (RVP.'PRESSION') LC ;
  807. S2 = S2 + (RVP.'GRADP') ;
  808. FINS ;
  809. *
  810. *
  811. * 4) Calcul de M*D^(-1) puis Ct*P pour l'algorithme de projection (RV.'PROJ')
  812. * ************************************************************************
  813. *
  814. SI TESTPRJ ;
  815. DT = (RV.'PASDETPS'.'DELTAT') * (RV.'ALFA') ;
  816. RVP.'DELTAT' = DT ;
  817. *
  818. *
  819. * ------------------------------------
  820. * Calcul du produit M*D^(-1) si requis
  821. * ------------------------------------
  822. *
  823. SI TMDM1 ;
  824. SI ((EXIS RV 'MDM1') ET (NON CALPRE)) ;
  825. MDM1 = RV.'MDM1' ;
  826. SINON ;
  827. SI (IIMPR > 0) ;
  828. MESS 'On calcule M*D^(-1)' ;
  829. FINS ;
  830. *
  831. STN MATN = KOPS 'MATRIK' ;
  832. ROW = 0. ;
  833. REPE BLOCPJ1 (DIME (RVP.'LISTOPER')) ;
  834. NOMPER1 = EXTR &BLOCPJ1 (RVP.'LISTOPER') ;
  835. NOMTAB1 = CHAI &BLOCPJ1 NOMPER1 ;
  836. SI (EGA NOMPER1 'KBBT') ;
  837. NOMIV1 = EXTR (RVP.NOMTAB1.'LISTINCO') 1 ;
  838. *
  839. REPE BLOCPJ2 (DIME (RV.'LISTOPER')) ;
  840. NOMPER2 = EXTR &BLOCPJ2 (RV.'LISTOPER') ;
  841. NOMTAB2 = CHAI &BLOCPJ2 NOMPER2 ;
  842. SI (EGA NOMPER2 'DFDT') ;
  843. NOMIV2 = EXTR (RV.NOMTAB2.'LISTINCO') 1 ;
  844. *
  845. SI (EGA NOMIV2 NOMIV1) ;
  846. SI (EGA TYPROJ 'PENA') ;
  847. RWI = (RV.NOMTAB2).ARG1 ;
  848. TRWI = TYPE RWI ;
  849. SI (EGA TRWI 'MOT') ;
  850. RWI = RV.'INCO'.RWI ;
  851. TRWI = TYPE RWI ;
  852. FINS ;
  853. SI (IIMPR > 0) ;
  854. MESS 'Type du coefficient RWI = ' TRWI ;
  855. FINS ;
  856. *
  857. SI (EGA TRWI 'FLOTTANT') ;
  858. ROW = MAXI (PROG RWI ROW) ;
  859. FINS ;
  860. SI ((EGA TRWI 'CHPOINT') OU (EGA TRWI 'MCHAML')) ;
  861. ROW = MAXI (PROG (MAXI RWI) ROW) ;
  862. FINS ;
  863. SINON ;
  864. ROW = 1. ;
  865. FINS ;
  866. * MESS 'ROW = ' ROW ;
  867. *
  868. MSI MAI = DFDT (RV.NOMTAB2) ;
  869. MATN = MATN ET MAI ;
  870. * DOMZP = RV.NOMTAB2.'DOMZ' ;
  871. * STI = KCHT DOMZP 'VECT' 'SOMMET' 'COMP' LC VUNI ;
  872. FINS ;
  873. FINS ;
  874. FIN BLOCPJ2 ;
  875. FINS ;
  876. FIN BLOCPJ1 ;
  877. *
  878. SI (EGA TYPROJ 'PENA') ;
  879. MDM1 = 1. ;
  880. SINON ;
  881. MDM1 = KMF MATN IDIAG ;
  882. FINS ;
  883. *
  884. SI (EGA ISCHT 1) ;
  885. MDM1 = MDM1 * ((DT*2.)/3.) ;
  886. SINON ;
  887. MDM1 = MDM1 * DT ;
  888. FINS ;
  889. *
  890. RV.'MDM1' = MDM1 ;
  891. FINS ;
  892. FINS ;
  893. *
  894. *
  895. * ---------------------------
  896. * Calcul du produit Ct*P[n-1]
  897. * ---------------------------
  898. *
  899. TVNP = RVP.'TVNP' ;
  900.  
  901. SI (EXIS (RV.'INCO') 'PRESSION') ;
  902. PPI = RV.'INCO'.'PRESSION' ;
  903. *
  904. * TPNM2 Elimination of end of step velocity (2*P[n] - P[n-1])
  905. SI (EXIS (RV.'INCO') 'PNM2') ;
  906. PPI = (2 * PPI) - (RV.'INCO'.'PNM2') ;
  907. FINS ;
  908. *
  909. * TGRAD=VRAI => Formulation en u*grad(p)
  910. * TGRAD=FAUX => Formulation en p*div(u)
  911. SI TGRAD ;
  912. CPRE = KMF (RVP.'MATG') PPI 'TRAN' ;
  913. SINON ;
  914. CPRE = KMF (RVP.'MATC') PPI 'TRAN' ;
  915. FINS ;
  916. *
  917. * Prise en compte des C.L. (opérateur VNIMP avec CENTREP0/CENTREP1)
  918. SI TVNP ;
  919. CXRE = KMF (RVP.'MBTR') PPI 'TRAN' ;
  920. CPRE = CPRE ET CXRE ;
  921. FINS ;
  922.  
  923. * Consistence selon Gresho
  924. SI TMDM1 ;
  925. GRADPRES = KOPS MDM1 '*' CPRE ;
  926. * Consistence selon Guermond
  927. SINON ;
  928. GRADPRES = CPRE ;
  929. FINS ;
  930. OUBL CPRE ;
  931. RV.'INCO'.'GRADPRES' = NOMC GRADPRES LCU LC ;
  932. *
  933. SI TGRAD ;
  934. S2 = SF - (RV.'INCO'.'GRADPRES') + ST ;
  935. SINON ;
  936. S2 = SF + (RV.'INCO'.'GRADPRES') + ST ;
  937. FINS ;
  938. *
  939. FINS ;
  940. FINS ;
  941. *
  942. SI (EXIS RV 'CLIM') ;
  943. S1 = RV.'CLIM' ;
  944. SINON ;
  945. S1 = ' ' ;
  946. FINS ;
  947. RV.'S2' = S2 ;
  948. *
  949. *
  950. * ===========================================
  951. * RÉSOLUTION DU SYSTÈME D'ÉQUATIONS DE Q.D.M.
  952. * ===========================================
  953. *
  954. SI ((NON TESTPRJ) OU (EGA MDFDT 0)) ;
  955. *
  956. RV.'METHINV'.'XINIT' = RV.'resmn' ;
  957. *
  958. * Méthode de projection incrémentale
  959. * **********************************
  960. SI (EXIS RV 'GPROJ') ;
  961. SI (IIMPR > 0) ;
  962. MESS 'Résolution QDM - méthode de projection incrémentale' ;
  963. FINS ;
  964. RES = GRESP MA1 S1 S2 RV ;
  965. *
  966. * Système monolithique d'équations (U-P)
  967. * **************************************
  968. SINON ;
  969. SI (IIMPR > 0) ;
  970. MESS 'Résolution QDM - système monolithique (U-P)' ;
  971. FINS ;
  972. *
  973. SI ((RV.'METHINV'.'TYPINV' EGA 1) ET TRESOU) ;
  974. SI (IIMPR > 0) ;
  975. MESS 'ILU-Crout => on utilise RESOU' ;
  976. FINS ;
  977. *
  978. MA1 = KOPS 'RIMA' MA1 'NSYM' ;
  979. INCPRES = EXTR RV.'LISTINCO' 1 ;
  980. OPTI 'INCO' INCPRES INCPRES ;
  981. *
  982. * RIGIDITE/CHPOINT associés aux C.L.
  983. SI (EXIS RV 'trico1') ;
  984. TRICO = RV.'trico1' ;
  985. SINON ;
  986. INC = EXTR (RV.'CLIM') 'COMP' ;
  987. TRICO = TABL 'ESCLAVE' ;
  988. RV.'trico1' = TRICO ;
  989. NCMP = DIME INC ;
  990. REPE BINC NCMP ;
  991. ICO = EXTR &BINC INC ;
  992. RCLIM = EXCO (RV.'CLIM') ICO ICO ;
  993. MIP = EXTR RCLIM 'MAILLAGE' ;
  994. TRICO.&BINC = BLOQ MIP ICO ;
  995. FIN BINC ;
  996. FINS ;
  997. S1 = VIDE 'CHPOINT/DISCRET' ;
  998. REPE BINC NCMP ;
  999. ICO = EXTR &BINC INC ;
  1000. RCLIM = EXCO (RV.'CLIM') ICO ICO ;
  1001. RIGP = TRICO.&BINC ;
  1002. S1 = S1 ET (DEPI RIGP RCLIM) ;
  1003. FIN BINC ;
  1004. *
  1005. RIGP = ET TRICO ;
  1006. MA1 = MA1 ET RIGP ;
  1007. FINS ;
  1008. *
  1009. RES = KRESP MA1 'TYPI' (RV.'METHINV')
  1010. 'CLIM' S1
  1011. 'SMBR' S2
  1012. 'IMPR' IMPKRES ;
  1013. FINS ;
  1014. FINS ;
  1015. *
  1016. * Méthode de projection standard
  1017. * ******************************
  1018. SI (TESTPRJ ET (NEG MDFDT 0)) ;
  1019. SI (IIMPR > 0) ;
  1020. MESS 'Résolution QDM - méthode de projection standard' ;
  1021. FINS ;
  1022. *
  1023. * LPART = EXTR MA1 'COMP' ;
  1024. LPART = KOPS 'EXTRCOUP' MA1 ;
  1025. NBPART = DIME LPART ;
  1026. * MESS 'Liste des composantes' ;
  1027. * LIST LPART ;
  1028. *
  1029. $MA1 = TABL 'ESCLAVE' ;
  1030. $S1 = TABL 'ESCLAVE' ;
  1031. $S2 = TABL 'ESCLAVE' ;
  1032. $RESMN = TABL 'ESCLAVE' ;
  1033. $TAB1 = TABL 'ESCLAVE' ;
  1034. *
  1035. SI (EXIS RV '$trico') ;
  1036. $TRICO = RV.'$trico' ;
  1037. SINON ;
  1038. $TRICO = TABL 'ESCLAVE' ;
  1039. RV.'$trico' = $TRICO ;
  1040. REPE BCLCOM NBPART ;
  1041. $TRICO.&BCLCOM = TABL 'ESCLAVE' ;
  1042. FIN BCLCOM ;
  1043. FINS ;
  1044. *
  1045. * On boucle sur les partitions du domaine (déjà déterminées plus haut)
  1046. REPE BCLCOM NBPART ;
  1047. * NMC = EXTR LPART &BCLCOM ;
  1048. NMC = LPART.&BCLCOM ;
  1049. * MESS 'NMC=' NMC &BCLCOM ;
  1050.  
  1051. * MA1I = EXTR MA1 (MOTS NMC) (MOTS NMC) ;
  1052. MA1I = KOPS 'EXTRINCO' NMC NMC MA1 ;
  1053. $MA1.&BCLCOM = MA1I ;
  1054. * $S1.&BCLCOM = EXCO S1 NMC NMC ;
  1055. * $S2.&BCLCOM = EXCO S2 NMC NMC ;
  1056. * RESMNI = EXCO (RV.'resmn') NMC 'NOID' NMC ;
  1057. $S1.&BCLCOM = EXCO S1 NMC 'NOID' ;
  1058. $S2.&BCLCOM = EXCO S2 NMC 'NOID' ;
  1059. RESMNI = EXCO (RV.'resmn') NMC 'NOID' ;
  1060. TAB1 = COPI RV.'METHINV' ;
  1061. * TAB1 = RV.'METHINV' ;
  1062. TAB1.'XINIT' = RESMNI ;
  1063. $TAB1.&BCLCOM = TAB1 ;
  1064. *
  1065. SI (IIMPR > 0) ;
  1066. MESS 'Mise à jour du préconditionnement' ;
  1067. FINS ;
  1068. SI (TAB1.'TYPINV' >EG 2) ;
  1069. LWORD1 LWORD2 = KOPS 'EXTRNINC' MA1I ;
  1070. WORD1 = EXTR LWORD1 1 ;
  1071. *
  1072. SI (NON (EXIS RV 'TABRES')) ;
  1073. RV.'TABRES' = TABL ;
  1074. FINS ;
  1075. TABRES = RV.'TABRES' ;
  1076. *
  1077. SI (TAB1.'CALPREC') ;
  1078. SI (NON (EXIS TABRES WORD1)) ;
  1079. TABRES.WORD1 = TABL ;
  1080. FINS ;
  1081. TABRES.WORD1.'MATASS' = MA1I ;
  1082. TABRES.WORD1.'MAPREC' = MA1I ;
  1083. FINS ;
  1084. TAB1.'MATASS' = TABRES.WORD1.'MATASS' ;
  1085. TAB1.'MAPREC' = TABRES.WORD1.'MAPREC' ;
  1086. FINS ;
  1087. FIN BCLCOM ;
  1088. *
  1089. * RÉSOLUTION PARALLÈLE
  1090. SI TKPR ;
  1091. SI IMPARA ;
  1092. MESS 'Les résolutions sont traitées en parallèle' ;
  1093. FINS ;
  1094. *
  1095. RES = ASSI 'TOUS' KRES $MA1 'TYPI' $TAB1
  1096. 'CLIM' $S1
  1097. 'SMBR' $S2
  1098. 'IMPR' IMPKRES ;
  1099. RES = ETG RES ;
  1100. *
  1101. * RÉSOLUTION SÉQUENTIELLE
  1102. SINON ;
  1103. SI IMPARA ;
  1104. MESS 'Les résolutions sont traitées séquentiellement' ;
  1105. FINS ;
  1106. *
  1107. RES MM1 = KOPS 'MATRIK' ;
  1108. $TRICO = RV.'$trico' ;
  1109. REPE BCLCOM NBPART ;
  1110. SI (($TAB1.&BCLCOM.'TYPINV' EGA 1) ET TRESOU) ;
  1111. MESS 'ILU-Crout => on utilise RESOU' ;
  1112. *
  1113. MATCOM = $MA1.&BCLCOM ;
  1114. MATCOM = KOPS 'RIMA' MATCOM 'NSYM' ;
  1115. *
  1116. * RIGIDITE/CHPOINT associés aux C.L.
  1117. SI (NON (EXIS ($TRICO.&BCLCOM) &BINC)) ;
  1118. INC = EXTR ($S1.&BCLCOM) 'COMP' ;
  1119. NCMP = DIME INC ;
  1120. REPE BINC NCMP ;
  1121. ICO = EXTR &BINC INC ;
  1122. OPTI 'INCO' ICO ICO ;
  1123. RCLIM = EXCO ICO ($S1.&BCLCOM) ICO ;
  1124. MIP = EXTR RCLIM 'MAILLAGE' ;
  1125. $TRICO.&BCLCOM.(&BINC) = BLOQ MIP ICO ;
  1126. FIN BINC ;
  1127. FINS ;
  1128. S1 = VIDE 'CHPOINT'/'DISCRET' ;
  1129. REPE BINC NCMP ;
  1130. ICO = EXTR &BINC INC ;
  1131. RCLIM = EXCO ICO ($S1.&BCLCOM) ICO ;
  1132. RIGP = $TRICO.(&BCLCOM).&BINC ;
  1133. S1 = S1 ET (DEPI RIGP RCLIM) ;
  1134. FIN BINC ;
  1135. *
  1136. RIGP = ETG ($TRICO.&BCLCOM) ;
  1137. MATCOM = MATCOM ET RIGP ;
  1138. *
  1139. RES1 = RESO MATCOM ($S2.&BCLCOM ET S1) ;
  1140. RES1 = RES1 'CHAN' 'ATTRIBUT' 'NATURE' 'DISCRET' ;
  1141.  
  1142. SINON ;
  1143. RES1 = KRES ($MA1.&BCLCOM) 'TYPI' ($TAB1.&BCLCOM)
  1144. 'CLIM' ($S1.&BCLCOM)
  1145. 'SMBR' ($S2.&BCLCOM)
  1146. 'IMPR' IMPKRES ;
  1147. FINS ;
  1148. RES = RES ET RES1 ;
  1149. FIN BCLCOM ;
  1150. FINS ;
  1151. FINS ;
  1152. *
  1153. *
  1154. * ===============================================
  1155. * ÉTAPE DE PROJECTION STANDARD : CALCUL DE C*Û[n]
  1156. * ===============================================
  1157. *
  1158. SI TESTPRJ ;
  1159. *
  1160. CUN = KMF (RVP.'MATC') RES ;
  1161. SI TVNP ;
  1162. CXN = KMF (RVP.'MBTR') RES ;
  1163. CUN = CUN ET CXN ;
  1164. FINS ;
  1165. *
  1166. * On calcule les seconds membres de l'équation de pression
  1167. * s'ils existent (opérateurs FIMP)
  1168. REPE BLOCPJ (DIME (RVP.'LISTOPER')) ;
  1169. NOMPER = EXTR &BLOCPJ (RVP.'LISTOPER') ;
  1170. NOMTAB = CHAI &BLOCPJ NOMPER ;
  1171. SI (EGA NOMPER 'FIMP') ;
  1172. MSI MAI = (TEXTE NOMPER) (RVP.NOMTAB) ;
  1173. CUN = CUN ET MSI ;
  1174. FINS ;
  1175. FIN BLOCPJ ;
  1176. *
  1177. CUN = CUN * ROW * (-1./DT) ;
  1178. *
  1179. * Calcul de Lamdba[n]
  1180. RVP.'METHINV'.'XINIT' = RV.'resmnp' ;
  1181. SL = KRESP MATPR 'TYPI' (RVP.'METHINV')
  1182. 'CLIM' (RVP.'CLIM' )
  1183. 'SMBR' CUN
  1184. 'IMPR' IMPKRES ;
  1185. OUBL CUN ;
  1186. RV.'resmnp' = SL ;
  1187. *
  1188. CTL = KMF (RVP.'MATC') SL 'TRAN' ;
  1189. SI TVNP ;
  1190. CXL = KMF (RVP.'MBTR') SL 'TRAN' ;
  1191. CTL = CTL + CXL ;
  1192. FINS ;
  1193. *
  1194. * Elimination of end of step velocity
  1195. * (si 'PNM2' existe, on ne corrige pas)
  1196. SI (NON TPNM2) ;
  1197. SI (NON (EXIS (RV.'INCO') 'PNM2')) ;
  1198. A = NOMC LCU LC (KOPS IDIGV '*' CTL) ;
  1199. RES = RES + ((1./ROW)*A*DT) ;
  1200. FINS ;
  1201. FINS ;
  1202. OUBL CTL ;
  1203. *
  1204. FINS ;
  1205. *
  1206. *
  1207. * ===================
  1208. * AVANCEMENT EN TEMPS
  1209. * ===================
  1210. *
  1211. * SI (TESTPRJ ET (EGA MDFDT 0)) ;
  1212. * IMPTCRR = IIMPR ;
  1213. * FINS ;
  1214. SI (EGA MDFDT 0) ;
  1215. IMPTCRR = IIMPR ;
  1216. FINS ;
  1217. IMPKRES = 0 ;
  1218. *
  1219. * MESS 'On passe dans TCRR ' ;
  1220. LCMP = EXTR RES 'COMP' ;
  1221. SI (EXIS LCMP 'LX') ;
  1222. * AG = EXTR (EXCO RES 'LX') 'MAILLAGE' ;
  1223. * AR = EXTR RES 'MAILLAGE' ;
  1224. * ATN = DIFF AR AG ;
  1225. * RES = REDU RES ATN ;
  1226. RES = ENLE RES 'LX' ;
  1227. FINS ;
  1228. *
  1229. EPS = TCRR RES OMEG (RV.'INCO') 'IMPR' IMPTCRR ;
  1230. *
  1231. RV.'METHINV'.'CALPREC' = FAUX ;
  1232. SI (TESTPRJ) ;
  1233. RVP.'METHINV'.'CALPREC'= FAUX ;
  1234. FINS ;
  1235. *
  1236. RV.'resmn' = RES ;
  1237. *
  1238. MENA ;
  1239. *
  1240. SI (RV.'STOPITER') ;
  1241. RV.'STOPITER' = FAUX ;
  1242. QUIT BLOCI ;
  1243. FINS ;
  1244. *
  1245. FIN BLOCI ;
  1246. * ******************************************
  1247. * ******************************************
  1248. * FIN DE LA BOUCLE INTERNE (= NON-LINEARITE)
  1249. * ******************************************
  1250. * ******************************************
  1251. *
  1252. *
  1253. * MESS ' ITMAX= ' ITMAX ' MDFDT= ' MDFDT ;
  1254. *
  1255. IRT=0 ;
  1256. SI (EGA ITMAX 0) ;
  1257. IRT = TCNM RV 'NOUP' ;
  1258. SINON ;
  1259. IRT = TCNM RV ;
  1260. FINS ;
  1261. *
  1262. * ------------------------------------------------------------
  1263. * Avancement en temps pour l'algorithme de projection standard
  1264. * => P[n+1] = P[n] + Lambda
  1265. * ------------------------------------------------------------
  1266. SI TESTPRJ ;
  1267. SI (NON (EXIS (RV.'INCO') 'PRESSION')) ;
  1268. RV.'INCO'.'PRESSION' = SL ;
  1269. SI TPNM2 ;
  1270. RV.'INCO'.'PNM2' = SL ;
  1271. FINS ;
  1272. SINON ;
  1273. PNM1 = RV.'INCO'.'PRESSION' ;
  1274. SI (EGA ISCHT 1) ;
  1275. PN = PNM1 + (1.5*SL) ;
  1276. SINON ;
  1277. PN = PNM1 + SL ;
  1278. FINS ;
  1279.  
  1280. SI TPNM2 ;
  1281. RV.'INCO'.'PNM2' = PNM1 ;
  1282. FINS ;
  1283. RV.'INCO'.'PRESSION' = PN ;
  1284. FINS ;
  1285.  
  1286. SI ((EGA TYPROJ 'PSCT') OU (EGA TYPROJ 'PENA')) ;
  1287. RV.'INCO'.'PRESSION' = SL ;
  1288. FINS ;
  1289.  
  1290. OUBL SL ;
  1291. FINS ;
  1292. *
  1293. * ----------------------------------------
  1294. * Traitement spécifique (non documenté...)
  1295. * ----------------------------------------
  1296. SI TESTRAN ;
  1297. RV.'CO'.'VITESSE' = KOPS (RV.'INCO'.NOMVI) '-' (RV.'SEDIM') ;
  1298. K = 'ABS' (RV.'INCO'.'KN') ;
  1299. E = 'ABS' (RV.'INCO'.'EN') ;
  1300. K = KOPS (KOPS K '*' K) '/' E ;
  1301. DIF = KOPS (KOPS K '*' 0.09) '+' (RV.'COEF') ;
  1302. RV.'CO'.'DIFFU' = NOEL RVPD DIF ;
  1303. RV.'CO'.'TEMPERA' = RV.'INCO'.'CN' ;
  1304. FINS ;
  1305. *
  1306. *
  1307. *
  1308. RV.'METHINV'.'CALPREC' = FAUX ;
  1309. SI (TESTPRJ) ;
  1310. RVP.'METHINV'.'CALPREC'= FAUX ;
  1311. FINS ;
  1312. *
  1313. MENA ;
  1314. *
  1315. SI (EGA IRT 1) ;
  1316. MESS 'TEMPS FINAL ATTEINT : ' (RV.'PASDETPS'.'TPS') ;
  1317. QUIT BLOC1 ;
  1318. FINS ;
  1319. *
  1320. SI (RV.'STOPPDT') ;
  1321. RV.'STOPPDT' = FAUX ;
  1322. MESS 'ARRET DEMANDÉ AU TEMPS : ' (RV.'PASDETPS'.'TPS') ;
  1323. QUIT BLOC1 ;
  1324. FINS ;
  1325. *
  1326. FIN BLOC1 ;
  1327. * *****************************************
  1328. * *****************************************
  1329. * FIN DE LA BOUCLE EXTERNE (= PAS DE TEMPS)
  1330. * *****************************************
  1331. * *****************************************
  1332. *
  1333. *
  1334. SI TESTPR ;
  1335. RVP.'PRESSION' = KCHT RVPD 'SCAL' 'CENTRE' (RVP.'PRESSION') ;
  1336. RVP.'PN' = ELNO RVPD (RVP.'PRESSION') ;
  1337. FINS ;
  1338. *
  1339. * OUBLI EVENTUEL des matrices présentes dans la table RV et qui prennent de
  1340. * la place dans les sauvegardes
  1341. *
  1342. 'SI' ('EXIS' rv 'DETMAT') ;
  1343. ldetmat = rv . 'DETMAT' ;
  1344. 'SINO' ;
  1345. ldetmat = faux ;
  1346. 'FINS' ;
  1347.  
  1348. 'SI' ldetmat ;
  1349. REPIX rv ;
  1350. 'FINS' ;
  1351.  
  1352.  
  1353.  
  1354. *
  1355. FINP ;
  1356. *
  1357.  
  1358.  

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