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    OF166741  26/05/11    21:15:08     12538          
  2.  
  3. *=======================================================================
  4. *             Continuation par pseudo longueur d'arc
  5. *=======================================================================
  6.  
  7.       SUBROUTINE HBMCON(KTKAM,KTQ,KTFEX,KTPAS,KTLIAA,KTEMP,KTLIAB,
  8.      &                  KTPHI,KCPR,KOCLFA,KOCLB1,NHBM,NFFT,KPARNUM,
  9.      &                  KSORT,NSPAS,ITER)
  10.  
  11.       IMPLICIT INTEGER (I-N)
  12.       IMPLICIT REAL*8 (A-H,O-Z)
  13.  
  14. -INC PPARAM
  15. -INC CCOPTIO
  16.  
  17. -INC TMDYNC
  18.  
  19.       SEGMENT mwork
  20.         REAL*8 dw0
  21.         REAL*8 FTEST(4), FTEST0(4)
  22.         REAL*8 Rco(NT,NT),fvec(NT),XAUX(NFFT)
  23.       ENDSEGMENT
  24.  
  25.       INTEGER NSPAS
  26.  
  27.       INTEGER IREDU, METSTB, CPOSRE, IREMAX,INFO
  28.       PARAMETER(IREMAX=6)
  29.       LOGICAL CHECK,CHECKBIF, ZINIT,ZSTAB
  30.       CHARACTER*8 FLAG
  31.  
  32.       REAL*8 ZERO,ONE,TWO
  33.       PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0)
  34.  
  35. *     .. External Functions ..
  36.       REAL*8   DDOT, DNRM2
  37.       EXTERNAL DDOT, DNRM2
  38.  
  39. *   Variables generalisees
  40.       MTQ = KTQ
  41. *   Partie lineaire
  42.       MTKAM = KTKAM
  43. *   Parametres numeriques
  44.       PARNUM = KPARNUM
  45. *   Reste des segments
  46.       MTPHI  = KTPHI
  47.       MTFEX  = KTFEX
  48.       MTPAS  = KTPAS
  49.       MTLIAA = KTLIAA
  50.       MTLIAB = KTLIAB
  51.       LOCLFA = KOCLFA
  52.       LOCLB1 = KOCLB1
  53.       MTEMP  = KTEMP
  54. *   Tableau des resultats
  55.       PSORT  = KSORT
  56. *
  57.       NT = Q1(/1)
  58.       NDDL = NT/(2*NHBM+1)
  59.       PDT = 1./NFFT
  60.       PSORT.CBIF = 0
  61.  
  62.       SEGINI,MWORK
  63.  
  64. *=======================================================================
  65. *===   INITIALISATION
  66. *=======================================================================
  67.  
  68.        II=1
  69.        CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  70.        IF (BAL.EQ.1) THEN
  71.          DO I=1,NT
  72.            Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  73.          ENDDO
  74.        END IF
  75.        CALL COPYMAT(NT,RX,Rco)
  76.        CALL DGESV(NT,1,Rco,NT,IPIV,Rw,NT,INFO)
  77.        DO I=1,NT
  78.           dX(I) = -Rw(I)
  79.        ENDDO
  80. c        WRITE(*,*) '>>> dX=',(dX(IOU),IOU=1,NT)
  81.        a = dnrm2(NT,dX,1)
  82.        dw = ONE/sqrt(a*a+ONE)
  83.        IF (ISENS.LT.ZERO) THEN
  84.          dw = -dw
  85.        ENDIF
  86.        t0(NT+1)=dw
  87.        DO J=1,NT
  88.           dX(J) = dX(J)*dw
  89.           t0(J) = dX(J)
  90.        ENDDO
  91.  
  92. *      Sauvegarde des donnees de sortie
  93.        DO I=1,NT
  94.          QSAVE(I,II)=Q1(I)
  95.        ENDDO
  96.        WSAVE(II)=OMEG
  97.  
  98. *      Stabilite de la solution initiale
  99.        ZINIT = .TRUE.
  100. *      CALL STBHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,t0,ZINIT,
  101. *     &                   LSAVE(1,1,II),ZSAVE(II),FTEST,FLAG,ZSTAB)
  102.        CALL HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,ONE,ONE,ZINIT,
  103.      &                        LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB)
  104.        ZINIT = .FALSE.
  105.  
  106. *      Message relatif au pas #0 (solution initiale)
  107.       IF(IIMPI.GE.1) THEN
  108.          WRITE(IOIMP,667)
  109.          WRITE(IOIMP,660)
  110.          WRITE(IOIMP,667)
  111.          r_z=dnrm2(NT,Q1,1)
  112. c          WRITE(IOIMP,666) (II-1),ITER,OMEG,r_z,ZSTAB
  113.          WRITE(IOIMP,666) (II-1),ITER,OMEG,r_z,CPOSRE
  114.          IF (IIMPI.GE.2) WRITE(IOIMP,667)
  115.       ENDIF
  116.  
  117. *=======================================================================
  118. *      1ere prediction
  119. *=======================================================================
  120.        OMEG = OMEG + DS*dw
  121.        DO J=1,NT
  122.          Q1(J) = Q1(J)+DS*dX(J)
  123.        ENDDO
  124. c        WRITE(*,*) '>>> DS=',DS
  125. c        WRITE(*,*) '>>> Q1^pred=',(Q1(IOU),IOU=1,NT)
  126.        IREDU = 0
  127.        II = 2
  128.        METSTB = 0
  129.  
  130. *=======================================================================
  131. *======================= Boucle sur la frequence =======================
  132. *=======================================================================
  133.        DO WHILE (II .LE. NBPAS)
  134.  
  135. *       === Correction =================================================
  136.         CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
  137.      &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'CF',ITER)
  138.  
  139. *       -----------------------------
  140. *       -----------------------------
  141. *       ---- si NEWT a converge, ----
  142. *       -----------------------------
  143. *       -----------------------------
  144.         IF (.NOT.CHECK) THEN
  145.  
  146. *          Sauvegarde des donnees de sortie
  147. *          Frequence, coeffs de Fourier, Norme des coeffs de Fourier
  148.            DO I=1,NT
  149.              QSAVE(I,II)=Q1(I)
  150.            ENDDO
  151.            WSAVE(II)=OMEG
  152. *
  153. *          === Prediction ==============================================
  154. *          Pas tangent a la courbe
  155.            CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  156.            IF (BAL.EQ.1) THEN
  157.              DO I=1,NT
  158.                Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  159.              ENDDO
  160.            END IF
  161.            DO KK= 1,NT
  162.              Rw2(KK) = Rw(KK)
  163.            ENDDO
  164.            CALL COPYMAT(NT,RX,Rco)
  165.            CALL DGESV(NT,1,Rco,NT,IPIV,Rw,NT,INFO)
  166.            DO J=1,NT
  167.               dX(J) = -Rw(J)
  168.               tp(J) = dX(J)
  169.            ENDDO
  170.            tp(NT+1) = ONE
  171.            a=dnrm2(NT,dX,1)
  172.            dw0 = dw
  173.            dw = ONE/sqrt(ONE+a**2)
  174. c             WRITE(*,*) 'dw=',dw
  175. c          signe de dw = | fourni par l'utilisateur si 1er pas
  176. c                        | tel que l'on continu colineairement au pas precedent
  177.            IF (II.EQ.2) THEN
  178.              IF (ISENS.LT.ZERO) THEN
  179.                dw = -dw
  180.              ENDIF
  181.            ELSE
  182.              aux = DDOT(NT+1,t0,1,tp,1)
  183.              dw = (SIGN(dw,aux))
  184.            ENDIF
  185.            DO K = 1,NT
  186.               dX(K) = dX(K)*dw
  187.               tp(K) = dX(K)
  188.            ENDDO
  189.            tp(NT+1) = dw
  190.  
  191. *          === Stabilite pour la solution convergee ====================
  192. *          via Methode de Hill
  193.            IF (METSTB.EQ.0) THEN
  194. *             FTEST0 = FTEST
  195. *             CALL STBHILL(NT,NHBM,NDDL,OMEG,Q1,Rco,XM,XASM,t0,ZINIT,
  196. *     &                    LSAVE(1,1,II),ZSAVE(II),FTEST,FLAG,ZSTAB)
  197.              DO J=1,NT
  198.              DO I=1,NT
  199.                  RX(I,J) = ZZ(I,J)-JAC(I,J)
  200.              ENDDO
  201.              ENDDO
  202.              CALL HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,dw0,dw,ZINIT,
  203.      &                        LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB)
  204. *          via Matrice de Monodromie
  205.            ELSE
  206. *             CALL STBMONO(NT,NDDL,NFFT,NHBM,MTQ,MTKAM,MTPHI,MTLIAA,
  207. *    &          MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,FLAG)
  208.            ENDIF
  209. *          Sauvegarde des donnees de sortie: LSAVE et ZSTAB rempli par STBHILL
  210.  
  211. *          WRITE(*,*) 'Pas ',II,', w=',OMEG,'STAB =',ZSTAB
  212.  
  213. *          === Message relatif au pas #II ==============================
  214.            IF(IIMPI.GE.1) THEN
  215.               r_z=dnrm2(NT,Q1,1)
  216. c               WRITE(IOIMP,666) (II-1),ITER,OMEG,r_z,ZSTAB
  217.               WRITE(IOIMP,666) (II-1),ITER,OMEG,r_z,CPOSRE
  218.               IF (IIMPI.GE.2) WRITE(IOIMP,667)
  219.            ENDIF
  220. *
  221. *          === Bifurcation? ============================================
  222.            IF ((FLAG.NE.'S')) THEN
  223.              WRITE(IOIMP,*) 'Bifurcation detectee! Type :',FLAG
  224. cbp            IF (FLAG.NE.'N') THEN
  225.                CALL HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,
  226.      &         MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECKBIF,FLAG,
  227.      &         PARNUM,PSORT)
  228.                IF (.NOT.CHECKBIF) THEN
  229.                    WRITE(IOIMP,*) 'Bifurcation localisee :)'
  230.                ELSE
  231.                    WRITE(IOIMP,*) 'Bifurcation non localisee :('
  232.                ENDIF
  233. cbp            ENDIF
  234.            ENDIF
  235.  
  236. *          === Tests d'arret ===========================================
  237. *          On verifie la condition d'arret par rapport a OMEG
  238. *           IF ((OMEG.GE.PARFIN).OR.(OMEG.LE.PARINI)) THEN
  239.            IF ((OMEG-PARFIN)*(PARFIN-PARINI).GE.0) THEN
  240.                WRITE(IOIMP,*) 'Fin de la continuation apres',II,' pas.'
  241.                IF (II.LT.NBPAS) THEN
  242. *                  WRITE(*,*) 'NPAS sera egal a:',II
  243.                   NT1 = QSAVE(/1)
  244.                   NA1 = LSAVE(/2)/2
  245. *                  NPAS = QSAVE(/2)
  246.                   NBIFU = QBIFU(/2)
  247. *                  write(*,*) NT1,NA1,NPAS,NBIFU
  248.                   NPAS = II
  249.                   SEGADJ, PSORT
  250.                ENDIF
  251.                KSORT = PSORT
  252.                RETURN
  253.            ENDIF
  254.  
  255. *          === prediction + donnees pour le prochain pas ===============
  256. *          ajustement automatique de la longueur du pas
  257. c            if(II.ge.62) write(*,*) '>>>t0=',(t0(iou),iou=1,NT+1)
  258. c            if(II.ge.62) write(*,*) '>>>tp=',(tp(iou),iou=1,NT+1)
  259.            CALL ADPAS(DS,DSMIN,DSMAX,ITER,ITERMOY,ANGMIN,ANGMAX,
  260.      &                NT+1,t0,tp)
  261. c            WRITE(*,*) '>>> DS=',DS
  262. *          stockage des donnees de fin de pas utiles (dans *OLD)
  263. *          pas predictif : Q1, OMEG
  264.            DO I=1,NT+1
  265.               t0(I) = tp(I)
  266.            ENDDO
  267. c            DSREDU=MAX(0.25*DS,DSMIN)
  268.            DO KK = 1,NT
  269. c c             ici, on anticipe la reduction du pas predicteur
  270. c c             en maintenant la direction dX fixe
  271. c               QOLD(KK) = Q1(KK)+DSREDU*dX(KK)
  272. c             et ici, on laisse la possibilite a DX d'evoluer
  273.               QOLD(KK) = Q1(KK)
  274.               Q1(KK)   = Q1(KK)+DS*dX(KK)
  275.            ENDDO
  276. c          idem pour OMEGOLD que QOLD
  277. c            OMEGOLD = OMEG + DSREDU*dw
  278.            OMEGOLD = OMEG
  279.            OMEG    = OMEG + DS*dw
  280.            IREDU = 0
  281.            ISTRA2=0
  282. cdebug
  283. c            if(II.ge.62.and.II.le.68) then
  284. c            write(*,*) 'prediction Q=',(Q1(iou),iou=1,NT)
  285. c            write(*,*) 'DS, OMEG=',DS,OMEG
  286. c            endif
  287.  
  288. *       -----------------------------------
  289. *       -----------------------------------
  290. *       ---- si NEWT n'a pas converge, ----
  291. *       -----------------------------------
  292. *       -----------------------------------
  293.         ELSE
  294.  
  295. *          === Message relatif au pas non converge #II =================
  296.            IF(IIMPI.GE.1) THEN
  297.               r_z=dnrm2(NT,Q1,1)
  298.               WRITE(IOIMP,665) (II-1),ITER,OMEG,r_z
  299.               IF (IIMPI.GE.2) WRITE(IOIMP,667)
  300.            ENDIF
  301.  
  302. *          === strategie 2 (non-convergence alors qu'on a deja =========
  303. *          ===              atteint DSMIN) : on quitte !       =========
  304.            IF (DS.LE.DSMIN) THEN
  305.              DS=DSMIN
  306.              DO I=1,NT
  307.                Q1(I) = QOLD(I)
  308.              ENDDO
  309.              OMEG = OMEGOLD
  310. c            et on arrete
  311.              IREDU=IREMAX
  312.              II = II-1
  313.  
  314. *          === strategie 1 : Recommencer le pas en diminuant DS ========
  315. c          test sur II>1 car pour l'instant on n'a pas branche la recup d'une solution initiale
  316.            ELSEIF (II.GT.1) THEN
  317. c c            en maintenant la direction dX fixe
  318. c c            (attention : DS doit etre = DSREDU !!!)
  319. c              DS = MAX(0.25*DS,DSMIN)
  320. c              DO I=1,NT
  321. c                Q1(I) = QOLD(I)
  322. c              ENDDO
  323. c              OMEG = OMEGOLD
  324. *            mise a jour de la prediction avec nouveau DS (dX et dw n'ont pas change)
  325.              DS = MAX(0.25D0*DS,DSMIN)
  326.              DO I=1,NT
  327.                Q1(I) = QOLD(I)+DS*dX(I)
  328.              ENDDO
  329.              OMEG = OMEGOLD + DS*dw
  330.              II = II-1
  331.              IREDU = IREDU + 1
  332. *            message
  333.              IF (IIMPI.GE.2) WRITE(IOIMP,668) IREDU,DS
  334.  
  335.            ENDIF
  336.  
  337.  
  338.  888       continue
  339. c          apres reduction de 0.25**10 = 1.E-6 , echec !
  340. c                             0.25**5  = 1.E-3
  341.            IF (IREDU.GE.IREMAX) THEN
  342. *             message interne
  343.               IF (IIMPI.GE.1) WRITE(IOIMP,669) DS
  344. c             Pas de convergence apres %i1 iterations. L'execution continue
  345.               INTERR(1)=ITERMAX
  346.               CALL ERREUR(151)
  347. *             ajustement des tableaux de sorties
  348.               IF (II.LT.NBPAS) THEN
  349.                   NT1 = QSAVE(/1)
  350.                   NA1 = LSAVE(/2)/2
  351. *                  NPAS = QSAVE(/2)
  352.                   NBIFU = QBIFU(/2)
  353.                   NPAS = II
  354.                   SEGADJ, PSORT
  355.               ENDIF
  356.               KSORT = PSORT
  357.               RETURN
  358.            ENDIF
  359.  
  360.         ENDIF
  361. *       ----------------------------------------------
  362. *       ----------------------------------------------
  363. *       ---- fin distinction NEWT converge ou pas ----
  364. *       ----------------------------------------------
  365. *       ----------------------------------------------
  366.  
  367. *        incrementation du compteur
  368.          II = II+1
  369.  
  370.       ENDDO
  371. *=======================================================================
  372. *       fin de la boucle sur la frequence
  373. *=======================================================================
  374.  
  375. *   message
  376. c      IF(IIMPI.GE.1) WRITE(IOIMP,667)
  377. c      WRITE(IOIMP,*) NBPAS, 'pas de continuation realises!'
  378.  
  379.       KSORT = PSORT
  380.  
  381.       SEGSUP,MWORK
  382.  
  383. *=======================================================================
  384. *      Mise en forme des messages
  385. *=======================================================================
  386. *   Dernier message pour fermer le tableau dans le cas iimpi=1
  387.       IF (IIMPI.EQ.1) WRITE(IOIMP,667)
  388.  
  389. c  660   FORMAT(' |  Pas  | Iter |     w     |     Q     | Stab |')
  390.  660   FORMAT(' |  Pas  | Iter |     w     |     Q     | Unst |')
  391.  665   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' |   -  |')
  392. c  666   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ',L3,'  |')
  393.  666   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ',I3,'  |')
  394.  667   FORMAT(' +-------+------+-----------+-----------+------+')
  395.  668   FORMAT(' + cont : ds=ds_initial/4**',I2,'=',F13.9)
  396.  669   FORMAT(' + cont : pas de convergence avec ds=dsmin=',F13.9)
  397.  
  398.       RETURN
  399.       END
  400.  
  401.  
  402.  

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