Télécharger darcysat.procedur

Retour à la liste

Numérotation des lignes :

  1. * DARCYSAT PROCEDUR JC220346 16/07/08 21:15:01 9008
  2. * GBM REGLER LE PROBLEME DES TRACES - VF EFMH
  3.  
  4. * DARCYSAT PROCEDUR MAUGIS 02/12/19 21:15:02 4527
  5. ***********************************************************************
  6. 'DEBPROC' DARCYSAT SATUR*'TABLE' ;
  7.  
  8. NDIME = 'VALE' DIME ;
  9. * MESSAGE
  10. * Par défaut, on affiche beaucoup d'information
  11. 'SI' ('EXISTE' SATUR 'MESSAGES' ) ;
  12. DEBUG = SATUR.'MESSAGES' > 0 ;
  13. 'SINON' ;
  14. DEBUG = VRAI ;
  15. 'FINSI' ;
  16.  
  17. *
  18. ********************************************************************
  19. *********************************************************************
  20. *********************************************************************
  21. * LECTURES *
  22. *********************************************************************
  23. *********************************************************************
  24. *********************************************************************
  25.  
  26. *
  27. * CONDITIONS INITIALES OU DU DERNIER PAS CALCULE
  28. *
  29. TPSINI CHRG TRH FLU0 LAST1 = SATUTILS LICODINI SATUR ;
  30. *
  31. * DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES
  32. *
  33. MDMH MCHYB VHYB MAILSOM MAILCENT ABMC = SATUTILS LIMODELE SATUR ;
  34. *
  35. * LOI_SATURATION
  36. *
  37. LOIS = SATUTILS LILOISAT SATUR ;
  38. *
  39. * LOI_PERMEABILITE
  40. *
  41. LOIP = SATUTILS LILOIPER SATUR ;
  42. *
  43. * param physiques
  44. *
  45. COFEMMAG L_GRAV RHOWG DAXE CONVH = SATUTILS LIPHYSIK SATUR ;
  46. *
  47. * RECUPERATION DES DONNEES NUMERIQUES
  48. *
  49. NIVEAU TETA NPAS0 NITER0 ERR0 CFL0 ITMAXI OKPENAL
  50. cofpenal npenaldt cofdiv nmaxdt = SATUTILS LIPARANU SATUR ;
  51.  
  52. NIVEAU = 'MOT' NIVEAU;
  53.  
  54. * TEMPS_CALCULES
  55. isauv = LAST1 ;
  56. DCAL LTCALCUL DTAUTO ICAL tmin NPAS0 DELTAT TFINAL isor tpsor dsor
  57. = SATUTILS LITABTPS SATUR TPSINI debug NPAS0 ;
  58. *
  59. * initialisations pour TRANGEOL
  60. *
  61. CHCLIM GEOL1 TABMODI = SATUTILS YNYTYAL SATUR TETA ;
  62. *
  63. * integrale des chargements
  64. *
  65. startflu startmix startsou FLUIMP FLUMIX TERSOU = SATUTILS YNTGFLUX
  66. SATUR ;
  67.  
  68.  
  69.  
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77. *******************************************************************
  78. * PARAMETRES DIVERS PRECALCULES
  79. *******************************************************************
  80. * Coordonnées à prendre en compte pour la gravité
  81. * (si AXE_G est nul, ZGRAV aussi, et donc pas d'effet de gravité).
  82. XAXE = 'COORDONNEE' 1 DAXE;
  83. YAXE = 'COORDONNEE' 2 DAXE;
  84. 'SI' (NDIME 'EGA' 2) ;
  85. XCO YCO = 'COOR' MAILSOM ;
  86. ZFF = 'KCHT' MDMH 'SCAL' 'SOMMET' 'COMP' 'SCAL'
  87. ((XCO * XAXE) + (YCO * YAXE));
  88. XCO YCO = 'COOR' MAILCENT ;
  89. ZCC = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL'
  90. ((XCO * XAXE) + (YCO * YAXE)) ;
  91. 'SINON' ;
  92. ZAXE = 'COORDONNEE' 3 DAXE;
  93. XCO YCO ZCO = 'COOR' MAILSOM ;
  94. ZFF = 'KCHT' MDMH 'SCAL' 'SOMMET' 'COMP' 'SCAL'
  95. ((XCO * XAXE) + (YCO * YAXE) + (ZCO * ZAXE));
  96. XCO YCO ZCO = 'COOR' MAILCENT ;
  97. ZCC = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL'
  98. ((XCO * XAXE) + (YCO * YAXE) + (ZCO * ZAXE)) ;
  99. 'FINSI' ;
  100.  
  101.  
  102. *--------------------------------------------------------------------*
  103. * RAPPEL DES DONNEES ENTREES DANS LA PROCEDURE *
  104. *--------------------------------------------------------------------*
  105.  
  106. SATUTILS rekapitu SATUR L_GRAV rhowg DAXE convh cfl0 teta LAST1 debug ;
  107.  
  108.  
  109. ***********************************************************************
  110. ***********************************************************************
  111. ***********************************************************************
  112. * RESOLUTION *
  113. ***********************************************************************
  114. ***********************************************************************
  115. ***********************************************************************
  116.  
  117.  
  118. *
  119. *--------------------------------------------------------------------*
  120. * BOUCLE RESOLVANT LE SYSTEME POUR CHAQUE PAS DE TEMPS *
  121. *--------------------------------------------------------------------*
  122. *
  123. * Initialisations :
  124. * =================
  125. * gbm a bugué sur PENAL à FAUX
  126. PENAL = OKPENAL ;
  127. TPS = TPSINI + DELTAT ;
  128. *
  129. *-- Paramètres physiques
  130.  
  131. *- Pression tronquée à utiliser dans les lois de comportement :
  132.  
  133. * tronkature pression et chgt 'unités
  134. PNS = SATUTILS TRONKP MDMH CHRG TRH L_GRAV NIVEAU RHOWG CONVH ;
  135. PNSC = SATUTILS TRONKP MDMH CHRG TRH L_GRAV 'CENTRE' RHOWG CONVH ;
  136.  
  137. *- calcul teneur en eau, saturation et premier terme capacité capillaire
  138. * Attention, la teneur calculée utilise une porosité constante au cours
  139. * du pas de temps. GBM refait à l'identique dans boucle plus loin.
  140. * GBM - ATTENTION PEUT ETRE INUTIL ICI SAUF INITIALISATION FLUX
  141. 'SI' ( 'EXISTE' SATUR 'CONSERVATIF' ) ;
  142. CONSER = SATUR .'CONSERVATIF' ;
  143. 'FINSI' ;
  144.  
  145. 'SI' ('EGA' CONSER VRAI) ;
  146. nbds SAT TENN CAPA PORO1 = SATUTILS KALSAT MDMH LOIS 'CENTRE' PNSC ;
  147. nbds SAT TEN CAPAD PORO1 = SATUTILS KALSAT MDMH LOIS NIVEAU PNS ;
  148. 'SINON' ;
  149.  
  150. nbds SAT TEN CAPA PORO1 = SATUTILS KALSAT MDMH LOIS NIVEAU PNS ;
  151. 'FINSI';
  152.  
  153.  
  154.  
  155. *- pas de temps initial automatique - rendre calcul CFL optionnel
  156. VVOL = VHYB * PORO1 ;
  157. DELTAT DTI = SATUTILS YNYTDT MDMH FLU0 MCHYB VVOL CFL0 debug DTAUTO
  158. DELTAT ;
  159.  
  160. *- Valeur des variables au pas de temps précédent
  161. * GBM PAS TOUTES UTILES - A REVOIR
  162. TRHANC = TRH ;
  163. CHRGANC = CHRG ;
  164. FLUANC = FLU0 ;
  165. PNSANC = PNS ;
  166. PNSCANC = PNSC ;
  167. TENANC = TEN ;
  168. SATANC = SAT ;
  169. TPSANC = TPSINI ;
  170. PERFANC = 0.D0 ;
  171.  
  172. *- Valeur des variables à l'itération précédente
  173. * pour la trace de charge servant au calcul du résidu,
  174. * on met n'importe quoi qui ait la bonne structure.
  175. * Ici, TRH n'a pas de multiplicateur de lagrange, donc l'estimation
  176. * de FLRES sera foireuse au tout prsavreemier calcul du résidu.
  177. TRH2N = TRH ;
  178. CAPAN = 0.D0 ;
  179. PERFN = 0.D0 ;
  180.  
  181. NOMESPL = 'H' ;
  182.  
  183. *
  184. * initialisation du terme source intégral.
  185. *
  186.  
  187. 'SI' ( 'EXISTE' SATUR 'SOURCE' ) ;
  188. TERSC2M1 = 'NOMC' (NOMESPL) ('TIRE' TERSOU TPSINI) ;
  189. 'FINSI' ;
  190.  
  191. 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ;
  192. FLUIMPM1 = 'TIRE' FLUIMP TPSINI ;
  193. 'FINSI' ;
  194.  
  195. 'SI' ('EXISTE' SATUR 'FLUMIXTE') ;
  196. FLUMMPM1 = 'TIRE' FLUMIX TPSINI ;
  197. 'FINSI' ;
  198.  
  199. FLU1 = FLU0 ;
  200. GEOL1 . 'CONCENTRATION' = CHRG ;
  201.  
  202. 'SI' ('EGA' CONSER VRAI) ;
  203. TENNC = TENN ;
  204. 'FINSI' ;
  205.  
  206. ************** NPX***************************************************
  207. *
  208. * Sauvegarde
  209. * ==========
  210.  
  211. MESS 'SAUVEGARDE INITIALE---------------------------------------' ;
  212.  
  213.  
  214.  
  215. TMP TMP2 = SATUTILS SAVRESU SATUR TPSOR ISOR TPSANC TPS -1 DELTAT
  216. TRHANC TRH CHRGANC CHRG TENANC TEN PNSCANC PNSC SATANC SAT FLUANC
  217. FLU0 NIVEAU DSOR 0 MDMH ;
  218.  
  219. ***************FIN NPX***********************************************
  220.  
  221. *======================================================== transitoire
  222.  
  223. 'SI' ('EGA' NPAS0 0) ; NPAS0 = -1 ; 'FINSI' ;
  224. 'REPETER' TRANSI NPAS0 ;
  225.  
  226. *====================================================================
  227. *
  228.  
  229. IPAS = ICAL + &TRANSI - 1 ;
  230.  
  231.  
  232. *-----------------------------------------------------------------
  233.  
  234. 'REPETER' PENALDT NPENALDT ;
  235.  
  236. CHRG = CHRGANC ;
  237.  
  238. TPS = TPSANC + DELTAT ;
  239.  
  240.  
  241. *------------------------------------------------------------------
  242. * Affichage information si debug vrai
  243. 'MESSAGE' 'deltat dti' deltat dti;
  244. SATUTILS AFFICH &PENALDT &TRANSI LAST1 TPS TPSANC DELTAT DTI
  245. debug ;
  246.  
  247. ***************** INITIALISATION PAS DE TPS **************************
  248. *
  249. *- Incorporation des CLs
  250. *
  251. 'SI' ('EXISTE' SATUR 'TRACE_IMPOSE') ;
  252. CHARIMPO = 'TIRE' SATUR . 'TRACE_IMPOSE' TPS ;
  253. CHCLIM . 'DIRICHLET' = 'NOMC' NOMESPL CHARIMPO ;
  254. 'FINSI' ;
  255.  
  256.  
  257. 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ;
  258. 'SI' (tpsanc '<EG' startflu) ;
  259. FLUIMPO = 'TIRE' FLUIMP TPS ;
  260. 'SINON' ;
  261. FLUIMPO = 'COPIER' FLUIMPM1 ;
  262. 'FINSI' ;
  263. FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET' ;
  264. 'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF')
  265. 'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH')
  266. 'ET' ('EGA' teta 1.D0))) ;
  267. * On écrit les flux sous forme intégrale sauf en explicite et
  268. * kranck-Nickholson pour EFMH. On les différentie ici
  269. FLUMP = (FLUIMPO '-' FLUIMPM1) '/' DELTAT ;
  270. 'SINON' ;
  271. FLUMP = 'COPIER' FLUIMPO ;
  272. 'FINSI' ;
  273. CHCLIM . 'NEUMANN' = 'NOMC' NOMESPL FLUMP ;
  274. 'FINSI' ;
  275.  
  276. 'SI' ('EXISTE' SATUR 'FLUMIXTE') ;
  277. 'SI' (tpsanc '<EG' startmix) ;
  278. FLUMMPO = 'TIRE' FLUMIX TPS ;
  279. 'SINON' ;
  280. FLUMMPO = 'COPIER' FLUMMPM1 ;
  281. 'FINSI' ;
  282. FLUMMPO = 'CHAN' 'ATTRIBUT' FLUMMPO 'NATURE' 'DISCRET' ;
  283. 'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF')
  284. 'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH')
  285. 'ET' ('EGA' teta 1.D0))) ;
  286. * On écrit les flux sous forme intégrale sauf en explicite et
  287. * kranck-Nickholson pour EFMH. On les différentie ici
  288. FLUMMP = ( FLUMMPO '-' FLUMMPM1) '/' DELTAT ;
  289. 'SINON' ;
  290. FLUMMP = 'COPIER' FLUMMPO ;
  291. 'FINSI' ;
  292. CHCLIM . 'FLUMIXTE' = TABLE ;
  293. CHCLIM . 'FLUMIXTE' . 'VAL' = 'NOMC' NOMESPL FLUMMP ;
  294. CHCLIM . 'FLUMIXTE' . 'COEFA' = SATUR . 'FLUMIXTE' . 'MIXCOFA' ;
  295. * GBM rappel il y aura des modifs à faire car HYDRAU
  296. CHCLIM . 'FLUMIXTE' . 'COEFB' = SATUR . 'FLUMIXTE' . 'MIXCOFB' ;
  297. 'FINSI' ;
  298.  
  299. GEOL1 . 'CLIMITES' = CHCLIM ;
  300.  
  301. *
  302. *- Calcul de la contribution des termes sources
  303. *
  304. TERSCE = 'MANU' 'CHPO' MAILCENT 1 'SOUR' 0. ;
  305. TERSCE = 'NOMC' NOMESPL TERSCE ;
  306.  
  307. * Terme source propre :
  308. 'SI' ( 'EXISTE' SATUR 'SOURCE' ) ;
  309. 'SI' (tpsans <EG startsou) ;
  310. * la source varie encore et est non nulle
  311. TERSC2 = 'TIRE' TERSOU TPS ;
  312. 'SINON' ;
  313. * la source est nulle on garde l'intégrale précédente
  314. TERSC2 = TERSC2M1 ;
  315. 'FINSI' ;
  316. TERSC2 = 'NOMC' NOMESPL TERSC2 ;
  317. TERSCV = (TERSC2 '-' TERSC2M1) '/' DELTAT ;
  318. TERSCE = TERSCV + TERSCE ;
  319. 'FINSI' ;
  320.  
  321. * Le terme source physique est tersce.
  322.  
  323. GEOL1 . 'DELTAT' = DELTAT ;
  324. GEOL1 . 'METHINV' = SATUR . 'METHINV' ;
  325.  
  326. *----------------------------------------------------- non linéaire
  327.  
  328. 'REPETER' LINEAR itmaxi ;
  329. 'SI' ('EGA' CONSER VRAI) ;
  330. GEOL1 . 'CONCENTRATION' = CHRG ;
  331. 'FINSI' ;
  332.  
  333. *__________________________________________________________________
  334.  
  335. * Préparation relaxation - GBM tester efficacité
  336. * ======================
  337. 'SI' ('EGA' &LINEAR 1) ;
  338. * pas de relaxation à la première itération
  339. cofrelax = 1.D0 ;
  340. 'SINON' ;
  341. 'SI' (&LINEAR '<EG' (itmaxi/2)) ;
  342. * on relaxe à partir du deuxième pas de temps
  343. cofrelax = SATUR.'SOUS_RELAXATION' ;
  344. 'SINON' ;
  345. * on relaxe plus si on dépasse itmaxi/2
  346. * GBM GBM GBMGBM 0.5 normalement
  347. cofrelax = 0.5D0 ;
  348. 'FINSI' ;
  349. 'FINSI' ;
  350.  
  351. * Calcul des nouveaux paramètres
  352. * ==============================
  353.  
  354. *-- Pression tronquée à utiliser dans les lois de comportement:
  355. * P1 pression pascal, PNS pression tronquée
  356. PNS = SATUTILS TRONKP MDMH CHRG TRH L_GRAV NIVEAU RHOWG CONVH ;
  357. PNSC = SATUTILS TRONKP MDMH CHRG TRH L_GRAV 'CENTRE' RHOWG CONVH ;
  358.  
  359.  
  360. * GBM reprendre calcul sat, ten et capa
  361. * mettre un flag pour calculer CAPA ou pas
  362. * retirer PNSANC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  363. * TENANC
  364.  
  365. 'SI' ('EGA' CONSER VRAI) ;
  366. nbds SATC TENC CAPA = SATUTILS KALSOT MDMH LOIS 'CENTRE' PNSC ;
  367. dum1 SAT TEN CAPAD = SATUTILS KALSOT MDMH LOIS NIVEAU PNS ;
  368. 'SINON' ;
  369. nbds SATC TEN CAPA = SATUTILS KALSOT MDMH LOIS 'CENTRE' PNSC ;
  370. dum1 SAT TEN CAPA = SATUTILS KALSOT MDMH LOIS NIVEAU PNS ;
  371. 'FINSI' ;
  372.  
  373.  
  374.  
  375. *-- Calcul perméabilité
  376. * PERF = 0. ;
  377. 'REPETER' BB ('DIME' LOIP.'INDEX') ;
  378. NOMT = 'CHAINE' LOIP.'INDEX'.&BB ;
  379. LOI1 = LOIP . NOMT ;
  380. MOD1 = LOI1.'MODELE' ;
  381. NOMKR = 'MOT' ('TEXTE' ('CHAINE' LOI1.'NOM_PROCEDURE')) ;
  382.  
  383. 'SI' (('EGA' LOI1.'SOUSTYPE' ('CHAINE' PUISSANCE) ) 'OU'
  384. ('EGA' LOI1.'SOUSTYPE' ('CHAINE' PERSONNELLE)) 'OU'
  385. ('EGA' LOI1.'SOUSTYPE' ('CHAINE' MUALEM)) 'OU'
  386. ('EGA' LOI1.'SOUSTYPE' ('CHAINE' BURDINE)) 'OU'
  387. ('EGA' LOI1.'SOUSTYPE' ('CHAINE' BROOKS_COREY)) 'OU'
  388. ('EGA' LOI1.'SOUSTYPE' ('CHAINE' MUALEM_BURDINE))) ;
  389. * réduction de la saturation au sous-domaine (local)
  390. SATL = 'REDU' SAT ('DOMA' MOD1 NIVEAU) ;
  391. PERFL = ('TEXTE' NOMKR) LOI1 SATL ;
  392. PERFL = RECENTRE PERFL MOD1 NIVEAU ;
  393. 'SINON' ;
  394. PNL = 'REDU' PNS ('DOMA' MOD1 NIVEAU) ;
  395. PERFL = ('TEXTE' NOMKR) LOI1 PNL ;
  396. PERFL = RECENTRE PERFL MOD1 NIVEAU ;
  397. 'FINSI' ;
  398.  
  399. * 'SI' ('EGA' NIVEAU 'FACE') ;
  400. * LOI1.'ABMC' = 'ABS' ('DOMA' MOD1 'ORIENTAT') ;
  401. * 'FINSI' ;
  402.  
  403.  
  404. * concaténation perméabilité
  405. LCOMP = 'EXTR' PERFL 'COMP' ;
  406. NCOMP = 'DIME' LCOMP ;
  407. 'REPETER' CONSPER NCOMP ;
  408. J = &CONSPER ;
  409. COMPI = 'EXTR' LCOMP J ;
  410. PERFINI = 'MANU' 'CHPO' ('DOMA' MDMH 'CENTRE') 1 COMPI 0.
  411. 'NATURE' 'DISCRET' ;
  412. PERFI = 'EXCO' PERFL COMPI COMPI ;
  413. PERFI = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' COMPI PERFINI PERFI;
  414. *
  415. 'SI' (&BB 'EGA' 1) ;
  416. 'SI' ('EGA' J 1) ;
  417. PERFT = PERFI ;
  418. 'SINON' ;
  419. PERFT = PERFT 'ET' PERFI ;
  420. 'FINSI' ;
  421. 'SINON' ;
  422. 'SI' ('EGA' J 1) ;
  423. PERFT = PERFT 'ET' PERFI ;
  424. 'SINON' ;
  425. PERFT = PERFT 'ET' PERFI ;
  426. 'FINSI' ;
  427. 'FINSI' ;
  428. *
  429. 'FIN' CONSPER ;
  430. *
  431. 'OUBLIER' satl ;
  432. 'OUBLIER' perl;
  433.  
  434. 'FIN' BB ;
  435.  
  436. *-- coefficient d'emmagasinement si milieu saturé
  437. * GBM ATTENTION ICI BIZARRE
  438. * GBM MET A 0 CAR DISCONTINUITE DU TERME DEVANT DT, 0 AVANT
  439. * SAT ET POSITIF APRES. DOIT ETRE MIS DANS LES LOIS DE SAT
  440. * ET EVOLUER CONTINUEMENT.
  441. EMMAG = ('MASQUE' ('NOMC' SATC 'SCAL') 'EGSUP' 0.99 ) * COFEMMAG;
  442.  
  443. * EMMAG = NOMC SCAL (SATC * COFEMMAG) ;
  444. EMMAG = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL' EMMAG ;
  445. * 'MESSAGE' 'maxi min emmag' ('MINIMUM' emmag) ('MAXIMUM' emmag);
  446.  
  447.  
  448. * Construction système matriciel
  449. * ==============================
  450.  
  451. *-- Penalisation -
  452. * GBM ERREUR sur NOMC 'SOUR'
  453. 'SI' PENAL ;
  454. TERSC1 = TERSCE + ('NOMC' NOMESPL
  455. ((CHRG - CHRGANC) * VHYB / DELTAT * cofpenal)) ;
  456. 'SINON' ;
  457. TERSC1 = TERSCE ;
  458. 'FINSI' ;
  459.  
  460.  
  461. 'SI' ('EGA' CONSER VRAI) ;
  462. TERSC1 = TERSC1 - ('NOMC' NOMESPL
  463. ((TENC - TENNC) * VHYB / DELTAT )) ;
  464. 'FINSI' ;
  465. GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSC1 ;
  466.  
  467. * relaxation
  468.  
  469.  
  470. * 'MESSAGE' 'coefrelaaaaaaaaaaaaaaaa' cofrelax;
  471. 'SI' ('NEG' cofrelax 1.D0 1.D-14) ;
  472. * on relaxe CAPA et PERF
  473. CAPAR = (CAPA '*' cofrelax)
  474. '+' (CAPAN '*' (1. - cofrelax)) ;
  475. PERFR = ( PERFT '*' cofrelax)
  476. '+' (PERFN '*' (1. - cofrelax)) ;
  477. 'SINON' ;
  478. CAPAR = CAPA ;
  479. PERFR = PERFT ;
  480. 'FINSI' ;
  481.  
  482.  
  483.  
  484.  
  485.  
  486. *-- coef devant la derivée temps
  487. COFDT = CAPAR '+' EMMAG ;
  488.  
  489. * permeabilité
  490. * LAM = ('NOMC' PERFR 'K' NATURE DIFFUS) ;
  491.  
  492. * GBM - penalisation un peu particuliere ! ecrire ce que cela
  493. * veut dire
  494. 'SI' PENAL ;
  495. COFDT = COFDT '+' COFPENAL ;
  496. 'FINSI' ;
  497.  
  498. * initialisation TRANSGEOL - gbm bien verif les tabmodi
  499. GEOL1 . 'POROSITE' = COFDT ;
  500. TABMODI . 'POROSITE' = VRAI ;
  501. GEOL1 . 'DIFFUSIVITE' = PERFR ;
  502. TABMODI . 'DIFFUSIV' = VRAI ;
  503.  
  504. * Résolution :
  505. * ============
  506.  
  507. * GBM appel trangeol - tester peut etre aussi premier penal
  508. 'SI' ( (&TRANSI 'EGA' 1)) ;
  509. GEOLPF1 GEOLPF2 = TRANGEOL MDMH GEOL1 ;
  510. 'SINON' ;
  511. GEOLPF1 GEOLPF2 = TRANGEOL MDMH GEOL1 GEOL2 ;
  512. 'FINSI' ;
  513. CHRG = GEOLPF1 . 'CONCENTRATION' ;
  514. FLU1 = GEOLPF1 . 'FLUXDIFF' ;
  515. 'SI' ('EGA' SATUR . 'TYPDISCRETISATION' EFMH) ;
  516. TRH2 = GEOLPF2 . 'TRACE_CONC' ;
  517. 'SINON' ;
  518. * GBM §§§§§§§§§§§§§§§§§§§§§§§
  519. TRH2 = 0.D0 * FLU1;
  520. 'FINSI' ;
  521.  
  522. * GBM - a clarifier
  523. TRH = TRH2 ;
  524.  
  525. * FIN DE BOUCLE :
  526. * =================
  527.  
  528.  
  529.  
  530. * test de convergence sur le résidu
  531. * =================================
  532. * A-t-on convergé à l'itération précédente (avec les nouveaux paramètres) ?
  533. * flux correspondant au résidu système résolu sans relaxation
  534.  
  535. **
  536. ** GBM reflechir comment améliorer la conservation
  537.  
  538. 'SI' (&LINEAR 'EGA' 1) ;
  539. OKCONV = FAUX ;
  540. maxer = 1.D0;
  541. 'SINON' ;
  542.  
  543. RES1 = 'MAXIMUM' ('RESULT' ('ABS' (TEN - TENN))) ;
  544. RES2 =
  545. ((0.1D0 * ('MAXIMUM' ('RESULT' ('ABS' (TEN)))))
  546. + ( ('MAXIMUM' ('RESULT' ('ABS' (TEN '-' TENANC))))));
  547.  
  548.  
  549.  
  550.  
  551. RES3 = 'MAXIMUM' ('RESULT' ('ABS' (PERFR - PERFN))) ;
  552. RES4 = ('MAXIMUM' ('RESULT' ('ABS' ((perfr '-' perfn)))))
  553. '/' ('MAXIMUM' ('RESULT' PERFR)) ;
  554. RES3 = RES3 '/' (
  555. 'MAXIMUM' ('RESULT' ('ABS' (PERFR - PERFANC))) + 1.D-100) ;
  556.  
  557.  
  558.  
  559.  
  560. * 'MESSAGE' 'minmax res1 res2' res1 res2 ('MAXIMUM'
  561. * (TEN - TENANC)) ;
  562. MAXER = RES1 '/' RES2 ;
  563. MAXER = MAXER '+' RES3 ;
  564.  
  565. * on prend en compte la taille du pas de tps comparé à CFL
  566. 'SI' (&TRANSI 'EGA' 1) ;
  567. * Pour le premier pas de temps les flux ne sont pas connus
  568. * donc le pas de temps CFL est mal évalué. On ne l'inclue
  569. * pas dans le crtiere de convergence
  570. OKCONV = MAXER < (ERR0) ;
  571. 'SINON' ;
  572. OKCONV = MAXER < (ERR0 * DTI '/' DELTAT) ;
  573. * 'MESSAGE' 'error4 err3' res4 res3 ;
  574. OKCONV = OKCONV 'ET' (RES4 < ( ERR0)) ;
  575. 'FINSI' ;
  576.  
  577. 'SI' (DTAUTO) ;
  578. 'SI' (&LINEAR 'EGA' 2) ;
  579. * on impose de boucler au moins 2 fois en pas de tps auto
  580. OKCONV = FAUX ;
  581. 'FINSI' ;
  582. 'FINSI' ;
  583. 'FINSI' ;
  584.  
  585.  
  586. * sorties textes
  587. * ==============
  588. *
  589.  
  590. 'SI' debug ;
  591. * nb d'espaces avant le nb d'itérations, pour faire joli
  592. * moins joli si &LINEAR > 9999.
  593. IL = &LINEAR - 1 ;
  594. TXT = '| ' ;
  595. 'SI' (IL < 1000) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ;
  596. 'SI' (IL < 100) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ;
  597. 'SI' (IL < 10) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ;
  598. * nb d'espaces avant le nb d'éléments désaturés, pour faire joli
  599. * moins joli si nbds > 9999.
  600. esp = ' ' ;
  601. 'SI' (nbds < 1000) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ;
  602. 'SI' (nbds < 100) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ;
  603. 'SI' (nbds < 10) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ;
  604.  
  605. * affichage des paramètres pour chaque itération
  606. * GBM change les maxi affichés - pas maxi TP1 mais maxi CHRG
  607. 'MESS' ('CHAINE' TXT IL ' |' MAXER ' | ' ('MINI' CHRG)
  608. ' | ' ('MAXI' CHRG) ' |' ('MAXI' CAPA)
  609. ' |' ('MINI' PERFR) ' | ' esp nbds ' | ') ;
  610. 'SI' OKCONV ;
  611. 'MESS' ('CHAINE' ' -------------------------------'
  612. '-------------------------------'
  613. '-------------------------------') ;
  614. 'FINSI' ;
  615. 'FINSI' ;
  616.  
  617. * tests sortie
  618. * ============
  619.  
  620.  
  621. *--- sortie si convergence
  622. 'SI' OKCONV ;
  623. * GBM - NON ON NE PEUT PAS SORTIR SI ON n'A PAS FAIT DE CALCUL
  624. * Le test précédent n'a pas de sens au premier indice.
  625. 'QUITTER' LINEAR ;
  626. 'FINSI' ;
  627.  
  628.  
  629. * petit message pour prévenir qu'on a changé la relaxation
  630. 'SI' (DEBUG 'ET' ('EGA' &LINEAR (itmaxi/2))) ;
  631. mess 'On accroît la sous-relaxation à 0.5' ;
  632. 'FINSI' ;
  633.  
  634.  
  635. *-- Sauvegarde des valeurs du pas de l'itération précédente
  636. TRH2N = TRH2 ;
  637. CAPAN = CAPAR ;
  638. PERFN = PERFR ;
  639. CHRGN = CHRG ;
  640. TENN = TEN ;
  641.  
  642.  
  643. *_____________________________________________________ non linéaire
  644.  
  645. 'FIN' LINEAR ;
  646. *__________________________________________________________________
  647.  
  648. 'SI' ('EGA' CONSER VRAI) ;
  649. TENNC = TENC ;
  650. 'FINSI' ;
  651.  
  652. * GBM - gerer les tabmodi à faux !!!!!!!!!!!!!!!!!!!!!!!!!
  653.  
  654. 'SI' ('NON' OKCONV) ;
  655. * message, adaptation du pas de tps ou penalisation
  656. DELTAT CFL0 = SATUTILS TESARRET &PENALDT NPENALDT LAST1 &LINEAR
  657. DELTAT maxer OKPENAL
  658. cofpenal CFL0 COFDIV debug
  659. SATUR tps ;
  660. 'SINON' ;
  661. * on a convergé
  662. * PENAL = FAUX ;
  663. 'QUITTER' PENALDT ;
  664. 'FINSI' ;
  665.  
  666.  
  667. *--------------------------------------------------------- artifice
  668.  
  669. 'FIN' PENALDT ;
  670.  
  671. *------------------------------------------------------------------
  672.  
  673.  
  674. FLU0 = FLU1 ;
  675. LAST1 = LAST1 + 1 ;
  676.  
  677. * Sauvegarde
  678. * ==========
  679.  
  680. ISOR ISAUV = SATUTILS SAVRESU SATUR TPSOR ISOR TPSANC TPS ISAUV DELTAT
  681. TRHANC TRH CHRGANC CHRG TENANC TEN PNSCANC PNSC SATANC SAT FLUANC
  682. FLU0 NIVEAU DSOR &LINEAR MDMH ;
  683.  
  684.  
  685. 'MESSAGE' ' ' ;
  686. 'MESSAGE' ' ' ;
  687. 'MESSAGE' 'RESIDU EN TEMPS '
  688. ((tfinal '-' tpsini) *
  689. ('MAXIMUM' (('RESULT' ('ABS' (TEN '-' TENANC))) '/'
  690. (deltat * ('RESULT' ('ABS' TEN)))))) ;
  691. 'MESSAGE' ' ' ;
  692. 'MESSAGE' ' ' ;
  693.  
  694.  
  695.  
  696. * Temps caractéristique et nouveau pas de temps :
  697. * ===============================================
  698.  
  699. *-- pas de temps idéal lié à la cinétique d'hydratation :
  700. * nb de face avec des flux pas trop faibles (> 1.D-9 x maximum)
  701.  
  702.  
  703.  
  704. *-- Nouveau pas de temps :
  705. 'SI' ('EXISTE' SATUR 'TEMPS_CALCULES') ;
  706. 'SI' (ICAL < DCAL) ;
  707. ICAL = ICAL + 1 ;
  708. DELTAT = ('EXTR' LTCALCUL ICAL) - TPS ;
  709. * GBM verifier que DTAUTO = FAUX dans ce cas sinon écrasé apres.
  710. 'FINSI' ;
  711. 'SINON' ;
  712. * quand on a donné SATUR.'DT_INITIAL', le pas de temps n'est
  713. * imposé qu'au départ. Après, il est automatique.
  714. DTAUTO = VRAI ;
  715. 'FINSI' ;
  716. 'SI' DTAUTO ;
  717. * pas de temps automatique
  718. * on modifie le CFL de façon à approcher nb d'itérations visé.
  719. 'MESSAGE' 'avant' CFL0;
  720. CFL0 = ((((Niter0 '/' 1.)) / &LINEAR) ** 0.5) ;
  721.  
  722. 'MESSAGE' 'apres' CFL0;
  723. SATUR.'CFL' = CFL0 ;
  724. 'FINSI' ;
  725.  
  726. VVOL = VHYB * PORO1 ;
  727. 'MESS' 'mini max capa' ('MINIMUM' (CAPA))
  728. ('MAXIMUM' (CAPA)) ;
  729. *'MESS' 'mini max perm' ('MINIMUM' (PERFR))
  730. * ('MAXIMUM' (PERFR)) ;
  731. DELTAT DTI = SATUTILS YNYTDT MDMH FLU0 MCHYB VVOL CFL0 debug
  732. DTAUTO DELTAT TERSCE ;
  733.  
  734. *-- sortie en cas de temps limite dépassé
  735. 'SI' ( TPS '>EG' TFINAL ) ;
  736. 'MENAGE' ;
  737. 'QUITTER' TRANSI ;
  738. 'FINSI' ;
  739.  
  740. * Préparation pas de temps suivant :
  741. * ==================================
  742.  
  743. * Sauvegarde des valeurs du pas de temps précédent
  744. TRHANC = TRH ;
  745. CHRGANC = CHRG ;
  746. FLUANC = FLU0 ;
  747. PNSANC = PNS ;
  748. PNSCANC = PNSC ;
  749. TENANC = TEN ;
  750. SATANC = SAT ;
  751. TPSANC = TPS ;
  752. PERFANC = PERFR ;
  753. * GBM REGARDE
  754. 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ;
  755. FLUIMPM1 = FLUIMPO ;
  756. 'FINSI' ;
  757. 'SI' ('EXISTE' SATUR 'FLUMIXTE') ;
  758. FLUMMPM1 = FLUMMPO ;
  759. 'FINSI' ;
  760. 'SI' ('EXISTE' SATUR 'SOURCE') ;
  761. TERSC2M1 = TERSC2 ;
  762. 'FINSI' ;
  763.  
  764. GEOL1 = 'TABLE' GEOLPF1 ;
  765. GEOL2 = 'TABLE' GEOLPF2 ;
  766.  
  767.  
  768.  
  769. *================================================= boucle transitoire
  770.  
  771. 'FIN' TRANSI ;
  772.  
  773. *====================================================================
  774. *
  775. 'FINPROC' ;
  776.  
  777.  
  778.  
  779.  
  780.  

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