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    CB215821  23/01/25    21:15:19     11573          
  2.  
  3.       SUBROUTINE HBMCO3(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,TVN
  19.        PARAMETER(IREMAX=6)
  20.        segment mwork
  21.          REAL*8 Rco(NT,NT)
  22.          REAL*8 fvec(NT),Vii(NT+2),QOLD1(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 mais ici energie) 
  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,INDNNM
  157.          LOGICAL  JANAL
  158.          REAL*8   E0
  159.       ENDSEGMENT
  160. *
  161. *     Segment des resultats:
  162. *     ---------------------
  163.       SEGMENT  PSORT
  164.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  165.         REAL*8 VSAVE(NPAS)
  166.         LOGICAL ZSAVE(NPAS)
  167.         CHARACTER*2 TYPBIF(NBIFU)
  168.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  169.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  170.         INTEGER CBIF
  171.       ENDSEGMENT
  172. *     QSAVE(i,j) = Q harmonique i au pas j
  173. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  174. *     ZSAVE(j)   = stabilite au j-eme pas
  175. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  176. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  177. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  178. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  179. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  180. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  181.  
  182. *     Segment des tableaux de travail:
  183. *     -------------------------------
  184.       SEGMENT  MTEMP
  185.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  186.          REAL*8 T02(NT1+2), TP2(NT1+2)
  187.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  188.          REAL*8 res
  189.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  190.          REAL*8 QOLD(NT1),OMEGOLD
  191.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  192.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  193.       ENDSEGMENT
  194. *     Jacobiennes augmentees 
  195. *     Ja : [ RX  Rw ; dX dw]
  196. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  197.  
  198. *      SEGMENT  NNNN
  199. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  200. *      ENDSEGMENT
  201.  
  202. *************************** fin TMDYNC.INC *****************************
  203.  
  204.        REAL*8 ZERO,ONE,TWO
  205.        PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0)
  206.  
  207. *     .. External Functions ..
  208.        REAL*8 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. *   On introduit un terme additionnel aux equations d'equilibre:
  238. *     
  239. *      R(Q,w,a) = Z(w)*Q - FNL(Q) + a*(LoIn)*Q = 0
  240. *
  241. *   Ce dernier etant dissipatif, on a un cycle si et seulement si a=0.
  242. *   Le vecteur vitesse, (LoIn)*Q, est le vecteur propre associe a la 
  243. *   valeur propre triviale de Rx (du a l'invariance par translation). 
  244. *   Ainsi, son inclusion ne modifie pas le spectre de cette matrice.
  245. *   Le probleme est desormais bien pose et on continue les solutions par 
  246. *   rapport a l'energie totale du systeme.
  247. *-----------------------------------------------------------------------
  248.        II=1
  249.  
  250.        CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.false.)
  251.        CALL COPYMAT(NT,RX,Rco)
  252.        CALL DGESV(NT,1,Rco,NT,IPIV,Rw,NT,INFO)
  253.        DO I=1,NT
  254.           dX(I) = -Rw(I)
  255.        ENDDO 
  256.        a = dnrm2(NT,dX,1)
  257.        dw = ONE/sqrt(a**2+ONE)  
  258.        DO J=1,NT
  259.           dX(J) = dX(J)*dw
  260.           t0(J) = dX(J)
  261.        ENDDO
  262.        t0(NT+1) = dw
  263.  
  264. *      Sauvegarde des donnees de sortie
  265.        DO I=1,NT
  266.          QSAVE(I,II)=Q1(I)
  267.        ENDDO
  268.        WSAVE(II)=OMEG  
  269.  
  270. *-----------------------------------------------------------------------         
  271. * A IMPLEMENTER
  272. *      Stabilite de la solution initiale
  273. *       ZINIT = .TRUE.
  274. *       CALL HBMHNNM(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,ONE,ONE,ZINIT,
  275. *     &                        LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB)
  276. *       ZINIT = .FALSE.
  277. *----------------------------------------------------------------------- 
  278.  
  279. *      Message relatif au pas #0 (solution initiale)
  280.        IF (IIMPI.GE.1) THEN
  281.          WRITE(IOIMP,667) 
  282.          WRITE(IOIMP,660)
  283.          WRITE(IOIMP,667) 
  284.          Q1NRM2=dnrm2(NT,Q1,1)
  285.          WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2
  286.          IF (IIMPI.GE.2) WRITE(IOIMP,667) 
  287.        ENDIF
  288.  
  289. *=======================================================================  
  290. *      1ere prediction
  291. *======================================================================= 
  292.        CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,Ra)   
  293.        DO I = 1,NT
  294.          DO J = 1,NT
  295.            Jaa(I,J) = RX(I,J)
  296.          ENDDO
  297.          Jaa(I,NT+1)  = Rw(I)
  298.          Jaa(I,NT+2)  = Ra(I)
  299.          Jaa(NT+2,I) = dX(I)
  300.          Jaa(NT+1,I) = -Ra(I)
  301.        ENDDO
  302.        Jaa(NT+1,NT+1) = ZERO
  303.        Jaa(NT+1,NT+2) = ZERO
  304.        Jaa(NT+2,NT+1) = dw
  305.        Jaa(NT+2,NT+2) = ZERO
  306.        DO I = 1,NT+2
  307.          Vii(I) = ZERO
  308.        ENDDO
  309.        Vii(NT+2) = ONE       
  310.        CALL DGESV(NT+2,1,Jaa,NT+2,IPIV3,Vii,NT+2,INFO)
  311.        tvn = dnrm2(NT+2,Vii,1)
  312.        DO I=1,NT
  313.           dX(I) = ISENS*Vii(I)/tvn
  314.        ENDDO 
  315.        dw = ISENS*Vii(NT+1)/tvn
  316. *       
  317.        DO I = 1,NT+1
  318.           tP(I) = Vii(I)
  319.        ENDDO
  320. *       
  321.        OMEG = OMEG + DS*dw
  322.        DO J=1,NT
  323.          QOLD(J) = Q1(J)
  324.            Q1(J) = Q1(J)+DS*dX(J)
  325.        ENDDO        
  326.  
  327.  
  328.        IREDU = 0
  329.        II = 2
  330.  
  331. *=======================================================================  
  332. *======================= Boucle sur la frequence =======================        
  333. *=======================================================================  
  334.        DO WHILE (II .LE. NBPAS)
  335. *       === Correction =================================================
  336.         CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
  337.      &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'CN',ITER)   
  338.  
  339. *       -----------------------------
  340. *       -----------------------------
  341. *       ---- si NEWT a converge, ----
  342. *       -----------------------------
  343. *       -----------------------------
  344.         IF (.NOT.CHECK) THEN
  345.  
  346. *          Sauvegarde des donnees de sortie
  347. *          Frequence, coeffs de Fourier, Norme des coeffs de Fourier           
  348.            DO I=1,NT
  349.              QSAVE(I,II)=Q1(I)
  350.            ENDDO
  351.            WSAVE(II)=OMEG           
  352. * 
  353. *          === Prediction ==============================================
  354. *          Pas tangent a la courbe 
  355.            DO I = 1,NT+2
  356.              Vii(I) = ZERO
  357.            ENDDO
  358.            Vii(NT+2) = ONE       
  359.            CALL DGESV(NT+2,1,Jaa,NT+2,IPIV3,Vii,NT+2,INFO)
  360.            tvn = dnrm2(NT+1,Vii,1)
  361.            DO I=1,NT
  362.              dX(I) = Vii(I)/tvn
  363.            ENDDO 
  364.            dw = Vii(NT+1)/tvn
  365.            DO I = 1,NT+1
  366.              tP(I) = Vii(I)/tvn
  367.            ENDDO
  368.  
  369. *-----------------------------------------------------------------------
  370. * A IMPLEMENTER
  371. *          === Stabilite pour la solution convergee ====================
  372. *          via Methode de Hill
  373. *           IF (METSTB.EQ.0) THEN
  374. *             FTEST0 = FTEST
  375. *             CALL STBHILL(NT,NHBM,NDDL,OMEG,Q1,Rco,XM,XASM,t0,ZINIT,
  376. *     &                    LSAVE(1,1,II),ZSAVE(II),FTEST,FLAG,ZSTAB)
  377. *             DO J=1,NT
  378. *             DO I=1,NT
  379. *                 RX(I,J) = ZZ(I,J)-JAC(I,J)
  380. *             ENDDO           
  381. *             ENDDO
  382. *             CALL HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,dw0,dw,ZINIT,
  383. *     &                        LSAVE(1,1,II),ZSAVE(II),CPOSRE,FLAG,ZSTAB)
  384. *          via Matrice de Monodromie    
  385. *           ELSE 
  386. *             CALL STBMONO(NT,NDDL,NFFT,NHBM,MTQ,MTKAM,MTPHI,MTLIAA,                   
  387. *    &          MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,FLAG)
  388. *           ENDIF           
  389. *          Sauvegarde des donnees de sortie: LSAVE et ZSTAB rempli par STBHILL
  390. *          WRITE(*,*) 'Pas ',II,', w=',OMEG,'STAB =',ZSTAB
  391. *-----------------------------------------------------------------------
  392.  
  393. *          === Message relatif au pas #II ==============================
  394.            IF (IIMPI.GE.1) THEN
  395.               Q1NRM2=dnrm2(NT,Q1,1)
  396. c               WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2,ZSTAB
  397.               WRITE(IOIMP,666) (II-1),ITER,OMEG,Q1NRM2
  398.               IF (IIMPI.GE.2) WRITE(IOIMP,667) 
  399.            ENDIF
  400. *  
  401. *-----------------------------------------------------------------------
  402. * A IMPLEMENTER
  403. *=== Bifurcation? ============================================
  404. *           IF ((FLAG.NE.'S')) THEN 
  405. *             WRITE(IOIMP,*) 'Bifurcation detectee! Type :',FLAG
  406. *               CALL HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,
  407. *     &         MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECKBIF,FLAG,
  408. *     &         PARNUM,PSORT)
  409. *               IF (.NOT.CHECKBIF) THEN
  410. *                   WRITE(IOIMP,*) 'Bifurcation localisee :)'
  411. *               ELSE
  412. *                   WRITE(IOIMP,*) 'Bifurcation non localisee :('
  413. *               ENDIF
  414. *           ENDIF
  415. *-----------------------------------------------------------------------
  416.  
  417. *          === Tests d'arret ===========================================        
  418. *          On verifie la condition d'arret par rapport a OMEG
  419. *           IF ((OMEG.GE.PARFIN).OR.(OMEG.LE.PARINI)) THEN
  420.            XPARA = (dnrm2(NT,Q1,1))**2
  421. cbp : incoherence avec la maniere de calculer l'energie dans hbmnewt?
  422.            IF ((XPARA-PARFIN)*(PARFIN-PARINI).GE.0) THEN
  423.                WRITE(IOIMP,*) 'Fin de la continuation apres',II,' pas.'
  424.                IF (II.LT.NBPAS) THEN
  425. *                  WRITE(*,*) 'NPAS sera egal a:',II
  426.                   NT1 = QSAVE(/1)
  427.                   NA1 = LSAVE(/2)/2
  428. *                  NPAS = QSAVE(/2)
  429.                   NBIFU = QBIFU(/2)
  430. *                  write(*,*) NT1,NA1,NPAS,NBIFU
  431.                   NPAS = II
  432.                   SEGADJ, PSORT
  433.                ENDIF
  434.                KSORT = PSORT
  435.                RETURN
  436.            ENDIF
  437.  
  438. *          === prediction + donnees pour le prochain pas ===============
  439. *          ajustement automatique de la longueur du pas
  440. c            if(II.ge.62) write(*,*) '>>>t0=',(t0(iou),iou=1,NT+1)
  441. c            if(II.ge.62) write(*,*) '>>>tp=',(tp(iou),iou=1,NT+1)
  442.            CALL ADPAS(DS,DSMIN,DSMAX,ITER,ITERMOY,ANGMIN,ANGMAX,
  443.      &                NT+1,t0,tp)
  444. c            WRITE(*,*) '>>> DS=',DS
  445. *          stockage des donnees de fin de pas utiles (dans *OLD)
  446. *          pas predictif : Q1, OMEG
  447.            DO I=1,NT+1
  448.               t0(I) = tp(I)
  449.            ENDDO
  450. c            DSREDU=MAX(0.25*DS,DSMIN)
  451.            DO KK = 1,NT
  452. c c             ici, on anticipe la reduction du pas predicteur 
  453. c c             en maintenant la direction dX fixe
  454. c               QOLD(KK) = Q1(KK)+DSREDU*dX(KK)
  455. c             et ici, on laisse la possibilite a DX d'evoluer
  456.               QOLD(KK) = Q1(KK)
  457.               Q1(KK)   = Q1(KK)+DS*dX(KK)
  458.            ENDDO
  459. c          idem pour OMEGOLD que QOLD
  460. c            OMEGOLD = OMEG + DSREDU*dw
  461.            OMEGOLD = OMEG
  462.            OMEG    = OMEG + DS*dw
  463.            IREDU = 0
  464.            ISTRA2=0
  465. cdebug           
  466. c            if(II.ge.62.and.II.le.68) then
  467. c            write(*,*) 'prediction Q=',(Q1(iou),iou=1,NT)
  468. c            write(*,*) 'DS, OMEG=',DS,OMEG
  469. c            endif
  470.  
  471. *       -----------------------------------
  472. *       -----------------------------------
  473. *       ---- si NEWT n'a pas converge, ----
  474. *       -----------------------------------
  475. *       -----------------------------------
  476.         ELSE
  477.  
  478. *          === Message relatif au pas non converge #II =================
  479.            IF (IIMPI.GE.1) THEN
  480.               Q1NRM2=dnrm2(NT,Q1,1)
  481.               WRITE(IOIMP,665) (II-1),ITER,OMEG,Q1NRM2
  482.               IF (IIMPI.GE.2) WRITE(IOIMP,667) 
  483.            ENDIF
  484.  
  485. *          === strategie 2 (non-convergence alors qu'on a deja ========= 
  486. *          ===              atteint DSMIN) : on quitte !       =========
  487.            IF (DS.LE.DSMIN) THEN
  488.              DS=DSMIN
  489.              DO I=1,NT
  490.                Q1(I) = QOLD1(I)
  491.              ENDDO
  492.              OMEG = OMEGOLD           
  493. c            et on arrete
  494.              IREDU=IREMAX
  495.              II = II-1
  496.  
  497. *          === strategie 1 : Recommencer le pas en diminuant DS ======== 
  498. c          test sur II>1 car pour l'instant on n'a pas branche la recup d'une solution initiale
  499.            ELSEIF (II.GT.1) THEN
  500. c c            en maintenant la direction dX fixe 
  501. c c            (attention : DS doit etre = DSREDU !!!)
  502. c              DS = MAX(0.25*DS,DSMIN)
  503. c              DO I=1,NT
  504. c                Q1(I) = QOLD(I)
  505. c              ENDDO
  506. c              OMEG = OMEGOLD
  507. *            mise a jour de la prediction avec nouveau DS (dX et dw n'ont pas change)
  508.              DS = MAX(0.25D0*DS,DSMIN)
  509.              DO I=1,NT
  510.                Q1(I) = QOLD(I)+DS*dX(I)
  511.              ENDDO
  512.              OMEG = OMEGOLD + DS*dw
  513.              II = II-1
  514.              IREDU = IREDU + 1
  515. *            message
  516.              IF (IIMPI.GE.2) WRITE(IOIMP,668) IREDU,DS
  517.  
  518.            ENDIF
  519.  
  520.  
  521.  888       continue                  
  522. c          apres reduction de 0.25**10 = 1.E-6 , echec !
  523. c                             0.25**5  = 1.E-3
  524.            IF (IREDU.GE.IREMAX) THEN        
  525. *             message interne
  526.               IF (IIMPI.GE.1) WRITE(IOIMP,669) DS
  527. c             Pas de convergence apres %i1 iterations. L'execution continue
  528.               INTERR(1)=ITERMAX              
  529.               CALL ERREUR(151)
  530. *             ajustement des tableaux de sorties
  531.               IF (II.LT.NBPAS) THEN
  532.                   NT1 = QSAVE(/1)
  533.                   NA1 = LSAVE(/2)/2
  534. *                  NPAS = QSAVE(/2)
  535.                   NBIFU = QBIFU(/2)
  536.                   NPAS = II
  537.                   SEGADJ, PSORT
  538.               ENDIF
  539.               KSORT = PSORT
  540.               RETURN
  541.            ENDIF
  542.  
  543.         ENDIF 
  544. *       ----------------------------------------------
  545. *       ----------------------------------------------
  546. *       ---- fin distinction NEWT converge ou pas ----  
  547. *       ----------------------------------------------
  548. *       ----------------------------------------------
  549.  
  550. *        incrementation du compteur
  551.          II = II+1
  552.  
  553.         ENDDO
  554. *=======================================================================  
  555. *       fin de la boucle sur l'energie        
  556. *=======================================================================  
  557.  
  558.  
  559. *       message
  560. c         IF(IIMPI.GE.1) WRITE(IOIMP,667) 
  561. c       WRITE(IOIMP,*) NBPAS, 'pas de continuation realises!' 
  562.         KSORT = PSORT
  563.  
  564.         segsup,MWORK
  565.  
  566. *=======================================================================
  567. *      Mise en forme des messages
  568. *=======================================================================  
  569. *      dernier message pour fermer le tableau dans le cas iimpi=1
  570.        IF (IIMPI.EQ.1) WRITE(IOIMP,667) 
  571.  
  572. c  660   FORMAT(' |  Pas  | Iter |     w     |     Q     | Stab |')
  573.  660   FORMAT(' |  Pas  | Iter |     w     |     Q     | ')
  574.  665   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ')
  575. c  666   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ',L3,'  |')
  576.  666   FORMAT(' | ',I5,' | ',I4,' | ',F9.5,' | ',F9.5,' | ')
  577.  667   FORMAT(' +-------+------+-----------+-----------+')
  578.  668   FORMAT(' + cont : ds=ds_initial/4**',I2,'=',F13.9)
  579.  669   FORMAT(' + cont : pas de convergence avec ds=dsmin=',F13.9)
  580.  
  581.  
  582.       END
  583.  
  584.  
  585.  
  586.  
  587.  
  588.  
  589.  

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