* fichier : sinebum_fmm4.dgibi
*************************************************************************
* *
* SINE-SHAPED BUMP *
* *
* CALCUL DE L'ECOULEMENT SUBSONIQUE ISENTROPIQUE STATIONNAIRE DANS *
* UN CANAL *
* *
* Methode implicite sans matrice pour les equations d'Euler *
* (bas Mach) *
* *
* BECCANTINI A., SFME/LTMF, DEC 2003 *
* *
* *
*************************************************************************
*************************************************************************
*************************************************************************
*********** SOLUTION OF THE EULER EQUATIONS ****************************
*********** FMM ****************************
*************************************************************************
*************************************************************************
*
*
* RVX . 'MODEL' : model object
* RVX . 'RN0' : density
* RVX . 'GN0' : qdm
* RVX . 'RET0' : total energy
*
* RVX . 'PGAS' : gas model
* RVX . 'PGAS' . 'GAMN' : gamma
*
* RVX . 'GRAVITY' : gravity
*
* RVX . 'LISTCONS' : name of the conservative variables
* RVX . 'LISTPRIM' : name of the primitive variables
*
* RVX . 'LISTERR' : name of the error variables
*
* RVX . 'METHOD' : numerical scheme
* RVX . 'CUTOFF' : cut off speed
*
* RVX . 'SPACEA' : space accuracy
* RVX . 'LIMITER' : limiter type
* RVX . 'TIMEA' : time accuracy
*
* RVX . 'T0' : initial time
* RVX . 'TFINAL' : final time
* RVX . 'DTPS' ('CFL') : time step or CFL number
*
* RVX . 'DCFL' : CFL number for dual time
*
* Error criteria for dual time loop:
* RVX . 'NDTITER' : number of iterations
* RVX . 'EPSDT' : relative error
*
* RVX . 'NJAC' : Jacobi iterations
*
'DEBPROC' PEUSM ;
'ARGUMENT' RVX*TABLE ;
*
'SI' ('EXISTE' RVX 'RESULTS') ;
'MESS' 'Table RESULTS already exists' ;
'SINON' ;
RVX . 'RESULTS' = 'TABLE' ;
RVX . 'RESULTS' . 'TPS' = RVX . 'T0' ;
RVX . 'RESULTS' . 'RN' = 'COPIER' (RVX . 'RN0') ;
RVX . 'RESULTS' . 'GN' = 'COPIER' (RVX . 'GN0') ;
RVX . 'RESULTS' . 'RET' = 'COPIER' (RVX . 'RET0') ;
RVX . 'RESULTS' . 'LISTLINF' = 'PROG' ;
RVX . 'RESULTS' . 'LISTITDT' = 'LECT' ;
RVX . 'RESULTS' . 'LISTITER' = 'LECT' ;
RVX . 'RESULTS' . 'NITER' = 0 ;
'FINSI' ;
MDOMINT = RVX . 'MODEL' ;
NELT = 'NBEL' ('DOMA' MDOMINT 'CENTRE') ;
GAMN = RVX . 'PGAS' . 'GAMN' ;
*
LISTINCO = RVX . 'LISTCONS' ;
LISTPRIM = RVX . 'LISTPRIM' ;
LISTERR = RVX . 'LISTERR' ;
*
* Upwind scheme
*
METO = RVX . 'METHOD' ;
* Space accuracy (1 or 2) and limiter
* Time accuracy (1 or 2)
ORDESP = RVX . 'SPACEA' ;
TYPELIM = RVX . 'LIMITER' ;
ORDTPS = RVX . 'TIMEA' ;
* Initial/final time
* Deltat or CFL
TPS = RVX . 'RESULTS' . 'TPS' ;
TFINAL = RVX . 'TFINAL' ;
'SI' ('EXISTE' RVX 'DTPS') ;
'SI' ('EXISTE' RVX 'CFL') ;
'MESSAGE' 'DTPS or CFL ???' ;
'ERREUR' 21 ;
'FINSI' ;
DTPS = RVX . 'DTPS' ;
'SINON' ;
CFL = RVX . 'CFL' ;
'FINSI' ;
* Dual time iterations
NDT = RVX . 'NDTITER' ;
* Relative error
EPSDT = RVX . 'EPSDT' ;
* Jacobi iterations
* NJAC = RVX . 'NJAC' ;
* Cut off speed
ICO = RVX . 'CUTOFF' ;
*
**** Conservative variables
*
*
MOT1 = 'EXTRAIRE' LISTINCO 1 ;
'SI' ('EGA' ('VALE' 'DIME') 2) ;
NOMMOM = 'EXTRAIRE' LISTINCO
('LECT' 2 3 ) ;
NOMVEL = 'MOTS' 'UX' 'UY' ;
MOT2 = 'EXTRAIRE' LISTINCO 4 ;
'SINON' ;
NOMMOM = 'EXTRAIRE' LISTINCO
('LECT' 2 3 4) ;
NOMVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
MOT2 = 'EXTRAIRE' LISTINCO 5 ;
'FINSI' ;
RN0 = 'REDU' (RVX . 'RESULTS' . 'RN') ('DOMA' MDOMINT 'CENTRE') ;
GN0 = 'REDU' (RVX . 'RESULTS' . 'GN') ('DOMA' MDOMINT 'CENTRE') ;
RET0 = 'REDU' (RVX . 'RESULTS' . 'RET') ('DOMA' MDOMINT 'CENTRE') ;
*
**** Primitive variables
*
VN0 PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
MOTRN = 'EXTRAIRE' LISTPRIM 1 ;
'SI' ('EGA' ('VALE' 'DIME') 2) ;
MOTVN = 'EXTRAIRE' LISTPRIM
('LECT' 2 3 ) ;
MOTPN = 'EXTRAIRE' LISTPRIM 4 ;
'SINON' ;
MOTVN = 'EXTRAIRE' LISTPRIM
('LECT' 2 3 4) ;
MOTPN = 'EXTRAIRE' LISTPRIM 5 ;
'FINSI' ;
*
**** Coeff to compute gradients (MCHSCA, MCHVEC)
* Boundary conditions have to be taken into account
*
RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
RN0 VN0 PN0 GAMN ;
MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
'SI' (MAILLIM 'EGA' 0) ;
GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN0 ;
GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
NOMVEL GN0 ;
'SINON' ;
SCALBC = 'MANUEL' 'CHPO' MAILLIM 1 'SCAL' 1.0 ;
'SI' (('VALEUR' 'DIME') 'EGA' 2) ;
VECTBC = 'MANUEL' 'CHPO' MAILLIM
2 'UX' 0.0 'UY' 0.0 ;
* Same name as NOMVEL
'SINON' ;
VECTBC = 'MANUEL' 'CHPO' MAILLIM
3 'UX' 0.0 'UY' 0.0 'UZ' ;
'FINSI' ;
GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN0
'CLIM' SCALBC ;
GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
NOMVEL VN0
'CLIM' VECTBC ;
'FINSI' ;
*
*
*** After each dual time loop, we could display
* the evolution of the error in the dual time loop
* the evolution of the dual time step (the safety factor)
*
LISTLINF = RVX . 'RESULTS' . 'LISTLINF' ;
LISTITDT = RVX . 'RESULTS' . 'LISTITDT' ;
LISTITER = RVX . 'RESULTS' . 'LISTITER' ;
*
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ;
'TEMPS' ;
************************************************************************
************************************************************************
**** Temporal loop *****************************************************
************************************************************************
************************************************************************
RN_N1M1 = 'COPIER' RN0 ;
GN_N1M1 = 'COPIER' GN0 ;
RET_N1M1 = 'COPIER' RET0 ;
AA = 'DIME' LISTITER ;
'SI' (AA > 0) ;
PTITER = 'EXTRAIRE' LISTITER AA ;
'SINON' ;
PTITER = 0 ;
'FINSI' ;
DUSDT = 0.0D0 ;
'REPETER' BLITER -1 ;
PTITER = PTITER '+' 1 ;
*
**** Personal procedure
*
PROCPT RVX ;
*
*
**** _N1M = (t^n,\tau^m)
* _N1M1 = (t^n,\tau^{m+1})
*
*
************************************************************************
****** Loop on dual time***********************************************
************************************************************************
*
*
*** DUSDT0 is the increment of DUSDT in the previous (physical) time
* iteration.
*
DUSDT0 = DUSDT ;
DUSDT = 0.0D0 ;
'REPETER' BLDT NDT ;
RN_N1M = RN_N1M1 ;
GN_N1M = GN_N1M1 ;
RET_N1M = RET_N1M1 ;
*
**** Primitive variables
*
VN_N1M PN_N1M = 'PRIM' 'PERFMONO' RN_N1M GN_N1M RET_N1M GAMN ;
*
*** Boundary conditions
*
RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
RN_N1M VN_N1M PN_N1M GAMN ;
*
****** First/second order reconstruction
*
'SI' (ORDESP 'EGA' 2) ;
GRADRN ALRN0 = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') RN_N1M
'CLIM' ('EXCO' MOTRN RCHLIM 'SCAL') 'GRADGEO' MCHSCA ;
GRADPN ALPN0 = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
('MOTS' 'SCAL') PN_N1M
'CLIM' ('EXCO' MOTPN RCHLIM 'SCAL') 'GRADGEO' MCHSCA ;
GRADVN ALVN0 = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' TYPELIM
('MOTS' 'UX' 'UY') VN_N1M
'CLIM' ('EXCO' MOTVN RCHLIM MOTVN) 'GRADGEO' MCHVEC ;
* 'SI' (&BLDT < NLCB) ;
* ALRN0 = 'COPIER' ALRN ;
* ALPN0 = 'COPIER' ALPN ;
* ALVN0 = 'COPIER' ALVN ;
* 'SINON' ;
* 'SI' (&BLDT 'EGA' NLCB) ;
* 'MESSAGE' ;
* 'MESSAGE' 'On gele les limiteurs!!!' ;
* 'MESSAGE' ;
* 'FINSI' ;
* 'FINSI' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
MDOMINT
RN_N1M GRADRN ALRN0
VN_N1M GRADVN ALVN0
PN_N1M GRADPN ALPN0
GAMN ;
'SINON' ;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
MDOMINT
RN_N1M VN_N1M PN_N1M GAMN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
MDOMINT LISTINCO ROF VITF PF GAMF MAILLIM
ICO ICO ;
RESIDU = RESIDU '+' RESLIM ;
*
**** La gravite
*
RESGRA = 'FIMP' 'VF' 'GRAVMONO' 'RESI' LISTINCO
RN_N1M GN_N1M1 (RVX . 'GRAVITY') ;
RESIDU = RESIDU '+' RESGRA ;
*
****** Residuum for dual tims stepping also involved the
* variation of the conserved variables with respect
* to time
'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
RESIDU = RESIDU '-' DUSDT ;
'SINON' ;
RESIDU = RESIDU '-' ((1.5 '*' DUSDT) '-' (0.5 '*' DUSDT0)) ;
'FINSI' ;
*
*** Time step at the first iteration/jacobi iteration
*
'SI' (&BLDT 'EGA' 1) ;
'SI' ('EXISTE' RVX 'CFL') ;
DTPS = (RVX . 'CFL') '*' 2.0D0 '*' DELTAT ;
'SINON' ;
DTPS = RVX . 'DTPS' ;
'FINSI' ;
DTPS = 'MINIMUM' ('PROG' DTPS ((TFINAL '-' TPS) '*' 1.001)) ;
TPS = TPS '+' DTPS ;
*
NJAC = 'ENTIER' ('MINIMUM' (RVX . 'NJACITER')) ;
NJAC0 = NJAC ;
'SINON' ;
NJAC = 'IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
(RVX . 'NJACLERR') (RVX . 'NJACITER') ;
NJAC = 'MAXIMUM' ('PROG' NJAC0 NJAC) ;
NJAC = 'ENTIER' NJAC ;
NJAC0 = NJAC ;
'FINSI' ;
*
*** JACOBI
*
'SI' ('NON' ('EXISTE' RVX 'TYPEJAC')) ;
LOGBF = (&BLDT '/' 2 '*' 2) 'EGA' &BLDT ;
'SI' LOGBF ;
RVX . 'TYPEJAC' = 'LJACOF' ;
'SINON' ;
RVX . 'TYPEJAC' = 'LJACOB' ;
'FINSI' ;
'FINSI' ;
*
**** CFL dual
*
'SI' (&BLDT 'EGA' 1) ;
SAFFACD = ('MINIMUM' (RVX . 'DCFL')) '*' 2 ;
'SINON' ;
SAFFACD = ('IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
(RVX . 'DCFLERR') (RVX . 'DCFL')) '*' 2 ;
'FINSI' ;
'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
* DUN IPRO = 'DETO' (RVX . 'TYPEJAC')
DUN IPRO = 'KONV' 'VF' 'PMON1FMM' (RVX . 'TYPEJAC')
LISTINCO MDOMINT RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
DTPS SAFFACD
NJAC 'CLIM' LISTPRIM RCHLIM (RN_N1M '*' 0) ;
'SINON' ;
* DUN IPRO = 'DETO' (RVX . 'TYPEJAC')
DUN IPRO = 'KONV' 'VF' 'PMON1FMM' (RVX . 'TYPEJAC')
LISTINCO MDOMINT RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
(DTPS '/' 1.5)
SAFFACD NJAC 'CLIM' LISTPRIM RCHLIM (RN_N1M '*' 0) ;
'FINSI' ;
* 'FINSI' ;
'SI' (IPRO 'NEG' 0) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dans FMM' ;
'MESSAGE' ;
'ERREUR' 21 ;
'FINSI' ;
*
**** We compute DUSDT for the future loop
*
DUSDT = DUSDT '+' (DUN '/' DTPS) ;
*
**** We evaluate the conservative variables at t^{n+1}, \tau^{m+1}
*
DRN = 'EXCO' MOT1 DUN 'SCAL' ;
DGN = 'EXCO' NOMMOM DUN NOMVEL ;
DRET = 'EXCO' MOT2 DUN 'SCAL' ;
RN_N1M1 = RN_N1M '+' DRN ;
GN_N1M1 = GN_N1M '+' DGN ;
RET_N1M1 = RET_N1M '+' DRET ;
ERRINF = 'MAXIMUM' DUN 'ABS' LISTERR ;
LISTLINF = LISTLINF 'ET' ('PROG' ERRINF) ;
LISTITDT = LISTITDT 'ET' ('LECT' &BLDT) ;
LISTITER = LISTITER 'ET' ('LECT' PTITER) ;
'SI' (RVX . 'RELERR') ;
* Relative error
'SI' (&BLDT 'EGA' 1) ;
ERRINF0 = ERRINF ;
'SINON' ;
'SI' (ERRINF < (EPSDT '*' ERRINF0)) ;
'QUITTER' BLDT ;
'FINSI' ;
'FINSI' ;
'SINON' ;
'SI' (ERRINF < EPSDT) ;
'QUITTER' BLDT ;
'FINSI' ;
'FINSI' ;
'SI' ((&BLDT 'EGA' 1) 'OU'
(((&BLDT '/' 20) '*' 20) 'EGA' &BLDT)) ;
'MESSAGE' ;
'MESSAGE'
('CHAINE' 'ITER =' PTITER ' TPS =' TPS
' DTITER =' &BLDT ' LINF =' ERRINF
' DCFL =' SAFFACD ' NJAC =' NJAC) ;
'MESSAGE' ;
'FINSI' ;
'FIN' BLDT ;
************************************************************************
****** End of the loop on dual time*************************************
************************************************************************
*
*** Update of RVX . 'RESULTS'
*
RVX . 'RESULTS' . 'RN' = RN_N1M1 ;
RVX . 'RESULTS' . 'GN' = GN_N1M1 ;
RVX . 'RESULTS' . 'RET' = RET_N1M1 ;
RVX . 'RESULTS' . 'LISTLINF' = LISTLINF ;
RVX . 'RESULTS' . 'LISTITDT' = LISTITDT ;
RVX . 'RESULTS' . 'LISTITER' = LISTITER ;
RVX . 'RESULTS' . 'TPS' = TPS ;
RVX . 'RESULTS' . 'NITER' = (RVX . 'RESULTS' . 'NITER') '+' 1 ;
*
'SI' (TPS '>EG' TFINAL) ;
'QUITTER' BLITER ;
'FINSI' ;
'FIN' BLITER ;
TCPU = 'TEMPS' 'NOEC' ;
RVX . 'RESULTS' . 'TCPU' = TCPU ;
'FINPROC' ;
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
TYEL = 'TRI3' ;
* TYEL = 'QUA4' ;
'OPTION' 'DIME' 2 'ELEM' TYEL 'ISOV' 'SULI'
'ECHO' 0 'TRAC' 'X' ;
* GRAPH = VRAI ;
GRAPH = FAUX ;
******************
******************
******************
**** MAILLAGE ****
******************
******************
******************
******************
*
*
RAF = 2 ;
NY = 5 '*' RAF ;
NX1 = 4 '*' RAF ;
NX2 = 2 '*' NX1 ;
NX3 = NX1 ;
NX = (NX1 '+' NX2 '+' NX3) ;
DX = (4.0 '/' NX) ;
A0 = -2.0 0.0 ;
A1 = -1.0 0.0 ;
A2 = 1.0 0.0 ;
A3 = 2.0 0.0 ;
A4 = 2.0 1.0 ;
A5 = -2.0 1.0 ;
*
**** LIGB
*
LIGB1 = A0 'DROIT' NX1 A1 ;
* LIGB2
xcel = ('COORDONNEE' 1 A1) '+' DX ;
ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
ACEL = xcel ycel ;
LIGB2 = A1 'DROIT' 1 ACEL ;
'REPETER' BL1 (NX2 '-' 2) ;
ACEL0 = ACEL ;
xcel = xcel '+' DX ;
ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
ACEL = xcel ycel ;
LIGB2 = LIGB2 'ET' (ACEL0 'DROIT' 1 ACEL) ;
'FIN' BL1;
LIGB2 = LIGB2 'ET' (ACEL 'DROIT' 1 A2) ;
LIGB3 = A2 'DROIT' NX3 A3 ;
LIGB = LIGB1 'ET' LIGB2 'ET' LIGB3 ;
*
**** LIGH
*
LIGH = A4 'DROIT' NX A5 ;
*
**** DOMINT
*
DOMINT = LIGH 'REGLER' NY ('INVERSE' LIGB) ;
* DOMINT = LIGB 'REGLER' NY ('INVERSE' LIGH) ;
LIGCON = 'CONTOUR' DOMINT ;
*
*** LIGG
*
LIGG = LIGCON 'ELEM' 'COMP' A5 A0 ;
*
**** LIGD
*
LIGD = LIGCON 'ELEM' 'COMP' A3 A4 ;
'SI' GRAPH ;
'TRACER' DOMINT 'TITRE'
('CHAINE' 'Domaine, nbel=' ('NBEL' DOMINT)) ;
'FINSI' ;
*
**** Creation of MODELS
*
MDOMINT = 'MODELISER' DOMINT 'EULER' ;
MLIGG = 'MODELISER' LIGG 'EULER' ;
MLIGD = 'MODELISER' LIGD 'EULER' ;
TDOMINT = 'DOMA' MDOMINT 'VF' ;
TLIGG = 'DOMA' MLIGG 'VF' ;
TLIGD = 'DOMA' MLIGD 'VF' ;
QDOMINT = TDOMINT . 'QUAF' ;
QLIGD = TLIGD . 'QUAF' ;
QLIGG = TLIGG . 'QUAF' ;
'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QLIGG ;
'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QLIGD ;
***************************************************************
***************************************************************
***************************************************************
************** Initial conditions *****************************
***************************************************************
***************************************************************
***************************************************************
GAMAIR = 1.4 ;
RO = 1.4 ;
P = 1.0 ;
U = 0.01 ;
UCO = 1.0D-3*U ;
***************************************************************
***************************************************************
***************************************************************
************** Boundary Conditions ****************************
***************************************************************
***************************************************************
***************************************************************
*
* Subsonic inlet.
* We impose the total enthalpy and the entropy.
* The velocity direction is supposed to be 90 degrees
*
HTIN = (GAMAIR '/' (GAMAIR '-' 1.0)) '*' P ;
ECIN = 0.5 '*' RO '*' U '*' U ;
HTIN = HTIN '+' ECIN ;
HTIN = HTIN '/' RO ;
SINL = P '/' (RO '**' GAMAIR) ;
CHPINL = 'MANUEL' 'CHPO' ('DOMA' MLIGG 'CENTRE') 2 'HT' HTIN
'S' sinl ;
*
* Subsonic outlet
* We impose the pressure
*
CHPOUT = 'MANUEL' 'CHPO' ('DOMA' MLIGD 'CENTRE') 1 'PN' P ;
***************************************************************
***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
***************************************************************
'DEBPROC' CONS ;
'ARGUMENT' RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
RVN = RN '*' VN ('MOTS' 'SCAL' 'SCAL') ('MOTS' 'UX' 'UY')
('MOTS' 'UX' 'UY') ;
CELL = 'PSCAL' RVN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
RECIN = 0.5 '*' CELL ;
REIN = PN '/' (GAMN '-' 1.0) ;
RETN = RECIN '+' REIN ;
DETR CELL ;
DETR RECIN ;
DETR REIN ;
'RESPRO' RVN RETN ;
'FINPROC' ;
***************************************************************
***************************************************************
***************************************************************
RN0 = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' ro
'NATURE' 'DISCRET' ;
VN0 = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2 'UX' u
'UY' 0.0 'NATURE' 'DISCRET' ;
PN0 = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' p
'NATURE' 'DISCRET' ;
GAMN = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' GAMAIR
'NATURE' 'DISCRET' ;
GN0 RET0 = CONS RN0 VN0 PN0 GAMN ;
VN1 PN1 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
ERRO = 'MAXIMUM' (PN1 '-' PN0) 'ABS' ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' 'Problem in the ic file!!!'
'ERREUR' 5 ;
'FINSI' ;
*
**** Plot of IC
*
MOD1 = 'MODELISER' ('DOMA' MDOMINT 'MAILLAGE' ) 'THERMIQUE';
'SI' GRAPH ;
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN0 ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN0 ;
CHM_VN = 'KCHA' MDOMINT 'CHAM' VN0 ;
'TRACER' CHM_RN MOD1
'TITR' ('CHAINE' 'RN at t=' 0.0);
'TRACER' CHM_PN MOD1
'TITR' ('CHAINE' 'PN at t=' 0.0);
'TRACER' CHM_VN MOD1
'TITR' ('CHAINE' 'MACHN at t=' 0.0);
'FINSI' ;
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
************ COMPUTATION OF THE SOLUTION *****************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
******PROCEDURE PROCPT **************************************************
*************************************************************************
*
*
'DEBPROC' PROCPT ;
'FINPROC' ;
*
*************************************************************************
******PROCEDURE PROLIM **************************************************
*************************************************************************
*
'DEBPROC' PROLIM ;
'ARGUMENT' RVX*'TABLE' MDOMINT*'MMODEL' LISTINCO*'LISTMOTS'
LISTPRIM*'LISTMOTS'
RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
*
*** Boundary conditions
* Subsonic inlet and subsonic outlet
*
RCHLIM1 RESLIM1 = 'KONV' 'VF' 'PERFMONO' 'CLIM' 'RESI'
MDOMINT (RVX . 'MODLI1') LISTINCO LISTPRIM
RN VN PN GAMN (RVX. 'CHPLI1') 'INSU' ;
RCHLIM2 RESLIM2 = 'KONV' 'VF' 'PERFMONO' 'CLIM' 'RESI'
MDOMINT (RVX . 'MODLI2') LISTINCO LISTPRIM
RN VN PN GAMN (RVX. 'CHPLI2') 'OUTP' ;
* RCHLIM = primitive variables on the boundary
* RESLIM = contribution of boundary conditions to the residuum
RESLIM = RESLIM1 '+' RESLIM2 ;
RCHLIM = RCHLIM1 '+' RCHLIM2 ;
'RESPRO' RCHLIM RESLIM ;
'FINPROC' ;
*************************************************************************
****** FIN PROCEDURE PROLIM *********************************************
*************************************************************************
**********************************************************
**********************************************************
**********************************************************
************** COMPUTATION OF THE SOLUTION ***************
**********************************************************
**********************************************************
**********************************************************
**********************************************************
*
RV = 'TABLE' ;
RV . 'MODEL' = MDOMINT ;
*
**** Conservative variables / primitive variables
*
MOTRN = 'RN' ;
MOTGNX = 'RUX' ;
MOTGNY = 'RUY' ;
MOTVNX = 'UX' ;
MOTVNY = 'UY' ;
MOTVN = 'MOTS' 'UX' 'UY' ;
NOMMOM = 'MOTS' MOTGNX MOTGNY ;
MOTRET = 'RETN' ;
MOTPN = 'PN' ;
*
RV . 'LISTCONS' = 'MOTS' MOTRN MOTGNX MOTGNY MOTRET ;
RV . 'LISTPRIM' = 'MOTS' MOTRN MOTVNX MOTVNY MOTPN ;
*
RV . 'RN0' = RN0 ;
RV . 'GN0' = GN0 ;
RV . 'RET0' = RET0 ;
*
**** Gas property/gravity
*
RV . 'PGAS' = 'TABLE' ;
RV . 'PGAS' . 'GAMN' = GAMN ;
RV . 'GRAVITY' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2 'UX' 0.0
'UY' 0.0 ;
*
* Table for BC
*
RV . 'PROLIM' = 'TABLE' ;
RV . 'PROLIM' . 'MODLI1' = MLIGG ;
RV . 'PROLIM' . 'MODLI2' = MLIGD ;
RV . 'PROLIM' . 'CHPLI1' = CHPINL ;
RV . 'PROLIM' . 'CHPLI2' = CHPOUT ;
*
**** Numerical parameters
*
*
* Variable to compute L2 error
*
RV . 'LISTERR' = 'MOTS' MOTGNX MOTGNY ;
* RV . 'LISTERR' = 'MOTS' MOTRN ;
* RV . 'LISTERR' = 'MOTS' MOTRET ;
*
* Upwind scheme
*
* RV . 'METHOD' = 'RUSANOLM' ;
* RV . 'METHOD' = 'HLLCLM' ;
* RV . 'METHOD' = 'ROELM' ;
RV . 'METHOD' = 'AUSMPUP' ;
RV . 'CUTOFF' = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' uco
'NATURE' 'DISCRET' ;
*
* Reconstruction/limiter
* Time accuracy (1 or 2)
* Iterations
* Final time
RV . 'SPACEA' = 2 ;
RV . 'LIMITER' = 'NOLIMITE' ;
* RV . 'LIMITER' = 'LIMITEUR' ;
RV . 'TIMEA' = 1 ;
*
**** Phisical time
*
RV . 'T0' = 0 ;
RV . 'TFINAL' = 1.0D6 ;
* RV . 'CFL' = 20. ;
RV . 'DTPS' = 1.0D16 ;
*
**** Dual time
*
* Safety factor for the dual time step
* Max. Dual time iterations
* Relative error
*
RV . 'DCFLERR' = 'PROG' 16. -32. ;
RV . 'DCFL' = 'PROG' 1.0D2 1.0D2 ;
RV . 'NDTITER' = 2000 ;
RV . 'RELERR' = FAUX ;
RV . 'EPSDT' = 1.0D-13 ;
**** Jacobi iterations
RV . 'TYPEJAC' = 'PJACO' ;
* RV . 'TYPEJAC' = 'LJACOF' ;
* RV . 'TYPEJAC' = 'LJACOB' ;
* RV . 'TYPEJAC' = 'LJACOFB' ;
RV . 'NJACITER' = 'PROG' 15 15 ;
RV . 'NJACLERR' = 'PROG' -1 -16 ;
*
**** Parameters for PROCPT
*
*
'TEMPS' 'ZERO' ;
PEUSM RV ;
'TEMPS' ;
* 'TEMPS' 'IMPR' ;
RN = RV . 'RESULTS' . 'RN' ;
GN = RV . 'RESULTS' . 'GN' ;
RET = RV . 'RESULTS' . 'RET' ;
LISTITER = RV . 'RESULTS' . 'LISTITER' ;
LISTITDT = RV . 'RESULTS' . 'LISTITDT' ;
LISTLINF = RV . 'RESULTS' . 'LISTLINF' ;
**********************************************************
**********************************************************
**********************************************************
************** PLOTS *************************************
**********************************************************
**********************************************************
**********************************************************
**********************************************************
*
*** Convergence evolution inside of each iteration
*
'SI' (VRAI) ;
NITERE = 'DIME' LISTITER ;
I1 = 1 ;
AA = 'PROG' ;
BB = 'PROG' ;
CC = 'PROG' ;
'REPETER' BLITER NITERE ;
I0 = I1 ;
I1 = 'EXTRAIRE' LISTITER &BLITER ;
'SI' (I1 'EGA' I0) ;
AA = AA 'ET' ('PROG' ('EXTRAIRE' LISTITDT &BLITER)) ;
BB = BB 'ET' ('PROG' ('EXTRAIRE' LISTLINF &BLITER)) ;
'SINON' ;
everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
(('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
('LOG' 10.)) ;
'SI' GRAPH ;
'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
(I1 '-' 1)) ;
'FINSI' ;
AA = 'PROG' ;
BB = 'PROG' ;
'FINSI' ;
'FIN' BLITER ;
everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
(('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
('LOG' 10.)) ;
'SI' GRAPH ;
'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
(I1)) ;
'FINSI' ;
'FINSI' ;
*
**** The 2D graphics
*
VN PN = 'PRIM' 'PERFMONO' RN GN RET GAMN 'TRICHE' ;
CN2 = GAMN '*' (PN '/' RN) ;
VN2 = 'PSCAL' VN VN MOTVN MOTVN ;
MACHN2 = VN2 '/' CN2 ;
MACHN = MACHN2 '**' 0.5 ;
HTN = (GAMN '/' (GAMN '-' 1.0)) '*' (PN '/' RN) ;
ECIN = 0.5 '*' ('PSCAL' VN VN MOTVN MOTVN) ;
HTN = HTN '+' ECIN ;
SN = PN '/' (RN '**' GAMAIR) ;
tps = RV . 'RESULTS' . 'TPS' ;
'SI' GRAPH ;
CHM_RN = 'KCHA' MDOMINT 'CHAM' RN ;
CHM_VN = 'KCHA' MDOMINT 'CHAM' VN ;
CHM_PN = 'KCHA' MDOMINT 'CHAM' PN ;
CHM_MN = 'KCHA' MDOMINT 'CHAM' MACHN ;
CHM_HTN = 'KCHA' MDOMINT 'CHAM' HTN ;
CHM_SN = 'KCHA' MDOMINT 'CHAM' SN ;
'TRAC' CHM_RN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'rho at t=' TPS) ;
'TRAC' CHM_VN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'v at t= ' TPS) ;
'TRAC' CHM_PN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'p at t= ' TPS) ;
'TRAC' CHM_MN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'Mach at t= ' TPS) ;
'TRAC' CHM_HTN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'ht at t= ' TPS ' hl =' HTL) ;
'TRAC' CHM_SN MDOMINT ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'sn at t= ' TPS ' hl =' SL) ;
'OPTION' 'ISOV' 'LIGN' ;
RNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' RN) ;
'TRACER' DOMINT RNV ('CONTOUR' DOMINT) 15 'TITRE' 'ro';
PNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' PN) ;
'TRACER' DOMINT PNV ('CONTOUR' DOMINT) 15 'TITRE' 'p' ;
MNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' MACHN) ;
'TRACER' DOMINT MNV ('CONTOUR' DOMINT) 15 'TITRE' 'Mach' ;
'FINSI' ;
* TEST
LERR = 'EXTRAIRE' EVERR 'ORDO' ;
AA = 10. '**' ('EXTRAIRE' LERR 1) ;
BB = 10. '**' ('EXTRAIRE' LERR ('DIME' LERR)) ;
'SI' ((BB '/' AA) '>' 1.0D-5) ;
'ERREUR' 'Probleme de convergence' ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer