Télécharger unsat_lindiriEFMH.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : unsat_lindiriEFMH.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ********************************************************************
  5. GRAPH = faux ;
  6. ******************** CAS TEST : unsat_lindiri.dgibi **********************
  7. *
  8. * Test de fonctionnement de DARCYSAT en 1D avec effet de gravité
  9. * en regime transitoire.
  10. *
  11. *-------------------------------------------------------------------
  12. *
  13. * A. GENTY - DM2S/SFME/LSET - 03/2009
  14. *
  15. * ==================================================================
  16. * Infiltration d'eau a pression d'eau imposee (h2<=0) depuis la
  17. * surface dans un milieu 1D non sature limite par une surface
  18. * inferieure a pression d'eau (h0<=0) imposee.
  19. *
  20. * h = h2 (charge H = h2 + L)
  21. * z=L ---------------------------
  22. * I I
  23. * I I
  24. * I I
  25. * I I
  26. * I hi = solution du regime I
  27. * I permanent avec h1 I
  28. * I impose en z=L I
  29. * I (charge Hi = hi + z) I
  30. * flux nul I I flux nul
  31. * I I
  32. * I I
  33. * I I
  34. * I I
  35. * I I
  36. * I I
  37. * z=0 ---------------------------
  38. * x=0 h = h0 (charge H = h0) x=x0
  39. *
  40. *
  41. * La solution analytique (en pression) du probleme est obtenue
  42. * par linearisation (Cf. Srivastava et Yeh, 1991, WRR, Vol 27,
  43. * N 5, pp 753-762) puis par separation des variables.
  44. * ==================================================================
  45. *
  46. * - conditions initiales :
  47. * > la pression initiale hi correspond a la solution du probleme
  48. * en regime permanent avec h=h1 impose en z=L.
  49. *
  50. * - conditions aux limites :
  51. * > la pression sur la frontiere inferieure du domaine est imposee
  52. * a h = h0 (charge H = h + z = h0)
  53. * > la pression sur la frontiere superieure du domaine est imposee
  54. * a h = h2 (charge H = h2 + L)
  55. * > les autres faces exterieures du domaine sont imposees a flux nul
  56. *
  57. * ===================================================================
  58. * Les options de modélisation declarées dans la table transmise à
  59. * la procédure DARCYSAT sont les suivantes :
  60. *
  61. * - les effets gravitationnels sont pris en compte :
  62. * FORCE_GRAVITE = 1.
  63. *
  64. * - Le calcul est réalisé avec l'option de calcul de pas de temps
  65. * fixe dt = 0.1 h sur une duree totale 8 h.
  66. *
  67. * - Les options numeriques testees sont les suivantes:
  68. * EFMH DECENTRE TYPINV = 3 PRECOND = 5
  69. *
  70. * - La précision de convergence demandée est de 1e-05 m
  71. *
  72. * ===================================================================
  73. * RESULTATS ATTENDUS
  74. *
  75. * Le champs de pression sera compare a la solution analytique a
  76. * 8 heures par l intermediare des erreurs L_infini et L_2.
  77. * Une erreur L2 inferieure a 1.2 10-2 est attendue.
  78.  
  79. * ===================================================================
  80. *
  81. 'OPTION' 'ECHO' 1 ;
  82. 'SAUTER' 'PAGE';
  83. *
  84. 'TITRE' 'Transitoire homogene 1D : Dirichlet' ;
  85. OPTI DIME 2 ELEM QUA4 ;
  86. OPTI ISOV SULI ;
  87. *================================================
  88. * Parametres du calcul
  89. *---------------------
  90. * longueur milieu (m)
  91. LONG1 = 1.D0 ;
  92. * largeur du milieu (m) (libre car pheno 1D)
  93. LARG1 = 1.D0 ;
  94. * alpha (m-1) pour la loi exponentielle
  95. ALPHA1 = 10.D0 ;
  96. * permeabilite a saturation (m/s)
  97. KSAT1 = 0.01D0 / 3600.D0 ;
  98. * porosite (-)
  99. PORO1 = 0.4D0 ;
  100. * teneur en eau residuelle (-)
  101. TRESID1 = 0.06D0 ;
  102. * pression imposee en bas (m)
  103. PSI0 = 0.D0 ;
  104. * pression imposee en haut initiale (m)
  105. PSI1 = -1.1D0 ;
  106. * pression imposee en haut t>0 (m)
  107. PSI2 = -1.7D0 ; ;
  108. *================================================
  109. *-------------------------------------------------------------------
  110. *--------------------- Création du maillage 2D ---------------------
  111. *
  112. *- Niveau de rafinement
  113. *- Nombre de mailles = 50 * NIV1
  114. *- Pas de temps = (Temps final) / (20 * (NIV1**2))
  115. NIV1 = 2 ;
  116. *- Densités de mailles
  117. NDX1 = 1 ;
  118. NDY1 = NIV1 * 50 ;
  119. *
  120. *- Cotes du domaine rectangulaire (x0 et z0)
  121. XCA1 = LARG1 ;
  122. YCA1 = LONG1 ;
  123. *
  124. *- Création des points du domaine
  125. A1 = 0. 0. ;
  126. A2 = XCA1 0. ;
  127. A3 = XCA1 YCA1 ;
  128. A4 = 0. YCA1 ;
  129. *
  130. *- Création des lignes
  131. LIN1 = DROIT NDX1 A1 A2 ;
  132. LIN2 = DROIT NDY1 A2 A3 ;
  133. LIN3 = DROIT NDX1 A3 A4 ;
  134. LIN4 = DROIT NDY1 A4 A1 ;
  135. *
  136. *- Faces droite et gauche
  137. LDROIT = LIN2 ;
  138. LGAUCH = LIN4 ;
  139. *
  140. *- Création du domaine
  141. MASSIF0 = DALLER LIN1 LIN2 LIN3 LIN4 PLAN ;
  142. LHAUT = LIN3 COUL BLEU ;
  143. LBAS = LIN1 COUL ROUG ;
  144. *
  145. *-Tracé du maillage avec faces sup et inf pour conditions limites
  146. ELIM 0.00001 (MASSIF0 ET LBAS ET LHAUT ET LGAUCH ET LDROIT) ;
  147. 'SI' GRAPH ;
  148. TRAC CACH (MASSIF0 ET LBAS ET LHAUT);
  149. 'FINSI' ;
  150. *
  151. ***************************************************************************
  152. ***************************************************************************
  153. *
  154. *- Création des maillages contenant tous les points
  155. QFTOT = 'CHANGER' MASSIF0 'QUAF' ;
  156. QFSORB = 'CHANGER' LBAS 'QUAF' ;
  157. QFSORH = 'CHANGER' LHAUT 'QUAF' ;
  158. *
  159. 'ELIMINATION' 0.00001 (QFTOT 'ET' QFSORB 'ET' QFSORH) ;
  160. *
  161. *-------------------------------------------------------------------
  162. *----------------------- Création du Modèle ------------------------
  163. *
  164. MODHYB = 'MODELE' QFTOT 'DARCY' 'ISOTROPE' ;
  165. MODSORB = 'MODELE' QFSORB 'DARCY' 'ISOTROPE' ;
  166. MODSORH = 'MODELE' QFSORH 'DARCY' 'ISOTROPE' ;
  167. CESORB = 'DOMA' MODSORB 'CENTRE' ;
  168. CESORH = 'DOMA' MODSORH 'CENTRE' ;
  169. HYCEN = 'DOMA' MODHYB 'CENTRE' ;
  170. HYFAC = 'DOMA' MODHYB 'FACE' ;
  171. HYSOM = 'DOMA' MODHYB 'SOMMET' ;
  172. *
  173. ZSS = COOR 2 HYSOM ;
  174. ZCC = COOR 2 HYCEN ;
  175. ZFF = COOR 2 HYFAC ;
  176. *-------------------------------------------------------------------
  177. *--------------- Conditions aux limites (blocages) ------------------
  178. *
  179. LICALC = 'PROG' 0.D0 1.D8 ;
  180. LIST1 = 'PROG' 2 * 1.D0 ;
  181. *
  182. *- surface inferieure
  183. *--------------------
  184. * Pression imposee P = HN0
  185. HN0 = PSI0 ;
  186. * On impose la charge H = P + z (z=0 en bas)
  187. ESORB = 'MANU' 'CHPO' CESORB 1 'TH' HN0 'NATURE' 'DISCRET';
  188. VALI0 = 'CHAR' 'THIM' (ESORB) ('EVOL' 'MANU' LICALC LIST1) ;
  189. *
  190. *- face superieure
  191. *-----------------
  192. *- Pression imposee initial P = HN1
  193. HN1 = PSI1 + LONG1 ;
  194. *- Pression imposee P = HN2
  195. HN2 = PSI2 + LONG1 ;
  196. HN2P = HN2 + LONG1 ;
  197. * On impose la charge H = P + z (z=L en haut)
  198. ESORH = 'MANU' 'CHPO' CESORH 1 'TH' HN2P 'NATURE' 'DISCRET';
  199. VALI2 = 'CHAR' 'THIM' (ESORH) ('EVOL' 'MANU' LICALC LIST1) ;
  200. *
  201. *--------------------------------------------------------------------
  202. *----------------- initialisation des inconnues ---------------------
  203. *
  204. *- charge initiale (en metres d'eau) dans le domaine (H = P + z)
  205. * Pression h
  206. ZAD1 = ZCC * ALPHA1 ;
  207. *
  208. L1 = ALPHA1 * LONG1 ;
  209. CST1 = EXP(ALPHA1 * HN0) ;
  210. CST2A = EXP(ALPHA1 * HN1) ;
  211. CST2B = EXP(ALPHA1 * HN2) ;
  212. CST3 = 1.0D0 - (EXP(-1.0D0*L1)) ;
  213. CST4 = (CST1 - CST2A) / CST3 ;
  214. CST5 = (CST1 - CST2B) / CST3 ;
  215. *
  216. EAP1 = ((EXP(-1. * ZAD1)) - 1.0D0) * CST4 + CST1 ;
  217. PM = (LOG(EAP1)) / ALPHA1 ;
  218. *
  219. * Traces de controle condition initiale en pression
  220. * valeur aux noeuds
  221. ZAD2 = ZSS * ALPHA1 ;
  222. EAP2 = ((EXP(-1. * ZAD2)) - 1.0D0) * CST4 + CST1 ;
  223. EAP3 = ((EXP(-1. * ZAD2)) - 1.0D0) * CST5 + CST1 ;
  224. PM2 = LOG(EAP2) * (1. / ALPHA1) ;
  225. PM3 = LOG(EAP3) * (1. / ALPHA1) ;
  226. *
  227. CHP2 = ELNO MODHYB PM ;
  228. 'SI' GRAPH ;
  229. TRAC CHP2 MASSIF0 ;
  230. 'FINSI' ;
  231. *
  232. CHP1 = KCHA MODHYB PM 'CHAM' ;
  233. 'SI' GRAPH ;
  234. TRAC MODHYB CHP1 ;
  235. 'FINSI' ;
  236. *
  237. PIN1 = REDU CHP2 LIN2 ;
  238. EVIN0 = EVOL BLEU 'CHPO' PIN1 LIN2 ;
  239. 'SI' GRAPH ;
  240. DESS EVIN0 ;
  241. 'FINSI' ;
  242. *
  243. PINI1 = REDU PM2 LIN2 ;
  244. PINI2 = REDU PM3 LIN2 ;
  245. EVIN1 = EVOL VERT 'CHPO' PINI1 LIN2 ;
  246. EVIN2 = EVOL ROUG 'CHPO' PINI2 LIN2 ;
  247. 'SI' GRAPH ;
  248. DESS (EVIN1 ET EVIN0 ET EVIN2) ;
  249. 'FINSI' ;
  250. *
  251. *
  252. * Charge H
  253. HM = PM '+' ZCC;
  254. *
  255. H0 = nomc (HM) 'H' 'NATURE' 'DISCRET';
  256. *
  257. *- Flux
  258. QFACE0 = MANU CHPO HYFAC 1 'FLUX' 0.D0 'NATURE' 'DISCRET' ;
  259. *
  260. *--------------------------------------------------------------------
  261. *----------------------------- Calcul -------------------------------
  262. *
  263. * ---------------------------
  264. * = Table DARCY_TRANSITOIRE =
  265. * ---------------------------
  266. *- initialisation table
  267. SATUR = 'TABLE' ;
  268. SATUR. 'TEMPS' = 'TABLE' ;
  269. SATUR. 'TRACE_CHARGE' = 'TABLE' ;
  270. SATUR. 'CHARGE' = 'TABLE' ;
  271. SATUR. 'FLUX' = 'TABLE' ;
  272. *
  273. *- données géométriques
  274. SATUR. 'SOUSTYPE' = 'DARCY_TRANSATUR' ;
  275. SATUR. 'MODELE' = MODHYB ;
  276. *
  277. *- instant initial
  278. SATUR. 'TEMPS' . 0 = 0. ;
  279. SATUR. 'CHARGE' . 0 = H0 ;
  280. SATUR. 'FLUX' . 0 = QFACE0 ;
  281. *
  282. *- flux impose
  283. *SATUR. 'FLUX_IMPOSE' = VALI1 ;
  284. ***********************************************************
  285. * SEULS GROS CHGTS GBM
  286. ***********************************************************
  287. *
  288. *- conditions aux limites et chargements
  289. * GBM - MODIFIE
  290. SATUR . 'TRACE_IMPOSE' = VALI0 ET VALI2 ;
  291. *
  292. * GBM - NOUVEAU
  293. SATUR . 'LUMP' = FAUX ;
  294. SATUR . 'TYPDISCRETISATION' = 'EFMH' ;
  295. SATUR . 'DECENTR' = VRAI ;
  296. *
  297. * gravite
  298. SATUR . 'FORCE_GRAVITE' = 1.D0 ;
  299. *
  300. ************************************************************
  301. *- données physiques
  302. ************************************************************
  303. *- Loi de permeabilite relative (exponentielle de h)
  304. LoiP = 'TABLE' 'EXPONENTIELLE';
  305. LoiP. 'ALPHA' = ALPHA1 ;
  306. LoiP. 'PERMSAT' = KSAT1 ;
  307. LoiP. 'COEF_N' = 1.0 ;
  308. LoiP. 'COEF_C' = 1.0 ;
  309. SATUR.'LOI_PERMEABILITE' = 'TABLE' LoiP ;
  310. *
  311. *****************************************************************
  312. *****************************************************************
  313. *- Loi de succion (exponentielle de h)
  314. LOIS = 'TABLE' 'EXPONENTIELLE' ;
  315. LOIS.'PORO' = PORO1 ;
  316. LOIS.'TERESIDU' = TRESID1 ;
  317. LOIS.'COEF_N' = 1.0 ;
  318. LOIS.'COEF_C' = 1.0 ;
  319. LOIS.'BHETA' = ALPHA1 ;
  320. SATUR . 'LOI_SATURATION' = 'TABLE' LOIS ;
  321. *
  322. *****************************************************************
  323. *****************************************************************
  324. *- données numériques
  325. *
  326. * Temps de calcul final 8 heures
  327. TF1 = 8.0D0 * 3600.D0 ;
  328. SATUR. 'TEMPS_FINAL' = TF1 ;
  329. SATUR. 'HOMOGENEISATION' = 'CHAINE' 'CENTRE' ;
  330. SATUR. 'RESIDU_MAX' = 1.D-5 ;
  331. SATUR. 'NITER' = 50 ;
  332. SATUR. 'MESSAGES' = 2 ;
  333. *
  334. * SOLVEUR Grad conj, ILU0
  335. TABRES = table METHINV;
  336. TABRES . 'TYPINV' = 3;
  337. TABRES . 'PRECOND' = 5;
  338. *
  339. SATUR . 'METHINV' = TABRES;
  340. SATUR . 'METHINV' . RESID = 1.D-25;
  341. *
  342. NBRPAS1 = (20. * (NIV1**2)) ;
  343. DELTAT1 = TF1 / NBRPAS1 ;
  344. SATUR . 'TEMPS_CALCULES' = ('PROG' 0. 'PAS' 1. NBRPAS1) * (DELTAT1) ;
  345. *
  346. *********************************************************************
  347. *------------------------- Exécution procédure ----------------------
  348. *********************************************************************
  349. *
  350. DARCYSAT SATUR ;
  351. *
  352. *********************************************************************
  353. *
  354. *-Legende des traces
  355. TABI1 = TABLE ;
  356. TABI1.'TITRE'= TABLE ;
  357. TABI1.1 = 'REGU MARQ CARR' ;
  358. TABI1.'TITRE' . 1 = MOT 'Initial 0 h' ;
  359. TABI1.2 = 'REGU MARQ TRIA' ;
  360. TABI1.'TITRE' . 2 = MOT 'Calcul 8 h' ;
  361. TABI1.3 = 'REGU MARQ LOSA' ;
  362. TABI1.'TITRE' . 3 = MOT 'Analytique 8 h' ;
  363. TABI1.4 = 'REGU MARQ TRIB' ;
  364. TABI1.'TITRE' . 4 = MOT 'Permanent' ;
  365. *
  366. *-------------------------------------------------------------------
  367. *-------------------------- Extration des charges ---------------------
  368. NDI1 = (('DIME' SATUR . TEMPS) '-' 1) ;
  369. * Dates de traces
  370. T44 = SATUR . TEMPS . NDI1 ;
  371. CALCTH4 = SATUR . CHARGE . NDI1 ;
  372. *
  373. *-------------------------------------------------------------------
  374. *------------------- Calcul et trace de la pression ----------------
  375. CPRE4 = CALCTH4 '-' ('NOMC' 'H' ZCC) ;
  376. *
  377. CALOO4 = KCHA MODHYB CPRE4 'CHAM';
  378. *
  379. SAVCAL = NOMC SCAL CPRE4 ;
  380. *
  381. * Trace sur elements
  382. 'SI' GRAPH ;
  383. TRAC MODHYB CALOO4 ;
  384. 'FINSI' ;
  385. *
  386. * Trace sur noeuds
  387. CALNO4 = ELNO MODHYB CPRE4 ;
  388. CONT1 = CONTOUR MASSIF0 ;
  389. *
  390. 'SI' GRAPH ;
  391. TRAC CALNO4 (MASSIF0) CONT1 ;
  392. 'FINSI' ;
  393. *
  394. * Trace 1D sur les noeuds
  395. VL4D = REDU CALNO4 LIN2 ;
  396. *
  397. EVOL4 = EVOL VERT 'CHPO' VL4D LIN2 ;
  398. *
  399. 'SI' GRAPH ;
  400. DESS (EVIN1 ET EVOL4) LEGE TABI1 ;
  401. DESS (EVIN1 ET EVOL4 ET EVIN2) ;
  402. 'FINSI' ;
  403. *
  404. *--------------------------------------------------------------------------
  405. * =========================================================================
  406. * Solution analytique Richards lineaire Dirichlet
  407. * u = exp(alpha.h)
  408. * u(Z,t) = us(Z) + Somme ( Bn . sin(n.Pi.Z/L) . exp(-Z/2) . exp(-p.p.t) )
  409. *--------------------------------------------------------------------------
  410. *
  411. * recuperation des coordonnees Z aux centres
  412. ZCCZ1 = COOR 2 HYCEN ;
  413. *
  414. * constantes utiles
  415. DATE1 = T44 ;
  416. TAD1 = (ALPHA1 * KSAT1 * DATE1) / (PORO1 - TRESID1) ;
  417. KZER1 = LoiP. 'PERMSAT' ;
  418. ALPHA1 = ABS((LoiP. 'ALPHA')) ;
  419. HZER1 = HN0 ;
  420. LNEW1 = ALPHA1 * LONG1 ;
  421. *
  422. * Z addimentionnel
  423. ZCC1 = ZCCZ1 * ALPHA1 ;
  424. *
  425. *----------------------------------------------------------------------------
  426. * Calcul du terme stationaire us(Z)
  427. * us(Z) = TT1(Z) = exp(alpha*h0) -
  428. * (1-exp(-z))*(exp(alpha*h0)-exp(alpha*h2))/(1-exp(-L))
  429. *----------------------------------------------------------------------------
  430. CC1 = EXP(ALPHA1*HN0) ;
  431. CC2 = EXP(ALPHA1*HN1) ;
  432. CC3 = EXP(ALPHA1*HN2) ;
  433. CC4 = EXP(-1.D0 * LNEW1) ;
  434. CUZ1 = (CC3 - CC1) / (1.D0 - CC4) ;
  435. TT1 = CC1 + (CUZ1 * (1.D0 - (EXP(-1.D0 * ZCC1)))) ;
  436. CSTB1 = (CC2 - CC3) / (1.D0 - CC4) ;
  437. BN1 = 2.D0 / LNEW1 ;
  438. *
  439. *----------------------------------------------------------------------------
  440. * calcul du terme Somme
  441. * Somme (Bn sin(n.Pi.z/L) exp(-z/2) exp(-p*p*t)
  442. * Remarque: l'indice n=0 n'est pas calcule dans la somme
  443. *----------------------------------------------------------------------------
  444. * constantes trigo pour calcul
  445. PI1 = PI ;
  446. CSD1 = 180.D0 / PI1 ;
  447. *
  448. * calcul de exp(-Z/2)
  449. TLZ1 = KOPS ZCC1 '*' -0.5D0 ;
  450. CSTE2 = EXP(TLZ1) ;
  451. *
  452. CSL1 = EXP(0.5D0 * LNEW1) ;
  453. CSL2 = (CC3 - CC2) ;
  454. NBOUC1 = 1000 ;
  455. NNN1 = 1.D0 ;
  456. SOM1 = 0.D0 ;
  457. *---------------------------------------------
  458. * calcul de la somme
  459. *---------------------------------------------
  460. REPETER BOU1 NBOUC1 ;
  461. *
  462. * Calcul des L_n et somme en sin(L_n.Z),
  463. * avec ici L_n = n.Pi/L
  464. *
  465. MESS 'Indice de INTG' NNN1 ;
  466. *******************************
  467. LM1 = NNN1 * PI1 / LNEW1 ;
  468. *
  469. *******************************
  470. * calcul de SI1 = sin(L_n Z)
  471. MU1Z = KOPS ZCC1 '*' LM1 ;
  472. * en degres
  473. MU2Z = MU1Z * CSD1 ;
  474. SI1 = SIN(MU2Z) ;
  475. ***************************************
  476. * calcul de exp(-t p**2)
  477. * avec p**2 = 1/4 + L_n**2
  478. LM2 = LM1 * LM1 ;
  479. PCAR1 = 0.25D0 + LM2 ;
  480. TEX1 = -1.D0 * TAD1 * PCAR1 ;
  481. NP1 = EXP(TEX1) ;
  482. ***************************************
  483. * calcul de Bn
  484. * avec Bn = -(2/L) . (exp(alpha.h1)-exp(alpha.h2)) .
  485. * ((-1)**n) . L_n . exp(L/2) / p**2
  486. PMU1 = (-1.D0) ** NNN1 ;
  487. BN2 = (PMU1 * LM1) / (0.25D0 + LM2) ;
  488. ***************************************
  489. * calcul produit total et somme
  490. COEF1 = SI1 * NP1 * BN2 ;
  491. SOM1 = SOM1 + COEF1 ;
  492. NNN1 = NNN1 + 1. ;
  493. FIN BOU1 ;
  494. ***************************************
  495. * fin de calcul du terme somme
  496. *-----------------------------------------------
  497. *
  498. * Calcul de u
  499. KTOT2 = BN1 * SOM1 * CSTE2 * CSL2 * CSL1 ;
  500. KTOT1 = KTOT2 + TT1 ;
  501. *
  502. **************************************************************
  503. *
  504. * calcul de h associe a u
  505. VPHK2 = LOG(KTOT1) ;
  506. VCH1 = VPHK2 / ALPHA1 ;
  507. * Trace h sur elements
  508. HFI1E = KCHA MODHYB VCH1 'CHAM';
  509. 'SI' GRAPH ;
  510. TRAC MODHYB HFI1E ;
  511. 'FINSI' ;
  512. * Trace h sur noeuds
  513. HFI1N = ELNO MODHYB VCH1 ;
  514. 'SI' GRAPH ;
  515. TRAC HFI1N (MASSIF0) CONT1 ;
  516. 'FINSI' ;
  517. *
  518. **************************************************************
  519. VL8D = REDU HFI1N LIN2 ;
  520. EVOL8 EVOL BLEU 'CHPO' VL8D LIN2 ;
  521. *
  522. 'SI' GRAPH ;
  523. DESS (EVIN1 ET EVOL4 ET EVOL8 ET EVIN2) LEGE TABI1 ;
  524. 'FINSI' ;
  525. **************************************************************
  526. *
  527. * Calcul de l'erreur relative par element (en pression)
  528. error = ((SAVCAL '-' VCH1) '/' VCH1) 'ABS' ;
  529. error = 'KCHA' modhyb error CHAM;
  530. 'SI' GRAPH ;
  531. 'TRACER' modhyb error;
  532. 'FINSI' ;
  533. *
  534. * Erreur Linfini
  535. NUME1 = ABS(SAVCAL '-' VCH1) ;
  536. DENU1 = ABS(VCH1) ;
  537. DENU2 = MAXI (DENU1) ;
  538. ERR1 = NUME1 / DENU2 ;
  539. ER1C = KCHA MODHYB ERR1 'CHAM';
  540. 'SI' GRAPH ;
  541. TRAC MODHYB ER1C ;
  542. 'FINSI' ;
  543. *
  544. ERINF1 = MAXI (ERR1) ;
  545. MESS 'ERREUR L INFINI' ERINF1;
  546. *
  547. **********************
  548. * Calcul de l'erreur L2 par element (en pression)
  549. ERRL21 = NUME1 ** 2.D0 ;
  550. ERRL22 = DENU1 ** 2.D0 ;
  551. ERRL2 = (ERRL21 '/' ERRL22) ** 0.5D0 ;
  552. ERRL3 = 'KCHA' MODHYB ERRL2 CHAM ;
  553. 'SI' GRAPH ;
  554. TRAC MODHYB ERRL3 ;
  555. 'FINSI' ;
  556. *
  557. * Calcul erreur L2
  558. SURF1 = 'DOMA' MODHYB 'VOLUME' ;
  559. ERNUM1 = ERRL21 * SURF1 ;
  560. ERDEN1 = ERRL22 * SURF1 ;
  561. INT1 = RESU ERNUM1 ;
  562. INT2 = RESU ERDEN1 ;
  563. ERL2 = (MAXI (INT1 / INT2)) ** 0.5 ;
  564. MESS 'ERREUR L2' ERL2 ;
  565. *
  566. **************************************************************
  567. * Test de non regression sur erreur L2
  568. ERMAX1 = 0.012 ;
  569. SI (ERL2 > ERMAX1) ;
  570. mess 'ERREUR ANORMALE' ERL2 '>' ERMAX1 ;
  571. 'ERREUR' 5;
  572. SINON ;
  573. 'ERREUR' 0;
  574. FINSI ;
  575. *
  576. *-------------------------------------------------------------------
  577. *
  578. FIN;
  579.  
  580.  
  581.  
  582.  
  583.  
  584.  
  585.  

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