Télécharger la plaquette Cast3M


fr en

Télécharger srivastava1VF.dgibi

Retour à la liste

Numérotation des lignes : 

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

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