* fichier : nlin_te_unstat.dgibi
*************************************************************************
*
* NOM : nlin_te_unstat.dgibi
*
* DESCRIPTION : We compute the flow governed by the Navier-Stokes
* equations, in a T jonction with high temperature
* difference.
* Both Boussinesq approximation and Low Mach number
* approximation are considered.
* See also technical report SFME/LTMF/RT/O5-067/A
* In this test case we simply perform 2 time iterations
* and we verify that two point-fix algorithms converge
* (to the same solution).
*
* 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 !
*************************************************************************
*
*
*************************************************
****** LIST OF SYMBOLS **************************
*************************************************
*
* \int = integral
* \div = divergence
* \der = derivative (\der{a}{b} = \frac{d a}{d b}
* \grad = gradient
* \cdot = scalar product
* ^ = power
*
* (a,b) = \int_{\Omega} a * b dV
*
* Np = test function for the divergence equations
* = shape function for the pressure
* Nt = test function for the temperature
* = shape function for the temperature equation
* Nv = test function for the speed equations (x,y)
* = shape function for the speed
*
*************************************************************************
*************************************************************************
******* 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
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 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 matflu
*
************************************************************************
*
* NOM : MATFLU
*
* DESCRIPTION : We compute the matrix arising from
*
* \int_{\delta \Omega} (coef q \cdot n) Nd dS
*
* where q is the primal variable ('UX', 'UY')
* Nd is the shape function of the dual variable
*
* 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 !
************************************************************************
*
*
**** We compute the matrix
*
* \int_{\delta \Omega} (q \cdot n) Nd dS
*
* N.B. Be carefull in specifying n (the normal at the surfaces)
* No singular points are aviable.
*
* Arguments
*
* _st = MAILLAGE, quaf mesh of the surface
* discp = MOT, name of the discretisation of the primal variable
* nomd = MOT, name of the dual variable
* discd = MOT, name of the discretisation of the dual variable
* disco = MOT, discretisation of the coefficient
* coef = coefficient
* discn = MOT, name of the discretisation of the normals
* nx = CHPOINT
* ny = CHPOINT
*
* nomp = name of the primal variable are 'UX' 'UY'
*
'DEBPROC' MATFLU ;
'ARGUMENT' _st*'MAILLAGE' ;
'ARGUMENT' discp*'MOT ' ;
'ARGUMENT' nomd*'MOT ' ;
'ARGUMENT' discd*'MOT ' ;
'ARGUMENT' discco*'MOT ' ;
'ARGUMENT' coef*'CHPOINT' ;
'ARGUMENT' discn*'MOT ' ;
'ARGUMENT' nx*'CHPOINT' ;
'ARGUMENT' ny*'CHPOINT' ;
*
idim = 'VALEUR' 'DIME' ;
*
* Matrix A for the operator NLIN
*
numop = 1 ;
numder = idim ;
numvar = 2 ;
* The variables involved in A are qx, qy
numdat = 3 ;
* coef, nx, ny
numcof = 3 ;
* coef, nx, ny ('IDEN')
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' = discp ;
A . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
A . 'VAR' . 2 . 'DISC' = discp ;
*
* coef
A . 'DAT' . 1 . 'NOMDDL' = 'EXTRAIRE' coef 'COMP';
A . 'DAT' . 1 . 'DISC' = discco ;
A . 'DAT' . 1 . 'VALEUR' = coef ;
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
* nx
A . 'DAT' . 2 . 'NOMDDL' = 'EXTRAIRE' nx 'COMP';
A . 'DAT' . 2 . 'DISC' = discn ;
A . 'DAT' . 2 . 'VALEUR' = nx ;
A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
* ny
A . 'DAT' . 3 . 'NOMDDL' = 'EXTRAIRE' ny 'COMP';
A . 'DAT' . 3 . 'DISC' = discn ;
A . 'DAT' . 3 . 'VALEUR' = ny ;
A . 'COF' . 3 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 3 . 'LDAT' = 'LECT' 3 ;
*
* Matrix B for the operator NLIN
*
numvar =1 ;
B = ININLIN numop numvar 0 0 numder ;
*
* For B no coefficients are involved
* but 1
*
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomd ;
B . 'VAR' . 1 . 'DISC' = discd ;
* (coef * nx * qx '+' coef * ny qy) Nd
A . 1 . 1 . 0 = 'LECT' 1 2;
A . 1 . 2 . 0 = 'LECT' 1 3;
B . 1 . 1 . 0 = 'LECT' ;
*
mat = 'NLIN' discg _st A B methgau ;
'RESPRO' mat ;
'FINPROC' ;
*
*ENDPROCEDUR matflu
*BEGINPROCEDUR chpflu
*
************************************************************************
*
* NOM : CHPFLU
*
* DESCRIPTION : We compute the CHPOIN
*
* -q = lambda (\grad(T)
*
* on the dof of the dual variable
* T = CHPOINT
*
* N.B. Be carefull in specifying n (the normal at the
* surfaces).
* No singular points are admitted.
*
* 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' CHPFLU ;
*
'ARGUMENT' _mt*'MAILLAGE' ;
'ARGUMENT' disct*'MOT ' ;
'ARGUMENT' T*'CHPOINT' ;
'ARGUMENT' discd*'MOT ' ;
'ARGUMENT' discla*'MOT ' ;
'ARGUMENT' lambda*'CHPOINT' ;
*
* _mt = MAILLAGE, solid quaf mesh
* disct = MOT, name of the discretisation of T
* T = CHPOINT (one component, the name of which is free)
* discd = MOT, name of the discretisation of the dual variable
* discla = MOT, name of the discretisation of lambda
* lambda = CHPOINT (one component, SCAL)
*
idim = 'VALEUR' 'DIME' ;
*
*
* Matrix A for the operator NLIN
*
numop = 2 ;
numder = idim ;
numvar = 1 ;
* One variable, T
numdat = 1 ;
* -lambda
numcof = 1 ;
* Coef: -lambda (three 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' = 'EXTR' TN 'COMP' ;
A . 'VAR' . 1 . 'DISC' = disct ;
* lambda
A . 'DAT' . 1 . 'NOMDDL' = ('EXTRAIRE' lambda 'COMP') ;
A . 'DAT' . 1 . 'DISC' = discla ;
A . 'DAT' . 1 . 'VALEUR' = lambda ;
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
*
* 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' = discd ;
B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
B . 'VAR' . 2 . 'DISC' = discd ;
*
* Contribution
*
* lambda \dep{T}{x} Nd
A . 1 . 1 . 1 = 'LECT' 1 ;
B . 1 . 1 . 0 = 'LECT' ;
* lambda \dep{T}{y} Nd
A . 2 . 1 . 2 = 'LECT' 1 ;
B . 2 . 2 . 0 = 'LECT' ;
*
mat = 'NLIN' discg _mt A B methgau ;
*
* We compute coef
*
matm1 = MATMAS _mt 'UX' 'UX' discd discd ;
matm2 = MATMAS _mt 'UY' 'UY' discd discd ;
matm = matm1 'ET' matm2 ;
mtotik = 'KOPS' 'RIMA' matm ;
smbt = mat '*' t ;
chlim pipi = 'KOPS' matrik ;
'MESSAGE' ;
coef = RESOUP FAUX mtotik mtotik smbt chlim 1.0D-12 200 'ITER' ;
*
'RESPRO' coef ;
'FINPROC' ;
*
*ENDPROCEDUR chpflu
*BEGINPROCEDUR matkt
*
************************************************************************
*
* NOM : MATKT
*
* DESCRIPTION : We compute the matrix
*
* ((u_x \dep{T}{x} '+' u_y \dep{T}{y}) , Nt)
*
* 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' MATKT ;
*
'ARGUMENT' _mt*'MAILLAGE' ;
'ARGUMENT' discv*'MOT ' ;
'ARGUMENT' ux*'CHPOINT' ;
'ARGUMENT' uy*'CHPOINT' ;
'ARGUMENT' nomt *'MOT ' ;
'ARGUMENT' disct*'MOT ' ;
*
* discv = discretisation type for u_x and u_y
* ux = chpoint ('SCAL', on discv)
* uy = chpoint ('SCAL', on discv)
* nomt = name of the primal and the dual variable t
* disct = discretisation type for t
*
idim = 'VALEUR' 'DIME' ;
*
numop = 1 ;
numder = idim ;
numvar = 1 ;
* One variable, T
numdat = 2 ;
* ux, uy
numcof = 2 ;
* Coef: ux uy (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' nomt ;
A . 'VAR' . 1 . 'DISC' = disct ;
* ux
A . 'DAT' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 1 . 'DISC' = discv ;
A . 'DAT' . 1 . 'VALEUR' = ux ;
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
* uy
A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 2 . 'DISC' = discv ;
A . 'DAT' . 2 . 'VALEUR' = uy ;
A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
*
* Matrix B for the operator NLIN
*
numvar = 1 ;
B = ININLIN numop numvar 0 0 numder ;
*
* For B, no coefficients are involved
* but 1
*
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' nomt ;
B . 'VAR' . 1 . 'DISC' = disct ;
*
* Contribution
*
* ux \dep{T}{x} Nd
A . 1 . 1 . 1 = 'LECT' 1 ;
B . 1 . 1 . 0 = 'LECT' ;
* uy \dep{T}{y} Nd
A . 1 . 1 . 2 = 'LECT' 2 ;
*
mat = 'NLIN' discg _mt A B methgau ;
'RESPRO' mat ;
'FINPROC' ;
*
*ENDPROCEDUR matkt
*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 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 mattau
*
************************************************************************
*
* NOM : MATTAU
*
* DESCRIPTION : We compute
*
* ( (-coef*div(tau))|x, Nv)
* ( (-coef*div(tau))|y, Nv)
*
* where
*
* Nv is the form function for the speed (u_x,u_y)^T
* = the test function for the speed equations
*
* tau_{xy} = mu (dep{u_x}{y} '+' dep{u_y}{x} '-'
* \delta_{x,y} div (u))
*
* These terms are integrated by part. Only the volumic contribution
* is computed.
*
*
* 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' MATTAU ;
*
'ARGUMENT' _mt*'MAILLAGE' ;
'ARGUMENT' discv*'MOT ' ;
'ARGUMENT' discmu*'MOT ' ;
'ARGUMENT' discco*'MOT ' ;
'ARGUMENT' mu*'CHPOINT' ;
'ARGUMENT' coef*'CHPOINT' ;
*
* _mt = surface QUAF mesh
* discv = discretization type of the speed
* discmu = discretization type of mu
* discco = discretization type of coef
* mu = CHPOINT ('SCAL')
* coef = CHPOINT ('SCAL')
*
* N.B. Names of the primal variable = names of the dual variable =
* 'UX', 'UY'
*
idim = 'VALEUR' 'DIME' ;
*
* Matrix A for the operator NLIN
*
numop = 4 ;
numder = idim ;
numvar = 2 ;
*
numdat = 3 ;
* The data: mu, 2./3., 4./3.
numcof = 3 ;
*
* One function: f(mu)= mu
*
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 ;
*
* mu
*
mmu = 'EXTRAIRE' mu 'COMP' ;
A . 'DAT' . 1 . 'NOMDDL' = mmu ;
A . 'DAT' . 1 . 'DISC' = discmu ;
A . 'DAT' . 1 . 'VALEUR' = mu ;
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
*
* -2/3
*
A . 'DAT' . 2 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 2 . 'DISC' = 'CSTE' ;
A . 'DAT' . 2 . 'VALEUR' = (-2./3.) ;
A . 'COF' . 2 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 2 . 'LDAT' = 'LECT' 2 ;
*
* 4/3
*
A . 'DAT' . 3 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'DAT' . 3 . 'DISC' = 'CSTE' ;
A . 'DAT' . 3 . 'VALEUR' = (4./3.) ;
A . 'COF' . 3 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 3 . 'LDAT' = 'LECT' 3 ;
*
* Matrix B for the operator NLIN
*
numdat = 1 ;
numcof = 3 ;
B = ININLIN numop numvar numdat numcof numder ;
*
B . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'UX' ;
B . 'VAR' . 1 . 'DISC' = discv ;
B . 'VAR' . 2 . 'NOMDDL' = 'MOTS' 'UY' ;
B . 'VAR' . 2 . 'DISC' = discv ;
*
B . 'DAT' . 1 . 'NOMDDL' = ('EXTRAIRE' coef 'COMP') ;
B . 'DAT' . 1 . 'DISC' = discco ;
B . 'DAT' . 1 . 'VALEUR' = coef ;
* coef, \dep{coef}{x}, \dep{coef}{y}
B . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
B . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
B . 'COF' . 2 . 'COMPOR' = 'D/DX1' ;
B . 'COF' . 2 . 'LDAT' = 'LECT' 1 ;
B . 'COF' . 3 . 'COMPOR' = 'D/DX2' ;
B . 'COF' . 3 . 'LDAT' = 'LECT' 1 ;
*
* Dual variable 'UX'
*
* \int
* (2 mu \dep{u_x}{x} coef \dep{N_j}{x} '+'
* 2 mu \dep{u_x}{x} N_j \dep{coef}{x} '+'
* -2/3 mu div(u) coef \dep{N_j}{x} '+'
* -2/3 mu div(u) N_j \dep{coef}{x} '+'
* mu \dep{u_x}{y} coef \dep{N_j}{y} '+'
* mu \dep{u_x}{y} N_j \dep{coef}{y} '+'
* mu \dep{u_y}{x} coef \dep{N_j}{y} '+'
* mu \dep{u_y}{x} N_j \dep{coef}{y}) =
*
* \int
* ((4/3 mu \dep{u_x}{x}) '-' (2/3 mu \dep{u_y}{y})) *
* (coef \dep{N_j}{x} '+' N_j \dep{coef}{x})
* '+'
* (mu \dep{u_x}{y}) '+' (mu \dep{u_y}{x}) *
* (coef \dep{N_j}{y} '+' N_j \dep{coef}{y})
*
A . 1 . 1 . 1 = 'LECT' 1 3 ;
A . 1 . 2 . 2 = 'LECT' 1 2 ;
B . 1 . 1 . 1 = 'LECT' 1 ;
B . 1 . 1 . 0 = 'LECT' 2 ;
*
A . 2 . 1 . 2 = 'LECT' 1 ;
A . 2 . 2 . 1 = 'LECT' 1 ;
B . 2 . 1 . 2 = 'LECT' 1 ;
B . 2 . 1 . 0 = 'LECT' 3 ;
*
* Dual variables 'UY'
*
* \int
* ((4/3 mu \dep{u_y}{y}) '-' (2/3 mu \dep{u_x}{x})) *
* (coef \dep{N_j}{y} '+' N_j \dep{coef}{y})
* '+'
* (mu \dep{u_y}{x}) '+' (mu \dep{u_x}{y}) *
* (coef \dep{N_j}{x} '+' N_j \dep{coef}{x})
*
*
A . 3 . 2 . 2 = 'LECT' 1 3 ;
A . 3 . 1 . 1 = 'LECT' 1 2 ;
B . 3 . 2 . 2 = 'LECT' 1 ;
B . 3 . 2 . 0 = 'LECT' 3 ;
*
A . 4 . 2 . 1 = 'LECT' 1 ;
A . 4 . 1 . 2 = 'LECT' 1 ;
B . 4 . 2 . 1 = 'LECT' 1 ;
B . 4 . 2 . 0 = 'LECT' 2 ;
*
mat = 'NLIN' discg _mt A B methgau ;
'RESPRO' mat ;
'FINPROC' ;
*
*ENDPROCEDUR mattau
*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). The surface integral must be computed elsewhere.
*
* 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
*
'DEBPROC' INTINVT ;
'ARGUMENT' _mt*'MAILLAGE' ;
'ARGUMENT' disct*'MOT ' ;
'ARGUMENT' tn*'CHPOINT ' ;
*
idim = 'VALEUR' 'DIME' ;
*
* Matrix A for the operator NLIN
*
numop = 1 ;
numder = idim ;
numvar = 1 ;
* The variable involved in A is the primal variable
numdat = 1 ;
numcof = 1 ;
*
discg = 'LINE' ;
* Geometrical discretisation
methgau = 'GAU7' ;
* Gauss points involved
*
A = ININLIN numop numvar numdat numcof numder ;
A . 'VAR' . 1 . 'NOMDDL' = 'MOTS' 'SCAL' ;
A . 'VAR' . 1 . 'DISC' = 'QUAF' ;
*
A . 'DAT' . 1 . 'NOMDDL' = ('EXTRAIRE' TN 'COMP') ;
A . 'DAT' . 1 . 'DISC' = disct ;
A . 'DAT' . 1 . 'VALEUR' = tn ;
A . 'COF' . 1 . 'COMPOR' = 'IDEN' ;
A . 'COF' . 1 . 'LDAT' = 'LECT' 1 ;
*
* 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' 'SCAL' ;
B . 'VAR' . 1 . 'DISC' = 'CSTE' ;
*
A . 1 . 1 . 0 = 'LECT' -1 ;
B . 1 . 1 . 0 = 'LECT' ;
*
mat = 'NLIN' discg _mt A B methgau ;
*
chpbid = 'MANUEL' 'CHPO' _mt 1 'SCAL' 1.0D0 ;
'RESPRO' ('SOMT' (mat '*' chpbid)) ;
'FINPROC' ;
*
*************************************************************************
*************************************************************************
******* END OF THE PERSONAL L PROCEDURES ********************************
*************************************************************************
*************************************************************************
*
************************************************************************
************** Begin of file *******************************************
************************************************************************
**
*************************************************
****** LIST OF SYMBOLS **************************
*************************************************
*
* \int = integral
* \div = divergence
* \der = derivative (\der{a}{b} = \frac{d a}{d b}
* \grad = gradient
* \cdot = scalar product
* ^ = power
*
* (a,b) = \int_{\Omega} a * b dV
*
* Np = test function for the divergence equations
* = form function for the pressure
* Nt = test function for the temperature
* = form function for the temperature equation
* Nv = test function for the speed equations (x,y)
* = form function for the speed
*
*
************************************************************************
************** Begin of file *******************************************
************************************************************************
NX = 6 ;
*
* IC incoming from other computations
*
* REPSAU = '/test4/becc/Dtmp/Dcaviteimplpv/' ;
REPSAU = './Dsauv/' ;
*
****************************************
**** Parameters for the computation ****
****************************************
*
* BDF1 or BDF2
*
LOGBDF1 = FAUX ;
* Tensorial/scalar/diagonal upwinding
LOGUP = FAUX ;
BETAUP = 0.0 ;
* BETAUP = 1 -> SU
* BETAUP = 0 -> No SU
LOGSCAL = FAUX ;
'MESSAGE' ;
'MESSAGE' 'Scalar numerical diffusion' ;
'MESSAGE' ;
LOGDIA = FAUX ;
* NB 'SI' LOGDIA, LOGSCAL = FAUX
'SI' LOGDIA ;
LOGSCAL = FAUX ;
'FINSI' ;
*
* Boussinesq or Low Mach?
*
LBOUS = FAUX ;
*
DISCT = 'QUAF' ;
DISCMU = DISCT ;
DISCLA = DISCT ;
*
* NB If DISCLA, DISCMU, DISCT are not 'QUAF'
* SUPG does not work
*
TFINAL = 0.25 ;
* NINT = 10 ;
EPSERR = 1.0D-20 ;
GRAPH = FAUX ;
*
MTRA = 'X' ;
NFREQ = 2 ;
* METINV = MOT 'ITER' ;
METINV = MOT 'DIRECT' ;
* PREC = 5 ;
PREC = 5 ;
NFREQP = 5 ;
ERROPRE = 1.0D-5 ;
* We compute the preconditioner each NFREQP-th iterations
* If the residial is lower than ERROPRE, we stop computing
* the preconditioner
* Number of linear iterations (for NX = 160)
NITERT = NX ;
* NITERU = NX '/' 1 ;
* NITERUT = NX '/' 4 ;
* NITERP = NX ;
* NITERT = NX ;
* NITERU = NX ;
* NITERUT = NX '*' 2;
* NITERP = NX '*' 4 ;
* CFLREF = 1.0D6 ;
* DT computed with respect to the reference speed
*
*
* Fixed options:
* - LINEar/bilinear continuos pressure
* - QUAFdratic speed
* - Name of speed (and name of speed duals variables) :
* 'UX', 'UY'
* - Name of temperature: T
* - Name of pressure: NOMPRE (Lagrange multiplicator)
*
*
'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'TRAC' MTRA 'ECHO' 1 ;
*
* Geometry dimensions
*
D1 = 1.0 ;
L1 = 12.0 '*' d1 ;
L2 = 9.0 '*' d1 ;
L3 = 7.0 '*' d1 ;
DX = d1 '/' NX ;
DXMAX = 4 '*' DX ;
*
* d1 used into Reynolds number
*
* Physical properties, reference/initial states, BC
*
R = 288. ;
gamma = 1.4 ;
g = 9.0 ;
*
PR = 0.7 ;
RE = 100. ;
TH = 267.8 ;
TC = 178.6 ;
mt = 1.0 ;
mH = 0.8 '*' mt;
mC = 0.2 '*' mt ;
*
cp = R * gamma '/' (gamma '-' 1.0) ;
*
T0 = 223.2 ;
P0 = 0.6429D5 ;
*
rho0 = P0 '/' (R '*' T0) ;
rhoC0 = P0 '/' (R '*' TC) ;
rhoH0 = P0 '/' (R '*' TH) ;
*
uC0 = mC '/' rhoC0 ;
uH0 = mH '/' rhoH0 ;
*
* mu (dynamic viscosity) and lambda (thermal diffusivity) at the
* reference state
*
mu0 = (mt '*' d1) '/' Re ;
lambda0 = mu0 '*' (cp '/' Pr) ;
*
***********************
**** Mesh *************
***********************
*
*
* P31 P32 P33 P34
*
*
* TUBE2 TUBE3 TUBE4
*
*
* P21 P22 P23 P24
*
*
*
* T
* U
* B
* 1
*
*
*
*
*
*
*
*
* P12 P13
*
*
*
*
P11 = 0.0 0.0 ;
P12 = L1 0.0 ;
P13 = (L1 '+' d1) 0.0 ;
*
P21 = 0.0 L3 ;
P22 = L1 L3 ;
P23 = (L1 '+' d1) L3 ;
P24 = (L1 '+' d1 '+' L2) L3 ;
*
P31 = 0.0 (L3 '+' d1) ;
P32 = L1 (L3 '+' d1) ;
P33 = (L1 '+' d1) (L3 '+' d1) ;
P34 = (L1 '+' d1 '+' L2) (L3 '+' d1) ;
*
* Lignes
*
P12P13 = P12 'DROIT' NX P13 ;
P22P23 = P22 'DROIT' NX P23 ;
P32P33 = P32 'DROIT' NX P33 ;
*
P12P22 = P12 'DROIT' P22 'DINI' DXMAX 'DFIN' DX ;
P13P23 = P12P22 'PLUS' (d1 0.0) ;
'ELIMINATION' (P13P23 'ET' P13 'ET' P23) (DX '/' 100.) ;
TUB1 = 'DALL' P12P13 P13P23 ('INVERSE' P22P23) ('INVERSE' P12P22)
'PLAN' ;
*
P24P34 = P24 'DROIT' NX P34 ;
P23P33 = P23 'DROIT' NX P33 ;
P22P32 = P22 'DROIT' NX P32 ;
P21P31 = P21 'DROIT' NX P31 ;
*
P24P23 = P24 'DROIT' P23 'DINI' DXMAX 'DFIN' DX ;
P34P33 = P24P23 'PLUS' (0.0 d1) ;
'ELIMINATION' (P34P33 'ET' P34 'ET' P33) (DX '/' 100.) ;
*
P21P22 = P21 'DROIT' P22 'DINI' DXMAX 'DFIN' DX ;
P31P32 = P21P22 'PLUS' (0.0 d1) ;
'ELIMINATION' (P31P32 'ET' P31 'ET' P32) (DX '/' 100.) ;
*
TUB2 = 'DALL' P21P22 P22P32 ('INVERSE' P31P32) ('INVERSE' P21P31)
'PLAN' ;
TUB3 = 'DALL' P22P23 P23P33 ('INVERSE' P32P33) ('INVERSE' P22P32)
'PLAN' ;
TUB4 = 'DALL' P24P34 P34P33 ('INVERSE' P23P33) ('INVERSE' P24P23)
'PLAN' ;
*
DOM1 = TUB1 'ET' TUB2 'ET' TUB3 'ET' TUB4 ;
*
'SI' GRAPH ;
'TRACER' DOM1 'TITRE' 'Mesh' ;
'FINSI' ;
*
* DOM1 and company are made of 'LINE' elts
*
QDOM1 = 'CHANGER' DOM1 'QUAF' ;
QP12P13 = 'CHANGER' P12P13 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP12P13 ;
QP22P23 = 'CHANGER' P22P23 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP22P23 ;
QP32P33 = 'CHANGER' P32P33 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP32P33 ;
*
QP12P22 = 'CHANGER' P12P22 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP12P22 ;
QP13P23 = 'CHANGER' P13P23 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP13P23 ;
*
QP24P34 = 'CHANGER' P24P34 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP24P34 ;
QP23P33 = 'CHANGER' P23P33 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP23P33 ;
QP22P32 = 'CHANGER' P22P32 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP22P32 ;
QP21P31 = 'CHANGER' P21P31 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP21P31 ;
*
QP24P23 = 'CHANGER' P24P23 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP24P23 ;
QP34P33 = 'CHANGER' P34P33 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP34P33 ;
*
QP21P22 = 'CHANGER' P21P22 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP21P22 ;
QP31P32 = 'CHANGER' P31P32 'QUAF' ;
'ELIMINATION' QDOM1 (DX '/' 100.) QP31P32 ;
*
* QDOM1 and friends are now QUAF
*
VOL = 'MESURE' QDOM1 ;
*
* Mesh dimension for upwinding
* HUP on DISCT (otherwise, it does not work)
*
HUP = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' ('MAXIMUM' ('PROG' DX DXMAX)) ;
*********************************************************
*** Dynamic viscosity and thermal diffusivity and the ***
*** unitary CHPOINT ***
*********************************************************
*
* CHMU and CHLAM and CHUN
*
'SI' ('EGA' DISCMU 'LINE') ;
CHMU = 'MANUEL' 'CHPO' DOM1 1 'SCAL' mu0 ;
'FINS' ;
'SI' ('EGA' DISCMU 'QUAF') ;
CHMU = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' mu0 ;
'FINS' ;
*
'SI' ('EGA' DISCLA 'LINE') ;
CHLAM = 'MANUEL' 'CHPO' DOM1 1 'SCAL' lambda0 ;
'FINS' ;
'SI' ('EGA' DISCLA 'QUAF') ;
CHLAM = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' lambda0 ;
'FINS' ;
*
'SI' ('EGA' DISCT 'LINE') ;
CHUN = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
'FINS' ;
'SI' ('EGA' DISCT 'QUAF') ;
CHUN = 'MANUEL' 'CHPO' QDOM1 1 'SCAL' 1.0 ;
'FINS' ;
**************
*** DT *******
**************
*
* DT = CFLREF '*' DX '/' (mt '/' rho0) ;
* NITER = 'ENTIER' (TFINAL '/' DT) ;
NITER = 2 ;
DT = TFINAL '/' NITER ;
*
*
******************************************
*** Fully-implicit (or projection) ? *****
******************************************
*
'REPETER' BLMET 2 ;
'SI' ('EGA' &BLMET 1) ;
LOGFI = VRAI ;
'SINON' ;
LOGFI = FAUX ;
'FINSI' ;
'SI' (LOGFI) ;
NOMPRE = 'LX' ;
NINT = 8 ;
'SINON' ;
NOMPRE = 'P' ;
NINT = 30 ;
'FINSI' ;
***************************************
********** IC / BC ********************
***************************************
*
PTHER = P0 ;
*
* IC/BC for the dynamic pressure (alaways 'LINE')
*
* Projection: imposition of PNLIM '+' tau.n on MPLIM (exit)
* imposition of \delta PNLIM = 0 on MPLIM (exit)
* imposition of the speed on the rest
* imposition of \grad \delta PNLIM . n = 0 on MPLIM
*
* Fully-implicit: the same (except thet there is not \delta PNLIM)
*
* We impose the pressure and a free boundary condition on a point
* belonging to the upper well.
MPLIM = QP22P32 'POIN' 'PROC'
((L1 '+' (d1 '/' 2)) (L3 '+' d1)) ;
MPLIM = ('MANU' 'POI1' MPLIM ) 'COULEUR' 'ROUG' ;
PNLIM = 'MANUEL' 'CHPO' MPLIM 1 NOMPRE 0.0 ;
* PNLIM not taken into account in the fully implicit case
PN = 'MANUEL' 'CHPO' DOM1 1 NOMPRE 0.0 ;
*
* BC/IC for the speed (always on 'QUAF')
*
QZEROV = ('DIFF' ('CHANGER' ('CONTOUR' QDOM1) 'POI1')
('CHANGER' 'POI1' (QP21P31 'ET' QP12P13 'ET' QP24P34)))
'COULEUR' 'VERT' ;
QZEROV = QZEROV 'ET' P21 'ET' P31 'ET' P24 'ET' P34 'ET' P12 'ET' P13;
QZEROX = (QZEROV 'ET' ('CHANGER' 'POI1' QP12P13)) 'COULEUR'
'VERT' ;
QZEROY = (QZEROV 'ET' ('CHANGER' 'POI1' QP24P34)
'ET' ('CHANGER' 'POI1' QP21P31)) 'COULEUR'
'VERT' ;
QZEROY = 'DIFF' QZEROY MPLIM ;
*
'SI' GRAPH ;
'TRACER' (QDOM1 'ET' (MPLIM 'COULEUR' ROUG))
'TITRE' 'Point for PLIM in projection.' ;
'TRACER' (QDOM1 'ET' (QZEROX 'COULEUR' ROUG))
'TITRE' 'Points in which ux is 0' ;
'TRACER' (QDOM1 'ET' (QZEROY 'COULEUR' ROUG))
'TITRE' 'Points in which uy is 0' ;
'FINSI' ;
*
* P12P13 is the cold leg
*
UN1 = 'MANUEL' 'CHPO' QZEROX 1 'UX' 0.0 ;
UN2 = 'MANUEL' 'CHPO' QZEROY 1 'UY' 0.0 ;
*
* UN1 and UN2 used to correct the BC on the speed during
* the internal iterations, as well as UNGEOMC and
* UNGEOMH
*
UNGEOMC = ('COORDONNEE' 1 QP12P13) '-' (L1 '+' (d1 '/' 2.)) ;
UNGEOMC = UNGEOMC '/' d1 ;
UNGEOMC = UNGEOMC * UNGEOMC ;
UNGEOMC = UNGEOMC '-' 0.25 ;
UNGEOMC = UNGEOMC '*' -6 ;
UN3 = UNGEOMC * uC0 ;
UN3 = 'NOMC' 'UY' UN3 ;
*
* P24P34 is the hot leg
*
UNGEOMH = ('COORDONNEE' 2 QP24P34) '-' (L3 '+' (d1 '/' 2.)) ;
UNGEOMH = UNGEOMH '/' d1 ;
UNGEOMH = UNGEOMH * UNGEOMH ;
UNGEOMH = UNGEOMH '-' 0.25 ;
UNGEOMH = UNGEOMH '*' 6 ;
UN4 = UNGEOMH * uH0 ;
UN4 = 'NOMC' 'UX' UN4 ;
*
* P21P31 is the outlet
*
UNGEOMO = ('COORDONNEE' 2 QP21P31) '-' (L3 '+' (d1 '/' 2.)) ;
UNGEOMO = UNGEOMO '/' d1 ;
UNGEOMO = UNGEOMO * UNGEOMO ;
UNGEOMO = UNGEOMO '-' 0.25 ;
UNGEOMO = UNGEOMO '*' 6 ;
*
* UN5 will be updated once the pressure derivative is
* known. For the moment UN5 = 0
*
UN5 = UNGEOMO '*' 0.0 ;
UN5 = 'NOMC' 'UX' UN5 ;
UNLIM = UN1 '+' UN2 '+' UN3 '+' UN4 '+' UN5 ;
*
UN = ('MANUEL' 'CHPO' QDOM1 2 'UX' 0.0 'UY' 0.0) '+' UNLIM ;
DUNLIM = 0.0 '*' UNLIM ;
*
* BC/IC for the temperature ('QUAF'/'LINE')
*
'SI' ('EGA' DISCT 'LINE') ;
TNLIM = 'MANUEL' 'CHPO' P12P13 1 'T' TC ;
TNLIM = TNLIM '+' ('MANUEL' 'CHPO' P24P34 1 'T' TH) ;
TN = ('MANUEL' 'CHPO' DOM1 1 'T' T0) ;
TN = TN '+' TNLIM '-' ('REDU' TN (P12P13 'ET' P24P34)) ;
DOMT = DOM1 ;
'FINS' ;
*
'SI' ('EGA' DISCT 'QUAF') ;
TNLIM = 'MANUEL' 'CHPO' QP12P13 1 'T' TC ;
TNLIM = TNLIM '+' ('MANUEL' 'CHPO' QP24P34 1 'T' TH) ;
TN = ('MANUEL' 'CHPO' QDOM1 1 'T' T0) ;
TN = TN '+' TNLIM '-' ('REDU' TN (QP12P13 'ET' QP24P34)) ;
DOMT = QDOM1 ;
'FINS' ;
*
* Plot of IC
*
'SI' GRAPH ;
VECN = 'VECTEUR' 0.5 UN ;
'TRACER' QDOM1 VECN ('CONTOUR' QDOM1) 'TITRE' 'Speed';
'TRACER' TN DOMT 'TITRE' 'Temperature' ;
'FINSI' ;
*
********************************
***** Constant Matrices ********
********************************
* (with respect to time)
*
* Mass matrix for the speed
*
* (ux,Nv)
* (uy,Nv)
*
MMASV = (MATMAS QDOM1 'UX' 'UX' 'QUAF' 'QUAF') 'ET'
(MATMAS QDOM1 'UY' 'UY' 'QUAF' 'QUAF') ;
MMASVSDT = MMASV '/' DT ;
MMASIK = 'KOPS' 'RIMA' MMASVSDT ;
*
*
* Mass matrix for the temperature
*
* (T,Nt)
*
*
MMAST = (MATMAS QDOM1 'T' 'T' DISCT DISCT) ;
MMASTSDT = MMAST '/' DT ;
*
* (\div u, Np)
*
MDIV = MATDIV QDOM1 NOMPRE 'LINE' 'QUAF' ;
*
* Matrix for the gravity source term in the speed equation
*
* (g (rho0/rho '-' 1), Nv) j
*
* where * rho0 = P0 '/' (R '*' T0)
* * (1/rho) = (R '*' TN) '/' PTHER
*
* N.B. rho0/rho is a CHPOINT with one component ('T'), with the same
* discretisation as TN (disct).
*
* nomp nomd discp discd
MSOUGRA = MATMAS QDOM1 'T' 'UY' DISCT 'QUAF' ;
*
* Matrix for DPTHSDT source term in the temperature equation
*
* ( (1/(rho cp)) \der{PTHER}{t}, Nt)
* where * Nt is the test function for the temperature
* * \der{PTHER}{t} ia a FLOTTANT
* * (1/(rho cp)) = (R '*' TN) '/' (cp PTHER)
* is a CHPOINT with one component ('T'),
* with the same discretisation as TN (disct)
*
MSOUTEM = MATMAS QDOM1 'T' 'T' DISCT DISCT ;
*
* In the divergence equation, we have
*
* \div(u) = (-1/(gamma * PTHER)) '*' \der{PTHER}{t}
* '+'
* ((gamma '-' 1) '/' (gamma '*' PTHER)) \div(lambda \grad(T))
*
* Discretisation is perform by the scalar product of the equation with
* Np
*
* Matrix to compute the contribution of \div(lambda \grad(T)) to
* the divergence of the speed.
*
* Surface contributions in QP12P13 (cold inlet), QP24P34 (cold inlet).
* At the oulet, we impone natural boundary conditions. Elsewhere, we
* have adiabatic surfaces
*
*
* \int_s lambda (grad T) Np \cdot n dS
* Here (lambda grad T) is the primal variable
*
*'DEBPROC' MATFLU ;
*'ARGUMENT' _st*'MAILLAGE' ;
*'ARGUMENT' discp*'MOT ' ;
*'ARGUMENT' nomd*'MOT ' ;
*'ARGUMENT' discd*'MOT ' ;
*'ARGUMENT' discco*'MOT ' ;
*'ARGUMENT' coef*'CHPOINT' ;
*'ARGUMENT' discn*'MOT ' ;
*'ARGUMENT' nx*'CHPOINT' ;
*'ARGUMENT' ny*'CHPOINT' ;
*
*
NORMX = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 0.0 ;
NORMY = 'MANUEL' 'CHPO' DOM1 1 'SCAL' -1.0 ;
COEFBID = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
MSMDCOLD = MATFLU QP12P13 'QUAF' NOMPRE 'LINE' 'LINE' COEFBID
'LINE' NORMX NORMY ;
NORMX = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
NORMY = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 0.0 ;
COEFBID = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
MSMDHOT = MATFLU QP24P34 'QUAF' NOMPRE 'LINE' 'LINE' COEFBID
'LINE' NORMX NORMY ;
*
* - (\int lambda (grad T) \cdot (grad Np) dV) =
* -1. '*' MATLAP
*
*'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' ;
*
*
COEFBID = 'MANUEL' 'CHPO' DOM1 'SCAL' 1.0 ;
MSMDVO = (MATLAP QDOM1 'T' DISCT NOMPRE 'LINE' DISCLA CHLAM
'LINE' COEFBID) '*' -1. ;
*
* Matrix to compute the contribution of \div(lambda \grad(T)) to
* the divergence condition. This matrix has as primal variable both
* t ('T') and grad T ('UX', 'UY')!
*
MSMD = MSMDCOLD 'ET' MSMDHOT 'ET' MSMDVO ;
*
* Matrix to compute the contribution of
*
* (-1/(gamma * PTHER)) '*' \der{PTHER}{t}
*
MSMDIV = MATMAS QDOM1 'SCAL' NOMPRE 'LINE' 'LINE' ;
*
* Test: DIV U = -u0 * LINJ
*
'SI' VRAI ;
uinj = 'SOMT' (MDIV '*' UN) ;
uinj1 = d1 * (Uc0 '+' Uh0) ;
ERRO = (UINJ1 '+' uinj) 'ABS' ;
'SI' (ERRO > 1.0D-10) ;
'ERREUR' 21 ;
'FINSI' ;
'FINS' ;
*
*************************************************************************
* Temporal loop *********************************************************
*************************************************************************
*
TAB1 = 'TABLE' ;
TAB1.'TITRE'= TABLE ;
TAB1.1='TIRR ';
TAB1.2='TIRC ';
TAB1.3='REGU ';
TAB1.'TITRE' . 1 = MOT 'Mass cons.';
TAB1.'TITRE' . 2 = MOT 'Div. cons.';
TAB1.'TITRE' . 3 = MOT 'Analyt.';
TITP = 'CHAINE' 'Pressure from BC.' ;
*
TPS = 0.0 ;
*
LIT = 'PROG' ;
LTPS = 'PROG' ;
LER = 'PROG' ;
LERT = 'PROG' ;
LERP = 'PROG' ;
LPTHER = 'PROG' ;
ERRO = 1.0D10 ;
ERROT = 1.0D10 ;
*
TABRESTN = 'TABLE' ;
TABRESTN . 'XINIT' = TN ;
TABRESUN = 'TABLE' ;
'SI' LOGFI ;
TABRESUN . 'XINIT' = UN '+' PN ;
'SINON' ;
TABRESUN . 'XINIT' = UN ;
TABRESDU = 'TABLE' ;
TABRESDU . 'XINIT' = UN '*' 0.0 ;
TABRESDP = 'TABLE' ;
TABRESDP . 'XINIT' = (PN '*' 0.0) ;
'FINSI' ;
*
INDITE = 0 ;
*
UNM = 'COPIER' UN ;
TNM = 'COPIER' TN ;
*
*****************************************
*********** TIME LOOP *******************
*****************************************
*
'REPETER' BL1 NITER ;
UNM2 = 'COPIER' UNM ;
TNM2 = 'COPIER' TNM ;
UNM = 'COPIER' UN ;
TNM = 'COPIER' TN ;
PTHERM2 = PTHERM ;
PTHERM = PTHER ;
*
**** Time step
*
TPS = TPS '+' DT ;
*
**** Pressure and its time variation
*
PTHER = P0 ;
'SI' (('EGA' &BL1 1) 'OU' LOGBDF1) ;
DPTHSDT = (PTHER '-' PTHERM) '/' DT ;
'SINON' ;
DPTHSDT = ((3 '*' (PTHER '-' PTHERM)) '-'
(PTHERM '-' PTHERM2)) '/' (2. '*' DT) ;
'FINSI' ;
*
*
**** Speeds vary during the time iterations
*
rhoC = PTHER '/' (R '*' TC) ;
rhoH = PTHER '/' (R '*' TH) ;
*
uC = mC '/' rhoC ;
uH = mH '/' rhoH ;
*
UN3 = UNGEOMC * uc ;
UN3 = 'NOMC' 'UY' UN3 ;
UN4 = UNGEOMH * uh ;
UN4 = 'NOMC' 'UX' UN4 ;
UNLIM = UN1 '+' UN2 '+' UN3 '+' UN4 ;
*
*** We have to correct the speed also at the outlet
*
UNLIM1 = 'REDU' UN ('EXTRAIRE' UNLIM 'MAILLAGE') ;
UN = UN '+' (UNLIM '-' UNLIM1) ;
*
**** INTERNAL ITERATIONS ****************
*
'REPETER' BLINT NINT ;
*
INDITE = INDITE '+' 1 ;
*
************************************************************************
**** Divergence constraint and correction of the speed at the outlet ***
************************************************************************
*
* We compute qdiv = (div(u),Npj), Npj being the form function
* for the pressure.
*
* First of all, let us evaluate
*
* coef (\div (lambda \grad T), Npj)
*
'SI' LBOUS ;
coef = 0.0 ;
'SINON' ;
coef = (gamma '-' 1) '/' (gamma '*' PTHER) ;
'FINSI' ;
*
* Surface integral contribution
* CHPFLU
* 'ARGUMENT' _mt*'MAILLAGE' ;
* 'ARGUMENT' disct*'MOT ' ;
* 'ARGUMENT' T*'CHPOINT' ;
* 'ARGUMENT' discd*'MOT ' ;
* 'ARGUMENT' discla*'MOT ' ;
* 'ARGUMENT' lambda*'CHPOINT' ;
*
lamgt = CHPFLU QDOM1 DISCT TN 'QUAF' DISCLA CHLAM ;
qdiv = MSMD '*' ((coef '*' TN) '+' (coef '*' lamgt)) ;
*
* Since we are in a open domain, we must know
* 1/(gamma P) (dP/dt)
*
chcoef = 'MANUEL' 'CHPO' DOM1 1 'SCAL'
(-1 '*' DPTHSDT '/' (gamma '*' PTHER)) ;
qdiv = qdiv '+' (MSMDIV '*' chcoef) ;
*
* Finally, we correct the speed at the outlet in order to satisfy
* \sum divu
*
cacca = 'SOMT' (qdiv '-' (MDIV '*' UN)) ;
cacca = cacca '/' d1 ;
* cacca is the increment of average speed which exit from QP12P13
UN5 = UNGEOMO * cacca ;
UN5 = 'NOMC' 'UX' UN5 ;
UN = UN '+' UN5 ;
UNLIM = UNLIM '+' UN5 ;
*
* Test: we check whether (MDIV '*' UN) = qdiv
*
cacca = 'SOMT' (qdiv '-' (MDIV '*' UN)) ;
'SI' (('ABS' cacca) > 1.0D-12) ;
VECN = 'VECTEUR' 0.5 UN ;
'TRACER' QDOM1 VECN ('CONTOUR' QDOM1) 'TITRE' 'Speed';
'ERREUR' 21 ;
'FINSI' ;
*
* We save the values of PTHER, UN, TN, PN to compute
* the convergence error
*
LOGER = ('EGA' &BLINT 1) 'OU'
('EGA' ((&BLINT '/' NFREQ) '*' NFREQ) &BLINT) ;
*
'SI' LOGER ;
UNL = 'COPIER' UN ;
PNL = 'COPIER' PN ;
TNL = 'COPIER' TN ;
'FINSI' ;
*
UNTX = ('EXCO' 'UX' UN) ;
UNTY = ('EXCO' 'UY' UN) ;
*
**********************************************************************
*** Tensorial numerical viscosity ************************************
**********************************************************************
*
'SI' LOGUP ;
UMOD = (UNTX '*' UNTX) '+' (UNTY '*' UNTY) ;
UMOD = UMOD '**' 0.5 ;
* We add epsilon to avoid UMOD = 0
UMODEP = UMOD '+' (1.0D-6 '*' u0) ;
*
* We compute \alpha = \frac{\lambda}{\rho c_p} =
* \frac{(\gamma '-' 1) T}{\gamma P} \lambda
*
*
ALPHA = ((gamma '-' 1) '/' (gamma '*' PTHER)) '*'
('NOMC' TN 'SCAL') '*' CHLAM ;
* The reciprocal of Peclet
USPE = ALPHA '/' (UMODEP '*' HUP) ;
* The reciprocal of Reynolds
USRE = Pr '*' USPE ;
*
GT = -1 '*' (USPE '-' 0.5) ;
GT = GT '*' ('MASQUE' GT 'SUPERIEUR' 0.0) ;
GT = GT '*' (HUP '/' UMODEP) ;
*
GV = -1 '*' (USRE '-' 0.5) ;
GV = GV '*' ('MASQUE' GV 'SUPERIEUR' 0.0) ;
GV = GV '*' (HUP '/' UMODEP) ;
*
* The tensorial matrix
*
NUXX = UNTX '*' UNTX ;
NUYY = UNTY '*' UNTY ;
NUXY = UNTX '*' UNTY ;
NUYX = NUXY ;
*
'SI' LOGSCAL ;
NU1 = (NUXX '+' NUXY) ;
NU2 = (NUXX '-' NUXY) ;
NU1 = NU1 '+' (NU2 'ABS') ;
NU = 0.5 '*' NU1 ;
NU1 = (NU '+' NUYX) ;
NU2 = (NU '-' NUYX) ;
NU1 = NU1 '+' (NU2 'ABS') ;
NU = 0.5 '*' NU1 ;
NU1 = (NU '+' NUYY) ;
NU2 = (NU '-' NUYY) ;
NU1 = NU1 '+' (NU2 'ABS') ;
NU = 0.5 '*' NU1 ;
NUXX = NU ;
NUXY = 0 '*' NU ;
NUYX = 0 '*' NU ;
NUYY = NU ;
'FINSI' ;
'SI' LOGDIA ;
NUXY = 0 '*' NUXY ;
NUYX = 0 '*' NUYX ;
'FINSI' ;
*
* Diffusion matrix for the temperature
*
MLAPTD = MATLAPT QDOM1 'T' DISCT 'T' DISCT 'QUAF'
(GT '*' NUXX) (GT '*' NUXY) (GT '*' NUXY) (GT '*' NUYY)
'QUAF' (CHUN '*' BETAUP) ;
MLAPVD = (MATLAPT QDOM1 'UX' 'QUAF' 'UX' 'QUAF' 'QUAF'
(GV '*' NUXX) (GV '*' NUXY) (GV '*' NUXY) (GV '*' NUYY)
'QUAF' (CHUN '*' BETAUP)) 'ET'
(MATLAPT QDOM1 'UY' 'QUAF' 'UY' 'QUAF' 'QUAF'
(GV '*' NUXX) (GV '*' NUXY) (GV '*' NUXY) (GV '*' NUYY)
'QUAF' (CHUN '*' BETAUP)) ;
'FINSI' ;
*
*******************
*** Temperature ***
*******************
*
* ((u \cdot \grad) T^{n+1}, Nt), LHS
*
MKONT = MATKT QDOM1 'QUAF' UNTX UNTY 'T' DISCT ;
*
* (-coef \div (alpha \grad T), Nt), LHS
* Contribution of the integral on the volume, LHS
* The contribution of the integral of surface have
* to be specified (if it is not 0)
*
* coef = 1 '/' (rho '*' cp) =
* = ((gamma '-' 1) '*' TN) '/' (gamma '*' PTHER)
*
* NB coef = CHPOINT with one component ('SCAL') on DISCT
*
DISCOF = DISCT ;
'SI' LBOUS ;
*
* In Boussinesq approx coef = 1 '/' (rho0 '*' cp)
*
coef = ((gamma '-' 1) '/' (gamma '*' p0)) '*' T0 ;
COEFL = 'MANUEL' 'CHPO' DOMT 1 'SCAL' coef ;
'SINON' ;
COEFL = ((gamma '-' 1) '/' (gamma '*' pther)) '*'
('NOMC' TN 'SCAL') ;
'FINSI' ;
MLAPT = (MATLAP QDOM1 'T' DISCT 'T' DISCT DISCLA CHLAM DISCOF
COEFL) ;
*
* Source term
*
* (coef * Tn,Lt) = coef (Tn,Lt) RHS
*
* where coef = ((gamma '-' 1) '*' DPTHSDT) '/' (gamma '*' PTHER)
* =FLOTTANT
*
DISCOF = DISCT ;
'SI' LBOUS ;
COEFL = 'MANUEL' 'CHPO' DOMT 1 'T' 0.0 ;
'SINON' ;
COEFL = (((gamma '-' 1) '*' DPTHSDT) '/' (gamma '*' PTHER))
'*' TN ;
'FINSI' ;
qtt = MSOUTEM '*' COEFL ;
'SI' (('EGA' &BL1 1) 'OU' LOGBDF1) ;
*
*** First time iteration: BDF1
*
* ((1 '/' dt) T^{n}, Nt), RHS
*
qtt = qtt '+' (MMASTSDT '*' TNM) ;
mtott = mmastsdt 'ET' mkont 'ET' mlapt ;
'SINON' ;
qtt = qtt '+' (MMASTSDT '*' ((2 '*' TNM) '-' (0.5 '*' TNM2))) ;
mtott = (1.5 '*' mmastsdt) 'ET' mkont 'ET' mlapt ;
'FINSI' ;
*
* Temperature: LHS, RHS:
*
'MESSAGE' ;
'MESSAGE' 'Computation of T' ;
'SI' LOGUP ;
mtott = mtott 'ET' mlaptd ;
'FINSI' ;
smbt = qtt ;
mtotik = 'KOPS' 'RIMA' mtott ;
LOGPREC = (&BLINT 'EGA' 1) 'OU' (INDITE 'EGA' 1) 'OU'
('EGA' METINV 'DIRECT') ;
'SI' LOGPREC ;
mprect = mtotik ;
'SINON' ;
'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
(ERROT > ERROPRE)) ;
mprect = mtotik ;
'MESSAGE' 'Computation of preconditioner' ;
'FINSI' ;
'FINSI' ;
TN = RESOUP FAUX mtotik mprect smbt (TNLIM) 1.0D-60 NITERT
PREC TABRESTN METINV ;
*
*************
*** Speed ***
*************
*
* ((1 '/' dt) ux^{n}, Nv), RHS
* ((1 '/' dt) uy^{n}, Nv), RHS
*
'SI' (('EGA' &BL1 1) 'OU' LOGBDF1) ;
*
*** First time iteration: BDF1
*
qtv = (MMASVSDT '*' UNM) ;
mtotv = MMASVSDT ;
'SINON' ;
qtv = (MMASVSDT '*' ((2 '*' UNM) '-' (0.5 '*' UNM2))) ;
mtotv = 1.5 '*' MMASVSDT ;
'FINSI' ;
*
*
* ((u \cdot \grad) ux^{n+1}, Nv), LHS
* ((u \cdot \grad) uy^{n+1}, Nv), LHS
*
MKONV = (MATKV QDOM1 'QUAF' 'QUAF' UNTX UNTY) ;
*
* (-coef \div (nu \tau), Nv), LHS
* contribution of the integral on the volume, LHS
* Contribution of the integral of surface have to be specified
* (if it is not 0).
*
*
* coef = 1/rho = (R '/' PTHER) '*' TN
*
DISCOF = DISCT ;
'SI' LBOUS ;
*
* coef = RT0/P0
*
coef = (R '/' P0) '*' T0 ;
COEFL = 'MANUEL' 'CHPO' DOMT 1 'SCAL' coef ;
'SINON' ;
COEFL = (R '/' PTHER) '*' ('NOMC' TN 'SCAL') ;
'FINSI' ;
MLAPV = MATTAU QDOM1 'QUAF' DISCMU DISCOF CHMU COEFL ;
*
* (coef \dep{p}{x}, Nv), LHS (RHS in projection)
* (coef \dep{p}{y}, Nv), LHS (RHS in projection)
* contribution of the integral on the volume, RHS
* Contribution of the integral of surface have to be specified (if it
* is not 0).
*
* NB coef of gradP and div tau are the same
* But gradP is on the RHS, tau in the LHS
*
MPRE = MATPRE QDOM1 NOMPRE 'LINE' 'QUAF' DISCOF COEFL ;
'SI' ('NON' LOGFI) ;
qtv = qtv '+' (MPRE '*' (-1 '*' PN)) ;
'FINSI' ;
*
* (g ((rho_0 '/' rho) '-' 1) , Nv) j, RHS
*
'SI' LBOUS ;
*
* (g (T/T0 '-' 1) , Nv) j, RHS
*
CHPO1 = ((TN '/' T0) '-' 1.) '*' g ;
'SINON' ;
CHPO1 = ((TN '*' (P0 '/' (T0'*' PTHER))) '-' 1.) '*' g ;
'FINSI' ;
qtv = qtv '+' (MSOUGRA '*' CHPO1) ;
*
* Speed: LHS, RHS
*
'SI' LOGFI ;
'MESSAGE' ;
'MESSAGE' 'Computation of U' ;
mtotv = mtotv 'ET' MKONV 'ET' MLAPV 'ET' MPRE 'ET' MDIV ;
'SI' LOGUP ;
mtotv = mtotv 'ET' mlapvd ;
'FINSI' ;
*
smbv = qtv '+' qdiv ;
mtotik = 'KOPS' 'RIMA' mtotv ;
* mtotik = 'ET' mtotik mdiv2 ;
'SI' LOGPREC ;
mprecufi = mtotik ;
'SINON' ;
'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
(ERRO > ERROPRE)) ;
mprecufi = mtotik ;
'MESSAGE' 'Computation of preconditioner' ;
'FINSI' ;
'FINSI' ;
'SI' ('EGA' ('TYPE' NITERU) ENTIER) ;
UNP = RESOUP FAUX mtotik mprecufi smbv UNLIM 1.0D-60 NITERU
PREC TABRESUN METINV ;
'SINON' ;
UNP = RESOUP FAUX mtotik mprecufi smbv UNLIM 1.0D-16 200
PREC TABRESUN METINV ;
'FINSI' ;
UN = 'EXCO' ('MOTS' 'UX' 'UY') UNP ('MOTS' 'UX' 'UY') ;
PN = 'EXCO' ('MOTS' NOMPRE) UNP ('MOTS' NOMPRE) ;
'SINON' ;
*
**** Speed without divergence constraint
*
'MESSAGE' ;
'MESSAGE' 'Computation of Util' ;
mtotv = mtotv 'ET' MKONV 'ET' MLAPV ;
'SI' LOGUP ;
mtotv = mtotv 'ET' mlapvd ;
'FINSI' ;
smbv = qtv ;
mtotik = 'KOPS' 'RIMA' mtotv ;
'SI' LOGPREC ;
mprecut = mtotik ;
'SINON' ;
'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
(ERRO > ERROPRE)) ;
mprecut = mtotik ;
'MESSAGE' 'Computation of preconditioner' ;
'FINSI' ;
'FINSI' ;
'SI' ('EGA' ('TYPE' NITERUT) ENTIER) ;
UNTIL = RESOUP FAUX mtotik mprecut smbv UNLIM 1.0D-60
NITERUT PREC TABRESUN METINV ;
'SINON' ;
UNTIL = RESOUP FAUX mtotik mprecut smbv UNLIM 1.0D-16
200 PREC TABRESUN METINV ;
'FINSI' ;
*
**** Dynamic pressure
*
* First of all, we have to compute the laplacian
* of the dynamic pressure increment
*
* (-coef \div coef1 \grad P, Np), LHS
*
* coef = -1
* coef1 = 1/rho = (R '/' PTHER) '*' TN
*
COEF = 'MANUEL' 'CHPO' DOM1 1 'SCAL' -1.0 ;
DISCOF = DISCT ;
'SI' LBOUS ;
*
* coef1 = RT0/P0
*
coef1 = (R '/' P0) '*' T0 ;
COEF1 = 'MANUEL' 'CHPO' DOMT 1 'SCAL' coef1 ;
'SINON' ;
COEF1 = (R '/' PTHER) '*' ('NOMC' TN 'SCAL') ;
'FINSI' ;
MLAPP = (MATLAP QDOM1 NOMPRE 'LINE' NOMPRE 'LINE' DISCOF COEF1
'LINE' COEF) ;
'MESSAGE' ;
'MESSAGE' 'Computation of deltaP' ;
qpv = MDIV '*' UNTIL ;
qpv = (qpv '-' qdiv) '/' DT ;
'SI' ('NON' (('EGA' &BL1 1) 'OU' LOGBDF1)) ;
*
********** BDF2
*
qpv = qpv '*' 1.5 ;
'FINSI' ;
*
mtotik = 'KOPS' 'RIMA' MLAPP ;
'SI' LOGPREC ;
mprecp = mtotik ;
'SINON' ;
'SI' (('EGA' ((INDITE '/' NFREQP) '*' NFREQP) INDITE) 'ET'
(ERRO > ERROPRE)) ;
mprecp = mtotik ;
'MESSAGE' 'Computation of preconditioner' ;
'FINSI' ;
'FINSI' ;
'SI' ('EGA' ('TYPE' NITERP) ENTIER) ;
DELTAP = RESOUP FAUX mtotik mprecp qpv PNLIM 1.0D-60 NITERP
PREC TABRESDP METINV ;
'SINON' ;
DELTAP = RESOUP FAUX mtotik mprecp qpv PNLIM 1.0D-16 200
PREC TABRESDP METINV ;
'FINSI' ;
PN = PN '+' DELTAP ;
*
**** Speed (with divergence constraint)
*
'MESSAGE' ;
'MESSAGE' 'Computation of deltaU' ;
qtv = MPRE '*' DELTAP ;
DELTAU = RESOUP FAUX MMASIK MMASIK qtv DUNLIM 1.0D-16 200
PREC TABRESDU METINV ;
UN = UNTIL '-' DELTAU ;
'FINSI' ;
*
***********************************************
**** Post-treatment and convergence check *****
***********************************************
*
'SI' LOGER ;
ERRO = 'MAXIMUM' (UN '-' UNL) 'ABS' ;
ERROT = 'MAXIMUM' (TN '-' TNL) 'ABS' ;
ERROP = 'MAXIMUM' (PN '-' PNL) 'ABS' ;
'MESSAGE'
'***********************************************************' ;
'MESSAGE' ('CHAINE' 'TEMPORAL ITER = ' &BL1
', TPS = ' TPS ', ERRO = ' ERRO
', INTIT = ' &BLINT) ;
'MESSAGE'
'***********************************************************' ;
LIT = LIT 'ET' ('PROG' INDITE) ;
L10 = 'LOG' 10. ;
LER = LER 'ET' ('PROG' (('LOG' (ERRO '+' 1.0D-20)) '/' L10)) ;
* EVER = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LER ;
L10 = 'LOG' 10. ;
LERT = LERT 'ET'
('PROG' (('LOG' (ERROT '+' 1.0D-20)) '/' L10));
LERP = LERP 'ET'
('PROG' (('LOG' (ERROP '+' 1.0D-20)) '/' L10));
EVERT = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LERT ;
'SI' GRAPH ;
TIT = 'CHAINE' 'TPS =' TPS ' , IT = ' &BL1 ' , INTIT = '
&BLINT ;
'SI' (((INDITE '/' (2 '*' NFREQ)) '*'
(2 '*' NFREQ)) 'EGA' INDITE) ;
'DESSIN' EVERT 'TITR'
('CHAINE' 'Convergence history (T). ' TIT)
'YBOR' -16. 3 'NCLK' ;
* 'DESSIN' EVER 'TITR' 'Convergence history (v)' 'NCLK' ;
'SINON' ;
VUN = 'VECTEUR' 0.1 UN ;
'TRACER' TN DOM1 VUN 'TITRE' TIT 'NCLK' ;
* EVP = 'EVOL' 'MANU' 'it' LTPS 'Pther' LPMA ;
* EVP2 = 'EVOL' 'MANU' 'it' LTPS 'Pther' LPCO ;
* 'DESSIN' (EVP 'ET' EVP2) 'TITR' TITP 'LEGE' TAB1 'NCLK' ;
'FINSI' ;
'FINSI' ;
'FINSI' ;
'SI' (ERROT < EPSERR) ;
'QUITTER' BLINT ;
'FINSI' ;
*
'FIN' BLINT ;
*
***** END OF INTERNAL ITERATIONS
*
LTPS = LTPS 'ET' ('PROG' TPS) ;
LPTHER = LPTHER 'ET' ('PROG' PTHER) ;
'FIN' BL1 ;
*
*
*************************************************************************
* End of temporal loop **************************************************
*************************************************************************
*
*
* 'OUBLIER' RIGIDITE object
*
'OUBLIER' MFLU ;
'OUBLIER' MTOTV ;
'OUBLIER' MPRE ;
'OUBLIER' MLAPV ;
'OUBLIER' MKONV ;
'OUBLIER' MTOTT ;
'OUBLIER' MLAPT ;
'OUBLIER' MKONT ;
'OUBLIER' MSMDIV ;
'OUBLIER' MSMD ;
'OUBLIER' MSMD5 ;
'OUBLIER' MSMD4 ;
'OUBLIER' MSMD3 ;
'OUBLIER' MSMD2 ;
'OUBLIER' MSMD1 ;
'OUBLIER' MSOUTEM ;
'OUBLIER' MSOUGRA ;
'OUBLIER' MDIV ;
'OUBLIER' MMASTSDT;
'OUBLIER' MMAST ;
'OUBLIER' MMASVSDT;
'OUBLIER' MMASV ;
'OUBLIER' MLAPTD ;
'OUBLIER' MLAPVD ;
*
* 'OUBLIER' MATRIK object
*
'OUBLIER' MTOTIK ;
'OUBLIER' MMASIK ;
'OUBLIER' MPRECT ;
'OUBLIER' MPRECUT ;
'OUBLIER' MPRECP ;
'OUBLIER' MPRECUFI ;
*
'MENAGE' ;
*
*
* 'SI' LBOUS ;
* 'OPTION' 'SAUVER' ('CHAINE' REPSAU 'bous'
* 'g' ('ENTIER' g) 'Re100' NX '.sauv') ;
* 'SINON' ;
* 'OPTION' 'SAUVER' ('CHAINE' REPSAU
* 'g' ('ENTIER' g) 'Re100' NX '.sauv') ;
* 'FINSI' ;
* 'SAUVER' ;
*
* Test on the convergence error
*
'SI' LOGFI ;
TNFI = 'COPIER' TN ;
aa = 'IPOL' (NINT '*' 0.9999) ('EXTR' evert 'ABSC')
('EXTR' evert 'ORDO') ;
'SI' (aa > -5.) ;
'MESSAGE' 'FI: convergerce error inacceptable' ;
'MESSAGE' ('CHAINE' 'erro = ' aa) ;
'ERREUR' 5 ;
'FINSI' ;
aa = 'IPOL' (2 '*' NINT '*' 0.9999) ('EXTR' evert 'ABSC')
('EXTR' evert 'ORDO') ;
'SI' (aa > -5.) ;
'MESSAGE' 'FI: convergerce error inacceptable' ;
'MESSAGE' ('CHAINE' 'erro = ' aa) ;
'ERREUR' 5 ;
'FINSI' ;
'SINON' ;
TNPRO = 'COPIER' TN ;
aa = 'IPOL' (NINT '*' 0.9999) ('EXTR' evert 'ABSC')
('EXTR' evert 'ORDO') ;
'SI' (aa > -4.) ;
'MESSAGE' 'PROJ: convergerce error inacceptable' ;
'MESSAGE' ('CHAINE' 'erro = ' aa) ;
'ERREUR' 5 ;
'FINSI' ;
aa = 'IPOL' (2 '*' NINT '*' 0.9999) ('EXTR' evert 'ABSC')
('EXTR' evert 'ORDO') ;
'SI' (aa > -4.) ;
'MESSAGE' 'PROJ: convergerce error inacceptable' ;
'MESSAGE' ('CHAINE' 'erro = ' aa) ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
'FIN' BLMET ;
*
* Test on the results obatined via the projection and the fully implicit
* approach.
*
ERRO = 'MAXI' ((TNPRO - TNFI) '/' T0) 'ABS' ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' 'FI and PROJ give different results!' ;
'MESSAGE' ('CHAINE' 'erro = ' erro) ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer