1 : * EXEC PROCEDUR MAGN 11/09/12 21:16:08 7117 2 : *X EXEC (Procedure) 3 : * Procedure EXEC 4 : * 5 : * Objet : Execute un algorithme décrit dans une table RV 6 : * de type EQEX. 7 : * Cette table est créée par l'opérateur EQEX 8 : * 9 : * Syntaxe : EXEC RV ; 10 : * 11 : *********************************************************************** 12 : * VERSION : ???? 13 : * HISTORIQUE : 20/12/99: gounand 14 : * Rajout de la gestion de la matrice servant à l'assemblage 15 : * (rv . 'METHINV' . 'MATASS') 16 : * 17 : * HISTORIQUE : 12/05/06 : ajout STOPITER, NUITER et STOPPDT 18 : * HISTORIQUE : 21/12/07 : ajout projection algébrique incrémentale 19 : * (cf. GRESP) 20 : * HISTORIQUE : 21 : ************************************************************************ 22 : * DISCPRES = (d'apres KMIC) 23 : * KPRE=3 pression P0 KPRE=4 pression P1 KPRE=2 cas macro 1ère génération 24 : * KPRE=5 pression MSOMMET 25 : ************************************************************************ 26 : 'DEBPROC' EXEC ; 27 : 'ARGUMENT' rv*'TABLE ' nparti/'ENTIER' lopp/'LISTMOTS'; 28 : * 29 : * 30 : * 31 : * Préparation pour paralléliser le comportement 32 : * 33 : Si(EXIST RV 'KPME'); 34 : KPME=rv.'KPME'; 35 : Sinon; 36 : KPME=FAUX; 37 : Finsi; 38 : 39 : si (exist nparti);npart=nparti;sinon;npart=valeur assis;finsi; 40 : Si(NON (KPME));npart=1;Finsi; 41 : 42 : * On partitionne en fonction du nombre de proc ou plus voir moins 43 : * si précisé dans l'instruction exec rv n ; 44 : * Si Npart1); 85 : $mtx= part npart $mti 'NOOPT'; 86 : Sinon ; 87 : $mtx= $mti; 88 : Finsi; 89 : TBPART.$mti=$mtx; 90 : Finsi ; 91 : 92 : Si (Npart > 1); 93 : 94 : repeter bloc1 npart; 95 : tnsi=copie (rv.notable); 96 : tnsi.'DOMZ'=$mti; 97 : tnsi= enlev tnsi 'TDOMZ'; 98 : fin bloc1; 99 : 100 : TBPART.notable=$tns; 101 : *Sinon; 102 : *TBPART.notable=notable; 103 : Finsi; 104 : 105 : Finsi; 106 : fin BPART; 107 : 108 : Finsi; 109 : 110 : 111 : TBPART=rv.'TBPART'; 112 : * 113 : * Fin Préparation pour paralléliser le comportement 114 : * 115 : * 116 : * 117 : 118 : *Logique d'arret 119 : 'SI' ('NON' ('EXISTE' rv 'STOPITER')) ; 120 : rv . 'STOPITER' = FAUX ; 121 : 'FINSI' ; 122 : 'SI' ('NON' ('EXISTE' rv 'STOPPDT')) ; 123 : rv . 'STOPPDT' = FAUX ; 124 : 'FINSI' ; 125 : 126 : * Logique pilotant le recalcul de la matrice de pression 127 : CALPRE=FAUX ; 128 : Si ('EXIST' rv 'CALPRE') ; 129 : vertytab rv 'CALPRE' 'LOGIQUE' ; 130 : CALPRE=rv.'CALPRE' ; 131 : Finsi ; 132 : 133 : Si(NON('EXIST' rv 'XEQUA')); 134 : rv.'XEQUA'= FAUX ; 135 : Finsi ; 136 : 137 : 138 : 'SI' ('EGA' (rv . 'NAVISTOK') 0) ; 139 : EXAC rv ; 140 : 'QUITTER' EXEC ; 141 : 'FINSI' ; 142 : 143 : *******PROJ************************************************************* 144 : *******PROJ Trois possibilités ***************************************** 145 : *******PROJ TMDM1 (vrai)_> Correction Gresho NON 146 : *******PROJ TGRAD (vrai)-> Formulation Gradient NON 147 : *******PROJ TPNM2 (vrai)-> elimination end of step velocity OUI 148 : *******PROJ TMDM1 = VRAI TGRAD = FAUX TPNM2 = FAUX ancien algo 149 : TYPROJ='VPI1' ; 150 : ROW = 1. ; 151 : 'SI' ('EXIST' RV 'TYPROJ') ; 152 : TYPROJ=RV.'TYPROJ' ; 153 : 'SI' (non (exist (MOTS 'PSCT' 'VPI1' 'VPI2' 'PENA') TYPROJ)); 154 : Mess '*********************************************************' ; 155 : Mess ' ERREUR ERREUR ERREUR ERREUR ERREUR ERREUR ' ; 156 : Mess ' ' ; 157 : Mess 'Le mot ' TYPROJ ' n existe pas dans la liste.' ; 158 : Mess 'Les méthodes de projection autorisées sont :' ; 159 : Mess 'VPI1, VPI2, PSCT et PENA' ; 160 : Mess '*********************************************************' ; 161 : erreur 21 ; 162 : quitter EXEC ; 163 : 'FINSI' ; 164 : 'FINSI' ; 165 : 166 : 'SI' ('EGA' TYPROJ 'PSCT') ; 167 : TGRAD = FAUX ; 168 : TMDM1 = FAUX ; 169 : TPNM2 = FAUX ; 170 : 'FINSI' ; 171 : 'SI' ('EGA' TYPROJ 'VPI1') ; 172 : TGRAD = FAUX ; 173 : TMDM1 = VRAI ; 174 : TPNM2 = FAUX ; 175 : 'FINSI' ; 176 : 'SI' ('EGA' TYPROJ 'VPI2') ; 177 : TGRAD = FAUX ; 178 : TMDM1 = FAUX ; 179 : TPNM2 = VRAI ; 180 : 'FINSI' ; 181 : 'SI' ('EGA' TYPROJ 'PENA') ; 182 : TGRAD = FAUX ; 183 : TMDM1 = FAUX ; 184 : TPNM2 = FAUX ; 185 : 'FINSI' ; 186 : 187 : *******PROJ Quatres possibilités *************************************** 188 : *******PROJ************************************************************* 189 : 190 : nomvi=rv . 'NOMVI' ; 191 : 192 : 'SI' ('EGA' ('VALEUR' 'DIME') 2) ; 193 : vnul=0.D0 0.D0 ; 194 : vuni=1.D0 1.D0 ; 195 : cnvi1= 'MOT' ('TEXTE' ('CHAINE' 1 nomvi)) ; 196 : cnvi2= 'MOT' ('TEXTE' ('CHAINE' 2 nomvi)) ; 197 : lc= mots cnvi1 cnvi2 ; 198 : lcu=mots 'UX' 'UY' ; 199 : 'SINON' ; 200 : vnul=0.D0 0.D0 0.D0 ; 201 : vuni=1.D0 1.D0 1.D0 ; 202 : cnvi1= 'MOT' ('TEXTE' ('CHAINE' 1 nomvi)) ; 203 : cnvi2= 'MOT' ('TEXTE' ('CHAINE' 2 nomvi)) ; 204 : cnvi3= 'MOT' ('TEXTE' ('CHAINE' 3 nomvi)) ; 205 : lc= mots cnvi1 cnvi2 cnvi3 ; 206 : lcu=mots 'UX' 'UY' 'UZ' ; 207 : 'FINSI' ; 208 : 'SI' ('NON' ('EXISTE' rv 'OMEGA')) ; 209 : omeg=1.D0 ; 210 : 'SINON' ; 211 : omeg=rv . 'OMEGA' ; 212 : 'FINSI' ; 213 : 214 : testpr ='EXISTE' rv 'PRESSION' ; 215 : testprj ='EXISTE' rv 'PROJ' ; 216 : testran=testpr 'ET' ('EXISTE' rv 'CO') ; 217 : 218 : si testpr ; rvp = rv.'PRESSION' ; 219 : achp matpr= 'KOPS' 'MATRIK' ; 220 : rvp.'MATP'= matpr; 221 : finsi ; 222 : si testprj; rvp = rv.'PROJ' ; Finsi ; 223 : 224 : 'SI' (rv.'IMPR' >EG 1) ; 225 : si testpr ; 226 : mess 'Algorithme semi explicite (ANCIEN)'; 227 : mess '=================================='; 228 : Finsi ; 229 : si testprj; 230 : mess 'Algorithme de Projection'; 231 : mess '========================'; 232 : finsi ; 233 : si ( non (testpr ou testprj)) ; 234 : mess 'Algorithme standard implicite ou explicite '; 235 : mess '==========================================='; 236 : finsi ; 237 : 'FINSI' ; 238 : 239 : 'SI' ('NON' ('EXISTE' rv 'HIST')) ; 240 : rv . 'HIST' = 'TABLE' ; 241 : 'FINSI' ; 242 : 243 : 244 : ITMA=(rv . 'ITMA') ; 245 : IMPTCRR=0 ; 246 : IMPKRES=0 ; 247 : 'SI' ('<EG' ITMA 1) ; 248 : ITMA=1 ; 249 : IMPTCRR=1 ; 250 : 'FINSI' ; 251 : Si (testprj ) ; IMPTCRR=1 ; finsi ; 252 : Si ( non (testpr ou testprj)) ; IMPTCRR=1 ; finsi ; 253 : 254 : 255 : * Gestion de la matrice de préconditionnement 256 : * 257 : * calprec : doit-on recalculer le préconditionneur 258 : * fcprect : fréquence de recalcul du préc. en fn du pas de temps 259 : * fcpreci : fréquence de recalcul du préc. dans la boucle 260 : * d'itérations internes pour les non-linéarités 261 : * fcprectp : idem pour la matrice de pression 262 : * fcprecip : idem pour la matrice de pression 263 : * resmn : le résidu au pas de temps ou à l'itération interne précédente 264 : * 265 : calprec = VRAI ; 266 : fcprect = rv. 'METHINV' . 'FCPRECT' ; 267 : fcpreci = rv. 'METHINV' . 'FCPRECI' ; 268 : resmn maprec = 'KOPS' 'MATRIK' ; 269 : resmn1 maprec1 = 'KOPS' 'MATRIK' ; 270 : resmn2 maprec2 = 'KOPS' 'MATRIK' ; 271 : 'SI'(exist rv 'resmn');resmn=rv.'resmn' ; finsi ; 272 : 'SI' ('NON' ('EXIST' rv.'METHINV' 'CALPREC')) ; 273 : rv.'METHINV'.'CALPREC'= VRAI ; 274 : 'FINSI' ; 275 : 'SI' testprj ; 276 : fcprectp = rvp. 'METHINV' . 'FCPRECT' ; 277 : fcprecip = rvp. 'METHINV' . 'FCPRECI' ; 278 : 'SI' ('NON' ('EXIST' rvp.'METHINV' 'CALPREC')) ; 279 : rvp.'METHINV'.'CALPREC'= VRAI ; 280 : 'FINSI' ; 281 : 'FINSI' ; 282 : 283 : 'REPETER' bloc1 ITMA ; 284 : 'SI' ('MULT' &bloc1 fcprect) ; 285 : rv.'METHINV'.'CALPREC' = VRAI ; 286 : 'FINSI' ; 287 : 'SI' (testprj) ; 288 : 'SI' ('MULT' &bloc1 fcprectp) ; 289 : rvp.'METHINV'.'CALPREC' = VRAI ; 290 : 'FINSI' ; 291 : 'FINSI' ; 292 : testp1 = FAUX ; 293 : si (testpr ou testprj) ; 294 : Si CALPRE ; testp1=VRAI ; Finsi ; 295 : Si (non (exist rvp 'MATC')) ; testp1=VRAI ; Finsi ; 296 : finsi ; 297 : 298 : 'REPETER' bloci (rv . 'NITER') ; 299 : rv . 'NUITER' = &bloci ; 300 : 'SI' ('MULT' &bloci fcpreci) ; 301 : rv.'METHINV'.'CALPREC' = VRAI ; 302 : 'FINSI' ; 303 : 304 : 'SI' (testprj) ; 305 : 'SI' ('MULT' &bloci fcprecip) ; 306 : rvp.'METHINV'.'CALPREC' = VRAI ; 307 : 'FINSI' ; 308 : 'FINSI' ; 309 : st mat = 'KOPS' 'MATRIK' ; 310 : sf mau = 'KOPS' 'MATRIK' ; 311 : 312 : mdfdt = 0 ; 313 : 'REPETER' bloc2 ('DIME' (rv . 'LISTOPER')) ; 314 : nomper = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ; 315 : notable= 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ; 316 : 317 : * mess 'Operateur Bloc2 mdfdt ? ' nomper ; 318 : 319 : rvn=rv . notable; 320 : rvna=TBPART. notable; 321 : TASSI=FAUX; 322 : Si('EXIST' rvna 'SOUSTYPE'); 323 : Si('EGA' (rvna.'SOUSTYPE') 'ESCLAVE'); 324 : Si(Npart > 1); TASSI=VRAI; 325 : finsi; 326 : finsi; 327 : finsi; 328 : mdfdt = mdfdt + rvn . 'KOPT' . 'KFORM' ; 329 : 330 : si (ega nomper 'DFDT '); 331 : ISCHT= (rvn.kopt.'ISCHT') ; 332 : Si TASSI ; 333 : msi mai=ASSI TOUS ('TEXTE' nomper) (rvna) ; 334 : msi=et msi; 335 : mai=et mai; 336 : Sinon; 337 : msi mai=('TEXTE' nomper) (rvn) ; 338 : Finsi; 339 : mat = mat 'ET' mai ; 340 : st = st 'ET' msi ; 341 : 342 : sinon ; 343 : 344 : * mess ' 1OPER =' nomper; 345 : Si TASSI ; 346 : msi mai=ASSI TOUS ('TEXTE' nomper) (rvna) ; 347 : msi=et msi; 348 : mai=et mai; 349 : Sinon; 350 : msi mai=('TEXTE' nomper) (rvn) ; 351 : Finsi; 352 : mau = mau 'ET' mai ; 353 : sf = sf 'ET' msi ; 354 : 355 : finsi ; 356 : 357 : 'FIN' bloc2 ; 358 : 359 : s2 = sf et st ; 360 : ma1 = mau 'ET' mat ; 361 : 362 : ********** Traitement PRESSION **************************** 363 : *mess 'Traitement PRESSION'; 364 : 365 : 'SI' (exist rv 'rvpd') ; 366 : rvpd = rv.'rvpd' ; 367 : IDigv=rv.'IDigv'; 368 : Digv=rv.'Digv'; 369 : matpr= rvp.'MATP' ; 370 : 'FINSI' ; 371 : 372 : */1 CALCUL de CMCT pour PRESSION 373 : 'SI' (testpr et testp1) ; 374 : 375 : rvpd = (rvp.'DOMAINE') ; 376 : Diago = doma rvpd 'XXDIAGSI' ; 377 : si ( ega ('VALEUR' 'DIME') 2) ; 378 : Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) ; 379 : sinon ; 380 : Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) et 381 : ( exco Diago 'SCAL' cnvi3 ) ; 382 : finsi ; 383 : Digv = kcht rvpd vect sommet comp lc Digv ; 384 : 385 : 'SI' ('EXISTE' rv 'CLIM') ; 386 : rvp . 'CLIM' = rv . 'CLIM' ; 387 : Digv = 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1; 388 : 'FINSI' ; 389 : 390 : IDigv= 'INVERSE' digv ; 391 : 392 : 'SI' ('NON' ('EXISTE' rvp 'DIAGV')) ; 393 : rvp . 'DIAGV' = digv ; 394 : 'FINSI' ; 395 : 396 : rvp . 'MATC' = 'KMAB' rvp ; 397 : rvp . 'PRESSION' = 'KCHT' rvpd 'SCAL' 398 : 'CENTRE' 0.D0 ; 399 : rvp . 'GRADP' = 'KCHT' rvpd 'VECT' 400 : 'SOMMET' vnul ; 401 : 402 : achp matpr= 'KOPS' 'MATRIK' ; 403 : rvp . 'MATP' = matpr ; 404 : rv.'rvpd'=rvpd; 405 : rv.'IDigv'=IDigv; 406 : rv.'Digv'=Digv; 407 : 408 : 'FINSI' ; 409 : 410 : */3 CALCUL de CMCT pour PROJ 411 : 'SI' (testprj et testp1) ; 412 : 413 : Diago mma = 'KOPS' 'MATRIK' ; 414 : 415 : idfdt = 0 ; 416 : 'REPETER' blocj ('DIME' (rv . 'LISTOPER')) ; 417 : nomper = 'EXTRAIRE' &blocj (rv . 'LISTOPER') ; 418 : notable= 'MOT' ('TEXTE' ('CHAINE' &blocj nomper)) ; 419 : * mess 'Operateur-> ' nomper ; 420 : si ((ega nomper 'DFDT ') et 421 : ('EXIST' (rv.notable) 'LISTINCO')); 422 : si ('EXIST' (rv. notable.'LISTINCO') nomvi); 423 : idfdt=idfdt + 1 ; 424 : Diago=Diago et (doma (rv. notable .'DOMZ') 'XXDIAGSI'); 425 : rvpd=rv. notable . 'DOMZ'; 426 : finsi ; 427 : finsi ; 428 : 'FIN' blocj ; 429 : si (ega idfdt 0) ; mess ' Pas de DFDT ?? ' ; erreur 21 ; finsi ; 430 : 431 : si ( ega ('VALEUR' 'DIME') 2) ; 432 : Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) ; 433 : sinon ; 434 : Digv= ( exco Diago 'SCAL' cnvi1 ) et ( exco Diago 'SCAL' cnvi2 ) et 435 : ( exco Diago 'SCAL' cnvi3 ) ; 436 : finsi ; 437 : * Digv = kcht rvpd vect sommet comp lc (Digv*Ro) ; 438 : Digv = kcht rvpd vect sommet comp lc Digv ; 439 : * Matrice diagonale ne contenant pas les conditions aux limites : Diag 440 : Diag = copier Digv; 441 : IDiag = 'INVERSE' Diag ; 442 : 443 : Digv = 'KOPS' Digv 'CLIM' (rv . 'CLIM') -1; 444 : 445 : IDigv= 'INVERSE' digv ; 446 : 447 : 'SI' ('NON' ('EXISTE' rvp 'DIAGV')) ; 448 : rvp . 'DIAGV' = digv ; 449 : 'FINSI' ; 450 : 451 : rvp.'INCO'=rv.'INCO' ; 452 : 453 : sp map = 'KOPS' 'MATRIK' ; 454 : sr mar = 'KOPS' 'MATRIK' ; 455 : svnpc mvnpc = 'KOPS' 'MATRIK' ; 456 : 'REPETER' blocpj ('DIME' (rvp . 'LISTOPER')) ; 457 : nomper = 'EXTRAIRE' &blocpj (rvp . 'LISTOPER') ; 458 : notable= 'MOT' ('TEXTE' ('CHAINE' &blocpj nomper)) ; 459 : * mess 'Procedure PROJ Operateur : ' nomper ; 460 : 461 : si (EGA nomper 'KBBT'); 462 : rvp . notable . 'KOPT' . 'IKOMP' = 1 ; 463 : finsi ; 464 : 465 : msi mai=('TEXTE' nomper) (rvp . notable) ; 466 : 467 : *******PROJ************************************************************* 468 : *******PROJ TGRAD calcul maig (Formulation en Gradient) DEBUT*** ******* 469 : Si TGRAD ; 470 : mess 'Procedure PROJCT (Gradient) Operateur : ' nomper ; 471 : rvp . notable . 'KOPT' . 'IKOMP' = 0 ; 472 : msig maig=('TEXTE' nomper) (rvp . notable) ; 473 : Finsi ; 474 : *******PROJ -> MATG FIN***** 475 : *******PROJ///////////////////////////////////////////////////////////// 476 : 477 : TVNP=FAUX; 478 : TVNPC=FAUX; 479 : Si (('EGA' nomper 'VNIMP ') et (non (ega (rvp.'DISCPRES') 5))); 480 : TVNP=VRAI ; 481 : mar = mar 'ET' mai ; 482 : sr = sr 'ET' msi ; 483 : Finsi ; 484 : 485 : * mess 'nomper=' nomper (rvp.'DISCPRES'); 486 : Si (('EGA' nomper 'VNIMP ') et (ega (rvp.'DISCPRES') 5)); 487 : TVNPC=VRAI; 488 : mvnpc=mvnpc et mai ; 489 : svnpc=svnpc et msi ; 490 : Finsi ; 491 : 492 : Si (non ('EGA' nomper 'VNIMP ')); 493 : map = map 'ET' mai ; 494 : sp = sp 'ET' msi ; 495 : Finsi ; 496 : 497 : 'FIN' blocpj ; 498 : 499 : matpc mac = kops 'CMCTSPLT' map ; 500 : *******PROJ************************************************************* 501 : *******PROJ TGRAD (vrai Formulation en Gradient) CMCTSPLT DEBUT*** ******* 502 : Si TGRAD ; 503 : matpcg macg = kops 'CMCTSPLT' maig ; 504 : Finsi ; 505 : *******PROJ -> MATG FIN***** 506 : *******PROJ///////////////////////////////////////////////////////////// 507 : 508 : Si TVNPC; mac=mac et mvnpc; Finsi; 509 : rvp . 'MATC' = mac ; 510 : 511 : *******PROJ************************************************************* 512 : *******PROJ TGRAD (vrai Formulation en Gradient) -> MATG DEBUT*** ******* 513 : SI TGRAD ; rvp . 'MATG' = macg; FINSI ; 514 : *******PROJ -> MATG FIN***** 515 : *******PROJ///////////////////////////////////////////////////////////// 516 : 517 : scr mcr = 'KOPS' 'MATRIK' ; 518 : Si TVNP ; 519 : rvp . 'MBTR' = mar ; 520 : Dunit=Idigv; 521 : crt= kops 'CMCT' mac mar (Dunit) ; 522 : ctr= kops 'CMCT' mar mac (Dunit) ; 523 : rrt= kops 'CMCT' mar mar (Dunit) ; 524 : mcr= mcr et crt et ctr et rrt ; 525 : Finsi ; 526 : rvp . 'TVNP' = TVNP ; 527 : 528 : Si( ega (rvp.'DISCPRES') 5 ) ; 529 : mess ' Cas des pressions continues' ; 530 : matpr = matpc et mcr ; 531 : Sinon ; 532 : mess ' Cas des pressions discontinues' ; 533 : 534 : matpr = kops 'CMCT' mac mac IDigv ; 535 : matpr = matpc et matpr et mcr ; 536 : 537 : Finsi ; 538 : 539 : 540 : rvp . 'MATP' = matpr ; 541 : rv.'rvpd'=rvpd; 542 : rv.'IDigv'=IDigv; 543 : rv.'Digv'=Digv; 544 : 'FINSI' ; 545 : ************* 546 : ** Fin calcul CMCT pour testp1 et testprj 547 : ****************************************** 548 : * t 549 : *************** Calcul de C p Pression 550 : 'SI' testpr ; 551 : 552 : dt = (rv . 'PASDETPS' . 'DELTAT') '*' (rv . 'ALFA') ; 553 : rvp . 'DELTAT' = dt ; 554 : f = 'COPIER' s2 ; 555 : u = rv . 'INCO' . nomvi ; 556 : 557 : lc = 'EXTRAIRE' digv 'COMP' ; 558 : fu = 'KCHT' rvpd 'VECT' 'SOMMET' 559 : 'COMP' lc ('EXCO' f lc) ; 560 : 561 : 'SI' ('EXISTE' rv 'CLIM') ; 562 : dti = -1.D0 '/' dt ; 563 : dm1f= 'KOPS' (dt '*' ('KOPS' fu '*' IDigv)) 564 : 'CLIM' (rv . 'CLIM') 3 ; 565 : dm1f= dti * dm1f ; 566 : 'SINON' ; 567 : dm1f= (-1.D0) '*' ('KOPS' fu '*' IDigv) ; 568 : 'FINSI' ; 569 : 570 : rvp . 'PRESSION' = 'KMF' (rvp . 'MATC') dm1f ; 571 : 572 : 'KRES' rvp (rvp . 'PRESSION') 573 : 'BETA' (rvp . 'KBETA') (rvp . 'BETA') 574 : 'PIMP' (rvp . 'KPIMP') (rvp . 'PIMP') ; 575 : rvp . 'GRADP' = 'KMTP' 1 (rvp . 'MATC') 576 : (rvp . 'PRESSION') lc ; 577 : s2 = s2 + (rvp . 'GRADP') ; 578 : rv.'INCO'.'PRESSION'=rvp . 'PRESSION'; 579 : 580 : 'FINSI' ; 581 : 582 : *******PROJ************************************************************* 583 : *******PROJ MDM1 Ici on calcule M D-1 pour ensuite calculer Ctp DEBUT*** ******* 584 : *******PROJ on ne calcule M D-1 que si (rv 'MDM1') n'existe pas DEBUT*** ******* 585 : *******PROJ ou CALPRE VRAI DEBUT*** ******* 586 : 'SI' testprj ; 587 : 588 : dt = (rv . 'PASDETPS' . 'DELTAT') '*' (rv . 'ALFA') ; 589 : rvp . 'DELTAT' = dt ; 590 : 591 : ** Produit M D-1 592 : 'SI' TMDM1 ; 593 : 'SI' (('EXIST' rv 'MDM1') et (NON CALPRE)) ; 594 : MDM1=rv.'MDM1' ; 595 : 'SINON' ; 596 : mess ' On calcule MD-1 ' ; 597 : stn matn= 'KOPS' 'MATRIK' ; 598 : 599 : ROW = 0.; 600 : 'REPETER' blocpj1 ('DIME' (rvp . 'LISTOPER')) ; 601 : nomper1 = 'EXTRAIRE' &blocpj1 (rvp . 'LISTOPER') ; 602 : notable1= 'MOT' ('TEXTE' ('CHAINE' &blocpj1 nomper1)) ; 603 : * mess 'Operateur ' nomper1 ; 604 : 605 : si (ega nomper1 'KBBT '); 606 : nomiv1= extr (rvp . notable1 . 'LISTINCO') 1 ; 607 : 608 : 'REPETER' blocpj2 ('DIME' (rv . 'LISTOPER')) ; 609 : nomper2 = 'EXTRAIRE' &blocpj2 (rv . 'LISTOPER') ; 610 : notable2= 'MOT' ('TEXTE' ('CHAINE' &blocpj2 nomper2)) ; 611 : * mess 'Operateur ' nomper2 ; 612 : 613 : si (ega nomper2 'DFDT '); 614 : nomiv2= extr (rv . notable2 . 'LISTINCO') 1 ; 615 : si ('EGA' nomiv2 nomiv1) ; 616 : 617 : Si (EGA TYPROJ 'PENA'); 618 : rwi=(rv . notable2) . ARG1; 619 : trwi = type rwi ; mess ' type du coeff rwi ' trwi; 620 : 621 : Si (EGA trwi 'FLOTTANT'); 622 : ROW = maxi (prog rwi ROW); 623 : Finsi ; 624 : Si ((EGA trwi 'CHPOIN') ou (EGA trwi 'MCHAML')) ; 625 : ROW = maxi (prog (maxi rwi) ROW); 626 : Finsi ; 627 : Si (EGA (type rwi) 'MOT'); 628 : rwi = rv . 'INCO' . rwi ; 629 : Si (EGA (type rwi) 'FLOTTANT'); 630 : ROW = maxi (prog rwi ROW); 631 : Finsi ; 632 : Si ((EGA trwi 'CHPOIN') ou (EGA trwi 'MCHAML')) ; 633 : ROW = maxi (prog (maxi rwi) ROW); 634 : Finsi ; 635 : Finsi ; 636 : Sinon ; 637 : ROW = 1. ; 638 : Finsi ; 639 : * mess 'ROW = ' ROW ; 640 : 641 : domzp=rv . notable2 . 'DOMZ' ; 642 : msi mai='DFDT' (rv . notable2) ; 643 : matn = matn et mai ; 644 : sti = kcht domzp vect sommet comp lc vuni; 645 : finsi ; 646 : finsi ; 647 : 'FIN' blocpj2 ; 648 : finsi ; 649 : 'FIN' blocpj1 ; 650 : 651 : Si (EGA TYPROJ 'PENA'); 652 : MDM1=1. ; 653 : SINON ; 654 : MDM1='KMF' matn IDiag ; 655 : FINSI ; 656 : 657 : Si(EGA ISCHT 1); 658 : MDM1=MDM1 * ((dt*2.)/3.) ; 659 : Sinon ; 660 : MDM1=MDM1 * dt ; 661 : Finsi ; 662 : 663 : rv.'MDM1' = MDM1 ; 664 : 665 : *******PROJ MDM1 -> rv.'MDM1' FIN***** 666 : *******PROJ///////////////////////////////////////////////////////////// 667 : 'FINSI'; 668 : 'FINSI'; 669 : *******PROJ************************************************************* 670 : *******PROJ t n-1 DEBUT*** ******* 671 : *******PROJ Calcul de C P DEBUT*** ******* 672 : 673 : TVNP=rvp.'TVNP' ; 674 : 675 : 'SI' ('EXIST' (rv.'INCO') 'PRESSION') ; 676 : PPI = rv.'INCO'.'PRESSION' ; 677 : 'SI' ('EXIST' (rv.'INCO') 'PNM2') ; 678 : * TPNM2 elimination of end of step velocity (2 Pn - Pn-1) 679 : PPI = 2*(PPI) - rv.'INCO'.'PNM2' ; 680 : 'FINSI' ; 681 : 682 : *Formulation en u grad p 683 : Si (Exist rvp 'MATG'); 684 : cpre = 'KMF' (rvp . 'MATG') PPI 'TRAN' ; 685 : *Formulation en p div u 686 : Sinon; 687 : cpre = 'KMF' (rvp . 'MATC') PPI 'TRAN' ; 688 : Finsi ; 689 : 690 : Si TVNP ; 691 : cxre = 'KMF' (rvp . 'MBTR') PPI 'TRAN' ; 692 : cpre = cpre et cxre ; 693 : Finsi ; 694 : 695 : *******PROJ MDM1 696 : Si TMDM1 ; 697 : * consistence selon (Gresho) 698 : gradpres = 'KOPS' MDM1 '*' cpre ; 699 : Sinon ; 700 : * consistence selon (Guermond) 701 : gradpres = cpre ; 702 : Finsi ; 703 : 704 : oublier cpre ; 705 : rv.'INCO'.'GRADPRES' = 'NOMC' lcu lc gradpres; 706 : 707 : *Formulation en u grad p 708 : Si (Exist rvp 'MATG'); 709 : s2 = sf - (rv.'INCO'.'GRADPRES') + st ; 710 : Sinon ; 711 : *Formulation en p div u 712 : s2 = sf + (rv.'INCO'.'GRADPRES') + st ; 713 : Finsi ; 714 : 715 : 'FINSI' ; 716 : 717 : *******PROJ t n-1 FIN***** 718 : *******PROJ MDM1 '*' C P -> rv.'INCO'.'GRADPRES' FIN***** 719 : *******PROJ t n-1 FIN***** 720 : *******PROJ si (p div u) -> S2 = F + C P + st FIN***** 721 : *******PROJ t n-1 FIN***** 722 : *******PROJ si (u grad p) -> S2 = F - C P + st FIN***** 723 : *******PROJ///////////////////////////////////////////////////////////// 724 : 'FINSI' ; 725 : 726 : ********** FIN Traitement PRESSION ************************ 727 : * mess ' FIN Traitement PRESSION' ; 728 : 729 : 'SI' ('EXISTE' rv 'CLIM') ; 730 : s1 = rv . 'CLIM' ; 731 : 'SINON' ; 732 : s1=' ' ; 733 : 'FINSI' ; 734 : rv . 'S2' = s2 ; 735 : 736 : ******************************************************************* 737 : * Résolution hors QDM méthode de projection 738 : 'SI'((NON testprj) ou (EGA mdfdt 0)); 739 : *mess '* Résolution hors QDM méthode de projection '; 740 : 741 : rv . 'METHINV' . 'XINIT' = resmn ; 742 : 743 : 'SI' ('EXISTE' rv 'GPROJ') ; 744 : res = GRESP ma1 s1 s2 rv ; 745 : 'SINON' ; 746 : res = KRESP ma1 'TYPI' (rv . 'METHINV') 747 : 'CLIM' s1 748 : 'SMBR' s2 749 : 'IMPR' IMPKRES ; 750 : 'FINSI' ; 751 : 752 : 'FINSI' ; 753 : 754 : 'SI' (testprj et (non(EGA mdfdt 0))); 755 : *******PROJ************************************************************* 756 : *******PROJ Résolution QDM DEBUT*** ******* 757 : *******PROJ (On éclate les résolutions) DEBUT*** ******* 758 : 759 : * Résolution QDM méthode de projection (On éclate les résolutions) 760 : *mess ' Résolution QDM méthode de projection'; 761 : 762 : lpart = KOPS 'EXTRCOUP' ma1 ; 763 : nbpart = dime lpart ; 764 : 765 : $ma1 =table 'ESCLAVE'; 766 : $s1 =table 'ESCLAVE'; 767 : $s2 =table 'ESCLAVE'; 768 : $resmn=table 'ESCLAVE'; 769 : $tab1 =table 'ESCLAVE'; 770 : 771 : repeter Bclcom nbpart ; 772 : nmc=lpart.&Bclcom; 773 : * mess ' Liste des composantes'; 774 : * list nmc ; 775 : MA1i = kops 'EXTRINCO' nmc nmc ma1 ; 776 : $ma1 .&Bclcom = MA1i ; 777 : $s1 .&Bclcom = exco s1 nmc 'NOID' ; 778 : $s2 .&Bclcom = exco s2 nmc 'NOID' ; 779 : resmni = exco resmn nmc 'NOID' ; 780 : TAB1 = copie rv . 'METHINV' ; 781 : * TAB1 = rv . 'METHINV' ; 782 : TAB1 . 'XINIT' = resmni ; 783 : $tab1 .&Bclcom = TAB1 ; 784 : 785 : *------------------------------------------------------------------------- 786 : *mess ' Mise à jour du préconditionnement'; 787 : 'SI' ('>EG' TAB1 . 'TYPINV' 2) ; 788 : lword1 lword2 = 'KOPS' 'EXTRNINC' MA1i ; 789 : word1 = 'EXTR' lword1 1 ; 790 : 791 : 'SI' ('NON' ('EXIS' rv 'TABRES')) ; 792 : rv . 'TABRES' = 'TABLE' ; 793 : 'FINSI' ; 794 : TABRES=rv.'TABRES'; 795 : 796 : 'SI' (TAB1 . 'CALPREC') ; 797 : 'SI' ('NON' ('EXIS' TABRES word1)) ; 798 : TABRES . word1 = 'TABLE' ; 799 : TABRES . word1 . 'MATASS' = MA1i ; 800 : TABRES . word1 . 'MAPREC' = MA1i ; 801 : 'MESS' 'On recalcule le preconditionneur'; 802 : 'SINON' ; 803 : TABRES . word1 . 'MATASS' = MA1i ; 804 : TABRES . word1 . 'MAPREC' = MA1i ; 805 : 'MESS' 'On recalcule le preconditionneur'; 806 : 'FINSI' ; 807 : 'FINSI' ; 808 : TAB1 . 'MATASS' = TABRES.word1.'MATASS' ; 809 : TAB1 . 'MAPREC' = TABRES.word1.'MAPREC' ; 810 : 'FINSI' ; 811 : *------------------------------------------------------------------------- 812 : Fin Bclcom ; 813 : 814 : Si(EXIST RV 'KPR'); 815 : KPR=rv.'KPR'; 816 : Sinon; 817 : KPR=FAUX; 818 : Finsi; 819 : 820 : Si KPR ; 821 : mess ' On parallélise les résolutions '; 822 : res = ASSI TOUS KRES $ma1 'TYPI' $tab1 823 : 'CLIM' $s1 824 : 'SMBR' $s2 825 : 'IMPR' IMPKRES ; 826 : res=et res ; 827 : Sinon ; 828 : mess ' Les résolutions sont traitées séquentiellement'; 829 : res mm1= kops 'MATRIK' ; 830 : repeter Bclcom nbpart ; 831 : resi= KRES ($ma1.&bclcom) 'TYPI' ($tab1.&bclcom) 832 : 'CLIM' ($s1.&bclcom) 833 : 'SMBR' ($s2.&bclcom) 834 : 'IMPR' IMPKRES ; 835 : res=res et resi ; 836 : Fin Bclcom ; 837 : Finsi; 838 : 839 : *mess ' Fin Résolution QDM méthode de projection'; 840 : * Fin Résolution QDM méthode de projection 841 : *------------------------------------------------------------------------------- 842 : *******PROJ Résolution QDM -> res FIN***** 843 : *******PROJ (On éclate les résolutions) FIN***** 844 : *******PROJ///////////////////////////////////////////////////////////// 845 : 'FINSI' ; 846 : 847 : 'SI' testprj ; 848 : *******PROJ************************************************************* 849 : *******PROJ ETAPE DE PROJECTION DEBUT*** ******* 850 : *******PROJ on calcul cun (alias c U tilde) DEBUT*** ******* 851 : * mess 'ETAPE DE PROJECTION ' ; 852 : * _n 853 : * C U 854 : cun = ('KMF' (rvp . 'MATC') res) ; 855 : Si TVNP ; 856 : cxn = ('KMF' (rvp . 'MBTR') res) ; 857 : cun = cun et cxn ; 858 : Finsi ; 859 : 860 : * mess ' On calcule les seconds membres de l équation de pression ' ; 861 : * mess ' s ils existent (opérateurs FIMP) ' ; 862 : 863 : 'REPETER' blocpj ('DIME' (rvp . 'LISTOPER')) ; 864 : nomper = 'EXTRAIRE' &blocpj (rvp . 'LISTOPER') ; 865 : notable= 'MOT' ('TEXTE' ('CHAINE' &blocpj nomper)) ; 866 : si (ega nomper 'FIMP '); 867 : * mess 'Seconds membres de l équation de pression Opérateur ' nomper ; 868 : msi mai=('TEXTE' nomper) (rvp . notable) ; 869 : cun = cun et msi ; 870 : finsi ; 871 : 'FIN' blocpj ; 872 : 873 : cun = cun * ROW * (-1./dt) ; 874 : 875 : rvp . 'METHINV' . 'XINIT' = resmn2 ; 876 : 877 : sl= KRESP matpr 'TYPI' (rvp . 'METHINV') 878 : 'CLIM' (rvp . CLIM ) 879 : 'SMBR' cun 880 : 'IMPR' IMPKRES ; oublier cun ; 881 : 882 : oublier resmn2 ; resmn2 = sl ; 883 : ctl = 'KMF' (rvp . 'MATC') sl 'TRAN' ; 884 : Si TVNP ; 885 : cxl = 'KMF' (rvp . 'MBTR') sl 'TRAN' ; 886 : ctl = ctl + cxl ; 887 : Finsi ; 888 : 889 : * elimination of end of step velocity (Si exist PNM2) (on ne corrige pas) 890 : 'SI'(NON TPNM2); 891 : 'SI' (NON ('EXIST' (rv.'INCO') 'PNM2')) ; 892 : a=nomc lcu lc ('KOPS' IDigv '*' ctl) ; 893 : res = res + ((1./ROW)*a*dt) ; 894 : 'FINSI' ; 895 : 'FINSI' ; 896 : 897 : oublier ctl ; 898 : 899 : *******PROJ -> (alias c U tilde) FIN***** 900 : *******PROJ -> sl (alias lambda) FIN***** 901 : *******PROJ -> on corrige u FIN***** 902 : *******PROJ///////////////////////////////////////////////////////////// 903 : 'FINSI' ; 904 : 905 : ***** Avancement en temps 906 : *? Si (testprj et ('EGA' mdfdt 0)) ; IMPTCRR=RV.'IMPR' ; finsi ; 907 : Si ('EGA' mdfdt 0) ; IMPTCRR=RV.'IMPR' ; finsi ; 908 : IMPKRES=0 ; 909 : * mess 'On passe dans TCRR ' ; 910 : eps = 'TCRR' res omeg (rv . 'INCO') 'IMPR' IMPTCRR ; 911 : 912 : 'OUBLIER' resmn ; 913 : rv.'METHINV'.'CALPREC' = FAUX ; 914 : 'SI' (testprj) ; 915 : rvp.'METHINV'.'CALPREC'= FAUX ; 916 : 'FINSI' ; 917 : resmn = res ; 918 : 'MENAGE' ; 919 : * 920 : 'SI' (rv . 'STOPITER') ; 921 : rv . 'STOPITER' = FAUX ; 922 : 'QUITTER' bloci ; 923 : 'FINSI' ; 924 : 'FIN' bloci ; 925 : irt=0 ; 926 : 927 : * mess ' ITMA= ' (rv . 'ITMA') ' mdfdt= ' mdfdt ; 928 : 'SI' ('EGA' (rv . 'ITMA') 0) ; 929 : * mess ' ON PASSE PLUS DS TCNM>>>>>>>>>>>>>>' ; 930 : irt = 'TCNM' rv 'NOUP'; 931 : 'SINON' ; 932 : * mess ' ON PASSE DS TCNM>>>>' (rv . 'ITMA') mdfdt; 933 : irt = 'TCNM' rv ; 934 : 'FINSI' ; 935 : 936 : ***** Avancement en temps algorithme de projection 937 : *******PROJ************************************************************* 938 : *******PROJ P(n+1) = P(n) + sl DEBUT*** ******* 939 : * Algorithme de projection P(n+1) = P(n) + sl 940 : 'SI' testprj ; 941 : 942 : 'SI' (NON ('EXIST' (rv.'INCO') 'PRESSION')); 943 : rv.'INCO'.'PRESSION' = sl; 944 : Si TPNM2 ; rv.'INCO'.'PNM2' = sl ; Finsi ; 945 : 'SINON' ; 946 : PNM1=rv.'INCO'.'PRESSION'; 947 : si ( EGA ISCHT 1) ; 948 : PN=PNM1 + (1.5*sl) ; 949 : sinon ; 950 : PN=PNM1 + sl ; 951 : finsi ; 952 : 953 : Si TPNM2 ; rv.'INCO'.'PNM2'=PNM1 ; Finsi ; 954 : rv.'INCO'.'PRESSION'=PN ; 955 : 'FINSI' ; 956 : 957 : Si ((EGA TYPROJ 'PSCT') ou (EGA TYPROJ 'PENA')) ; 958 : rv.'INCO'.'PRESSION' = sl; 959 : Finsi ; 960 : 961 : 'OUBLIER' sl ; 962 : 'FINSI' ; 963 : *******PROJ P(n+1) = P(n) + sl FIN***** 964 : *******PROJ///////////////////////////////////////////////////////////// 965 : ***** Avancement en temps FIN 966 : 967 : 968 : 'SI' testran ; 969 : rv . 'CO' . 'VITESSE' ='KOPS' (rv . 'INCO' . nomvi) 970 : '-' (rv . 'SEDIM') ; 971 : k = 'ABS' (rv . 'INCO' . 'KN') ; 972 : e = 'ABS' (rv . 'INCO' . 'EN') ; 973 : k = 'KOPS' ('KOPS' k '*' k) '/' e ; 974 : dif = 'KOPS' ('KOPS' k '*' 0.09) '+' (rv . 'COEF') ; 975 : rv . 'CO' . 'DIFFU' = 'NOEL' rvpd dif ; 976 : rv . 'CO' . 'TEMPERA' = rv . 'INCO' . 'CN' ; 977 : 'FINSI' ; 978 : rv.'METHINV'.'CALPREC' = FAUX ; 979 : 'SI' (testprj) ; 980 : rvp.'METHINV'.'CALPREC'= FAUX ; 981 : 'FINSI' ; 982 : 'MENAGE' ; 983 : 984 : 'SI' ('EGA' irt 1) ; 985 : 'MESSAGE' ' Temps final atteint : ' 986 : (rv . 'PASDETPS' . 'TPS') ; 987 : 'QUITTER' bloc1 ; 988 : 'FINSI' ; 989 : 'SI' (rv . 'STOPPDT') ; 990 : rv . 'STOPPDT' = FAUX ; 991 : 'MESSAGE' ' Arret demandé au temps : ' (rv . 'PASDETPS' . 'TPS') ; 992 : 'QUITTER' bloc1 ; 993 : 'FINSI' ; 994 : 'FIN' bloc1 ; 995 : 996 : 'SI' testpr ; 997 : rvp . 'PRESSION' = 'KCHT' rvpd 'SCAL' 'CENTRE' 998 : (rvp . 'PRESSION') ; 999 : rvp . 'PN' = 'ELNO' rvpd 1000 : (rvp . 'PRESSION') ; 1001 : 'FINSI' ; 1002 : ************************ E X E C ************************************ 1003 : 'FINPROC' ; 1004 : 1005 : 1006 : 1007 : 1008 : 1009 :
© Cast3M 2003 - All rights reserved.
Disclaimer