C HBMCO3 SOURCE OF166741 26/05/11 21:15:08 12538 SUBROUTINE HBMCO3(KTKAM,KTQ,KTFEX,KTPAS,KTLIAA,KTEMP,KTLIAB, & KTPHI,KCPR,KOCLFA,KOCLB1,NHBM,NFFT,KPARNUM, & KSORT,NSPAS,ITER) *======================================================================= * Continuation par pseudo longueur d'arc *======================================================================= IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC TMDYNC SEGMENT Mwork REAL*8 FTEST(4), FTEST0(4) REAL*8 XAUX(NFFT), fvec(NT),Vii(NT+2),QOLD1(NT) REAL*8 Rco(NT,NT) ENDSEGMENT INTEGER SPAS, IREDU, METSTB, CPOSRE, IREMAX,INFO LOGICAL CHECK,CHECKBIF, ZINIT,ZSTAB CHARACTER*8 FLAG PARAMETER(IREMAX=6) REAL*8 ZERO,ONE,TWO PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0) REAL*8 Q1NRM2,TVN C A implementer REAL*8 dw0 * .. External Functions .. REAL*8 DNRM2 * Variables generalisees MTQ = KTQ * Partie lineaire MTKAM = KTKAM * Parametres numeriques PARNUM = KPARNUM * Reste des segments MTPHI = KTPHI MTFEX = KTFEX MTPAS = KTPAS MTLIAA = KTLIAA MTLIAB = KTLIAB LOCLFA = KOCLFA LOCLB1 = KOCLB1 MTEMP = KTEMP * Tableau des resultats PSORT = KSORT * NT = Q1(/1) NDDL = NT/(2*NHBM+1) PDT = ONE / NFFT PSORT.CBIF = 0 segini,MWORK *======================================================================= *=== INITIALISATION *======================================================================= * On introduit un terme additionnel aux equations d'equilibre: * * R(Q,w,a) = Z(w)*Q - FNL(Q) + a*(LoIn)*Q = 0 * * Ce dernier etant dissipatif, on a un cycle si et seulement si a=0. * Le vecteur vitesse, (LoIn)*Q, est le vecteur propre associe a la * valeur propre triviale de Rx (du a l'invariance par translation). * Ainsi, son inclusion ne modifie pas le spectre de cette matrice. * Le probleme est desormais bien pose et on continue les solutions par * rapport a l'energie totale du systeme. *----------------------------------------------------------------------- II=1 CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.false.) CALL COPYMAT(NT,RX,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,Rw,NT,INFO) DO I=1,NT dX(I) = -Rw(I) ENDDO a = dnrm2(NT,dX,1) dw = ONE/sqrt(a**2+ONE) DO J=1,NT dX(J) = dX(J)*dw t0(J) = dX(J) ENDDO t0(NT+1) = dw * Sauvegarde des donnees de sortie DO I=1,NT QSAVE(I,II)=Q1(I) ENDDO WSAVE(II)=OMEG *----------------------------------------------------------------------- * A IMPLEMENTER * Stabilite de la solution initiale * ZINIT = .TRUE. * CALL HBMHNNM(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,ONE,ONE,ZINIT, * & LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB) * ZINIT = .FALSE. *----------------------------------------------------------------------- * Message relatif au pas #0 (solution initiale) IF (IIMPI.GE.1) THEN WRITE(IOIMP,667) WRITE(IOIMP,660) WRITE(IOIMP,667) Q1NRM2=dnrm2(NT,Q1,1) WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2 IF (IIMPI.GE.2) WRITE(IOIMP,667) ENDIF *======================================================================= * 1ere prediction *======================================================================= CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,Ra) DO I = 1,NT DO J = 1,NT Jaa(I,J) = RX(I,J) ENDDO Jaa(I,NT+1) = Rw(I) Jaa(I,NT+2) = Ra(I) Jaa(NT+2,I) = dX(I) Jaa(NT+1,I) = -Ra(I) ENDDO Jaa(NT+1,NT+1) = ZERO Jaa(NT+1,NT+2) = ZERO Jaa(NT+2,NT+1) = dw Jaa(NT+2,NT+2) = ZERO DO I = 1,NT+2 Vii(I) = ZERO ENDDO Vii(NT+2) = ONE CALL DGESV(NT+2,1,Jaa,NT+2,IPIV3,Vii,NT+2,INFO) tvn = ISENS * ONE / dnrm2(NT+2,Vii,1) DO I=1,NT dX(I) = Vii(I) * tvn ENDDO dw = Vii(NT+1) * tvn * DO I = 1,NT+1 tP(I) = Vii(I) ENDDO * OMEG = OMEG + DS*dw DO J=1,NT QOLD(J) = Q1(J) Q1(J) = Q1(J)+DS*dX(J) ENDDO IREDU = 0 II = 2 *======================================================================= *======================= Boucle sur la frequence ======================= *======================================================================= DO WHILE (II .LE. NBPAS) * === Correction ================================================= CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM, & MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'CN',ITER) * ----------------------------- * ----------------------------- * ---- si NEWT a converge, ---- * ----------------------------- * ----------------------------- IF (.NOT.CHECK) THEN * Sauvegarde des donnees de sortie * Frequence, coeffs de Fourier, Norme des coeffs de Fourier DO I=1,NT QSAVE(I,II)=Q1(I) ENDDO WSAVE(II)=OMEG * * === Prediction ============================================== * Pas tangent a la courbe DO I = 1,NT+2 Vii(I) = ZERO ENDDO Vii(NT+2) = ONE CALL DGESV(NT+2,1,Jaa,NT+2,IPIV3,Vii,NT+2,INFO) tvn = ONE / dnrm2(NT+1,Vii,1) DO I=1,NT dX(I) = Vii(I) * tvn ENDDO dw = Vii(NT+1) * tvn DO I = 1,NT+1 tP(I) = Vii(I) * tvn ENDDO *----------------------------------------------------------------------- * A IMPLEMENTER * === Stabilite pour la solution convergee ==================== * via Methode de Hill * IF (METSTB.EQ.0) THEN * FTEST0 = FTEST * CALL STBHILL(NT,NHBM,NDDL,OMEG,Q1,Rco,XM,XASM,t0,ZINIT, * & LSAVE(1,1,II),ZSAVE(II),FTEST,FLAG,ZSTAB) * DO J=1,NT * DO I=1,NT * RX(I,J) = ZZ(I,J)-JAC(I,J) * ENDDO * ENDDO * CALL HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,dw0,dw,ZINIT, * & LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB) * via Matrice de Monodromie * ELSE * CALL STBMONO(NT,NDDL,NFFT,NHBM,MTQ,MTKAM,MTPHI,MTLIAA, * & MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,FLAG) * ENDIF * Sauvegarde des donnees de sortie: LSAVE et ZSTAB rempli par STBHILL * WRITE(*,*) 'Pas ',II,', w=',OMEG,'STAB =',ZSTAB *----------------------------------------------------------------------- * === Message relatif au pas #II ============================== IF (IIMPI.GE.1) THEN Q1NRM2=dnrm2(NT,Q1,1) c WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2,ZSTAB WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2 IF (IIMPI.GE.2) WRITE(IOIMP,667) ENDIF * *----------------------------------------------------------------------- * A IMPLEMENTER *=== Bifurcation? ============================================ * IF ((FLAG.NE.'S')) THEN * WRITE(IOIMP,*) 'Bifurcation detectee! Type :',FLAG * CALL HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP, * & MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECKBIF,FLAG, * & PARNUM,PSORT) * IF (.NOT.CHECKBIF) THEN * WRITE(IOIMP,*) 'Bifurcation localisee :)' * ELSE * WRITE(IOIMP,*) 'Bifurcation non localisee :(' * ENDIF * ENDIF *----------------------------------------------------------------------- * === Tests d'arret =========================================== * On verifie la condition d'arret par rapport a OMEG * IF ((OMEG.GE.PARFIN).OR.(OMEG.LE.PARINI)) THEN XPARA = (dnrm2(NT,Q1,1))**2 cbp : incoherence avec la maniere de calculer l'energie dans hbmnewt? IF ((XPARA-PARFIN)*(PARFIN-PARINI).GE.0) THEN WRITE(IOIMP,*) 'Fin de la continuation apres',II,' pas.' IF (II.LT.NBPAS) THEN * WRITE(*,*) 'NPAS sera egal a:',II NT1 = QSAVE(/1) NA1 = LSAVE(/2)/2 * NPAS = QSAVE(/2) NBIFU = QBIFU(/2) * write(*,*) NT1,NA1,NPAS,NBIFU NPAS = II SEGADJ, PSORT ENDIF KSORT = PSORT RETURN ENDIF * === prediction + donnees pour le prochain pas =============== * ajustement automatique de la longueur du pas c if(II.ge.62) write(*,*) '>>>t0=',(t0(iou),iou=1,NT+1) c if(II.ge.62) write(*,*) '>>>tp=',(tp(iou),iou=1,NT+1) CALL ADPAS(DS,DSMIN,DSMAX,ITER,ITERMOY,ANGMIN,ANGMAX, & NT+1,t0,tp) c WRITE(*,*) '>>> DS=',DS * stockage des donnees de fin de pas utiles (dans *OLD) * pas predictif : Q1, OMEG DO I=1,NT+1 t0(I) = tp(I) ENDDO c DSREDU=MAX(0.25*DS,DSMIN) DO KK = 1,NT c c ici, on anticipe la reduction du pas predicteur c c en maintenant la direction dX fixe c QOLD(KK) = Q1(KK)+DSREDU*dX(KK) c et ici, on laisse la possibilite a DX d'evoluer QOLD(KK) = Q1(KK) Q1(KK) = Q1(KK)+DS*dX(KK) ENDDO c idem pour OMEGOLD que QOLD c OMEGOLD = OMEG + DSREDU*dw OMEGOLD = OMEG OMEG = OMEG + DS*dw IREDU = 0 ISTRA2=0 cdebug c if(II.ge.62.and.II.le.68) then c write(*,*) 'prediction Q=',(Q1(iou),iou=1,NT) c write(*,*) 'DS, OMEG=',DS,OMEG c endif * ----------------------------------- * ----------------------------------- * ---- si NEWT n'a pas converge, ---- * ----------------------------------- * ----------------------------------- ELSE * === Message relatif au pas non converge #II ================= IF (IIMPI.GE.1) THEN Q1NRM2=dnrm2(NT,Q1,1) WRITE(IOIMP,665) (II-1),ITER,OMEG,Q1NRM2 IF (IIMPI.GE.2) WRITE(IOIMP,667) ENDIF * === strategie 2 (non-convergence alors qu'on a deja ========= * === atteint DSMIN) : on quitte ! ========= IF (DS.LE.DSMIN) THEN DS=DSMIN DO I=1,NT Q1(I) = QOLD1(I) ENDDO OMEG = OMEGOLD c et on arrete IREDU=IREMAX II = II-1 * === strategie 1 : Recommencer le pas en diminuant DS ======== c test sur II>1 car pour l'instant on n'a pas branche la recup d'une solution initiale ELSEIF (II.GT.1) THEN c c en maintenant la direction dX fixe c c (attention : DS doit etre = DSREDU !!!) c DS = MAX(0.25*DS,DSMIN) c DO I=1,NT c Q1(I) = QOLD(I) c ENDDO c OMEG = OMEGOLD * mise a jour de la prediction avec nouveau DS (dX et dw n'ont pas change) DS = MAX(0.25D0*DS,DSMIN) DO I=1,NT Q1(I) = QOLD(I)+DS*dX(I) ENDDO OMEG = OMEGOLD + DS*dw II = II-1 IREDU = IREDU + 1 * message IF (IIMPI.GE.2) WRITE(IOIMP,668) IREDU,DS ENDIF 888 continue c apres reduction de 0.25**10 = 1.E-6 , echec ! c 0.25**5 = 1.E-3 IF (IREDU.GE.IREMAX) THEN * message interne IF (IIMPI.GE.1) WRITE(IOIMP,669) DS c Pas de convergence apres %i1 iterations. L'execution continue INTERR(1)=ITERMAX CALL ERREUR(151) * ajustement des tableaux de sorties IF (II.LT.NBPAS) THEN NT1 = QSAVE(/1) NA1 = LSAVE(/2)/2 * NPAS = QSAVE(/2) NBIFU = QBIFU(/2) NPAS = II SEGADJ, PSORT ENDIF KSORT = PSORT RETURN ENDIF ENDIF * ---------------------------------------------- * ---------------------------------------------- * ---- fin distinction NEWT converge ou pas ---- * ---------------------------------------------- * ---------------------------------------------- * incrementation du compteur II = II+1 ENDDO *======================================================================= * fin de la boucle sur l'energie *======================================================================= * message c IF(IIMPI.GE.1) WRITE(IOIMP,667) c WRITE(IOIMP,*) NBPAS, 'pas de continuation realises!' KSORT = PSORT SEGSUP,MWORK *======================================================================= * Mise en forme des messages *======================================================================= * dernier message pour fermer le tableau dans le cas iimpi=1 IF (IIMPI.EQ.1) WRITE(IOIMP,667) c 660 FORMAT(' | Pas | Iter | w | Q | Stab |') 660 FORMAT(' | Pas | Iter | w | Q | ') 665 FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ') c 666 FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ',L3,' |') 666 FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ') 667 FORMAT(' +-------+------+-----------+-----------+') 668 FORMAT(' + cont : ds=ds_initial/4**',I2,'=',F13.9) 669 FORMAT(' + cont : pas de convergence avec ds=dsmin=',F13.9) RETURN END