*****************************************************
************************************************************************
* Section : Fluides Thermique
************************************************************************
* fichier :  vahldavis.dgibi                        *
** modifie le 15/06/2014 passage EQPR -> EQEX       *
*****************************************************
COMPLET = FAUX ;
SI ( COMPLET ) ;
      d1 = 0.02  ; 
      d2 = 0.1;
      NITER = 5000 ;
      CFL = 1.0 ;
SINON ;
      d1 = 0.05 ;
      d2 = 0.5 ;
      NITER = 600 ;
      CFL = 1.0 ;
FINSI ;
GRAPH = FAUX ;
KPRESS='CENTRE';
DISCR ='MACRO';
BETA=1.;


***************************************************************
*               CAVITE CARREE - VAHL-DAVIS                    *
* A. Chene/H. Paillere TTMF Aout 1997                         *
***************************************************************

******************************
*ESTIMATION DE LA CONVERGENCE*
******************************

DEBPROC CALCUL;
ARGU RVX*'TABLE';

RV = RVX.'EQEX';

DD = RV.PASDETPS.'NUPASDT';
NN = DD/5;
LO = (DD-(5*NN)) EGA 0;

SI (LO);

UN   = RV.INCO.'UN';
UNM1 = RV.INCO.'UNM1';

unx = kcht $MT scal sommet (exco 'UX' un);
unm1x = kcht $MT scal sommet (exco 'UX' unm1);
uny = kcht $MT scal sommet (exco 'UY' un);
unm1y = kcht $MT scal sommet (exco 'UY' unm1);

ERRX = KOPS unx '-' unm1x;
ERRY = KOPS uny '-' unm1y;
ELIX = MAXI ERRX 'ABS';
ELIY = MAXI ERRY 'ABS';
ELIX = (LOG (ELIX + 1.0E-10))/(LOG 10.);
ELIY = (LOG (ELIY + 1.0E-10))/(LOG 10.);
MESSAGE 'ITER' RV.PASDETPS.'NUPASDT'   'ERREUR LINF' ELIX ELIY;
RV.INCO.'UNM1' = KCHT $MT 'VECT' 'SOMMET' (RV.INCO.'UN');
IT = PROG RV.PASDETPS.'NUPASDT';
ER = PROG ELIY;
RV.INCO.'IT' = (RV.INCO.'IT') ET IT;
RV.INCO.'ER' = (RV.INCO.'ER') ET ER;

FINSI;
as2 ama1  = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1 ;

************
* MAILLAGE *
************

OPTI DIME 2;
OPTI ELEM QUA8;

p1  = 0.  0.;
p15 = 0.5 0.;
p2  = 1.  0.;
p25 = 1.  0.5;
p3  = 1.  1.;
p35 = 0.5 1.;
p4  = 0.  1.;
p45 = 0.  0.5;

bas  = p1 d dini d1 dfin d2 p15 d dini d2 dfin d1 p2;
cdro = p2 d dini d1 dfin d2 p25 d dini d2 dfin d1 p3;
haut = p3 d dini d1 dfin d2 p35 d dini d2 dfin d1 p4;
cgau = p4 d dini d1 dfin d2 p45 d dini d2 dfin d1 p1;

cnt = bas et cdro et haut et cgau;
mt = bas cdro haut cgau daller;

********
* MODE *
********

Mmt= chan mt quaf ;
$mt = MODE Mmt 'NAVIER_STOKES' DISCR ;
doma $mt 'IMPR' ;
mt = doma $mt 'MAILLAGE' ;

**************
* PARAMETRES *
**************

Pr = 0.71;
Ra = 1.e6;
Gr = Ra/Pr;
NU = 1/(Gr**0.5);
ALF= NU/Pr;
gb = 0 -1;
uref = 1;

***********************
* CREATION DES TABLES *
***********************

RV = EQEX $MT 'ITMA' NITER 'ALFA' CFL
 'ZONE' $MT 'OPER' CALCUL
 'OPTI' 'SUPG' 
 'ZONE' $MT 'OPER' 'NS' NU GB 'TN' 0.5 'INCO' 'UN'
 'OPTI' 'SUPG' 'MMPG'
 'ZONE' $MT 'OPER' 'TSCAL' ALF 'UN' 0. 'INCO' 'TN'
 'OPTI' 'CENTREE'
  ZONE  $MT 'OPER' 'DFDT' 1. 'UN' 'DELTAT' INCO 'UN'
  ZONE  $MT 'OPER' 'DFDT' 1. 'TN' 'DELTAT' INCO 'TN'
;

RV = EQEX RV
 'CLIM' 'UN' 'UIMP' cnt 0.
        'UN' 'VIMP' cnt 0.
        'TN' 'TIMP' cgau 1.
        'TN' 'TIMP' cdro 0.;

RVP = EQEX 'OPTI' 'EF' KPRESS
 'ZONE' $MT  OPER  KBBT  -1. beta  INCO 'UN' 'PRES'
;

    rvp.'METHINV'.TYPINV=1 ;
    rvp.'METHINV'.IMPINV=0 ;
    rvp.'METHINV'.NITMAX=300;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-8 ;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;

  RV.'PROJ'     =RVP ;

 RV.INCO = TABLE INCO;
 RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.);
 RV.INCO.'PRES'= kcht $MT SCAL KPRESS 0.;
 corx = coor 1 mt;
 RV.INCO.'TN' = kcht $MT SCAL SOMMET (1.-corx);

 RV.INCO.'UNM1' = kcht $MT VECT SOMMET (1.E-5 1.E-5);
 RV.INCO.'IT' = PROG 1;
 RV.INCO.'ER' = PROG 0.;

EXEC RV;

*************
* RESULTATS *
*************

Mcdro=chan cdro quaf ;
elim (mt et cdro) 1.e-3 ;
$cdro = mode Mcdro 'NAVIER_STOKES' DISCR;
DY = DOMA $cdro 'VOLUME';
DYT = SOMT DY;
elim (Mmt et Mcdro)1.e-3;
  
GRADT = KOPS RV.INCO.'TN' 'GRAD' $MT;
DTDX = KCHT $MT 'SCAL' 'CENTRE' (EXCO 'UX' GRADT);
DTDX = ELNO $MT DTDX;

DTDXp = KCHT $CDRO 'SCAL' 'SOMMET' DTDX;
DTDXp = KOPS DTDXp '*' (-1.); 

EVOL1 = EVOL 'CHPO' DTDXp SCAL CDRO;
EVOL1 = EVOL1 COUL ROUG ;

DTDXpc = NOEL $CDRO DTDXp;
num = KOPS DTDXpc '*' DY;
num = SOMT NUM;
num = (-1)*NUM/DYT;

MESSAGE 'NUSSELT MOYEN' NUM;
MESSAGE 'NUSSELT MAX' (MAXI DTDXpc);
MESSAGE 'NUSSELT MIN' (MINI DTDXpc);

SI ( (MAXI DTDXp) < 12.00 ) ;
        ERREUR 5 ;
FINSI ;
SI ( ABS (NUM) < 7.55 ) ;
        ERREUR 5 ;
FINSI ;
SI ( (MINI RV.INCO.'ER') > -4.95 ) ;
        ERREUR 5 ;
FINSI ;

SI ( GRAPH ) ;
OPTI ISOV SULI;
trace mt 'TITR' 'MAILLAGE' ;
trac RV.INCO.'TN' mt (cont mt) 14 'TITR' 'TEMPERATURE' ;

unch = vect RV.INCO.'UN' 1. UX UY JAUNE;
trace unch mt (cont mt) 'TITR' 'CHAMP DE VITESSE' ;

EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
     (RV.INCO.'ER') ;
dess EVOL4 'XBOR' 0. 10000. 'YBOR' -10.0 0.0 
        'TITR' 'CONVERGENCE VERS LE STATIONNAIRE' ;

TAB1 = TABLE ;
TAB1.1 = 'MARQ LOSA';
DESS EVOL1 'XBOR' 0. 1. 'GRIL' TAB1 'TITR' 'NUSSELT A LA PAROI' ;

un = RV.INCO.'UN';
un2=KOPS un '*' (Pr*(Gr**0.5));
sw = (-1.) * (KOPS un 'ROT' $mt) ;
rk = EQEX $mt 'OPTI' 'EF' 'IMPL'
ZONE  $mt OPER LAPN 1.    INCO 'PSI'
ZONE  $mt OPER FIMP sw    INCO 'PSI'
'CLIM' 'PSI' 'TIMP' (cgau et haut et bas et cdro) 0.;
rk.'INCO'=table 'INCO' ;
rk.'INCO'.'PSI'=kcht $mt scal sommet 0. ;
EXEC rk ;
psi=rk.'INCO'.'PSI';
psi2=kops psi '*' (Pr*(Gr**0.5)) ;
LISTE (MAXI PSI) ; 
LISTE (MINI PSI) ;
trac psi mt CNT 14 'TITR' 'FONCTION DE COURANT' ;
FINSI ;

FIN ;