* fichier : nlin_cavity_HP.dgibi
*************************************************************************
*
* NOM : nlin_cavity_HP.dgibi
*
* DESCRIPTION : We compute the flow governed by the Navier-Stokes
* equations, in a square cavity with high temperature
* difference.
* Both Boussinesq approximation and Low Mach number
* approximation are considered.
* See also technical report SFME/LTMF/RT/O4-031/A
*
* 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
*
* PMASS = thermodyn. pressure via the integration of the EOF (by
* imposing that the integral of the density keeps contant)
*
* PCOMP = thermodyn. pressure via the divergence constraint.
*
* PTHER = PMASS or PCOMP
*
*
*************************************************************************
*************************************************************************
******* 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 *******************************************
************************************************************************
*
NXS2 = 10 ;
'REPETER' BLPAR 2 ;
*
* IC incoming from other computations?
*
* REPSAU = '/test4/becc/Dtmp/Dcaviteimplpv/' ;
* REPSAU = './' ;
VARREST = FAUX ;
'SI' VARREST ;
'OPTION' 'RESTITUER' ('CHAINE' 'file' NXS2 '.sauv') ;
* 'OPTION' 'RESTITUER' ('CHAINE' 'resu' NXS2 '.sauv') ;
'RESTITUER' ;
* We give PMASS, PCOMP, UN, PN, TN
*
* Modif. to use old results
*
'SI' ('NEG' ('TYPE' PMASS) 'FLOTTANT') ;
PMASS = PTHER ;
PMASSM = PTHER ;
PCOMP = PTHER1;
PCOMPM = PTHER1;
'FINSI' ;
'FINSI' ;
*
****************************************
**** Parameters for the computation ****
****************************************
*
* Fully-implicit (or projection) ?
LOGFI = VRAI ;
'SI' (LOGFI) ;
NOMPRE = 'LX' ;
'SINON' ;
NOMPRE = 'P' ;
'FINSI' ;
'SI' VARREST ;
*
* Modif. to use old results
*
PN = 'NOMC' NOMPRE PN ;
*
'FINSI' ;
*
* Pressure from the compatibility condition.
*
LOGCOMP = FAUX ;
*
*
* Boussinesq or Low Mach?
*
'SI' ('EGA' &BLPAR 1) ;
LBOUS = FAUX ;
'SINON' ;
LBOUS = VRAI ;
'FINSI' ;
* \lambda and \mu as function of T?
VARPRO = FAUX ;
*
DISCT = 'QUAF' ;
DISCMU = DISCT ;
DISCLA = DISCT ;
NITMAX = 1200 ;
GRAPH = FAUX ;
MTRA = 'X' ;
NFREQ = 10 ;
* METINV = MOT 'ITER' ;
METINV = MOT 'DIRECT' ;
* PREC = 5 ;
PREC = 5 ;
NFREQP = 40 ;
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 NXS2 = 160)
NITERT = NXS2 '/' 2 ;
NITERU = NXS2 '/' 1 ;
* NITERUT = NXS2 '/' 4 ;
* NITERP = NXS2 ;
* NITERT = NXS2 ;
* NITERU = NXS2 ;
* NITERUT = NXS2 '*' 2;
* NITERP = NXS2 '*' 4 ;
CFLREF = 50. ;
* If I take a constant DT, I can gain CPU time in computing the (real)
* speed -> DT computed with respect to the reference speed
* NB. Even with a variable DT, we have to compute an initial DT in order
* to avoid DT=\infty with U=0.
*
*
* Fixed options:
* - LINEar/bilinear continuus 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 ;
*
* Physical properties
*
R = 287. ;
gamma = 1.4 ;
g = 9.81 ;
PR = 0.71 ;
*
cp = R * gamma '/' (gamma '-' 1.0) ;
*
* Reference/initial state
*
T0 = 600. ;
P0 = 101325. ;
rho0 = P0 '/' (R '*' T0) ;
*
* mu (dynamic viscosity) and lambda (thermal diffusivity) at the
* reference state
*
'SI' varpro ;
esse = 110.5 ;
tstar = 273. ;
mustar = 1.68D-5 ;
'FINSI' ;
mu0 = 2.95456D-5 ;
lambda0 = mu0 '*' (cp '/' Pr) ;
*
* Hot/cold temperature
*
* epsilon = 2 \Delta T /T_m
* T_h = T_m (1 '+' epsilon)
* T_c = T_m (1 '-' epsilon)
*
epsilon = 0.6 ;
TH = T0 '*' (1 '+' epsilon) ;
TC = T0 '*' (1 '-' epsilon) ;
* L as function of Rayleight (Grashof) number
RAY = 1.0D6 ;
GR = RAY '/' PR ;
* GR = UREF^2 * L^2 * rho0^2 '/' mu0
UREF2L2 = (mu0 '*' mu0 '*' GR) '/' (rho0 * rho0);
UREF2SL = g '*' (TH '-' TC) '/' T0 ;
L3 = UREF2L2 '/' UREF2SL ;
L = L3 '**' (1. '/' 3.) ;
*
* Test
*
UREF = ((g '*' L '*' (TH '-' TC) '/' T0)) '**' 0.5 ;
RE = rho0 '*' UREF '*' L '/' mu0 ;
GR1 = RE '*' RE ;
'SI' (((GR '-' GR1) '/' GR) > 1.0D-4) ;
'ERREUR' 21 ;
'FINSI' ;
*
***********************
**** Mesh *************
***********************
*
* A4 A3
*
*
*
*
* A1 A2
*
A1 = 0.0 0.0 ;
A12 = (L '/' 2.) 0.0 ;
A2 = L 0.0 ;
A23 = L (L '/' 2.) ;
A3 = L L ;
A34 = (L '/' 2.) L ;
A4 = 0.0 L ;
A41 = 0.0 (L '/' 2.) ;
*
DX = L '/' (2 '*' NXS2) ;
DXMINC = DX '/' 4 ;
DXMINH = DX '/' 2. ;
DXMAX = DX '*' 2. ;
DY = DX ;
'SI' LBOUS ;
DXMINH = DX '/' 3 ;
DXMINC = DX '/' 3 ;
'FINSI' ;
*
A1A2 = (A1 'DROIT' A12 'DINI' DXMINC 'DFIN' DXMAX) 'ET'
(A12 'DROIT' A2 'DINI' DXMAX 'DFIN' DXMINH) ;
A2A3 = (A2 'DROIT' A23 'DINI' DXMINH 'DFIN' DXMAX) 'ET'
(A23 'DROIT' A3 'DINI' DXMAX 'DFIN' DXMINH) ;
A3A4 = (A3 'DROIT' A34 'DINI' DXMINH 'DFIN' DXMAX) 'ET'
(A34 'DROIT' A4 'DINI' DXMAX 'DFIN' DXMINC) ;
A4A1 = (A4 'DROIT' A41 'DINI' DXMINH 'DFIN' DXMAX) 'ET'
(A41 'DROIT' A1 'DINI' DXMAX 'DFIN' DXMINH) ;
* DOM1 = 'SURFACE' (A1A2 'ET' A2A3 'ET' A3A4 'ET' A4A1) 'PLAN' ;
DOM1 = 'DALL' A1A2 A2A3 A3A4 A4A1 'PLAN' ;
'SI' GRAPH ;
'TRACER' DOM1 'TITRE' 'Mesh' ;
'FINSI' ;
*
* DOM1, A1A2, A2A3, A3A4, A4A1 are made of 'LINE' elts
*
QDOM1 = 'CHANGER' DOM1 'QUAF' ;
QA4A1 = 'CHANGER' A4A1 'QUAF' ;
QA1A2 = 'CHANGER' A1A2 'QUAF' ;
QA2A3 = 'CHANGER' A2A3 'QUAF' ;
QA3A4 = 'CHANGER' A3A4 'QUAF' ;
'ELIMINATION' (QA4A1 'ET' QDOM1) (DXMINC '/' 10) ;
'ELIMINATION' (QA1A2 'ET' QDOM1) (DXMINC '/' 10) ;
'ELIMINATION' (QA2A3 'ET' QDOM1) (DXMINC '/' 10) ;
'ELIMINATION' (QA3A4 'ET' QDOM1) (DXMINC '/' 10) ;
*
* QDOM1, QA1A2, QA2A3, QA3A4, QA4A1 are made of 'QUAF' elts
*
'SI' VARREST ;
*
* Compatibility between the mesh of the IC and this mesh
*
NN1 = 'NBNO' QDOM1 ;
AA = 'EXTRAIRE' UN 'MAILLAGE' ;
BB = QDOM1 'ET' AA ;
NN2 = 'NBNO' BB ;
'ELIMINATION' AA (DXMINC '/' 10.) QDOM1 ;
NN3 = 'NBNO' (QDOM1 'ET' AA) ;
'SI' ((NN3 'NEG' NN1) 'OU' (NN2 'NEG' (2 '*' NN1))) ;
'ERREUR' 21 ;
'FINSI' ;
AA = 'EXTRAIRE' PN 'MAILLAGE' ;
'ELIMINATION' AA (DXMINC '/' 10.) DOM1 ;
'FINSI' ;
**************
*** DT *******
**************
*
*
DT = CFLREF '*' DX '/' UREF ;
*
***************************************
********** IC / BC ********************
***************************************
*
* IC for thermodyn. pressure
*
'SI' ('NON' VARREST) ;
PMASS = P0 ;
PCOMP = P0 ;
PTHER = P0 ;
'FINSI' ;
*
* IC/BC for the dynamic pressure (alaways 'LINE')
*
* Projection: imposition of PNLIM on a point
* imposition of the speed on the contour
*
* Fully-implicit: imposition of the speed on the contour but in one
* point
*
AA = A1A2 'POIN' 'PROC' ((L '/' 2) 0.0) ;
MPLIM = 'MANUEL' 'POI1' AA ;
PNLIM = 'MANUEL' 'CHPO' MPLIM 1 NOMPRE 0.0 ;
* PNLIM not taken into account in the fully implicit case
'SI' ('NON' VARREST) ;
PN = 'MANUEL' 'CHPO' DOM1 1 NOMPRE 0.0 ;
'FINSI' ;
*
* BC/IC for the speed (always on 'QUAF')
*
* BC on the normal speed
*
'SI' LOGFI ;
LIMDIRN = 'DIFF' ('CHANGER' ('CONTOUR' QDOM1) 'POI1') MPLIM ;
'SINON' ;
LIMDIRN = 'CHANGER' ('CONTOUR' QDOM1) 'POI1' ;
'FINSI' ;
LIMDIRT = ('CHANGER' ('CONTOUR' QDOM1) 'POI1') ;
'SI' GRAPH ;
'TRACER' (QDOM1 'ET' (MPLIM 'COULEUR' ROUG))
'TITRE' 'Point for PLIM' ;
'TRACER' (QDOM1 'ET' (LIMDIRT 'COULEUR' ROUG))
'TITRE' 'Points in which ux is 0' ;
'TRACER' (QDOM1 'ET' (LIMDIRN 'COULEUR' ROUG))
'TITRE' 'Points in which uy is 0' ;
'FINSI' ;
UN1 = 'MANUEL' 'CHPO' LIMDIRT 1 'UX' 0.0 ;
UN2 = 'MANUEL' 'CHPO' LIMDIRN 1 'UY' 0.0 ;
UNLIM = UN1 '+' UN2 ;
'SI' VARREST ;
UNCAN = ('REDU' UN ('CONTOUR' QDOM1)) '*' -1 ;
UN = UN '+' UNLIM '+' UNCAN ;
'SINON' ;
UN = ('MANUEL' 'CHPO' QDOM1 2 'UX' 0.0 'UY' 0.0) '+' UNLIM ;
'FINSI' ;
*
* BC/IC for the temperature ('QUAF'/'LINE')
*
'SI' ('EGA' DISCT 'LINE') ;
TN1 = 'MANUEL' 'CHPO' A4A1 1 'T' Tc ;
TN2 = 'MANUEL' 'CHPO' A2A3 1 'T' Th ;
TNLIM = TN1 '+' TN2 ;
'SI' ('NON' VARREST) ;
TN = ('MANUEL' 'CHPO' DOM1 1 'T' T0) ;
'FINSI' ;
TN = TN '+' TNLIM '-' ('REDU' TN (A4A1 'ET' A2A3)) ;
DOMT = DOM1 ;
'FINS' ;
*
'SI' ('EGA' DISCT 'QUAF') ;
TN1 = 'MANUEL' 'CHPO' QA4A1 1 'T' Tc ;
TN2 = 'MANUEL' 'CHPO' QA2A3 1 'T' Th ;
TNLIM = TN1 '+' TN2 ;
'SI' ('NON' VARREST) ;
TN = ('MANUEL' 'CHPO' QDOM1 1 'T' T0) ;
'FINSI' ;
TN = TN '+' TNLIM '-' ('REDU' TN (QA4A1 'ET' QA2A3)) ;
DOMT = QDOM1 ;
'FINS' ;
*
*************************************************
*** Dynamic viscosity and thermal diffusivity ***
*************************************************
*
* CHMU and CHLAM
*
'SI' VARPRO ;
TNSCA = 'NOMC' TN 'SCAL' ;
CHTSTA = (TNSCA '*' 0.0) '+' tstar ;
CHESSE = (TNSCA '*' 0.0) '+' esse ;
CHTSPE = CHTSTA '+' CHESSE ;
CHMU = mustar '*' (((TNSCA '/' TSTAR) '**' 1.5)
'*' (CHTSPE '/' (TNSCA '+' CHESSE))) ;
CHLAM = CHMU '*' (cp '/' PR) ;
'SINON' ;
'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' ;
'FINSI' ;
*
* Plot of IC
*
'SI' GRAPH ;
VECN = 'VECTEUR' 0.1 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
*
*
* A4 A3
*
*
*
*
* A1 A2
*
* \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' ;
*
NXA4A1 = 'MANUEL' 'CHPO' DOM1 1 'SCAL' -1.0 ;
NYA4A1 = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 0.0 ;
COEFBID = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
MSMD41 = MATFLU QA4A1 'QUAF' NOMPRE 'LINE' 'LINE' COEFBID
'LINE' NXA4A1 NYA4A1 ;
*
*
NXA2A3 = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 1.0 ;
NYA2A3 = 'MANUEL' 'CHPO' DOM1 1 'SCAL' 0.0 ;
MSMD23 = MATFLU QA2A3 'QUAF' NOMPRE 'LINE' 'LINE' COEFBID
'LINE' NXA2A3 NYA2A3 ;
*
*
* - (\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' ;
*
*
*
MSMD1 = (MATLAP QDOM1 'T' DISCT NOMPRE 'LINE' DISCLA CHLAM
'LINE' COEFBID) '*' -1. ;
MSMD = MSMD1 'ET' MSMD23 'ET' MSMD41 ;
*
* Matrix to compute the contribution of
*
* (-1/(gamma * PTHER)) '*' \der{PTHER}{t}
*
* Matrix to compute the contribution of \div(lambda \grad(T)) to
*
MSMDIV = MATMAS QDOM1 'SCAL' NOMPRE 'LINE' 'LINE' ;
*************************************************************************
* Temporal loop *********************************************************
*************************************************************************
*
TAB1 = 'TABLE' ;
TAB1.'TITRE'= TABLE ;
TAB1.1='TIRR ';
TAB1.'TITRE' . 1 = MOT 'Mass cons.';
TAB1.'TITRE' . 2 = MOT 'Div. cons.';
TAB1.2='MARQ CROI';
'SI' LOGCOMP ;
TITP = 'CHAINE' 'Pressure from div cons.' ;
'SINON' ;
TITP = 'CHAINE' 'Pressure from mass cons.' ;
'FINSI' ;
*
TPS = 0.0 ;
*
LIT = 'PROG' ;
LIT1 = 'PROG' ;
LER = 'PROG' ;
LERT = 'PROG' ;
LERP = 'PROG' ;
LPCO = 'PROG' ;
LPMA = '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' ;
*
'REPETER' BL1 NITMAX ;
LIT1 = LIT1 'ET' ('PROG' &BL1) ;
PNM = 'COPIER' PN ;
UNM = 'COPIER' UN ;
TNM = 'COPIER' TN ;
PMASSM = PMASS ;
PCOMPM = PCOMP ;
*
*** Stabilization attempt for the convective term
*
UNTX = ('EXCO' 'UX' UNM) ;
UNTY = ('EXCO' 'UY' UNM) ;
*
**** Time step
*
TPS = TPS '+' DT ;
*****************************************
***** Computation of the transport ******
***** properties ************************
*****************************************
*
'SI' VARPRO ;
TNSCA = 'NOMC' TN 'SCAL' ;
CHMU = mustar '*' (((TNSCA '/' TSTAR) '**' 1.5)
'*' (CHTSPE '/' (TNSCA '+' CHESSE))) ;
CHLAM = CHMU '*' (cp '/' PR) ;
'FINSI' ;
*
*****************************************
**** Divergence constraint and PCOMP ***
*****************************************
*
* We compute qdiv = (div(u),Npj), Npj being the shape 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' ;
* Volume 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)) ;
*
* Surfaces integral contribution
*
*
* We add all the terms.
* Since (div u, 1) = 0
*
coef1 = ('SOMT' qdiv) ;
coef1 = coef1 '/' (l '*' l) ;
chcoef = 'MANUEL' 'CHPO' DOM1 1 'SCAL' (coef1 '*' -1) ;
qdiv = qdiv '+' (MSMDIV '*' chcoef) ;
*
* coef1 thus computed is equal to
*
* coef1 = \frac{\gamma-1}{\gamma PCOMP V} \times
* \sum_j (\div (\lambda \grad T), Npj)
*
* Let us now compute DPCOMSDT equal to
*
* DPCOMSDT = \frac{\gamma-1}{V} \sum_j (\div (\lambda \grad T) , Npj)
* = \gamma PCOMPM coef1
*
*
DPCOMSDT = coef1 '*' (gamma * PCOMPM) ;
PCOMP = PCOMPM '+' (DPCOMSDT '*' DT) ;
LPCO = LPCO 'ET' ('PROG' PCOMP) ;
*
* NB in Boussinesq, DPCOMSDT = 0.0, PCOMP = PCOMPM = P0
*
*********************************************
*** PMASS, via the integration of the EOS ***
*********************************************
*
* Average temperature on each cell
* discpri discdu
INTUST = INTINVT QDOM1 DISCT TN ;
PMASS = (R '*' l '*' l '*' rho0) '/' INTUST ;
DPMASSDT = (PMASS '-' PMASSM) '/' DT ;
'SI' LBOUS;
PMASS = PMASSM ;
DPMASSDT = 0.0 ;
'FINSI' ;
LPMA = LPMA 'ET' ('PROG' PMASS) ;
*
*************
*** PTHER ***
*************
*
'SI' LOGCOMP ;
PTHER = PCOMP ;
DPTHSDT = DPCOMSDT ;
'SINON' ;
PTHER = PMASS ;
DPTHSDT = DPMASSDT ;
'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) '*' TNM) '/' (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' TNM '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
*
COEFL = (((gamma '-' 1) '*' DPTHSDT) '/' (gamma '*' PTHER))
'*' TNM ;
qtt = MSOUTEM '*' COEFL ;
*
* ((1 '/' dt) T^{n}, Nt), RHS
*
qtt = qtt '+' (MMASTSDT '*' TNM) ;
*
* Temperature: LHS, RHS:
*
'MESSAGE' ;
'MESSAGE' 'Computation of T' ;
mtott = mmastsdt 'ET' mkont 'ET' mlapt ;
smbt = qtt ;
mtotik = 'KOPS' 'RIMA' mtott ;
'SI' ((&BL1 'EGA' 1) 'OU' ('EGA' METINV 'DIRECT')) ;
mprect = mtotik ;
'SINON' ;
'SI' (('EGA' ((&BL1 '/' NFREQP) '*' NFREQP) &BL1) '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
*
qtv = (MMASVSDT '*' UNM) ;
*
* ((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 = MMASVSDT 'ET' MKONV 'ET' MLAPV 'ET' MPRE 'ET' MDIV ;
* mtotv = MMASVSDT 'ET' MKONV 'ET' MLAPV ;
* mdiv2 = 'KOPS' 'RIMA' mdiv ;
* mdiv2 = 'KOPS' 'RELA' mdiv2 'SIMP' ;
smbv = qtv '+' qdiv ;
mtotik = 'KOPS' 'RIMA' mtotv ;
* mtotik = 'ET' mtotik mdiv2 ;
'SI' ((&BL1 'EGA' 1) 'OU' ('EGA' METINV 'DIRECT')) ;
mprecufi = mtotik ;
'SINON' ;
'SI' (('EGA' ((&BL1 '/' NFREQP) '*' NFREQP) &BL1) '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 = mmasvsdt 'ET' MKONV 'ET' MLAPV ;
smbv = qtv ;
mtotik = 'KOPS' 'RIMA' mtotv ;
'SI' ((&BL1 'EGA' 1) 'OU' ('EGA' METINV 'DIRECT')) ;
mprecut = mtotik ;
'SINON' ;
'SI' (('EGA' ((&BL1 '/' NFREQP) '*' NFREQP) &BL1) '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 ;
mtotik = 'KOPS' 'RIMA' MLAPP ;
'SI' ((&BL1 'EGA' 1) 'OU' ('EGA' METINV 'DIRECT')) ;
mprecp = mtotik ;
'SINON' ;
'SI' (('EGA' ((&BL1 '/' NFREQP) '*' NFREQP) &BL1) '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 UNLIM 1.0D-16 200
PREC TABRESDU METINV ;
UN = UNTIL '-' DELTAU ;
'FINSI' ;
*
***********************************************
**** Post-treatment and convergence check *****
***********************************************
*
'SI' (((&BL1 '/' NFREQ) '*' NFREQ) 'EGA' &BL1) ;
ERRO = 'MAXIMUM' (UN '-' UNM) 'ABS' ;
ERROT = 'MAXIMUM' (TN '-' TNM) 'ABS' ;
ERROP = 'MAXIMUM' (PN '-' PNM) 'ABS' ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'TEMPORAL ITER = ' &BL1
', TPS = ' TPS ', ERRO = '
ERRO) ;
'MESSAGE' ;
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 ;
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 ;
'SI' (((&BL1 '/' (2 '*' NFREQ)) '*'
(2 '*' NFREQ)) 'EGA' &BL1) ;
'DESSIN' EVERT 'TITR' 'Convergence history (T)' 'YBOR'
-7. 1 'NCLK' ;
* 'DESSIN' EVER 'TITR' 'Convergence history (v)' 'NCLK' ;
'SINON' ;
VUN = 'VECTEUR' UN ;
'TRACER' TN DOM1 VUN 'TITRE'
('CHAINE' 'TPS =' TPS) 'NCLK' ;
* EVP = 'EVOL' 'MANU' 'it' LIT1 'Pther' LPMA ;
* EVP2 = 'EVOL' 'MANU' 'it' LIT1 'Pther' LPCO ;
* 'DESSIN' (EVP 'ET' EVP2) 'TITR' TITP 'LEGE' TAB1 'NCLK' ;
'FINSI' ;
'FINSI' ;
'SI' ((ERROT < 5.0D-2) 'ET' (ERRO < 1.0D-4)) ;
'QUITTER' BL1 ;
'FINSI' ;
'FINSI' ;
*
'FIN' BL1 ;
*
***** Post treatment
*
EVER = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LER ;
'SI' GRAPH ;
'DESSIN' EVER 'TITR' 'Convergence history (v)' ;
'DESSIN' EVERT 'TITR' 'Convergence history (T)' ;
'FINSI' ;
EVERP = 'EVOL' 'MANU' 'it' LIT 'Log(Err)' LERP ;
'SI' GRAPH ;
'DESSIN' EVERP 'TITR' 'Convergence history (P)' ;
'FINSI' ;
*
EVP = 'EVOL' 'MANU' 'it' LIT1 'Pther' LPMA ;
EVP2 = 'EVOL' 'MANU' 'it' LIT1 'Pther (comp. constr.)' LPCO ;
'SI' GRAPH ;
'DESSIN' (EVP 'ET' EVP2) 'TITR' 'Pressure evolution' 'LEGE' TAB1 ;
'FINSI' ;
UX = 'EXCO' UN 'UX' ;
UY = 'EXCO' UN 'UY' ;
VECN = 'VECTEUR' UN ;
'SI' GRAPH ;
'TRACER' TN DOM1 'TITRE' 'T' ;
'TRACER' UX DOM1 VECN 'TITR' 'ux' ;
'TRACER' UY DOM1 VECN 'TITR' 'uy' ;
'FINSI' ;
*
**** Thermal flux and nusselt
*
lamgt = CHPFLU QDOM1 DISCT TN 'QUAF' DISCLA CHLAM ;
NORMX = 'MANUEL' 'CHPO' QDOM1 1 'UX' 1.0 ;
NORMY = 'MANUEL' 'CHPO' QDOM1 1 'UY' 0.0 ;
QQ = 'PSCAL' lamgt (normx '+' normy) ('MOTS' 'UX' 'UY')
('MOTS' 'UX' 'UY') ;
QQ41 = 'REDU' QQ QA4A1 ;
QQ23 = 'REDU' QQ QA2A3 ;
EVQLEFT = 'EVOL' 'CHPO' QQ41 QA4A1 ;
EVQRIGH = 'EVOL' 'CHPO' (QQ23) QA2A3 ;
*
coef = lambda0 '*' (Th '-' Tc) '/' l ;
NULEFT = EVQLEFT '/' coef ;
NURIGH = EVQRIGH '/' coef ;
*
TAB1 = 'TABLE' ;
TAB1 . 1 = 'MOT' 'TIRR ';
TAB1 . 2 = 'MOT' 'MARQ CROI NOLI';
TAB1 . 'TITRE' = 'TABLE' ;
TAB1 . 'TITRE' . 1 = 'MOT' 'Left' ;
TAB1 . 'TITRE' . 2 = 'MOT' 'Right' ;
'SI' GRAPH ;
'DESSIN' (NULEFT 'ET' NURIGH) 'TITRE' 'Nusselt' 'GRILL'
LEGE TAB1 ;
'FINSI' ;
*
aa = 'PRIM' NULEFT ;
aa = 'IPOL' (0.99999 '*' l) ('EXTR' aa 'ABSC') ('EXTR' aa 'ORDO') ;
aa = aa '/' l ;
'SI' GRAPH ;
'DESSIN' NULEFT 'TITRE' ('CHAINE' 'Average Nusselt = ' aa) 'GRILL'
'MIMA' ;
'FINSI' ;
*
bb = 'PRIM' NURIGH ;
bb = 'IPOL' (0.99999 '*' l) ('EXTR' bb 'ABSC') ('EXTR' bb 'ORDO') ;
bb = bb '/' l ;
'SI' GRAPH ;
'DESSIN' NURIGH 'TITRE' ('CHAINE' 'Average Nusselt = ' bb)
'MIMA' ;
'FINSI' ;
*
**** Temperature
*
'SI' GRAPH ;
'OPTION' 'ISOV' 'LIGN' ;
'TRACER' DOMT TN ('CONTOUR' DOMT) 'TITR' 'Temperature' ;
'FINSI' ;
**
** 'OUBLIER' RIGIDITE object
**
* 'OUBLIER' MFLU ;
* 'OUBLIER' MTOTV ;
* 'OUBLIER' MPRE ;
* 'OUBLIER' MLAPV ;
* 'OUBLIER' MKONV ;
* 'OUBLIER' MTOTT ;
* 'OUBLIER' MLAPT ;
* 'OUBLIER' MKONT ;
* 'OUBLIER' MSMDIV ;
* 'OUBLIER' MSMD ;
* 'OUBLIER' MSMD23 ;
* 'OUBLIER' MSMD41 ;
* 'OUBLIER' MSMD1 ;
* 'OUBLIER' MSOUTEM ;
* 'OUBLIER' MSOUGRA ;
* 'OUBLIER' MDIV ;
* 'OUBLIER' MMASTSDT;
* 'OUBLIER' MMAST ;
* 'OUBLIER' MMASVSDT;
* 'OUBLIER' MMASV ;
**
** 'OUBLIER' MATRIK object
**
* 'OUBLIER' MTOTIK ;
* 'OUBLIER' MMASIK ;
* 'OUBLIER' MPRECT ;
* 'OUBLIER' MPRECUT ;
* 'OUBLIER' MPRECP ;
* 'OUBLIER' MPRECUFI ;
*
* 'MENAGE' ;
*
* 'OPTION' 'SAUVER' ('CHAINE' REPSAU
* 'resu' NXS2 '.sauv') ;
* 'SAUVER' ;
*
*
* Tests
*
'SI' (&BL1 > 200) ;
'MESSAGE' 'Convergence not reached.' ;
'ERREUR' 5 ;
'FINSI' ;
'SI' LBOUS ;
aa = extr NULEFT 'ORDO' ;
bb = extr NURIGH 'ORDO' ;
erro = 'MAXIMUM' (aa - bb) 'ABS' ;
'SI' (erro > 1.0D-8) ;
'MESSAGE' ;
'MESSAGE' 'Problem in Boussinesq approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
maxaa = 'MAXIMUM' AA ;
erro = ('ABS' (maxaa '-' 17.9)) '/' 17.9 ;
'SI' (erro > 1.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in Boussinesq approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
'SINON' ;
erro = ('ABS' (PTHER '-' 86734.)) '/' 86734. ;
'SI' (erro > 1.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
* Values on the cold wall (left, A4A1)
* Values on the hot wall (right, A2A3)
aa = extr NULEFT 'ORDO' ;
bb = extr NURIGH 'ORDO' ;
maxaa = 'MAXIMUM' AA ;
minaa = 'MINIMUM' AA ;
maxbb = 'MAXIMUM' BB ;
minbb = 'MINIMUM' BB ;
erro = 'ABS' ((maxaa '-' 16.36)) '/' 16.36 ;
'SI' (erro > 5.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
erro = 'ABS' ((maxbb '-' 19.60)) '/' 19.60 ;
'SI' (erro > 5.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
erro = 'ABS' ((minaa '-' 0.854)) '/' 16.36 ;
'SI' (erro > 5.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
erro = 'ABS' ((minbb '-' 1.073)) '/' 19.60 ;
'SI' (erro > 5.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
aa = 'PRIM' NULEFT ;
aa = 'IPOL' (0.99999 '*' l) ('EXTR' aa 'ABSC') ('EXTR' aa 'ORDO') ;
aa = aa '/' l ;
erro = ('ABS' (aa '-' 8.86)) '/' 8.86 ;
'SI' (erro > 5.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
aa = 'PRIM' NURIGH ;
aa = 'IPOL' (0.99999 '*' l) ('EXTR' aa 'ABSC') ('EXTR' aa 'ORDO') ;
aa = aa '/' l ;
erro = ('ABS' (aa '-' 8.86)) '/' 8.86 ;
'SI' (erro > 5.0D-2) ;
'MESSAGE' ;
'MESSAGE' 'Problem in LM approx.' ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
*
'FIN' BLPAR ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer