* fichier :  convnonlin1.dgibi
************************************************************************
* Section : Fluides Transitoire
************************************************************************

COMPLET = FAUX ;
GRAPH = FAUX ;
NX = 10 ;
NY = 10 ;
CFL = 0.7 ;
'SI' ( COMPLET ) ;
      TFIN = 10. ;
'SINON' ;
      TFIN = 1. ;
'FINSI' ;

***********************************************************
* EQUATION DE CONVECTION NON-LINEAIRE Ut + div (F(U)) = 0 *
* AVEC F(U) = U*U/2 1_x + ln U 1_y    U > 0               *
* RESOLUE SOUS FORME NON-CONSERVATIVE                     *
*       dU/dt + U dU/dx + 1/U dU/dy = 0                   *
* AVEC OPTION SUPGDC (par défaut)                         *
* SOLUTION ANALYTIQUE Uex(x,y,t) = (x-t)/(y-t)            *
* Probleme hyperbolique avec conditions aux limites       *
* variant dans le temps                                   *
* H. PAILLERE AVRIL 1997                                  *
***********************************************************

*******************************************************************
* PROCEDURE POUR CALCULER CHAMP DE VITESSE, METTRE A JOUR LES C.L *
* TESTER LA CONVERGENCE VERS L'ETAT STATIONNAIRE ET L'ERREUR      *
* ENTRE SOLUTION NUMERIQUE ET SOLUTION ANALYTIQUE                 *
*******************************************************************

'DEBPROC' CALCUL ;
'ARGU' RVX*'TABLE' ;

RV = RVX.'EQEX' ;
CN = RV.INCO.'CN' ;

* calcul du champ de vitesse 

VX1 = CN ;
VX1 = 'NOMC' 'UX' VX1 ;
VY1 = 'KOPS' 1. '/' CN ;
VY1 = 'NOMC' 'UY' VY1 ;
VV='KCHT' $DOMTOT 'VECT' 'SOMMET' 'COMP' (MOTS 'UX' 'UY')
        (VX1 'ET' VY1) ;
RV.INCO.'VITESSE' = VV ;

* mise a jour des conditions aux limites

CL_CN = 'EXCO' 'CN' RV.'CLIM' 'CN' ;
UEXAC = 'KOPS' ('KOPS' XX '-' RV.PASDETPS.'TPS') '/'
             ('KOPS' YY '-' RV.PASDETPS.'TPS') ;
CL_NOUV = 'KCHT' $CL 'SCAL' 'SOMMET' UEXAC ;
CL_NOUV = 'NOMC' 'CN' CL_NOUV 'NATU' 'DISCRET' ;

RV.'CLIM' = 0.D0 * CL_CN + CL_NOUV ;

* calcul de l'erreur de convergence vers l'etat stationnaire
* et de l'erreur entre solution numerique et solution analytique

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

ERR = 'KOPS' (RV.INCO.'CN') '-' (RV.INCO.'CNM1') ;
ELI = 'MAXI' ERR 'ABS' ;
ELI = ('LOG' (ELI + 1.0E-20))/(LOG 10.) ;
'MESSAGE' 'TEMPS = ' RV.PASDETPS.'TPS' 
        'ITER ' RV.PASDETPS.'NUPASDT' '   ERREUR LINF (dUdt) ' ELI ;
IT = 'PROG' RV.PASDETPS.'NUPASDT' ;
ER = 'PROG' ELI ;
RV.INCO.'IT' = (RV.INCO.'IT') 'ET' IT ;
RV.INCO.'ERCONV' = (RV.INCO.'ERCONV') 'ET' ER ;

UEXAC = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' UEXAC ;
UEXAC = 'NOMC' 'CN' UEXAC 'NATU' 'DISCRET' ;
ERRP = 'KOPS' UEXAC '-' RV.INCO.'CN' ;
ERRP = 'KOPS' ERRP '*' ERRP ;
ERR2 = 0. ;
'REPETER' BLOC1 (NBNO DOMTOT) ;
P1 = (DOMA $DOMTOT SOMMET) 'POIN' &BLOC1 ;
ERR1 = 'EXTR' ERRP 'SCAL' P1 ;
ERR2 = ERR2 + ERR1 ;
'FIN' BLOC1 ;
ERR2 = ERR2/(NBNO DOMTOT) ;
ERR2 = ERR2 ** 0.5 ;
ER = 'PROG' ERR2 ;
RV.INCO.'ERSOLU' = (RV.INCO.'ERSOLU') 'ET' ER ;
'MESSAGE' 'ERREUR |Uex(x,y,t) - Unum(x,y,t)| EN NORME L2 = '  ERR2 ;

'FINSI' ;

RV.INCO.'CNM1' = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'CN') ;

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


'OPTI' 'DIME' 2 ;
'OPTI' 'ELEM' 'QUA4' ;

*************
* DIFFUSION *
*************
DIF = 1.0E-10 ;

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

A1 = -2.0 -2.0 ;
A2 = -1.0 -2.0 ;
A3 = -1.0 -1.0 ;
A4 = -2.0 -1.0 ;

FBAS = A1 'DROI' NX A2 ;
FDRO = A2 'DROI' NY A3 ;
FHAU = A3 'DROI' NX A4 ;
FGAU = A4 'DROI' NY A1 ;

DOMTOT = 'DALL' FBAS FDRO FHAU FGAU 'PLAN' ;

************************************
* CREATION DU MODELE NAVIER_STOKES *
************************************

GEOCL = FBAS 'ET' FGAU ;
'ELIM' GEOCL 1.E-3 ;

MDOMTOT = CHAN DOMTOT QUAF ;
$DOMTOT = MODE MDOMTOT 'NAVIER_STOKES' LINE ; DOMA $DOMTOT 'IMPR' ;

MGEOCL = CHAN GEOCL QUAF ;
ELIM MDOMTOT MGEOCL 1.E-3 ;
$CL     = MODE MGEOCL 'NAVIER_STOKES' LINE ;

*********************
* SOLUTION INITIALE *
*********************

XX YY = 'COOR' (DOMA $DOMTOT MAILLAGE) ;
SOLINI = 'KOPS' XX '/' YY ;
SOLCL = 'KCHT' $CL 'SCAL' 'SOMMET' SOLINI ;

***********************
* TABLE DE RESOLUTION *
***********************

RV = 'EQEX' $DOMTOT 'ITMA' 1000 'ALFA' CFL 'TFINAL' (TFIN/2.)
        'ZONE' $DOMTOT 
        'OPER' CALCUL 
        'OPTI' 'SUPGDC'
        'OPER' TSCAL DIF 'VITESSE' 0.0 'INCO' 'CN'
        'OPTI' 'CENTREE'
        'OPER' DFDT  1. 'CN' 'DELTAT'       'INCO' 'CN'
        'CLIM' 'CN' 'TIMP' GEOCL SOLCL ;

RV.INCO = 'TABLE' 'INCO' ;
RV.INCO.'CN' = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' SOLINI ;
RV.INCO.'VITESSE'= CHVIT ;

RV.INCO.'CNM1' = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' SOLINI ;
RV.INCO.'IT' = 'PROG' 1 ;
RV.INCO.'ERCONV' = 'PROG' 0. ;
RV.INCO.'ERSOLU' = 'PROG' 0. ;

* ON EXECTUTE JUSQU'A TFIN/2. *

EXEC RV ;

* ON TESTE L'OPTION SUPG (IDCEN=3) DE TFIN/2. A TFIN *

RV.'2TSCAL'.'KOPT'.'IDCEN' = 3 ;
RV.'TFINAL' = TFIN ;
EXEC RV ;

'SI' ( GRAPH ) ;
'OPTI' 'ISOV' 'SULI' ;
'TRACE' RV.INCO.'CN' DOMTOT ('CONT' DOMTOT) 20 ;
EVOL4 = 'EVOL' 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
           (RV.INCO.'ERCONV') ;
'DESS' EVOL4 'XBOR' 0. 2000. 'YBOR' -20.0 0.0 ;
'FINSI' ;

'SI' ( 'MINI' ( RV.INCO.'ERCONV' ) > -2.0 ) ;
        'ERREUR' 5 ;
'FINSI' ;

'SI' ( 'MAXI' ( RV.INCO.'ERSOLU' ) > 0.09 ) ;
        'ERREUR' 5 ;
'FINSI' ;

'FIN' ;