1 : * GRESP PROCEDUR GOUNAND 12/02/01 21:15:14 7266 2 : ************************************************************************ 3 : * NOM : GRESP 4 : * DESCRIPTION : Résout un système par une méthode de projection 5 : * algébrique incrémentale 6 : * 7 : * LANGAGE : GIBIANE-CAST3M 8 : * AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) 9 : * mél : gounand@semt2.smts.cea.fr 10 : ********************************************************************** 11 : * VERSION : v1, 22/11/2007, version initiale 12 : * HISTORIQUE : v1, 22/11/2007, création 13 : * HISTORIQUE : 21/12/2011 : correction d'un bug dans le 14 : * préconditionnement, nettoyage KOPS 15 : * HISTORIQUE : 16 : ************************************************************************ 17 : * 18 : 'DEBPROC' GRESP ; 19 : 'ARGUMENT' mat*'MATRIK' ; 20 : 'ARGUMENT' ccli*'CHPOINT' ; 21 : 'ARGUMENT' smbr*'CHPOINT' ; 22 : 'ARGUMENT' rv*'TABLE' ; 23 : * 24 : * Nom des inconnues : vitesses et pressions 25 : * Dans rvg . 'METHINV', on stocke l'éventuelle table d'inversion 26 : * pour le laplacien de pression 27 : * 28 : debug = FAUX ; 29 : rvg = rv . 'GPROJ' ; 30 : rvm = rv . 'METHINV' ; 31 : 'SI' ('EXISTE' rvg 'dblproj') ; 32 : dblproj = rvg . 'dblproj' ; 33 : 'SINON' ; 34 : dblproj = VRAI ; 35 : 'FINSI' ; 36 : * 37 : 'SI' ('EXISTE' rvg 'NOPREC') ; 38 : noprec = rvg . 'NOPREC' ; 39 : 'SINON' ; 40 : noprec = faux ; 41 : 'FINSI' ; 42 : 'SI' debug ; 43 : 'MESSAGE' ('CHAINE' 'noprec=' noprec) ; 44 : 'FINSI' ; 45 : * 46 : 'SI' ('NON' noprec) ; 47 : 'SI' ('NON' ('EXISTE' rvg 'preco')) ; 48 : rvg . 'preco' = 'TABLE' ; 49 : 'FINSI' ; 50 : rvgp = rvg . 'preco' ; 51 : * Précision relative utilisée pour le préconditionnement du calcul de la 52 : * matrice de pression 53 : precrel = 1.D-12 ; 54 : * precrel = 1.D-2 ; 55 : 'FINSI' ; 56 : * Projection incrémentale donc pas besoin de XINIT 57 : 'SI' ('EXIS' rvm 'XINIT' ) ; 58 : 'OUBLIER' rvm 'XINIT' ; 59 : 'FINSI' ; 60 : *fd 'OUBLIER' rvm 'XINIT' ; 61 : * On a mis chaine sinon plantage qd la pression s'appelle 'PRES' ; 62 : ngvit = 'CHAINE' rvg . 'NOMVIT' ; 63 : ngpre = 'CHAINE' rvg . 'NOMPRES' ; 64 : 'SI' ('NON' ('EXISTE' rvg 'METHINV')) ; 65 : rvg . 'METHINV' = rvm ; 66 : 'FINSI' ; 67 : rvgm = rvg . 'METHINV' ; 68 : * 69 : dim = 'VALEUR' 'DIME' ; 70 : *fd 71 : * Si l'utilisateur a transmis les noms de composantes on ne fait rien. 72 : * Sinon (il utilise EXEC) on construit la liste nivit = 'MOTS' '1UN' '2UN' 73 : 'SI' ('EXIS' rvg 'COMPVIT') ; 74 : nivit = rvg . 'COMPVIT' ; 75 : nivit2 = nivit ; 76 : nipre = rvg . 'COMPPRES' ; 77 : nipre2 = nipre ; 78 : 'SINON' ; 79 : nivit = 'MOTS' ; 80 : 'REPETER' idim dim ; 81 : niviti = 'CHAINE' &idim ngvit ; 82 : nivit = 'ET' nivit ('MOTS' niviti) ; 83 : 'FIN' idim ; 84 : nipre = 'MOTS' ngpre ; 85 : nivit2 = 'EXTRAIRE' ('MOTS' 'UX' 'UY' 'UZ') ('LECT' 1 'PAS' 1 dim) ; 86 : nipre2 = 'MOTS' 'SCAL' ; 87 : 'FINSI' ; 88 : * Construit la liste des noms d'inconnues ni vitesse ni pression 89 : nitout = 'EXTR' mat 'COMP' ; 90 : niaut = 'MOTS' ; 91 : 'REPETER' itout ('DIME' nitout) ; 92 : ninc = 'EXTRAIRE' nitout &itout ; 93 : lexis = FAUX ; 94 : 'REPETER' ivit ('DIME' nivit) ; 95 : lexis = 'OU' lexis ('EGA' ninc ('EXTRAIRE' nivit &ivit)) ; 96 : 'FIN' ivit ; 97 : 'REPETER' ipre ('DIME' nipre) ; 98 : lexis = 'OU' lexis ('EGA' ninc ('EXTRAIRE' nipre &ipre)) ; 99 : 'FIN' ipre ; 100 : 'SI' ('NON' lexis) ; 101 : niaut = 'ET' niaut ('MOTS' ninc) ; 102 : 'FINSI' ; 103 : 'FIN' itout ; 104 : * 105 : * Initialisation de la solution du système total 106 : * 107 : sol = 0. ; 108 : * 109 : * On résout d'abord la partie ni vitesse ni pression 110 : * 111 : 'SI' ('>' ('DIME' niaut) 0) ; 112 : mtot = 'EXTRAIRE' mat niaut niaut ; 113 : ftot = 'EXCO' niaut smbr niaut 'NOID' ; 114 : cltot = 'EXCO' niaut ccli niaut 'NOID' ; 115 : rvm . 'MATASS' =mtot ; 116 : rvm . 'MAPREC' =mtot ; 117 : ct = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm ; 118 : sol = '+' ct sol ; 119 : 'FINSI' ; 120 : ******************************************************************** 121 : * * 122 : * On résout la partie vitesse-pression par projection incrémentale * 123 : * * 124 : ******************************************************************** 125 : nivp = 'ET' nivit nipre ; 126 : * Réduction des arguments de la proc à la vitesse-pression 127 : *mat = 'KOPS' 'EXTRINCO' mat nivp nivp ; 128 : mat = 'EXTRAIRE' mat nivp nivp ; 129 : ccli = 'EXCO' nivp ccli nivp 'NOID' ; 130 : smbr = 'EXCO' nivp smbr nivp 'NOID' ; 131 : * 132 : un = rv . 'INCO' . ngvit ; 133 : pn = rv . 'INCO' . ngpre ; 134 : un2 = 'NOMC' nivit2 nivit un ; 135 : pn2 = 'NOMC' nipre2 nipre pn ; 136 : incovp = un2 '+' pn2 ; 137 : vpetit = -1.D200 ; 138 : * Comme on va faire une méthode incrémentale, il faut imposer 139 : * les conditions aux limites de Dirichlet sur l'inconnue 140 : ccli1 = 'MASQUE' ccli 'SUPERIEUR' vpetit ; 141 : ccli0 = '-' ccli1 1.D0 ; 142 : niccli = 'EXTRAIRE' ccli 'COMP' ; 143 : dcli = '-' ccli ('*' incovp ccli1 niccli niccli niccli) ; 144 : incovp = '+' incovp dcli ; 145 : * 146 : * Matrice : partie diagonale en vitesse et contrainte 147 : * on est obligé de reconstruire la transposée car KBBT 148 : * ne stocke que B et pas Bt 149 : * 150 : mkvit = 'EXTRAIRE' mat nivit nivit ; 151 : mkcontr = 'EXTRAIRE' mat nivit nipre ; 152 : mkcontrt = 'KOPS' 'TRANSPOS' mkcontr ; 153 : Kt = mkvit 'ET' mkcontr 'ET' mkcontrt ; 154 : * 155 : * Construction du résidu * -1 156 : * 157 : residu = '-' smbr ('*' Kt incovp) ; 158 : * 159 : * Calcul de la matrice de pression et de la matrice masse diagonalisée 160 : * 161 : * On lumpe la partie diagonale et on l'inverse 162 : chpo1 = 'MASQUE' un2 'SUPERIEUR' vpetit ; 163 : mkvitd = '*' mkvit chpo1 ; 164 : * 165 : * Choix 1 : si la matrice masse lumpée 166 : * a tous ses termes positifs, on en prend l'inverse 167 : * Avantage : simple, autorise le préconditionnement 168 : * car ne change pas à chaque itération 169 : * Inconvénient : ne marche pas en axi quadratique sur l'axe 170 : * Choix 2 : on prend la diagonale de la matrice totale 171 : * on la met à l'échelle pour avoir la même masse que mkvitd 172 : * Avantage : marche en axi, sur maillage déformés 173 : * Inconvénient : change à chaque itération quand le problème 174 : * est non linéaire (=> recalcul de la matrice de pression), 175 : * pas forcément positif. 176 : * Choix 3 : valeur absolue du précédent + petite tolérance 177 : * Avantage : marche peut-être 178 : * Inconvénient : change à chaque itération. 179 : * 180 : * 181 : mmi = 'MINIMUM' mkvitd ; 182 : 'SI' ('>' mmi 1.D-200) ; 183 : imkvitd1 = 'INVERSE' mkvitd ; 184 : 'SI' debug ; 'MESSAGE' 'Choix 1 masse lumpée' ; 'FINSI' ; 185 : 'SINON' ; 186 : mkdia = 'EXTR' mkvit 'DIAG' ; 187 : mmi = 'MINIMUM' mkdia ; 188 : 'SI' ('>' mmi 1.D-200) ; 189 : 'SI' debug ; 'MESSAGE' 'Choix 2 masse lumpée' ; 'FINSI' ; 190 : 'SINON' ; 191 : mkdia = 'ABS' mkdia ; 192 : mma = 'MAXIMUM' mkdia ; 193 : mkdia = '+' mkdia ('*' mma 1.D-8) ; 194 : 'SI' debug ; 'MESSAGE' 'Choix 3 masse lumpée' ; 'FINSI' ; 195 : 'FINSI' ; 196 : masvitd = 'MAXIMUM' ('RESULT' mkvitd) ; 197 : * 'MESSAGE' ('CHAINE' 'totmas=' masvitd) ; 198 : masdia = 'MAXIMUM' ('RESULT' mkdia) ; 199 : imkvitd1 = '*' ('INVERSE' mkdia) ('/' masdia masvitd) ; 200 : 'FINSI' ; 201 : * On surcharge par 0. l'inverse de la diagonale là où il y a des clims 202 : * de Dirichlet 203 : imkvitd = '-' imkvitd1 ('*' imkvitd1 ccli1 niccli niccli niccli) ; 204 : * 205 : * Attention ! On préconditionne éventuellement le calcul de la matrice 206 : * de pression 207 : * 208 : precok = FAUX ; 209 : *'MESSAGE' ('CHAINE' 'precok = ' precok) ; 210 : 'SI' ('NON' noprec) ; 211 : exis1 = 'EXISTE' rvgp 'mklapphi' ; 212 : exis2 = 'EXISTE' rvgp 'imkvitd' ; 213 : 'SI' (exis1 'ET' exis2) ; 214 : imkvd = rvgp . 'imkvitd' ; 215 : mkl = rvgp . 'mklapphi' ; 216 : * On vérifie l'égalité des mkvitd à un facteur constant près 217 : ech = 'MAXIMUM' imkvitd 'ABS' ; 218 : echp = 'MAXIMUM' imkvd 'ABS' ; 219 : alfa = '/' ech echp ; 220 : * 'MESSAGE' ('CHAINE' 'alfa = ' alfa) ; 221 : dimkv = '-' imkvitd ('*' imkvd alfa) ; 222 : ndimkv = 'EXTRAIRE' dimkv 'COMP' ; 223 : inimkv = 'INVERSE' ('+' imkvitd ccli1) ; 224 : dimkvs = '*' dimkv inimkv ndimkv ndimkv ndimkv ; 225 : mdimkvs = 'MAXIMUM' dimkvs 'ABS' ; 226 : egaimkv = 'EGA' mdimkvs 0. precrel ; 227 : * 228 : 'SI' debug ; 229 : 'MESSAGE' ('CHAINE' 'egaimkv = ' egaimkv 230 : ' a ' ('*' mdimkvs 100.) ' %' ) ; 231 : 'FINSI' ; 232 : precok = egaimkv ; 233 : 'FINSI' ; 234 : 'FINSI' ; 235 : * 236 : chdbg = 'CHAINE' 'Matrice de pression MP : ' ; 237 : 'SI' precok ; 238 : 'SI' ('NON' ('EGA' alfa 1.D0 precrel)) ; 239 : chdbg = 'CHAINE' chdbg 'MP(i) = MP(i-1) * ' alfa ; 240 : * mklapphi = '*' (rvgp . 'mklapphi') alfa ; 241 : * imkvitd = '*' (rvgp . 'imkvitd') alfa ; 242 : mklapphi = rvgp . 'mklapphi' ; 243 : imkvitd = rvgp . 'imkvitd' ; 244 : 'SINON' ; 245 : chdbg = 'CHAINE' chdbg 'MP(i) = MP(i-1)' ; 246 : mklapphi = rvgp . 'mklapphi' ; 247 : imkvitd = rvgp . 'imkvitd' ; 248 : 'FINSI' ; 249 : 'SINON' ; 250 : * Dans le cas général où mkcontr et mkcontrt sont différents, 251 : * il faudrait mettre : 252 : *mklapphi = 'KOPS' 'CMCT' mkcontr imchd (kops mkcontrt 'TRANSPOS') ; 253 : chdbg = 'CHAINE' chdbg 'calcul...' ; 254 : mklapphi = 'KOPS' 'CMCT' mkcontr imkvitd mkcontr ; 255 : alfa = 1.D0 ; 256 : 'FINSI' ; 257 : 'SI' debug ; 'MESSAGE' chdbg ; 'FINSI' ; 258 : 'SI' ('NON' noprec) ; 259 : rvgp . 'mklapphi' = mklapphi ; 260 : rvgp . 'imkvitd' = imkvitd ; 261 : 'FINSI' ; 262 : *'MESSAGE' 'Fin du calcul de la matrice de pression' ; 263 : * 264 : * Calcul éventuel d'une estimation de la pression (double projection) 265 : * 266 : *'MESSAGE' ('CHAINE' 'dblproj=' dblproj) ; 267 : 'SI' dblproj ; 268 : desqdm = 'EXCO' nivit residu nivit 'NOID' ; 269 : idesq = '*' desqdm imkvitd nivit nivit nivit ; 270 : * 1) Celui-ci se compense avec le 2) 271 : * donc pas nécessaire 272 : * 'SI' ('NEG' alfa 1.D0 1.D-13) ; 273 : * idesq = '*' idesq alfa ; 274 : * 'FINSI' ; 275 : didesq = '*' mkcontr idesq ; 276 : mtot = mklapphi ; 277 : *dbg dmtot = 'EXTRAIRE' mklapphi 'DIAG' ; 278 : *dbg HCHPO dmtot 'dmtot' ; 279 : ftot = didesq ; 280 : cltot = 'EXCO' nipre ccli0 nipre 'NOID' ; 281 : rvgm .'MATASS' =mtot ; 282 : rvgm . 'MAPREC' =mtot ; 283 : *dbg dphi1 tt = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm 284 : *dbg 'IMPR' 3 'LTIME' VRAI ; 285 : dphi1 = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm ; 286 : * 2) Se compense avec le 1) 287 : * donc pas nécessaire 288 : * 'SI' ('NEG' alfa 1.D0 1.D-13) ; 289 : * dphi1 = '/' dphi1 alfa ; 290 : * 'FINSI' ; 291 : * HCHPO dphi1 'dphi1' ; 292 : gdphi1 = '*' mkcontrt dphi1 ; 293 : residu = residu '-' gdphi1 ; 294 : 'FINSI' ; 295 : * 296 : * Calcul de l'incrément de vitesse dv* à partir de la qdm 297 : * On splitte les équations par direction 298 : * 299 : dvitstar = 0. ; 300 : ncomp = 'DIME' nivit ; 301 : 'REPETER' iicomp ncomp ; 302 : icomp = &iicomp ; 303 : licomp = 'LECT' icomp ; 304 : mcomp = 'EXTRAIRE' nivit licomp ; 305 : mtot = 'EXTRAIRE' mkvit mcomp mcomp ; 306 : ftot = 'EXCO' mcomp residu mcomp 'NOID' ; 307 : cltot = 'EXCO' mcomp ccli0 mcomp 'NOID' ; 308 : rvm . 'MATASS' =mtot ; 309 : rvm . 'MAPREC' =mtot ; 310 : *dbg ct tt = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm 311 : *dbg 'IMPR' 3 'LTIME' VRAI ; 312 : ct = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm ; 313 : dvitstar = '+' dvitstar ct ; 314 : 'FIN' iicomp ; 315 : * HCHPO dvitstar 'dvitstar' ; 316 : * 317 : * Calcul de l'incrément de pression dphi pour vérifier la conservation 318 : * de la masse 319 : * 320 : * Pas besoin de mettre dvitstar à 0 sur les clims Dirichlet car c'est 321 : * déjà le cas ! 322 : smb1 = '*' mkcontr dvitstar ; 323 : mtot = mklapphi ; 324 : ftot = '-' smb1 ('EXCO' nipre residu nipre 'NOID') ; 325 : cltot = 'EXCO' nipre ccli0 nipre 'NOID' ; 326 : rvgm .'MATASS' =mtot ; 327 : rvgm . 'MAPREC' =mtot ; 328 : dphi = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm ; 329 : 'SI' ('NEG' alfa 1.D0 1.D-13) ; 330 : dphi = '/' dphi alfa ; 331 : 'FINSI' ; 332 : *HCHPO dphi 'dphi' ; 333 : * 334 : * Calcul de l'incrément final de vitesse dvit 335 : * 336 : gdphi = '*' mkcontrt dphi ; 337 : * On surcharge par 0. là où il y a des clims de Dirichlet 338 : * Pas besoin !!! vu que imkvitd est là ! 339 : *!! gdphi = '-' gdphi ('*' gdphi ccli1 niccli niccli niccli) ; 340 : gdphi = '*' gdphi imkvitd nivit nivit nivit ; 341 : 'SI' ('NEG' alfa 1.D0 1.D-13) ; 342 : gdphi = '*' gdphi alfa ; 343 : 'FINSI' ; 344 : dvit = dvitstar '-' gdphi ; 345 : *HCHPO dvit 'dvit' ; 346 : * 347 : 'SI' dblproj ; 348 : dphi = dphi '+' dphi1 ; 349 : 'FINSI' ; 350 : incrfin = dvit '+' dphi ; 351 : solvp = incovp '+' incrfin ; 352 : * 353 : sol = '+' solvp sol ; 354 : * 355 : 'RESPRO' sol ; 356 : * 357 : * End of procedure file GRESP 358 : * 359 : 'FINPROC' ; 360 : 361 : 362 : 363 : 364 : 365 : 366 : 367 :
© Cast3M 2003 - All rights reserved.
Disclaimer