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***************************************************
  • *
  • * Sauvegarde
  • * ==========
  •  
  • MESS 'SAUVEGARDE INITIALE---------------------------------------' ;
  •  
  •  
  •  
  • TMP TMP2 = SATUTILS SAVRESU SATUR TPSOR ISOR TPSANC TPS -1 DELTAT
  • TRHANC TRH CHRGANC CHRG TENANC TEN PNSCANC PNSC SATANC SAT FLUANC
  • FLU0 NIVEAU DSOR 0 MDMH ;
  •  
  • ***************FIN NPX***********************************************
  •  
  • *======================================================== transitoire
  •  
  • 'SI' ('EGA' NPAS0 0) ; NPAS0 = -1 ; 'FINSI' ;
  • 'REPETER' TRANSI NPAS0 ;
  •  
  • *====================================================================
  • *
  •  
  • IPAS = ICAL + &TRANSI - 1 ;
  •  
  •  
  • *-----------------------------------------------------------------
  •  
  • 'REPETER' PENALDT NPENALDT ;
  •  
  • CHRG = CHRGANC ;
  •  
  • TPS = TPSANC + DELTAT ;
  •  
  •  
  • *------------------------------------------------------------------
  • * Affichage information si debug vrai
  • 'MESSAGE' 'deltat dti' deltat dti;
  • SATUTILS AFFICH &PENALDT &TRANSI LAST1 TPS TPSANC DELTAT DTI
  • debug ;
  •  
  • ***************** INITIALISATION PAS DE TPS **************************
  • *
  • *- Incorporation des CLs
  • *
  • 'SI' ('EXISTE' SATUR 'TRACE_IMPOSE') ;
  • CHARIMPO = 'TIRE' SATUR . 'TRACE_IMPOSE' TPS ;
  • CHCLIM . 'DIRICHLET' = 'NOMC' NOMESPL CHARIMPO ;
  • 'FINSI' ;
  •  
  •  
  • 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ;
  • 'SI' (tpsanc '<EG' startflu) ;
  • FLUIMPO = 'TIRE' FLUIMP TPS ;
  • 'SINON' ;
  • FLUIMPO = 'COPIER' FLUIMPM1 ;
  • 'FINSI' ;
  • FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET' ;
  • 'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF')
  • 'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH')
  • 'ET' ('EGA' teta 1.D0))) ;
  • * On écrit les flux sous forme intégrale sauf en explicite et
  • * kranck-Nickholson pour EFMH. On les différentie ici
  • FLUMP = (FLUIMPO '-' FLUIMPM1) '/' DELTAT ;
  • 'SINON' ;
  • FLUMP = 'COPIER' FLUIMPO ;
  • 'FINSI' ;
  • CHCLIM . 'NEUMANN' = 'NOMC' NOMESPL FLUMP ;
  • 'FINSI' ;
  •  
  • 'SI' ('EXISTE' SATUR 'FLUMIXTE') ;
  • 'SI' (tpsanc '<EG' startmix) ;
  • FLUMMPO = 'TIRE' FLUMIX TPS ;
  • 'SINON' ;
  • FLUMMPO = 'COPIER' FLUMMPM1 ;
  • 'FINSI' ;
  • FLUMMPO = 'CHAN' 'ATTRIBUT' FLUMMPO 'NATURE' 'DISCRET' ;
  • 'SI' (('EGA' (SATUR . 'TYPDISCRETISATION') 'VF')
  • 'OU' (('EGA' (SATUR . 'TYPDISCRETISATION') 'EFMH')
  • 'ET' ('EGA' teta 1.D0))) ;
  • * On écrit les flux sous forme intégrale sauf en explicite et
  • * kranck-Nickholson pour EFMH. On les différentie ici
  • FLUMMP = ( FLUMMPO '-' FLUMMPM1) '/' DELTAT ;
  • 'SINON' ;
  • FLUMMP = 'COPIER' FLUMMPO ;
  • 'FINSI' ;
  • CHCLIM . 'FLUMIXTE' = TABLE ;
  • CHCLIM . 'FLUMIXTE' . 'VAL' = 'NOMC' NOMESPL FLUMMP ;
  • CHCLIM . 'FLUMIXTE' . 'COEFA' = SATUR . 'FLUMIXTE' . 'MIXCOFA' ;
  • * GBM rappel il y aura des modifs à faire car HYDRAU
  • CHCLIM . 'FLUMIXTE' . 'COEFB' = SATUR . 'FLUMIXTE' . 'MIXCOFB' ;
  • 'FINSI' ;
  •  
  • GEOL1 . 'CLIMITES' = CHCLIM ;
  •  
  • *
  • *- Calcul de la contribution des termes sources
  • *
  • TERSCE = 'MANU' 'CHPO' MAILCENT 1 'SOUR' 0. ;
  • TERSCE = 'NOMC' NOMESPL TERSCE ;
  •  
  • * Terme source propre :
  • 'SI' ( 'EXISTE' SATUR 'SOURCE' ) ;
  • 'SI' (tpsans <EG startsou) ;
  • * la source varie encore et est non nulle
  • TERSC2 = 'TIRE' TERSOU TPS ;
  • 'SINON' ;
  • * la source est nulle on garde l'intégrale précédente
  • TERSC2 = TERSC2M1 ;
  • 'FINSI' ;
  • TERSC2 = 'NOMC' NOMESPL TERSC2 ;
  • TERSCV = (TERSC2 '-' TERSC2M1) '/' DELTAT ;
  • TERSCE = TERSCV + TERSCE ;
  • 'FINSI' ;
  •  
  • * Le terme source physique est tersce.
  •  
  • GEOL1 . 'DELTAT' = DELTAT ;
  • GEOL1 . 'METHINV' = SATUR . 'METHINV' ;
  •  
  • *----------------------------------------------------- non linéaire
  •  
  • 'REPETER' LINEAR itmaxi ;
  • 'SI' ('EGA' CONSER VRAI) ;
  • GEOL1 . 'CONCENTRATION' = CHRG ;
  • 'FINSI' ;
  •  
  • *__________________________________________________________________
  •  
  • * Préparation relaxation - GBM tester efficacité
  • * ======================
  • 'SI' ('EGA' &LINEAR 1) ;
  • * pas de relaxation à la première itération
  • cofrelax = 1.D0 ;
  • 'SINON' ;
  • 'SI' (&LINEAR '<EG' (itmaxi/2)) ;
  • * on relaxe à partir du deuxième pas de temps
  • cofrelax = SATUR.'SOUS_RELAXATION' ;
  • 'SINON' ;
  • * on relaxe plus si on dépasse itmaxi/2
  • * GBM GBM GBMGBM 0.5 normalement
  • cofrelax = 0.5D0 ;
  • 'FINSI' ;
  • 'FINSI' ;
  •  
  • * Calcul des nouveaux paramètres
  • * ==============================
  •  
  • *-- Pression tronquée à utiliser dans les lois de comportement:
  • * P1 pression pascal, PNS pression tronquée
  • PNS = SATUTILS TRONKP MDMH CHRG TRH L_GRAV NIVEAU RHOWG CONVH ;
  • PNSC = SATUTILS TRONKP MDMH CHRG TRH L_GRAV 'CENTRE' RHOWG CONVH ;
  •  
  •  
  • * GBM reprendre calcul sat, ten et capa
  • * mettre un flag pour calculer CAPA ou pas
  • * retirer PNSANC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  • * TENANC
  •  
  • 'SI' ('EGA' CONSER VRAI) ;
  • nbds SATC TENC CAPA = SATUTILS KALSOT MDMH LOIS 'CENTRE' PNSC ;
  • dum1 SAT TEN CAPAD = SATUTILS KALSOT MDMH LOIS NIVEAU PNS ;
  • 'SINON' ;
  • nbds SATC TEN CAPA = SATUTILS KALSOT MDMH LOIS 'CENTRE' PNSC ;
  • dum1 SAT TEN CAPA = SATUTILS KALSOT MDMH LOIS NIVEAU PNS ;
  • 'FINSI' ;
  •  
  •  
  •  
  • *-- Calcul perméabilité
  • * PERF = 0. ;
  • 'REPETER' BB ('DIME' LOIP.'INDEX') ;
  • NOMT = 'CHAINE' LOIP.'INDEX'.&BB ;
  • LOI1 = LOIP . NOMT ;
  • MOD1 = LOI1.'MODELE' ;
  • NOMKR = 'MOT' ('TEXTE' ('CHAINE' LOI1.'NOM_PROCEDURE')) ;
  •  
  • 'SI' (('EGA' LOI1.'SOUSTYPE' ('CHAINE' PUISSANCE) ) 'OU'
  • ('EGA' LOI1.'SOUSTYPE' ('CHAINE' PERSONNELLE)) 'OU'
  • ('EGA' LOI1.'SOUSTYPE' ('CHAINE' MUALEM)) 'OU'
  • ('EGA' LOI1.'SOUSTYPE' ('CHAINE' BURDINE)) 'OU'
  • ('EGA' LOI1.'SOUSTYPE' ('CHAINE' BROOKS_COREY)) 'OU'
  • ('EGA' LOI1.'SOUSTYPE' ('CHAINE' MUALEM_BURDINE))) ;
  • * réduction de la saturation au sous-domaine (local)
  • SATL = 'REDU' SAT ('DOMA' MOD1 NIVEAU) ;
  • PERFL = ('TEXTE' NOMKR) LOI1 SATL ;
  • PERFL = RECENTRE PERFL MOD1 NIVEAU ;
  • 'SINON' ;
  • PNL = 'REDU' PNS ('DOMA' MOD1 NIVEAU) ;
  • PERFL = ('TEXTE' NOMKR) LOI1 PNL ;
  • PERFL = RECENTRE PERFL MOD1 NIVEAU ;
  • 'FINSI' ;
  •  
  • * 'SI' ('EGA' NIVEAU 'FACE') ;
  • * LOI1.'ABMC' = 'ABS' ('DOMA' MOD1 'ORIENTAT') ;
  • * 'FINSI' ;
  •  
  •  
  • * concaténation perméabilité
  • LCOMP = 'EXTR' PERFL 'COMP' ;
  • NCOMP = 'DIME' LCOMP ;
  • 'REPETER' CONSPER NCOMP ;
  • J = &CONSPER ;
  • COMPI = 'EXTR' LCOMP J ;
  • PERFINI = 'MANU' 'CHPO' ('DOMA' MDMH 'CENTRE') 1 COMPI 0.
  • 'NATURE' 'DISCRET' ;
  • PERFI = 'EXCO' PERFL COMPI COMPI ;
  • PERFI = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' COMPI PERFINI PERFI;
  • *
  • 'SI' (&BB 'EGA' 1) ;
  • 'SI' ('EGA' J 1) ;
  • PERFT = PERFI ;
  • 'SINON' ;
  • PERFT = PERFT 'ET' PERFI ;
  • 'FINSI' ;
  • 'SINON' ;
  • 'SI' ('EGA' J 1) ;
  • PERFT = PERFT 'ET' PERFI ;
  • 'SINON' ;
  • PERFT = PERFT 'ET' PERFI ;
  • 'FINSI' ;
  • 'FINSI' ;
  • *
  • 'FIN' CONSPER ;
  • *
  • 'OUBLIER' satl ;
  • 'OUBLIER' perl;
  •  
  • 'FIN' BB ;
  •  
  • *-- coefficient d'emmagasinement si milieu saturé
  • * GBM ATTENTION ICI BIZARRE
  • * GBM MET A 0 CAR DISCONTINUITE DU TERME DEVANT DT, 0 AVANT
  • * SAT ET POSITIF APRES. DOIT ETRE MIS DANS LES LOIS DE SAT
  • * ET EVOLUER CONTINUEMENT.
  • EMMAG = ('MASQUE' ('NOMC' SATC 'SCAL') 'EGSUP' 0.99 ) * COFEMMAG;
  •  
  • * EMMAG = NOMC SCAL (SATC * COFEMMAG) ;
  • EMMAG = 'KCHT' MDMH 'SCAL' 'CENTRE' 'COMP' 'SCAL' EMMAG ;
  • * 'MESSAGE' 'maxi min emmag' ('MINIMUM' emmag) ('MAXIMUM' emmag);
  •  
  •  
  • * Construction système matriciel
  • * ==============================
  •  
  • *-- Penalisation -
  • * GBM ERREUR sur NOMC 'SOUR'
  • 'SI' PENAL ;
  • TERSC1 = TERSCE + ('NOMC' NOMESPL
  • ((CHRG - CHRGANC) * VHYB / DELTAT * cofpenal)) ;
  • 'SINON' ;
  • TERSC1 = TERSCE ;
  • 'FINSI' ;
  •  
  •  
  • 'SI' ('EGA' CONSER VRAI) ;
  • TERSC1 = TERSC1 - ('NOMC' NOMESPL
  • ((TENC - TENNC) * VHYB / DELTAT )) ;
  • 'FINSI' ;
  • GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSC1 ;
  •  
  • * relaxation
  •  
  •  
  • * 'MESSAGE' 'coefrelaaaaaaaaaaaaaaaa' cofrelax;
  • 'SI' ('NEG' cofrelax 1.D0 1.D-14) ;
  • * on relaxe CAPA et PERF
  • CAPAR = (CAPA '*' cofrelax)
  • '+' (CAPAN '*' (1. - cofrelax)) ;
  • PERFR = ( PERFT '*' cofrelax)
  • '+' (PERFN '*' (1. - cofrelax)) ;
  • 'SINON' ;
  • CAPAR = CAPA ;
  • PERFR = PERFT ;
  • 'FINSI' ;
  •  
  •  
  •  
  •  
  •  
  • *-- coef devant la derivée temps
  • COFDT = CAPAR '+' EMMAG ;
  •  
  • * permeabilité
  • * LAM = ('NOMC' PERFR 'K' NATURE DIFFUS) ;
  •  
  • * GBM - penalisation un peu particuliere ! ecrire ce que cela
  • * veut dire
  • 'SI' PENAL ;
  • COFDT = COFDT '+' COFPENAL ;
  • 'FINSI' ;
  •  
  • * initialisation TRANSGEOL - gbm bien verif les tabmodi
  • GEOL1 . 'POROSITE' = COFDT ;
  • TABMODI . 'POROSITE' = VRAI ;
  • GEOL1 . 'DIFFUSIVITE' = PERFR ;
  • TABMODI . 'DIFFUSIV' = VRAI ;
  •  
  • * Résolution :
  • * ============
  •  
  • * GBM appel trangeol - tester peut etre aussi premier penal
  • 'SI' ( (&TRANSI 'EGA' 1)) ;
  • GEOLPF1 GEOLPF2 = TRANGEOL MDMH GEOL1 ;
  • 'SINON' ;
  • GEOLPF1 GEOLPF2 = TRANGEOL MDMH GEOL1 GEOL2 ;
  • 'FINSI' ;
  • CHRG = GEOLPF1 . 'CONCENTRATION' ;
  • FLU1 = GEOLPF1 . 'FLUXDIFF' ;
  • 'SI' ('EGA' SATUR . 'TYPDISCRETISATION' EFMH) ;
  • TRH2 = GEOLPF2 . 'TRACE_CONC' ;
  • 'SINON' ;
  • * GBM §§§§§§§§§§§§§§§§§§§§§§§
  • TRH2 = 0.D0 * FLU1;
  • 'FINSI' ;
  •  
  • * GBM - a clarifier
  • TRH = TRH2 ;
  •  
  • * FIN DE BOUCLE :
  • * =================
  •  
  •  
  •  
  • * test de convergence sur le résidu
  • * =================================
  • * A-t-on convergé à l'itération précédente (avec les nouveaux paramètres) ?
  • * flux correspondant au résidu système résolu sans relaxation
  •  
  • **
  • ** GBM reflechir comment améliorer la conservation
  •  
  • 'SI' (&LINEAR 'EGA' 1) ;
  • OKCONV = FAUX ;
  • maxer = 1.D0;
  • 'SINON' ;
  •  
  • RES1 = 'MAXIMUM' ('RESULT' ('ABS' (TEN - TENN))) ;
  • RES2 =
  • ((0.1D0 * ('MAXIMUM' ('RESULT' ('ABS' (TEN)))))
  • + ( ('MAXIMUM' ('RESULT' ('ABS' (TEN '-' TENANC))))));
  •  
  •  
  •  
  •  
  • RES3 = 'MAXIMUM' ('RESULT' ('ABS' (PERFR - PERFN))) ;
  • RES4 = ('MAXIMUM' ('RESULT' ('ABS' ((perfr '-' perfn)))))
  • '/' ('MAXIMUM' ('RESULT' PERFR)) ;
  • RES3 = RES3 '/' (
  • 'MAXIMUM' ('RESULT' ('ABS' (PERFR - PERFANC))) + 1.D-100) ;
  •  
  •  
  •  
  •  
  • * 'MESSAGE' 'minmax res1 res2' res1 res2 ('MAXIMUM'
  • * (TEN - TENANC)) ;
  • MAXER = RES1 '/' RES2 ;
  • MAXER = MAXER '+' RES3 ;
  •  
  • * on prend en compte la taille du pas de tps comparé à CFL
  • 'SI' (&TRANSI 'EGA' 1) ;
  • * Pour le premier pas de temps les flux ne sont pas connus
  • * donc le pas de temps CFL est mal évalué. On ne l'inclue
  • * pas dans le crtiere de convergence
  • OKCONV = MAXER < (ERR0) ;
  • 'SINON' ;
  • OKCONV = MAXER < (ERR0 * DTI '/' DELTAT) ;
  • * 'MESSAGE' 'error4 err3' res4 res3 ;
  • OKCONV = OKCONV 'ET' (RES4 < ( ERR0)) ;
  • 'FINSI' ;
  •  
  • 'SI' (DTAUTO) ;
  • 'SI' (&LINEAR 'EGA' 2) ;
  • * on impose de boucler au moins 2 fois en pas de tps auto
  • OKCONV = FAUX ;
  • 'FINSI' ;
  • 'FINSI' ;
  • 'FINSI' ;
  •  
  •  
  • * sorties textes
  • * ==============
  • *
  •  
  • 'SI' debug ;
  • * nb d'espaces avant le nb d'itérations, pour faire joli
  • * moins joli si &LINEAR > 9999.
  • IL = &LINEAR - 1 ;
  • TXT = '| ' ;
  • 'SI' (IL < 1000) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ;
  • 'SI' (IL < 100) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ;
  • 'SI' (IL < 10) ; txt = 'CHAINE' txt ' ' ; 'FINSI' ;
  • * nb d'espaces avant le nb d'éléments désaturés, pour faire joli
  • * moins joli si nbds > 9999.
  • esp = ' ' ;
  • 'SI' (nbds < 1000) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ;
  • 'SI' (nbds < 100) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ;
  • 'SI' (nbds < 10) ; esp = 'CHAINE' esp ' ' ; 'FINSI' ;
  •  
  • * affichage des paramètres pour chaque itération
  • * GBM change les maxi affichés - pas maxi TP1 mais maxi CHRG
  • 'MESS' ('CHAINE' TXT IL ' |' MAXER ' | ' ('MINI' CHRG)
  • ' | ' ('MAXI' CHRG) ' |' ('MAXI' CAPA)
  • ' |' ('MINI' PERFR) ' | ' esp nbds ' | ') ;
  • 'SI' OKCONV ;
  • 'MESS' ('CHAINE' ' -------------------------------'
  • '-------------------------------'
  • '-------------------------------') ;
  • 'FINSI' ;
  • 'FINSI' ;
  •  
  • * tests sortie
  • * ============
  •  
  •  
  • *--- sortie si convergence
  • 'SI' OKCONV ;
  • * GBM - NON ON NE PEUT PAS SORTIR SI ON n'A PAS FAIT DE CALCUL
  • * Le test précédent n'a pas de sens au premier indice.
  • 'QUITTER' LINEAR ;
  • 'FINSI' ;
  •  
  •  
  • * petit message pour prévenir qu'on a changé la relaxation
  • 'SI' (DEBUG 'ET' ('EGA' &LINEAR (itmaxi/2))) ;
  • mess 'On accroît la sous-relaxation à 0.5' ;
  • 'FINSI' ;
  •  
  •  
  • *-- Sauvegarde des valeurs du pas de l'itération précédente
  • TRH2N = TRH2 ;
  • CAPAN = CAPAR ;
  • PERFN = PERFR ;
  • CHRGN = CHRG ;
  • TENN = TEN ;
  •  
  •  
  • *_____________________________________________________ non linéaire
  •  
  • 'FIN' LINEAR ;
  • *__________________________________________________________________
  •  
  • 'SI' ('EGA' CONSER VRAI) ;
  • TENNC = TENC ;
  • 'FINSI' ;
  •  
  • * GBM - gerer les tabmodi à faux !!!!!!!!!!!!!!!!!!!!!!!!!
  •  
  • 'SI' ('NON' OKCONV) ;
  • * message, adaptation du pas de tps ou penalisation
  • DELTAT CFL0 = SATUTILS TESARRET &PENALDT NPENALDT LAST1 &LINEAR
  • DELTAT maxer OKPENAL
  • cofpenal CFL0 COFDIV debug
  • SATUR tps ;
  • 'SINON' ;
  • * on a convergé
  • * PENAL = FAUX ;
  • 'QUITTER' PENALDT ;
  • 'FINSI' ;
  •  
  •  
  • *--------------------------------------------------------- artifice
  •  
  • 'FIN' PENALDT ;
  •  
  • *------------------------------------------------------------------
  •  
  •  
  • FLU0 = FLU1 ;
  • LAST1 = LAST1 + 1 ;
  •  
  • * Sauvegarde
  • * ==========
  •  
  • ISOR ISAUV = SATUTILS SAVRESU SATUR TPSOR ISOR TPSANC TPS ISAUV DELTAT
  • TRHANC TRH CHRGANC CHRG TENANC TEN PNSCANC PNSC SATANC SAT FLUANC
  • FLU0 NIVEAU DSOR &LINEAR MDMH ;
  •  
  •  
  • 'MESSAGE' ' ' ;
  • 'MESSAGE' ' ' ;
  • 'MESSAGE' 'RESIDU EN TEMPS '
  • ((tfinal '-' tpsini) *
  • ('MAXIMUM' (('RESULT' ('ABS' (TEN '-' TENANC))) '/'
  • (deltat * ('RESULT' ('ABS' TEN)))))) ;
  • 'MESSAGE' ' ' ;
  • 'MESSAGE' ' ' ;
  •  
  •  
  •  
  • * Temps caractéristique et nouveau pas de temps :
  • * ===============================================
  •  
  • *-- pas de temps idéal lié à la cinétique d'hydratation :
  • * nb de face avec des flux pas trop faibles (> 1.D-9 x maximum)
  •  
  •  
  •  
  • *-- Nouveau pas de temps :
  • 'SI' ('EXISTE' SATUR 'TEMPS_CALCULES') ;
  • 'SI' (ICAL < DCAL) ;
  • ICAL = ICAL + 1 ;
  • DELTAT = ('EXTR' LTCALCUL ICAL) - TPS ;
  • * GBM verifier que DTAUTO = FAUX dans ce cas sinon écrasé apres.
  • 'FINSI' ;
  • 'SINON' ;
  • * quand on a donné SATUR.'DT_INITIAL', le pas de temps n'est
  • * imposé qu'au départ. Après, il est automatique.
  • DTAUTO = VRAI ;
  • 'FINSI' ;
  • 'SI' DTAUTO ;
  • * pas de temps automatique
  • * on modifie le CFL de façon à approcher nb d'itérations visé.
  • 'MESSAGE' 'avant' CFL0;
  • CFL0 = ((((Niter0 '/' 1.)) / &LINEAR) ** 0.5) ;
  •  
  • 'MESSAGE' 'apres' CFL0;
  • SATUR.'CFL' = CFL0 ;
  • 'FINSI' ;
  •  
  • VVOL = VHYB * PORO1 ;
  • 'MESS' 'mini max capa' ('MINIMUM' (CAPA))
  • ('MAXIMUM' (CAPA)) ;
  • *'MESS' 'mini max perm' ('MINIMUM' (PERFR))
  • * ('MAXIMUM' (PERFR)) ;
  • DELTAT DTI = SATUTILS YNYTDT MDMH FLU0 MCHYB VVOL CFL0 debug
  • DTAUTO DELTAT TERSCE ;
  •  
  • *-- sortie en cas de temps limite dépassé
  • 'SI' ( TPS '>EG' TFINAL ) ;
  • 'MENAGE' ;
  • 'QUITTER' TRANSI ;
  • 'FINSI' ;
  •  
  • * Préparation pas de temps suivant :
  • * ==================================
  •  
  • * Sauvegarde des valeurs du pas de temps précédent
  • TRHANC = TRH ;
  • CHRGANC = CHRG ;
  • FLUANC = FLU0 ;
  • PNSANC = PNS ;
  • PNSCANC = PNSC ;
  • TENANC = TEN ;
  • SATANC = SAT ;
  • TPSANC = TPS ;
  • PERFANC = PERFR ;
  • * GBM REGARDE
  • 'SI' ('EXISTE' SATUR 'FLUX_IMPOSE') ;
  • FLUIMPM1 = FLUIMPO ;
  • 'FINSI' ;
  • 'SI' ('EXISTE' SATUR 'FLUMIXTE') ;
  • FLUMMPM1 = FLUMMPO ;
  • 'FINSI' ;
  • 'SI' ('EXISTE' SATUR 'SOURCE') ;
  • TERSC2M1 = TERSC2 ;
  • 'FINSI' ;
  •  
  • GEOL1 = 'TABLE' GEOLPF1 ;
  • GEOL2 = 'TABLE' GEOLPF2 ;
  •  
  •  
  •  
  • *================================================= boucle transitoire
  •  
  • 'FIN' TRANSI ;
  •  
  • *====================================================================
  • *
  • 'FINPROC' ;
  •  
  •  
  •  
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales