Télécharger la plaquette Cast3M


fr en

Télécharger dvispw.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  dvispw.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. * CAVITE CARREE ­ VAHL­DAVIS Methode de projection implicite *
  5. * Malvina Renesson Aout 1999 * 
  6. * Version utilisant l'option EQUA de EQEX
  7. * le terme de flottabilite est traite dans la procedure CALCUL
  8. * La boucle en temps est exterieure a EXEC
  9. *************************************************************** 
  10. *$$$$ EXEC
  11. * EXEC      PROCEDUR  MAGN      03/03/31    21:15:04     4631           
  12. *X EXEC (Procedure)
  13. *   Procedure EXEC
  14. *
  15. *   Objet   :  Execute un algorithme décrit dans une table RV
  16. *              de type EQEX.
  17. *              Cette table est créée par l'opérateur EQEX
  18. *
  19. *   Syntaxe :  EXEC RV ;
  20. *
  21. ***********************************************************************
  22. * VERSION    : ????
  23. * HISTORIQUE : 20/12/99: gounand
  24. * Rajout de la gestion de la matrice servant à l'assemblage
  25. * (rv . 'METHINV' . 'MATASS') 
  26. *
  27. * HISTORIQUE : 10/11/00: magnaud
  28. * Reecriture de exec pour tirer parti du // (operateur ASSI et RESO)
  29. * Pour l'instant cette procedure est testee experimentalement dans
  30. * la directive EQUA de EXEC
  31. *
  32. * HISTORIQUE : 04/04/01: magnaud
  33. * Calcul et impression tout les FIDT cycles de la norme Linf du residu
  34. *
  35. * HISTORIQUE :
  36. * HISTORIQUE :
  37. ************************************************************************
  38. * DISCPRES = (d'apres KMIC)
  39. * KPRE=3 pression P0   KPRE=4 pression P1  KPRE=2 cas macro 1ère génération
  40. * KPRE=5 pression MSOMMET
  41. ************************************************************************
  42. 'DEBPROC'  EXEC ;
  43. 'ARGUMENT' rv*'TABLE   ' ;
  44.  
  45. * Logique pilotant le recalcul de la matrice de pression
  46.  CALPRE=FAUX                                                   ;
  47.  Si ('EXIST' rv 'CALPRE')                                      ;
  48.   vertytab rv 'CALPRE' 'LOGIQUE'                               ;
  49.   CALPRE=rv.'CALPRE'                                           ;
  50.  Finsi                                                         ;
  51.  
  52.  Si(NON('EXIST' rv 'XEQUA'));
  53.  rv.'XEQUA'= FAUX ;
  54.  Finsi ;
  55.  
  56. ********************************************************************************
  57. ********************************************************************************
  58. ****** Directive EQUA **********************************************************
  59. ********************************************************************************
  60. ********************************************************************************
  61.  
  62.  Si (rv.'XEQUA');
  63.  
  64.     IMPKRES=0 ;
  65.   Si ('EXIST' rv 'NUMASSI'); numassi=rv.'NUMASSI';
  66.   Sinon;                     numassi=1     ;
  67.   Finsi;
  68.  
  69.  
  70.   Xdfdt  = FAUX ;
  71.  
  72.   ITMA=rv . 'ITMA';
  73.  'SI' ('<EG' ITMA 1) ;
  74.   ITMA=1 ;
  75.  'FINSI'  ;
  76.  
  77.   NITER=rv . 'NITER';
  78.  'SI' ('<EG' NITER 1) ;
  79.   NITER=1 ;
  80.  'FINSI'  ;
  81.  
  82.   OMEGA=rv . 'OMEGA';
  83.  'SI' ('<EG' NITER 1) ;
  84.   OMEGA=1.;
  85.  'FINSI'  ;
  86.  
  87.  NUPDT=rv.'PASDETPS'.'NUPASDT';
  88.  
  89.  'SI'((NON (rv.'TRAN')) et ('>' ITMA 1) );
  90.  Mess ' Incoherence : Algorithme Permanent et ITMA >1';
  91.  QUITTER EXEC ;
  92.  'FINSI' ;
  93.  
  94.  'SI'((rv.'TRAN') et (NON (rv.'PROJ')) et ('EGA' NUPDT 1));
  95.  Mess ' Algorithme transitoire standard implicite ou explicite ';
  96.  mess '========================================================';
  97.  'FINSI' ;
  98.  
  99.  'SI'((rv.'PROJ') et ('EGA' NUPDT 1));
  100.   mess 'Algorithme de Projection';
  101.   mess '========================';
  102.  'FINSI' ;
  103.  
  104.   'SI' ( non (rv.'TRAN'));
  105.   mess 'Algorithme de recherche d un permanent ';
  106.   mess '=======================================';
  107.  'FINSI' ;
  108.  
  109.  
  110. 'REPETER' BCLTPS ITMA ;
  111.   si (rv.'PROJ');
  112.    testp1 = non ( exist rv 'MATC') ;
  113.   sinon ;
  114.    testp1 = FAUX  ;
  115.   finsi ;
  116.  'REPETER' BCLINT NITER ;
  117.  Si (non (rv.'TRAN'));
  118.    mess 'ITERATION ' &BCLINT;
  119.  Finsi ;
  120.  
  121.   Si (rv.'XRIG');
  122.    MRIG='RGDT    ';
  123.    sf mau  = assi numassi 'KOPS' 'RIGIDITE' ;
  124.    sp map  = assi numassi 'KOPS' 'RIGIDITE' ;
  125.    sr mar  = assi numassi 'KOPS' 'RIGIDITE' ;
  126.   Sinon ;
  127.    MRIG=TEXT '    ';
  128.    sf mau  = assi numassi 'KOPS' 'MATRIK' ;
  129.    sp map  = assi numassi 'KOPS' 'MATRIK' ;
  130.    sr mar  = assi numassi 'KOPS' 'MATRIK' ;
  131.   Finsi ;
  132.  
  133.   NRES=0    ;
  134.  
  135. *****************************************************
  136. *  ._._.                  \|/
  137. *  ._._.   CALCUL MATRICES ELEMENTAIRES
  138. *  ._._.                  /|\
  139. *****************************************************
  140.  
  141.  
  142.    'REPETER' bloc2 ('DIME' (rv . 'LISTOPER')) ;
  143.       nomper = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ;
  144. *     mess ' Début BLOC2 nomper=' nomper ;
  145.       notable= 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ;
  146.        Si('>EG' (rv.'IMPR') 2); mess 'Operateur ' nomper ; Finsi ;
  147.       td=rv.notable ;
  148.       tdk=td.'KOPT' ;
  149.       iarg= rv.notable.'IARG' ;
  150.    arg1 = TEXT '   ' ; arg2 = TEXT '   ' ; arg3 = TEXT '   ' ;
  151.    arg4 = TEXT '   ' ; arg5 = TEXT '   ' ; arg6 = TEXT '   ' ;
  152.    arg7 = TEXT '   ' ; arg8 = TEXT '   ' ; arg9 = TEXT '   ' ;
  153.  
  154.    Si ('>EG' iarg 1 );
  155.       Si ( EGA ('TYPE' (td.'ARG1')) 'MOT     ' ) ;
  156.       arg1=rv.'INCO'.(td.'ARG1') ;
  157.       Sinon ;
  158.       arg1=(td.'ARG1') ;
  159.       Finsi ;
  160.    Finsi ;
  161.  
  162.    Si ('>EG' iarg 2 );
  163.       Si ( EGA ('TYPE' (td.'ARG2')) 'MOT     ' ) ;
  164.       arg2=rv.'INCO'.(td.'ARG2') ;
  165.       Sinon ;
  166.       arg2=(td.'ARG2') ;
  167.       Finsi ;
  168.    Finsi ;
  169.  
  170.    Si ('>EG' iarg 3 );
  171.       Si ( EGA ('TYPE' (td.'ARG3')) 'MOT     ' ) ;
  172.       arg3=rv.'INCO'.(td.'ARG3') ;
  173.       Sinon ;
  174.       arg3=(td.'ARG3') ;
  175.       Finsi ;
  176.    Finsi ;
  177.  
  178.    Si ('>EG' iarg 4 );
  179.       Si ( EGA ('TYPE' (td.'ARG4')) 'MOT     ' ) ;
  180.       arg4=rv.'INCO'.(td.'ARG4') ;
  181.       Sinon ;
  182.       arg4=(td.'ARG4') ;
  183.       Finsi ;
  184.    Finsi ;
  185.  
  186. * Traitement du cas transitoire
  187.    XTRA =CHAI 'PERM' ; DT   = TEXT '   ' ;  SCHE=TEXT '   ' ;
  188.    inc1 = TEXT '   ' ; inc2 = TEXT '   ' ;
  189.    imc1 = TEXT '   ' ; imc2 = TEXT '   ' ;
  190.    BETAT= TEXT '   ' ;
  191.    XBDF2 = FAUX ;
  192.  
  193.    XOT=(('EGA' nomper 'NS      ') ou ('EGA' nomper 'TSCA    ')
  194.      ou ('EGA' nomper 'DFDT    ') ou ('EGA' nomper 'TSCAL   '));
  195.  
  196.    Si (rv.'PROJ' et ('EGA' nomper 'KBBT    ')) ;
  197.     XTRA=CHAI 'PROJ' ;
  198.    Finsi ;
  199.  
  200.    Si (rv.'TRAN') ;
  201.      Xdfdt  = VRAI ;
  202.      XTRA=CHAI 'TRANS' ;
  203.      Si (rv.'PROJ') ; XTRA=CHAI 'PROJ' ; Finsi ;
  204.      DT=(rv.'DELTAT') ;
  205.      Si ( EGA ('TYPE' DT) 'MOT     ' ) ;
  206.        DT=rv.'INCO'.DT ;
  207.      Finsi ;
  208.      SCHE = rv.'SCHEMAT';
  209.      Si (EGA SCHE 'SEMI');
  210.        BETAT=(rv.'Betat') ;
  211.      Finsi ;
  212.      Si (EGA SCHE 'BDF2');
  213.        XBDF2 = VRAI ;
  214.      Finsi ;
  215.    Finsi ;
  216.  
  217.      nbic=dime td.'LISTINCO';
  218.  
  219.      Si ('>EG' nbic 1 );
  220.        ninc1=extr td.'NUMEINCO' 1 ;
  221.        a= chai 'FORMAT' '(A3,I1)' 'INC' ninc1 ;
  222.        inc1=rv.'INCO'.(rv.a);
  223.          Si (XBDF2) ;
  224.           a= chai 'FORMAT' '(A3,I1)' 'IMC' ninc1 ;
  225.           imc1=rv.'INCO'.(rv.a);
  226.          Finsi;
  227.      Finsi ;
  228.  
  229.      Si ('>EG' nbic 2 );
  230.        ninc2=extr td.'NUMEINCO' 2 ;
  231.        a= chai 'FORMAT' '(A3,I1)' 'INC' ninc2 ;
  232.        inc2=rv.'INCO'.(rv.a);
  233.          Si (XBDF2) ;
  234.           a= chai 'FORMAT' '(A3,I1)' 'IMC' ninc2 ;
  235.           imc2=rv.'INCO'.(rv.a);
  236.          Finsi;
  237.      Finsi ;
  238.  
  239.  
  240.    MM=TEXT '    ';
  241.    DCEN='????' ;
  242.    Si(EGA (tdk.'MTRMASS') 2 );MM='MMDIAGO';Finsi;
  243.    Si(EGA (tdk.'IDCEN') 1 );DCEN='CENTREE';Finsi;
  244.    Si(EGA (tdk.'IDCEN') 2 );DCEN='SUPGDC ';Finsi;
  245.    Si(EGA (tdk.'IDCEN') 3 );DCEN='SUPG   ';Finsi;
  246.    Si(EGA (tdk.'IKOMP') 0 );MC=TEXT '    ';Finsi;
  247.    Si(EGA (tdk.'IKOMP') 1 );MC='CONS';Finsi;
  248.  
  249.    INEFMD= tdk.'INEFMD' ;
  250.    tdz=TEXT '   '  ;
  251.    Si (EGA nomper 'VNIMP   ');
  252.    tdz=CHAI 'V' ;
  253.    nomper= text 'NS';
  254.    Finsi ;
  255.  
  256. * temp place ;
  257. * mess ' Avant ' nomper  ' tdz ' tdz ' DT=' DT ' XTRA=' ; list xtra ;
  258. * mess SCHE ; mess BETAT;
  259.    msi mai =assi  numassi  nomper tdz
  260.      (td.'DOMZ') (td.'LISTINCO') (td.'TYPEINCO') (rv.'TYPPRESS')
  261.      INEFMD MRIG MM 'CONS' 'STABP' (tdk.'STABP') 'CMD' (tdk.'CMD') DCEN
  262.      XTRA DT SCHE BETAT
  263.     'INCO' inc1 imc1 inc2 imc2
  264.    (td.'IARG') arg1 arg2 arg3 arg4 ;
  265. * temp place ;
  266. * mess ' Apres ' nomper;
  267.  
  268.        Si ((rv.'PROJ') et ('EGA' nomper 'KBBT    '));
  269.         msp =assi  numassi  msi ;
  270.         mac =assi  numassi  mai ;
  271.         'ITER' bloc2 ;
  272.        Finsi ;
  273.  
  274.        Si ((rv.'PROJ') et ('EGA' nomper 'DFDT    '));
  275.         mam =assi  numassi  mai ;
  276.        Finsi ;
  277.  
  278. *      Si ((rv.'PROJ') et ('EGA' nomper 'TOIMP   '));
  279. *       list msi ;
  280. *      Finsi ;
  281.  
  282.        Si ((rv.'PROJ') et ('EGA' nomper 'NS      '));
  283.      XTRA=CHAI 'TRANS' ;
  284.      DT=(rv.'DELTAT') ;
  285.      Si ( EGA ('TYPE' DT) 'MOT     ' ) ;
  286.        DT=rv.'INCO'.DT ;
  287.      Finsi ;
  288.      SCHE = rv.'SCHEMAT';
  289.      Si (EGA SCHE 'SEMI');
  290.        BETAT=(rv.'Betat') ;
  291.      Finsi ;
  292.      Si (EGA SCHE 'BDF2');
  293.        XBDF2 = VRAI ;
  294.      Finsi ;
  295. *       mess ' Appel DFDT DT=' DT ' XTRA=' ; list xtra ;
  296.         msm mam =assi  numassi  DFDT
  297.         (td.'DOMZ') (td.'LISTINCO') (td.'TYPEINCO') (rv.'TYPPRESS')
  298.         INEFMD MRIG MM 'CONS' 'CMD' (tdk.'CMD') DCEN
  299.         XTRA DT SCHE BETAT
  300.        'INCO' inc1 imc1 inc2 imc2
  301.         1 arg1 ;
  302.         oubli msm;
  303.        Finsi ;
  304.  
  305.     mau = assi numassi mau 'ET' mai ;
  306.     sf  = assi numassi sf  'ET' msi ;
  307.  
  308.  
  309.      rv.'TVNP'=FAUX;
  310.      rv.'TVNPC'=FAUX;
  311.  
  312.    Si ('EGA' nomper  'VNIMP   ');                          *
  313.    mess ' Cas oper VNIMP' ;
  314.     Si (ega (rv.'TYPPRESS') 'MSOMMET');
  315.      rv.'TVNPC'=VRAI;
  316.      mvnpc=mai ;
  317.      svnpc=msi ;
  318.     Sinon ;
  319.      rv.'TVNP'=VRAI ;
  320.      mar = assi numassi mar 'ET' mai ;
  321.      sr  = assi numassi sr  'ET' msi ;
  322.     Finsi ;
  323.    Sinon ;
  324.  
  325.    Finsi ;
  326.  
  327.  'FIN' bloc2 ;
  328. *mess 'Fin Bloc2 ' ;
  329.  
  330.  'SI'((rv.'TRAN') et (NON Xdfdt));
  331.  Mess ' Incoherence : Algorithme Transitoire '
  332.       ' et pas d increment en temps';
  333.  Mess ' Il manque un operateur DFDT NS ou TSCA ';
  334.  QUITTER EXEC ;
  335.  'FINSI' ;
  336.  
  337.  NRES=NRES+1;
  338. *mess 'NRES=' NRES;
  339.  TABRES=chai 'FORMAT' '(A6,I1)' 'TABRES' NRES;
  340.  Si (NON('EXIST' rv TABRES)); rv.TABRES=TABLE ; Finsi ;
  341.  Si (NON ('EGA' ('TYPE' rv.TABRES) 'TABLE   '));
  342.   rv.TABRES=TABLE ;
  343.  Finsi ;
  344.  
  345.  rv.TABRES.'mau' =mau     ;
  346.  rv.TABRES.'sf'  =sf      ;
  347.  rv.TABRES.'clim'=rv.'CLIM' ;
  348.  rv.TABRES.'METHINV'=rv.'METHINV' ;
  349.  
  350. *****************************************************
  351. *              __      |\
  352. *   \         /_/o     | |   \|/
  353. *    \=>      \_\/      \|  --O--  PROJECTION
  354. *             /|\      /\    /|\
  355. *************/*|*\***********************************
  356.  
  357.  'SI' (rv.'PROJ');
  358.     SCHE = rv.'SCHEMAT';
  359.     cnvi1=chai 'FORMAT' '(I1,A3)' 1 (extr rv.'LISTINCO' 1) ;
  360.     cnvi2=chai 'FORMAT' '(I1,A3)' 2 (extr rv.'LISTINCO' 1) ;
  361.     cnvi3=chai 'FORMAT' '(I1,A3)' 3 (extr rv.'LISTINCO' 1) ;
  362.     cnppp=chai 'FORMAT' '(I1,A3)' 3 (extr rv.'LISTINCO' 2) ;
  363.   'SI' (testp1) ;
  364.     mess ' CALCUL de CMCT ' ;
  365.     Diag = doma (rv.'MODELE') 'XXDIAGSI' ;
  366.     Digv= (nomc Diag cnvi1) et (nomc Diag cnvi2) ;
  367.      Si ('EGA' (VALE DIME) 3);
  368.       Digv= Digv et (nomc Diag cnvi3);
  369.      Finsi ;
  370.  
  371. * Diag : Matrice diagonale de base ('SCAL')
  372. * Digv : Matrice diagonale de base ('1U' '2U' '3U')
  373. * Digvl: Matrice diagonale contenant les conditions aux limites
  374.  
  375.     IDiag = 'INVERSE' Digv ;
  376.  
  377.    Digvl= 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1;
  378.    IDigvl= 'INVERSE' Digvl ;
  379.  
  380. *
  381. * Calcul matrices MD-1
  382. *
  383.     MDM1='KMF' mam IDiag ;
  384.      Si ('EGA' SCHE 'BDF2'); MDM1=MDM1 * (((rv.'DELTAT')*2.)/3.) ;
  385.      Sinon ; MDM1=MDM1 * (rv.'DELTAT') ;
  386.      Finsi ;
  387. *
  388. * Calcul matrices C
  389. *
  390. *   msp mac =assi  numassi  'KBBT'
  391. *   (rv.'MODELE') (rv.'LISTINCO') (rv.'TYPEINCO') (rv.'TYPPRESS')
  392. *   MRIG  (1) (1.) ;
  393.  
  394.     matpc mac  = kops 'CMCTSPLT' mac ;
  395.  
  396.      Si ('EGA' (rv.'TYPPRESS') 'MSOMMET');
  397.       mess ' Cas des pressions continues' ;
  398.       incp= mots (extr rv.'LISTINCO' 2) ;
  399.       msp matpr =assi  numassi  'LAPN'
  400.       (rv.'MODELE') incp (mots 'TEMPERAT')
  401.       MRIG 'LINE' (1) (1.) ;
  402.  
  403.      Sinon ;
  404.       mess ' Cas des pressions discontinues' ;
  405.  
  406.       Si (rv.'TVNPC'); mac=mac et mvnpc; Finsi;
  407.        rv . 'MATC' = mac ;
  408. *  list mac 4 ;
  409.        scr mcr  = 'KOPS' 'MATRIK' ;
  410.         Si (rv.'TVNP') ;
  411.          rv . 'MBTR' = mar ;
  412.          Dunit=Idigv;
  413.          crt= kops 'CMCT' mac mar (Dunit)  ;
  414.          ctr= kops 'CMCT' mar mac (Dunit)  ;
  415.          rrt= kops 'CMCT' mar mar (Dunit)  ;
  416.          mcr= mcr et crt et ctr et rrt ;
  417.         Finsi ;
  418.  
  419. *     list mac 3 ;list IDigvl;
  420.        matpr = kops 'CMCT' mac mac IDigvl;
  421. *   mess 'APRES CMCT' ; list matpr 3 ;
  422.        matpr = matpc et matpr et mcr ;
  423.  
  424.      Finsi ;
  425.  
  426.        rv.'IDigvl'= IDigvl ;
  427.        rv.'MATP'  = matpr ;
  428.        rv.'MATC'  = mac   ;
  429.        rv.'MDM1'  = MDM1  ;
  430.  
  431.   'FINSI' ;
  432.  
  433.    IDigvl= rv.'IDigvl';
  434.    matpr = rv.'MATP' ;
  435.    mac   = rv.'MATC' ;
  436.    MDM1  = rv.'MDM1' ;
  437.  
  438.  
  439.    'SI' ('EXIST' (rv.'INCO') 'PINM1') ;
  440. **
  441. ** mess '           t  n-1' ;
  442. ** mess 'Calcul de C  P'    ;
  443. **
  444.  
  445.      pinm1= rv.'INCO'.'PINM1';
  446. *   cpi2=(extr rv.'LISTINCO' 2) ;
  447. *   pinm1= nomc (rv.'INCO'.(rv.'INC2')) cpi2;
  448.     cpre = 'KMF' (rv . 'MATC') pinm1 'TRAN' ;
  449.      Si (rv.'TVNP') ;
  450.       cxre = 'KMF' (rv . 'MBTR') pinm1 'TRAN' ;
  451.       cpre = cpre et cxre ;
  452.      Finsi ;
  453.  
  454.      Si ('EGA' (rv.'TYPPRESS') 'MSOMMET');
  455.       gradpres = 'KOPS' MDM1 '*' cpre ;
  456.      Sinon ;
  457.       gradpres = 'KOPS' MDM1 '*' cpre ;
  458.      Finsi ;
  459.       gradpres = nomc (mots 'UX' 'UY' 'UZ')
  460.                    (mots cnvi1 cnvi2 cnvi3) gradpres;
  461.  
  462.       oublier cpre ;
  463.       rv.'INCO'.'GRADPRES' = gradpres ;
  464.       sf = sf - (rv.'INCO'.'GRADPRES') ;
  465.       rv.'TABRES1'.'sf'=sf ;
  466.  
  467.    'FINSI' ;
  468.     clim = exco rv.'CLIM' (mots cnvi1 cnvi2 cnvi3) 'NOID' ;
  469.     rv.'TABRES1'.'clim'=clim ;
  470.  
  471. **********  FIN Traitement PRESSION ************************
  472. *   mess ' FIN Traitement PRESSION' ;
  473.  
  474.  'FINSI' ;
  475. *******************************************************
  476.  
  477.  
  478. *****************************************************
  479. *   \___\  /
  480. *        \/                     RESOLUTION
  481. *         \               \|/
  482. *************************--O--***********************
  483.     res mau  = assi numassi 'KOPS' 'MATRIK' ;
  484.     Si (rv.XRIG);
  485.     res = chan res 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
  486.     Finsi ;
  487. Repeter Bclres NRES ;
  488. *   mess ' Resolution ' NRES ;
  489.     TABRES=rv.(chai 'FORMAT' '(A6,I1)' 'TABRES' &Bclres);
  490.     MTIV = TABRES.'METHINV';
  491.     mau  = TABRES.'mau' ;
  492.     sf   = TABRES.'sf'  ;
  493.     clim = TABRES.'clim';
  494. *........................................................
  495.  Si (rv.'XRIG') ;
  496.     gclim=extr clim 'MAILLAGE' ;
  497.     rclim=bloq 'T' gclim ;
  498.     dclim=depi rclim (nomc 'Q' clim);
  499.     mau = mau et rclim ;
  500.     sf  = sf  et dclim ;
  501. *   list mau; list sf ;
  502.     res1 = 'RESO' mau sf ;
  503.     res = ASSI numassi res1 et res ;
  504.  Sinon ;
  505.      MATASS= mau ;
  506.      TYRENU=MTIV.'TYRENU';PCMLAG=MTIV.'PCMLAG';KTYPI =MTIV.'TYPINV';
  507.      XINIT=TEXT '   ' ;DI2 = TEXT '   ' ;DI3 = TEXT '   ' ;
  508.      DI4 = TEXT '   ' ;DI5 = TEXT '   ' ;
  509.      MAPREC=TEXT '   ';KPREC =TEXT '   ';PR3 = TEXT '   ' ;
  510.      PR4 = TEXT '   ' ;PR5 = TEXT '   ' ;
  511.    Si (('>EG' KTYPI 2) ET ('<EG' KTYPI 5));
  512.     XINIT mar  = assi numassi 'KOPS' 'MATRIK' ;
  513.     DI2=MTIV.'NITMAX';DI3=MTIV.'RESID';DI4=MTIV.'BCGSBTOL';
  514.      Si ('EGA' KTYPI 5);DI5=MTIV.'GMRESTRT';Finsi ;
  515.     MAPREC= mau ;
  516.     KPREC=MTIV.'PRECOND';
  517.      Si ('EGA' KPREC 4);PR3=MTIV.'MILURELX';Finsi ;
  518.      Si (('EGA' KPREC 5) OU ('EGA' KPREC 6));
  519.      PR4=MTIV.'ILUTLFIL';PR5=MTIV.'ILUTDTOL';
  520.      Finsi ;
  521.    Finsi ;
  522. *     list mau 3 ;
  523. *     list sf ;
  524.       res1 = ASSI numassi 'KRES'  'LL' mau
  525.      'TYPI' MTIV KTYPI 'CLIM' clim 'SMBR' sf
  526.      'IMPR' IMPKRES MATASS TYRENU PCMLAG
  527.      XINIT DI2 DI3 DI4 DI5
  528.      MAPREC KPREC PR3 PR4 PR5;
  529.       res = ASSI numassi res1 et res ;
  530.  Finsi;
  531. *........................................................
  532. Fin Bclres ;
  533.  
  534.    rv.inco.'Res'=res ;
  535.  
  536.  'SI' (rv.'PROJ');
  537. **
  538. *      mess 'RESOLUTION : ETAPE DE PROJECTION ' ;
  539. **   _n
  540. ** C U
  541.       cun = ('KMF' (rv . 'MATC') res)*(1./(rv.'DELTAT')) ;
  542.       Si (rv.'TVNP');
  543.       cxn = ('KMF' (rv . 'MBTR') res)*(1./(rv.'DELTAT')) ;
  544.       cun = cun et cxn ;
  545.       Finsi ;
  546.       clim=exco (rv.'CLIM') (mots (extr rv.'LISTINCO' 2)) 'NOID';
  547. *     MTIV=rv.'METHINV2' ;
  548.  
  549. *........................................................
  550.  Si (rv.'XRIG') ;
  551.     gclim=extr clim 'MAILLAGE' ;
  552.     rclim=bloq 'T' gclim ;
  553.     dclim=depi rclim (nomc 'Q' clim);
  554.     matpr = matpr et rclim ;
  555.     cun   = cun   et dclim ;
  556.     sl  = 'RESO' matpr cun ;
  557.  Sinon ;
  558.      MATASS= matpr ;
  559.      TYRENU=MTIV.'TYRENU';PCMLAG=MTIV.'PCMLAG';KTYPI =MTIV.'TYPINV';
  560.      XINIT=TEXT '   ' ;DI2 = TEXT '   ' ;DI3 = TEXT '   ' ;
  561.      DI4 = TEXT '   ' ;DI5 = TEXT '   ' ;
  562.      MAPREC=TEXT '   ';KPREC =TEXT '   ';PR3 = TEXT '   ' ;
  563.      PR4 = TEXT '   ' ;PR5 = TEXT '   ' ;
  564.    Si (('>EG' KTYPI 2) ET ('<EG' KTYPI 5));
  565.     XINIT mar  = assi numassi 'KOPS' 'MATRIK' ;
  566.     DI2=MTIV.'NITMAX';DI3=MTIV.'RESID';DI4=MTIV.'BCGSBTOL';
  567.      Si ('EGA' KTYPI 5);DI5=MTIV.'GMRESTRT';Finsi ;
  568.     MAPREC= matpr;
  569.     KPREC=MTIV.'PRECOND';
  570.      Si ('EGA' KPREC 4);PR3=MTIV.'MILURELX';Finsi ;
  571.      Si (('EGA' KPREC 5) OU ('EGA' KPREC 6));
  572.      PR4=MTIV.'ILUTLFIL';PR5=MTIV.'ILUTDTOL';
  573.      Finsi ;
  574.    Finsi ;
  575.       sl = ASSI numassi 'KRES'  'LL' matpr
  576.      'TYPI' MTIV KTYPI 'CLIM' clim 'SMBR' cun
  577.      'IMPR' IMPKRES MATASS TYRENU PCMLAG
  578.      XINIT DI2 DI3 DI4 DI5
  579.      MAPREC KPREC PR3 PR4 PR5;
  580.  Finsi;
  581. *........................................................
  582.       ctl = 'KMF' (rv . 'MATC') sl 'TRAN' ;
  583.       Si(rv.'TVNP') ;
  584.       cxl = 'KMF' (rv . 'MBTR') sl 'TRAN' ;
  585.       ctl = ctl + cxl ;
  586.       Finsi ;
  587.       a=KOPS IDigvl '*' ctl ;
  588.       au = nomc (mots 'UX' 'UY' 'UZ')(mots cnvi1 cnvi2 cnvi3) a ;
  589.       oublier ctl ;
  590.       res  =  res - (au*(rv.'DELTAT')) ;
  591.  'FINSI' ;
  592.  
  593.  
  594. *************************--O--***********************
  595. *   \                     /|\
  596. *   /\____                      RELAXATION
  597. *  /  \   \
  598. *****************************************************
  599.    nbic=dime rv.'LISTINCO';
  600.    omega=rv.'OMEGA';
  601.    NUPDT=rv.'PASDETPS'.'NUPASDT';
  602.    LRA= CHAI 'Pas de temps ' NUPDT ' ERREUR LINF : ' ;
  603.      'SI' ('MULT' NUPDT (rv.'FIDT'));
  604.        Si('EXIST' rv.'PASDETPS' 'NUPDT');
  605.         rv.'PASDETPS'.'NUPDT'=rv.'PASDETPS'.'NUPDT' et (PROG NUPDT) ;
  606.        Sinon ;
  607.         rv.'PASDETPS'.'NUPDT'= (PROG NUPDT) ;
  608.        Finsi ;
  609.      'FINSI' ;
  610.    Repeter Relax nbic;
  611.  
  612.     a= chai 'FORMAT' '(A3,I1)' 'INC' &relax ;
  613.  
  614.     Si (XBDF2)  ;
  615.     b= chai 'FORMAT' '(A3,I1)' 'IMC' &relax ;
  616.     rv.'INCO'.(rv.b)=rv.'INCO'.(rv.a);
  617. *   mess ' On fait ' (rv.a) ' -> ' (rv.b) ;
  618.     Finsi ;
  619.  
  620.     nc =extr rv.'LISTINCO' &relax;
  621.     nc1=mots nc     ;
  622.     nc2=mots 'SCAL' ;
  623.     ty=extr rv.'TYPEINCO' &relax;
  624. *   mess ' Composante ' nc ' Type inco ' ty ;
  625.     Si (('EGA' ty 'VITESSE ') et ('EGA' ('VALEUR' 'DIME') 2));
  626.     nc2=mots 'UX' 'UY' ;
  627.     a1= chai 'FORMAT' '(I1,A3)' 1 nc ;
  628.     a2= chai 'FORMAT' '(I1,A3)' 2 nc ;
  629.     nc1=mots a1 a2 ;
  630.     Finsi ;
  631.     Si (('EGA' ty 'VITESSE ') et ('EGA' ('VALEUR' 'DIME') 3));
  632.     nc2=mots 'UX' 'UY' 'UZ' ;
  633.     a1= chai 'FORMAT' '(I1,A3)' 1 nc ;
  634.     a2= chai 'FORMAT' '(I1,A3)' 2 nc ;
  635.     a3= chai 'FORMAT' '(I1,A3)' 3 nc ;
  636.     nc1=mots a1 a2 a3 ;
  637.     Finsi ;
  638.     ri= exco nc1 res nc2 'NOID';
  639.     nbc=dime (extr ri 'COMP') ;
  640.     Si('EGA' nbc 0);
  641.      Si('EGA' NUPDT 1);
  642.       Mess ' ATTENTION La (ou les) composante(s) suivantes ne sont pas'
  643.       ' calculees' ;
  644.       list nc1 ;
  645.      Finsi ;
  646.     Sinon ;
  647.      'SI' ('MULT' NUPDT (rv.'FIDT'));
  648.      eri=ri - rv.'INCO'.(rv.a);
  649.      ELI=MAXI ERI 'ABS';
  650.      ELI=(LOG (ELI + 1.0E-10))/(LOG 10.);
  651.      LRA= CHAI LRA '  ' (rv.a) ' ' ELI ;
  652.      ER = PROG ELI;
  653.      ERA=CHAI 'ERLI' (rv.a);
  654.       Si('EXIST' rv.'PASDETPS' ERA);
  655.        rv.'PASDETPS'.ERA=rv.'PASDETPS'.ERA et ER ;
  656.       Sinon ;
  657.        rv.'PASDETPS'.ERA= ER ;
  658.       Finsi ;
  659.      'FINSI' ;
  660.  
  661.     rv.'INCO'.(rv.a)=(omega*ri) + ((1. - omega)*rv.'INCO'.(rv.a)) ;
  662.     Finsi ;
  663.  
  664.    Fin Relax ;
  665.  
  666.      'SI' ('MULT' NUPDT (rv.'FIDT'));
  667.       MESS LRA ;
  668.      'FINSI' ;
  669.  
  670.  Fin BCLINT ;
  671.  
  672. *****  Avancement en temps
  673.  rv.inco.'Resnm'=res ;
  674.  rv.'PASDETPS'.'NUPASDT'=rv.'PASDETPS'.'NUPASDT' + 1 ;
  675. * Algorithme de projection  P(n+1) = P(n) + (sl / dt)
  676.     'SI' (rv.'PROJ');
  677.  
  678.     'SI' ('EXIST' (rv.'INCO') 'PINM1') ;
  679.  
  680. *      PNM1=(-1.)*rv.'INCO'.'PRESSION' ;
  681.        PNM1=(-1.)*rv.'INCO'.(rv.'INC2') ;
  682.        si (XBDF2) ;
  683.        PN=PNM1 + (1.5*sl) ;
  684.        sinon ;
  685.        PN=PNM1 + (sl) ;
  686.        finsi ;
  687.  
  688. *      rv.'INCO'.'PRESSION'=(-1.)*PN ;
  689.        rv.'INCO'.(rv.'INC2')=(-1.)*PN ;
  690.        rv.'INCO'.'PINM1'=PN ;
  691.      'SINON' ;
  692. *      rv.'INCO'.'PRESSION' = (-1.)*(sl) ;
  693.        rv.'INCO'.(rv.'INC2')= (-1.)*(sl) ;
  694.        rv.'INCO'.'PINM1'=(2.*sl);
  695.     'FINSI' ;
  696.     'OUBLIER' sl ;
  697.    'FINSI' ;
  698. *****  Avancement en temps  FIN
  699. Fin BCLTPS ;
  700.  
  701.  QUITTER EXEC ;
  702.  Finsi ;
  703. ********************************************************************************
  704. ********************************************************************************
  705. ****** Fin directive EQUA ******************************************************
  706. ********************************************************************************
  707. ********************************************************************************
  708.  
  709. 'SI' ('EGA'  (rv . 'NAVISTOK') 0) ;
  710.   EXAC rv ;
  711. 'QUITTER' EXEC ;
  712. 'FINSI'  ;
  713.  
  714. *******PROJ*************************************************************
  715. *******PROJ Trois possibilités *****************************************
  716. *******PROJ TMDM1 (vrai)_> Correction Gresho                  NON
  717. *******PROJ TGRAD (vrai)-> Formulation Gradient               NON
  718. *******PROJ TPNM2 (vrai)-> elimination end of step velocity   OUI
  719. *******PROJ TMDM1 = VRAI TGRAD = FAUX TPNM2 = FAUX  ancien algo
  720.          TYPROJ='VPI1'              ;
  721.        'SI' ('EXIST' RV 'TYPROJ')   ;
  722.          TYPROJ=RV.'TYPROJ'         ;
  723.           'SI' (non (exist (MOTS 'PSCT' 'VPI1' 'VPI2') TYPROJ));
  724.  Mess '*********************************************************' ;
  725.  Mess ' ERREUR   ERREUR   ERREUR   ERREUR   ERREUR   ERREUR     ' ;
  726.  Mess '                                                         ' ;
  727.              Mess 'Le mot ' TYPROJ ' n existe pas dans la liste.';
  728.              Mess 'Les méthodes de projection autorisées sont :';
  729.              Mess 'VPI1 et VPI2 '                               ;
  730.  Mess '*********************************************************' ;
  731.              erreur 21              ;
  732.              quitter EXEC           ;
  733.           'FINSI'                   ;
  734.        'FINSI'                      ;
  735.  
  736.        'SI' ('EGA' TYPROJ 'PSCT')   ;
  737.             TGRAD = FAUX ;
  738.             TMDM1 = FAUX ;
  739.             TPNM2 = FAUX ;
  740.        'FINSI'                      ;
  741.        'SI' ('EGA' TYPROJ 'VPI1')   ;
  742.             TGRAD = FAUX ;
  743.             TMDM1 = VRAI ;
  744.             TPNM2 = FAUX ;
  745.        'FINSI'                      ;
  746.        'SI' ('EGA' TYPROJ 'VPI2')   ;
  747.             TGRAD = FAUX ;
  748.             TMDM1 = FAUX ;
  749.             TPNM2 = VRAI ;
  750.        'FINSI'                      ;
  751.  
  752. *******PROJ Trois possibilités *****************************************
  753. *******PROJ*************************************************************
  754.  
  755.   nomvi=rv . 'NOMVI'  ;
  756.  
  757. 'SI'  ('EGA' ('VALEUR' 'DIME') 2) ;
  758.   vnul=0.D0 0.D0 ;
  759.   vuni=1.D0 1.D0 ;
  760.   cnvi1= 'MOT' ('TEXTE' ('CHAINE' 1 nomvi)) ;
  761.   cnvi2= 'MOT' ('TEXTE' ('CHAINE' 2 nomvi)) ;
  762.   lc= mots cnvi1 cnvi2 ;
  763.   lcu=mots 'UX' 'UY' ;
  764. 'SINON' ;
  765.   vnul=0.D0 0.D0 0.D0 ;
  766.   vuni=1.D0 1.D0 1.D0 ;
  767.   cnvi1= 'MOT' ('TEXTE' ('CHAINE' 1 nomvi)) ;
  768.   cnvi2= 'MOT' ('TEXTE' ('CHAINE' 2 nomvi)) ;
  769.   cnvi3= 'MOT' ('TEXTE' ('CHAINE' 3 nomvi)) ;
  770.   lc= mots cnvi1 cnvi2 cnvi3 ;
  771.   lcu=mots 'UX' 'UY' 'UZ' ;
  772. 'FINSI' ;
  773. 'SI' ('NON' ('EXISTE' rv 'OMEGA')) ;
  774.   omeg=1.D0 ;
  775. 'SINON' ;
  776.   omeg=rv . 'OMEGA' ;  
  777. 'FINSI' ;
  778.  
  779. testpr ='EXISTE' rv 'PRESSION' ;
  780. testprj ='EXISTE' rv 'PROJ' ;
  781. testran=testpr 'ET' ('EXISTE' rv 'CO') ;
  782.  
  783.   si testpr ;  rvp  = rv.'PRESSION' ;
  784.   achp  matpr= 'KOPS' 'MATRIK' ;
  785.   rvp.'MATP'= matpr;
  786.   finsi ;
  787.   si testprj;  rvp  = rv.'PROJ' ; Finsi ;
  788.  
  789.   'SI' (rv.'IMPR' >EG 1) ;
  790.   si testpr ;
  791.   mess 'Algorithme semi explicite (ANCIEN)'; 
  792.   mess '==================================';
  793.   Finsi ;
  794.   si testprj;
  795.   mess 'Algorithme de Projection';
  796.   mess '========================';
  797.   finsi ;
  798.   si ( non (testpr ou testprj)) ;
  799.   mess 'Algorithme standard implicite ou explicite ';
  800.   mess '===========================================';
  801.   finsi ;
  802.   'FINSI' ;
  803.  
  804. 'SI' ('NON' ('EXISTE' rv 'HIST')) ;
  805.   rv . 'HIST' = 'TABLE' ;
  806. 'FINSI' ;
  807.  
  808.  
  809. ITMA=(rv . 'ITMA') ;
  810. IMPTCRR=0 ;
  811. IMPKRES=0 ;
  812. 'SI' ('<EG' ITMA 1) ;
  813.   ITMA=1 ;
  814.   IMPTCRR=1 ;
  815. 'FINSI'  ;
  816.  Si (testprj ) ; IMPTCRR=1         ;  finsi ;
  817.  Si ( non (testpr ou testprj)) ; IMPTCRR=1 ; finsi ;
  818.  
  819.  
  820. * Gestion de la matrice de préconditionnement
  821. *
  822. * calass  : doit-on recalculer l'assemblage
  823. * calprec : doit-on recalculer le préconditionneur
  824. * calassp : idem pour la matrice de pression
  825. * calprecp: idem pour la matrice de pression
  826. * fcprect : fréquence de recalcul du préc. en fn du pas de temps
  827. * fcpreci : fréquence de recalcul du préc. dans la boucle
  828. *           d'itérations internes pour les non-linéarités
  829. * resmn : le résidu au pas de temps ou à l'itération interne précédente
  830. *
  831. calprec = VRAI ;
  832. fcprect = rv. 'METHINV' . 'FCPRECT' ;
  833. fcpreci = rv. 'METHINV' . 'FCPRECI' ;
  834. resmn maprec = 'KOPS' 'MATRIK' ;
  835. resmn1 maprec1 = 'KOPS' 'MATRIK' ;
  836. resmn2 maprec2 = 'KOPS' 'MATRIK' ;
  837. 'SI'(exist rv 'resmn');resmn=rv.'resmn' ; finsi ;
  838. 'SI'(non (exist rv 'calprec'));rv.'calprec'=VRAI;finsi;
  839. 'SI'(non (exist rv 'calass' ));rv.'calass' =VRAI;calass =VRAI;finsi;
  840. 'REPETER' bloc1 ITMA ;
  841. * mess ' DEBUT BLOC1 ';
  842.  
  843.    testp1 = FAUX  ;
  844.   si (testpr ou testprj) ;
  845.    Si CALPRE  ; testp1=VRAI ; Finsi ;
  846.    Si (non (exist rvp 'MATC')) ; testp1=VRAI ; Finsi ;
  847.   finsi ;
  848.  
  849.   'SI' ('MULT' (&bloc1 '+' 1) fcprect) ;
  850.     rv.'calprec' = VRAI ;
  851.   'FINSI' ;
  852.   'REPETER' bloci (rv . 'NITER') ;
  853.     'SI' ('MULT' (&bloci '+' 1) fcpreci) ;
  854.     rv.'calprec' = VRAI ;
  855.     'FINSI' ;
  856.  
  857.   'SI' (testprj) ;
  858.     calprecp=VRAI ;
  859.     calassp =VRAI ;
  860.     si (exist rvp 'calprecp' ) ; calprecp=rvp.'calprecp';
  861.     maprec1=rvp.'maprec1' ;finsi ;
  862.     si (exist rvp 'calassp' ) ; calassp=rvp.'calassp';
  863.     matass1=rvp.'matass1' ;finsi ;
  864.     fcprectp = rvp . 'METHINV' . 'FCPRECT' ;
  865.     'SI' ('MULT' (&bloc1 '+' 1) fcprectp) ;
  866.     calprecp= VRAI ;
  867.     'FINSI' ;
  868.   'FINSI' ;
  869.  
  870.  
  871.  
  872.     st mat  = 'KOPS' 'MATRIK' ;
  873.     sf mau  = 'KOPS' 'MATRIK' ;
  874.  
  875.  
  876.       mdfdt = 0 ;
  877.     'REPETER' bloc2 ('DIME' (rv . 'LISTOPER')) ;
  878.       nomper = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ;
  879.       notable= 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ;
  880. *     mess 'Operateur Bloc2 mdfdt ? ' nomper ;
  881.       mdfdt = mdfdt + rv. notable . 'KOPT' . 'KFORM' ;
  882.  
  883.       si (ega nomper 'DFDT    ');
  884.       ISCHT= (rv . notable.kopt.'ISCHT') ;
  885.       msi mai=('TEXTE' nomper) (rv . notable) ;
  886.       mat = mat 'ET' mai ;
  887.       st  = st  'ET' msi ;
  888.  
  889.       sinon ;
  890.  
  891. *     mess ' 1OPER =' nomper;
  892.       msi mai=('TEXTE' nomper) (rv . notable) ;
  893.       mau = mau 'ET' mai ;
  894.       sf  = sf  'ET' msi ;
  895.  
  896.       finsi ;
  897.  
  898.     'FIN' bloc2 ;
  899.  
  900.       s2 = sf et st ;
  901.       ma1 = mau 'ET' mat ;
  902.  
  903. **********  Traitement PRESSION ****************************
  904. *mess 'Traitement PRESSION';
  905.  
  906. 'SI' (exist rv 'rvpd') ;
  907.   rvpd = rv.'rvpd' ;
  908.   IDigv=rv.'IDigv';
  909.    Digv=rv.'Digv';
  910.   matpr= rvp.'MATP' ;
  911. 'FINSI' ;
  912.  
  913. */1    CALCUL de CMCT  pour PRESSION
  914. 'SI' (testpr et testp1) ;
  915.  
  916.   rvpd = (rvp.'DOMAINE') ;
  917.   Diago = doma rvpd 'XXDIAGSI' ;
  918.   si ( ega ('VALEUR' 'DIME') 2) ;
  919.   Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) ;
  920.   sinon ;
  921.   Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) et
  922.         ( exco Diago 'SCAL' cnvi3 ) ;
  923.   finsi ;
  924.   Digv = kcht rvpd vect sommet comp lc Digv ;
  925.  
  926.     'SI' ('EXISTE' rv 'CLIM') ;
  927.     rvp . 'CLIM' = rv . 'CLIM' ;
  928.     Digv = 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1;
  929.     'FINSI' ;
  930.  
  931.   IDigv= 'INVERSE' digv ;
  932.  
  933.     'SI' ('NON' ('EXISTE' rvp 'DIAGV')) ;
  934.     rvp . 'DIAGV' = digv ;
  935.     'FINSI' ;
  936.  
  937.   rvp . 'MATC' = 'KMAB' rvp ;
  938.   rvp . 'PRESSION' = 'KCHT' rvpd 'SCAL'
  939.                      'CENTRE' 0.D0 ;
  940.   rvp . 'GRADP'    = 'KCHT' rvpd 'VECT'
  941.                      'SOMMET' vnul ;
  942.  
  943.   achp  matpr= 'KOPS' 'MATRIK' ;
  944.   rvp . 'MATP' = matpr ;
  945.   rv.'rvpd'=rvpd;
  946.   rv.'IDigv'=IDigv;
  947.   rv.'Digv'=Digv;
  948.  
  949. 'FINSI' ;
  950.  
  951. */3   CALCUL de CMCT  pour PROJ
  952. 'SI' (testprj et testp1) ;
  953.  
  954.     Diago mma = 'KOPS' 'MATRIK' ;
  955.  
  956.       idfdt = 0 ;
  957.     'REPETER' blocj ('DIME' (rv . 'LISTOPER')) ;
  958.       nomper = 'EXTRAIRE' &blocj (rv . 'LISTOPER') ;
  959.       notable= 'MOT' ('TEXTE' ('CHAINE' &blocj nomper)) ;
  960. *?    mess 'Operateur ' nomper ;
  961.       si (ega nomper 'DFDT    ');
  962.       idfdt=idfdt + 1 ;
  963.       Diago = Diago et (doma (rv. notable . 'DOMZ') 'XXDIAGSI') ;
  964.       rvpd=rv. notable . 'DOMZ';
  965.       finsi ;
  966.     'FIN' blocj ;
  967.     si (ega idfdt 0) ; mess ' Pas de DFDT ?? ' ; erreur 21 ; finsi ;
  968.  
  969.   si ( ega ('VALEUR' 'DIME') 2) ;
  970.   Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) ;
  971.   sinon ;
  972.   Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) et
  973.         ( exco Diago 'SCAL' cnvi3 ) ;
  974.   finsi ;
  975.   Digv = kcht rvpd vect sommet comp lc Digv  ;
  976. * Matrice diagonale ne contenant pas les conditions aux limites : Diag
  977.      Diag  = copier Digv;
  978.      IDiag = 'INVERSE' Diag ;
  979.  
  980.   Digv = 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1;
  981.  
  982.   IDigv= 'INVERSE' digv ;
  983.  
  984.     'SI' ('NON' ('EXISTE' rvp 'DIAGV')) ;
  985.     rvp . 'DIAGV' = digv ;
  986.     'FINSI' ;
  987.  
  988.   rvp.'INCO'=rv.'INCO' ;
  989.  
  990.     sp map  = 'KOPS' 'MATRIK' ;
  991.     sr mar  = 'KOPS' 'MATRIK' ;
  992.     svnpc mvnpc = 'KOPS' 'MATRIK' ;
  993.     'REPETER' blocpj ('DIME' (rvp . 'LISTOPER')) ;
  994.       nomper = 'EXTRAIRE' &blocpj (rvp . 'LISTOPER') ;
  995.       notable= 'MOT' ('TEXTE' ('CHAINE' &blocpj nomper)) ;
  996.     mess 'Procedure PROJ Operateur : ' nomper ;
  997.  
  998.       si (EGA nomper 'KBBT');
  999.       rvp . notable . 'KOPT' . 'IKOMP' = 1 ;
  1000.       finsi ;
  1001.  
  1002.       msi mai=('TEXTE' nomper) (rvp . notable) ;
  1003.  
  1004. *******PROJ*************************************************************
  1005. *******PROJ TGRAD  calcul maig  (Formulation en Gradient)       DEBUT*** *******
  1006.   Si TGRAD ;
  1007.     mess 'Procedure PROJCT (Gradient) Operateur : ' nomper ;
  1008.       rvp . notable . 'KOPT' . 'IKOMP' = 0 ;
  1009.       msig maig=('TEXTE' nomper) (rvp . notable) ;
  1010.   Finsi ;
  1011. *******PROJ                                             -> MATG FIN*****
  1012. *******PROJ/////////////////////////////////////////////////////////////
  1013.  
  1014.       TVNP=FAUX;
  1015.       TVNPC=FAUX;
  1016.      Si (('EGA' nomper  'VNIMP   ') et (non (ega (rvp.'DISCPRES') 5)));
  1017.       TVNP=VRAI ;
  1018.       mar = mar 'ET' mai ;
  1019.       sr  = sr  'ET' msi ;
  1020.      Finsi ;
  1021.  
  1022. *     mess 'nomper=' nomper (rvp.'DISCPRES');
  1023.      Si (('EGA' nomper  'VNIMP   ') et (ega (rvp.'DISCPRES') 5));
  1024.       TVNPC=VRAI;
  1025.       mvnpc=mvnpc et mai ;
  1026.       svnpc=svnpc et msi ;
  1027.      Finsi ;
  1028.  
  1029.      Si (non ('EGA' nomper  'VNIMP   '));
  1030.       map = map 'ET' mai ;
  1031.       sp  = sp  'ET' msi ;
  1032.      Finsi ;
  1033.  
  1034.     'FIN' blocpj ;
  1035.  
  1036.   matpc mac  = kops 'CMCTSPLT' map ;
  1037. *******PROJ*************************************************************
  1038. *******PROJ TGRAD  (vrai  Formulation en Gradient) CMCTSPLT     DEBUT*** *******
  1039.   Si TGRAD ;
  1040.   matpcg macg  = kops 'CMCTSPLT' maig ;
  1041.   Finsi ;
  1042. *******PROJ                                             -> MATG FIN*****
  1043. *******PROJ/////////////////////////////////////////////////////////////
  1044.  
  1045.     Si TVNPC; mac=mac et mvnpc; Finsi;
  1046.   rvp . 'MATC' = mac ;
  1047.  
  1048. *******PROJ*************************************************************
  1049. *******PROJ TGRAD  (vrai  Formulation en Gradient) -> MATG      DEBUT*** *******
  1050.     SI TGRAD ;   rvp . 'MATG' = macg;  FINSI ;
  1051. *******PROJ                                             -> MATG FIN*****
  1052. *******PROJ/////////////////////////////////////////////////////////////
  1053.  
  1054.   scr mcr  = 'KOPS' 'MATRIK' ;
  1055.   Si TVNP ;
  1056.   rvp . 'MBTR' = mar ;
  1057.   Dunit=Idigv;
  1058.   crt= kops 'CMCT' mac mar (Dunit)  ;
  1059.   ctr= kops 'CMCT' mar mac (Dunit)  ;
  1060.   rrt= kops 'CMCT' mar mar (Dunit)  ;
  1061.   mcr= mcr et crt et ctr et rrt ;
  1062.   Finsi ;
  1063.   rvp . 'TVNP' = TVNP ;
  1064.  
  1065.   Si( ega (rvp.'DISCPRES') 5 ) ;
  1066.   mess ' Cas des pressions continues' ;
  1067.   matpr = matpc et mcr ;
  1068.   Sinon ;
  1069.   mess ' Cas des pressions discontinues' ;
  1070.  
  1071.   matpr = kops 'CMCT' mac mac IDigv ;
  1072.   matpr = matpc et matpr et mcr ;
  1073.  
  1074.   Finsi ;
  1075.  
  1076.  
  1077.   rvp . 'MATP' = matpr ;
  1078.   rv.'rvpd'=rvpd;
  1079.   rv.'IDigv'=IDigv;
  1080.   rv.'Digv'=Digv;
  1081. 'FINSI' ;
  1082. *************
  1083. ** Fin calcul CMCT pour testp1 et testprj
  1084. ******************************************
  1085. *                           t
  1086. ***************  Calcul de C p   Pression
  1087. 'SI' testpr ;
  1088.  
  1089.       dt = (rv . 'PASDETPS' . 'DELTAT') '*' (rv . 'ALFA') ;
  1090.       rvp . 'DELTAT' = dt ;
  1091.       f = 'COPIER' s2 ;
  1092.       u = rv . 'INCO' . nomvi ;
  1093.  
  1094.       lc = 'EXTRAIRE' digv  'COMP' ;
  1095.       fu = 'KCHT'     rvpd 'VECT' 'SOMMET'
  1096.                       'COMP' lc ('EXCO' f lc) ;
  1097.  
  1098.       'SI' ('EXISTE' rv 'CLIM') ;
  1099.         dti = -1.D0 '/' dt ;
  1100.         dm1f= 'KOPS' (dt '*' ('KOPS'  fu '*' IDigv))
  1101.                      'CLIM' (rv . 'CLIM') 3 ;
  1102.         dm1f= dti * dm1f ;
  1103.       'SINON'  ;
  1104.         dm1f= (-1.D0) '*' ('KOPS' fu '*' IDigv)  ;
  1105.       'FINSI'  ;
  1106.  
  1107.       rvp . 'PRESSION' = 'KMF' (rvp . 'MATC') dm1f ;
  1108.  
  1109.       'KRES' rvp (rvp . 'PRESSION')
  1110.              'BETA' (rvp . 'KBETA') (rvp . 'BETA')
  1111.              'PIMP' (rvp . 'KPIMP') (rvp . 'PIMP') ;
  1112.       rvp . 'GRADP' = 'KMTP' 1 (rvp . 'MATC')
  1113.                                (rvp . 'PRESSION') lc ;
  1114.       s2 = s2 + (rvp . 'GRADP') ;
  1115.       rv.'INCO'.'PRESSION'=rvp . 'PRESSION';
  1116.  
  1117. 'FINSI' ;
  1118.  
  1119. *******PROJ*************************************************************
  1120. *******PROJ MDM1 Ici on calcule M D-1 pour ensuite calculer Ctp DEBUT*** *******
  1121. *******PROJ on ne calcule M D-1 que si (rv 'MDM1') n'existe pas DEBUT*** *******
  1122. *******PROJ ou CALPRE VRAI                                      DEBUT*** *******
  1123. 'SI' testprj ;
  1124.  
  1125.       dt = (rv . 'PASDETPS' . 'DELTAT') '*' (rv . 'ALFA') ;
  1126.       rvp . 'DELTAT' = dt ;
  1127.  
  1128. ** Produit M D-1
  1129.  'SI' TMDM1                                  ;
  1130.   'SI' (('EXIST' rv 'MDM1') et (NON CALPRE)) ;
  1131.      MDM1=rv.'MDM1' ;
  1132.   'SINON' ;
  1133.      mess ' On calcule MD-1 ' ;
  1134.      stn matn= 'KOPS' 'MATRIK' ;
  1135.  
  1136.     'REPETER' blocpj1 ('DIME' (rvp . 'LISTOPER')) ;
  1137.       nomper1 = 'EXTRAIRE' &blocpj1 (rvp . 'LISTOPER') ;
  1138.       notable1= 'MOT' ('TEXTE' ('CHAINE' &blocpj1 nomper1)) ;
  1139. *    mess 'Operateur ' nomper1 ;
  1140.  
  1141.       si (ega nomper1 'KBBT    ');
  1142.       nomiv1= extr  (rvp . notable1 . 'LISTINCO')  1 ;
  1143.  
  1144.        'REPETER' blocpj2 ('DIME' (rv . 'LISTOPER')) ;
  1145.          nomper2 = 'EXTRAIRE' &blocpj2 (rv . 'LISTOPER') ;
  1146.          notable2= 'MOT' ('TEXTE' ('CHAINE' &blocpj2 nomper2)) ;
  1147. *        mess 'Operateur ' nomper2 ;
  1148.  
  1149.          si (ega nomper2 'DFDT    ');
  1150.           nomiv2= extr  (rv . notable2 . 'LISTINCO')  1 ;
  1151.            si ('EGA' nomiv2 nomiv1) ;
  1152.  
  1153.             domzp=rv . notable2 . 'DOMZ' ;
  1154.             msi mai='DFDT' (rv . notable2) ;
  1155.             matn = matn et mai ;
  1156.             sti = kcht domzp vect sommet comp lc vuni;
  1157.            finsi ;
  1158.          finsi ;
  1159.        'FIN' blocpj2 ;
  1160.       finsi ;
  1161.     'FIN' blocpj1 ;
  1162.  
  1163.   MDM1='KMF' matn IDiag ;
  1164.  
  1165.     Si(EGA ISCHT 1);
  1166.     MDM1=MDM1 * ((dt*2.)/3.) ;
  1167.     Sinon ;
  1168.     MDM1=MDM1 * dt ;
  1169.     Finsi ;
  1170.  
  1171.   rv.'MDM1' = MDM1 ;
  1172.  
  1173. *******PROJ                                   MDM1 -> rv.'MDM1' FIN*****
  1174. *******PROJ/////////////////////////////////////////////////////////////
  1175.   'FINSI';
  1176.  'FINSI';
  1177. *******PROJ*************************************************************
  1178. *******PROJ            t  n-1                                   DEBUT*** *******
  1179. *******PROJ Calcul de C  P                                      DEBUT*** *******
  1180.  
  1181.   TVNP=rvp.'TVNP' ;
  1182.  
  1183.   'SI' ('EXIST' (rv.'INCO') 'PRESSION') ;
  1184.        PPI = rv.'INCO'.'PRESSION'       ;
  1185.      'SI' ('EXIST' (rv.'INCO') 'PNM2')  ;
  1186. * TPNM2 elimination of end of step velocity (2 Pn - Pn-1)
  1187.        PPI = 2*(PPI) - rv.'INCO'.'PNM2' ;
  1188.      'FINSI'                            ;
  1189.  
  1190. *Formulation en u grad p
  1191.       Si (Exist rvp 'MATG');
  1192.       cpre = 'KMF' (rvp . 'MATG') PPI                    'TRAN' ;
  1193. *Formulation en p div u
  1194.       Sinon;
  1195.       cpre = 'KMF' (rvp . 'MATC') PPI                    'TRAN' ;
  1196.       Finsi ;
  1197.  
  1198.       Si TVNP ;
  1199.       cxre = 'KMF' (rvp . 'MBTR') PPI                    'TRAN' ;
  1200.       cpre = cpre et cxre ;
  1201.       Finsi ;
  1202.  
  1203. *******PROJ MDM1
  1204.       Si TMDM1 ;
  1205. * consistence selon (Gresho)
  1206.       gradpres = 'KOPS' MDM1 '*' cpre ;
  1207.       Sinon ;
  1208. * consistence selon (Guermond)
  1209.       gradpres =                 cpre     ;
  1210.       Finsi ;
  1211.  
  1212.       oublier cpre ;
  1213.       rv.'INCO'.'GRADPRES' = 'NOMC' lcu lc gradpres;
  1214.  
  1215. *Formulation en u grad p
  1216.       Si (Exist rvp 'MATG');
  1217.       s2 = sf - (rv.'INCO'.'GRADPRES') + st ;
  1218.       Sinon ;
  1219. *Formulation en p div u
  1220.       s2 = sf + (rv.'INCO'.'GRADPRES') + st     ;
  1221.       Finsi ;
  1222.  
  1223.   'FINSI' ;
  1224.  
  1225. *******PROJ                  t  n-1                             FIN*****
  1226. *******PROJ        MDM1 '*' C  P        -> rv.'INCO'.'GRADPRES' FIN*****
  1227. *******PROJ                                          t  n-1     FIN*****
  1228. *******PROJ         si (p div u)        -> S2 = F + C  P  + st  FIN*****
  1229. *******PROJ                                          t  n-1     FIN*****
  1230. *******PROJ         si (u grad p)       -> S2 = F - C  P  + st  FIN*****
  1231. *******PROJ/////////////////////////////////////////////////////////////
  1232. 'FINSI' ;
  1233.  
  1234. **********  FIN Traitement PRESSION ************************
  1235. *   mess ' FIN Traitement PRESSION' ;
  1236.  
  1237.     'SI' ('EXISTE' rv 'CLIM') ;
  1238.       s1 = rv . 'CLIM' ;
  1239.     'SINON' ;
  1240.       s1=' ' ;
  1241.     'FINSI'  ;
  1242.     rv . 'S2' = s2 ;
  1243.  
  1244. *******************************************************************
  1245. * Résolution hors QDM méthode de projection
  1246.  'SI'((NON testprj) ou (EGA mdfdt 0));
  1247. *mess '* Résolution hors QDM méthode de projection ';
  1248.  
  1249. *   mess ' ====================> calprec '; list calprec ;
  1250.     'SI' rv.'calass' ;
  1251.      rv . 'METHINV' . 'MATASS' =ma1 ;
  1252.      rv.'calass'=FAUX ;
  1253.     'FINSI' ;
  1254.  
  1255.     'SI' rv.'calprec' ;
  1256.      rv . 'METHINV' . 'MAPREC' = ma1    ;
  1257.      rv.'calprec'=FAUX ;
  1258.     'FINSI' ;
  1259.  
  1260.     rv . 'METHINV' . 'XINIT'  = resmn ;
  1261.  
  1262.     res = 'KRES'  ma1 'TYPI' (rv . 'METHINV')
  1263.                       'CLIM' s1
  1264.                       'SMBR' s2
  1265.                       'IMPR' IMPKRES ;
  1266.  
  1267.  'FINSI' ;
  1268.  
  1269. 'SI' (testprj et (non(EGA mdfdt 0)));
  1270. *******PROJ*************************************************************
  1271. *******PROJ  Résolution QDM                                     DEBUT*** *******
  1272. *******PROJ  (On éclate les résolutions)                        DEBUT*** *******
  1273.  
  1274. * Résolution QDM méthode de projection (On éclate les résolutions)
  1275. *mess ' Résolution QDM méthode de projection';
  1276.    res mm1= kops 'MATRIK' ;
  1277.    liscom=extr s1 'COMP';
  1278.    nbcom=dime liscom ;
  1279. *   MEss ' NBCOM=' nbcom ;
  1280.  
  1281.    Si(NON(EXIST RV 'TABRES')); RV.'TABRES'=TABLE ; Finsi ;
  1282.    Si(NON(EGA ('TYPE' RV.'TABRES') 'TABLE')); RV.'TABRES'=TABLE;Finsi;
  1283.  
  1284.    repeter Bclcom nbcom  ;
  1285.    nmc=extr liscom &Bclcom;
  1286. *   Mess ' Composante : ' nmc ;
  1287.     ma1i=kops 'EXTRCOMP' nmc 'SCAL' ma1;
  1288.     s1i    = exco s1    nmc 'NOID' ;
  1289.     s2i    = exco s2    nmc 'NOID' ;
  1290.     resmni = exco resmn nmc 'NOID' ;
  1291.  
  1292.     MATASSI =chai 'FORMAT' '(A6,I1)' 'MATASS' &Bclcom;
  1293.     MAPRECI =chai 'FORMAT' '(A6,I1)' 'MAPREC' &Bclcom;
  1294.  
  1295.     'SI' (rv.'calass' ou (NON ('EXIST' (RV.'TABRES') MATASSI))) ;
  1296.      rv . 'TABRES' . MATASSI = ma1i;
  1297.      rv.'calass'=FAUX ;
  1298.     'FINSI' ;
  1299.  
  1300.     'SI' (rv.'calprec' ou (NON ('EXIST' (RV.'TABRES') MAPRECI)));
  1301.      mess 'On recalcule le preconditionneur ';
  1302.      rv . 'TABRES' . MAPRECI = ma1i;
  1303.      rv.'calprec'=FAUX ;
  1304.     'FINSI' ;
  1305.  
  1306.     rv . 'METHINV' . 'MATASS' = rv . 'TABRES' . MATASSI ;
  1307.     rv . 'METHINV' . 'MAPREC' = rv . 'TABRES' . MAPRECI ;
  1308.     rv . 'METHINV' . 'XINIT'  = resmni ;
  1309.  
  1310.  
  1311.     resi = 'KRES'  ma1i 'TYPI' (rv . 'METHINV')
  1312.                       'CLIM' s1i
  1313.                       'SMBR' s2i
  1314.                       'IMPR' IMPKRES ;
  1315.     resi = nomc resi nmc ;
  1316.     res=res et resi ;
  1317.     Fin Bclcom ;
  1318.     oubli ma1 ;
  1319. *mess ' Fin Résolution QDM méthode de projection';
  1320. * Fin Résolution QDM méthode de projection
  1321. *******PROJ  Résolution QDM                              -> res FIN*****
  1322. *******PROJ  (On éclate les résolutions)                        FIN*****
  1323. *******PROJ/////////////////////////////////////////////////////////////
  1324. 'FINSI' ;
  1325.  
  1326. 'SI' testprj ;
  1327. *******PROJ*************************************************************
  1328. *******PROJ  ETAPE DE PROJECTION                                DEBUT*** *******
  1329. *******PROJ  on calcul cun (alias c U tilde)                    DEBUT*** *******
  1330. *     mess 'ETAPE DE PROJECTION ' ;
  1331. *   _n
  1332. * C U
  1333.       cun = ('KMF' (rvp . 'MATC') res) ;
  1334.       Si TVNP ;
  1335.       cxn = ('KMF' (rvp . 'MBTR') res) ;
  1336.       cun = cun et cxn ;
  1337.       Finsi ;
  1338.  
  1339. *    mess ' On calcule les seconds membres de l équation de pression ' ;
  1340. *    mess ' s ils existent (opérateurs FIMP)                         ' ;
  1341.  
  1342.     'REPETER' blocpj ('DIME' (rvp . 'LISTOPER')) ;
  1343.       nomper = 'EXTRAIRE' &blocpj (rvp . 'LISTOPER') ;
  1344.       notable= 'MOT' ('TEXTE' ('CHAINE' &blocpj nomper)) ;
  1345.       si (ega nomper 'FIMP    ');
  1346. *   mess 'Seconds membres de l équation de pression Opérateur ' nomper ;
  1347.       msi mai=('TEXTE' nomper) (rvp . notable) ;
  1348.       cun = cun et msi ;
  1349.       finsi ;
  1350.     'FIN' blocpj ;
  1351.  
  1352.       cun = cun * (-1./dt) ;
  1353.  
  1354. *   mess ' ====================> calprecp '; list calprecp ; 
  1355.       'SI' calassp;
  1356.       matass1=matpr ;
  1357.       rvp.'matass1'=matass1 ;
  1358.       'FINSI' ;
  1359.  
  1360.      'SI' calprecp                          ;
  1361.       'SI' (RV.IMPR >EG 1)                  ;
  1362.       'MESS' ' On preconditionne CMCT '     ;
  1363.       'FINSI'                               ;
  1364.        maprec1 = matpr ;
  1365.        rvp.'maprec1'=maprec1 ;
  1366.       'FINSI' ;
  1367.       rvp . 'METHINV' . 'MATASS' = matass1 ;
  1368.       rvp . 'METHINV' . 'MAPREC' = maprec1 ;
  1369.       rvp . 'METHINV' . 'XINIT'  = resmn2  ;
  1370.  
  1371.                     sl= 'KRES' matpr 'TYPI' (rvp . 'METHINV')
  1372.                         'CLIM' (rvp . CLIM )
  1373.                         'SMBR' cun
  1374.                         'IMPR' IMPKRES  ;  oublier cun ;
  1375.  
  1376.       oublier resmn2 ; resmn2 = sl ;
  1377.       ctl = 'KMF' (rvp . 'MATC') sl 'TRAN' ;
  1378.       Si TVNP ;
  1379.       cxl = 'KMF' (rvp . 'MBTR') sl 'TRAN' ;
  1380.       ctl = ctl + cxl ;
  1381.       Finsi ;
  1382.  
  1383. * elimination of end of step velocity (Si exist PNM2) (on ne corrige pas)
  1384.      'SI'(NON(TPNM2));
  1385.      'SI' (NON ('EXIST' (rv.'INCO') 'PNM2')) ;
  1386.       a=nomc lcu lc ('KOPS' IDigv '*' ctl) ;
  1387.       res  =  res + (a*dt) ;
  1388.      'FINSI' ;
  1389.      'FINSI' ;
  1390.  
  1391.       oublier ctl ;
  1392.  
  1393. *******PROJ                                -> (alias c U tilde) FIN*****
  1394. *******PROJ                                -> sl (alias lambda) FIN*****
  1395. *******PROJ                                -> on corrige u      FIN*****
  1396. *******PROJ/////////////////////////////////////////////////////////////
  1397. 'FINSI' ;
  1398.  
  1399. *****  Avancement en temps
  1400. *?  Si (testprj et ('EGA' mdfdt 0)) ; IMPTCRR=RV.'IMPR' ;  finsi ;
  1401.  Si ('EGA' mdfdt 0) ; IMPTCRR=RV.'IMPR' ;  finsi ;
  1402.     IMPKRES=0 ;
  1403. *   mess 'On passe dans TCRR ' ;
  1404.     eps = 'TCRR' res omeg (rv . 'INCO') 'IMPR' IMPTCRR ;
  1405.  
  1406.     'OUBLIER' resmn ;
  1407.     resmn = res ;
  1408.     'MENAGE' ;
  1409.   'FIN' bloci ;
  1410.   irt=0 ;
  1411.  
  1412. * mess ' ITMA= ' (rv . 'ITMA') ' mdfdt= ' mdfdt ;
  1413.   'SI' ('EGA' (rv . 'ITMA') 0) ;
  1414. *   mess ' ON PASSE PLUS DS TCNM>>>>>>>>>>>>>>' ;
  1415.      irt = 'TCNM' rv 'NOUP';
  1416.   'SINON'                                         ;
  1417. *   mess ' ON PASSE DS TCNM>>>>' (rv . 'ITMA') mdfdt;
  1418.      irt = 'TCNM' rv ;
  1419.   'FINSI' ;
  1420.  
  1421. *****  Avancement en temps  algorithme de projection
  1422.  
  1423.      'SI' (testprj) ;
  1424.        'SI' calprecp;
  1425.          calprecp= FAUX ;
  1426.        'SINON' ;
  1427.        'SI' ('NON' calassp) ;
  1428.                 'OUBLIER' matpr ;
  1429.        'FINSI' ;         
  1430.        'FINSI' ;
  1431.        calassp=FAUX ;
  1432.        rvp.'calassp'=calassp ;
  1433.        rvp.'calprecp'=calprecp ;
  1434.      'FINSI' ;
  1435.  
  1436. *******PROJ*************************************************************
  1437. *******PROJ  P(n+1) = P(n) + sl                                 DEBUT*** *******
  1438. * Algorithme de projection  P(n+1) = P(n) + sl
  1439.    'SI' testprj ;
  1440.  
  1441.      'SI' (NON ('EXIST' (rv.'INCO') 'PRESSION'));
  1442.          rv.'INCO'.'PRESSION' = sl;
  1443.          Si TPNM2 ; rv.'INCO'.'PNM2' = sl ; Finsi ;
  1444.      'SINON' ;
  1445.          PNM1=rv.'INCO'.'PRESSION';
  1446.         si ( EGA ISCHT 1) ;
  1447.          PN=PNM1 + (1.5*sl) ;
  1448.         sinon ;
  1449.          PN=PNM1 + sl ;
  1450.         finsi ;
  1451.  
  1452.        Si TPNM2 ; rv.'INCO'.'PNM2'=PNM1 ;  Finsi ;
  1453.        rv.'INCO'.'PRESSION'=PN   ;
  1454.      'FINSI' ;
  1455.  
  1456.      Si (EGA TYPROJ 'PSCT') ;
  1457.      rv.'INCO'.'PRESSION' = sl;
  1458.      Finsi ;
  1459.  
  1460.        'OUBLIER' sl ;
  1461.    'FINSI' ;
  1462. *******PROJ  P(n+1) = P(n) + sl                                 FIN*****
  1463. *******PROJ/////////////////////////////////////////////////////////////
  1464. *****  Avancement en temps  FIN
  1465.  
  1466.  
  1467.   'SI' testran ;
  1468.     rv . 'CO' . 'VITESSE' ='KOPS' (rv . 'INCO' . nomvi)
  1469.                              '-'  (rv . 'SEDIM') ;
  1470.     k = 'ABS' (rv . 'INCO' . 'KN') ; 
  1471.     e = 'ABS' (rv . 'INCO' . 'EN') ;
  1472.     k = 'KOPS' ('KOPS' k '*' k) '/' e ;
  1473.     dif = 'KOPS' ('KOPS' k '*' 0.09) '+' (rv . 'COEF') ;
  1474.     rv . 'CO' . 'DIFFU'   = 'NOEL' rvpd dif ;
  1475.     rv . 'CO' . 'TEMPERA' = rv . 'INCO' . 'CN' ;
  1476.   'FINSI' ;
  1477.   'MENAGE' ;
  1478.  
  1479.   'SI' ('EGA' irt 1) ;
  1480.     'MESSAGE' ' Temps final atteint : '
  1481.               (rv . 'PASDETPS' . 'TPS') ;
  1482.     'QUITTER' bloc1 ;
  1483.   'FINSI' ;
  1484. 'FIN' bloc1 ;
  1485.  
  1486. 'SI' testpr ;
  1487.   rvp . 'PRESSION' = 'KCHT' rvpd 'SCAL' 'CENTRE'
  1488.                             (rvp . 'PRESSION') ;
  1489.   rvp . 'PN'       = 'ELNO' rvpd
  1490.                             (rvp . 'PRESSION') ;
  1491. 'FINSI' ;
  1492. ************************    E X E C ************************************
  1493. 'FINPROC' ;
  1494.  
  1495.  
  1496. GRAPH=VRAI ;
  1497. GRAPH=FAUX ;
  1498. COMPLET=VRAI;
  1499. COMPLET=FAUX ;
  1500.  
  1501. Si (NON Complet);
  1502. itmax= 80 ;
  1503. DISCR='LINE' ;
  1504. d1 = 0.04; d2 = 0.12; 
  1505. DT=.5;
  1506. Ra = 1.e6 ;
  1507. Sinon ;
  1508. *Ra = 1.e5 ;
  1509. *itmax= 600 ;
  1510. Ra = 1.e6 ;
  1511. itmax= 600 ;
  1512. DISCR='LINE' ;
  1513. d1 = 0.01; d2 = 0.07;
  1514. DT=.5;
  1515. Finsi ;
  1516.  
  1517. KPRESS='MSOMMET' ;
  1518. KSUPG ='CENTREE' ;
  1519.  
  1520. ****************************** 
  1521. *ESTIMATION DE LA CONVERGENCE* 
  1522. ****************************** 
  1523. DEBPROC CALCUL; 
  1524. ARGU RV*'TABLE' NUPASDT*ENTIER;
  1525. * Actualisation du terme GB
  1526. TN=rv.INCO.'TN';
  1527. mt=doma $mt maillage ;
  1528.  rv.INCO.'GB'= (manu 'CHPO' mt 1 'UX' 0. 'NATURE' 'DIFFUS') et
  1529.                (nomc 'UY' (TN - 0.5) 'NATURE' 'DIFFUS' );
  1530.  'FINPROC' ;
  1531.  
  1532. ************ 
  1533. * MAILLAGE * 
  1534. ************ 
  1536. OPTI ELEM QUA8;
  1537. crit = 0.0000001; 
  1538. xxll = 1.;
  1539. p1 = 0. 0.; 
  1540. p15 = (xxll/2.) 0.; 
  1541. p2 = xxll 0.; 
  1542. p25 = xxll 0.5; 
  1543. p3 = xxll 1.; 
  1544. p35 = (xxll/2.) 1.; 
  1545. p4 = 0. 1.; 
  1546. p45 = 0. 0.5; 
  1547.  
  1548. bas1 = p1 d dini d1 dfin d2 p15; 
  1549. bas2 = p15 d dini d2 dfin d1 p2;
  1550. bas = (bas1 et bas2);elim crit bas; 
  1551. cdro = p2 d dini d1 dfin d2 p25 d dini d2 dfin d1 p3; 
  1552. haut = p3 d dini d1 dfin d2 p35 d dini d2 dfin d1 p4; 
  1553. cgau = p4 d dini d1 dfin d2 p45 d dini d2 dfin d1 p1; 
  1554. coupx = bas plus (0. 0.5) ;
  1555. coupy=cgau plus (0.5 0.) ;
  1556. cnt = bas et cdro et haut et cgau; 
  1557. mt = bas cdro haut cgau daller; 
  1558. bord = cont mt; 
  1559. *mt = surf bord; 
  1560. elim (mt et coupx) 1.e-5;
  1561. Si GRAPH ; trace mt ; Finsi ;
  1562. *option donn 5;
  1563. orienter mt; 
  1564. ******** 
  1565. * DOMA * 
  1566. ******** 
  1567.  
  1568. Mmt = chan QUAF mt ;
  1569. $mt = MODE Mmt 'NAVIER_STOKES' DISCR ;
  1570. mt = doma $mt maillage;
  1571.  
  1572. Mcgau=chan cgau QUAF ; 
  1573. $cgau=mode Mcgau 'NAVIER_STOKES' DISCR ;
  1574.  
  1575. Mcdro=chan cdro QUAF ; 
  1576. $cdro=mode Mcdro 'NAVIER_STOKES' DISCR ;
  1577.  
  1578. Mcnt=chan cnt QUAF ; 
  1579. $cnt=mode Mcnt 'NAVIER_STOKES' DISCR ;
  1580.  
  1581. Mcoupx = chan coupx 'QUAF' ;
  1582. $coupx=mode Mcoupx 'NAVIER_STOKES' DISCR ;
  1583.  
  1584. Mcoupy = chan coupy 'QUAF' ;
  1585. $coupy=mode Mcoupy 'NAVIER_STOKES' DISCR ;
  1586.  
  1587. Elim (Mmt et Mcgau et Mcdro et Mcnt et Mcoupx
  1588. et Mcoupy) 1.e-5 ;
  1589. cgau=doma $cgau maillage ;
  1590. cdro=doma $cdro maillage ;
  1591. cnt=doma $cnt maillage ;
  1592. coupx=doma $coupx maillage ;
  1593. coupy=doma $coupy maillage ;
  1594.  
  1595. *option donn 5; 
  1596. ************** 
  1597. * PARAMETRES * 
  1598. ************** 
  1599. Pr = 0.71; 
  1600. Gr = Ra/Pr; 
  1601. NU = 1/(Gr**0.5); 
  1602. ALF= NU/Pr; 
  1603. gb = kcht $mt vect sommet (0. -1.) ;
  1604. uref = 0.2 ;
  1605.  
  1606. *********************** 
  1607. * CREATION DES TABLES *
  1608. *********************** 
  1609. ****
  1610. * La cavité est fermée il faut imposer la pression en 1 point !
  1611. prep1=doma $mt kpress;
  1612. bcp=elem  prep1 POI1 (lect 1 ) ;
  1613.  
  1614. RV = EQEX EQUA  'U' 'P' $mt 'VITESSE' 'PRESSION' KPRESS
  1615.  'PROJ' DT  'EUL_IMPL' 'UN' 'PN'
  1616.  'OMEGA' 1.    'NITER' 1  ITMA 1 'FIDT' 1
  1617.  'OPTI' 'EF' 'IMPL' KSUPG 'MMDIAGO'
  1618.  'ZONE' $MT 'OPER' 'NS'   1. NU 'UN' 'GB'    'INCO' 'U'
  1619.  'ZONE' $MT 'OPER' 'KBBT' -1.                'INCO' 'U' 'P'
  1620.  'CLIM' 'U' 'UIMP' cnt 0. 'U' 'VIMP' cnt 0. 'P' TIMP bcp 0.
  1621. ;
  1622.  
  1623.  rv.'METHINV'.TYPINV=3 ;
  1624.  rv.'METHINV'.IMPINV=0 ;
  1625.  rv.'METHINV'.NITMAX=400;
  1626.  rv.'METHINV'.PRECOND=3 ;
  1627.  rv.'METHINV'.RESID  =1.e-10;
  1628.  rv.'METHINV' . 'FCPRECT'=1   ;
  1629.  rv.'METHINV' . 'FCPRECI'=1   ;
  1630.  
  1631.  
  1632. RT = EQEX EQUA  'T' $mt 'TEMPERAT'
  1633.  'TRAN' DT  'EUL_IMPL' 'TN'
  1634. *'TRAN' DT  'BDF2' 'TN' 'TNM'
  1635.  'OMEGA' 1.    'NITER' 1  ITMA 1   'FIDT' 1
  1636.  'OPTI' 'EF' 'IMPL' KSUPG 'MMDIAGO'
  1637.  'ZONE' $MT  'OPER' 'TSCAL' 1. ALF 'UN' 0. 'INCO' 'T'
  1638.  'CLIM' 'T' 'TIMP'  cgau   1. 'T' 'TIMP' cdro 0.
  1639. ; 
  1640.  rt.'METHINV'.TYPINV=3 ;
  1641.  rt.'METHINV'.IMPINV=0 ;
  1642.  rt.'METHINV'.NITMAX=400;
  1643.  rt.'METHINV'.PRECOND=3 ;
  1644.  rt.'METHINV'.RESID  =1.e-10;
  1645.  rt.'METHINV' . 'FCPRECT'=1   ;
  1646.  rt.'METHINV' . 'FCPRECI'=1   ;
  1647.  
  1648. ****
  1649.  
  1650. RV.INCO = TABLE INCO;
  1651. RT.INCO = RV.INCO   ;
  1652. RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.); 
  1653. corsx=1. - (coor 1 (doma $mt sommet)) ;
  1654. RV.INCO.'TN' = kcht $MT SCAL SOMMET corsx;
  1655. RV.INCO.'TNM' = kcht $MT SCAL SOMMET 0.5;
  1656. RV.INCO.'PN' = kcht $MT SCAL KPRESS 0. ;
  1657. RV.INCO.'UNM1' = kcht $MT VECT SOMMET (1.E-5 1.E-5);
  1658. RV.INCO.'IT' = PROG 1; 
  1659. RV.INCO.'ER' = PROG 0.; 
  1660. *RV.INCO.'GB' =  gb ;
  1661. tn=RV.INCO.'TN' ;
  1662. rv.INCO.'GB'= (manu 'CHPO' mt 1 'UX' 0. 'NATURE' 'DIFFUS') et
  1663.               (nomc 'UY' ((TN - 0.5)) 'NATURE' 'DIFFUS' );
  1664.  
  1665. rv.'FIDT'=5 ;
  1666. rt.'FIDT'=5 ;
  1667. Repeter Bclt itmax ;
  1668.  
  1669. calcul RV &bclt ;
  1670. EXEC RV;
  1671. EXEC RT ;
  1672.  
  1673. Fin Bclt ;
  1674.  
  1675. ************* 
  1676. * RESULTATS * 
  1677. ************* 
  1678. *option sauv 'SAV';
  1679. *sauv Pr Ra Gr NU ALF gb uref RV;
  1680. *************
  1681. OPTI ISOV SULI; 
  1682.  
  1683.  
  1684. un=(RV.INCO.'UN') * (Pr*(Gr**0.5)) ;
  1685. unch = vect un (0.3/(Pr*(Gr**0.5))) UX UY JAUNE; 
  1686.  
  1687. u1y = exco un 'UX' ;
  1688. u2x = exco un 'UY' ;
  1689. evu=evol 'CHPO'  u2x  coupx ;
  1690. evv=evol 'CHPO'  u1y  coupy ;
  1691.  
  1692. TN= RV.INCO.'TN' ;
  1693. GRADT = KOPS TN 'GRAD' $MT;
  1694. DTDX = KCHT $MT 'SCAL' 'CENTRE' (EXCO 'UX' GRADT); 
  1695. DTDX = (-1)*(ELNO $MT DTDX); 
  1696.  
  1697. EVOLd = EVOL 'CHPO' DTDX SCAL cdro ;
  1698. EVOLg = EVOL 'CHPO' DTDX SCAL cgau ;
  1699. evol1=evold et evolg ;
  1700. LISTE EVOL1; 
  1701. TAB1 = TABLE; 
  1702. TAB1.1 = 'MARQ LOSA';
  1703. TAB1.2 = 'MARQ CROI'; 
  1704.  
  1705. Si GRAPH ; 
  1706. DESS EVOL1 'XBOR' 0. 1. 'GRIL' TAB1 'TITRE' ' Nusselt' ; 
  1707. trac RV.INCO.'TN' mt (cont mt) 14 'TITRE' ' Isothermes'  ;  
  1708. trace unch mt (cont mt) 'TITRE' ' Vitesses' ; 
  1709. dess evu  'TITRE' 'Coupe ox à y=1/2 (uy)' ;
  1710. dess evv  'TITRE' 'Coupe oy à x=1/2 (ux)' ;
  1711.  
  1712. TAB1.'TITRE'=TABLE ;
  1713. TAB1.'TITRE' . 1 = mot 'Vitesse '     ;
  1714. TAB1.'TITRE' . 2 = mot 'Temperature '     ;
  1715. ERUN=rv.'PASDETPS'.'ERLIUN' ;
  1716. ERTN=rt.'PASDETPS'.'ERLITN' ;
  1717. EVTP=rt.'PASDETPS'.'NUPDT'  ;
  1718. EVU=EVOL 'MANU' 'Pas de Tps' EVTP 'LOG--E--inf' ERUN;
  1719. EVT=EVOL 'MANU' 'Pas de Tps' EVTP 'LOG--E--inf' ERTN;
  1720. EVOLR=EVU et EVT ;
  1721. dess EVOLR 'XBOR' 0. itmax 'YBOR' -10.0 0.0 LEGE TAB1
  1722. 'TITRE' 'Erreurs L-infini U et T ';
  1723. Finsi;
  1724.  
  1725. Nusd = (redu DTDX cdro);
  1726. Nusg = (redu DTDX cgau);
  1727.  
  1728. Numg = somt (Nusg*(doma $cgau 'XXDIAGSI')) ;
  1729. Numd = somt (Nusd*(doma $cdro 'XXDIAGSI')) ;
  1730.  
  1731. un2=KOPS un '*' (Pr*(Gr**0.5)); 
  1732. sw = KOPS un 'ROT' $mt ; 
  1733. rk = EQEX $mt 'OPTI' 'EF' 'IMPL' 
  1734. ZONE $mt OPER LAPN -1. INCO 'PSI' 
  1735. ZONE $mt OPER FIMP sw INCO 'PSI' 
  1736. 'CLIM' 'PSI' 'TIMP' cnt 0.; 
  1737. rk.'INCO'=table 'INCO' ; 
  1738. rk.'INCO'.'PSI'=kcht $mt scal sommet 0. ; 
  1739. EXec rk ;
  1740. psi=rk.'INCO'.'PSI'; 
  1741. psi2=kops psi '*' (Pr*(Gr**0.5)); 
  1742.  
  1743. Si GRAPH ; 
  1744. trac psi mt CNT 14 'TITRE' 'Fonction de courant '; 
  1745. Finsi;
  1746.  
  1747. opti 'ECHO' 0 ;
  1748. MESS '----------- Ra = ' Ra '------------------------------------';
  1749. MESS '                     NUSSELT MOYEN   NUSSELT MAX   NUSSELT MIN' ;
  1750. MESS 'CAST3M : Gauche     ' NUMG '  ' (MAXI Nusg) '  ' (MINI Nusg)    ;
  1751. MESS 'CAST3M : Droite     ' NUMD '  ' (MAXI Nusd) '  ' (MINI Nusd)    ;
  1752. Si (EGA Ra 1.e5);
  1753. MESS 'Vahl Davis : Ra 1.e6    4.519              7.717         0.729';
  1754. Finsi ;
  1755. Si (EGA Ra 1.e6);
  1756. MESS 'Vahl Davis : Ra 1.e6    8.8               17.925         0.989';
  1757. Finsi ;
  1758. MESS '---------------------------------------------------------------';
  1759. MESS '                         Umax            Vmax   ' ;
  1760. MESS 'CAST3M      :          '
  1761.  (MAXI (exco un 'UX')) (MAXI (exco un 'UY'));
  1762. Si (EGA Ra 1.e5);
  1763. MESS 'Vahl Davis  : Ra 1.e6    34.73            68.59  ' ;
  1764. Finsi ;
  1765. Si (EGA Ra 1.e6);
  1766. MESS 'Vahl Davis  : Ra 1.e6    64.63           219.36  ' ;
  1767. Finsi ;
  1768. MESS '---------------------------------------------------------------';
  1769. MESS '                          Psi max         ' ;
  1770. MESS 'CAST3M      :           '  (MAXI PSI)       ;
  1771. Si (EGA Ra 1.e5);
  1772. MESS 'Vahl Davis  : Ra 1.e6      9.612          ' ;
  1773. Finsi;
  1774. Si (EGA Ra 1.e6);
  1775. MESS 'Vahl Davis  : Ra 1.e6     16.75           ' ;
  1776. Finsi;
  1777. MESS '---------------------------------------------------------------';
  1778. opti 'ECHO' 1 ; 
  1779.  
  1780. Si (NON COMPLET) ;
  1781.  
  1782. ern=abs( NUMD - 7.635)/7.635 ;
  1783. ernm=abs((MAXI Nusg) - 12.37)/12.37;
  1784. eru=abs((MAXI (exco un 'UY')) - 252.5)/252.5;
  1785. erp=abs((MAXI PSI) - 18.74)/18.74;
  1786. mess 'ern = ' ern ' ernm = ' ernm ' eru = ' eru  ' erp = ' erp ;
  1787. Si (ern  > 0.002); erreur 5 ; Finsi ;
  1788. Si (ernm > 0.002); erreur 5 ; Finsi ;
  1789. Si (eru  > 0.002); erreur 5 ; Finsi ;
  1790. Si (erp  > 0.002); erreur 5 ; Finsi ;
  1791.  
  1792. Finsi ;
  1793.  
  1794.  FIN ;
  1795.  
  1796.  
  1797.  
  1798.  
  1799.  
  1800.  
  1801.  
  1802.  
  1803.  
  1804.  
  1805.  
  1806.  
  1807.  
  1808.  

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