Télécharger la plaquette Cast3M


fr en

Télécharger nlin_te_unstat.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  nlin_te_unstat.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *************************************************************************
  5. *
  6. * NOM         : nlin_te_unstat.dgibi
  7. *
  8. * DESCRIPTION : We compute the flow governed by the Navier-Stokes
  9. *               equations, in a T jonction with high temperature
  10. *               difference.
  11. *               Both Boussinesq approximation and Low Mach number
  12. *               approximation are considered.
  13. *               See also technical report SFME/LTMF/RT/O5-067/A
  14. *               In this test case we simply perform 2 time iterations
  15. *               and we verify that two point-fix algorithms converge
  16. *               (to the same solution).
  17. *
  18. * LANGAGE     : GIBIANE-CAST3M
  19. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  20. *               mél : beccantini@cea.fr
  21. *************************************************************************
  22. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  23. * en cas de modification de ce sous-programme afin de faciliter
  24. * la maintenance !
  25. *************************************************************************
  26. *
  27. *
  28. *************************************************
  29. ****** LIST OF SYMBOLS **************************
  30. *************************************************
  31. *
  32. * \int  = integral
  33. * \div  = divergence
  34. * \der  = derivative (\der{a}{b} = \frac{d a}{d b} 
  35. * \grad = gradient
  36. * \cdot = scalar product
  37. * ^     = power
  38. *
  39. * (a,b) = \int_{\Omega} a * b dV
  40. *
  41. *  Np   = test function for the divergence equations
  42. *       = shape function for the pressure
  43. *  Nt   = test function for the temperature
  44. *       = shape function for the temperature equation 
  45. *  Nv   = test function for the speed equations (x,y)
  46. *       = shape function for the speed
  47. *
  48. *************************************************************************
  49. *************************************************************************
  50. ******* PERSONAL PROCEDURES *********************************************
  51. *************************************************************************
  52. *************************************************************************
  53. *
  54. *BEGINPROCEDUR matmas
  55. *
  56. ************************************************************************
  57. *
  58. * NOM         : MATMAS
  59. *
  60. * DESCRIPTION : We compute 
  61. *
  62. *               (Np , Nd )
  63. *
  64. * LANGAGE     : GIBIANE-CAST3M
  65. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  66. *               mél : beccantini@cea.fr
  67. **********************************************************************
  68. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  69. * en cas de modification de ce sous-programme afin de faciliter
  70. * la maintenance !
  71. ************************************************************************
  72. *
  73. *
  74. 'DEBPROC' MATMAS ;
  75. *
  76. 'ARGUMENT'  _mt*'MAILLAGE'    ;
  77. 'ARGUMENT'  nomp*'MOT     '  ;
  78. 'ARGUMENT'  nomd*'MOT     '  ;
  79. 'ARGUMENT'  discp*'MOT     '  ;
  80. 'ARGUMENT'  discd*'MOT     ' ;
  81. *
  82. *  Arguments
  83. *  _mt   = MAILLAGE, quaf mesh
  84. *  nomp  = MOT, name of the primal variable
  85. *  nomd  = MOT, name of the dual variable
  86. *  discp = MOT, name of the discretization of the primal variable
  87. *  discd = MOT, name of the discretization of the dual variable
  88. *
  89. idim   = 'VALEUR' 'DIME' ;
  90. *
  91. *   Matrix A for the operator NLIN
  92. *
  93. numop  = 1 ;
  94. numder = idim ;
  95. numvar = 1 ;
  96. * The variable involved in A is the primal variable
  97. numdat = 0 ;
  98. numcof = 0 ;
  99. *
  100. discg = 'LINE' ;
  101. * Geometrical discretisation
  102. methgau = 'GAU7' ;
  103. * Gauss points involved
  104. *
  105. A = ININLIN numop numvar numdat numcof numder ;   
  106. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomp ;
  107. A . 'VAR' . 1 . 'DISC'   = discp ;
  108. *
  109. *   Matrix B for the operator NLIN
  110. *
  111. B = ININLIN numop numvar 0 0 numder ;
  112. * For B the coef no coefficients are involved
  113. * but 1
  114. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
  115. B . 'VAR' . 1 . 'DISC'   = discd ;
  116. *
  117. A . 1 . 1 . 0 = 'LECT' ;
  118. B . 1 . 1 . 0 = 'LECT'  ;
  119. *
  120. mat = 'NLIN' discg _mt A B methgau ;
  121. 'RESPRO' mat ;
  122. 'FINPROC' ;
  123. *
  124. *ENDPROCEDUR matmas
  125. *BEGINPROCEDUR matdiv
  126. *
  127. ************************************************************************
  128. *
  129. * NOM         : MATDIV
  130. *
  131. * DESCRIPTION : We compute the matrix arising from the computation of
  132. *               the scalar product
  133. *
  134. *               (div(u),Np)
  135. *
  136. *               where Np = shape function for pressure
  137. *
  138. *
  139. * LANGAGE     : GIBIANE-CAST3M
  140. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  141. *               mél : beccantini@cea.fr
  142. **********************************************************************
  143. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  144. * en cas de modification de ce sous-programme afin de faciliter
  145. * la maintenance !
  146. ************************************************************************
  147. *
  148. *
  149. 'DEBPROC' MATDIV ;
  150. *
  151. 'ARGUMENT'   _mt*'MAILLAGE' ;
  152. 'ARGUMENT'  nompre*'MOT     ' ;
  153. 'ARGUMENT'  discp*'MOT     ' ;
  154. 'ARGUMENT'  discv*'MOT     ' ;
  155. *
  156. *  _mt   = solid QUAF mesh
  157. *  nompre = name of the pressure
  158. *  discp = discretization type for the pressure
  159. *          (for instance 'LINE')
  160. *  discv = discretization type for the speed
  161. *          (for instance 'QUAF')
  162. *
  163. idim   = 'VALEUR' 'DIME' ;
  164. *
  165. *    Two contributions
  166. *
  167. *     \dep{u_x}{x} Np
  168. *   + \dep{u_y}{y} Np
  169. *
  170. *   Matrix A for the operator NLIN
  171. *
  172. numop  = 1 ;
  173. numder = idim ;
  174. numvar = 2 ;
  175. * The variables involved in A are u_x, u_y
  176. numdat = 0 ;
  177. * Zero data
  178. numcof = 0 ;
  179. * Zero functions (but 1)
  180. discg = 'LINE' ;
  181. * Geometrical discretisation
  182. methgau = 'GAU7' ;
  183. * Gauss points involved
  184. *
  185. A = ININLIN numop numvar numdat numcof numder ;   
  186. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  187. A . 'VAR' . 1 . 'DISC'   = discv ;
  188. A . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  189. A . 'VAR' . 2 . 'DISC'   = discv ;
  190. *
  191. *   Matrix B for the operator NLIN
  192. *
  193. numvar = 1 ;
  194. B = ININLIN numop numvar 0 0 numder ;
  195. * For B the coef no coefficients are involved
  196. * but 1
  197. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nompre ;
  198. B . 'VAR' . 1 . 'DISC'   = discp ;
  199. *
  200. * Contribution
  201. *
  202. * 1 \dep{u_x}{x} Np
  203. A . 1 . 1 . 1 = 'LECT'  ;
  204. B . 1 . 1 . 0 = 'LECT'  ;
  205. * 1 \dep{u_y}{y} Np
  206. A . 1 . 2 . 2 = 'LECT'  ;
  207. *
  208. mdiv = 'NLIN' discg _mt A B  methgau ;
  209. 'RESPRO' mdiv ;
  210. 'FINPROC' ;
  211. *
  212. *ENDPROCEDUR matdiv
  213. *BEGINPROCEDUR matflu
  214. *
  215. ************************************************************************
  216. *
  217. * NOM         : MATFLU
  218. *
  219. * DESCRIPTION : We compute the matrix arising from
  220. *
  221. *               \int_{\delta \Omega} (coef q \cdot n) Nd dS
  222. *
  223. *               where q is the primal variable ('UX', 'UY')
  224. *               Nd is the shape function of the dual variable
  225. *
  226. * LANGAGE     : GIBIANE-CAST3M
  227. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  228. *               mél : beccantini@cea.fr
  229. **********************************************************************
  230. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  231. * en cas de modification de ce sous-programme afin de faciliter
  232. * la maintenance !
  233. ************************************************************************
  234. *
  235. *
  236. **** We compute the matrix
  237. *
  238. *    \int_{\delta \Omega} (q \cdot n) Nd dS
  239. *
  240. *    N.B. Be carefull in specifying n (the normal at the surfaces)
  241. *         No singular points are aviable.
  242. *
  243. *    Arguments
  244. *
  245. *    _st    = MAILLAGE, quaf mesh of the surface
  246. *    discp  = MOT, name of the discretisation of the primal variable
  247. *    nomd   = MOT, name of the dual variable
  248. *    discd  = MOT, name of the discretisation of the dual variable
  249. *    disco  = MOT, discretisation of the coefficient
  250. *    coef   = coefficient
  251. *    discn  = MOT, name of the discretisation of the normals
  252. *    nx     = CHPOINT 
  253. *    ny     = CHPOINT 
  254. *
  255. *    nomp   = name of the primal variable are 'UX' 'UY'
  256. *
  257. 'DEBPROC' MATFLU ;
  258. 'ARGUMENT'  _st*'MAILLAGE' ;
  259. 'ARGUMENT'  discp*'MOT     ' ;
  260. 'ARGUMENT'  nomd*'MOT     ' ;
  261. 'ARGUMENT'  discd*'MOT     ' ;
  262. 'ARGUMENT'  discco*'MOT     ' ;
  263. 'ARGUMENT'  coef*'CHPOINT' ;
  264. 'ARGUMENT'  discn*'MOT     ' ;
  265. 'ARGUMENT'  nx*'CHPOINT' ;
  266. 'ARGUMENT'  ny*'CHPOINT' ;
  267. *
  268. idim   = 'VALEUR' 'DIME' ;
  269. *
  270. *   Matrix A for the operator NLIN
  271. *
  272. numop  = 1 ;
  273. numder = idim ;
  274. numvar = 2 ;
  275. * The variables involved in A are qx, qy
  276. numdat = 3 ;
  277. * coef, nx, ny
  278. numcof = 3 ;
  279. * coef, nx, ny ('IDEN')
  280. discg = 'LINE' ;
  281. * Geometrical discretisation
  282. methgau = 'GAU7' ;
  283. * Gauss points involved
  284. *
  285. A = ININLIN numop numvar numdat numcof numder ;   
  286. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  287. A . 'VAR' . 1 . 'DISC'   = discp ;
  288. A . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  289. A . 'VAR' . 2 . 'DISC'   = discp ;
  290. *
  291. * coef
  292. A . 'DAT' . 1 . 'NOMDDL' = 'EXTRAIRE' coef 'COMP';
  293. A . 'DAT' . 1 . 'DISC' =  discco ;
  294. A . 'DAT' . 1 . 'VALEUR' = coef ;
  295. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  296. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  297. * nx 
  298. A . 'DAT' . 2 . 'NOMDDL' = 'EXTRAIRE' nx 'COMP';
  299. A . 'DAT' . 2 . 'DISC' =  discn ;
  300. A . 'DAT' . 2 . 'VALEUR' = nx ;
  301. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  302. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  303. * ny
  304. A . 'DAT' . 3 . 'NOMDDL' = 'EXTRAIRE' ny 'COMP';
  305. A . 'DAT' . 3 . 'DISC' =  discn ;
  306. A . 'DAT' . 3 . 'VALEUR' = ny ;
  307. A . 'COF' . 3 . 'COMPOR' = 'IDEN' ;
  308. A . 'COF' . 3 . 'LDAT' = 'LECT' 3 ;
  309. *
  310. *   Matrix B for the operator NLIN
  311. *
  312. numvar =1 ;
  313. B = ININLIN numop numvar 0 0 numder ;
  314. *
  315. * For B no coefficients are involved
  316. * but 1
  317. *
  318. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
  319. B . 'VAR' . 1 . 'DISC'   = discd ;
  320. * (coef * nx * qx '+' coef * ny qy) Nd
  321. A . 1 . 1 . 0 = 'LECT' 1 2;
  322. A . 1 . 2 . 0 = 'LECT' 1 3;
  323. B . 1 . 1 . 0 = 'LECT'  ;
  324. *
  325. mat = 'NLIN' discg _st A B methgau ;
  326. 'RESPRO' mat ;
  327. 'FINPROC' ;
  328. *
  329. *ENDPROCEDUR matflu
  330. *BEGINPROCEDUR chpflu
  331. *
  332. ************************************************************************
  333. *
  334. * NOM         : CHPFLU
  335. *
  336. * DESCRIPTION : We compute the CHPOIN
  337. *               
  338. *               -q = lambda (\grad(T) 
  339. *
  340. *               on the dof of the dual variable
  341. *               T  = CHPOINT
  342. *
  343. *               N.B. Be carefull in specifying n (the normal at the
  344. *                    surfaces).
  345. *                    No singular points are admitted.
  346. *
  347. * LANGAGE     : GIBIANE-CAST3M
  348. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  349. *               mél : beccantini@cea.fr
  350. **********************************************************************
  351. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  352. * en cas de modification de ce sous-programme afin de faciliter
  353. * la maintenance !
  354. ************************************************************************
  355. *
  356. *
  357. 'DEBPROC' CHPFLU ;
  358. *
  359. 'ARGUMENT'  _mt*'MAILLAGE' ;
  360. 'ARGUMENT'  disct*'MOT     ' ;
  361. 'ARGUMENT'  T*'CHPOINT'     ;
  362. 'ARGUMENT'  discd*'MOT     ' ;
  363. 'ARGUMENT'  discla*'MOT     ' ;
  364. 'ARGUMENT'  lambda*'CHPOINT' ;
  365. *
  366. *  _mt    = MAILLAGE, solid quaf mesh
  367. *  disct  = MOT, name of the discretisation of T
  368. *  T      = CHPOINT (one component, the name of which is free)
  369. *  discd  = MOT, name of the discretisation of the dual variable
  370. *  discla = MOT, name of the discretisation of lambda
  371. *  lambda = CHPOINT (one component, SCAL)
  372. *
  373. idim   = 'VALEUR' 'DIME' ;
  374. *
  375. *
  376. *  Matrix A for the operator NLIN
  377. *
  378. numop  = 2 ;
  379. numder = idim ;
  380. numvar = 1 ;
  381. * One variable, T
  382. numdat = 1 ;
  383. * -lambda
  384. numcof = 1 ;
  385. * Coef: -lambda (three identity functions)
  386. discg = 'LINE' ;
  387. * Geometrical discretisation
  388. methgau = 'GAU7' ;
  389. * Gauss points involved
  390. *
  391. *   Matrix A for the operator NLIN
  392. *
  393. A = ININLIN numop numvar numdat numcof numder ;   
  394. A . 'VAR' . 1 . 'NOMDDL' = 'EXTR' TN 'COMP' ;
  395. A . 'VAR' . 1 . 'DISC'   = disct ;
  396. * lambda
  397. A . 'DAT' . 1 . 'NOMDDL' = ('EXTRAIRE' lambda 'COMP')  ;
  398. A . 'DAT' . 1 . 'DISC' =  discla ;
  399. A . 'DAT' . 1 . 'VALEUR' = lambda ;
  400. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  401. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  402. * 
  403. *   Matrix B for the operator NLIN
  404. *
  405. numvar = 2 ;
  406. B = ININLIN numop numvar 0 0 numder ;
  407. *
  408. * For B, no coefficients are involved
  409. * but 1
  410. *
  411. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  412. B . 'VAR' . 1 . 'DISC'   = discd ;
  413. B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  414. B . 'VAR' . 2 . 'DISC'   = discd ;
  415. *
  416. * Contribution
  417. *
  418. * lambda \dep{T}{x} Nd
  419. A . 1 . 1 . 1 = 'LECT' 1  ;
  420. B . 1 . 1 . 0 = 'LECT'  ;
  421. * lambda \dep{T}{y}  Nd
  422. A . 2 . 1 . 2 = 'LECT' 1  ;
  423. B . 2 . 2 . 0 = 'LECT'  ;
  424. *
  425. mat = 'NLIN' discg _mt A B  methgau ;
  426. *
  427. * We compute coef
  428. *
  429.  matm1 = MATMAS _mt 'UX' 'UX' discd discd ;
  430.  matm2 = MATMAS _mt 'UY' 'UY' discd discd ;
  431.  matm = matm1 'ET' matm2 ;
  432.  mtotik = 'KOPS' 'RIMA' matm ;
  433.  smbt = mat '*' t ;
  434.  chlim pipi = 'KOPS' matrik ;
  435.  'MESSAGE' ;
  436.  coef = RESOUP FAUX mtotik mtotik smbt chlim 1.0D-12 200 'ITER' ;
  437. *
  438. 'RESPRO' coef ;
  439. 'FINPROC' ;
  440. *
  441. *ENDPROCEDUR chpflu
  442. *BEGINPROCEDUR matkt
  443. *
  444. ************************************************************************
  445. *
  446. * NOM         : MATKT
  447. *
  448. * DESCRIPTION : We compute the matrix
  449. *               
  450. *               ((u_x \dep{T}{x} '+' u_y \dep{T}{y}) , Nt)
  451. *
  452. * LANGAGE     : GIBIANE-CAST3M
  453. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  454. *               mél : beccantini@cea.fr
  455. **********************************************************************
  456. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  457. * en cas de modification de ce sous-programme afin de faciliter
  458. * la maintenance !
  459. ************************************************************************
  460. *
  461. *
  462. 'DEBPROC' MATKT ;
  463. *
  464. 'ARGUMENT'   _mt*'MAILLAGE' ;
  465. 'ARGUMENT'  discv*'MOT     ' ;
  466. 'ARGUMENT'    ux*'CHPOINT' ;
  467. 'ARGUMENT'    uy*'CHPOINT' ;
  468. 'ARGUMENT'  nomt *'MOT     ' ;
  469. 'ARGUMENT'  disct*'MOT     ' ;
  470. *
  471. *    discv = discretisation type for u_x and u_y
  472. *    ux    = chpoint ('SCAL', on discv)
  473. *    uy    = chpoint ('SCAL', on discv)
  474. *    nomt  = name of the primal and the dual variable t
  475. *    disct = discretisation type for t
  476. *
  477. idim   = 'VALEUR' 'DIME' ;
  478. *
  479. numop  = 1 ;
  480. numder = idim ;
  481. numvar = 1 ;
  482. * One variable, T
  483. numdat = 2 ;
  484. * ux, uy
  485. numcof = 2 ;
  486. * Coef: ux  uy (two identity functions)
  487. discg = 'LINE' ;
  488. * Geometrical discretisation
  489. methgau = 'GAU7' ;
  490. * Gauss points involved
  491. *
  492. *   Matrix A for the operator NLIN
  493. *
  494. A = ININLIN numop numvar numdat numcof numder ;   
  495. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomt ;
  496. A . 'VAR' . 1 . 'DISC'   = disct ;
  497. * ux
  498. A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  499. A . 'DAT' . 1 . 'DISC' =  discv ;
  500. A . 'DAT' . 1 . 'VALEUR' = ux ;
  501. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  502. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  503. * uy 
  504. A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  505. A . 'DAT' . 2 . 'DISC' =  discv ;
  506. A . 'DAT' . 2 . 'VALEUR' = uy ;
  507. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  508. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  509. * 
  510. *   Matrix B for the operator NLIN
  511. *
  512. numvar = 1 ;
  513. B = ININLIN numop numvar 0 0 numder ;
  514. *
  515. * For B, no coefficients are involved
  516. * but 1
  517. *
  518. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomt ;
  519. B . 'VAR' . 1 . 'DISC'   = disct ;
  520. *
  521. * Contribution
  522. *
  523. * ux \dep{T}{x} Nd
  524. A . 1 . 1 . 1 = 'LECT' 1 ;
  525. B . 1 . 1 . 0 = 'LECT'   ;
  526. * uy \dep{T}{y}  Nd
  527. A . 1 . 1 . 2 = 'LECT' 2 ;
  528. *
  529. mat = 'NLIN' discg _mt A B  methgau ;
  530. 'RESPRO' mat ;
  531. 'FINPROC' ;
  532. *
  533. *ENDPROCEDUR matkt
  534. *BEGINPROCEDUR matlap
  535. *
  536. ************************************************************************
  537. *
  538. * NOM         : MATLAP
  539. *
  540. * DESCRIPTION : We compute the matrix arising from the computation of
  541. *               the scalar product
  542. *
  543. *               ( -(coef \div( \alpha \grad T)) , Nd )
  544. *
  545. *               where Nd = test function for the dual variable.
  546. *
  547. *               -\int_{\Omega} coef \div( \alpha \grad T) Nd dV =
  548. *                  -\int_{\delta \Omega}
  549. *                     coef \alpha (\grad T \cdot n) Nd dS '+'
  550. *                  +\int_{\Omega}
  551. *                     coef \alpha (\grad T \cdot \grad Nd) dV '+'
  552. *                  +\int_{\Omega}
  553. *                     Nd \alpha (\grad T \cdot \grad coef) dV
  554. *
  555. *               The matrix is issue from the 2nd and 3rd contribution
  556. *               (volume integrals).
  557. *
  558. *
  559. * LANGAGE     : GIBIANE-CAST3M
  560. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  561. *               mél : beccantini@cea.fr
  562. **********************************************************************
  563. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  564. * en cas de modification de ce sous-programme afin de faciliter
  565. * la maintenance !
  566. ************************************************************************
  567. *
  568. *
  569. 'DEBPROC' MATLAP ;
  570. *
  571. 'ARGUMENT'  _mt*'MAILLAGE'    ;
  572. 'ARGUMENT'   nomt*'MOT     '  ;
  573. 'ARGUMENT'  disct*'MOT     '  ;
  574. 'ARGUMENT'   nomd*'MOT     '  ;
  575. 'ARGUMENT'  discd*'MOT     '  ;
  576. 'ARGUMENT'  discal*'MOT     ' ;
  577. 'ARGUMENT'  alpha*'CHPOINT'   ;
  578. 'ARGUMENT'  discof*'MOT     ' ;
  579. 'ARGUMENT'  coef*'CHPOINT'    ;
  580. *
  581. *    Arguments
  582. *
  583. *    _mt    = surface QUAF mesh
  584. *    nomt   = name of the primal variable T
  585. *    disct  = discretization type of T
  586. *    nomd   = name of the dual variable
  587. *    discd  = discretization type of the dual variable
  588. *    discal = discretization type of alpha
  589. *    discof = discretization type of coef
  590. *    alpha  = CHPOINT ('SCAL')
  591. *    coef   = CHPOINT ('SCAL')
  592. *
  593. idim   = 'VALEUR' 'DIME' ;
  594. *
  595. *    Four contributions
  596. *
  597. *    (coef alpha \dep{T}{x} \dep{Nd}{x}
  598. *   + coef alpha \dep{T}{y} \dep{Nd}{y}
  599. *   + Nd alpha \dep{T}{x} \dep{coef}{x}
  600. *   + Nd alpha \dep{T}{y} \dep{coef}{y})
  601. *
  602. *   Matrix A for the operator NLIN
  603. *
  604. numop  = 4 ;
  605. numder = idim ;
  606. numvar = 1 ;
  607. * The variable involved in A is the primal variable T
  608. numdat = 2 ;
  609. * Two data: alpha, coef
  610. numcof = 4 ;
  611. * Four functions: f(alpha) = alpha
  612. *                 f(coef)  = coef
  613. *                 f(coef) = \dep{coef}{x}
  614. *                 f(coef) = \dep{coef}{y}
  615. discg = 'LINE' ;
  616. * Geometrical discretisation
  617. methgau = 'GAU7' ;
  618. * Gauss points involved
  619. *
  620. A = ININLIN numop numvar numdat numcof numder ;   
  621. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomt ;
  622. A . 'VAR' . 1 . 'DISC'   = disct ;
  623. * alpha
  624. A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  625. A . 'DAT' . 1 . 'DISC' = discal ;
  626. A . 'DAT' . 1 . 'VALEUR' = alpha ;
  627. * coef
  628. A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  629. A . 'DAT' . 2 . 'DISC' = discof ;
  630. A . 'DAT' . 2 . 'VALEUR' = coef ;
  631. * Function alpha
  632. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  633. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  634. * Function coef
  635. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  636. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  637. * Function \dep{coef}{x}
  638. A . 'COF' . 3 . 'COMPOR' = 'D/DX1' ;
  639. A . 'COF' . 3 . 'LDAT' = 'LECT' 2  ;
  640. * Function \dep{coef}{y}
  641. A . 'COF' . 4 . 'COMPOR' = 'D/DX2' ;
  642. A . 'COF' . 4 . 'LDAT' = 'LECT' 2 ;
  643. *
  644. *   Matrix B for the operator NLIN
  645. *
  646. B = ININLIN numop numvar 0 0 numder ;
  647. * For B the coef no coefficients are involved
  648. * but 1
  649. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
  650. B . 'VAR' . 1 . 'DISC'   = discd ;
  651. *
  652. * Contribution
  653. *
  654. * coef alpha \dep{T}{x} \dep{Nd}{x}
  655. A . 1 . 1 . 1 = 'LECT' 1 2 ;
  656. B . 1 . 1 . 1 = 'LECT'  ;
  657. * coef alpha \dep{T}{y} \dep{Nd}{y}
  658. A . 2 . 1 . 2 = 'LECT' 1 2 ;
  659. B . 2 . 1 . 2 = 'LECT'  ;
  660. * Nd alpha \dep{T}{x} \dep{coef}{x}
  661. A . 3 . 1 . 1 = 'LECT' 1 3 ;
  662. B . 3 . 1 . 0 = 'LECT'  ;
  663. * Nd alpha \dep{T}{y} \dep{coef}{y}
  664. A . 4 . 1 . 2 = 'LECT' 1 4 ;
  665. B . 4 . 1 . 0 = 'LECT'  ;
  666. mlapn = 'NLIN' discg _mt A B  methgau ;
  667. 'RESPRO' mlapn ;
  668. 'FINPROC' ;
  669. *ENDPROCEDUR matlap
  670. *BEGINPROCEDUR matkv
  671. *
  672. ************************************************************************
  673. *
  674. * NOM         : MATKV
  675. *
  676. * DESCRIPTION : We compute the matrix
  677. *
  678. *               ((cx \dep{u_x}{x} '+' cy \dep{u_x}{y}), Nvx)
  679. *               ((cx \dep{u_y}{x} '+' cy \dep{u_y}{y}), Nvy)
  680. *
  681. *               where Nv = shape function for (u_x,u_y)^T
  682. *
  683. *               Names of primal variables = names of the dual variables
  684. *                = 'UX' 'UY'
  685. *
  686. * LANGAGE     : GIBIANE-CAST3M
  687. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  688. *               mél : beccantini@cea.fr
  689. **********************************************************************
  690. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  691. * en cas de modification de ce sous-programme afin de faciliter
  692. * la maintenance !
  693. ************************************************************************
  694. *
  695. 'DEBPROC' MATKV ;
  696. 'ARGUMENT'   _mt*'MAILLAGE' ;
  697. 'ARGUMENT' discv*'MOT     ' ;
  698. 'ARGUMENT' discc*'MOT     ' ;
  699. 'ARGUMENT'     cx*'CHPOINT' ;
  700. 'ARGUMENT'     cy*'CHPOINT' ;
  701. *
  702. *  _mt   = solid quaf mesh
  703. *  discv = discretization type ux, uy
  704. *  discc = discretization type for cx, cy
  705. *  cx    = CHPOINT ('SCAL')
  706. *  cy    = CHPOINT ('SCAL')
  707. *
  708. idim   = 'VALEUR' 'DIME' ;
  709. *
  710. numop  = 2 ;
  711. numder = idim ;
  712. numvar = 2 ;
  713. * Two variables, ux, uy
  714. numdat = 2 ;
  715. * cx, cy
  716. numcof = 2 ;
  717. * Coef: cx  cy (two identity functions)
  718. discg = 'LINE' ;
  719. * Geometrical discretisation
  720. methgau = 'GAU7' ;
  721. * Gauss points involved
  722. *
  723. *   Matrix A for the operator NLIN
  724. *
  725. A = ININLIN numop numvar numdat numcof numder ;   
  726. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  727. A . 'VAR' . 1 . 'DISC'   = discv  ;
  728. A . 'VAR' . 2 .  'NOMDDL' = 'MOTS' 'UY' ;
  729. A . 'VAR' . 2 .  'DISC'   = discv  ;
  730. * cx
  731. A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  732. A . 'DAT' . 1 . 'DISC' =  discc ;
  733. A . 'DAT' . 1 . 'VALEUR' = cx ;
  734. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  735. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  736. * cy 
  737. A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  738. A . 'DAT' . 2 . 'DISC' =  discc ;
  739. A . 'DAT' . 2 . 'VALEUR' = cy ;
  740. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  741. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  742. * 
  743. *   Matrix B for the operator NLIN
  744. *
  745. numvar = 2 ;
  746. B = ININLIN numop numvar 0 0 numder ;
  747. *
  748. * For B, no coefficients are involved
  749. * but 1
  750. *
  751. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  752. B . 'VAR' . 1 . 'DISC'   = discv ;
  753. B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  754. B . 'VAR' . 2 . 'DISC'   = discv ;
  755. *
  756. * Contribution
  757. *
  758. * cx \dep{ux}{x} Nvx
  759. A . 1 . 1 . 1 = 'LECT' 1 ;
  760. B . 1 . 1 . 0 = 'LECT'   ;
  761. * cy \dep{ux}{y}  Nvx
  762. A . 1 . 1 . 2 = 'LECT' 2 ;
  763. * cx \dep{uy}{x} Nvy
  764. A . 2 . 2 . 1 = 'LECT' 1 ;
  765. B . 2 . 2 . 0 = 'LECT'   ;
  766. * cy \dep{uy}{y}  Nd
  767. A . 2 . 2 . 2 = 'LECT' 2 ;
  768. *
  769. mat = 'NLIN' discg _mt A B  methgau ;
  770. 'RESPRO' mat ;
  771. 'FINPROC' ;
  772. *
  773. *ENDPROCEDUR matkv
  774. *BEGINPROCEDUR mattau
  775. *
  776. ************************************************************************
  777. *
  778. * NOM         : MATTAU
  779. *
  780. * DESCRIPTION : We compute
  781. *
  782. *               ( (-coef*div(tau))|x, Nv)
  783. *               ( (-coef*div(tau))|y, Nv)
  784. *
  785. *    where
  786. *
  787. *    Nv is the form function for the speed (u_x,u_y)^T
  788. *          =     the test function for the speed equations   
  789. *
  790. *    tau_{xy} = mu (dep{u_x}{y} '+' dep{u_y}{x} '-'
  791. *                   \delta_{x,y} div (u))
  792. *
  793. *    These terms are integrated by part. Only the volumic contribution
  794. *    is computed. 
  795. *
  796. *
  797. * LANGAGE     : GIBIANE-CAST3M
  798. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  799. *               mél : beccantini@cea.fr
  800. *
  801. **********************************************************************
  802. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  803. * en cas de modification de ce sous-programme afin de faciliter
  804. * la maintenance !
  805. ************************************************************************
  806. *
  807. 'DEBPROC' MATTAU ;
  808. *
  809. 'ARGUMENT'     _mt*'MAILLAGE' ;
  810. 'ARGUMENT'   discv*'MOT     ' ;
  811. 'ARGUMENT'  discmu*'MOT     ' ;
  812. 'ARGUMENT'  discco*'MOT     ' ;
  813. 'ARGUMENT'      mu*'CHPOINT' ;
  814. 'ARGUMENT'    coef*'CHPOINT' ;
  815. *
  816. *  _mt    = surface QUAF mesh
  817. *  discv  = discretization type of the speed
  818. *  discmu = discretization type of mu
  819. *  discco = discretization type of coef
  820. *  mu     = CHPOINT ('SCAL')
  821. *  coef   = CHPOINT ('SCAL')
  822. *
  823. *    N.B. Names of the primal variable = names of the dual variable =
  824. *         'UX', 'UY'
  825. *
  826. idim   = 'VALEUR' 'DIME' ;
  827. *
  828. *  Matrix A for the operator NLIN
  829. *
  830. numop  = 4 ;
  831. numder = idim ;
  832. numvar = 2 ;
  833. * 
  834. numdat = 3 ;
  835. * The data: mu, 2./3., 4./3.
  836. numcof = 3 ;
  837. *
  838. * One function: f(mu)= mu
  839. *
  840. discg = 'LINE' ;
  841. * Geometrical discretisation
  842. methgau = 'GAU7' ;
  843. * Gauss points involved
  844. *
  845. A = ININLIN numop numvar numdat numcof numder ;   
  846. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  847. A . 'VAR' . 1 . 'DISC'   = discv ;
  848. A . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  849. A . 'VAR' . 2 . 'DISC'   = discv ;
  850. *
  851. * mu
  852. *
  853. mmu = 'EXTRAIRE' mu 'COMP' ;
  854. A . 'DAT' . 1 . 'NOMDDL' = mmu ;
  855. A . 'DAT' . 1 . 'DISC' =  discmu ;
  856. A . 'DAT' . 1 . 'VALEUR' = mu ;
  857. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  858. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  859. *
  860. * -2/3
  861. *
  862. A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  863. A . 'DAT' . 2 . 'DISC' =  'CSTE' ;
  864. A . 'DAT' . 2 . 'VALEUR' = (-2./3.) ;
  865. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  866. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  867. *
  868. * 4/3
  869. *
  870. A . 'DAT' . 3 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  871. A . 'DAT' . 3 . 'DISC' =  'CSTE' ;
  872. A . 'DAT' . 3 . 'VALEUR' = (4./3.) ;
  873. A . 'COF' . 3 . 'COMPOR' = 'IDEN' ;
  874. A . 'COF' . 3 . 'LDAT' = 'LECT' 3 ;
  875. *
  876. *   Matrix B for the operator NLIN
  877. *
  878. numdat = 1 ;
  879. numcof = 3 ;
  880. B = ININLIN numop numvar numdat numcof numder ;
  881. *
  882. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  883. B . 'VAR' . 1 . 'DISC'   = discv ;
  884. B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  885. B . 'VAR' . 2 . 'DISC'   = discv ;
  886. *
  887. B . 'DAT' . 1 . 'NOMDDL' = ('EXTRAIRE' coef 'COMP') ;
  888. B . 'DAT' . 1 . 'DISC' =  discco ;
  889. B . 'DAT' . 1 . 'VALEUR' = coef ;
  890. * coef, \dep{coef}{x}, \dep{coef}{y}
  891. B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  892. B . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  893. B . 'COF' . 2 . 'COMPOR' = 'D/DX1' ;
  894. B . 'COF' . 2 . 'LDAT' = 'LECT' 1 ;
  895. B . 'COF' . 3 . 'COMPOR' = 'D/DX2' ;
  896. B . 'COF' . 3 . 'LDAT' = 'LECT' 1 ;
  897. *
  898. *   Dual variable 'UX'
  899. *
  900. *   \int
  901. *   (2 mu \dep{u_x}{x} coef \dep{N_j}{x} '+'
  902. *    2 mu \dep{u_x}{x} N_j \dep{coef}{x} '+'
  903. *    -2/3 mu div(u)  coef \dep{N_j}{x} '+'
  904. *    -2/3 mu div(u)  N_j \dep{coef}{x} '+' 
  905. *     mu \dep{u_x}{y} coef \dep{N_j}{y} '+'
  906. *     mu \dep{u_x}{y} N_j \dep{coef}{y} '+'
  907. *     mu \dep{u_y}{x} coef \dep{N_j}{y} '+'
  908. *     mu \dep{u_y}{x} N_j \dep{coef}{y}) =
  909. *
  910. *   \int
  911. *   ((4/3 mu \dep{u_x}{x}) '-' (2/3 mu \dep{u_y}{y})) *
  912. *        (coef \dep{N_j}{x} '+' N_j \dep{coef}{x})
  913. *   '+'
  914. *   (mu \dep{u_x}{y}) '+' (mu \dep{u_y}{x}) *
  915. *        (coef \dep{N_j}{y} '+' N_j \dep{coef}{y})
  916. *
  917. A . 1 . 1 . 1 = 'LECT' 1 3  ;
  918. A . 1 . 2 . 2 = 'LECT' 1 2  ;
  919. B . 1 . 1 . 1 = 'LECT' 1 ;
  920. B . 1 . 1 . 0 = 'LECT' 2 ;
  921. *
  922. A . 2 . 1 . 2 = 'LECT' 1 ;
  923. A . 2 . 2 . 1 = 'LECT' 1 ;
  924. B . 2 . 1 . 2 = 'LECT' 1 ;
  925. B . 2 . 1 . 0 = 'LECT' 3 ;
  926. *
  927. *   Dual variables 'UY'   
  928. *
  929. *   \int
  930. *   ((4/3 mu \dep{u_y}{y}) '-' (2/3 mu \dep{u_x}{x})) *
  931. *        (coef \dep{N_j}{y} '+' N_j \dep{coef}{y})
  932. *   '+'
  933. *   (mu \dep{u_y}{x}) '+' (mu \dep{u_x}{y}) *
  934. *        (coef \dep{N_j}{x} '+' N_j \dep{coef}{x})
  935. *    
  936. *
  937. A . 3 . 2 . 2 = 'LECT' 1 3 ;
  938. A . 3 . 1 . 1 = 'LECT' 1 2 ;
  939. B . 3 . 2 . 2 = 'LECT' 1 ;
  940. B . 3 . 2 . 0 = 'LECT' 3 ;
  941. *
  942. A . 4 . 2 . 1 = 'LECT' 1 ;
  943. A . 4 . 1 . 2 = 'LECT' 1 ;
  944. B . 4 . 2 . 1 = 'LECT' 1 ;
  945. B . 4 . 2 . 0 = 'LECT' 2 ;
  946. *
  947. mat = 'NLIN' discg _mt A B methgau ;
  948. 'RESPRO' mat ;
  949. 'FINPROC' ;
  950. *
  951. *ENDPROCEDUR mattau
  952. *BEGINPROCEDUR matpre
  953. *
  954. ************************************************************************
  955. *
  956. * NOM         : MATPRE
  957. *
  958. * DESCRIPTION : We compute the integral of volume of
  959. *
  960. *              (coef \dep{p}{x} , Nv)
  961. *              (coef \dep{p}{y} , Nv) 
  962. *
  963. *              where Nv = test function for u_x = test function for u_y.
  964. *
  965. *              \int_{\Omega} coef  \dep{p}{x} Nv dV =
  966. *              \int_{\delta \Omega} coef p Nv (i \cdot n) dS '+'
  967. *              -\int_{\Omega} coef p \dep{nv}{x} dV '+'
  968. *              -\int_{\Omega} Nv p \dep{coef}{x} dV
  969. *
  970. *              Here we only compute the 2nd and 3rd contribution (volume
  971. *              integrals). The surface integral must be computed elsewhere.
  972. *
  973. * LANGAGE     : GIBIANE-CAST3M
  974. * AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  975. *               mél : beccantini@cea.fr
  976. *
  977. **********************************************************************
  978. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  979. * en cas de modification de ce sous-programme afin de faciliter
  980. * la maintenance !
  981. ************************************************************************
  982. *
  983. *
  984. 'DEBPROC' MATPRE ;
  985. *
  986. 'ARGUMENT'   _mt*'MAILLAGE'   ;
  987. 'ARGUMENT'  nompre*'MOT     ' ;
  988. 'ARGUMENT'  discp*'MOT     '  ;
  989. 'ARGUMENT'  discv*'MOT     '  ;
  990. 'ARGUMENT'  discc*'MOT     '  ;
  991. 'ARGUMENT'  coef*'CHPOINT '   ;
  992. *
  993. * _mt    = surface QUAF mesh
  994. * nompre = name of the pressure (usuallt 'LX')
  995. * discp  = discretization type of p (usually LINE)
  996. * discv  = discretization type of v (usually QUAF)
  997. * discc = discretization type of coef
  998. *  coef   = CHPOINT ('SCAL')
  999. *
  1000. idim   = 'VALEUR' 'DIME' ;
  1001. *
  1002. *  Matrix A for the operator NLIN
  1003. *
  1004. numop  = 4 ;
  1005. numder = idim ;
  1006. numvar = 1 ;
  1007. * The variable involved in A is the primal p
  1008. numdat = 1 ;
  1009. * The data: coef
  1010. numcof = 3 ;
  1011. * Three functions: f(coef)= coef
  1012. *                  f(coef)=\dep{coef}{x}
  1013. *                  f(coef)=\dep{coef}{y}
  1014. *
  1015. discg = 'LINE' ;
  1016. * Geometrical discretisation
  1017. methgau = 'GAU7' ;
  1018. * Gauss points involved
  1019. *
  1020. A = ININLIN numop numvar numdat numcof numder ;   
  1021. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nompre ;
  1022. A . 'VAR' . 1 . 'DISC'   = discp ;
  1023. *
  1024. * coef
  1025. *
  1026. mcoef = 'EXTRAIRE' coef 'COMP' ;
  1027. A . 'DAT' . 1 . 'NOMDDL' = mcoef ;
  1028. A . 'DAT' . 1 . 'DISC' =  discc ;
  1029. A . 'DAT' . 1 . 'VALEUR' = coef ;
  1030. * coef, \dep{coef}{x}, \dep{coef}{y}
  1031. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  1032. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  1033. A . 'COF' . 2 . 'COMPOR' = 'D/DX1' ;
  1034. A . 'COF' . 2 . 'LDAT' = 'LECT' 1 ;
  1035. A . 'COF' . 3 . 'COMPOR' = 'D/DX2' ;
  1036. A . 'COF' . 3 . 'LDAT' = 'LECT' 1 ;
  1037. *
  1038. *   Matrix B for the operator NLIN
  1039. *
  1040. numvar = 2 ;
  1041. *
  1042. * ux, uy
  1043. *
  1044. numdat = 1 ;
  1045. numcof = 1 ;
  1046. B = ININLIN numop numvar numdat numcof numder ;
  1047. *
  1048. * For B the coef no coefficients are involved
  1049. * but -1
  1050. *
  1051. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  1052. B . 'VAR' . 1 . 'DISC'   = discv ;
  1053. B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  1054. B . 'VAR' . 2 . 'DISC'   = discv ;
  1055. * -1
  1056. B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  1057. B . 'DAT' . 1 . 'DISC' =  'CSTE' ;
  1058. B . 'DAT' . 1 . 'VALEUR' = -1.0 ;
  1059. B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  1060. B . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  1061. *
  1062. *   -\int_{\Omega} coef p \dep{Nv}{x} dV
  1063. *
  1064. A . 1 . 1 . 0 = 'LECT' 1 ;
  1065. B . 1 . 1 . 1 = 'LECT' 1 ;
  1066. *
  1067. *   -\int_{\Omega} Nv p \dep{coef}{x} dV
  1068. *
  1069. A . 2 . 1 . 0 = 'LECT' 2 ;
  1070. B . 2 . 1 . 0 = 'LECT' 1 ;
  1071. *
  1072. *   -\int_{\Omega} coef p \dep{Nv}{y} dV (second dual variable)
  1073. *
  1074. A . 3 . 1 . 0 = 'LECT' 1 ;
  1075. B . 3 . 2 . 2 = 'LECT' 1 ;
  1076. *
  1077. *   -\int_{\Omega} Nv p \dep{coef}{y} dV (second dual variable)
  1078. *
  1079. A . 4 . 1 . 0 = 'LECT' 3 ;
  1080. B . 4 . 2 . 0 = 'LECT' 1 ;
  1081. *
  1082. mat = 'NLIN' discg _mt A B methgau ;
  1083. 'RESPRO' mat ;
  1084. 'FINPROC' ;
  1085. *
  1086. *ENDPROCEDUR matpre
  1087. *BEGINPROCEDUR resoup
  1088. **********************************************************
  1089. **** Resolution of a linear system ***********************
  1090. **********************************************************
  1091. *
  1092. 'DEBPROC' RESOUP ;
  1093.    'ARGUMENT' LOGTPS*'LOGIQUE ' ;
  1094.    'ARGUMENT' mat*'MATRIK  ' ;
  1095.    'ARGUMENT' matpre/'MATRIK  ' ;
  1096.    'ARGUMENT' smb*'CHPOINT ' ;
  1097.    'ARGUMENT' ccl*'CHPOINT ' ;
  1098.    'ARGUMENT' res*'FLOTTANT' ;
  1099.    'ARGUMENT' nit*'ENTIER'   ;
  1100.    'ARGUMENT' pre/'ENTIER'   ;
  1101.    'ARGUMENT' gggtre/'TABLE   ' ;
  1102.    'ARGUMENT' gggtcv/'LOGIQUE ' ;
  1103.    'ARGUMENT' typslv/'MOT     ' ;
  1104. *
  1105. *  Det. tt s.t.
  1106. *
  1107. *  mat . tt = smb
  1108. *  tt(boundary) = ccl
  1109. *
  1110. *  
  1111. *  mat = matrix to "inverse"
  1112. *  matpre = matrix for which preconditioner has been computed (optional)
  1113. *  smb = right hand side
  1114. *  ccl = imposed boundary conditions
  1115. *  res = rvk . 'RESID'
  1116. *  nit = max number of linear iteration (for an iterative solver)
  1117. *  pre = type of preconditioning matrix (for an iterative solver)
  1118. *  gggtre = table of resolution to store 'XINIT'
  1119. *  gggtcv = graphic of the history for the iterative solver (optional FAUX)
  1120. *  typslv = 'DIRECT'/'ITER' (optional, ITER)
  1121. *
  1122.    'SI' ('NON' ('EXISTE' typslv)) ;
  1123.       tslv = VRAI ;
  1124.    'SINON' ;      
  1125.       'SI' ('EGA' typslv 'DIRECT') ;
  1126.          tslv = FAUX ;
  1127.       'SINON' ;
  1128.          'SI' ('EGA' typslv 'ITER') ;
  1129.             tslv = VRAI ;
  1130.          'SINON' ;
  1131.             'ERREUR' 'Solveur DIRECT ou ITERatif ?' ;
  1132.          'FINSI' ;
  1133.       'FINSI' ;
  1134.    'FINSI' ;
  1135.     gggniv = 1 ;
  1136.    'SI' ('NON' ('EXISTE' pre)) ;
  1137.       pre = 5 ;
  1138.    'FINSI' ;
  1139.    'SI' ('OU' ('NON' ('EXISTE' gggtcv)) ('NON' tslv)) ;
  1140.       gggtcv = FAUX ;
  1141.    'FINSI' ;      
  1142.    rv = 'EQEX' ;
  1143.    rvk = rv . 'METHINV' ;
  1144.    'SI' ('EXISTE' matpre) ;
  1145.       rvk . 'MATASS' = matpre ;
  1146.       rvk . 'MAPREC' = matpre ;
  1147.    'SINON' ;      
  1148.       rvk . 'MATASS' = mat ;
  1149.       rvk . 'MAPREC' = mat ;
  1150.    'FINSI' ;      
  1151.    'SI' tslv ;
  1152.       rvk . 'TYPINV' = 3 ;
  1153.    'SINON' ;      
  1154.       rvk . 'TYPINV' = 1 ;
  1155.    'FINSI' ;      
  1156.    rvk . 'SCALING' = 1 ;
  1157. *  Scaling factor   
  1158.    rvk . 'OUBMAT' = 1 ;
  1159. *  oubmat = 1 -> We eliminate the elementary matrix
  1160. *  'ILUTPPIV' -> 1 We always search for a new pivot
  1161. *             -> 0 We do not search for a new pivot
  1162.    rvk . 'ILUTPPIV' = 0.1D0 ;
  1163. *   rvk . 'ILUTPPIV' = 0.01D0 ;
  1164. *   rvk . 'ILUTPPIV' = 0.00001D0 ;
  1165. *   rvk . 'ILUTPPIV' = 0.D0 ;
  1166.    rvk . 'IMPINV' = gggniv ;
  1167.    rvk . 'TYRENU' = 'SLOA' ;
  1168.    rvk . 'PCMLAG' = 'APR2' ;
  1169.    rvk . 'NITMAX' = nit ;
  1170.    rvk . 'RESID' = res ;
  1171.    rvk . 'PRECOND' = pre ;
  1172.    rvk . 'BCGSBTOL' = 1.D-200 ;
  1173.    rvk . 'ILUTLFIL' = 1.5 ;
  1174. *
  1175.    'SI' ('EXISTE' gggtre) ;
  1176.       'SI' ('EXISTE' gggtre 'XINIT') ;
  1177.          rvk . 'XINIT' = gggtre . 'XINIT' ;
  1178.       'FINSI' ;
  1179.    'FINSI' ;      
  1180.    'TEMPS' 'ZERO' ;
  1181.    tt = 'KRES' mat 'TYPI' rvk 'SMBR' smb 'CLIM' ccl ;
  1182.    TABTPS = TEMP 'NOEC';
  1183.    tcpu = TABTPS.'TEMPS_CPU'.'INITIAL';
  1184.    tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
  1185.    'SI' LOGTPS ;
  1186.       'MESSAGE' ('CHAINE' 'tcpu (s) = ' tcpus) ;
  1187.    'FINSI' ;
  1188.    'SI' ('ET' gggtcv ('EXISTE' rvk 'CONVINV')) ;
  1189.       lcvg = rvk . 'CONVINV' ;
  1190.       nit = 'DIME' lcvg ;
  1191.       'SI' ('>' nit 1) ;
  1192.          lit = 'PROG' 1.D0 PAS 1.D0 NPAS ('-' nit 1) ;
  1193.          evit = 'EVOL' 'MANU' 'iter' lit 'log residu'
  1194.           ('/' ('LOG' lcvg) ('LOG' 10.D0)) ;
  1195.          'DESSIN' evit 'TITR' titglob ;
  1196.       'FINSI' ;         
  1197.    'FINSI' ;
  1198.    'SI' ('EXISTE' gggtre) ;
  1199.       gggtre . 'XINIT' = tt ;
  1200.    'FINSI' ;      
  1201.    'RESPRO' tt ;
  1202. *
  1203. * End of procedure file RESOU
  1204. *
  1205. 'FINPROC' ;
  1206. *
  1207. **********************************************************
  1208. **** End of the Resolution of a linear system ************
  1209. **********************************************************
  1210. *ENDPROCEDUR resoup
  1211. *
  1212.  'DEBPROC' INTINVT ;
  1213.  'ARGUMENT'   _mt*'MAILLAGE' ;
  1214.  'ARGUMENT' disct*'MOT     ' ;
  1215.  'ARGUMENT'    tn*'CHPOINT ' ;
  1216. *
  1217. idim   = 'VALEUR' 'DIME' ;
  1218. *
  1219. *   Matrix A for the operator NLIN
  1220. *
  1221.  numop  = 1 ;
  1222.  numder = idim ;
  1223.  numvar = 1 ;
  1224. * The variable involved in A is the primal variable
  1225.  numdat = 1 ;
  1226.  numcof = 1 ;
  1227. *
  1228.  discg = 'LINE' ;
  1229. * Geometrical discretisation
  1230.  methgau = 'GAU7' ;
  1231. * Gauss points involved
  1232. *
  1233. A = ININLIN numop numvar numdat numcof numder ;
  1234. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  1235. A . 'VAR' . 1 . 'DISC'   = 'QUAF' ;
  1236. *
  1237. A . 'DAT' . 1 . 'NOMDDL' = ('EXTRAIRE' TN 'COMP') ;
  1238. A . 'DAT' . 1 . 'DISC'   =  disct ;
  1239. A . 'DAT' . 1 . 'VALEUR' =  tn ;
  1240. A . 'COF' . 1 . 'COMPOR' =  'IDEN' ;
  1241. A . 'COF' . 1 . 'LDAT'   =  'LECT' 1 ;
  1242. *
  1243. *   Matrix B for the operator NLIN
  1244. *
  1245. B = ININLIN numop numvar 0 0 numder ;
  1246. * For B the coef no coefficients are involved
  1247. * but 1
  1248. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  1249. B . 'VAR' . 1 . 'DISC'   = 'CSTE' ;
  1250. *
  1251. A . 1 . 1 . 0 = 'LECT' -1 ;
  1252. B . 1 . 1 . 0 = 'LECT'    ;
  1253. *
  1254. mat = 'NLIN' discg _mt A B methgau ;
  1255. *
  1256. chpbid = 'MANUEL' 'CHPO' _mt 1 'SCAL' 1.0D0 ;
  1257. 'RESPRO' ('SOMT' (mat '*' chpbid)) ;
  1258. 'FINPROC' ;
  1259. * 
  1260. *************************************************************************
  1261. *************************************************************************
  1262. ******* END OF THE PERSONAL L PROCEDURES ********************************
  1263. *************************************************************************
  1264. *************************************************************************
  1265. *
  1266. ************************************************************************
  1267. ************** Begin of file *******************************************
  1268. ************************************************************************
  1269. **
  1270. *************************************************
  1271. ****** LIST OF SYMBOLS **************************
  1272. *************************************************
  1273. *
  1274. * \int  = integral
  1275. * \div  = divergence
  1276. * \der  = derivative (\der{a}{b} = \frac{d a}{d b} 
  1277. * \grad = gradient
  1278. * \cdot = scalar product
  1279. * ^     = power
  1280. *
  1281. * (a,b) = \int_{\Omega} a * b dV
  1282. *
  1283. *  Np   = test function for the divergence equations
  1284. *       = form function for the pressure
  1285. *  Nt   = test function for the temperature
  1286. *       = form function for the temperature equation 
  1287. *  Nv   = test function for the speed equations (x,y)
  1288. *       = form function for the speed
  1289. *
  1290. *    
  1291. ************************************************************************
  1292. ************** Begin of file *******************************************
  1293. ************************************************************************
  1294.  NX = 6 ;
  1295. *
  1296. * IC incoming from other computations
  1297. *
  1298. * REPSAU = '/test4/becc/Dtmp/Dcaviteimplpv/' ;
  1299.  REPSAU = './Dsauv/' ;
  1300. *
  1301. ****************************************
  1302. **** Parameters for the computation ****
  1303. ****************************************
  1304. *
  1305. * BDF1 or BDF2
  1306. *
  1307.  LOGBDF1 = FAUX ;
  1308. * Tensorial/scalar/diagonal upwinding
  1309.  LOGUP = FAUX ;
  1310.  BETAUP = 0.0 ;
  1311. * BETAUP = 1 -> SU
  1312. * BETAUP = 0 -> No SU
  1313.  LOGSCAL = FAUX ;
  1314.  'MESSAGE' ;
  1315.  'MESSAGE' 'Scalar numerical diffusion' ;
  1316.  'MESSAGE' ;
  1317.  LOGDIA = FAUX ;
  1318. * NB 'SI' LOGDIA, LOGSCAL = FAUX
  1319.  'SI' LOGDIA ;
  1320.     LOGSCAL = FAUX ;
  1321.  'FINSI' ;
  1322. * 
  1323. * Boussinesq or Low Mach?
  1324. *
  1325.  LBOUS = FAUX ;
  1326. * 
  1327.  DISCT = 'QUAF' ;
  1328.  DISCMU = DISCT ;
  1329.  DISCLA = DISCT ;
  1330. * 
  1331. * NB If DISCLA, DISCMU, DISCT are not 'QUAF'
  1332. *    SUPG does not work
  1333. *
  1334.  TFINAL = 0.25 ;
  1335. * NINT = 10 ;
  1336.  EPSERR = 1.0D-20 ;
  1337.  GRAPH = FAUX ;
  1338. * 
  1339.  MTRA = 'X' ;
  1340.  NFREQ = 2 ;
  1341. * METINV = MOT 'ITER' ;
  1342.  METINV = MOT 'DIRECT' ;
  1343. * PREC   = 5 ;
  1344.  PREC = 5 ;
  1345.  NFREQP = 5 ;
  1346.  ERROPRE = 1.0D-5 ;
  1347. * We compute the preconditioner each NFREQP-th iterations
  1348. * If the residial is lower than ERROPRE, we stop computing
  1349. * the preconditioner
  1350. * Number of linear iterations (for NX = 160)
  1351.  NITERT = NX  ;
  1352. * NITERU = NX '/' 1 ;
  1353. * NITERUT = NX '/' 4 ;
  1354. * NITERP = NX ;
  1355. * NITERT = NX ;
  1356. * NITERU = NX ;
  1357. * NITERUT = NX  '*' 2;
  1358. * NITERP = NX '*' 4 ;
  1359. * CFLREF = 1.0D6 ;
  1360. * DT computed with respect to the reference speed
  1361. *
  1362. *
  1363. * Fixed options:
  1364. * - LINEar/bilinear continuos pressure 
  1365. * - QUAFdratic speed
  1366. * - Name of speed (and name of speed duals variables) :
  1367. *   'UX', 'UY'
  1368. * - Name of temperature: T
  1369. * - Name of pressure: NOMPRE (Lagrange multiplicator)
  1370. *
  1371. *
  1372.  'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'TRAC' MTRA 'ECHO' 1 ;
  1373. *
  1374. * Geometry dimensions
  1375. *
  1376.  D1   = 1.0 ;
  1377.  L1   = 12.0 '*' d1 ;
  1378.  L2   = 9.0 '*' d1 ;
  1379.  L3   = 7.0 '*' d1 ;
  1380.  DX     = d1 '/' NX ;
  1381.  DXMAX  = 4 '*' DX ;
  1382. * 
  1383. * d1 used into Reynolds number 
  1384. *
  1385. * Physical properties, reference/initial states, BC
  1386. *
  1387.  R = 288. ;
  1388.  gamma = 1.4 ;
  1389.  g = 9.0 ;
  1390. *
  1391.  PR = 0.7 ;
  1392.  RE = 100. ;
  1393.  TH = 267.8 ;
  1394.  TC = 178.6 ;
  1395.  mt = 1.0 ;
  1396.  mH = 0.8 '*' mt;
  1397.  mC = 0.2 '*' mt ;
  1398. *
  1399.  cp = R * gamma '/' (gamma '-' 1.0) ; 
  1400. *
  1401.  T0 = 223.2 ;
  1402.  P0 = 0.6429D5 ;
  1403. *
  1404.  rho0 = P0 '/' (R '*' T0) ;
  1405.  rhoC0 = P0 '/' (R '*' TC) ;
  1406.  rhoH0 = P0 '/' (R '*' TH) ;
  1407. *
  1408.  uC0 = mC '/' rhoC0 ;
  1409.  uH0 = mH '/' rhoH0 ;
  1410. *
  1411. * mu (dynamic viscosity) and lambda (thermal diffusivity) at the
  1412. * reference state 
  1413. *
  1414.  mu0 = (mt '*' d1) '/' Re ;
  1415.  lambda0 = mu0 '*' (cp '/' Pr) ;
  1416. * 
  1417. ***********************
  1418. **** Mesh *************
  1419. ***********************
  1420. *
  1421. *
  1422. * P31             P32            P33              P34
  1423. *
  1424. *
  1425. *       TUBE2            TUBE3          TUBE4
  1426. *
  1427. *
  1428. * P21             P22            P23              P24
  1429. *
  1430. *
  1431. *
  1432. *                          T
  1433. *                          U
  1434. *                          B
  1435. *                          1
  1436. *
  1437. *
  1438. *
  1439. *
  1440. *
  1441. *
  1442. *
  1443. *
  1444. *                 P12            P13
  1445. *
  1446. *
  1447. * 
  1448. * 
  1449.  P11 = 0.0 0.0 ;
  1450.  P12 = L1  0.0 ;
  1451.  P13 = (L1 '+' d1) 0.0 ;
  1452. *
  1453.  P21 = 0.0 L3 ;
  1454.  P22 = L1  L3 ;
  1455.  P23 = (L1 '+' d1) L3 ;
  1456.  P24 = (L1 '+' d1 '+' L2) L3 ;
  1457. * 
  1458.  P31 = 0.0 (L3 '+' d1) ;
  1459.  P32 = L1  (L3 '+' d1) ;
  1460.  P33 = (L1 '+' d1) (L3 '+' d1) ;
  1461.  P34 = (L1 '+' d1 '+' L2) (L3 '+' d1) ;
  1462. *
  1463. * Lignes
  1464. *
  1465.  P12P13 = P12 'DROIT' NX P13 ;
  1466.  P22P23 = P22 'DROIT' NX P23 ;
  1467.  P32P33 = P32 'DROIT' NX P33 ;
  1468. *
  1469.  P12P22 = P12 'DROIT' P22 'DINI' DXMAX 'DFIN' DX ;
  1470.  P13P23 = P12P22 'PLUS' (d1 0.0) ;
  1471.  'ELIMINATION' (P13P23 'ET' P13 'ET' P23) (DX '/' 100.) ;
  1472.  TUB1 = 'DALL' P12P13 P13P23 ('INVERSE' P22P23) ('INVERSE' P12P22)
  1473.     'PLAN' ;
  1474. *
  1475.  P24P34 = P24 'DROIT' NX P34 ;
  1476.  P23P33 = P23 'DROIT' NX P33 ;
  1477.  P22P32 = P22 'DROIT' NX P32 ;
  1478.  P21P31 = P21 'DROIT' NX P31 ;
  1479. *
  1480.  P24P23 = P24 'DROIT' P23 'DINI' DXMAX 'DFIN' DX ;
  1481.  P34P33 = P24P23 'PLUS' (0.0 d1) ;
  1482.  'ELIMINATION' (P34P33 'ET' P34 'ET' P33) (DX '/' 100.) ;
  1483. *
  1484.  P21P22 = P21 'DROIT' P22 'DINI' DXMAX 'DFIN' DX ;
  1485.  P31P32 = P21P22 'PLUS' (0.0 d1) ;
  1486.  'ELIMINATION' (P31P32 'ET' P31 'ET' P32) (DX '/' 100.) ;
  1487. * 
  1488.  TUB2 = 'DALL' P21P22 P22P32 ('INVERSE' P31P32) ('INVERSE' P21P31)
  1489.     'PLAN' ;
  1490.  TUB3 = 'DALL' P22P23 P23P33 ('INVERSE' P32P33) ('INVERSE' P22P32)
  1491.     'PLAN' ;
  1492.  TUB4 = 'DALL' P24P34 P34P33 ('INVERSE' P23P33) ('INVERSE' P24P23)
  1493.     'PLAN' ;
  1494. *   
  1495.  DOM1 = TUB1 'ET' TUB2 'ET' TUB3 'ET' TUB4 ;
  1496. * 
  1497.  'SI' GRAPH ;
  1498.     'TRACER' DOM1 'TITRE' 'Mesh' ;
  1499.  'FINSI' ;
  1500. *
  1501. * DOM1 and company are made of 'LINE' elts
  1502. *
  1503.  QDOM1 = 'CHANGER' DOM1 'QUAF' ;
  1504.  QP12P13 = 'CHANGER' P12P13 'QUAF' ;
  1505.  'ELIMINATION' QDOM1 (DX '/' 100.) QP12P13 ;
  1506.  QP22P23 = 'CHANGER' P22P23 'QUAF' ;
  1507.  'ELIMINATION' QDOM1 (DX '/' 100.) QP22P23 ;
  1508.  QP32P33 = 'CHANGER' P32P33 'QUAF' ;
  1509.  'ELIMINATION' QDOM1 (DX '/' 100.) QP32P33 ;
  1510. *
  1511.  QP12P22 = 'CHANGER' P12P22 'QUAF' ;
  1512.  'ELIMINATION' QDOM1 (DX '/' 100.) QP12P22 ;
  1513.  QP13P23 = 'CHANGER' P13P23 'QUAF' ;
  1514.  'ELIMINATION' QDOM1 (DX '/' 100.) QP13P23 ;
  1515. *
  1516.  QP24P34 = 'CHANGER' P24P34 'QUAF' ;
  1517.  'ELIMINATION' QDOM1 (DX '/' 100.) QP24P34 ;
  1518.  QP23P33 = 'CHANGER' P23P33 'QUAF' ;
  1519.  'ELIMINATION' QDOM1 (DX '/' 100.) QP23P33 ;
  1520.  QP22P32 = 'CHANGER' P22P32 'QUAF' ;
  1521.  'ELIMINATION' QDOM1 (DX '/' 100.) QP22P32 ;
  1522.  QP21P31 = 'CHANGER' P21P31 'QUAF' ;
  1523.  'ELIMINATION' QDOM1 (DX '/' 100.) QP21P31 ;
  1524. *
  1525.  QP24P23 = 'CHANGER' P24P23 'QUAF' ;
  1526.  'ELIMINATION' QDOM1 (DX '/' 100.) QP24P23 ;
  1527.  QP34P33 = 'CHANGER' P34P33 'QUAF' ;
  1528.  'ELIMINATION' QDOM1 (DX '/' 100.) QP34P33 ;
  1529. *
  1530.  QP21P22 = 'CHANGER' P21P22 'QUAF' ;
  1531.  'ELIMINATION' QDOM1 (DX '/' 100.) QP21P22 ;
  1532.  QP31P32 = 'CHANGER' P31P32 'QUAF' ;
  1533.  'ELIMINATION' QDOM1 (DX '/' 100.) QP31P32 ;
  1534. * 
  1535. * QDOM1 and friends are now QUAF
  1536. *
  1537.  VOL = 'MESURE' QDOM1 ;
  1538. *
  1539. * Mesh dimension for upwinding
  1540. * HUP on DISCT (otherwise, it does not work)
  1541. *
  1542.  HUP = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' ('MAXIMUM' ('PROG' DX DXMAX)) ;
  1543. *********************************************************
  1544. *** Dynamic viscosity and thermal diffusivity and the ***
  1545. *** unitary CHPOINT                                   ***
  1546. *********************************************************
  1547. *
  1548. * CHMU and CHLAM and CHUN
  1549. *
  1550.  'SI' ('EGA' DISCMU 'LINE') ;
  1551.     CHMU = 'MANUEL' 'CHPO' DOM1 1 'SCAL' mu0 ;
  1552.  'FINS' ;
  1553.  'SI' ('EGA' DISCMU 'QUAF') ;
  1554.     CHMU = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' mu0 ;
  1555.  'FINS' ;
  1556. *
  1557.  'SI' ('EGA' DISCLA 'LINE') ;
  1558.     CHLAM = 'MANUEL' 'CHPO' DOM1 1 'SCAL' lambda0 ;
  1559.  'FINS' ;
  1560.  'SI' ('EGA' DISCLA 'QUAF') ;
  1561.     CHLAM = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' lambda0 ;
  1562.  'FINS' ;
  1563. *
  1564.  'SI' ('EGA' DISCT  'LINE') ;
  1565.     CHUN  = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
  1566.  'FINS' ;
  1567.  'SI' ('EGA' DISCT  'QUAF') ;
  1568.     CHUN  = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' 1.0 ;
  1569.  'FINS' ;
  1570. **************
  1571. *** DT *******
  1572. **************
  1573. *
  1574. * DT = CFLREF '*' DX '/' (mt '/' rho0) ;
  1575. * NITER = 'ENTIER' (TFINAL '/' DT) ;
  1576.  NITER = 2 ;
  1577.  DT = TFINAL '/' NITER ;
  1578. *
  1579. *
  1580. ******************************************
  1581. *** Fully-implicit (or projection) ? *****
  1582. ******************************************
  1583. *
  1584.  'REPETER' BLMET 2 ;
  1585.  'SI' ('EGA' &BLMET 1) ;
  1586.     LOGFI = VRAI ;
  1587.  'SINON' ;
  1588.     LOGFI = FAUX ;
  1589.  'FINSI' ;
  1590.  'SI' (LOGFI) ;
  1591.     NOMPRE = 'LX' ;
  1592.     NINT = 8 ;
  1593.  'SINON' ;
  1594.     NOMPRE = 'P' ;
  1595.     NINT = 30 ;
  1596.  'FINSI' ;
  1597.  
  1598. ***************************************
  1599. ********** IC / BC ********************
  1600. ***************************************
  1601. *
  1602.  PTHER = P0 ;
  1603. *
  1604. * IC/BC for the dynamic pressure (alaways 'LINE')
  1605. *
  1606. * Projection: imposition of PNLIM '+' tau.n on MPLIM (exit)
  1607. *             imposition of \delta PNLIM = 0 on MPLIM (exit)
  1608. *             imposition of the speed on the rest
  1609. *             imposition of \grad \delta PNLIM . n = 0 on MPLIM
  1610. *
  1611. * Fully-implicit: the same (except thet there is not \delta PNLIM)
  1612. *
  1613. * We impose the pressure and a free boundary condition on a point
  1614. * belonging to the upper well.
  1615.  MPLIM = QP22P32 'POIN'  'PROC'
  1616.    ((L1 '+' (d1 '/' 2)) (L3 '+' d1)) ;
  1617.  MPLIM = ('MANU' 'POI1' MPLIM ) 'COULEUR' 'ROUG' ;
  1618.  PNLIM = 'MANUEL' 'CHPO' MPLIM 1 NOMPRE 0.0 ;
  1619. * PNLIM not taken into account in the fully implicit case 
  1620.  PN = 'MANUEL' 'CHPO' DOM1 1 NOMPRE 0.0 ;
  1621. *
  1622. * BC/IC for the speed (always on 'QUAF')
  1623. * 
  1624.  QZEROV  = ('DIFF' ('CHANGER' ('CONTOUR' QDOM1) 'POI1') 
  1625.     ('CHANGER' 'POI1' (QP21P31 'ET' QP12P13 'ET' QP24P34)))
  1626.     'COULEUR' 'VERT' ;
  1627.  QZEROV = QZEROV 'ET' P21 'ET' P31 'ET' P24 'ET' P34 'ET' P12 'ET' P13; 
  1628.  QZEROX = (QZEROV 'ET' ('CHANGER' 'POI1' QP12P13)) 'COULEUR'
  1629.     'VERT' ;
  1630.  QZEROY = (QZEROV 'ET' ('CHANGER' 'POI1' QP24P34)
  1631.   'ET' ('CHANGER' 'POI1' QP21P31)) 'COULEUR'
  1632.     'VERT' ;
  1633.  QZEROY = 'DIFF' QZEROY MPLIM ;
  1634. *
  1635.  'SI' GRAPH ;
  1636.     'TRACER' (QDOM1 'ET' (MPLIM 'COULEUR' ROUG))
  1637.        'TITRE' 'Point for PLIM in projection.' ;
  1638.     'TRACER' (QDOM1 'ET' (QZEROX 'COULEUR' ROUG))
  1639.        'TITRE' 'Points in which ux is 0' ;
  1640.     'TRACER' (QDOM1 'ET' (QZEROY 'COULEUR' ROUG))
  1641.        'TITRE' 'Points in which uy is 0' ;
  1642.  'FINSI' ;
  1643. *
  1644. * P12P13 is the cold leg
  1645. *
  1646. * JC220346 => removes duplicate nodes before MANU CHPO (& hope it's OK)
  1647. QZEROX = 'CHANGER' 'POI1' QZEROX ;
  1648. QZEROY = 'CHANGER' 'POI1' QZEROY ;
  1649.  
  1650.  UN1 = 'MANUEL' 'CHPO' QZEROX 1 'UX' 0.0 ;
  1651.  UN2 = 'MANUEL' 'CHPO' QZEROY 1 'UY' 0.0 ;
  1652. *
  1653. * UN1 and UN2 used to correct the BC on the speed during
  1654. * the internal iterations, as well as UNGEOMC and
  1655. * UNGEOMH
  1656. * 
  1657.  UNGEOMC = ('COORDONNEE' 1 QP12P13) '-' (L1 '+' (d1 '/' 2.)) ;
  1658.  UNGEOMC = UNGEOMC '/' d1 ;
  1659.  UNGEOMC = UNGEOMC * UNGEOMC ;
  1660.  UNGEOMC = UNGEOMC '-' 0.25 ;
  1661.  UNGEOMC = UNGEOMC '*' -6 ;
  1662.  UN3 = UNGEOMC * uC0 ;
  1663.  UN3 = 'NOMC' 'UY' UN3 ;
  1664. *
  1665. * P24P34 is the hot leg
  1666. *
  1667.  UNGEOMH = ('COORDONNEE' 2 QP24P34) '-' (L3 '+' (d1 '/' 2.)) ;
  1668.  UNGEOMH = UNGEOMH '/' d1 ;
  1669.  UNGEOMH = UNGEOMH * UNGEOMH ;
  1670.  UNGEOMH = UNGEOMH '-' 0.25 ;
  1671.  UNGEOMH = UNGEOMH '*' 6 ;
  1672.  UN4 = UNGEOMH * uH0 ;
  1673.  UN4 = 'NOMC' 'UX' UN4 ;
  1674. *
  1675. * P21P31 is the outlet
  1676. *
  1677.  UNGEOMO = ('COORDONNEE' 2 QP21P31) '-' (L3 '+' (d1 '/' 2.)) ;
  1678.  UNGEOMO = UNGEOMO '/' d1 ;
  1679.  UNGEOMO = UNGEOMO * UNGEOMO ;
  1680.  UNGEOMO = UNGEOMO '-' 0.25 ;
  1681.  UNGEOMO = UNGEOMO '*' 6 ;
  1682. * 
  1683. * UN5 will be updated once the pressure derivative is 
  1684. * known. For the moment UN5 = 0
  1685. *
  1686.  UN5 = UNGEOMO '*' 0.0 ;
  1687.  UN5 = 'NOMC' 'UX' UN5 ;
  1688.  UNLIM = UN1 '+' UN2 '+' UN3 '+' UN4 '+' UN5 ;
  1689. * 
  1690.  UN = ('MANUEL' 'CHPO' QDOM1 2 'UX' 0.0 'UY' 0.0) '+' UNLIM ;
  1691.  DUNLIM = 0.0 '*' UNLIM ;
  1692. *
  1693. * BC/IC for the temperature ('QUAF'/'LINE')
  1694. *
  1695.  'SI' ('EGA' DISCT 'LINE') ;
  1696.     TNLIM = 'MANUEL' 'CHPO' P12P13 1 'T' TC ;
  1697.     TNLIM = TNLIM '+' ('MANUEL' 'CHPO' P24P34 1 'T' TH) ;
  1698.     TN = ('MANUEL' 'CHPO' DOM1 1 'T' T0) ;
  1699.     TN = TN '+' TNLIM '-' ('REDU' TN (P12P13 'ET' P24P34)) ;
  1700.     DOMT = DOM1 ;
  1701.  'FINS' ;
  1702. *
  1703.  'SI' ('EGA' DISCT 'QUAF') ;
  1704.     TNLIM = 'MANUEL' 'CHPO' QP12P13 1 'T' TC ;
  1705.     TNLIM = TNLIM '+' ('MANUEL' 'CHPO' QP24P34 1 'T' TH) ;
  1706.     TN = ('MANUEL' 'CHPO' QDOM1 1 'T' T0) ;
  1707.     TN = TN '+' TNLIM '-' ('REDU' TN (QP12P13 'ET' QP24P34)) ;
  1708.     DOMT = QDOM1 ;
  1709.  'FINS' ;
  1710. *
  1711. * Plot of IC
  1712. *
  1713.  'SI' GRAPH ;
  1714.     VECN = 'VECTEUR' 0.5 UN ;
  1715.     'TRACER' QDOM1 VECN ('CONTOUR' QDOM1) 'TITRE' 'Speed';
  1716.     'TRACER' TN DOMT 'TITRE' 'Temperature' ;
  1717.  'FINSI' ;
  1718. * 
  1719. ********************************
  1720. ***** Constant Matrices ********
  1721. ********************************
  1722. * (with respect to time)
  1723. *
  1724. * Mass matrix for the speed
  1725. *
  1726. * (ux,Nv)
  1727. * (uy,Nv)
  1728. *             
  1729.  MMASV = (MATMAS QDOM1 'UX' 'UX' 'QUAF' 'QUAF') 'ET'
  1730.          (MATMAS QDOM1 'UY' 'UY' 'QUAF' 'QUAF') ;
  1731.  MMASVSDT = MMASV '/' DT ;
  1732.  MMASIK = 'KOPS' 'RIMA' MMASVSDT ;
  1733. * 
  1734. *
  1735. * Mass matrix for the temperature
  1736. *
  1737. * (T,Nt)
  1738. *
  1739. *             
  1740.  MMAST = (MATMAS QDOM1 'T' 'T' DISCT DISCT) ;
  1741.  MMASTSDT = MMAST '/' DT ;
  1742. *
  1743. * (\div u, Np) 
  1744. *        
  1745.  MDIV = MATDIV QDOM1 NOMPRE 'LINE' 'QUAF' ;
  1746. *
  1747. * Matrix for the gravity source term in the speed equation
  1748. *
  1749. * (g (rho0/rho '-' 1), Nv) j
  1750. *
  1751. * where *  rho0    = P0 '/' (R '*' T0)
  1752. *       * (1/rho)  = (R '*' TN) '/' PTHER 
  1753. *
  1754. * N.B. rho0/rho is a CHPOINT with one component ('T'), with the same
  1755. *      discretisation as TN (disct).
  1756. *
  1757. *                     nomp nomd discp discd 
  1758.  MSOUGRA = MATMAS QDOM1 'T' 'UY' DISCT 'QUAF' ;
  1759. *
  1760. * Matrix for DPTHSDT source term in the temperature equation
  1761. *
  1762. * ( (1/(rho cp)) \der{PTHER}{t}, Nt)
  1763. * where * Nt is the test function for the temperature
  1764. *       * \der{PTHER}{t} ia a FLOTTANT
  1765. *       * (1/(rho cp)) = (R '*' TN) '/' (cp PTHER)
  1766. *         is a CHPOINT with one component ('T'),
  1767. *         with the same discretisation as TN (disct)
  1768. * 
  1769.  MSOUTEM = MATMAS QDOM1 'T' 'T' DISCT DISCT ; 
  1770. *
  1771. * In the divergence equation, we have
  1772. *
  1773. * \div(u) = (-1/(gamma * PTHER)) '*' \der{PTHER}{t}
  1774. *           '+'
  1775. *           ((gamma '-' 1) '/' (gamma '*' PTHER)) \div(lambda \grad(T))
  1776. *
  1777. * Discretisation is perform by the scalar product of the equation with
  1778. * Np
  1779. *
  1780. * Matrix to compute the contribution of \div(lambda \grad(T)) to
  1781. * the divergence of the speed.
  1782. *
  1783. * Surface contributions in QP12P13 (cold inlet), QP24P34 (cold inlet).
  1784. * At the oulet, we impone natural boundary conditions. Elsewhere, we
  1785. * have adiabatic surfaces
  1786. *
  1787. *
  1788. * \int_s lambda (grad T) Np \cdot n dS
  1789. * Here (lambda grad T) is the primal variable 
  1790. *
  1791. *'DEBPROC' MATFLU ;
  1792. *'ARGUMENT'  _st*'MAILLAGE' ;
  1793. *'ARGUMENT'  discp*'MOT     ' ;
  1794. *'ARGUMENT'  nomd*'MOT     ' ;
  1795. *'ARGUMENT'  discd*'MOT     ' ;
  1796. *'ARGUMENT'  discco*'MOT     ' ;
  1797. *'ARGUMENT'  coef*'CHPOINT' ;
  1798. *'ARGUMENT'  discn*'MOT     ' ;
  1799. *'ARGUMENT'  nx*'CHPOINT' ;
  1800. *'ARGUMENT'  ny*'CHPOINT' ;
  1801. *
  1802. *
  1803.  NORMX = 'MANUEL' 'CHPO' DOM1 1 'SCAL'  0.0 ;
  1804.  NORMY = 'MANUEL' 'CHPO' DOM1 1 'SCAL' -1.0 ;
  1805.  COEFBID = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
  1806.  MSMDCOLD = MATFLU QP12P13 'QUAF' NOMPRE 'LINE' 'LINE' COEFBID
  1807.           'LINE' NORMX NORMY ;
  1808.  
  1809.  NORMX = 'MANUEL' 'CHPO' DOM1 1 'SCAL'   1.0 ;
  1810.  NORMY = 'MANUEL' 'CHPO' DOM1 1 'SCAL'   0.0 ;
  1811.  COEFBID = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
  1812.  MSMDHOT = MATFLU QP24P34 'QUAF' NOMPRE 'LINE' 'LINE' COEFBID
  1813.           'LINE' NORMX NORMY ;
  1814. *
  1815. * - (\int lambda (grad T) \cdot (grad Np) dV) =
  1816. *   -1. '*' MATLAP
  1817. * 
  1818. *'DEBPROC' MATLAP ;
  1819. *'ARGUMENT'  _mt*'MAILLAGE'    ;
  1820. *'ARGUMENT'   nomt*'MOT     '  ;
  1821. *'ARGUMENT'  disct*'MOT     '  ;
  1822. *'ARGUMENT'   nomd*'MOT     '  ;
  1823. *'ARGUMENT'  discd*'MOT     '  ;
  1824. *'ARGUMENT'  discal*'MOT     ' ;
  1825. *'ARGUMENT'  alpha*'CHPOINT'   ;
  1826. *'ARGUMENT'  discof*'MOT     ' ;
  1827. *'ARGUMENT'  coef*'CHPOINT'    ;
  1828. *
  1829. *
  1830.  COEFBID = 'MANUEL' 'CHPO' DOM1 'SCAL' 1.0 ;
  1831.  MSMDVO = (MATLAP   QDOM1 'T' DISCT NOMPRE 'LINE' DISCLA CHLAM
  1832.          'LINE' COEFBID) '*' -1. ;
  1833. *
  1834. * Matrix to compute the contribution of \div(lambda \grad(T)) to
  1835. * the divergence condition. This matrix has as primal variable both
  1836. * t ('T') and grad T ('UX', 'UY')!
  1837. *
  1838.  MSMD = MSMDCOLD 'ET' MSMDHOT 'ET' MSMDVO ;
  1839. *
  1840. * Matrix to compute the contribution of
  1841. *
  1842. * (-1/(gamma * PTHER)) '*' \der{PTHER}{t}
  1843. *
  1844.  MSMDIV = MATMAS QDOM1 'SCAL' NOMPRE 'LINE' 'LINE' ;
  1845. *
  1846. * Test: DIV U = -u0 * LINJ
  1847. *
  1848.  'SI' VRAI ;
  1849.      uinj = 'SOMT' (MDIV '*' UN) ;
  1850.      uinj1 = d1 * (Uc0 '+' Uh0) ;
  1851.      ERRO = (UINJ1 '+' uinj) 'ABS' ;
  1852.      'SI' (ERRO > 1.0D-10) ;
  1853.          'ERREUR' 21 ;
  1854.      'FINSI' ;
  1855.  'FINS' ;
  1856. *
  1857. *************************************************************************
  1858. * Temporal loop *********************************************************
  1859. *************************************************************************
  1860. *
  1861.  TAB1 = 'TABLE' ;
  1862.  TAB1.'TITRE'= TABLE ;
  1863.  TAB1.1='TIRR ';
  1864.  TAB1.2='TIRC ';
  1865.  TAB1.3='REGU ';
  1866.  TAB1.'TITRE' . 1  = MOT 'Mass cons.';
  1867.  TAB1.'TITRE' . 2  = MOT 'Div. cons.';
  1868.  TAB1.'TITRE' . 3  = MOT 'Analyt.';
  1869.  TITP = 'CHAINE' 'Pressure from BC.' ;
  1870. *
  1871.  TPS = 0.0 ;
  1872. * 
  1873.  LIT   = 'PROG' ;
  1874.  LTPS  = 'PROG' ;
  1875.  LER   = 'PROG' ;
  1876.  LERT  = 'PROG' ;
  1877.  LERP  = 'PROG' ;
  1878.  LPTHER  = 'PROG' ;
  1879.  ERRO = 1.0D10 ;
  1880.  ERROT = 1.0D10 ;
  1881. *
  1882.  TABRESTN = 'TABLE' ;
  1883.  TABRESTN . 'XINIT' = TN ;
  1884.  TABRESUN = 'TABLE' ;
  1885.  'SI' LOGFI ;
  1886.     TABRESUN . 'XINIT' = UN '+' PN ;
  1887.  'SINON' ;
  1888.     TABRESUN . 'XINIT' = UN ;
  1889.     TABRESDU = 'TABLE' ;
  1890.     TABRESDU . 'XINIT' = UN '*' 0.0 ;
  1891.     TABRESDP = 'TABLE' ;
  1892.     TABRESDP . 'XINIT' = (PN '*' 0.0) ;
  1893.  'FINSI' ;    
  1894. *
  1895.  INDITE = 0 ;
  1896. *
  1897.  UNM = 'COPIER' UN ;
  1898.  TNM = 'COPIER' TN ;
  1899. *
  1900. *****************************************
  1901. *********** TIME LOOP *******************
  1902. *****************************************
  1903. * 
  1904.  'REPETER' BL1 NITER ;
  1905.  UNM2 = 'COPIER' UNM ;
  1906.  TNM2 = 'COPIER' TNM ; 
  1907.  UNM = 'COPIER' UN ;
  1908.  TNM = 'COPIER' TN ;
  1909.  PTHERM2 = PTHERM ;
  1910.  PTHERM  = PTHER ;
  1911. *
  1912. **** Time step
  1913. *
  1914.  TPS = TPS '+' DT ;
  1915. *
  1916. **** Pressure and its time variation
  1917. *
  1918.  PTHER = P0 ;
  1919.  'SI' (('EGA' &BL1 1) 'OU' LOGBDF1) ;
  1920.     DPTHSDT = (PTHER '-' PTHERM) '/' DT ;
  1921.  'SINON' ;
  1922.     DPTHSDT = ((3 '*' (PTHER '-' PTHERM)) '-'
  1923.       (PTHERM '-' PTHERM2)) '/' (2. '*' DT) ;
  1924.  'FINSI' ;
  1925. *
  1926. *
  1927. **** Speeds vary during the time iterations
  1928. *
  1929.     rhoC = PTHER '/' (R '*' TC) ;
  1930.     rhoH = PTHER '/' (R '*' TH) ;
  1931. *
  1932.     uC = mC '/' rhoC ;
  1933.     uH = mH '/' rhoH ;
  1934. *
  1935.     UN3 = UNGEOMC * uc ;
  1936.     UN3 = 'NOMC' 'UY' UN3 ;
  1937.     UN4 = UNGEOMH * uh ;
  1938.     UN4 = 'NOMC' 'UX' UN4 ;
  1939.     UNLIM = UN1 '+' UN2 '+' UN3 '+' UN4 ;
  1940. *    
  1941. *** We have to correct the speed also at the outlet
  1942. *
  1943.     UNLIM1 = 'REDU' UN ('EXTRAIRE' UNLIM 'MAILLAGE') ;
  1944.     UN = UN '+' (UNLIM '-' UNLIM1) ;
  1945. * 
  1946. **** INTERNAL ITERATIONS ****************
  1947. *
  1948.  'REPETER' BLINT  NINT ;
  1949. *      
  1950.     INDITE = INDITE '+' 1 ;
  1951. *
  1952. ************************************************************************
  1953. **** Divergence constraint and correction of the speed at the outlet ***
  1954. ************************************************************************
  1955. *
  1956. *   We compute qdiv = (div(u),Npj), Npj being the form function
  1957. *   for the pressure.
  1958. *
  1959. *   First of all, let us evaluate
  1960. *
  1961. *   coef (\div (lambda \grad T), Npj)
  1962. *
  1963.     'SI' LBOUS ;
  1964.        coef = 0.0 ;
  1965.     'SINON' ;
  1966.        coef = (gamma '-' 1) '/' (gamma '*' PTHER) ;
  1967.     'FINSI' ;
  1968. *    
  1969. *   Surface integral contribution
  1970. *   CHPFLU
  1971. *   'ARGUMENT'  _mt*'MAILLAGE' ;
  1972. *   'ARGUMENT'  disct*'MOT     ' ;
  1973. *   'ARGUMENT'  T*'CHPOINT'     ;
  1974. *   'ARGUMENT'  discd*'MOT     ' ;
  1975. *   'ARGUMENT'  discla*'MOT     ' ;
  1976. *   'ARGUMENT'  lambda*'CHPOINT' ;
  1977. *
  1978.     lamgt = CHPFLU QDOM1 DISCT TN 'QUAF' DISCLA CHLAM ;
  1979.     qdiv = MSMD '*' ((coef '*' TN) '+' (coef '*' lamgt)) ;
  1980. *
  1981. *   Since we are in a open domain, we must know
  1982. *   1/(gamma P) (dP/dt)
  1983. *
  1984.     chcoef = 'MANUEL' 'CHPO' DOM1 1 'SCAL'
  1985.        (-1 '*' DPTHSDT '/' (gamma '*' PTHER)) ;
  1986.     qdiv = qdiv '+' (MSMDIV '*' chcoef) ;
  1987. *
  1988. *   Finally, we correct the speed at the outlet in order to satisfy
  1989. *   \sum divu
  1990. *
  1991.     cacca = 'SOMT' (qdiv '-' (MDIV '*' UN)) ;
  1992.     cacca = cacca '/' d1 ;
  1993. *   cacca is the increment of average speed which exit from QP12P13
  1994.     UN5 = UNGEOMO * cacca  ;
  1995.     UN5 = 'NOMC' 'UX' UN5 ;
  1996.     UN = UN '+' UN5 ;
  1997.     UNLIM = UNLIM '+' UN5 ;
  1998. *
  1999. *   Test: we check whether (MDIV '*' UN) = qdiv
  2000. *
  2001.     cacca = 'SOMT' (qdiv '-' (MDIV '*' UN)) ;
  2002.     'SI' (('ABS' cacca) > 1.0D-12) ;
  2003.        VECN = 'VECTEUR' 0.5 UN ;
  2004.        'TRACER' QDOM1 VECN ('CONTOUR' QDOM1) 'TITRE' 'Speed';
  2005.        'ERREUR' 21 ;
  2006.     'FINSI' ;
  2007. *
  2008. *   We save the values of PTHER, UN, TN, PN to compute
  2009. *   the convergence error
  2010. *
  2011.     LOGER = ('EGA' &BLINT 1) 'OU'
  2012.           ('EGA' ((&BLINT '/' NFREQ) '*' NFREQ) &BLINT) ;
  2013. *          
  2014.     'SI' LOGER ;
  2015.        UNL = 'COPIER' UN ;
  2016.        PNL = 'COPIER' PN ;
  2017.        TNL = 'COPIER' TN ;
  2018.     'FINSI' ;
  2019. *
  2020.     UNTX = ('EXCO' 'UX' UN) ;
  2021.     UNTY = ('EXCO' 'UY' UN) ;
  2022. *    
  2023. **********************************************************************
  2024. *** Tensorial numerical viscosity ************************************
  2025. **********************************************************************
  2026. *
  2027.     'SI' LOGUP ;
  2028.        UMOD = (UNTX '*' UNTX) '+' (UNTY '*' UNTY) ;
  2029.        UMOD = UMOD '**' 0.5 ;
  2030. *      We add epsilon to avoid UMOD = 0    
  2031.        UMODEP = UMOD '+' (1.0D-6 '*' u0) ;
  2032. *
  2033. *   We compute \alpha = \frac{\lambda}{\rho c_p} =
  2034. *              \frac{(\gamma '-' 1) T}{\gamma P} \lambda
  2035. *
  2036. *
  2037.        ALPHA = ((gamma '-' 1) '/' (gamma '*' PTHER)) '*'
  2038.           ('NOMC' TN 'SCAL') '*' CHLAM ;
  2039. *   The reciprocal of Peclet
  2040.        USPE = ALPHA '/' (UMODEP '*' HUP) ;
  2041. *   The reciprocal of Reynolds
  2042.        USRE = Pr '*' USPE ;
  2043. *
  2044.        GT = -1 '*' (USPE '-' 0.5) ;
  2045.        GT = GT '*' ('MASQUE' GT 'SUPERIEUR' 0.0) ;
  2046.        GT = GT '*' (HUP '/' UMODEP) ;
  2047. *
  2048.        GV = -1 '*' (USRE '-' 0.5) ;
  2049.        GV = GV '*' ('MASQUE' GV 'SUPERIEUR' 0.0) ;
  2050.        GV = GV '*' (HUP '/' UMODEP) ;
  2051. *
  2052. *      The tensorial matrix
  2053. *
  2054.        NUXX = UNTX '*' UNTX ;
  2055.        NUYY = UNTY '*' UNTY ;
  2056.        NUXY = UNTX '*' UNTY ;
  2057.        NUYX = NUXY ;
  2058. *       
  2059.        'SI' LOGSCAL ;
  2060.           NU1 = (NUXX '+' NUXY) ;
  2061.           NU2 = (NUXX '-' NUXY) ;
  2062.           NU1 = NU1 '+' (NU2 'ABS') ;
  2063.           NU  = 0.5 '*' NU1 ;
  2064.           NU1 = (NU '+' NUYX) ;
  2065.           NU2 = (NU '-' NUYX) ;
  2066.           NU1 = NU1 '+' (NU2 'ABS') ;
  2067.           NU  = 0.5 '*' NU1 ;
  2068.           NU1 = (NU '+' NUYY) ;
  2069.           NU2 = (NU '-' NUYY) ;
  2070.           NU1 = NU1 '+' (NU2 'ABS') ;
  2071.           NU  = 0.5 '*' NU1 ;
  2072.           NUXX = NU ;
  2073.           NUXY = 0 '*' NU ;
  2074.           NUYX = 0 '*' NU ;
  2075.           NUYY = NU ;
  2076.        'FINSI' ;
  2077.        'SI' LOGDIA ;
  2078.           NUXY = 0 '*' NUXY ;
  2079.           NUYX = 0 '*' NUYX ;
  2080.        'FINSI' ;
  2081. *    
  2082. *   Diffusion matrix for the temperature
  2083. *
  2084.        MLAPTD = MATLAPT QDOM1 'T' DISCT 'T' DISCT 'QUAF'
  2085.           (GT '*' NUXX) (GT '*' NUXY)  (GT '*' NUXY) (GT '*' NUYY)
  2086.           'QUAF' (CHUN '*' BETAUP) ;
  2087.        MLAPVD = (MATLAPT QDOM1 'UX' 'QUAF' 'UX' 'QUAF' 'QUAF'
  2088.           (GV '*' NUXX) (GV '*' NUXY)  (GV '*' NUXY) (GV '*' NUYY)
  2089.           'QUAF' (CHUN '*' BETAUP)) 'ET'
  2090.           (MATLAPT QDOM1 'UY' 'QUAF' 'UY' 'QUAF' 'QUAF'
  2091.           (GV '*' NUXX) (GV '*' NUXY)  (GV '*' NUXY) (GV '*' NUYY)
  2092.           'QUAF' (CHUN '*' BETAUP)) ;
  2093.     'FINSI' ;
  2094. *   
  2095. *******************
  2096. *** Temperature ***
  2097. *******************
  2098. *
  2099. *   ((u \cdot \grad) T^{n+1}, Nt), LHS
  2100. *
  2101.     MKONT = MATKT QDOM1 'QUAF' UNTX UNTY 'T' DISCT ;
  2102. *
  2103. *   (-coef \div (alpha \grad T), Nt), LHS
  2104. *   Contribution of the integral on the volume, LHS 
  2105. *   The contribution of the integral of surface have
  2106. *   to be specified (if it is not 0)
  2107. *
  2108. *   coef = 1 '/' (rho '*' cp) =
  2109. *        = ((gamma '-' 1) '*' TN) '/' (gamma '*' PTHER)
  2110. *
  2111. *   NB coef = CHPOINT with one component ('SCAL') on DISCT
  2112. *
  2113.     DISCOF = DISCT ;
  2114.     'SI' LBOUS ;
  2115. *
  2116. *   In Boussinesq approx coef = 1 '/' (rho0 '*' cp)
  2117. *
  2118.        coef = ((gamma '-' 1) '/' (gamma '*' p0)) '*' T0 ;
  2119.        COEFL = 'MANUEL' 'CHPO' DOMT 1 'SCAL' coef ;
  2120.     'SINON' ;
  2121.        COEFL = ((gamma '-' 1) '/' (gamma '*' pther)) '*'
  2122.                ('NOMC' TN 'SCAL') ;
  2123.     'FINSI' ;
  2124.     MLAPT = (MATLAP QDOM1 'T' DISCT 'T' DISCT DISCLA CHLAM DISCOF
  2125.              COEFL) ;
  2126. *
  2127. *   Source term
  2128. *
  2129. *   (coef * Tn,Lt) = coef (Tn,Lt) RHS
  2130. *    
  2131. *   where coef = ((gamma '-' 1) '*' DPTHSDT) '/' (gamma '*' PTHER)
  2132. *              =FLOTTANT
  2133. *
  2134.     DISCOF = DISCT ;
  2135.     'SI' LBOUS ;
  2136.        COEFL = 'MANUEL' 'CHPO' DOMT 1 'T' 0.0 ;
  2137.     'SINON' ;
  2138.        COEFL = (((gamma '-' 1) '*' DPTHSDT) '/' (gamma '*' PTHER))
  2139.            '*' TN ;
  2140.     'FINSI' ;
  2141.     qtt   = MSOUTEM '*' COEFL ;
  2142.     'SI' (('EGA' &BL1 1) 'OU' LOGBDF1) ;
  2143. *
  2144. *** First time iteration: BDF1     
  2145. *
  2146. *   ((1 '/' dt) T^{n}, Nt), RHS
  2147. *
  2148.        qtt  = qtt '+' (MMASTSDT '*' TNM) ;
  2149.        mtott = mmastsdt 'ET' mkont 'ET' mlapt ;
  2150.     'SINON' ;
  2151.        qtt = qtt '+' (MMASTSDT '*' ((2 '*' TNM) '-' (0.5 '*' TNM2))) ;
  2152.        mtott = (1.5 '*' mmastsdt) 'ET' mkont 'ET' mlapt ;
  2153.     'FINSI' ;
  2154. *
  2155. *   Temperature: LHS, RHS:
  2156. *
  2157.     'MESSAGE' ;
  2158.     'MESSAGE' 'Computation of T' ;
  2159.     'SI' LOGUP ;
  2160.        mtott = mtott 'ET' mlaptd ;
  2161.     'FINSI' ;
  2162.     smbt = qtt ;    
  2163.     mtotik = 'KOPS' 'RIMA' mtott ;
  2164.     LOGPREC = (&BLINT 'EGA' 1) 'OU' (INDITE 'EGA' 1) 'OU'
  2165.        ('EGA' METINV 'DIRECT') ;
  2166.     'SI' LOGPREC ;
  2167.        mprect = mtotik ;
  2168.     'SINON' ;
  2169.        'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
  2170.              (ERROT > ERROPRE)) ;
  2171.           mprect = mtotik ;
  2172.           'MESSAGE' 'Computation of preconditioner' ;
  2173.        'FINSI' ;
  2174.     'FINSI' ;
  2175.     TN = RESOUP FAUX mtotik mprect smbt (TNLIM) 1.0D-60 NITERT
  2176.          PREC TABRESTN METINV ;
  2177. *
  2178. *************
  2179. *** Speed ***
  2180. *************
  2181. *    
  2182. *   ((1 '/' dt) ux^{n}, Nv), RHS
  2183. *   ((1 '/' dt) uy^{n}, Nv), RHS
  2184. *
  2185.     'SI' (('EGA' &BL1 1) 'OU' LOGBDF1) ;
  2186. *
  2187. *** First time iteration: BDF1     
  2188. *
  2189.        qtv =  (MMASVSDT '*' UNM) ;
  2190.        mtotv = MMASVSDT ;
  2191.     'SINON' ;
  2192.        qtv =  (MMASVSDT '*' ((2 '*' UNM) '-' (0.5 '*' UNM2))) ;
  2193.        mtotv = 1.5 '*' MMASVSDT ;
  2194.     'FINSI' ;
  2195. *
  2196. *
  2197. *   ((u \cdot \grad) ux^{n+1}, Nv), LHS
  2198. *   ((u \cdot \grad) uy^{n+1}, Nv), LHS
  2199. *
  2200.     MKONV = (MATKV QDOM1 'QUAF' 'QUAF' UNTX UNTY) ;
  2201. *    
  2202. *   (-coef \div (nu \tau), Nv), LHS
  2203. *   contribution of the integral on the volume, LHS 
  2204. *   Contribution of the integral of surface have to be specified
  2205. *   (if it is not 0).
  2206. *
  2207. *
  2208. *   coef =  1/rho = (R '/' PTHER) '*' TN
  2209. *
  2210.     DISCOF = DISCT ;
  2211.     'SI' LBOUS ;
  2212. *
  2213. *      coef = RT0/P0
  2214. *
  2215.        coef = (R '/' P0) '*' T0 ;
  2216.        COEFL = 'MANUEL' 'CHPO' DOMT 1 'SCAL' coef ;
  2217.     'SINON' ;
  2218.        COEFL = (R '/' PTHER) '*' ('NOMC' TN 'SCAL') ;
  2219.     'FINSI' ;
  2220.     MLAPV = MATTAU QDOM1 'QUAF' DISCMU DISCOF CHMU COEFL ;
  2221. *
  2222. *   (coef \dep{p}{x}, Nv), LHS  (RHS in projection)
  2223. *   (coef \dep{p}{y}, Nv), LHS  (RHS in projection)
  2224. *   contribution of the integral on the volume, RHS 
  2225. *   Contribution of the integral of surface have to be specified (if it
  2226. *   is not 0).
  2227. *
  2228. *   NB coef of gradP and div tau are the same
  2229. *   But gradP is on the RHS, tau in the LHS
  2230. *
  2231.     MPRE = MATPRE QDOM1 NOMPRE 'LINE' 'QUAF' DISCOF COEFL ;
  2232.     'SI' ('NON' LOGFI) ;
  2233.        qtv = qtv '+' (MPRE '*' (-1 '*' PN)) ;
  2234.     'FINSI' ;
  2235. *
  2236. *   (g ((rho_0 '/' rho) '-' 1) , Nv) j, RHS
  2237. *    
  2238.     'SI' LBOUS ;
  2239. *
  2240. *       (g (T/T0 '-' 1) , Nv) j, RHS
  2241. *
  2242.         CHPO1 = ((TN '/' T0) '-' 1.) '*' g ;
  2243.     'SINON' ;
  2244.         CHPO1 = ((TN '*' (P0 '/' (T0'*' PTHER))) '-' 1.) '*' g ;
  2245.     'FINSI' ;
  2246.     qtv = qtv '+' (MSOUGRA '*' CHPO1) ;
  2247. *
  2248. *   Speed: LHS, RHS
  2249. *
  2250.     'SI' LOGFI ;
  2251.        'MESSAGE' ;
  2252.        'MESSAGE' 'Computation of U' ;
  2253.        mtotv = mtotv 'ET' MKONV 'ET' MLAPV 'ET' MPRE 'ET' MDIV ;
  2254.        'SI' LOGUP ;
  2255.           mtotv = mtotv 'ET' mlapvd ;
  2256.        'FINSI' ;
  2257. *       
  2258.        smbv  = qtv '+' qdiv ;
  2259.        mtotik = 'KOPS' 'RIMA' mtotv ;
  2260. *       mtotik = 'ET' mtotik mdiv2 ;
  2261.        'SI' LOGPREC ;
  2262.           mprecufi = mtotik ;
  2263.        'SINON' ;
  2264.           'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
  2265.                 (ERRO > ERROPRE)) ;
  2266.              mprecufi = mtotik ;
  2267.              'MESSAGE' 'Computation of preconditioner' ;
  2268.           'FINSI' ;
  2269.        'FINSI' ;
  2270.        'SI' ('EGA' ('TYPE' NITERU) ENTIER) ;
  2271.           UNP = RESOUP FAUX mtotik mprecufi smbv UNLIM 1.0D-60 NITERU
  2272.                 PREC TABRESUN METINV ;
  2273.        'SINON' ;
  2274.           UNP = RESOUP FAUX mtotik mprecufi smbv UNLIM 1.0D-16 200
  2275.                 PREC TABRESUN METINV ;
  2276.        'FINSI' ;
  2277.        UN = 'EXCO' ('MOTS' 'UX' 'UY') UNP ('MOTS' 'UX' 'UY') ;
  2278.        PN = 'EXCO' ('MOTS' NOMPRE) UNP ('MOTS' NOMPRE) ;
  2279.     'SINON' ;
  2280. *
  2281. ****   Speed without divergence constraint
  2282. *    
  2283.        'MESSAGE' ;
  2284.        'MESSAGE' 'Computation of Util' ;
  2285.        mtotv = mtotv 'ET' MKONV 'ET' MLAPV ;
  2286.        'SI' LOGUP ;
  2287.           mtotv = mtotv 'ET' mlapvd ;
  2288.        'FINSI' ;
  2289.        smbv  = qtv ;
  2290.        mtotik = 'KOPS' 'RIMA' mtotv ;
  2291.        'SI' LOGPREC ;
  2292.           mprecut = mtotik ;
  2293.        'SINON' ;
  2294.           'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
  2295.                 (ERRO > ERROPRE)) ;
  2296.              mprecut = mtotik ;
  2297.              'MESSAGE' 'Computation of preconditioner' ;
  2298.           'FINSI' ;
  2299.        'FINSI' ;
  2300.        'SI' ('EGA' ('TYPE' NITERUT) ENTIER) ;
  2301.           UNTIL = RESOUP FAUX mtotik mprecut smbv UNLIM 1.0D-60
  2302.              NITERUT PREC TABRESUN METINV ;
  2303.        'SINON' ;
  2304.           UNTIL = RESOUP FAUX mtotik mprecut smbv UNLIM 1.0D-16
  2305.              200 PREC TABRESUN METINV ;
  2306.        'FINSI' ;
  2307. *
  2308. ****   Dynamic pressure
  2309. *
  2310. *      First of all, we have to compute the laplacian
  2311. *      of the dynamic pressure increment
  2312. *
  2313. *      (-coef \div coef1 \grad P, Np), LHS
  2314. *
  2315. *      coef =  -1 
  2316. *      coef1 =  1/rho = (R '/' PTHER) '*' TN
  2317. *
  2318.        COEF = 'MANUEL' 'CHPO' DOM1 1 'SCAL' -1.0 ;
  2319.        DISCOF = DISCT ;
  2320.        'SI' LBOUS ;
  2321. *
  2322. *         coef1 = RT0/P0
  2323. *
  2324.           coef1 = (R '/' P0) '*' T0 ;
  2325.           COEF1 = 'MANUEL' 'CHPO' DOMT 1 'SCAL' coef1 ;
  2326.        'SINON' ;
  2327.           COEF1 = (R '/' PTHER) '*' ('NOMC' TN 'SCAL') ;
  2328.        'FINSI' ;
  2329.        MLAPP = (MATLAP QDOM1 NOMPRE 'LINE' NOMPRE 'LINE' DISCOF COEF1
  2330.                 'LINE'  COEF) ;
  2331.        'MESSAGE' ;
  2332.        'MESSAGE' 'Computation of deltaP' ;
  2333.        qpv = MDIV '*' UNTIL ;
  2334.        qpv = (qpv '-' qdiv) '/' DT ;
  2335.        'SI' ('NON' (('EGA' &BL1 1) 'OU' LOGBDF1)) ;
  2336. *
  2337. ********** BDF2
  2338. *
  2339.           qpv = qpv '*' 1.5 ;
  2340.        'FINSI' ;
  2341. *       
  2342.        mtotik = 'KOPS' 'RIMA' MLAPP ;
  2343.        'SI' LOGPREC ;
  2344.           mprecp = mtotik ;
  2345.        'SINON' ;
  2346.           'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
  2347.                 (ERRO > ERROPRE)) ;
  2348.              mprecp = mtotik ;
  2349.              'MESSAGE' 'Computation of preconditioner' ;
  2350.           'FINSI' ;
  2351.        'FINSI' ;
  2352.        'SI' ('EGA' ('TYPE' NITERP) ENTIER) ;
  2353.           DELTAP = RESOUP FAUX mtotik mprecp qpv PNLIM 1.0D-60 NITERP
  2354.               PREC TABRESDP METINV ;
  2355.        'SINON' ;
  2356.           DELTAP = RESOUP FAUX mtotik mprecp qpv PNLIM 1.0D-16 200 
  2357.               PREC TABRESDP METINV ;
  2358.        'FINSI' ;
  2359.        PN = PN '+' DELTAP ;
  2360. *
  2361. ****   Speed (with divergence constraint)
  2362. *
  2363.        'MESSAGE' ;
  2364.        'MESSAGE' 'Computation of deltaU' ;
  2365.        qtv = MPRE '*' DELTAP ;
  2366.        DELTAU = RESOUP FAUX  MMASIK MMASIK qtv DUNLIM 1.0D-16 200
  2367.             PREC TABRESDU METINV ;
  2368.        UN = UNTIL '-' DELTAU ;
  2369.     'FINSI' ;
  2370. *    
  2371. ***********************************************
  2372. **** Post-treatment and convergence check *****
  2373. ***********************************************
  2374. *
  2375.     'SI' LOGER ;
  2376.  
  2377.  
  2378.         ERRO  = 'MAXIMUM' (UN '-' UNL) 'ABS' ;
  2379.         ERROT = 'MAXIMUM' (TN '-' TNL) 'ABS' ;
  2380.         ERROP = 'MAXIMUM' (PN '-' PNL) 'ABS' ;
  2381.         'MESSAGE' 
  2382.  '***********************************************************' ;
  2383.         'MESSAGE' ('CHAINE' 'TEMPORAL ITER = ' &BL1
  2384.                    ', TPS = ' TPS  ',  ERRO = '  ERRO
  2385.                    ', INTIT = ' &BLINT) ;
  2386.         'MESSAGE'
  2387.  '***********************************************************' ;
  2388.         LIT = LIT 'ET' ('PROG' INDITE) ;
  2389.         L10 = 'LOG' 10. ;
  2390.         LER = LER 'ET' ('PROG' (('LOG' (ERRO '+' 1.0D-20)) '/' L10)) ;
  2391. *        EVER = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LER ;
  2392.         L10 = 'LOG' 10. ;
  2393.         LERT = LERT 'ET'
  2394.                ('PROG' (('LOG' (ERROT '+' 1.0D-20)) '/' L10));
  2395.         LERP = LERP 'ET'
  2396.                ('PROG' (('LOG' (ERROP '+' 1.0D-20)) '/' L10));
  2397.         EVERT = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LERT ;
  2398.         'SI' GRAPH ;
  2399.            TIT = 'CHAINE' 'TPS =' TPS ' , IT = ' &BL1 ' , INTIT = '
  2400.                  &BLINT ;
  2401.            'SI' (((INDITE '/' (2 '*' NFREQ)) '*'
  2402.                   (2 '*' NFREQ)) 'EGA' INDITE) ;
  2403.              'DESSIN' EVERT 'TITR'
  2404.                ('CHAINE' 'Convergence history (T). ' TIT)
  2405.                'YBOR' -16. 3  'NCLK' ;
  2406. *            'DESSIN' EVER 'TITR' 'Convergence history (v)' 'NCLK' ;
  2407.            'SINON' ;
  2408.               VUN = 'VECTEUR' 0.1 UN ;
  2409.               'TRACER' TN DOM1 VUN 'TITRE' TIT 'NCLK' ;
  2410. *              EVP  = 'EVOL' 'MANU' 'it' LTPS 'Pther' LPMA ;
  2411. *              EVP2 = 'EVOL' 'MANU' 'it' LTPS 'Pther' LPCO ;
  2412. *              'DESSIN' (EVP 'ET' EVP2) 'TITR' TITP 'LEGE' TAB1 'NCLK' ;
  2413.            'FINSI' ;
  2414.         'FINSI' ;
  2415.     'FINSI' ;
  2416.     'SI' (ERROT < EPSERR) ;
  2417.        'QUITTER' BLINT ;
  2418.     'FINSI' ;
  2419. *    
  2420.  'FIN' BLINT ;
  2421. *
  2422. ***** END OF INTERNAL ITERATIONS
  2423. *
  2424.  LTPS   = LTPS 'ET' ('PROG' TPS) ;
  2425.  LPTHER = LPTHER 'ET' ('PROG' PTHER) ;
  2426.  'FIN' BL1 ;
  2427. *
  2428. *
  2429. *************************************************************************
  2430. * End of temporal loop **************************************************
  2431. *************************************************************************
  2432. *
  2433. *
  2434. * 'OUBLIER' RIGIDITE object
  2435. *
  2436.  'OUBLIER' MFLU ;
  2437.  'OUBLIER' MTOTV   ;
  2438.  'OUBLIER' MPRE    ;
  2439.  'OUBLIER' MLAPV   ;
  2440.  'OUBLIER' MKONV   ;
  2441.  'OUBLIER' MTOTT   ;
  2442.  'OUBLIER' MLAPT   ;
  2443.  'OUBLIER' MKONT   ;
  2444.  'OUBLIER' MSMDIV  ;
  2445.  'OUBLIER' MSMD    ;
  2446.  'OUBLIER' MSMD5   ;
  2447.  'OUBLIER' MSMD4   ;
  2448.  'OUBLIER' MSMD3   ;
  2449.  'OUBLIER' MSMD2   ;
  2450.  'OUBLIER' MSMD1   ;
  2451.  'OUBLIER' MSOUTEM ;
  2452.  'OUBLIER' MSOUGRA ;
  2453.  'OUBLIER' MDIV    ;
  2454.  'OUBLIER' MMASTSDT;
  2455.  'OUBLIER' MMAST   ;
  2456.  'OUBLIER' MMASVSDT;
  2457.  'OUBLIER' MMASV   ;
  2458.  'OUBLIER' MLAPTD ;
  2459.  'OUBLIER' MLAPVD ;
  2460. *
  2461. * 'OUBLIER' MATRIK object
  2462. * 
  2463.  'OUBLIER' MTOTIK ;
  2464.  'OUBLIER' MMASIK ;
  2465.  'OUBLIER' MPRECT ;
  2466.  'OUBLIER' MPRECUT ;
  2467.  'OUBLIER' MPRECP ;
  2468.  'OUBLIER' MPRECUFI ;
  2469. *
  2470.  'MENAGE' ;
  2471. * 
  2472. * 
  2473. * 'SI' LBOUS ; 
  2474. *    'OPTION' 'SAUVER' ('CHAINE' REPSAU 'bous'
  2475. *                             'g' ('ENTIER' g) 'Re100' NX '.sauv') ;
  2476. * 'SINON' ;
  2477. *    'OPTION' 'SAUVER' ('CHAINE' REPSAU
  2478. *                             'g' ('ENTIER' g) 'Re100' NX '.sauv') ;
  2479. * 'FINSI' ;
  2480. * 'SAUVER' ;
  2481. *
  2482. * Test on the convergence error
  2483. *
  2484.  'SI' LOGFI ;
  2485.     TNFI = 'COPIER' TN ;
  2486.     aa = 'IPOL' (NINT '*' 0.9999) ('EXTR' evert 'ABSC')
  2487.        ('EXTR' evert 'ORDO') ;
  2488.     'SI' (aa > -5.) ;
  2489.        'MESSAGE' 'FI: convergerce error inacceptable' ;
  2490.        'MESSAGE' ('CHAINE' 'erro = ' aa) ;
  2491.        'ERREUR' 5 ;
  2492.     'FINSI' ;
  2493.     aa = 'IPOL' (2 '*' NINT '*' 0.9999) ('EXTR' evert 'ABSC')
  2494.        ('EXTR' evert 'ORDO') ; 
  2495.     'SI' (aa > -5.) ;
  2496.        'MESSAGE' 'FI: convergerce error inacceptable' ;
  2497.        'MESSAGE' ('CHAINE' 'erro = ' aa) ;
  2498.        'ERREUR' 5 ;
  2499.     'FINSI' ;
  2500.  'SINON' ;
  2501.     TNPRO = 'COPIER' TN ;
  2502.     aa = 'IPOL' (NINT '*' 0.9999) ('EXTR' evert 'ABSC')
  2503.        ('EXTR' evert 'ORDO') ;
  2504.     'SI' (aa > -4.) ;
  2505.        'MESSAGE' 'PROJ: convergerce error inacceptable' ;
  2506.        'MESSAGE' ('CHAINE' 'erro = ' aa) ;
  2507.        'ERREUR' 5 ;
  2508.     'FINSI' ;
  2509.     aa = 'IPOL' (2 '*' NINT '*' 0.9999) ('EXTR' evert 'ABSC')
  2510.        ('EXTR' evert 'ORDO') ; 
  2511.     'SI' (aa > -4.) ;
  2512.        'MESSAGE' 'PROJ: convergerce error inacceptable' ;
  2513.        'MESSAGE' ('CHAINE' 'erro = ' aa) ;
  2514.        'ERREUR' 5 ;
  2515.     'FINSI' ;
  2516.  'FINSI' ;
  2517.  
  2518.  'FIN' BLMET ;
  2519. *
  2520. * Test on the results obatined via the projection and the fully implicit
  2521. * approach.
  2522. *
  2523.  ERRO = 'MAXI' ((TNPRO - TNFI) '/' T0) 'ABS' ;
  2524.  'SI' (ERRO > 1.0D-6) ;
  2525.     'MESSAGE' 'FI and PROJ give different results!' ;
  2526.     'MESSAGE' ('CHAINE' 'erro = ' erro) ;
  2527.  'FINSI' ;
  2528.  'FIN' ; 
  2529.  
  2530.  
  2531.  
  2532.  
  2533.  
  2534.  
  2535.  
  2536.  
  2537.  

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