Télécharger la plaquette Cast3M


fr en

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.  

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