* fichier :  dvisp.dgibi
* CAVITE CARREE ­ VAHL­DAVIS Methode de projection implicite *
* Malvina Renesson Aout 1999 *
***************************************************************
GRAPH=VRAI ;
GRAPH=FAUX ;
COMPLET=VRAI;
 COMPLET=FAUX ;
TCPT   = FAUX;
TKPR   = VRAI;
TRESOU = FAUX;
IMPARA = FAUX;

Si (NON Complet);
 itmax= 80 ;
 itmax= 10 ;
DISCR='LINE' ;
*DISCR='MACRO';
*DISCR='QUAF';
d1 = 0.04; d2 = 0.12;
d1 = 0.06; d2 = 0.18;
*d1 = 0.01; d2 = 0.03;
** d1 = 0.2 ; d2 = 0.3 ;
 d1 = 0.01 ; d2 = 0.02 ;
*d1 = 0.005; d2 = 0.01 ;
DT=.25;
Ra = 1.e6 ;
Sinon ;
*Ra = 1.e5 ;
*itmax= 600 ;
Ra = 1.e6 ;
itmax= 600 ;
DISCR='LINE' ;
d1 = 0.01; d2 = 0.07;
DT=.25;

Finsi ;

KSUPG ='CENTREE' ;
KPRESS='MSOMMET' ;
*KPRESS='CENTREP1' ;
*KPRESS='CENTRE' ;

******************************
*ESTIMATION DE LA CONVERGENCE*
******************************
DEBPROC CALCUL;
ARGU RVX*'TABLE' MDZ/'MMODEL' MDZ/'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;
crit = 0.0000001;
xxll = 1.;
p1 = 0. 0.;
p15 = (xxll/2.) 0.;
p2 = xxll 0.;
p25 = xxll 0.5;
p3 = xxll 1.;
p35 = (xxll/2.) 1.;
p4 = 0. 1.;
p45 = 0. 0.5;

bas1 = p1 d dini d1 dfin d2 p15;
bas2 = p15 d dini d2 dfin d1 p2;
bas = (bas1 et bas2);elim crit bas;
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;
coupx = bas plus (0. 0.5) ;
coupy=cgau plus (0.5 0.) ;
cnt = bas et cdro et haut et cgau;
mt = bas cdro haut cgau daller;
bord = cont mt;
*mt = surf bord;
elim (mt et coupx) 1.e-5;
Si GRAPH ; trace mt ; Finsi ;
*option donn 5;
orienter mt;
********
* DOMA *
********

Mmt = chan QUAF mt ;
$mt = MODE Mmt 'NAVIER_STOKES' DISCR ;
mt = doma $mt maillage;

Mcgau=chan cgau QUAF ;
$cgau=mode Mcgau 'NAVIER_STOKES' DISCR ;

Mcdro=chan cdro QUAF ;
$cdro=mode Mcdro 'NAVIER_STOKES' DISCR ;

Mcnt=chan cnt QUAF ;
$cnt=mode Mcnt 'NAVIER_STOKES' DISCR ;

Mcoupx = chan coupx 'QUAF' ;
$coupx=mode Mcoupx 'NAVIER_STOKES' DISCR ;

Mcoupy = chan coupy 'QUAF' ;
$coupy=mode Mcoupy 'NAVIER_STOKES' DISCR ;

Elim (Mmt et Mcgau et Mcdro et Mcnt et Mcoupx
et Mcoupy) 1.e-5 ;
cgau=doma $cgau maillage ;
cdro=doma $cdro maillage ;
cnt=doma $cnt maillage ;
coupx=doma $coupx maillage ;
coupy=doma $coupy maillage ;

*option donn 5;
**************
* PARAMETRES *
**************
Pr = 0.71;
Gr = Ra/Pr;
NU = 1/(Gr**0.5);
ALF= NU/Pr;
gb = kcht $mt vect centre (0. -1.);
uref = 0.2 ;

***********************
* CREATION DES TABLES *
***********************
RV = EQEX  'OMEGA' 1.    'NITER' 1  ITMA itmax   'FIDT' 1
  'OPTI' 'EF' 'IMPL' KSUPG
'ZONE' $MT 'OPER' CALCUL
'ZONE' $MT 'OPER' 'NS' 1. 'UN'   NU 'GB' 'TN' 0.5 'INCO' 'UN'
'ZONE' $MT 'OPER' 'TSCAL' 1. 'UN'   ALF 0. 'INCO' 'TN'
'OPTI'  'CENTREE'
 ZONE  $mt    OPER  DFDT 1. 'UN' DT  INCO 'UN'
 ZONE  $mt    OPER  DFDT 1. 'TN' DT  INCO 'TN'
;
****

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

 rv.'METHINV'.TYPINV=1 ;
 rv.'METHINV'.IMPINV=0 ;
 rv.'METHINV'.NITMAX=400;
 rv.'METHINV'.PRECOND=3 ;
 rv.'METHINV'.RESID  =1.e-10;
 rv. 'METHINV' . 'FCPRECT'=1   ;
 rv. 'METHINV' . 'FCPRECI'=1   ;


****
* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt kpress;
bcp=elem  prep1 POI1 (lect 1 ) ;

RVP = EQEX
'OPTI' 'EF' KPRESS
ZONE  $mt    OPER  KBBT  -1.  INCO 'UN' 'PRES'
CLIM PRES TIMP bcp 0. ;

    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.INCO = TABLE INCO;
RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.);
mt=doma $mt 'MAILLAGE';
cx=coor 1 mt;
cy=coor 2 mt;
un1=(nomc cx 'UY') - (nomc cy 'UX');
RV.INCO.'UN1'= un1;   
ung= vect un1 0.5 ux uy jaune ;
*trace ung mt;
*RV.INCO.'TN' = kcht $MT SCAL SOMMET (1.-corsx);
RV.INCO.'TN' = kcht $MT SCAL SOMMET 0.5 ;
RV.INCO.'PRES' = kcht $MT SCAL KPRESS 0. ;
RV.INCO.'UNM1' = kcht $MT VECT SOMMET (1.E-5 1.E-5);
RV.INCO.'IT' = PROG 1;
RV.INCO.'ER' = PROG 0.;
RV.INCO.'GB' =  gb ;
RV.'PROJ' = RVP;


 RV.'TCPT'  = TCPT  ;
 RV.'TKPR'  = TKPR  ;
 RV.'TRESOU'= TRESOU;
 RV.'IMPARA'= IMPARA;

 rv.'ITMA'=1;
temps;
 REPETER BLOC itmax;
EXEC RV    ;
 FIN BLOC;

temps;

*************
* RESULTATS *
*************
*option sauv 'SAV';
*sauv Pr Ra Gr NU ALF gb uref RV;
*************
OPTI ISOV SULI;


un=(RV.INCO.'UN') * (Pr*(Gr**0.5)) ;
unch = vect un (0.3/(Pr*(Gr**0.5))) UX UY JAUNE;

u1y = exco un 'UX' ;
u2x = exco un 'UY' ;
evu=evol 'CHPO'  u2x  coupx ;
evv=evol 'CHPO'  u1y  coupy ;

EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG---E---inf'
(RV.INCO.'ER');

GRADT = KOPS RV.INCO.'TN' 'GRAD' $MT;
DTDX = KCHT $MT 'SCAL' 'CENTRE' (EXCO 'UX' GRADT);
DTDX = (-1)*(ELNO $MT DTDX);

EVOLd = EVOL 'CHPO' DTDX SCAL cdro ;
EVOLg = EVOL 'CHPO' DTDX SCAL cgau ;
evol1=evold et evolg ;
LISTE EVOL1;
TAB1 = TABLE;
TAB1.1 = 'MARQ LOSA';
TAB1.2 = 'MARQ CROI';

Si GRAPH ;
DESS EVOL1 'XBOR' 0. 1. 'GRIL' TAB1 'TITRE' ' Nusselt' ;
trac RV.INCO.'TN' mt (cont mt) 14 'TITRE' ' Isothermes'  ;
trace unch mt (cont mt) 'TITRE' ' Vitesses' ;
dess evu  'TITRE' 'Coupe ox à y=1/2 (uy)' ;
dess evv  'TITRE' 'Coupe oy à x=1/2 (ux)' ;
dess EVOL4 'XBOR' 0. itmax 'YBOR' -10.0 0.0;
Finsi;

Nusd = (redu DTDX cdro);
Nusg = (redu DTDX cgau);

Numg = somt (Nusg*(doma $cgau 'XXDIAGSI')) ;
Numd = somt (Nusd*(doma $cdro 'XXDIAGSI')) ;

un2=KOPS un '*' (Pr*(Gr**0.5));
sw = KOPS un 'ROT' $mt ;
rk = EQEX $mt 'OPTI' 'EF' 'IMPL'
ZONE $mt OPER LAPN -1. INCO 'PSI'
*ZONE $mt OPER KONV 1. 'UN' 1. INCO 'PSI'
 ZONE $mt OPER FIMP sw INCO 'PSI'
'CLIM' 'PSI' 'TIMP' cnt 0.;
rk.'INCO'=table 'INCO' ;
rk.'INCO'.'PSI'=kcht $mt scal sommet 0. ;
 rk.'METHINV'.TYPINV=1 ;


RK.'TCPT'=TCPT;
EXec rk ;
psi=rk.'INCO'.'PSI';
psi2=kops psi '*' (Pr*(Gr**0.5));

Si GRAPH ;
trac psi mt CNT 14 'TITRE' 'Fonction de courant ';
Finsi;

opti 'ECHO' 0 ;
MESS '----------- Ra = ' Ra '------------------------------------';
MESS '                     NUSSELT MOYEN   NUSSELT MAX   NUSSELT MIN' ;
MESS 'CAST3M : Gauche     ' NUMG '  ' (MAXI Nusg) '  ' (MINI Nusg)    ;
MESS 'CAST3M : Droite     ' NUMD '  ' (MAXI Nusd) '  ' (MINI Nusd)    ;
Si (EGA Ra 1.e5);
MESS 'Vahl Davis : Ra 1.e6    4.519              7.717         0.729';
Finsi ;
Si (EGA Ra 1.e6);
MESS 'Vahl Davis : Ra 1.e6    8.8               17.925         0.989';
Finsi ;
MESS '---------------------------------------------------------------';
MESS '                         Umax            Vmax   ' ;
MESS 'CAST3M      :          '
 (MAXI (exco un 'UX')) (MAXI (exco un 'UY'));
Si (EGA Ra 1.e5);
MESS 'Vahl Davis  : Ra 1.e6    34.73            68.59  ' ;
Finsi ;
Si (EGA Ra 1.e6);
MESS 'Vahl Davis  : Ra 1.e6    64.63           219.36  ' ;
Finsi ;
MESS '---------------------------------------------------------------';
MESS '                          Psi max         ' ;
MESS 'CAST3M      :           '  (MAXI PSI)       ;
Si (EGA Ra 1.e5);
MESS 'Vahl Davis  : Ra 1.e6      9.612          ' ;
Finsi;
Si (EGA Ra 1.e6);
MESS 'Vahl Davis  : Ra 1.e6     16.75           ' ;
Finsi;
MESS '---------------------------------------------------------------';
opti 'ECHO' 1 ;

FIN ;
opti donn 5;
Si (NON COMPLET) ;

ern=abs( NUMD - 7.25)/7.25 ; ernm=abs((MAXI Nusg) - 11.5)/11.5 ;
eru=abs((MAXI (exco un 'UY')) - 277.)/277.;
erp=abs((MAXI PSI) - 25.3)/25.3;
mess 'ern = ' ern ' ernm = ' ernm ' eru = ' eru  ' erp = ' erp ;
Si (ern  > 0.005); erreur 5 ; Finsi ;
Si (ernm > 0.005); erreur 5 ; Finsi ;
Si (eru  > 0.005); erreur 5 ; Finsi ;
Si (erp  > 0.005); erreur 5 ; Finsi ;

Finsi ;

FIN ;