Télécharger nlin_cavity.dgibi

Retour à la liste

* fichier :  nlin_cavity.dgibi
*************************************************************************
*
* NOM         : nlin_cavity.dgibi
*
* DESCRIPTION : We compute the flow governed by the incompressible
*               Navier-Stokes equations in a lid driven cavity.
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               mél : beccantini@cea.fr
*************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
*************************************************************************
*************************************************************************
*************************************************************************
******* PERSONAL PROCEDURES *********************************************
*************************************************************************
*************************************************************************
*
*BEGINPROCEDUR matmas
*
*************************************************************************
*
* NOM         : MATMAS
*
* DESCRIPTION : We compute 
*
*               (Np , Nd )
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               mél : beccantini@cea.fr
*************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
*************************************************************************
*
*
'DEBPROC' MATMAS ;
*
'ARGUMENT'  _mt*'MAILLAGE'    ;
'ARGUMENT'  nomp*'MOT     '  ;
'ARGUMENT'  nomd*'MOT     '  ;
'ARGUMENT'  discp*'MOT     '  ;
'ARGUMENT'  discd*'MOT     ' ;
*
*  Arguments
*  _mt   = MAILLAGE, quaf mesh
*  nomp  = MOT, name of the primal variable
*  nomd  = MOT, name of the dual variable
*  discp = MOT, name of the discretization of the primal variable
*  discd = MOT, name of the discretization of the dual variable
*
idim   = 'VALEUR' 'DIME' ;
*
*   Matrix A for the operator NLIN
*
numop  = 1 ;
numder = idim ;
numvar = 1 ;
* The variable involved in A is the primal variable T
numdat = 0 ;
numcof = 0 ;
*
discg = 'LINE' ;
* Geometrical discretisation
methgau = 'GAU7' ;
* Gauss points involved
*
A = ININLIN numop numvar numdat numcof numder ;   
A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomp ;
A . 'VAR' . 1 . 'DISC'   = discp ;
*
*   Matrix B for the operator NLIN
*
B = ININLIN numop numvar 0 0 numder ;
* For B the coef no coefficients are involved
* but 1
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
B . 'VAR' . 1 . 'DISC'   = discd ;
*
A . 1 . 1 . 0 = 'LECT' ;
B . 1 . 1 . 0 = 'LECT'  ;
*
mat = 'NLIN' discg _mt A B methgau ;
'RESPRO' mat ;
'FINPROC' ;
*
*ENDPROCEDUR matmas
*BEGINPROCEDUR matkv
*
*************************************************************************
*
* NOM         : MATKV
*
* DESCRIPTION : We compute the matrix
*
*               ((cx \dep{u_x}{x} '+' cy \dep{u_x}{y}), Nvx)
*               ((cx \dep{u_y}{x} '+' cy \dep{u_y}{y}), Nvy)
*
*               where Nv = shape function for (u_x,u_y)^T
*
*               Names of primal variables = names of the dual variables
*                = 'UX' 'UY'
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               mél : beccantini@cea.fr
*************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
*************************************************************************
*
'DEBPROC' MATKV ;
'ARGUMENT'   _mt*'MAILLAGE' ;
'ARGUMENT' discv*'MOT     ' ;
'ARGUMENT' discc*'MOT     ' ;
'ARGUMENT'     cx*'CHPOINT' ;
'ARGUMENT'     cy*'CHPOINT' ;
*
*  _mt   = solid quaf mesh
*  discv = discretization type ux, uy
*  discc = discretization type for cx, cy
*  cx    = CHPOINT ('SCAL')
*  cy    = CHPOINT ('SCAL')
*
idim   = 'VALEUR' 'DIME' ;
*
numop  = 2 ;
numder = idim ;
numvar = 2 ;
* Two variables, ux, uy
numdat = 2 ;
* cx, cy
numcof = 2 ;
* Coef: cx  cy (two identity functions)
discg = 'LINE' ;
* Geometrical discretisation
methgau = 'GAU7' ;
* Gauss points involved
*
*   Matrix A for the operator NLIN
*
A = ININLIN numop numvar numdat numcof numder ;   
A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
A . 'VAR' . 1 . 'DISC'   = discv  ;
A . 'VAR' . 2 .  'NOMDDL' = 'MOTS' 'UY' ;
A . 'VAR' . 2 .  'DISC'   = discv  ;
* cx
A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 1 . 'DISC' =  discc ;
A . 'DAT' . 1 . 'VALEUR' = cx ;
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
* cy 
A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 2 . 'DISC' =  discc ;
A . 'DAT' . 2 . 'VALEUR' = cy ;
A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
* 
*   Matrix B for the operator NLIN
*
numvar = 2 ;
B = ININLIN numop numvar 0 0 numder ;
*
* For B, no coefficients are involved
* but 1
*
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
B . 'VAR' . 1 . 'DISC'   = discv ;
B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
B . 'VAR' . 2 . 'DISC'   = discv ;
*
* Contribution
*
* cx \dep{ux}{x} Nvx
A . 1 . 1 . 1 = 'LECT' 1 ;
B . 1 . 1 . 0 = 'LECT'   ;
* cy \dep{ux}{y}  Nvx
A . 1 . 1 . 2 = 'LECT' 2 ;
* cx \dep{uy}{x} Nvy
A . 2 . 2 . 1 = 'LECT' 1 ;
B . 2 . 2 . 0 = 'LECT'   ;
* cy \dep{uy}{y}  Nd
A . 2 . 2 . 2 = 'LECT' 2 ;
*
mat = 'NLIN' discg _mt A B  methgau ;
'RESPRO' mat ;
'FINPROC' ;
*
*ENDPROCEDUR matkv
*BEGINPROCEDUR matlap
*
*************************************************************************
*
* NOM         : MATLAP
*
* DESCRIPTION : We compute the matrix arising from the computation of
*               the scalar product
*
*               ( -(coef \div( \alpha \grad T)) , Nd )
*
*               where Nd = test function for the dual variable.
*
*               -\int_{\Omega} coef \div( \alpha \grad T) Nd dV =
*                  -\int_{\delta \Omega}
*                     coef \alpha (\grad T \cdot n) Nd dS '+'
*                  +\int_{\Omega}
*                     coef \alpha (\grad T \cdot \grad Nd) dV '+'
*                  +\int_{\Omega}
*                     Nd \alpha (\grad T \cdot \grad coef) dV
*
*               The matrix is issue from the 2nd and 3rd contribution
*               (volume integrals).
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               mél : beccantini@cea.fr
*************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
*************************************************************************
*
*
'DEBPROC' MATLAP ;
*
'ARGUMENT'  _mt*'MAILLAGE'    ;
'ARGUMENT'   nomt*'MOT     '  ;
'ARGUMENT'  disct*'MOT     '  ;
'ARGUMENT'   nomd*'MOT     '  ;
'ARGUMENT'  discd*'MOT     '  ;
'ARGUMENT'  discal*'MOT     ' ;
'ARGUMENT'  alpha*'CHPOINT'   ;
'ARGUMENT'  discof*'MOT     ' ;
'ARGUMENT'  coef*'CHPOINT'    ;
*
*    Arguments
*
*    _mt    = surface QUAF mesh
*    nomt   = name of the primal variable T
*    disct  = discretization type of T
*    nomd   = name of the dual variable
*    discd  = discretization type of the dual variable
*    discal = discretization type of alpha
*    discof = discretization type of coef
*    alpha  = CHPOINT ('SCAL')
*    coef   = CHPOINT ('SCAL')
*
idim   = 'VALEUR' 'DIME' ;
*
*    Four contributions
*
*    (coef alpha \dep{T}{x} \dep{Nd}{x}
*   + coef alpha \dep{T}{y} \dep{Nd}{y}
*   + Nd alpha \dep{T}{x} \dep{coef}{x}
*   + Nd alpha \dep{T}{y} \dep{coef}{y})
*
*   Matrix A for the operator NLIN
*
numop  = 4 ;
numder = idim ;
numvar = 1 ;
* The variable involved in A is the primal variable T
numdat = 2 ;
* Two data: alpha, coef
numcof = 4 ;
* Four functions: f(alpha) = alpha
*                 f(coef)  = coef
*                 f(coef) = \dep{coef}{x}
*                 f(coef) = \dep{coef}{y}
discg = 'LINE' ;
* Geometrical discretisation
methgau = 'GAU7' ;
* Gauss points involved
*
A = ININLIN numop numvar numdat numcof numder ;   
A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomt ;
A . 'VAR' . 1 . 'DISC'   = disct ;
* alpha
A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 1 . 'DISC' = discal ;
A . 'DAT' . 1 . 'VALEUR' = alpha ;
* coef
A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 2 . 'DISC' = discof ;
A . 'DAT' . 2 . 'VALEUR' = coef ;
* Function alpha
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
* Function coef
A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
* Function \dep{coef}{x}
A . 'COF' . 3 . 'COMPOR' = 'D/DX1' ;
A . 'COF' . 3 . 'LDAT' = 'LECT' 2  ;
* Function \dep{coef}{y}
A . 'COF' . 4 . 'COMPOR' = 'D/DX2' ;
A . 'COF' . 4 . 'LDAT' = 'LECT' 2 ;
*
*   Matrix B for the operator NLIN
*
B = ININLIN numop numvar 0 0 numder ;
* For B the coef no coefficients are involved
* but 1
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
B . 'VAR' . 1 . 'DISC'   = discd ;
*
* Contribution
*
* coef alpha \dep{T}{x} \dep{Nd}{x}
A . 1 . 1 . 1 = 'LECT' 1 2 ;
B . 1 . 1 . 1 = 'LECT'  ;
* coef alpha \dep{T}{y} \dep{Nd}{y}
A . 2 . 1 . 2 = 'LECT' 1 2 ;
B . 2 . 1 . 2 = 'LECT'  ;
* Nd alpha \dep{T}{x} \dep{coef}{x}
A . 3 . 1 . 1 = 'LECT' 1 3 ;
B . 3 . 1 . 0 = 'LECT'  ;
* Nd alpha \dep{T}{y} \dep{coef}{y}
A . 4 . 1 . 2 = 'LECT' 1 4 ;
B . 4 . 1 . 0 = 'LECT'  ;
mlapn = 'NLIN' discg _mt A B  methgau ;
'RESPRO' mlapn ;
'FINPROC' ;
*ENDPROCEDUR matlap
*BEGINPROCEDUR matdiv
*
*************************************************************************
*
* NOM         : MATDIV
*
* DESCRIPTION : We compute the matrix arising from the computation of
*               the scalar product
*
*               (div(u),Np)
*
*               where Np = shape function for pressure
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               mél : beccantini@cea.fr
*************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
*************************************************************************
*
*
'DEBPROC' MATDIV ;
*
'ARGUMENT'   _mt*'MAILLAGE' ;
'ARGUMENT'  nompre*'MOT     ' ;
'ARGUMENT'  discp*'MOT     ' ;
'ARGUMENT'  discv*'MOT     ' ;
*
*  _mt   = solid QUAF mesh
*  nompre = name of the pressure
*  discp = discretization type for the pressure
*          (for instance 'LINE')
*  discv = discretization type for the speed
*          (for instance 'QUAF')
*
idim   = 'VALEUR' 'DIME' ;
*
*    Two contributions
*
*     \dep{u_x}{x} Np
*   + \dep{u_y}{y} Np
*
*   Matrix A for the operator NLIN
*
numop  = 1 ;
numder = idim ;
numvar = 2 ;
* The variables involved in A are u_x, u_y
numdat = 0 ;
* Zero data
numcof = 0 ;
* Zero functions (but 1)
discg = 'LINE' ;
* Geometrical discretisation
methgau = 'GAU7' ;
* Gauss points involved
*
A = ININLIN numop numvar numdat numcof numder ;   
A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
A . 'VAR' . 1 . 'DISC'   = discv ;
A . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
A . 'VAR' . 2 . 'DISC'   = discv ;
*
*   Matrix B for the operator NLIN
*
numvar = 1 ;
B = ININLIN numop numvar 0 0 numder ;
* For B the coef no coefficients are involved
* but 1
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nompre ;
B . 'VAR' . 1 . 'DISC'   = discp ;
*
* Contribution
*
* 1 \dep{u_x}{x} Np
A . 1 . 1 . 1 = 'LECT'  ;
B . 1 . 1 . 0 = 'LECT'  ;
* 1 \dep{u_y}{y} Np
A . 1 . 2 . 2 = 'LECT'  ;
*
mdiv = 'NLIN' discg _mt A B  methgau ;
'RESPRO' mdiv ;
'FINPROC' ;
*
*ENDPROCEDUR matdiv
*BEGINPROCEDUR matpre
*
*************************************************************************
*
* NOM         : MATPRE
*
* DESCRIPTION : We compute the integral of volume of
*
*              (coef \dep{p}{x} , Nv) 
*              (coef \dep{p}{y} , Nv) 
*
*              where Nv = test function for u_x = test function for u_y.
*
*              \int_{\Omega} coef  \dep{p}{x} Nv dV =
*              \int_{\delta \Omega} coef p Nv (i \cdot n) dS '+'
*              -\int_{\Omega} coef p \dep{nv}{x} dV '+'
*              -\int_{\Omega} Nv p \dep{coef}{x} dV
*
*              Here we only compute the 2nd and 3rd contribution (volume
*              integrals). 
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
*               mél : beccantini@cea.fr
*
*************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
*************************************************************************
*
*
'DEBPROC' MATPRE ;
*
'ARGUMENT'   _mt*'MAILLAGE'   ;
'ARGUMENT'  nompre*'MOT     ' ;
'ARGUMENT'  discp*'MOT     '  ;
'ARGUMENT'  discv*'MOT     '  ;
'ARGUMENT'  discc*'MOT     '  ;
'ARGUMENT'  coef*'CHPOINT '   ;
*
* _mt    = surface QUAF mesh
* nompre = name of the pressure (usuallt 'LX')
* discp  = discretization type of p (usually LINE)
* discv  = discretization type of v (usually QUAF)
* discc = discretization type of coef
*  coef   = CHPOINT ('SCAL')
*
idim   = 'VALEUR' 'DIME' ;
*
*  Matrix A for the operator NLIN
*
numop  = 4 ;
numder = idim ;
numvar = 1 ;
* The variable involved in A is the primal p
numdat = 1 ;
* The data: coef
numcof = 3 ;
* Three functions: f(coef)= coef
*                  f(coef)=\dep{coef}{x}
*                  f(coef)=\dep{coef}{y}
*
discg = 'LINE' ;
* Geometrical discretisation
methgau = 'GAU7' ;
* Gauss points involved
*
A = ININLIN numop numvar numdat numcof numder ;   
A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nompre ;
A . 'VAR' . 1 . 'DISC'   = discp ;
*
* coef
*
mcoef = 'EXTRAIRE' coef 'COMP' ;
A . 'DAT' . 1 . 'NOMDDL' = mcoef ;
A . 'DAT' . 1 . 'DISC' =  discc ;
A . 'DAT' . 1 . 'VALEUR' = coef ;
* coef, \dep{coef}{x}, \dep{coef}{y}
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
A . 'COF' . 2 . 'COMPOR' = 'D/DX1' ;
A . 'COF' . 2 . 'LDAT' = 'LECT' 1 ;
A . 'COF' . 3 . 'COMPOR' = 'D/DX2' ;
A . 'COF' . 3 . 'LDAT' = 'LECT' 1 ;
*
*   Matrix B for the operator NLIN
*
numvar = 2 ;
*
* ux, uy
*
numdat = 1 ;
numcof = 1 ;
B = ININLIN numop numvar numdat numcof numder ;
*
* For B the coef no coefficients are involved
* but -1
*
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
B . 'VAR' . 1 . 'DISC'   = discv ;
B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
B . 'VAR' . 2 . 'DISC'   = discv ;
* -1
B . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
B . 'DAT' . 1 . 'DISC' =  'CSTE' ;
B . 'DAT' . 1 . 'VALEUR' = -1.0 ;
B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
B . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
*
*   -\int_{\Omega} coef p \dep{Nv}{x} dV
*
A . 1 . 1 . 0 = 'LECT' 1 ;
B . 1 . 1 . 1 = 'LECT' 1 ;
*
*   -\int_{\Omega} Nv p \dep{coef}{x} dV
*
A . 2 . 1 . 0 = 'LECT' 2 ;
B . 2 . 1 . 0 = 'LECT' 1 ;
*
*   -\int_{\Omega} coef p \dep{Nv}{y} dV (second dual variable)
*
A . 3 . 1 . 0 = 'LECT' 1 ;
B . 3 . 2 . 2 = 'LECT' 1 ;
*
*   -\int_{\Omega} Nv p \dep{coef}{y} dV (second dual variable)
*
A . 4 . 1 . 0 = 'LECT' 3 ;
B . 4 . 2 . 0 = 'LECT' 1 ;
*
mat = 'NLIN' discg _mt A B methgau ;
'RESPRO' mat ;
'FINPROC' ;
*
*ENDPROCEDUR matpre
*BEGINPROCEDUR resoup
*************************************************************************
**** Resolution of a linear system **************************************
*************************************************************************
*
'DEBPROC' RESOUP ;
   'ARGUMENT' LOGTPS*'LOGIQUE ' ;
   'ARGUMENT' mat*'MATRIK  ' ;
   'ARGUMENT' matpre/'MATRIK  ' ;
   'ARGUMENT' smb*'CHPOINT ' ;
   'ARGUMENT' ccl*'CHPOINT ' ;
   'ARGUMENT' res*'FLOTTANT' ;
   'ARGUMENT' nit*'ENTIER'   ;
   'ARGUMENT' pre/'ENTIER'   ;
   'ARGUMENT' gggtre/'TABLE   ' ;
   'ARGUMENT' gggtcv/'LOGIQUE ' ;
   'ARGUMENT' typslv/'MOT     ' ;
*
*  Det. tt s.t.
*
*  mat . tt = smb
*  tt(boundary) = ccl
*
*  
*  mat = matrix to "inverse"
*  matpre = matrix for which preconditioner has been computed (optional)
*  smb = right hand side
*  ccl = imposed boundary conditions
*  res = rvk . 'RESID'
*  nit = max number of linear iteration (for an iterative solver)
*  pre = type of preconditioning matrix (for an iterative solver)
*  gggtre = table of resolution to store 'XINIT'
*  gggtcv = graphic of the history for the iterative solver
*           (optional FAUX)
*  typslv = 'DIRECT'/'ITER' (optional, ITER)
*
   'SI' ('NON' ('EXISTE' typslv)) ;
      tslv = VRAI ;
   'SINON' ;      
      'SI' ('EGA' typslv 'DIRECT') ;
         tslv = FAUX ;
      'SINON' ;
         'SI' ('EGA' typslv 'ITER') ;
            tslv = VRAI ;
         'SINON' ;
            'ERREUR' 'Solveur DIRECT ou ITERatif ?' ;
         'FINSI' ;
      'FINSI' ;
   'FINSI' ;
    gggniv = 1 ;
   'SI' ('NON' ('EXISTE' pre)) ;
      pre = 5 ;
   'FINSI' ;
   'SI' ('OU' ('NON' ('EXISTE' gggtcv)) ('NON' tslv)) ;
      gggtcv = FAUX ;
   'FINSI' ;      
   rv = 'EQEX' ;
   rvk = rv . 'METHINV' ;
   'SI' ('EXISTE' matpre) ;
      rvk . 'MATASS' = matpre ;
      rvk . 'MAPREC' = matpre ;
   'SINON' ;      
      rvk . 'MATASS' = mat ;
      rvk . 'MAPREC' = mat ;
   'FINSI' ;      
   'SI' tslv ;
      rvk . 'TYPINV' = 3 ;
   'SINON' ;      
      rvk . 'TYPINV' = 1 ;
   'FINSI' ;      
   rvk . 'SCALING' = 1 ;
*  Scaling factor   
   rvk . 'OUBMAT' = 1 ;
*  oubmat = 1 -> We eliminate the elementary matrix
*  'ILUTPPIV' -> 1 We always search for a new pivot
*             -> 0 We do not search for a new pivot
   rvk . 'ILUTPPIV' = 0.1D0 ;
*   rvk . 'ILUTPPIV' = 0.01D0 ;
*   rvk . 'ILUTPPIV' = 0.00001D0 ;
*   rvk . 'ILUTPPIV' = 0.D0 ;
   rvk . 'IMPINV' = gggniv ;
   rvk . 'TYRENU' = 'SLOA' ;
   rvk . 'PCMLAG' = 'APR2' ;
   rvk . 'NITMAX' = nit ;
   rvk . 'RESID' = res ;
   rvk . 'PRECOND' = pre ;
   rvk . 'BCGSBTOL' = 1.D-200 ;
   rvk . 'ILUTLFIL' = 1.5 ;
*
   'SI' ('EXISTE' gggtre) ;
      'SI' ('EXISTE' gggtre 'XINIT') ;
         rvk . 'XINIT' = gggtre . 'XINIT' ;
      'FINSI' ;
   'FINSI' ;      
   'TEMPS' 'ZERO' ;
   tt = 'KRES' mat 'TYPI' rvk 'SMBR' smb 'CLIM' ccl ;
   TABTPS = TEMP 'NOEC';
   tcpu = TABTPS.'TEMPS_CPU'.'INITIAL';
   tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
   'SI' LOGTPS ;
      'MESSAGE' ('CHAINE' 'tcpu (s) = ' tcpus) ;
   'FINSI' ;
   'SI' ('ET' gggtcv ('EXISTE' rvk 'CONVINV')) ;
      lcvg = rvk . 'CONVINV' ;
      nit = 'DIME' lcvg ;
      'SI' ('>' nit 1) ;
         lit = 'PROG' 1.D0 PAS 1.D0 NPAS ('-' nit 1) ;
         evit = 'EVOL' 'MANU' 'iter' lit 'log residu'
          ('/' ('LOG' lcvg) ('LOG' 10.D0)) ;
         'DESSIN' evit 'TITR' titglob ;
      'FINSI' ;         
   'FINSI' ;
   'SI' ('EXISTE' gggtre) ;
      gggtre . 'XINIT' = tt ;
   'FINSI' ;      
   'RESPRO' tt ;
*
* End of procedure file RESOU
*
'FINPROC' ;
*
*************************************************************************
**** End of the Resolution of a linear system ***************************
*************************************************************************
*ENDPROCEDUR resoup
*************************************************************************
*************************************************************************
******* END OF PERSONAL PROCEDURES **************************************
*************************************************************************
*************************************************************************
*
 'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'TRAC' 'X' ;
 GRAPH = FAUX ;
* 
* Linear/bilinear continuus pressure 
* QUAFdratic speed   
*
**** Mesh
*
*
* A4         A3
*
*
*
*
* A1         A2
*
 L  = 1.0 ;
 A1 = 0.0 0.0 ;
 A2 = L   0.0 ;
 A3 = L   L ;
 A4 = 0.0 L ;
*
 DX = L '/' (2 '*' 10) ;
 DY = DX ;
*
 A1A2 = A1 'DROIT' A2 'DINI' DX 'DFIN' DX ;
 A2A3 = A2 'DROIT' A3 'DINI' DY 'DFIN' DY ;
 A3A4 = A3 'DROIT' A4 'DINI' DX 'DFIN' DX ;
 A4A1 = A4 'DROIT' A1 'DINI' DY 'DFIN' DY ;
*
 DOM1 = 'SURFACE' (A1A2 'ET' A2A3 'ET' A3A4 'ET' A4A1) 'PLAN' ;
*DOM1 = 'CHANGER' DOM1 'QUAF' ;
 'SI' GRAPH ;
    'TRACER' DOM1 'TITRE' 'Mesh' ;
 'FINSI' ;
 QDOM1 = 'CHANGER' DOM1 'QUAF' ;
 QA4A1 = 'CHANGER' A4A1 'QUAF' ;
 QA1A2 = 'CHANGER' A1A2 'QUAF' ;
 QA2A3 = 'CHANGER' A2A3 'QUAF' ;
 QA3A4 = 'CHANGER' A3A4 'QUAF' ;
 'ELIMINATION' (QA4A1 'ET' QDOM1) 1.0D-4 ;
 'ELIMINATION' (QA1A2 'ET' QDOM1) 1.0D-4 ;
 'ELIMINATION' (QA2A3 'ET' QDOM1) 1.0D-4 ;
 'ELIMINATION' (QA3A4 'ET' QDOM1) 1.0D-4 ;
*
* Boundary conditions
*
 LIMDIR   = 'CHANGER' ('CONTOUR' QDOM1) 'POI1' ;
 LIMDIR1  = 'DIFF' ('CHANGER' 'POI1' QA3A4) ('MANUEL' 'POI1' A3) ;
 LIMDIR1  = 'DIFF' LIMDIR1 ('MANUEL' 'POI1' A4) ;
 LIMDIR2  = 'DIFF' LIMDIR LIMDIR1 ;
 UN1 = 'MANUEL' 'CHPO' LIMDIR1 2 'UX' 1.0 'UY' 0.0 ;
 UN2 = 'MANUEL' 'CHPO' LIMDIR2 2 'UX' 0.0 'UY' 0.0 ;
 UNLIM = UN1 '+' UN2 ;
 Re = 1000. ;
 PNLIM = 'MANUEL' 'CHPO' ('MANUEL' 'POI1' A4) 1 'LX' 0.0 ;
 UNPLIM = UNLIM '+' PNLIM ;
*
* IC
*
 UN = ('MANUEL' 'CHPO' QDOM1 2 'UX' 0.0 'UY' 0.0) '+' (UN1 '+' UN2) ;
*
 VECN = 'VECTEUR' 0.1 UN ;
 'SI' GRAPH ;
    'TRACER' DOM1 VECN ;
 'FINSI' ;
*
* Time step
*
 DT = 1.0D2 '*' DX '/' 1.0 ;
*
 NU = 'MANUEL' 'CHPO' DOM1 1 'SCAL' (1.0 '/' Re) ;
*
* Matrices
*
 MMAS1 = MATMAS QDOM1 'UX' 'UX' 'QUAF' 'QUAF' ;
 MMAS2 = MATMAS QDOM1 'UY' 'UY' 'QUAF' 'QUAF' ;
* MMASP = (MATMAS QDOM1 'LX' 'LX' 'LINE' 'LINE') '*' 0.0 ;
* Mass matrix for the pressure to help the linear solver to converge.
* From "Gounand's RECIPIES"
* MMAS = MMAS1 'ET' MMAS2 'ET' MMASP;
 MMAS = MMAS1 'ET' MMAS2 ;
 MKON = MATKV QDOM1 'QUAF' 'QUAF' ('EXCO' 'UX' UN) ('EXCO' 'UY' UN) ;
 CHONE = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
 MLAP1 = MATLAP  QDOM1 'UX' 'QUAF' 'UX' 'QUAF' 'LINE' NU 'LINE' CHONE ;
 MLAP2 = MATLAP  QDOM1 'UY' 'QUAF' 'UY' 'QUAF' 'LINE' NU 'LINE' CHONE ;
 MLAP = MLAP1 'ET' MLAP2 ;
 MDIV = MATDIV QDOM1 'LX' 'LINE' 'QUAF' ;
 MPRE = MATPRE QDOM1 'LX' 'LINE' 'QUAF' 'LINE' CHONE ;
*
* Total matrix
*
 MMASSDT = MMAS '/' DT ;
 mconst = MMASSDT 'ET' MLAP ;
 mconst = mconst 'ET' MPRE ;
 mconst = mconst 'ET' MDIV ;
 mtot = mconst 'ET' MKON ;
*
* Resolution
*
 TPS = 0.0 ;
 UN = 'COPIER' UN ;
 UNM = 'COPIER' UN ;
 LIT = 'PROG' ;
 LER = 'PROG' ;
 ERRO = 1.0D10 ;
* 
 TABRESUN = 'TABLE' ;
 UNP =  UN '+' ('MANUEL' 'CHPO' DOM1 1 'LX' 0.0) ;
 NITERU = 100 ;
 PREC = 3 ;
 METINV = MOT 'DIRECT' ;
* 
 'REPETER' BL1 200  ;
    TPS = TPS '+' DT ;
    qv  = MMASSDT '*' UN ;
    mtotik = 'KOPS' 'RIMA' mtot ;
    TABRESUN . 'XINIT' = UNP ;
    UNP  = RESOUP FAUX mtotik mtotik qv UNPLIM 1.0D-12 NITERU
               PREC TABRESUN METINV ;
    UN = 'EXCO' ('MOTS' 'UX' 'UY') UNP ('MOTS' 'UX' 'UY') ;
    'SI' (ERRO < 1.0D-5) ;
       'QUITTER' BL1 ;
    'FINSI' ;
    'SI' (((&BL1 '/' 5) '*' 5) 'EGA' &BL1) ;
        ERRO = 'MAXIMUM' (UN '-' UNM) 'ABS' ;
        'MESSAGE' ('CHAINE' 'ITER = ' &BL1 ',  TPS = ' TPS  ',  ERRO = '
           ERRO) ;
        UNM = 'COPIER' UN ;
        LIT = LIT 'ET' ('PROG' &BL1) ;
        L10 = 'LOG' 10. ;
        LER = LER 'ET' ('PROG' (('LOG' (ERRO '+' 1.0D-20)) '/' L10))  ;
        EVER = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LER ;
        'SI' GRAPH ;
           'DESSIN' EVER 'TITR' 'Convergence history' 'NCLK' ;
        'FINSI' ;
    'FINSI' ;
*
****** We update the convective matrix
*
    MKON = MATKV QDOM1 'QUAF' 'QUAF' ('EXCO' 'UX' UN) ('EXCO' 'UY' UN) ;
    MTOT = MCONST 'ET' MKON ;
 'FIN' BL1 ;
*
 UX = 'EXCO' UN 'UX' ;
 UY = 'EXCO' UN 'UY' ;
*
***** Post treatment
*
 VECN = 'VECTEUR' UN ;
 'SI' GRAPH ;
    'TRACER' UX DOM1 VECN 'TITR' 'ux' ;
    'TRACER' UY DOM1 VECN 'TITR' 'uy' ;
 'FINSI' ;
*
 QXMED = QA4A1 'PLUS' ((L '/' 2.) 0.0) ;
 QYMED = QA1A2 'PLUS' (0.0 (L '/' 2.)) ;
 'ELIMINATION' QDOM1 (DX '/' 100.) QXMED ;
 'ELIMINATION' QDOM1 (DX '/' 100.) QYMED ;
*
**** Evolution objects
*
 LUX = 'EXTRAIRE' ('EVOL' 'CHPO' ('EXCO' 'UX' UN) QXMED)  'ORDO' ;
 LY  = 'EXTRAIRE' ('EVOL' 'CHPO' ('COORDONNEE' 2 QXMED) QXMED) 'ORDO' ;
 UX_Y = 'EVOL' 'MANU' 'ux' LUX 'y' LY ;
* 
* 
 LUY = 'EXTRAIRE' ('EVOL' 'CHPO' ('EXCO' 'UY' UN) QYMED)  'ORDO' ;
 LX  = 'EXTRAIRE' ('EVOL' 'CHPO' ('COORDONNEE' 1 QYMED) QYMED) 'ORDO' ;
 UY_X = 'EVOL' 'MANU' 'x' LX 'uy' LUY ;
* 
*
*
 TAB1 = 'TABLE' ;
 TAB1 . 'TITRE'= 'TABLE' ;
 TAB1 . 2 = 'MARQ CROI NOLI';
 TAB1 . 'TITRE' . 2  = 'Reference (Su)' ;
 TAB1 . 1 = 'REGU' ;
 TAB1 . 'TITRE' . 1  = 'Numerical res.' ;
*
* Solution of Su
*
 LXSU = 'PROG' 0   0.0625  0.0703   0.0781  0.0938  0.1563  0.2266  
     0.2344   0.5  0.8047  0.8594   0.9063  0.9453  0.9531  0.9609  
     0.9688   1 ;
 LUY_X_SU = 'PROG' 0   0.2763  0.2918   0.3053  0.3282  0.3711  0.3279 
      0.3193   0.0243 -0.317 -0.4245 -0.5182 -0.3972 -0.3421 -0.2816 
      -0.2175   0 ;
 EVUYX = 'EVOL' 'MANU' 'x' LXSU 'uy' LUY_X_SU ;
*
 LYSU = 'PROG' 0. 0.0547 0.0625 0.0703 0.1016 0.1719 0.2812 0.4531
     0.5  0.6172  0.7344  0.8516  0.9531  0.9609  0.9688  0.9766  1 ;
 LUX_Y_SU = 'PROG' 0 -0.1788 -0.1999 -0.2204 -0.2972 -0.3834  -0.2759   
    -0.1058 -0.0605 0.0564 0.1857 0.3316 0.466 0.5109 0.5743 0.6582 1 ;
 EVUXY = 'EVOL' 'MANU' 'ux' LUX_Y_SU 'y' LYSU ;
*
 'SI' GRAPH ;
    'DESSIN' (UX_Y 'ET' EVUXY) 'TITRE'
       ('CHAINE' 'ux at y = ' (L '/' 2.)) 'LEGE' 
       TAB1 ;
    'DESSIN' (UY_X 'ET' EVUYX) 'TITRE'
       ('CHAINE' 'uy at x = ' (L '/' 2.)) 'LEGE' 
       TAB1 ;
 'FINSI' ;
*
* We check whether the convergence is reached
*
 NITER = &BL1 ;
*
 'SI' (NITER > 100) ;
     'MESSAGE' 'Convergence not reached' ;
     'ERREUR' 5 ;
 'FINSI' ;
*
* Difference between the solution of Su and the one here obtained
*
 LUX_Y_HE = 'IPOL' LYSU LY LUX ;
 ERRO = (LUX_Y_HE '-' LUX_Y_SU) 'ABS' ;
 EVERX =  'EVOL'  'MANU' 'y' lysu 'UX'  ERRO ;
 LUY_X_HE = 'IPOL' LXSU LX LUY ;
 ERRO = (LUY_X_HE '-' LUY_X_SU) 'ABS' ;
 EVERY =  'EVOL'  'MANU' 'x' lysu 'UY'  ERRO ;
 'SI' GRAPH ;
    'DESSIN' EVERX 'TITRE' 'Error on ux'  ;
    'DESSIN' EVERY 'TITRE' 'Error on uy'  ;
 'FINSI' ;
 aax = 'MAXIMUM' ('EXTRAIRE' ('PRIM'  EVERX) 'ORDO') ;
 aay = 'MAXIMUM' ('EXTRAIRE' ('PRIM'  EVERY) 'ORDO') ;
*
 'SI' (aax > 1.0D-2) ;
     'ERREUR' 5 ;
 'FINSI' ;
 'SI' (aay > 1.0D-2) ;
     'ERREUR' 5 ;
 'FINSI' ;
*
 'FIN' ;
*
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales