* fichier :  gridturb_expl.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                            *
**********************************************************

COMPLET = FAUX ;
GRAPH = FAUX ;

SI ( COMPLET ) ;
        D = 4 ;
        NITER = 1000 ;
SINON ;
        D = 1 ;
        NITER = 150 ;
FINSI ;

NX = D*10 ;
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.1 ;
A4 = 0.0 0.1 ;

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 ;
ELIM (DOMTOT et A4A1 et A2A3) 1.E-3 ;

$DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' 'MACRO' ;
DOMA $DOMTOT IMPR ;
$INLET  = 'MODE' A4A1 'NAVIER_STOKES' 'MACRO' ;
 INLET  = 'DOMA' $INLET 'MAILLAGE' ;
$OUTLET = 'MODE' A2A3 'NAVIER_STOKES' 'MACRO' ;

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

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

Uinlet = 1.0 ;
Kinlet = 0.01 ;
Einlet = 0.01 ;

Nu = 1.E-10 ;
tos = 0. 0. ;

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

C2 = 1.92 ;
X1 = COOR 1 (DOMA $DOMTOT SOMMET) ;
ZZZ = KOPS X1 '*' (EINLET/KINLET*(C2-1.0)) ;
ZZZ = KOPS ZZZ '+' 1.0 ;
KAN = KOPS ZZZ '**' (1./(1.-C2)) ;
KAN = KOPS KAN '*' KINLET ;

EAN = KOPS ZZZ '**' (C2/(1.-C2)) ;
EAN = KOPS EAN '*' EINLET ;

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

RV = EQEX $DOMTOT 'ITMA' NITER 'ALFA' 0.8
        'ZONE' $DOMTOT 'OPER' CALCUL
        'OPTI' 'SUPGDC'
        'ZONE' $DOMTOT 'OPER' 'NSKE' NU 'NUT' 'INCO' 'UN' 'KN' 'EN' 
        'OPTI' 'CENTREE'
        'ZONE' $OUTLET 'OPER' 'TOIM' tos 'INCO' 'UN'
        'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
        'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
        'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN';

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

RVP = EQPR $DOMTOT KTYPI 1
        ZONE $DOMTOT 'OPER' 'PRESSION' 0. ;

RV.INCO = TABLE INCO ;
RV.INCO.'UN' = KCHT $DOMTOT VECT SOMMET (1. 0.) ;
RV.INCO.'KN' = KCHT $DOMTOT SCAL SOMMET Kinlet ;
RV.INCO.'EN' = KCHT $DOMTOT SCAL SOMMET Einlet ;
RV.INCO.'ENM1' = KCHT $DOMTOT SCAL SOMMET Einlet ;
RV.INCO.'KNM1' = KCHT $DOMTOT SCAL SOMMET Kinlet ;
RV.INCO.'NUT' = KCHT $DOMTOT SCAL CENTRE 1.E-10 ;
RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERK' = PROG 0.;
RV.INCO.'ERE' = PROG 0.;
RV.PRESSION = RVP ;

EXEC RV ;

SI ( GRAPH ) ;
TRACE DOMTOT ;
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') ;
TAB1 = TABLE ;
TAB1.1 = 'MARQ CROI ' ;
TAB1.2 = 'MARQ TRIA ' ;
TAB1.'TITRE' = TABLE ;
TAB1.'TITRE'. 1 = MOT 'CONVERGENCE_K' ;
TAB1.'TITRE'. 2 = MOT 'CONVERGENCE_E' ;
DESS (EVOLK ET EVOLE) 'XBOR' 0. (2.*NITER) 'YBOR' -15.0 0.0 
                LEGE TAB1 ;
        
TRACE RV.INCO.'KN' DOMTOT CNT 'TITR' 'KN' 20 ;
TRACE RV.INCO.'EN' DOMTOT CNT 'TITR' 'EN' 20 ;

EVOLKN = EVOL (ROUG) 'CHPO' RV.INCO.'KN' 'SCAL' A1A2 ;
EVOLKA = EVOL (BLEU) 'CHPO' KAN 'SCAL' A1A2 ;
EVOLEN = EVOL (ROUG) 'CHPO' RV.INCO.'EN' 'SCAL' A1A2 ;
EVOLEA = EVOL (BLEU) 'CHPO' EAN 'SCAL' A1A2 ;

TAB1 = TABLE ;
TAB1.1 = 'MARQ CROI ' ;
TAB1.2 = ' ' ;
TAB1.'TITRE' = TABLE ;
TAB1.'TITRE'. 1 = MOT 'NUMERIQUE' ;
TAB1.'TITRE'. 2 = 'MOT' 'ANALYTIQUE' ;

DESS (EVOLKN ET EVOLKA) LEGE TAB1 'TITR' 'Grid Turbulence K' ;
DESS (EVOLEN ET EVOLEA) LEGE TAB1 'TITR' 'Grid Turbulence E' ;
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 ;

* TESTS D'ERREUR *

SI ( (MINI RV.INCO.'ERE') > -5.0 ) ;
        ERREUR 5 ;
FINSI ;
SI ( (MINI RV.INCO.'ERK') > -5.0 ) ;
        ERREUR 5 ;
FINSI ;

SI ( ERR2K > 1.5D-2 ) ;
        ERREUR 5 ;
FINSI ;

SI ( ERR2E > 3.0D-2 ) ;
        ERREUR 5 ;
FINSI ;

FIN ;