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

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