Télécharger la plaquette Cast3M


fr en

Télécharger hbmco3.eso

Retour à la liste

Numérotation des lignes : 

hbmco3

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

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