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

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