* fichier :  gridturb.dgibi
************************************************************************
* Section : Fluides Transitoire
************************************************************************
**********************************************************
* GRID TURBULENCE : convection of homogeneous turbulence *
* Analysis of the K-Epsilon TURBULENCE MODEL             *
* Mohammadi/Pironneau p. 74     (Wiley)                  *
* H. PAILLERE/TTMF/AVRIL 1997 (à l'origine 1/2 explicite)*
* Version implicite MARS 2000                            *
* Teste K-Epsilon en implicite méthode de projection     *
**********************************************************
* WARNING : Il faut que Nt/Nu > 1

*opti trace 'PSC';
GRAPH = VRAI ;
GRAPH = FAUX ;
KPRES='MSOMMET';
DISCR='LINE';
KSUPG='SUPGDC';
*
        D = 1 ;
        NITMA =200 ;
        LCFL = PROG 0.1 0.5 1. 10.;
        LCFL = PROG 1. 10. 50. 100.;
        lerkm= prog -1.5 -2. -4. -4. ;
        lerem= prog -1.5 -2. -4. -4. ;

NX = D*40 ;
NY = D ;

opti dime 2 ;
opti elem qua8 ;
opti isov suli ;

******************************************************************
* PROCEDURE POUR ESTIMER LA CONVERGENCE VERS L'ETAT STATIONNAIRE *
******************************************************************

DEBPROC CALCUL ;
ARGU RVX*'TABLE' ;

RV = RVX.'EQEX' ;

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

EN   = RV.INCO.'EN' ;
ENM1 = RV.INCO.'ENM1' ;
KN   = RV.INCO.'KN' ;
KNM1 = RV.INCO.'KNM1' ;

ERRE = KOPS EN '-' ENM1 ;
ELIE = MAXI ERRE 'ABS' ;
ELIE = (LOG (ELIE + 1.0E-20))/(LOG 20.) ;
ERRK = KOPS KN '-' KNM1 ;
ELIK = MAXI ERRK 'ABS' ;
ELIK = (LOG (ELIK + 1.0E-20))/(LOG 20.) ;
MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' '   ERREUR LINF ' ELIK ELIE ;
RV.INCO.'ENM1' = KCHT $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'EN') ;
RV.INCO.'KNM1' = KCHT $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'KN') ;
IT = PROG RV.PASDETPS.'NUPASDT' ;
ERK = PROG ELIK ;
ERE = PROG ELIE ;
RV.INCO.'IT' = (RV.INCO.'IT') ET IT ;
RV.INCO.'ERK' = (RV.INCO.'ERK') ET ERK ;
RV.INCO.'ERE' = (RV.INCO.'ERE') ET ERE ;
FINSI ;
 as2 ama1  = 'KOPS' 'MATRIK' ;
 FINPROC as2 ama1 ;

***********************
* GEOMETRY OF PROBLEM *
***********************

A1 = 0.0 0.0 ;
A2 = 1.0 0.0 ;
A3 = 1.0 0.02;
A4 = 0.0 0.02;

A1A2 = A1 'DROI' NX A2 ;
A2A3 = A2 'DROI' NY A3 ;
A3A4 = A3 'DROI' NX A4 ;
A4A1 = A4 'DROI' NY A1 ;

DOMTOT = 'DALLER' A1A2 A2A3 A3A4 A4A1 'PLAN' ;
DOMTOT = CHAN DOMTOT QUAF ;
A4A1   = CHAN A4A1   QUAF ;
A2A3   = CHAN A2A3   QUAF ;
A1A2   = CHAN A1A2   QUAF ;
A3A4   = CHAN A3A4   QUAF ;
ELIM (DOMTOT et A4A1 et A2A3 et A3A4 et A1A2) 1.E-9 ;
MCOTE  = A1A2 et A3A4;

$DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' DISCR ;
DOMA $DOMTOT IMPR ;
$INLET  = 'MODE' A4A1 'NAVIER_STOKES' DISCR ;
 INLET  = 'DOMA' $INLET 'MAILLAGE' ;
$OUTLET = 'MODE' A2A3 'NAVIER_STOKES' DISCR ;
 OUTLET = 'DOMA' $OUTLET 'MAILLAGE' ;
 poutlet= 'DOMA' $OUTLET 'MSOMMET'  ;
$COTE   = 'MODE' MCOTE 'NAVIER_STOKES' DISCR ;
COTE    = DOMA $COTE MAILLAGE ;

DOMTOT = DOMA $DOMTOT MAILLAGE ;
CNT = CONT DOMTOT ;

*****************
* PHYSICAL DATA *
*****************

Uinlet = 1.71;
Kinlet = 0.013;
Einlet = 0.029;
deltax = (coor 1 A2)/nx ;
Nu = 1.E-3 ;

*****************************************
* ANALYTICAL SOLUTION FOR K(x) AND E(x) *
*****************************************

C2 = 1.92 ;
X1 = COOR 1 (DOMA $DOMTOT SOMMET) ;
TINLET = KINLET/EINLET;
TAN=(((C2 - 1.)*X1) + (UINLET*TINLET))/UINLET;
M=1. - C2 ;
KAN=KINLET * (( 1. - ((M*X1)/UINLET/TINLET))**(1./M));
EAN=KAN*(INVE TAN) ;


********************
* SET OF EQUATIONS *
********************

RV = EQEX $DOMTOT 'ITMA' NITMA  'OMEGA' 1. 'NITER' 1 'FIDT' 20
        'ZONE' $DOMTOT 'OPER' CALCUL
        'OPTI' KSUPG 'IMPL' EF
        'ZONE' $DOMTOT 'OPER' 'NS' 'Ro' 'UN' 'MUF' 'INCO' 'UN'
 'ZONE' $DOMTOT 'OPER' 'KEPSILON' 'Ro' 'UN' 'Mu' 'DT' 'INCO' 'KN' 'EN'
 OPTI EFM1 'CENTREE'
 'ZONE' $DOMTOT 'OPER' 'DFDT' 'Ro' 'UN' 'DT'            'INCO' 'UN'
         ;

RV = EQEX RV
        'CLIM' 'UN' 'UIMP' INLET Uinlet
               'UN' 'VIMP' INLET 0.
               'UN' 'VIMP' COTE  0.
        'CLIM' 'KN' 'TIMP' INLET Kinlet
               'EN' 'TIMP' INLET Einlet
 ;

RVP = EQEX 'OPTI' 'EF' KPRES
        ZONE $DOMTOT 'OPER' KBBT  (-1.)     INCO 'UN' 'PRES'
 CLIM 'PRES' TIMP poutlet 0.
;

rvp.'METHINV' . 'FCPRECT'=300 ;



RV.INCO = TABLE INCO ;
RV.INCO.'UN' = KCHT $DOMTOT VECT SOMMET (uinlet 0.) ;
RV.INCO.'KN' = KCHT $DOMTOT SCAL SOMMET Kinlet;
RV.INCO.'EN' = KCHT $DOMTOT SCAL SOMMET 1.e-3   ;
RV.INCO.'ENM1' = KCHT $DOMTOT SCAL SOMMET Einlet ;
RV.INCO.'KNM1' = KCHT $DOMTOT SCAL SOMMET Kinlet ;
RV.INCO.'MUF' = KCHT $DOMTOT SCAL SOMMET 1.E-10 ;
RV.INCO.'Ro' = KCHT $DOMTOT SCAL SOMMET 1. ;
RV.INCO.'Mu' = KCHT $DOMTOT SCAL SOMMET Nu ;
RV.INCO.'PRES' = KCHT $DOMTOT SCAL KPRES 0. ;
RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERK' = PROG 0.;
RV.INCO.'ERE' = PROG 0.;
RV.INCO.'EAN'=EAN;
RV.INCO.'KAN'=KAN;

 RV.'PROJ' = RVP ;

EVOLKA = EVOL (BLEU) 'CHPO' KAN 'SCAL' A1A2 ;
EVOLEA = EVOL (BLEU) 'CHPO' EAN 'SCAL' A1A2 ;

NCFL = DIME LCFL ;
REPETER BCFL NCFL;
CFL = extr LCFL &BCFL ;
dt = CFL*deltax/Uinlet ;
mess ' CFL=' CFL ' deltax ' deltax ' dt = ' dt ;
RV.INCO.'DT' = dt;

RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERK' = PROG 0.;
RV.INCO.'ERE' = PROG 0.;
RV.'PASDETPS'.'NUPASDT'= 1;

EXEC RV ;

EVOLKA = EVOLKA et (EVOL (ROUG) 'CHPO' RV.INCO.'KN' 'SCAL' A1A2) ;
EVOLEA = EVOLEA et (EVOL (ROUG) 'CHPO' RV.INCO.'EN' 'SCAL' A1A2) ;

EVOLK = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
           (RV.INCO.'ERK') ;
EVOLE = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
           (RV.INCO.'ERE') ;
Si (EGA &BCFL 1);
EVOLER = EVOLK et EVOLE;
Sinon ;
EVOLER = EVOLK et EVOLE et EVOLER ;
Finsi ;
* TESTS D'ERREUR *

erkm= extr &bcfl lerkm;
erem= extr &bcfl lerem;
list rv.inco.'ERK' ;
mess ' Sur K' (MINI RV.INCO.'ERK') erkm;
list rv.inco.'ERE' ;
mess ' Sur E' (MINI RV.INCO.'ERE') erem;
SI ( (MINI RV.INCO.'ERE') > erem ) ;
*        ERREUR 5 ;
FINSI ;
SI ( (MINI RV.INCO.'ERK') > erkm ) ;
*        ERREUR 5 ;
FINSI ;
FIN BCFL ;

SI ( GRAPH ) ;
TRACE DOMTOT ;
TAB1 = TABLE ;
TAB1.1 = 'TIRC MARQ TRIA REGU' ;
TAB1.2 = 'TIRC MARQ TRIA REGU' ;
TAB1.3 = 'TIRL MARQ LOSA REGU' ;
TAB1.4 = 'TIRL MARQ LOSA REGU' ;
TAB1.5 = 'TIRM MARQ TRIB REGU' ;
TAB1.6 = 'TIRM MARQ TRIB REGU' ;
TAB1.7 = 'TIRR MARQ ETOI REGU' ;
TAB1.8 = 'TIRR MARQ ETOI REGU' ;

TAB1.'TITRE' = TABLE ;
REPETER BCL NCFL ;
TAB1.'TITRE'. (2*&bcl-1) = (CHAI 'Cv_ k : CFL=' (EXTR LCFL &bcl));
TAB1.'TITRE'. (2*&bcl)   = (CHAI 'Cv_ e : CFL=' (EXTR LCFL &bcl));
FIN BCL ;

titl = chai 'Convergence Linf : Cfl = ' cfl ;
DESS EVOLER 'XBOR' 0. (2.*NITMA) 'YBOR' -15.0 0.0 'GRIL'
 'TITR' titl  LEGE TAB1 ;
        
TRACE RV.INCO.'KN' DOMTOT CNT 'TITR' 'KN' 20 ;
TRACE RV.INCO.'EN' DOMTOT CNT 'TITR' 'EN' 20 ;


TAB1 = TABLE ;
TAB1.1 = 'MARQ CROI ' ;
TAB1.2 = 'TIRC MARQ TRIA REGU' ;
TAB1.3 = 'TIRL MARQ CARR REGU' ;
TAB1.4 = 'TIRM MARQ LOSA REGU' ;
TAB1.5 = 'TIRM MARQ ETOI REGU' ;
TAB1.'TITRE' = TABLE ;
TAB1.'TITRE'. 1 = 'MOT' 'k : ANALYTIQUE' ;
REPETER BCL NCFL ;
TAB1.'TITRE'. (1+&bcl) = (CHAI 'k : CFL=' (EXTR LCFL &bcl));
FIN BCL ;

mess ' CFL=' CFL ;
titk = chai 'Grid Turbulence K : Cfl = ' cfl ;
tite = chai 'Grid Turbulence E : Cfl = ' cfl ;
DESS EVOLKA LEGE TAB1 'TITR' titk 'GRIL';

TAB1.'TITRE'. 1 = 'MOT' 'epsilon ANALYTIQUE' ;
REPETER BCL NCFL ;
TAB1.'TITRE'. (1+&bcl) = (CHAI 'e : CFL=' (EXTR LCFL &bcl));
FIN BCL ;
DESS EVOLEA LEGE TAB1 'TITR' tite 'GRIL';
FINSI ;

ERRK = KOPS RV.INCO.'KN' '-' KAN ;
ERRK = KOPS ERRK '*' ERRK ;
ERR2 = 0. ;
REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE)) ;
P1 = (doma $DOMTOT SOMMET) POIN &BLOC1 ;
ERR1 = 'EXTR' ERRK 'SCAL' P1 ;
ERR2 = ERR2 + ERR1 ;
FIN BLOC1 ;
ERR2 = ERR2/(nbno (doma $DOMTOT MAILLAGE)) ;
ERR2 = ERR2 ** 0.5 ;
ERR2K = ERR2/Kinlet ;
MESSAGE 'ERREUR SUR K EN NORME L2 = ' ERR2K ;

ERRE = KOPS RV.INCO.'EN' '-' EAN ;
ERRE = KOPS ERRE '*' ERRE ;
ERR2 = 0. ;
REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE))  ;
P1 = (DOMA $DOMTOT SOMMET) POIN &BLOC1 ;
ERR1 = 'EXTR' ERRE 'SCAL' P1 ;
ERR2 = ERR2 + ERR1 ;
FIN BLOC1 ;
ERR2 = ERR2/(NBNO (DOMA $DOMTOT MAILLAGE)) ;
ERR2 = ERR2 ** 0.5 ;
ERR2E = ERR2/Einlet ;
MESSAGE 'ERREUR SUR E EN NORME L2 = ' ERR2E ;


FIN ;