1 : C TCRR SOURCE GOUNAND 12/12/20 21:15:24 7633 2 : SUBROUTINE TCRR 3 : C----------------------------------------------------------------------- 4 : C Test de convergence et mise à jour avec relaxation éventuelle des 5 : C CHPOINT contenu dans la table de sous type INCO. 6 : C----------------------------------------------------------------------- 7 : C 8 : C--------------------------- 9 : C Phrase d'appel (GIBIANE) : 10 : C--------------------------- 11 : C 12 : C FLOT1 = TCRR CHPO1 ( FLOT2 ) TAB1 ( 'IMPR' ENT1 ) ; 13 : C 14 : C------------------------ 15 : C Opérandes et résultat : 16 : C------------------------ 17 : C 18 : C E/ CHPO1 : CHPOINT contenant les nouvelles inconnues 19 : C E/ FLOT2 : REEL contenant le coefficient de relaxation (défaut 1.) 20 : C E/ ENT1 : ENTIER contenant la fréquence des impressions de controle 21 : C (suit le mot clé IMPR; si donné COMME SI ENT1=1) 22 : C E/S TAB1 : TABLE de sous type INCO contenant l'ensemble des champs 23 : C /S FLOT1 : REEL contenant l'erreur relative max 24 : C 25 : C---------------------------- 26 : C Indices de table modifiés : 27 : C---------------------------- 28 : C 29 : C Indices de la table INCO correspondant à une composante de CHPO1. 30 : C 31 : C-------------------------------- 32 : C Infos sur le calcul d'erreurs : 33 : C-------------------------------- 34 : C 35 : C Pour chaque composante, si CP est la nouvelle valeur à considérer 36 : C (i.e. après relaxation éventuelle) et CM l'ancienne valeur alors 37 : C on calcul MAX(CP) et ERCP=MAX(CP-CM)/MAX(CP), FLOT1 étant le maximum 38 : C des ERCP. 39 : C 40 : C----------------------------------------------------------------------- 41 : IMPLICIT INTEGER(I-N) 42 : IMPLICIT REAL*8 (A-H,O-Z) 43 : CHARACTER*8 NOMI,TYPE,TYP0,NOM 44 : CHARACTER*4 NOM4,NOMMAX 45 : LOGICAL LDUPL 46 : C 47 : -INC CCOPTIO 48 : -INC SMLENTI 49 : -INC SMTABLE 50 : -INC SMCHPOI 51 : -INC SMELEME 52 : C 53 : PARAMETER (NTB=1) 54 : DIMENSION KTAB(NTB) 55 : CHARACTER*8 LTAB(NTB) 56 : DATA LTAB /'INCO '/ 57 : C 58 : C- Initialisations 59 : C 60 : MLENTI = 0 61 : IPT1 = 0 62 : C 63 : C- Lecture de la table INCO 64 : C 65 : NTO = 1 65 : CALL LITABS(LTAB,KTAB,NTB,NTO,IRET) 66 : IF (IRET.EQ.0) RETURN 67 : C 68 : C- Lecture facultative du coefficient de relaxation 69 : C 69 : CALL LIRREE(OMEGA,0,IRET) 70 : IF (IRET.EQ.0) THEN 71 : OMEGA = 1.D0 72 : ENDIF 73 : OMEG1 = 1.D0 - OMEGA 74 : C 75 : C- Lecture facultative des fréquences d'impression de controle 76 : C 76 : CALL LIRCHA(NOM,0,IRET) 77 : IMPR = 0 78 : IF (IRET.EQ.0) THEN 79 : IMPR = 0 80 : ELSEIF (NOM.EQ.'IMPR ') THEN 80 : CALL LIRENT(IMPR,1,IRET) 81 : IF (IRET.EQ.0) RETURN 82 : ENDIF 83 : C 84 : C- Lecture du champoint contenant les nouveaux champs 85 : C 85 : CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) 86 : IF (IERR.NE.0) RETURN 87 : SEGACT MCHPOI 88 : NSOUPO = IPCHP(/1) 89 : C 90 : C- Index de la table INCO (cf opérateur INDEX) 91 : C 92 : MTAB1 = KTAB(1) 93 : SEGACT MTAB1 94 : NBINC1 = MTAB1.MLOTAB 94 : CALL ECROBJ('TABLE ',MTAB1) 94 : CALL INDETA 95 : TYPE='TABLE ' 95 : CALL LIROBJ(TYPE,MTAB2,1,IRET) 96 : SEGACT MTAB2 97 : C 98 : C-------------------------------- 99 : C= Calcul d'erreurs et relaxation 100 : C-------------------------------- 101 : C 102 : ERLM = 0.D0 103 : NOMMAX = ' ' 104 : IPTMAX = 0 105 : DO 70 K=1,NBINC1 105 : CALL ACEM(MTAB2,K,NOMI) 106 : C write(6,*)'NOMI=',NOMI 107 : TYPE=' ' 107 : CALL ACMO(MTAB1,NOMI,TYPE,MCHPO1) 108 : IF (TYPE.EQ.'CHPOINT ') THEN 109 : LDUPL=.FALSE. 110 : C 111 : C- LICHT active le MPOVAL du CHPOINT de pointeur MCHPO1 112 : C 113 : 113 : C CALL LICHT(MCHPO1,MPOVA1,TYP0,IGEOM) 114 : * gounand 07/12/2012 : dans TCRR, on change de stratégie et on crée un 115 : * chpoint tout neuf, cela permet d'éviter les appels à COPIER dans TCNM 116 : * et dans les procédures utilisateurs. On pourra aussi se référer au 117 : * pointeur pour préconditionner. 118 : SEGACT MCHPO1 119 : NSOUP1=MCHPO1.IPCHP(/1) 120 : IF(NSOUP1.EQ.0)GO TO 70 121 : DO 71 KNL=1,NSOUP1 122 : MSOUP1=MCHPO1.IPCHP(KNL) 123 : SEGACT MSOUP1 124 : MPOVA1=MSOUP1.IPOVAL 125 : IGEOM =MSOUP1.IGEOC 126 : IF(MPOVA1.EQ.0)GO TO 71 127 : * SEGACT,MPOVA1*MOD 128 : SEGACT,MPOVA1 129 : C 130 : C- KRIPAD donne le LISTENTI de correspondance 131 : C- entre les numérotations locale et globale 132 : C 133 : IF (IGEOM.NE.IPT1) THEN 134 : IF (IPT1.NE.0) SEGSUP MLENTI 135 : IPT1 = IGEOM 135 : CALL KRIPAD(IGEOM,MLENTI) 136 : 137 : ENDIF 138 : NPT = MPOVA1.VPOCHA(/1) 139 : NCI = MPOVA1.VPOCHA(/2) 140 : DO 60 I=1,NCI 141 : C 142 : C- CONVENTION sur les noms de composantes des champoints : 143 : C- CHPOINT scalaire -> indice de la table INCO 144 : C- CHPOINT vecteur -> rang de la composante + 3 premières lettres 145 : C- du nom de l'indice de la table INCO 146 : C 147 : IF (NCI.NE.1) THEN 148 : WRITE(NOM4,FMT='(I1,A3)')I,NOMI(1:3) 149 : ELSE 150 : NOM4 = NOMI(1:4) 151 : ENDIF 152 : C 153 : C- Recherche de la composante de nom NOM4. Si on la trouve, relaxation 154 : C- et calcul d'erreurs 155 : C 156 : DO 50 L=1,NSOUPO 157 : MSOUPO = IPCHP(L) 158 : SEGACT MSOUPO 159 : NC = NOCOMP(/2) 160 : DO 10 J=1,NC 161 : IF (NOCOMP(J).EQ.NOM4) GOTO 20 162 : 10 CONTINUE 163 : GOTO 40 164 : 20 CONTINUE 165 : * 166 : * Si on a trouvé une composante, on duplique MPOVA1 après le VERPAD 167 : * 168 : MELEME = IGEOC 168 : CALL VERPAD(MLENTI,MELEME,IRET) 169 : IF(IRET.NE.0)THEN 170 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 171 : MOTERR(1: 8) = 'INC '//NOM4 172 : MOTERR(9:16) = 'CHPOINT ' 172 : CALL ERREUR(788) 173 : RETURN 174 : ENDIF 175 : * duplication 176 : SEGINI,MPOVA2=MPOVA1 177 : SEGDES,MPOVA1 178 : MPOVA1=MPOVA2 179 : IF (.NOT.LDUPL) THEN 180 : *dbg WRITE(IOIMP,*) 'LDUPL NOMI=',NOMI 181 : SEGINI,MSOUP2=MSOUP1 182 : SEGDES,MSOUP1 183 : MSOUP1=MSOUP2 184 : SEGINI,MCHPO2=MCHPO1 185 : SEGDES,MCHPO1 186 : MCHPO1=MCHPO2 186 : CALL ECMO(MTAB1,NOMI,TYPE,MCHPO1) 187 : MCHPO1.IPCHP(KNL)=MSOUP1 188 : LDUPL=.TRUE. 189 : ENDIF 190 : MSOUP1.IPOVAL=MPOVA1 191 : * 192 : SEGACT MELEME 193 : MPOVAL = IPOVAL 194 : SEGACT MPOVAL 195 : NPT = VPOCHA(/1) 196 : VMAX = 0.D0 197 : ERMAX = 0.D0 198 : C write(6,*)' On relaxe NOMI NOM4 =',NOMI,NOM4 199 : DO 30 M=1,NPT 200 : MI = LECT(NUM(1,M)) 201 : IF(MI.EQ.0)GO TO 30 202 : A = OMEGA*VPOCHA(M,J) + OMEG1*MPOVA1.VPOCHA(MI 203 : $ ,I) 204 : AA = ABS(A) 205 : ER = A - MPOVA1.VPOCHA(MI,I) 206 : ER = ABS(ER) 207 : IF (ERMAX.LT.ER) THEN 208 : ERMAX = ER 209 : MER = MI 210 : ENDIF 211 : IF (AA.GT.VMAX) THEN 212 : VMAX = AA 213 : MMAX = MI 214 : ENDIF 215 : MPOVA1.VPOCHA(MI,I) = A 216 : 30 CONTINUE 217 : SEGDES MELEME 218 : SEGDES MPOVAL 219 : SEGDES MSOUPO 220 : ERL = ERMAX / (VMAX+1.D-30) 221 : IF (IMPR.GE.2) THEN 222 : WRITE(6,1000)NOM4,ERMAX,MER,ERL,VMAX,MMAX 223 : ENDIF 224 : IF (ERL.GT.ERLM) THEN 225 : ERLM = ERL 226 : NOMMAX = NOM4 227 : IPTMAX = MMAX 228 : ENDIF 229 : 40 CONTINUE 230 : 50 CONTINUE 231 : 60 CONTINUE 232 : SEGDES MPOVA1 233 : SEGDES MSOUP1 234 : 71 CONTINUE 235 : SEGDES MCHPO1 236 : ENDIF 237 : 70 CONTINUE 238 : C 239 : C- Affichage de l'erreur globale, ECRITURE et ménage 240 : C 241 : IF(IMPR.GE.1)WRITE(6,1010)NOMMAX,ERLM,IPTMAX 241 : CALL ECRREE(ERLM) 242 : IF (MLENTI.NE.0) SEGSUP MLENTI 243 : SEGDES MTAB1 244 : SEGDES MCHPOI 245 : SEGSUP MTAB2 246 : C 247 : RETURN 248 : 1000 FORMAT(' Comp : ',A4,' Err max : ',1PE8.1,' Pt ',I4,' Err % : ', 249 : & 1PE8.1,' Vmax : ',1PE8.1,' Pt ',I4) 250 : 1010 FORMAT(' Erreur relative maximale : sur la composante ', 251 : & A4,1PE8.1,' Pt ',I4) 252 : END 253 : 254 : 255 : 256 : 257 : 258 : 259 : 260 :
© Cast3M 2003 - All rights reserved.
Disclaimer