* fichier :  konv_scal_impl3d.dgibi
************************************************************************
* Section : Fluides Darcy
************************************************************************
***********************************************************
***********************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour le    ****
**** transport des scalaires                           **** 
**** OPERATEURS PRET, KONV                             ****
**** Implicit: calcul du jacobien du residu            ****
****                                                   ****
**** Methodes: UPWIND, CENTERED                        ****
****                                                   ****
**** A. BECCANTINI SFME/LMTF            DECEMBRE 2001  ****
***********************************************************
***********************************************************

 'OPTION'  'DIME' 3  'ELEM' CUB8  'ECHO'  0  'TRAC' 'X' ;


*****************************************************
*****************************************************
******** PROCEDURES *********************************
*****************************************************
*****************************************************

*
* Derivé partielle du residu en un point par rapport
* aux variable en un autre point
*

 'DEBPROC' JACOVA JACO*'MATRIK' $MODE*'MMODEL' LISTINCO*'LISTMOTS'
    PPRIM*'POINT' PDUAL*'POINT' MOTPRI*'MOT' MOTDUA*'MOT';
    
*   PPRIM = point ou est localisé la variable primale
*   PDUAL = point ou est localisé la variable duale
*   MOTPRI = nom de la composante concernante la variable primale
*   MOTDUA = nom de la composante concernante la variable duale

    ELT1 = 'MANUEL' 'POI1' PPRIM ;
    NDIM = 'DIME' LISTINCO ;
    CHPUN = 'MANUEL' 'CHPO' ELT1 1 MOTPRI 1.0 'NATURE' 'DISCRET' ;
    'REPETER' BL1 NDIM ;
       MOTCEL = 'EXTRAIRE' LISTINCO &BL1 ;
       CHPUN = CHPUN 'ET' ('MANUEL' 'CHPO' ('DOMA' $MODE 'CENTRE')
          1 MOTCEL 0.0 'NATURE' 'DISCRET') ;
    'FIN' BL1 ;
    D_DMOT = 'KOPS' JACO 'MULT' CHPUN ;
    SCAL = 'EXTRAIRE' D_DMOT PDUAL MOTDUA ;

 'FINPROC' SCAL  ;

 
 'DEBPROC' JACNUM $MODE*'MMODEL' NOMMET*'MOT' SN*'CHPOINT' UN*'CHPOINT'
    PPRIM*'POINT' PDUAL*'POINT'
    MOTPRI*'MOT' MOTDUA*'MOT' EPSILON*'FLOTTANT';

*   PPRIM = point ou est localisé la variable primale
*   PDUAL = point ou est localisé la variable duale
*   MOTPRI = nom de la composante concernante la variable primale
*   MOTDUA = nom de la composante concernante la variable duale

*   Le valeur dans l'état non-perturbé en PDUAL ;
    
    SNCEL   = 'COPIER'  SN ;
    
    SF =  'PRET' 'CLAUDEIS'  'FACE' 1  $MODE SN ;
    
    CHPRES0 DT = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'RESI' NOMMET
             $MODE SF UN ;

    VAL0 = 'EXTRAIRE' CHPRES0 PDUAL MOTDUA ;
    
*   On etabli la variable à perturber

    LISTINCO = 'EXTRAIRE' SN 'COMP' ;

    NDIM = 'DIME' LISTINCO ;
    'REPETER' BL1 NDIM ;
       MOTCEL = 'EXTRAIRE' LISTINCO &BL1 ;
       'SI' ('EGA' MOTCEL MOTPRI) ;
          ICEL = &BL1 ;
          'QUITTER' BL1 ;
       'FINSI' ;
    'FIN' BL1 ;
    
    'SI' (ICEL > NDIM) ;
       'MESSAGE' 'Procedure JACNUM' ;
       'MESSAGE' 'MOTPRI = ??? ';
       'ERREUR' 21 ;
    'FINSI' ;
    
    ELT1 = 'MANUEL' 'POI1' PPRIM ;
    
    SNCEL = ('MANUEL' 'CHPO' ELT1 1 MOTCEL EPSILON)
      '+' SN  ;

    SF =  'PRET' 'CLAUDEIS' 'FACE' 1  $MODE SNCEL ;
    
    CHPRES1 DT = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'RESI' NOMMET
             $MODE SF UN ;
      
    VAL1 = 'EXTRAIRE' CHPRES1 PDUAL MOTDUA ;

 'FINPROC' ((VAL1 '-' VAL0) '/' EPSILON)  ;

*****************************************************
*****************************************************
******** FIN PROCEDURES *****************************
*****************************************************
*****************************************************
*****************************************************

*
*** GRAPH
*

 GRAPH = FAUX ;
* GRAPH = VRAI ;

 ERRTOL = 1.0D-6 ;

***************************
***** DOMAINE SPATIAL  ****
***************************

***************************
***** DOMAINE SPATIAL  ****
***************************


 A0 = 0.0D0 0.0D0 0.0D0;
 A1 = 1.0D0 0.0D0 0.0D0;
 A2 = 1.0D0 1.0D0 0.0D0;
 A3 = 0.0D0 1.0D0 0.0D0;

 SUR1 = 'MANUEL' 'QUA4' A0 A1 A2 A3 ;
 DOM1 = SUR1 'VOLUME' 'TRAN' 1 (0.0 0.0 1.0) ;
 
 DOM2 = DOM1 'PLUS' (1.0 0.0 0.0) ;
 DOM3 = DOM1 'PLUS' (-1.0 0.0 0.0) ;
 DOM4 = DOM1 'PLUS' (0.0 1.0 0.0) ;
 DOM5 = DOM1 'PLUS' (0.0 -1.0 0.0) ;
 DOM6 = DOM1 'PLUS' (1.0 1.0 0.0) ;
 DOM7 = DOM1 'PLUS' (-1.0 -1.0 0.0) ;
 DOM8 = DOM1 'PLUS' (1.0 -1.0 0.0) ;
 DOM9 = DOM1 'PLUS' (-1.0 1.0 0.0) ;

 DOM10 = DOM1 'PLUS' (0.0 0.0 1.0) ;
 DOM11 = DOM1 'PLUS' (1.0 0.0 1.0) ;
 DOM12 = DOM1 'PLUS' (-1.0 0.0 1.0) ;
 DOM13 = DOM1 'PLUS' (0.0 1.0 1.0) ;
 DOM14 = DOM1 'PLUS' (0.0 -1.0 1.0) ;
 DOM15 = DOM1 'PLUS' (1.0 1.0 1.0) ;
 DOM16 = DOM1 'PLUS' (-1.0 -1.0 1.0) ;
 DOM17 = DOM1 'PLUS' (1.0 -1.0 1.0) ;
 DOM18 = DOM1 'PLUS' (-1.0 1.0 1.0) ;

 DOM19 = DOM1 'PLUS' (0.0 0.0 -1.0) ;
 DOM20 = DOM1 'PLUS' (1.0 0.0 -1.0) ;
 DOM21 = DOM1 'PLUS' (-1.0 0.0 -1.0) ;
 DOM22 = DOM1 'PLUS' (0.0 1.0 -1.0) ;
 DOM23 = DOM1 'PLUS' (0.0 -1.0 -1.0) ;
 DOM24 = DOM1 'PLUS' (1.0 1.0 -1.0) ;
 DOM25 = DOM1 'PLUS' (-1.0 -1.0 -1.0) ;
 DOM26 = DOM1 'PLUS' (1.0 -1.0 -1.0) ;
 DOM27 = DOM1 'PLUS' (-1.0 1.0 -1.0) ;

 DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 'ET' DOM4 'ET' DOM5
   'ET' DOM6 'ET' DOM7 'ET' DOM8 'ET' DOM9 'ET' DOM10
   'ET' DOM11 'ET' DOM12 'ET' DOM13 'ET' DOM14 'ET' DOM15
   'ET' DOM16 'ET' DOM17 'ET' DOM18 'ET' DOM19 'ET' DOM20
   'ET' DOM21 'ET' DOM22 'ET' DOM23 'ET' DOM24 'ET' DOM25
   'ET' DOM26 'ET' DOM27 ;

 'ELIMINATION' DOMTOT 0.0001 ;

*
**** Perturbation du domaine
*

 CHPBRU = ('NOMC' 'UX' ('BRUI' 'BLAN' 'UNIF' 0.0  0.4 DOMTOT)
            'NATU' 'DISCRET') 'ET'
            ('NOMC' 'UY' ('BRUI' 'BLAN' 'UNIF' 0.0  0.2 DOMTOT)
            'NATU' 'DISCRET') 'ET'
            ('NOMC' 'UZ' ('BRUI' 'BLAN' 'UNIF' 0.0  0.3 DOMTOT)
            'NATU' 'DISCRET');
            
 'FORME' CHPBRU ;

 $DOMTOT = 'MODELISER' DOMTOT 'EULER';
 $DOM6 = 'MODELISER'  DOM6 'EULER';
 $DOM9 = 'MODELISER'  DOM9 'EULER';
 TDOMTOT = 'DOMA'  $DOMTOT 'VF';
 TDOM6 = 'DOMA'  $DOM6 'VF';
 TDOM9 = 'DOMA'  $DOM9 'VF';
 MDOM6 = TDOM6 . 'QUAF' ;
 MDOM9 = TDOM9 . 'QUAF' ;
 MDOMTOT = TDOMTOT . 'QUAF' ;

 'ELIM' (MDOMTOT 'ET' MDOM6 'ET' MDOM9) 1.E-6 ;

 'SI' GRAPH;
      'TRACER' (('DOMA' $DOMTOT 'MAILLAGE' ) 'ET' 
      ('DOMA' $DOMTOT 'CENTRE')) 'TITRE' 'Domaine et centre' ;
 'FINSI' ;

*
***** Densité, pression, gamma ***********
*

 RN = ('BRUI' 'BLAN' 'UNIF' 1.11 0.5 ('DOMA' $DOMTOT 'CENTRE'))
       'NOMC' 'RN' 'NATU' 'DISCRET' ;
 PN = ('BRUI' 'BLAN' 'UNIF' 1234.1 800 ('DOMA' $DOMTOT 'CENTRE'))
       'NOMC' 'PN' 'NATU' 'DISCRET' ;
 GAMMAN = ('BRUI' 'BLAN' 'UNIF' 14.1 800 ('DOMA' $DOMTOT 'CENTRE'))
       'NOMC' 'GAMN' 'NATU' 'DISCRET' ;

 SN = RN 'ET' PN 'ET' GAMMAN ;

*
**** Vitesse aux faces
*

 UN = (('BRUI' 'BLAN' 'UNIF' 0.0 1.000 ('DOMA' $DOMTOT 'FACE'))
       'NOMC' 'UX' 'NATU' 'DISCRET') 'ET'
      (('BRUI' 'BLAN' 'UNIF' 0.0 1.0000 ('DOMA' $DOMTOT 'FACE'))
       'NOMC' 'UY' 'NATU' 'DISCRET') 'ET'
      (('BRUI' 'BLAN' 'UNIF' 0.0 1.0000 ('DOMA' $DOMTOT 'FACE'))
       'NOMC' 'UZ' 'NATU' 'DISCRET') ;
       
****************************************************
****************************************************
******* Calcul du jacobien et du residu   **********
****************************************************
****************************************************
*
* JACO est le jacobien
*
*


 'REPETER' BLMETO 2 ;

 'SI' (&BLMETO 'EGA' 1) ;
    METO = 'UPWIND' ;
 'FINSI' ;
 
 'SI' (&BLMETO 'EGA' 2) ;
    METO = 'CENTERED' ;
 'FINSI' ;
 
 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'METHODE = ' METO) ;
 'MESSAGE' ;
 
 SF =  'PRET' 'CLAUDEIS' 'FACE' 1  $DOMTOT SN ;

 DELTAS DT = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'RESI' METO
             $DOMTOT SF UN ;
             
 JACO =  'KONV' 'VF' 'CLAUDEIS' 'FACE' 'JACO' METO
             $DOMTOT SF UN ;


*****************************************************
******* TEST1 ***************************************
*****************************************************
*
* On compare le jacobien et la variation du residu
* en $DOM9 'CENTRE' par rapport à une variation
* infinitésimal en $DOM9 'CENTRE'
*

 PCEN9 = 'POIN' 1 ('DOMA' $DOM9 'CENTRE') ;

*
* Le jacobien exact. 
* DSiSj  = d(RES_Si)/dSj (variable primale en PCEN9, variable duale en PCEN9) ;
* ...

 LISTINCO = 'MOTS' 'RN' 'PN' 'GAMN' ;
 
 DS1S1   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN9 'RN' 'RN' ;
 DS1S3   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN9 'GAMN' 'RN' ;
 DS2S2   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN9 'PN' 'PN' ;
 DS3S3   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN9 'GAMN' 'GAMN' ;

 'SI' (('ABS' DS1S3) > ERRTOL) ;
    'MESSAGE' 'JACO CLAUDEIS = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
* Le jacobien numerique 


 DELTA = 1.0D-3 ;

 DS1S1N = JACNUM $DOMTOT METO SN UN  PCEN9 PCEN9
   'RN' 'RN' DELTA ;
 DS2S2N = JACNUM $DOMTOT METO SN UN  PCEN9 PCEN9
   'PN' 'PN' DELTA ;
 DS3S3N = JACNUM $DOMTOT METO SN UN  PCEN9 PCEN9
   'GAMN' 'GAMN' DELTA ;

 ERRO = 'MAXIMUM' ('PROG' (DS1S1 '-' DS1S1N)  (DS2S2 '-' DS2S2N)
    (DS3S3 '-' DS3S3N)) 'ABS' ;

 'SI' (ERRO > ERRTOL) ;
    'MESSAGE' 'JACO CLAUDEIS = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;

*****************************************************
*****************************************************
******* TEST2 ***************************************
*****************************************************
*****************************************************
*
* On compare le jacobien et la variation du residu
* en $DOM9 'CENTRE' par rapport à une variation
* infinitésimal en $DOM6 'CENTRE'
*

 PCEN6 = ('DOMA' $DOM6 'CENTRE') 'POIN' 1 ;
 
*
* Le jacobien exact. 
*

 DS1S1   = JACOVA JACO $DOMTOT LISTINCO PCEN6 PCEN9 'RN' 'RN' ;
 DS1S3   = JACOVA JACO $DOMTOT LISTINCO PCEN6 PCEN9 'GAMN' 'RN' ;
 DS2S2   = JACOVA JACO $DOMTOT LISTINCO PCEN6 PCEN9 'PN' 'PN' ;
 DS3S3   = JACOVA JACO $DOMTOT LISTINCO PCEN6 PCEN9 'GAMN' 'GAMN' ;

 'SI' (('ABS' DS1S3) > ERRTOL) ;
    'MESSAGE' 'JACO CLAUDEIS = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
* Le jacobien numerique 


 DELTA = 1.0D-3 ;

 DS1S1N = JACNUM $DOMTOT METO SN UN  PCEN6 PCEN9
   'RN' 'RN' DELTA ;
 DS2S2N = JACNUM $DOMTOT METO SN UN  PCEN6 PCEN9
   'PN' 'PN' DELTA ;
 DS3S3N = JACNUM $DOMTOT METO SN UN  PCEN6 PCEN9
   'GAMN' 'GAMN' DELTA ;

 ERRO = 'MAXIMUM' ('PROG' (DS1S1 '-' DS1S1N)  (DS2S2 '-' DS2S2N)
    (DS3S3 '-' DS3S3N)) 'ABS' ;

 'SI' (ERRO > ERRTOL) ;
    'MESSAGE' 'JACO CLAUDEIS = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;


*****************************************************
*****************************************************
******* TEST2 ***************************************
*****************************************************
*****************************************************
*
* On compare le jacobien et la variation du residu
* en $DOM6 'CENTRE' par rapport à une variation
* infinitésimal en $DOM9 'CENTRE'
*

*
* Le jacobien exact. 
*

 DS1S1   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN6 'RN' 'RN' ;
 DS1S3   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN6 'GAMN' 'RN' ;
 DS2S2   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN6 'PN' 'PN' ;
 DS3S3   = JACOVA JACO $DOMTOT LISTINCO PCEN9 PCEN6 'GAMN' 'GAMN' ;

 'SI' (('ABS' DS1S3) > ERRTOL) ;
    'MESSAGE' 'JACO CLAUDEIS = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
* Le jacobien numerique 


 DELTA = 1.0D-3 ;

 DS1S1N = JACNUM $DOMTOT METO SN UN  PCEN9 PCEN6
   'RN' 'RN' DELTA ;
 DS2S2N = JACNUM $DOMTOT METO SN UN  PCEN9 PCEN6
   'PN' 'PN' DELTA ;
 DS3S3N = JACNUM $DOMTOT METO SN UN  PCEN9 PCEN6
   'GAMN' 'GAMN' DELTA ;

 ERRO = 'MAXIMUM' ('PROG' (DS1S1 '-' DS1S1N)  (DS2S2 '-' DS2S2N)
    (DS3S3 '-' DS3S3N)) 'ABS' ;

 'SI' (ERRO > ERRTOL) ;
    'MESSAGE' 'JACO CLAUDEIS = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;

 'FIN' BLMETO ;
 
 'FIN' ;