Télécharger la plaquette Cast3M


fr en

Télécharger nlin_cavity_HP.dgibi

Retour à la liste

Numérotation des lignes : 

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

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