* 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 ;
tcpu = 'TEMPS' 'NOEC' ;
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 - All rights reserved.
Disclaimer