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