Télécharger la plaquette Cast3M


fr en

Télécharger hbmcon.eso

Retour à la liste

Numérotation des lignes : 

hbmcon

  1. C HBMCON    SOURCE    CB215821  23/01/25    21:15:20     11573          
  2.  
  3.        SUBROUTINE HBMCON(KTKAM,KTQ,KTFEX,KTPAS,KTLIAA,KTEMP,KTLIAB,KTPHI
  4.      &            ,KCPR,KOCLFA,KOCLB1,NHBM,NFFT,KPARNUM,KSORT,SPAS,ITER) 
  5.  
  6. *=======================================================================
  7. *             Continuation par pseudo longueur d'arc
  8. *======================================================================= 
  9.  
  10.        IMPLICIT INTEGER (I-N)
  11.        IMPLICIT REAL*8 (A-H,O-Z)
  12.  
  13.        INTEGER SPAS, IREDU, METSTB, CPOSRE, CBIF,IREMAX,INFO
  14.        LOGICAL CHECK,CHECKBIF
  15.        LOGICAL ZINIT,ZSTAB
  16.        CHARACTER*8 FLAG
  17.        REAL*8 FTEST(4), FTEST0(4)
  18.        REAL*8 XAUX(NFFT),nrmR,dw0,Q1NRM2
  19.        PARAMETER(IREMAX=6)
  20.        segment mwork
  21.          REAL*8 Rco(NT,NT)
  22.          REAL*8 fvec(NT)         
  23.        ENDSEGMENT
  24. *     
  25. -INC PPARAM    
  26. -INC CCOPTIO    
  27. *-INC TMDYNC.INC
  28. ************************** debut TMDYNC.INC ****************************
  29.  
  30. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  31. *    TODO   : a extraire dans un include des que stabilise
  32. *
  33. *     Segment des variables generalisees:
  34. *     -----------------------------------
  35.       SEGMENT MTQ
  36.          REAL*8 Q1(NT1)
  37.          REAL*8 OMEG,XPARA
  38.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  39.          REAL*8 dX(NT1), dw, dv
  40.       ENDSEGMENT
  41. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  42. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  43. *            avec q_i vecteur de dimension n ou n=nombre de modes
  44. *     OMEG : frequence fondamentale de l'approximation
  45. *     XPARA: parametre de continuation (par defaut la frequence) 
  46. *            \in [PARINI,PARFIN]
  47. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  48. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  49. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  50. *            lineaires et constantes
  51. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  52. *
  53. *
  54. *     Segment contenant les matrices XK, XASM et XM:
  55. *     ---------------------------------------------
  56.       SEGMENT MTKAM
  57.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  58.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  59. *         REAL*8 GAMFIN(NPC2,nl1)
  60.       ENDSEGMENT
  61. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  62. *     GAM et IGAM : matrices pour la FFT et son inverse
  63. *     GAMFIN      : 
  64. *
  65. *     Segment des deformees modales:
  66. *     ------------------------------
  67. *     (idem DYNE)
  68.       SEGMENT MTPHI
  69.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  70.          INTEGER IAROTA(NSB)
  71.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  72.       ENDSEGMENT
  73. *
  74. *     Segment descriptif des liaisons en base A:
  75. *     ------------------------------------------
  76. *     (idem DYNE)
  77.       SEGMENT MTLIAA
  78.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  79.          REAL*8  XPALA(NLIAA,NXPALA)
  80.       ENDSEGMENT
  81. *
  82. *     Segment descriptif des liaisons en base B:
  83. *     ------------------------------------------
  84. *     (idem DYNE)
  85.       SEGMENT MTLIAB
  86.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  87.          REAL*8  XPALB(NLIAB,NXPALB)
  88.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  89.       ENDSEGMENT
  90. *
  91. *     Segment representant les chargements exterieurs:
  92. *     -----------------------------------------------
  93.       SEGMENT MTFEX
  94.          REAL*8  FEXA(NT1)
  95.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  96.          INTEGER BAL
  97.       ENDSEGMENT
  98. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  99. *            Fourier et exprimes en base A
  100. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  101. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  102. *             compatibilite avec calcul des Fnl.
  103. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  104. *            (cad proportionnel a OMEG**2)
  105. *
  106. *     Segment "local" pour DEVLFA:
  107. *     ----------------------------
  108.       SEGMENT LOCLFA
  109.          REAL*8 FTEST(NA1,4)
  110.       ENDSEGMENT
  111. *      
  112. *     Segment "local" pour DEVLB1:
  113. *     ----------------------------
  114.       SEGMENT LOCLB1
  115.          REAL*8 FTEST2(NPLB,6)
  116.       ENDSEGMENT
  117. *      
  118. *     Segment contenant les variables au cours d un pas de temps:
  119. *     ----------------------------------------------------------
  120.       SEGMENT MTPAS
  121.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  122.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  123.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  124.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  125.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  126.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  127.       ENDSEGMENT
  128. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  129. *     XPTB       : deplacement du point d'une liaison en base B
  130. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  131. *     FEXB       : forces exterieures en base B (a priori uniquement 
  132. *                  pour les moments appliques aux rotations rigides ?)
  133. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  134. *                  maillage de contact dans certaines liaisons)
  135. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  136. *                       base A/B (= contributions a dFnl/dX)
  137. *
  138. *
  139. *     Segment des points de reference des modes (base A):
  140. *     --------------------------------------------------
  141.       SEGMENT MPREF
  142.          INTEGER IPOREF(NPREF)
  143.       ENDSEGMENT
  144. *
  145. *     Segment des points en base B:
  146. *     -----------------------------
  147.       SEGMENT NCPR(NBPTS)
  148. *     NCRP(#global) = #local dans XPTB (1er indice)      
  149. *
  150. *     Segment des parametres numeriques pour la continuation: 
  151. *     ------------------------------------------------------
  152.       SEGMENT PARNUM
  153.          CHARACTER*4 TYPS
  154.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  155.          REAL*8   PARINI,PARFIN
  156.          INTEGER  ITERMAX,NBPAS
  157.          LOGICAL  JANAL
  158.       ENDSEGMENT
  159. *
  160. *     Segment des resultats:
  161. *     ---------------------
  162.       SEGMENT  PSORT
  163.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  164.         REAL*8 VSAVE(NPAS)
  165.         LOGICAL ZSAVE(NPAS)
  166.         CHARACTER*2 TYPBIF(NBIFU)
  167.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  168.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  169.         INTEGER CBIF
  170.       ENDSEGMENT
  171. *     QSAVE(i,j) = Q harmonique i au pas j
  172. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  173. *     ZSAVE(j)   = stabilite au j-eme pas
  174. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  175. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  176. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  177. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  178. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  179. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  180.  
  181. *     Segment des tableaux de travail:
  182. *     -------------------------------
  183.       SEGMENT  MTEMP
  184.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  185.          REAL*8 T02(NT1+2), TP2(NT1+2)
  186.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  187.          REAL*8 res
  188.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  189.          REAL*8 QOLD(NT1),OMEGOLD
  190.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  191.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  192.       ENDSEGMENT
  193. *     Jacobiennes augmentees 
  194. *     Ja : [ RX  Rw ; dX dw]
  195. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  196.  
  197. *      SEGMENT  NNNN
  198. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  199. *      ENDSEGMENT
  200.  
  201. *************************** fin TMDYNC.INC *****************************
  202.  
  203.        REAL*8 ZERO,ONE,TWO
  204.        PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0)
  205.  
  206. *     .. External Functions ..
  207.        REAL*8   DDOT, DNRM2
  208.        EXTERNAL DDOT, DNRM2
  209.  
  210. *      Variables generalisees
  211.        MTQ = KTQ 
  212. *      Partie lineaire
  213.        MTKAM = KTKAM
  214. *      Parametres numeriques
  215.        PARNUM = KPARNUM 
  216. *      Reste des segments 
  217.        MTPHI  = KTPHI
  218.        MTFEX  = KTFEX
  219.        MTPAS  = KTPAS
  220.        MTLIAA = KTLIAA
  221.        MTLIAB = KTLIAB
  222.        LOCLFA = KOCLFA
  223.        LOCLB1 = KOCLB1     
  224.        MTEMP  = KTEMP
  225. *      Tableau des resultats
  226.        PSORT  = KSORT    
  227. *
  228.        NT = Q1(/1)
  229.        NDDL = NT/(2*NHBM+1)
  230.        PDT = 1./NFFT
  231.        CBIF = 0
  232.        segini,MWORK
  233.  
  234. *======================================================================= 
  235. *===   INITIALISATION 
  236. *======================================================================= 
  237.  
  238.        II=1
  239.        CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  240.        IF (BAL.EQ.1) THEN
  241.          DO I=1,NT
  242.            Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  243.          ENDDO
  244.        END IF
  245.        CALL COPYMAT(NT,RX,Rco)
  246.        CALL DGESV(NT,1,Rco,NT,IPIV,Rw,NT,INFO)
  247.        DO I=1,NT
  248.           dX(I) = -Rw(I)
  249.        ENDDO 
  250. c        WRITE(*,*) '>>> dX=',(dX(IOU),IOU=1,NT)
  251.        a = dnrm2(NT,dX,1)
  252.        dw = ONE/sqrt(a**2+ONE)
  253.        IF (ISENS.LT.ZERO) THEN
  254.          dw = -dw
  255.        ENDIF     
  256.        t0(NT+1)=dw
  257.        DO J=1,NT
  258.           dX(J) = dX(J)*dw
  259.           t0(J) = dX(J)
  260.        ENDDO
  261.  
  262. *      Sauvegarde des donnees de sortie
  263.        DO I=1,NT
  264.          QSAVE(I,II)=Q1(I)
  265.        ENDDO
  266.        WSAVE(II)=OMEG  
  267.  
  268. *      Stabilite de la solution initiale
  269.        ZINIT = .TRUE.
  270. *      CALL STBHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,t0,ZINIT,
  271. *     &                   LSAVE(1,1,II),ZSAVE(II),FTEST,FLAG,ZSTAB)
  272.        CALL HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,ONE,ONE,ZINIT,
  273.      &                        LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB)
  274.        ZINIT = .FALSE.
  275.  
  276. *      Message relatif au pas #0 (solution initiale)
  277.        IF(IIMPI.GE.1) THEN
  278.          WRITE(IOIMP,667) 
  279.          WRITE(IOIMP,660)
  280.          WRITE(IOIMP,667) 
  281.          Q1NRM2=dnrm2(NT,Q1,1)
  282. c          WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2,ZSTAB
  283.          WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2,CPOSRE
  284.          IF (IIMPI.GE.2) WRITE(IOIMP,667) 
  285.        ENDIF
  286.  
  287.  
  288. *=======================================================================  
  289. *      1ere prediction
  290. *=======================================================================  
  291.        OMEG = OMEG + DS*dw
  292.        DO J=1,NT
  293.          Q1(J) = Q1(J)+DS*dX(J)
  294.        ENDDO        
  295. c        WRITE(*,*) '>>> DS=',DS
  296. c        WRITE(*,*) '>>> Q1^pred=',(Q1(IOU),IOU=1,NT)
  297.        IREDU = 0
  298.        II = 2
  299.        METSTB = 0
  300.  
  301. *=======================================================================  
  302. *======================= Boucle sur la frequence =======================        
  303. *=======================================================================  
  304.        DO WHILE (II .LE. NBPAS)
  305.  
  306. *       === Correction =================================================
  307.         CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
  308.      &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'CF',ITER)   
  309.  
  310. *       -----------------------------
  311. *       -----------------------------
  312. *       ---- si NEWT a converge, ----
  313. *       -----------------------------
  314. *       -----------------------------
  315.         IF (.NOT.CHECK) THEN
  316.  
  317. *          Sauvegarde des donnees de sortie
  318. *          Frequence, coeffs de Fourier, Norme des coeffs de Fourier           
  319.            DO I=1,NT
  320.              QSAVE(I,II)=Q1(I)
  321.            ENDDO
  322.            WSAVE(II)=OMEG           
  323. * 
  324. *          === Prediction ==============================================
  325. *          Pas tangent a la courbe 
  326.            CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  327.            IF (BAL.EQ.1) THEN
  328.              DO I=1,NT
  329.                Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  330.              ENDDO
  331.            END IF
  332.            DO KK= 1,NT
  333.              Rw2(KK) = Rw(KK)
  334.            ENDDO
  335.            CALL COPYMAT(NT,RX,Rco)
  336.            CALL DGESV(NT,1,Rco,NT,IPIV,Rw,NT,INFO)
  337.            DO J=1,NT
  338.               dX(J) = -Rw(J)
  339.               tp(J) = dX(J)
  340.            ENDDO
  341.            tp(NT+1) = ONE
  342.            a=dnrm2(NT,dX,1)
  343.            dw0 = dw
  344.            dw = ONE/sqrt(ONE+a**2)
  345. c             WRITE(*,*) 'dw=',dw
  346. c          signe de dw = | fourni par l'utilisateur si 1er pas
  347. c                        | tel que l'on continu colineairement au pas precedent
  348.            IF (II.EQ.2) THEN
  349.              IF (ISENS.LT.ZERO) THEN
  350.                dw = -dw
  351.              ENDIF
  352.            ELSE 
  353.              aux = DDOT(NT+1,t0,1,tp,1)
  354.              dw = (SIGN(dw,aux))
  355.            ENDIF
  356.            DO K = 1,NT
  357.               dX(K) = dX(K)*dw
  358.               tp(K) = dX(K)
  359.            ENDDO 
  360.            tp(NT+1) = dw
  361.  
  362. *          === Stabilite pour la solution convergee ====================
  363. *          via Methode de Hill
  364.            IF (METSTB.EQ.0) THEN
  365. *             FTEST0 = FTEST
  366. *             CALL STBHILL(NT,NHBM,NDDL,OMEG,Q1,Rco,XM,XASM,t0,ZINIT,
  367. *     &                    LSAVE(1,1,II),ZSAVE(II),FTEST,FLAG,ZSTAB)
  368.              DO J=1,NT
  369.              DO I=1,NT
  370.                  RX(I,J) = ZZ(I,J)-JAC(I,J)
  371.              ENDDO           
  372.              ENDDO
  373.              CALL HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,dw0,dw,ZINIT,
  374.      &                        LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB)
  375. *          via Matrice de Monodromie    
  376.            ELSE 
  377. *             CALL STBMONO(NT,NDDL,NFFT,NHBM,MTQ,MTKAM,MTPHI,MTLIAA,                   
  378. *    &          MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,FLAG)
  379.            ENDIF           
  380. *          Sauvegarde des donnees de sortie: LSAVE et ZSTAB rempli par STBHILL
  381.  
  382. *          WRITE(*,*) 'Pas ',II,', w=',OMEG,'STAB =',ZSTAB
  383.  
  384. *          === Message relatif au pas #II ==============================
  385.            IF(IIMPI.GE.1) THEN
  386.               Q1NRM2=dnrm2(NT,Q1,1)
  387. c               WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2,ZSTAB
  388.               WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2,CPOSRE
  389.               IF (IIMPI.GE.2) WRITE(IOIMP,667) 
  390.            ENDIF
  391. *  
  392. *          === Bifurcation? ============================================
  393.            IF ((FLAG.NE.'S')) THEN 
  394.              WRITE(IOIMP,*) 'Bifurcation detectee! Type :',FLAG
  395. cbp            IF (FLAG.NE.'N') THEN
  396.                CALL HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,
  397.      &         MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECKBIF,FLAG,
  398.      &         PARNUM,PSORT)
  399.                IF (.NOT.CHECKBIF) THEN
  400.                    WRITE(IOIMP,*) 'Bifurcation localisee :)'
  401.                ELSE
  402.                    WRITE(IOIMP,*) 'Bifurcation non localisee :('
  403.                ENDIF
  404. cbp            ENDIF
  405.            ENDIF
  406.  
  407. *          === Tests d'arret ===========================================        
  408. *          On verifie la condition d'arret par rapport a OMEG
  409. *           IF ((OMEG.GE.PARFIN).OR.(OMEG.LE.PARINI)) THEN
  410.            IF ((OMEG-PARFIN)*(PARFIN-PARINI).GE.0) THEN
  411.                WRITE(IOIMP,*) 'Fin de la continuation apres',II,' pas.'
  412.                IF (II.LT.NBPAS) THEN
  413. *                  WRITE(*,*) 'NPAS sera egal a:',II
  414.                   NT1 = QSAVE(/1)
  415.                   NA1 = LSAVE(/2)/2
  416. *                  NPAS = QSAVE(/2)
  417.                   NBIFU = QBIFU(/2)
  418. *                  write(*,*) NT1,NA1,NPAS,NBIFU
  419.                   NPAS = II
  420.                   SEGADJ, PSORT
  421.                ENDIF
  422.                KSORT = PSORT
  423.                RETURN
  424.            ENDIF
  425.  
  426. *          === prediction + donnees pour le prochain pas ===============
  427. *          ajustement automatique de la longueur du pas
  428. c            if(II.ge.62) write(*,*) '>>>t0=',(t0(iou),iou=1,NT+1)
  429. c            if(II.ge.62) write(*,*) '>>>tp=',(tp(iou),iou=1,NT+1)
  430.            CALL ADPAS(DS,DSMIN,DSMAX,ITER,ITERMOY,ANGMIN,ANGMAX,
  431.      &                NT+1,t0,tp)
  432. c            WRITE(*,*) '>>> DS=',DS
  433. *          stockage des donnees de fin de pas utiles (dans *OLD)
  434. *          pas predictif : Q1, OMEG
  435.            DO I=1,NT+1
  436.               t0(I) = tp(I)
  437.            ENDDO
  438. c            DSREDU=MAX(0.25*DS,DSMIN)
  439.            DO KK = 1,NT
  440. c c             ici, on anticipe la reduction du pas predicteur 
  441. c c             en maintenant la direction dX fixe
  442. c               QOLD(KK) = Q1(KK)+DSREDU*dX(KK)
  443. c             et ici, on laisse la possibilite a DX d'evoluer
  444.               QOLD(KK) = Q1(KK)
  445.               Q1(KK)   = Q1(KK)+DS*dX(KK)
  446.            ENDDO
  447. c          idem pour OMEGOLD que QOLD
  448. c            OMEGOLD = OMEG + DSREDU*dw
  449.            OMEGOLD = OMEG
  450.            OMEG    = OMEG + DS*dw
  451.            IREDU = 0
  452.            ISTRA2=0
  453. cdebug           
  454. c            if(II.ge.62.and.II.le.68) then
  455. c            write(*,*) 'prediction Q=',(Q1(iou),iou=1,NT)
  456. c            write(*,*) 'DS, OMEG=',DS,OMEG
  457. c            endif
  458.  
  459. *       -----------------------------------
  460. *       -----------------------------------
  461. *       ---- si NEWT n'a pas converge, ----
  462. *       -----------------------------------
  463. *       -----------------------------------
  464.         ELSE
  465.  
  466. *          === Message relatif au pas non converge #II =================
  467.            IF(IIMPI.GE.1) THEN
  468.               Q1NRM2=dnrm2(NT,Q1,1)
  469.               WRITE(IOIMP,665) (II-1),ITER,OMEG,Q1NRM2
  470.               IF (IIMPI.GE.2) WRITE(IOIMP,667) 
  471.            ENDIF
  472.  
  473. *          === strategie 2 (non-convergence alors qu'on a deja ========= 
  474. *          ===              atteint DSMIN) : on quitte !       =========
  475.            IF (DS.LE.DSMIN) THEN
  476.              DS=DSMIN
  477.              DO I=1,NT
  478.                Q1(I) = QOLD(I)
  479.              ENDDO
  480.              OMEG = OMEGOLD           
  481. c            et on arrete
  482.              IREDU=IREMAX
  483.              II = II-1
  484.  
  485. *          === strategie 1 : Recommencer le pas en diminuant DS ======== 
  486. c          test sur II>1 car pour l'instant on n'a pas branche la recup d'une solution initiale
  487.            ELSEIF (II.GT.1) THEN
  488. c c            en maintenant la direction dX fixe 
  489. c c            (attention : DS doit etre = DSREDU !!!)
  490. c              DS = MAX(0.25*DS,DSMIN)
  491. c              DO I=1,NT
  492. c                Q1(I) = QOLD(I)
  493. c              ENDDO
  494. c              OMEG = OMEGOLD
  495. *            mise a jour de la prediction avec nouveau DS (dX et dw n'ont pas change)
  496.              DS = MAX(0.25D0*DS,DSMIN)
  497.              DO I=1,NT
  498.                Q1(I) = QOLD(I)+DS*dX(I)
  499.              ENDDO
  500.              OMEG = OMEGOLD + DS*dw
  501.              II = II-1
  502.              IREDU = IREDU + 1
  503. *            message
  504.              IF (IIMPI.GE.2) WRITE(IOIMP,668) IREDU,DS
  505.  
  506.            ENDIF
  507.  
  508.  
  509.  888       continue                  
  510. c          apres reduction de 0.25**10 = 1.E-6 , echec !
  511. c                             0.25**5  = 1.E-3
  512.            IF (IREDU.GE.IREMAX) THEN        
  513. *             message interne
  514.               IF (IIMPI.GE.1) WRITE(IOIMP,669) DS
  515. c             Pas de convergence apres %i1 iterations. L'execution continue
  516.               INTERR(1)=ITERMAX              
  517.               CALL ERREUR(151)
  518. *             ajustement des tableaux de sorties
  519.               IF (II.LT.NBPAS) THEN
  520.                   NT1 = QSAVE(/1)
  521.                   NA1 = LSAVE(/2)/2
  522. *                  NPAS = QSAVE(/2)
  523.                   NBIFU = QBIFU(/2)
  524.                   NPAS = II
  525.                   SEGADJ, PSORT
  526.               ENDIF
  527.               KSORT = PSORT
  528.               RETURN
  529.            ENDIF
  530.  
  531.  
  532.  
  533.         ENDIF 
  534. *       ----------------------------------------------
  535. *       ----------------------------------------------
  536. *       ---- fin distinction NEWT converge ou pas ----  
  537. *       ----------------------------------------------
  538. *       ----------------------------------------------
  539.  
  540. *        incrementation du compteur
  541.          II = II+1
  542.  
  543.         ENDDO
  544. *=======================================================================  
  545. *       fin de la boucle sur la frequence         
  546. *=======================================================================  
  547.  
  548.  
  549. *       message
  550. c         IF(IIMPI.GE.1) WRITE(IOIMP,667) 
  551. c       WRITE(IOIMP,*) NBPAS, 'pas de continuation realises!' 
  552.         KSORT = PSORT
  553.  
  554.         segsup,MWORK
  555.  
  556. *=======================================================================
  557. *      Mise en forme des messages
  558. *=======================================================================  
  559. *      dernier message pour fermer le tableau dans le cas iimpi=1
  560.        IF (IIMPI.EQ.1) WRITE(IOIMP,667) 
  561.  
  562. c  660   FORMAT(' |  Pas  | Iter |     w     |     Q     | Stab |')
  563.  660   FORMAT(' |  Pas  | Iter |     w     |     Q     | Unst |')
  564.  665   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' |   -  |')
  565. c  666   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ',L3,'  |')
  566.  666   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ',I3,'  |')
  567.  667   FORMAT(' +-------+------+-----------+-----------+------+')
  568.  668   FORMAT(' + cont : ds=ds_initial/4**',I2,'=',F13.9)
  569.  669   FORMAT(' + cont : pas de convergence avec ds=dsmin=',F13.9)
  570.  
  571.       END
  572.  
  573.  
  574.  
  575.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales