* fichier :  recircpp.dgibi
* tests ajustes par PV le 12/10/08 sur les valeurs obtenues en linux32 
opti dime 2 elem QUA8 ;

GRAPH=FAUX ;
*GRAPH=VRAI ;
nbit = 5 ;

*
* 1/  Exemple test : canal courbe
*                    Eclt laminaire incompressible
*                    On impose une vitesse a un bout
*                    V normale nulle sur les parois
* On indique comment imposer une Vitesse en entree d'un domaine
* (Vitesse normale non nulle et Vitesse tangentielle nulle !!)
* et comment imposer un glissement sur une paroi (Vitesse normale nulle)
*     Post  : debit sur le contour , Vitesses ,pression
*     Teste : vitesse normale imposee en 2D plan
*             algorithme de projection transitoire
*

DISCR='QUAF' ;
KPRES='MSOMMET' ;
R=4. ; n1=5 ; n2=20 ;
*R=4. ; n1=2 ; n2=5  ;
p1=0 0 ; p2= 1 0 ; pc= R 0 ;
entree = p1 d n1 p2 ;
sortie = entree tour pc (-60.) ;
sortie = inve sortie;
q2 = poin 1 sortie ;
q1 = sortie poin final ;
pp1=p1 c n2 pc q1 ;
pp2=p2 c n2 pc q2 ;
paroi= pp1 et pp2 ;
paroi = inve paroi ;

cnt= entree et pp2 et sortie et (inve pp1) ;
mt= surf cnt;
Mmt= chan quaf mt ;
Mcnt= chan quaf cnt ;
Mentree= chan quaf entree;
Msortie= chan quaf sortie;
Mparoi = chan quaf paroi ;

elim (Mmt et Mcnt et Mentree et Msortie et Mparoi) 1.e-5 ;
$mt = mode Mmt 'NAVIER_STOKES' DISCR ;
mt= doma $mt maillage ;
$entree= mode Mentree 'NAVIER_STOKES' DISCR ;
$sortie= mode Msortie 'NAVIER_STOKES' DISCR ;
$paroi = mode Mparoi  'NAVIER_STOKES' DISCR ;
entrep=doma $entree 'MSOMMET' ;
sortie=doma $sortie maillage ;

NU= 1.e-2 ;
U0=1. ;
DT=2. ;

* Comment imposer un debit d'entree
* En general il faut imposer la composante normale a
* la vitesse debitante et les composantes tangentielles a ZERO !
* On n'utilise donc pas VNIMP :

* On extrait le champ des normales de la surface d'entree
* en prenant garde a l'orientation des normales
nsor=doma $sortie 'NORMALEV';
* On extrait les composantes nx ny etc
nx= exco 'UX' nsor ;
ny= exco 'UY' nsor ;
* on cree les componsantes cartesiennes du champ de vitesse
* qui seront imposees en condition limite type Dirichlet
usx=u0*nx ;
usy=u0*ny ;

* Sur les parois la condition de glissement s'obtient
* effectivement en imposant une vitesse normale nulle a l'aide de VNIMP


rv = eqex 'ITMA' 1  'NITER' 1  'OMEGA' 1.
  OPTI EF 'CENTREE' 'IMPL' KPRES
  'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  'CLIM'
  'UN' UIMP sortie  usx  'UN' VIMP sortie usy
  ;

rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;

    RVP= EQEX 'OPTI' 'EF' KPRES
   'ZONE' $mt 'OPER' 'KBBT'  -1.  'INCO' 'UN' 'PRES'
   'ZONE' $paroi 'OPER' 'VNIMP' $mt 0. 'INCO' 'UN' 'PRES'
   'CLIM'
   'PRES' 'TIMP' entrep 0. ;
    ;

*  rvp. 'METHINV' . 'FCPRECT'=100 ;

    RV.'PROJ'= RVP ;


 exec rv ;

* affichage vitesse
un=rv.inco.'UN' ;
ung=vect un 0.1 ux uy jaune ;

Si GRAPH ;trace ung mt;Finsi ;
cnt1=chan cnt seg2 ;                                                    
Mcnt1= chan cnt1 quaf;                                                  
$cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
q=dbit (rv.inco.'UN') $cnt ;                                            
Mess ' Debit sur le contour ' q ;
 dq= (abs q) - 3.36856E-02 ;
 list dq ;
 si (dq  > 1.e-5 ) ; erreur 5 ; finsi ;

* affichage pression
* On extrait la bonne composante de la pression
* il y a aussi les multiplicateurs de Lagrange vitesse normale
pn=exco rv.inco.'PRESSION' 'PRES' ;

* On reaffecte les bon SPG au CHPOINT
pn= kcht $mt scal KPRES pn ;
mtp= doma $mt 'MMAIL' ;

Si GRAPH ;trace pn mtp;Finsi ;
mpn=(maxi pn);
dpn= abs (mpn - 1.7705);
mess ' MAXI PN = ' mpn 'dpn=' dpn ;
 si ( dpn > 1.e-4 ) ; erreur 5 ; finsi ;



*
*  2/ Exemple test : recirculation dans une cavite semi circulaire
*                    Eclt laminaire incompressible
*                    On impose une vitesse tg sur la surface libre
*                    V normale nulle sur le fond
*     Voir RECIRC.DGIBI pour l'algorithme Taylor Hood
*     Post  : debit sur le contour , fonction de courant , pression
*     Teste : vitesse normale imposee en 2D plan
*             algorithme de projection transitoire
*

p0=0. 0. ;
p1= -1. 0.;
p2= 1. 0. ;
p03= 0. -1. ;
p3= 1. -1. ;
p4= -1. -1. ;

n1=  5; n2= 5;
ab= p1 d n1 p2 ;
slib=  (chan ab POI1) ;
nab= nbel slib ;
slib= elem slib (lect 2 pas 1 (nab - 1));

bc= p2 c p0 n2 p03 c  p0 n2 p1;
*bc= p2 d n2 p3 d n2 p4 d n2 p1 ;

cnt=ab et bc ;
opti elem TRI6 ;
mt= cnt surf ;
Mmt=chan mt quaf ;
Mbc=chan bc quaf ;
elim (Mmt et Mbc) 1.e-5 ;

$mt = mode Mmt 'NAVIER_STOKES' DISCR;
$bc = mode Mbc 'NAVIER_STOKES' DISCR;

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

*NU= 1. / 61. ;
NU= 1.e-2 ;
U0=1. ;
DT=1. ;

rv = eqex 'ITMA' 10 'NITER' 1  'OMEGA' 1.
  OPTI EF 'CENTREE' 'IMPL' KPRES
  'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  'CLIM'
         'UN' UIMP slib 1. 'UN' VIMP ab 0.
  ;

rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
rv.'INCO'.'VN'=kcht $bc scal sommet 0. ;

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

*  rvp. 'METHINV' . 'FCPRECT'=100 ;

    RV.'PROJ'= RVP ;


exec rv ;

un = rv.inco.'UN' ;
ung= vect un 0.11  ux uy jaune ;

Si GRAPH ;trace ung Mmt ; Finsi ;

cnt1=chan cnt seg2 ;                                                    
Mcnt1= chan cnt1 quaf;                                                  
$cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
q=dbit (rv.inco.'UN') $cnt ;                                            
Mess ' Debit sur le contour ' q ;
 dq= (q - .22362) abs; list dq;
 si ( dq > 1.e-5 ) ; erreur 5 ; finsi ;

*
* Calcul de la fonction de courant
*

sw = 'KOPS' un 'ROT' $mt ;
rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
     'ZONE' $mt 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI'
     'ZONE' $mt 'OPER' 'FIMP' sw      'INCO' 'PSI'
     'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;

opti isov SULI ;
exec rk ;
psi=rk.inco.'PSI' ;
Si GRAPH ;trace psi Mmt cnt ; Finsi ;
mpsi=(mini psi);
dpsi= abs (mpsi + .17603) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
 si ( dpsi > 5.e-4 ) ; erreur 5 ; finsi ;

*
* Calcul de la pression
*

  pn=kcht $mt scal KPRES (exco rv.inco.'PRESSION' 'PRES') ;
  mtp=doma $mt 'MMAIL' ;
Si GRAPH ;  trace pn mtp ; Finsi ;
mpn=(maxi pn);
dpn= abs (mpn - .38065);
mess ' MAXI PN = ' mpn 'dpn=' dpn ;
 si ( dpn > 1.e-3 ) ; erreur 5 ; finsi ;


*
*  Troisieme cas : Meme probleme geometrie rectangulaire
*                  MACRO CENTRE
*
*
*
*

DISCR = 'QUAF';
KPRES= 'MSOMMET' ;

p0=0. 0. ;
p1= -1. 0.;
p2= 1. 0. ;
p03= 0. -1. ;
p3= 1. -1. ;
p4= -1. -1. ;

n1=  5 ; n2= 5 ;
ab= p1 d n1 p2 ;
slib=  (chan ab POI1) ;
nab= nbel slib ;
slib= elem slib (lect 2 pas 1 (nab - 1));

bc= p2 d n2 p3 d n2 p4 d n2 p1 ;

cnt=ab et bc ;

mt= cnt surf ;
Mmt=chan mt quaf ;
Mbc=chan bc quaf ;
elim (Mmt et Mbc) 1.e-5 ;

$mt = mode Mmt 'NAVIER_STOKES' DISCR;
$bc = mode Mbc 'NAVIER_STOKES' DISCR;

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

NU= 1.e-2 ;
U0=1. ;
DT=1. ;

rv = eqex 'ITMA' 10 'NITER' 1  'OMEGA' 1.
  OPTI EF 'CENTREE' 'IMPL' KPRES
  'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  'CLIM'
         'UN' UIMP slib 1. 'UN' VIMP ab 0.
  ;

rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
rv.'INCO'.'VN'=0. ;

    RVP= EQEX 'OPTI' 'EF' KPRES
   'ZONE' $mt 'OPER' 'KBBT'  -1.  'INCO' 'UN' 'PRES'
   'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
    CLIM PRES TIMP bcp 0. ;
    RV.'PROJ'= RVP ;

exec rv ;


un = rv.inco.'UN' ;
ung= vect un 0.11  ux uy jaune ;

Si GRAPH ;trace ung Mmt ; Finsi ;

cnt1=chan cnt seg2 ;                                                    
Mcnt1= chan cnt1 quaf;                                                  
$cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
q=dbit (rv.inco.'UN') $cnt ;                                            
Mess ' Debit sur le contour ' q ;
 dq=abs ( q - 2.78540E-03);
 si ( dq  > 1.e-3 ) ; erreur 5 ; finsi ;

sw = 'KOPS' un 'ROT' $mt ;
rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
     'ZONE' $mt 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI'
     'ZONE' $mt 'OPER' 'FIMP' sw      'INCO' 'PSI'
     'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;

opti isov SULI ;
exec rk ;
psi=rk.inco.'PSI' ;
mess ' MINI PSI = ' (mini psi) ;
Si GRAPH ;trace psi Mmt cnt ; Finsi ;
mpsi=(mini psi);
dpsi= abs (mpsi + .18561) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
 si ( dpsi > 1.e-3 ) ; erreur 5 ; finsi ;


*
* Calcul de la pression
*
  pn=kcht $mt scal KPRES (exco rv.inco.'PRESSION' 'PRES') ;
  mtp=doma $mt 'MMAIL';
Si GRAPH ;  trace pn mtp ; Finsi ;
mpn=(maxi pn);
dpn= abs (mpn - .43612);
mess ' MAXI PN = ' mpn 'dpn=' dpn ;
 si ( dpn > 0.1   ) ; erreur 5 ; finsi ;

FIN ;