Télécharger nlin_burger.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : nlin_burger.dgibi
  2. *************************************************************************
  3. *
  4. * NOM : nlin_burger.dgibi
  5. *
  6. * DESCRIPTION : We compute the 2D Burger problem.
  7. * Scalar numerical viscosity is added.
  8. *
  9. * LANGAGE : GIBIANE-CAST3M
  10. * AUTEUR : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  11. * mél : beccantini@cea.fr
  12. *************************************************************************
  13. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  14. * en cas de modification de ce sous-programme afin de faciliter
  15. * la maintenance !
  16. *************************************************************************
  17. *************************************************************************
  18. *************************************************************************
  19. ******* PERSONAL PROCEDURES *********************************************
  20. *************************************************************************
  21. *************************************************************************
  22. *
  23. *BEGINPROCEDUR matmas
  24. *
  25. *************************************************************************
  26. *
  27. * NOM : MATMAS
  28. *
  29. * DESCRIPTION : We compute
  30. *
  31. * (Np , Nd )
  32. *
  33. * LANGAGE : GIBIANE-CAST3M
  34. * AUTEUR : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  35. * mél : beccantini@cea.fr
  36. *************************************************************************
  37. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  38. * en cas de modification de ce sous-programme afin de faciliter
  39. * la maintenance !
  40. *************************************************************************
  41. *
  42. *
  43. 'DEBPROC' MATMAS ;
  44. *
  45. 'ARGUMENT' _mt*'MAILLAGE' ;
  46. 'ARGUMENT' nomp*'MOT ' ;
  47. 'ARGUMENT' nomd*'MOT ' ;
  48. 'ARGUMENT' discp*'MOT ' ;
  49. 'ARGUMENT' discd*'MOT ' ;
  50. *
  51. * Arguments
  52. * _mt = MAILLAGE, quaf mesh
  53. * nomp = MOT, name of the primal variable
  54. * nomd = MOT, name of the dual variable
  55. * discp = MOT, name of the discretization of the primal variable
  56. * discd = MOT, name of the discretization of the dual variable
  57. *
  58. idim = 'VALEUR' 'DIME' ;
  59. *
  60. * Matrix A for the operator NLIN
  61. *
  62. numop = 1 ;
  63. numder = idim ;
  64. numvar = 1 ;
  65. * The variable involved in A is the primal variable T
  66. numdat = 0 ;
  67. numcof = 0 ;
  68. *
  69. discg = 'LINE' ;
  70. * Geometrical discretisation
  71. methgau = 'GAU7' ;
  72. * Gauss points involved
  73. *
  74. A = ININLIN numop numvar numdat numcof numder ;
  75. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomp ;
  76. A . 'VAR' . 1 . 'DISC' = discp ;
  77. *
  78. * Matrix B for the operator NLIN
  79. *
  80. B = ININLIN numop numvar 0 0 numder ;
  81. * For B the coef no coefficients are involved
  82. * but 1
  83. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
  84. B . 'VAR' . 1 . 'DISC' = discd ;
  85. *
  86. A . 1 . 1 . 0 = 'LECT' ;
  87. B . 1 . 1 . 0 = 'LECT' ;
  88. *
  89. mat = 'NLIN' discg _mt A B methgau ;
  90. 'RESPRO' mat ;
  91. 'FINPROC' ;
  92. *
  93. *ENDPROCEDUR matmas
  94. *BEGINPROCEDUR matkv
  95. *
  96. *************************************************************************
  97. *
  98. * NOM : MATKV
  99. *
  100. * DESCRIPTION : We compute the matrix
  101. *
  102. * ((cx \dep{u_x}{x} '+' cy \dep{u_x}{y}), Nvx)
  103. * ((cx \dep{u_y}{x} '+' cy \dep{u_y}{y}), Nvy)
  104. *
  105. * where Nv = shape function for (u_x,u_y)^T
  106. *
  107. * Names of primal variables = names of the dual variables
  108. * = 'UX' 'UY'
  109. *
  110. * LANGAGE : GIBIANE-CAST3M
  111. * AUTEUR : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  112. * mél : beccantini@cea.fr
  113. *************************************************************************
  114. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  115. * en cas de modification de ce sous-programme afin de faciliter
  116. * la maintenance !
  117. *************************************************************************
  118. *
  119. 'DEBPROC' MATKV ;
  120. 'ARGUMENT' _mt*'MAILLAGE' ;
  121. 'ARGUMENT' discv*'MOT ' ;
  122. 'ARGUMENT' discc*'MOT ' ;
  123. 'ARGUMENT' cx*'CHPOINT' ;
  124. 'ARGUMENT' cy*'CHPOINT' ;
  125. *
  126. * _mt = solid quaf mesh
  127. * discv = discretization type ux, uy
  128. * discc = discretization type for cx, cy
  129. * cx = CHPOINT ('SCAL')
  130. * cy = CHPOINT ('SCAL')
  131. *
  132. idim = 'VALEUR' 'DIME' ;
  133. *
  134. numop = 2 ;
  135. numder = idim ;
  136. numvar = 2 ;
  137. * Two variables, ux, uy
  138. numdat = 2 ;
  139. * cx, cy
  140. numcof = 2 ;
  141. * Coef: cx cy (two identity functions)
  142. discg = 'LINE' ;
  143. * Geometrical discretisation
  144. methgau = 'GAU7' ;
  145. * Gauss points involved
  146. *
  147. * Matrix A for the operator NLIN
  148. *
  149. A = ININLIN numop numvar numdat numcof numder ;
  150. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  151. A . 'VAR' . 1 . 'DISC' = discv ;
  152. A . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  153. A . 'VAR' . 2 . 'DISC' = discv ;
  154. * cx
  155. A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  156. A . 'DAT' . 1 . 'DISC' = discc ;
  157. A . 'DAT' . 1 . 'VALEUR' = cx ;
  158. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  159. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  160. * cy
  161. A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  162. A . 'DAT' . 2 . 'DISC' = discc ;
  163. A . 'DAT' . 2 . 'VALEUR' = cy ;
  164. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  165. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  166. *
  167. * Matrix B for the operator NLIN
  168. *
  169. numvar = 2 ;
  170. B = ININLIN numop numvar 0 0 numder ;
  171. *
  172. * For B, no coefficients are involved
  173. * but 1
  174. *
  175. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
  176. B . 'VAR' . 1 . 'DISC' = discv ;
  177. B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
  178. B . 'VAR' . 2 . 'DISC' = discv ;
  179. *
  180. * Contribution
  181. *
  182. * cx \dep{ux}{x} Nvx
  183. A . 1 . 1 . 1 = 'LECT' 1 ;
  184. B . 1 . 1 . 0 = 'LECT' ;
  185. * cy \dep{ux}{y} Nvx
  186. A . 1 . 1 . 2 = 'LECT' 2 ;
  187. * cx \dep{uy}{x} Nvy
  188. A . 2 . 2 . 1 = 'LECT' 1 ;
  189. B . 2 . 2 . 0 = 'LECT' ;
  190. * cy \dep{uy}{y} Nd
  191. A . 2 . 2 . 2 = 'LECT' 2 ;
  192. *
  193. mat = 'NLIN' discg _mt A B methgau ;
  194. 'RESPRO' mat ;
  195. 'FINPROC' ;
  196. *
  197. *ENDPROCEDUR matkv
  198. *BEGINPROCEDUR matlap
  199. *
  200. *************************************************************************
  201. *
  202. * NOM : MATLAP
  203. *
  204. * DESCRIPTION : We compute the matrix arising from the computation of
  205. * the scalar product
  206. *
  207. * ( -(coef \div( \alpha \grad T)) , Nd )
  208. *
  209. * where Nd = test function for the dual variable.
  210. *
  211. * -\int_{\Omega} coef \div( \alpha \grad T) Nd dV =
  212. * -\int_{\delta \Omega}
  213. * coef \alpha (\grad T \cdot n) Nd dS '+'
  214. * +\int_{\Omega}
  215. * coef \alpha (\grad T \cdot \grad Nd) dV '+'
  216. * +\int_{\Omega}
  217. * Nd \alpha (\grad T \cdot \grad coef) dV
  218. *
  219. * The matrix is issue from the 2nd and 3rd contribution
  220. * (volume integrals).
  221. *
  222. *
  223. * LANGAGE : GIBIANE-CAST3M
  224. * AUTEUR : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  225. * mél : beccantini@cea.fr
  226. *************************************************************************
  227. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  228. * en cas de modification de ce sous-programme afin de faciliter
  229. * la maintenance !
  230. *************************************************************************
  231. *
  232. *
  233. 'DEBPROC' MATLAP ;
  234. *
  235. 'ARGUMENT' _mt*'MAILLAGE' ;
  236. 'ARGUMENT' nomt*'MOT ' ;
  237. 'ARGUMENT' disct*'MOT ' ;
  238. 'ARGUMENT' nomd*'MOT ' ;
  239. 'ARGUMENT' discd*'MOT ' ;
  240. 'ARGUMENT' discal*'MOT ' ;
  241. 'ARGUMENT' alpha*'CHPOINT' ;
  242. 'ARGUMENT' discof*'MOT ' ;
  243. 'ARGUMENT' coef*'CHPOINT' ;
  244. *
  245. * Arguments
  246. *
  247. * _mt = surface QUAF mesh
  248. * nomt = name of the primal variable T
  249. * disct = discretization type of T
  250. * nomd = name of the dual variable
  251. * discd = discretization type of the dual variable
  252. * discal = discretization type of alpha
  253. * discof = discretization type of coef
  254. * alpha = CHPOINT ('SCAL')
  255. * coef = CHPOINT ('SCAL')
  256. *
  257. idim = 'VALEUR' 'DIME' ;
  258. *
  259. * Four contributions
  260. *
  261. * (coef alpha \dep{T}{x} \dep{Nd}{x}
  262. * + coef alpha \dep{T}{y} \dep{Nd}{y}
  263. * + Nd alpha \dep{T}{x} \dep{coef}{x}
  264. * + Nd alpha \dep{T}{y} \dep{coef}{y})
  265. *
  266. * Matrix A for the operator NLIN
  267. *
  268. numop = 4 ;
  269. numder = idim ;
  270. numvar = 1 ;
  271. * The variable involved in A is the primal variable T
  272. numdat = 2 ;
  273. * Two data: alpha, coef
  274. numcof = 4 ;
  275. * Four functions: f(alpha) = alpha
  276. * f(coef) = coef
  277. * f(coef) = \dep{coef}{x}
  278. * f(coef) = \dep{coef}{y}
  279. discg = 'LINE' ;
  280. * Geometrical discretisation
  281. methgau = 'GAU7' ;
  282. * Gauss points involved
  283. *
  284. A = ININLIN numop numvar numdat numcof numder ;
  285. A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomt ;
  286. A . 'VAR' . 1 . 'DISC' = disct ;
  287. * alpha
  288. A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  289. A . 'DAT' . 1 . 'DISC' = discal ;
  290. A . 'DAT' . 1 . 'VALEUR' = alpha ;
  291. * coef
  292. A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
  293. A . 'DAT' . 2 . 'DISC' = discof ;
  294. A . 'DAT' . 2 . 'VALEUR' = coef ;
  295. * Function alpha
  296. A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
  297. A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
  298. * Function coef
  299. A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
  300. A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
  301. * Function \dep{coef}{x}
  302. A . 'COF' . 3 . 'COMPOR' = 'D/DX1' ;
  303. A . 'COF' . 3 . 'LDAT' = 'LECT' 2 ;
  304. * Function \dep{coef}{y}
  305. A . 'COF' . 4 . 'COMPOR' = 'D/DX2' ;
  306. A . 'COF' . 4 . 'LDAT' = 'LECT' 2 ;
  307. *
  308. * Matrix B for the operator NLIN
  309. *
  310. B = ININLIN numop numvar 0 0 numder ;
  311. * For B the coef no coefficients are involved
  312. * but 1
  313. B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
  314. B . 'VAR' . 1 . 'DISC' = discd ;
  315. *
  316. * Contribution
  317. *
  318. * coef alpha \dep{T}{x} \dep{Nd}{x}
  319. A . 1 . 1 . 1 = 'LECT' 1 2 ;
  320. B . 1 . 1 . 1 = 'LECT' ;
  321. * coef alpha \dep{T}{y} \dep{Nd}{y}
  322. A . 2 . 1 . 2 = 'LECT' 1 2 ;
  323. B . 2 . 1 . 2 = 'LECT' ;
  324. * Nd alpha \dep{T}{x} \dep{coef}{x}
  325. A . 3 . 1 . 1 = 'LECT' 1 3 ;
  326. B . 3 . 1 . 0 = 'LECT' ;
  327. * Nd alpha \dep{T}{y} \dep{coef}{y}
  328. A . 4 . 1 . 2 = 'LECT' 1 4 ;
  329. B . 4 . 1 . 0 = 'LECT' ;
  330. mlapn = 'NLIN' discg _mt A B methgau ;
  331. 'RESPRO' mlapn ;
  332. 'FINPROC' ;
  333. *ENDPROCEDUR matlap
  334. *BEGINPROCEDUR resoup
  335. *************************************************************************
  336. **** Resolution of a linear system **************************************
  337. *************************************************************************
  338. *
  339. 'DEBPROC' RESOUP ;
  340. 'ARGUMENT' LOGTPS*'LOGIQUE ' ;
  341. 'ARGUMENT' mat*'MATRIK ' ;
  342. 'ARGUMENT' matpre/'MATRIK ' ;
  343. 'ARGUMENT' smb*'CHPOINT ' ;
  344. 'ARGUMENT' ccl*'CHPOINT ' ;
  345. 'ARGUMENT' res*'FLOTTANT' ;
  346. 'ARGUMENT' nit*'ENTIER' ;
  347. 'ARGUMENT' pre/'ENTIER' ;
  348. 'ARGUMENT' gggtre/'TABLE ' ;
  349. 'ARGUMENT' gggtcv/'LOGIQUE ' ;
  350. 'ARGUMENT' typslv/'MOT ' ;
  351. *
  352. * Det. tt s.t.
  353. *
  354. * mat . tt = smb
  355. * tt(boundary) = ccl
  356. *
  357. *
  358. * mat = matrix to "inverse"
  359. * matpre = matrix for which preconditioner has been computed (optional)
  360. * smb = right hand side
  361. * ccl = imposed boundary conditions
  362. * res = rvk . 'RESID'
  363. * nit = max number of linear iteration (for an iterative solver)
  364. * pre = type of preconditioning matrix (for an iterative solver)
  365. * gggtre = table of resolution to store 'XINIT'
  366. * gggtcv = graphic of the history for the iterative solver
  367. * (optional FAUX)
  368. * typslv = 'DIRECT'/'ITER' (optional, ITER)
  369. *
  370. 'SI' ('NON' ('EXISTE' typslv)) ;
  371. tslv = VRAI ;
  372. 'SINON' ;
  373. 'SI' ('EGA' typslv 'DIRECT') ;
  374. tslv = FAUX ;
  375. 'SINON' ;
  376. 'SI' ('EGA' typslv 'ITER') ;
  377. tslv = VRAI ;
  378. 'SINON' ;
  379. 'ERREUR' 'Solveur DIRECT ou ITERatif ?' ;
  380. 'FINSI' ;
  381. 'FINSI' ;
  382. 'FINSI' ;
  383. gggniv = 1 ;
  384. 'SI' ('NON' ('EXISTE' pre)) ;
  385. pre = 5 ;
  386. 'FINSI' ;
  387. 'SI' ('OU' ('NON' ('EXISTE' gggtcv)) ('NON' tslv)) ;
  388. gggtcv = FAUX ;
  389. 'FINSI' ;
  390. rv = 'EQEX' ;
  391. rvk = rv . 'METHINV' ;
  392. 'SI' ('EXISTE' matpre) ;
  393. rvk . 'MATASS' = matpre ;
  394. rvk . 'MAPREC' = matpre ;
  395. 'SINON' ;
  396. rvk . 'MATASS' = mat ;
  397. rvk . 'MAPREC' = mat ;
  398. 'FINSI' ;
  399. 'SI' tslv ;
  400. rvk . 'TYPINV' = 3 ;
  401. 'SINON' ;
  402. rvk . 'TYPINV' = 1 ;
  403. 'FINSI' ;
  404. rvk . 'SCALING' = 1 ;
  405. * Scaling factor
  406. rvk . 'OUBMAT' = 1 ;
  407. * oubmat = 1 -> We eliminate the elementary matrix
  408. * 'ILUTPPIV' -> 1 We always search for a new pivot
  409. * -> 0 We do not search for a new pivot
  410. rvk . 'ILUTPPIV' = 0.1D0 ;
  411. * rvk . 'ILUTPPIV' = 0.01D0 ;
  412. * rvk . 'ILUTPPIV' = 0.00001D0 ;
  413. * rvk . 'ILUTPPIV' = 0.D0 ;
  414. rvk . 'IMPINV' = gggniv ;
  415. rvk . 'TYRENU' = 'SLOA' ;
  416. rvk . 'PCMLAG' = 'APR2' ;
  417. rvk . 'NITMAX' = nit ;
  418. rvk . 'RESID' = res ;
  419. rvk . 'PRECOND' = pre ;
  420. rvk . 'BCGSBTOL' = 1.D-200 ;
  421. rvk . 'ILUTLFIL' = 1.5 ;
  422. *
  423. 'SI' ('EXISTE' gggtre) ;
  424. 'SI' ('EXISTE' gggtre 'XINIT') ;
  425. rvk . 'XINIT' = gggtre . 'XINIT' ;
  426. 'FINSI' ;
  427. 'FINSI' ;
  428. 'TEMPS' 'ZERO' ;
  429. tt = 'KRES' mat 'TYPI' rvk 'SMBR' smb 'CLIM' ccl ;
  430. TABTPS = TEMP 'NOEC';
  431. tcpu = TABTPS.'TEMPS_CPU'.'INITIAL';
  432. tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
  433. 'SI' LOGTPS ;
  434. 'MESSAGE' ('CHAINE' 'tcpu (s) = ' tcpus) ;
  435. 'FINSI' ;
  436. 'SI' ('ET' gggtcv ('EXISTE' rvk 'CONVINV')) ;
  437. lcvg = rvk . 'CONVINV' ;
  438. nit = 'DIME' lcvg ;
  439. 'SI' ('>' nit 1) ;
  440. lit = 'PROG' 1.D0 PAS 1.D0 NPAS ('-' nit 1) ;
  441. evit = 'EVOL' 'MANU' 'iter' lit 'log residu'
  442. ('/' ('LOG' lcvg) ('LOG' 10.D0)) ;
  443. 'DESSIN' evit 'TITR' titglob ;
  444. 'FINSI' ;
  445. 'FINSI' ;
  446. 'SI' ('EXISTE' gggtre) ;
  447. gggtre . 'XINIT' = tt ;
  448. 'FINSI' ;
  449. 'RESPRO' tt ;
  450. *
  451. * End of procedure file RESOU
  452. *
  453. 'FINPROC' ;
  454. *
  455. *************************************************************************
  456. **** End of the Resolution of a linear system ***************************
  457. *************************************************************************
  458. *ENDPROCEDUR resoup
  459. *************************************************************************
  460. *************************************************************************
  461. ******* END OF PERSONAL PROCEDURES **************************************
  462. *************************************************************************
  463. *************************************************************************
  464. **
  465. *
  466. 'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
  467. GRAPH = FAUX ;
  468. *
  469. **** Mesh
  470. *
  471. *
  472. * A4 A3
  473. *
  474. *
  475. *
  476. *
  477. * A1 A2
  478. *
  479. L = 1.0 ;
  480. A1 = 0.0 0.0 ;
  481. A2 = L 0.0 ;
  482. A3 = L L ;
  483. A4 = 0.0 L ;
  484. *
  485. DX = L '/' 20 ;
  486. DY = DX ;
  487. *
  488. A1A2 = A1 'DROIT' A2 'DINI' DX 'DFIN' DX ;
  489. A2A3 = A2 'DROIT' A3 'DINI' DY 'DFIN' DY ;
  490. A3A4 = A3 'DROIT' A4 'DINI' DX 'DFIN' DX ;
  491. A4A1 = A4 'DROIT' A1 'DINI' DY 'DFIN' DY ;
  492. *
  493. DOM1 = 'SURFACE' (A1A2 'ET' A2A3 'ET' A3A4 'ET' A4A1) 'PLAN' ;
  494. PDOM1 = 'CHANGER' DOM1 'POI1' ;
  495. 'SI' GRAPH ;
  496. 'TRACER' DOM1 'TITRE' 'Mesh' ;
  497. 'FINSI' ;
  498. QDOM1 = 'CHANGER' DOM1 'QUAF' ;
  499. QA4A1 = 'CHANGER' A4A1 'QUAF' ;
  500. QA1A2 = 'CHANGER' A1A2 'QUAF' ;
  501. QA2A3 = 'CHANGER' A2A3 'QUAF' ;
  502. QA3A4 = 'CHANGER' A3A4 'QUAF' ;
  503. 'ELIMINATION' (QA4A1 'ET' QDOM1) 1.0D-4 ;
  504. 'ELIMINATION' (QA1A2 'ET' QDOM1) 1.0D-4 ;
  505. 'ELIMINATION' (QA2A3 'ET' QDOM1) 1.0D-4 ;
  506. 'ELIMINATION' (QA3A4 'ET' QDOM1) 1.0D-4 ;
  507. *
  508. * Boundary conditions
  509. *
  510. LIMDIR1 = 'DIFF' ('CHANGER' 'POI1' A4A1) ('MANUEL' 'POI1' A1) ;
  511. LIMDIR2 = 'DIFF' ('CHANGER' 'POI1' A1A2) ('MANUEL' 'POI1' A1) ;
  512. LIMDIR3 = ('MANUEL' 'POI1' A1) ;
  513. UN1 = 'MANUEL' 'CHPO' LIMDIR1 2 'UX' 1.0 'UY' 0.0 ;
  514. UN2 = 'MANUEL' 'CHPO' LIMDIR2 2 'UX' 0.0 'UY' 1.0 ;
  515. UN3 = 'MANUEL' 'CHPO' LIMDIR3 2 'UX' 1.0 'UY' 1.0 ;
  516. UNLIM = UN1 '+' UN2 '+' UN3 ;
  517. *
  518. * IC
  519. *
  520. UN = ('MANUEL' 'CHPO' DOM1 2 'UX' 0.0 'UY' 0.0) '+' UN1 '+'
  521. UN2 '+' UN3 ;
  522. *
  523. VECN = 'VECTEUR' UN ;
  524. 'SI' GRAPH ;
  525. 'TRACER' DOM1 VECN ;
  526. 'FINSI' ;
  527. *
  528. * Time step
  529. *
  530. DT = 1.0D16 '*' DX '/' 1.0 ;
  531. *
  532. NU = 'MANUEL' 'CHPO' DOM1 1 'SCAL' (0.25 '*' DX) ;
  533. CHONE = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
  534. *
  535. * Matrices
  536. *
  537. MMAS1 = MATMAS QDOM1 'UX' 'UX' 'LINE' 'LINE' ;
  538. MMAS2 = MATMAS QDOM1 'UY' 'UY' 'LINE' 'LINE' ;
  539. MMAS = MMAS1 'ET' MMAS2 ;
  540. MKON = MATKV QDOM1 'LINE' 'LINE' ('EXCO' 'UX' UN) ('EXCO' 'UY' UN) ;
  541. MLAP1 = MATLAP QDOM1 'UX' 'LINE' 'UX' 'LINE' 'LINE' NU 'LINE' CHONE ;
  542. MLAP2 = MATLAP QDOM1 'UY' 'LINE' 'UY' 'LINE' 'LINE' NU 'LINE' CHONE ;
  543. MLAP = MLAP1 'ET' MLAP2 ;
  544. *
  545. *
  546. * Total matrix
  547. *
  548. MMASSDT = MMAS '/' DT ;
  549. mconst = MMASSDT 'ET' MLAP ;
  550. mtot = mconst 'ET' MKON ;
  551. mtotm = 'KOPS' 'RIMA' mtot ;
  552. smbv pipi = 'KOPS' 'MATRIK' ;
  553. *
  554. * Resolution
  555. *
  556. * TAB1=TABLE;
  557. * TAB1.'TITRE'= TABLE ;
  558. * TAB1.1='TIRR ';
  559. * TAB1.'TITRE' . 1 = 'Exact' ;
  560. * TAB1.2='MARQ CROI';
  561. * TAB1.'TITRE' . 2 = 'Numerical' ;
  562. *
  563. TPS = 0.0 ;
  564. UN = 'COPIER' UN ;
  565. UNM = 'COPIER' UN ;
  566. LIT = 'PROG' ;
  567. LER = 'PROG' ;
  568. ERRO = 1.0D10 ;
  569. *
  570. TABRESUN = 'TABLE' ;
  571. TABRESUN . 'XINIT' = UN ;
  572. NITERU = 100 ;
  573. PREC = 5 ;
  574. METINV = MOT 'DIRECT' ;
  575. *
  576. 'REPETER' BL1 100 ;
  577. TPS = TPS '+' DT ;
  578. qv = MMASSDT '*' UN ;
  579. UN = RESOUP FAUX mtotm mtotm qv UNLIM 1.0D-60 NITERU
  580. PREC TABRESUN METINV ;
  581. 'SI' (ERRO < 1.0D-10) ;
  582. 'QUITTER' BL1 ;
  583. 'FINSI' ;
  584. 'SI' (((&BL1 '/' 5) '*' 5) 'EGA' &BL1) ;
  585. ERRO = 'MAXIMUM' (UN '-' UNM) 'ABS' ;
  586. 'MESSAGE' ('CHAINE' 'ITER = ' &BL1 ', TPS = ' TPS ', ERRO = '
  587. ERRO) ;
  588. UNM = 'COPIER' UN ;
  589. LIT = LIT 'ET' ('PROG' &BL1) ;
  590. L10 = 'LOG' 10. ;
  591. LER = LER 'ET' ('PROG' (('LOG' (ERRO '+' 1.0D-20)) '/' L10)) ;
  592. EVER = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LER ;
  593. 'SI' GRAPH ;
  594. 'DESSIN' EVER 'TITR' 'Convergence history' 'NCLK' ;
  595. 'FINSI' ;
  596. 'FINSI' ;
  597. *
  598. ****** We update the convective matrix
  599. *
  600. MKON = MATKV QDOM1 'LINE' 'LINE' ('EXCO' 'UX' UN) ('EXCO' 'UY' UN) ;
  601. MTOT = MCONST 'ET' MKON ;
  602. mtotm = 'KOPS' 'RIMA' mtot ;
  603. 'FIN' BL1 ;
  604. *
  605. UX = 'EXCO' UN 'UX' ;
  606. UY = 'EXCO' UN 'UY' ;
  607. *
  608. VECN = 'VECTEUR' UN ;
  609. 'SI' GRAPH ;
  610. 'TRACER' UX DOM1 VECN 'TITR' 'ux' ;
  611. 'TRACER' UY DOM1 VECN 'TITR' 'uy' ;
  612. 'FINSI' ;
  613. *
  614. * Test
  615. *
  616. NITER = &BL1 ;
  617. *
  618. 'SI' (NITER > 40) ;
  619. 'MESSAGE' 'Convergence not reached' ;
  620. 'ERREUR' 5 ;
  621. 'FINSI' ;
  622. *
  623. XX YY = 'COORDONNEE' DOM1 ;
  624. DOMH = (YY '-' XX) 'POIN' 'SUPERIEUR' 0.1 ;
  625. erro = 'MAXIMUM' (('EXCO' 'UX' ('REDU' UN DOMH)) '-' 1.0) 'ABS' ;
  626. 'SI' (erro > 1.0D-2) ;
  627. 'ERREUR' 5 ;
  628. 'FINSI' ;
  629. erro = 'MAXIMUM' (('EXCO' 'UY' ('REDU' UN DOMH))) 'ABS' ;
  630. 'SI' (erro > 1.0D-2) ;
  631. 'ERREUR' 5 ;
  632. 'FINSI' ;
  633. *
  634. 'FIN' ;
  635.  
  636.  
  637.  
  638.  
  639.  

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