Sources Esope | Cast3M

tcrr
C TCRR      SOURCE    CB215821  20/11/25    13:40:48     10792          
      SUBROUTINE TCRR
C-----------------------------------------------------------------------
C Test de convergence et mise à jour avec relaxation éventuelle des
C CHPOINT contenu dans la table de sous type INCO.
C-----------------------------------------------------------------------
C
C---------------------------
C Phrase d'appel (GIBIANE) :
C---------------------------
C
C FLOT1 = TCRR CHPO1 ( FLOT2 ) TAB1 ( 'IMPR' ENT1 ) ;
C
C------------------------
C Opérandes et résultat :
C------------------------
C
C E/  CHPO1 : CHPOINT contenant les nouvelles inconnues
C E/  FLOT2 : REEL contenant le coefficient de relaxation (défaut 1.)
C E/  ENT1  : ENTIER contenant la fréquence des impressions de controle
C             (suit le mot clé IMPR; si donné COMME SI ENT1=1)
C E/S TAB1  : TABLE de sous type INCO contenant l'ensemble des champs
C  /S FLOT1 : REEL contenant l'erreur relative max
C
C----------------------------
C Indices de table modifiés :
C----------------------------
C
C Indices de la table INCO correspondant à une composante de CHPO1.
C
C--------------------------------
C Infos sur le calcul d'erreurs :
C--------------------------------
C
C  Pour chaque composante, si CP est la nouvelle valeur à considérer
C  (i.e. après relaxation éventuelle) et CM l'ancienne valeur alors
C  on calcul MAX(CP) et ERCP=MAX(CP-CM)/MAX(CP), FLOT1 étant le maximum
C  des ERCP.
C
C-----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
-INC PPARAM
-INC CCOPTIO
-INC SMLENTI
-INC SMTABLE
-INC SMCHPOI
-INC SMELEME
C
      PARAMETER (NTB=1)
      DIMENSION   KTAB(NTB)
      CHARACTER*8 LTAB(NTB)
      DATA LTAB /'INCO    '/      
 
      CHARACTER*8 TYPE,TYP0,NOM
      CHARACTER*(LOCOMP) NOMI,MOCOMP,NOMMAX
      LOGICAL LDUPL
C
C- Initialisations
C
      MLENTI = 0
      IPT1   = 0
C
C- Lecture de la table INCO
C
      NTO = 1
      CALL LITABS(LTAB,KTAB,NTB,NTO,IRET)
      IF (IRET.EQ.0) RETURN
C
C- Lecture facultative du coefficient de relaxation
C
      CALL LIRREE(OMEGA,0,IRET)
      IF (IRET.EQ.0) THEN
         OMEGA = 1.D0
      ENDIF
      OMEG1 = 1.D0 - OMEGA
C
C- Lecture facultative des fréquences d'impression de controle
C
      CALL LIRCHA(NOM,0,IRET)
      IMPR = 0
      IF (IRET.EQ.0) THEN
         IMPR = 0
      ELSEIF (NOM.EQ.'IMPR    ') THEN
         CALL LIRENT(IMPR,1,IRET)
         IF (IRET.EQ.0) RETURN
      ENDIF
C
C- Lecture du champoint contenant les nouveaux champs
C
      CALL LIROBJ('CHPOINT',MCHPOI,1,IRET)
      IF (IERR.NE.0) RETURN
      SEGACT MCHPOI
      NSOUPO = IPCHP(/1)
C
C- Index de la table INCO (cf opérateur INDEX)
C
      MTAB1  = KTAB(1)
      SEGACT MTAB1
      NBINC1 = MTAB1.MLOTAB
      CALL ECROBJ('TABLE   ',MTAB1)
      CALL INDETA
      TYPE='TABLE   '
      CALL LIROBJ(TYPE,MTAB2,1,IRET)
      SEGACT MTAB2
C
C--------------------------------
C= Calcul d'erreurs et relaxation
C--------------------------------
C
      ERLM   = 0.D0
      NOMMAX = '    '
      IPTMAX = 0
      DO 70 K=1,NBINC1
         CALL ACEM(MTAB2,K,NOMI)
C     write(6,*)'NOMI=',NOMI
         TYPE=' '
         CALL ACMO(MTAB1,NOMI,TYPE,MCHPO1)
         IF (TYPE.EQ.'CHPOINT ') THEN
            LDUPL=.FALSE.
C
C- LICHT active le MPOVAL du CHPOINT de pointeur MCHPO1
C
 
C        CALL LICHT(MCHPO1,MPOVA1,TYP0,IGEOM)
* gounand 07/12/2012 : dans TCRR, on change de stratégie et on crée un
* chpoint tout neuf, cela permet d'éviter les appels à COPIER dans TCNM
* et dans les procédures utilisateurs. On pourra aussi se référer au
* pointeur pour préconditionner.
            SEGACT MCHPO1
            NSOUP1=MCHPO1.IPCHP(/1)
            IF(NSOUP1.EQ.0)GO TO 70
            DO 71 KNL=1,NSOUP1
               MSOUP1=MCHPO1.IPCHP(KNL)
               SEGACT MSOUP1
               MPOVA1=MSOUP1.IPOVAL
               IGEOM =MSOUP1.IGEOC
               IF(MPOVA1.EQ.0)GO TO 71
*               SEGACT,MPOVA1*MOD
               SEGACT,MPOVA1
C
C- KRIPAD donne le LISTENTI de correspondance
C- entre les numérotations locale et globale
C
               IF (IGEOM.NE.IPT1) THEN
                  IF (IPT1.NE.0) SEGSUP MLENTI
                  IPT1 = IGEOM
                  CALL KRIPAD(IGEOM,MLENTI)
 
               ENDIF
               NPT = MPOVA1.VPOCHA(/1)
               NCI = MPOVA1.VPOCHA(/2)
               DO 60 I=1,NCI
C
C- CONVENTION sur les noms de composantes des champoints :
C-   CHPOINT scalaire -> indice de la table INCO
C-   CHPOINT vecteur  -> rang de la composante + 3 premières lettres
C-                       du nom de l'indice de la table INCO
C
                  IF (NCI.NE.1) THEN
                     WRITE(MOCOMP,FMT='(I1)')I
                     MOCOMP=MOCOMP(1:1)//NOMI(1:LOCOMP-1)
                  ELSE
                     MOCOMP=NOMI
                  ENDIF
C
C- Recherche de la composante de nom MOCOMP. Si on la trouve, relaxation
C- et calcul d'erreurs
C
                  DO 50 L=1,NSOUPO
                     MSOUPO = IPCHP(L)
                     SEGACT MSOUPO
                     NC = NOCOMP(/2)
                     DO 10 J=1,NC
                        IF (NOCOMP(J).EQ.MOCOMP) GOTO 20
 10                  CONTINUE
                     GOTO 40
 20                  CONTINUE
*
* Si on a trouvé une composante, on duplique MPOVA1 après le VERPAD
*
                     MELEME = IGEOC
                     CALL VERPAD(MLENTI,MELEME,IRET)
                     IF(IRET.NE.0)THEN
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
                        MOTERR(1: 8) = 'INC '//MOCOMP
                        MOTERR(9:16) = 'CHPOINT '
                        CALL ERREUR(788)
                        RETURN
                     ENDIF
* duplication
                     SEGINI,MPOVA2=MPOVA1
                     MPOVA1=MPOVA2
                     IF (.NOT.LDUPL) THEN
*dbg                        WRITE(IOIMP,*) 'LDUPL NOMI=',NOMI
                        SEGINI,MSOUP2=MSOUP1
                        MSOUP1=MSOUP2
                        SEGINI,MCHPO2=MCHPO1
                        MCHPO1=MCHPO2
                        CALL ECMO(MTAB1,NOMI,TYPE,MCHPO1)
                        MCHPO1.IPCHP(KNL)=MSOUP1
                        LDUPL=.TRUE.
                     ENDIF
                     MSOUP1.IPOVAL=MPOVA1
*
                     SEGACT MELEME
                     MPOVAL = IPOVAL
                     SEGACT MPOVAL
                     NPT    = VPOCHA(/1)
                     VMAX   = 0.D0
                     ERMAX  = 0.D0
C       write(6,*)' On relaxe NOMI MOCOMP =',NOMI,MOCOMP
                     DO 30 M=1,NPT
                        MI = LECT(NUM(1,M))
                        IF(MI.EQ.0)GO TO 30
                        A  = OMEGA*VPOCHA(M,J) + OMEG1*MPOVA1.VPOCHA(MI
     $                       ,I)
                        AA = ABS(A)
                        ER = A - MPOVA1.VPOCHA(MI,I)
                        ER = ABS(ER)
                        IF (ERMAX.LT.ER) THEN
                           ERMAX = ER
                           MER   = MI
                        ENDIF
                        IF (AA.GT.VMAX) THEN
                           VMAX  = AA
                           MMAX  = MI
                        ENDIF
                        MPOVA1.VPOCHA(MI,I) = A
 30                  CONTINUE
                     ERL = ERMAX / (VMAX+1.D-30)
                     IF (IMPR.GE.2) THEN
                        WRITE(6,1000)MOCOMP,ERMAX,MER,ERL,VMAX,MMAX
                     ENDIF
                     IF (ERL.GT.ERLM) THEN
                        ERLM   = ERL
                        NOMMAX = MOCOMP
                        IPTMAX = MMAX
                     ENDIF
 40                  CONTINUE
 50               CONTINUE
 60            CONTINUE
 71         CONTINUE
         ENDIF
 70   CONTINUE
C
C- Affichage de l'erreur globale, ECRITURE et ménage
C
      IF(IMPR.GE.1)WRITE(6,1010)NOMMAX,ERLM,IPTMAX
      CALL ECRREE(ERLM)
      IF (MLENTI.NE.0) SEGSUP MLENTI
      SEGDES MTAB1
      SEGSUP MTAB2
C
      RETURN
 1000 FORMAT(' Comp : ',A8,' Err max : ',1PE8.1,' Pt ',I9,' Err % : ',
     &     1PE8.1,' Vmax : ',1PE8.1,' Pt ',I9)
 1010 FORMAT(' Erreur relative maximale : sur la composante ',
     &     A4,1PE8.1,' Pt ',I9)
      END