Télécharger la plaquette Cast3M


fr en

Télécharger hbmnewt.eso

Retour à la liste

Numérotation des lignes : 

hbmnewt

  1. C HBMNEWT   SOURCE    CB215821  23/01/25    21:15:22     11573          
  2.  
  3.       SUBROUTINE HBMNEWT(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,KPARNUM
  4.      &   ,KTLIAA,MTLIAB,KTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,TYPC,ITER)
  5. *
  6. *=======================================================================
  7. *
  8. *                     Methode de Newton-Raphson  
  9. *
  10. * Convergence vers une solution des equations d'equilibre (frequentiel).
  11. * TYPC(1) indique le type de calcul a realiser:
  12. *      = | 'I' : solution initiale 
  13. *        | 'C' : pas correcteur de la continuation 
  14. * TYPC(2) indique le type de calcul a realiser:
  15. *      = | 'F' : reponse Forcee
  16. *        | 'A' : reponse Autonome
  17. *        | 'N' : reponse libre : Non-linear normal modes
  18. *
  19. *=======================================================================     
  20.  
  21. *----- Declarations ----------------------------------------------------
  22.  
  23.        IMPLICIT INTEGER(I-N)
  24.        IMPLICIT REAL*8(A-H,O-Z)
  25.        REAL*8 fvec(NT),AVVP(NT),VCTCS(7),LAMBD(NDDL),AiDi(2),Di(2),bi(2)
  26.        REAL*8 J2(NT+1,NT+1), Ja2(NT+2,NT+2)
  27.        LOGICAL check,JANAL2
  28. c        CHARACTER*2 TYPC
  29.        INTEGER INFO
  30.        REAL*8 fref,fref0
  31.        REAL*8 gx(NT)
  32.  
  33. -INC SMLREEL
  34. -INC PPARAM
  35. -INC CCOPTIO
  36. *-INC TMDYNC.INC
  37. ************************** debut TMDYNC.INC ****************************
  38.  
  39. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  40. *    TODO   : a extraire dans un include des que stabilise
  41. *
  42. *     Segment des variables generalisees:
  43. *     -----------------------------------
  44.       SEGMENT MTQ
  45.          REAL*8 Q1(NT1)
  46.          REAL*8 OMEG,XPARA
  47.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  48.          REAL*8 dX(NT1), dw, dv
  49.       ENDSEGMENT
  50. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  51. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  52. *            avec q_i vecteur de dimension n ou n=nombre de modes
  53. *     OMEG : frequence fondamentale de l'approximation
  54. *     XPARA: parametre de continuation (par defaut la frequence) 
  55. *            \in [PARINI,PARFIN]
  56. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  57. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  58. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  59. *            lineaires et constantes
  60. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  61. *
  62. *
  63. *     Segment contenant les matrices XK, XASM et XM:
  64. *     ---------------------------------------------
  65.       SEGMENT MTKAM
  66.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  67.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  68. *         REAL*8 GAMFIN(NPC2,nl1)
  69.       ENDSEGMENT
  70. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  71. *     GAM et IGAM : matrices pour la FFT et son inverse
  72. *     GAMFIN      : 
  73. *
  74. *     Segment des deformees modales:
  75. *     ------------------------------
  76. *     (idem DYNE)
  77.       SEGMENT MTPHI
  78.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  79.          INTEGER IAROTA(NSB)
  80.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  81.       ENDSEGMENT
  82. *
  83. *     Segment descriptif des liaisons en base A:
  84. *     ------------------------------------------
  85. *     (idem DYNE)
  86.       SEGMENT MTLIAA
  87.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  88.          REAL*8  XPALA(NLIAA,NXPALA)
  89.       ENDSEGMENT
  90. *
  91. *     Segment descriptif des liaisons en base B:
  92. *     ------------------------------------------
  93. *     (idem DYNE)
  94.       SEGMENT MTLIAB
  95.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  96.          REAL*8  XPALB(NLIAB,NXPALB)
  97.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  98.       ENDSEGMENT
  99. *
  100. *     Segment representant les chargements exterieurs:
  101. *     -----------------------------------------------
  102.       SEGMENT MTFEX
  103.          REAL*8  FEXA(NT1)
  104.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  105.          INTEGER BAL
  106.       ENDSEGMENT
  107. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  108. *            Fourier et exprimes en base A
  109. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  110. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  111. *             compatibilite avec calcul des Fnl.
  112. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  113. *            (cad proportionnel a OMEG**2)
  114. *
  115. *     Segment "local" pour DEVLFA:
  116. *     ----------------------------
  117.       SEGMENT LOCLFA
  118.          REAL*8 FTEST(NA1,4)
  119.       ENDSEGMENT
  120. *      
  121. *     Segment "local" pour DEVLB1:
  122. *     ----------------------------
  123.       SEGMENT LOCLB1
  124.          REAL*8 FTEST2(NPLB,6)
  125.       ENDSEGMENT
  126. *      
  127. *     Segment contenant les variables au cours d un pas de temps:
  128. *     ----------------------------------------------------------
  129.       SEGMENT MTPAS
  130.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  131.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  132.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  133.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  134.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  135.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  136.       ENDSEGMENT
  137. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  138. *     XPTB       : deplacement du point d'une liaison en base B
  139. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  140. *     FEXB       : forces exterieures en base B (a priori uniquement 
  141. *                  pour les moments appliques aux rotations rigides ?)
  142. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  143. *                  maillage de contact dans certaines liaisons)
  144. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  145. *                       base A/B (= contributions a dFnl/dX)
  146. *
  147. *
  148. *     Segment des points de reference des modes (base A):
  149. *     --------------------------------------------------
  150.       SEGMENT MPREF
  151.          INTEGER IPOREF(NPREF)
  152.       ENDSEGMENT
  153. *
  154. *     Segment des points en base B:
  155. *     -----------------------------
  156.       SEGMENT NCPR(NBPTS)
  157. *     NCRP(#global) = #local dans XPTB (1er indice)      
  158. *
  159. *     Segment des parametres numeriques pour la continuation: 
  160. *     ------------------------------------------------------
  161.       SEGMENT PARNUM
  162.          CHARACTER*4 TYPS
  163.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  164.          REAL*8   PARINI,PARFIN
  165.          INTEGER  ITERMAX,NBPAS
  166.          LOGICAL  JANAL
  167.       ENDSEGMENT
  168. *
  169. *     Segment des resultats:
  170. *     ---------------------
  171.       SEGMENT  PSORT
  172.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  173.         REAL*8 VSAVE(NPAS)
  174.         LOGICAL ZSAVE(NPAS)
  175.         CHARACTER*2 TYPBIF(NBIFU)
  176.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  177.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  178.         INTEGER CBIF
  179.       ENDSEGMENT
  180. *     QSAVE(i,j) = Q harmonique i au pas j
  181. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  182. *     ZSAVE(j)   = stabilite au j-eme pas
  183. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  184. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  185. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  186. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  187. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  188. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  189.  
  190. *     Segment des tableaux de travail:
  191. *     -------------------------------
  192.       SEGMENT  MTEMP
  193.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  194.          REAL*8 T02(NT1+2), TP2(NT1+2)
  195.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  196.          REAL*8 res
  197.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  198.          REAL*8 QOLD(NT1),OMEGOLD
  199.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  200.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  201.       ENDSEGMENT
  202. *     Jacobiennes augmentees 
  203. *     Ja : [ RX  Rw ; dX dw]
  204. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  205.  
  206.  
  207. *************************** fin TMDYNC.INC *****************************
  208.  
  209.        CHARACTER*2 TYPC
  210.        REAL*8 ZERO,TWO,ONE
  211.        PARAMETER (ZERO=0.D0, TWO=2.D0,ONE=1.D0)
  212.        LOGICAL   NNM
  213. C Fonctions BLAS/LAPACK
  214.        REAL*8    DDOT, DNRM2
  215.        EXTERNAL  DDOT, DNRM2
  216.  
  217. c        NNM=(TYPC.EQ.'IN').OR.(TYPC.EQ.'CN')
  218.        NNM=TYPC(2:2).EQ.'N'
  219. *
  220. *      Segments
  221.        MTQ   = KTQ
  222.        MTKAM = KTKAM
  223.        MTEMP = KTEMP
  224.        PARNUM = KPARNUM
  225.        MTLIAA= KTLIAA 
  226. c        MTLIAB= KTLIAB
  227.        MTFEX = KTFEX  
  228.        ITER=-1
  229. c        RELAX=1.d0
  230.        JANAL2=JANAL
  231.  
  232. c *------ Test de convergence pour la solution ------------------
  233. c         test=0.D0
  234. c         DO i=1,NT
  235. c            IF(abs(fvec(i)).GT.test) test=abs(fvec(i))
  236. c         ENDDO
  237.  
  238.  
  239. *------ Calcul de la raideur dynamique tangente des efforts lineaires --
  240.  
  241. *       Matrice liee a la base modale (avec ou sans terme dissipatif selon NNM)
  242.         CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.not.NNM)   
  243.  
  244. *       cas Granger-Paidoussis : rajouter raideur de couplage
  245. c       aussi pour cas NNM ? a voir...
  246.         IF(.not.NNM) THEN
  247.           DO J = 1,IPALA(/1)
  248.             IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  249.              DO IJ = 1,7
  250.                VCTCS(IJ) = XPALA(J,IJ)
  251.              ENDDO
  252.              IF (TYPC.EQ.'CA') THEN 
  253.                VCTCS(4) = VV
  254.                XPALA(J,4) = VV
  255.              ENDIF
  256.              JG = NDDL
  257.              MLREE1 = IPALA(J,7)
  258. c              SEGACT, MLREE1
  259.              DO IJ = 1,JG
  260.                LAMBD(IJ) = MLREE1.PROG(IJ)
  261.              ENDDO
  262.              MLREE2 = IPALA(J,4)
  263.              MLREE3 = IPALA(J,8)
  264. c              SEGACT, MLREE2,MLREE3
  265. c *            Nombre de termes fixe, a generaliser si besoin.            
  266. c              JG = 2
  267. *            tentative de generalisation
  268.              JG = MLREE2.PROG(/1)
  269.              DO IJ = 1,JG
  270.                AiDi(IJ) = MLREE2.PROG(IJ)
  271.                Di(IJ)   = MLREE3.PROG(IJ)
  272.              ENDDO
  273. c              SEGDES,MLREE1,MLREE2,MLREE3     
  274.              CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)
  275.             ENDIF
  276.           ENDDO
  277.         ENDIF
  278.  
  279. *------ Calcul du residu initial ---------------------------------------
  280. *       (lie a la prediction par exemple)
  281.         its=0
  282.         CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  283.      &         MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,JANAL2,fvec,fref0,NNM)
  284.  
  285.         IF (IIMPI.GE.2) THEN
  286.           Q1NRM2=dnrm2(NT,Q1,1)
  287.           res=dnrm2(NT,fvec,1)/fref0
  288.           WRITE(IOIMP,990) its,OMEG,Q1NRM2,res,fref0,DS
  289. c           IF (IIMPI.GE.333) THEN
  290. c           WRITE(IOIMP,111) 'Q',its,(Q1(iou),iou=1,NT)
  291. c           WRITE(IOIMP,111) 'R',its,(fvec(iou),iou=1,NT)
  292. c c           do iou=1,NT
  293. c c             WRITE(IOIMP,111) 'ZZ_',iou,(ZZ(iou,jou),jou=1,NT)
  294. c c           enddo
  295. c           ENDIF
  296.         ENDIF
  297.  
  298.  
  299. *============================= Iterations ==============================      
  300.         DO its=1,ITERMAX 
  301.  
  302.  
  303. *    ----- DIFFICULTE A CONVERGER : ON RELAXE (test non concluant) ----
  304. c            IF (its .GT. 10) THEN
  305. c              RELAX=0.222d0             
  306. c            ENDIF
  307.  
  308. *    ----- CALCUL DE LA JACOBIENNE : dR/dX = Z(OMEGA) - dFNLdX ----
  309.            DO J=1,NT
  310.            DO I=1,NT
  311.                  RX(I,J) = ZZ(I,J)-JAC(I,J)
  312.            ENDDO 
  313.            ENDDO
  314. c            IF (IIMPI.GE.333)  THEN 
  315. c              DO iou=1,NT
  316. c                WRITE(IOIMP,111) 'RX_',iou,(RX(iou,jou),jou=1,NT)
  317. c              ENDDO
  318. c            ENDIF
  319.  
  320. *    ----- CAS 1: INITIALISATION (systeme force) ------------------
  321.            IF (TYPC.EQ.'IF') THEN
  322. *           Calcul de la correction dX, telle que: Rx*dX = R
  323. *           Appel a "DGESV" (LAPACK)
  324.             CALL DGESV(NT,1,RX,NT,IPIV,fvec,NT,INFO)
  325.             DO M = 1,NT
  326.                Q1(M) = Q1(M)-fvec(M)
  327.             ENDDO
  328. c             IF (IIMPI.GE.333)
  329. c      &      WRITE(IOIMP,111) 'dQ^',its,(fvec(iou),iou=1,NT)
  330. c           si on souhaite verifier la convergence en correction dX
  331. c           res=dnrm2(NT,fvec,1)/dnrm2(NT,Q1,1)               
  332.            ENDIF    
  333.  
  334. *    ----- CAS 2: CONTINUATION (systeme force) --------------------
  335.            IF (TYPC.EQ.'CF') THEN    
  336. *           Calcul de la correction dY, telle que: Ja*dY = {R,0}
  337. *           avec Ja la Jacobienne augmentee: Ja = {RX,Rw;dX,dw}
  338.             DO JJ=1,NT
  339.             DO II=1,NT
  340.                 Ja(II,JJ) = RX(II,JJ)
  341.             ENDDO
  342.             ENDDO
  343. *           Calcul de la derivee: R,w = dR/dw 
  344.             CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  345. *           cas Granger-Paidoussis : R,w = R,w - dFf/dw 
  346.             DO J = 1,IPALA(/1)
  347.               IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  348.                 VCTCS(4) = VV
  349.                 CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)     
  350.               ENDIF
  351.             ENDDO
  352. c           cas d'un balourd : force en OMEG**2 a deriver aussi
  353.             IF (BAL.EQ.1) THEN
  354.               DO I=1,NT
  355.                 Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  356.               ENDDO
  357.             ENDIF
  358.             DO KK = 1,NT
  359.                 Ja(KK,NT+1) = Rw(KK)
  360.                 Ja(NT+1,KK) = dX(KK)
  361. ctest -> nechange rien :  Ja(NT+1,KK) = DS*dX(KK)
  362. ctest                Ja(NT+1,KK) = (Q1(KK)-QOLD(KK))
  363.                 RHS(KK) = fvec(KK)
  364. cdebug
  365. c             if(OMEG.gt.0.153.and.OMEG.lt.0.170) 
  366. c      &      
  367. c               write(*,123) KK,(Ja(KK,iou),iou=1,6),RHS(KK)
  368.             ENDDO
  369.             Ja(NT+1,NT+1) = dw
  370. ctest -> nechange rien :      Ja(NT+1,NT+1) = DS*dw
  371. ctest            Ja(NT+1,NT+1) = OMEG-OMEGOLD
  372.             RHS(NT+1) = ZERO
  373. c             if(OMEG.gt.0.153.and.OMEG.lt.0.170) 
  374. c      &      
  375. c             write(*,123) NT+1,(Ja(NT+1,iou),iou=1,6),RHS(NT+1)
  376. *           Appel a "DGESV" (LAPACK)
  377.             CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO)
  378.             DO MM = 1,NT
  379. c                Q1(MM) = Q1(MM) - RELAX*RHS(MM)
  380.                Q1(MM) = Q1(MM) - RHS(MM)
  381.             ENDDO
  382. c             OMEG = OMEG - RELAX*RHS(NT+1)
  383.             OMEG = OMEG - RHS(NT+1)
  384.             CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  385. c           si on souhaite verifier la convergence en correction dX
  386. c             res=dnrm2(NT,RHS,1)/dnrm2(NT,Q1,1) 
  387. c             write(*,*) 'W,Q1=',OMEG,(Q1(iou),iou=1,NT)
  388.            ENDIF
  389.  
  390. *    ----- CAS 3: INITIALISATION (systeme autonome) --------------
  391.            IF (TYPC.EQ.'IA') THEN
  392. *          Calcul de la correction dY, telle que: Ja*dY = -[R;g]
  393. *          avec: Ja = [Rx Rw; gx 0]  
  394.             DO II=1,NT
  395.               DO JJ = 1,NT
  396.                 Ja(II,JJ) = RX(II,JJ)
  397.               ENDDO
  398.             ENDDO
  399.             VV = VCTCS(4)
  400. *           Calcul de la derivee: Rw = dR/dw 
  401.             CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)         
  402. *           cas Granger-Paidoussis : R,w = R,w - dFf/dw 
  403.             DO J = 1,IPALA(/1)
  404.               IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  405.                CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)     
  406.               ENDIF
  407.             ENDDO
  408.             DO KK = 1,NT
  409.               Ja(KK,NT+1) = Rw(KK)
  410.               RHS(KK) = fvec(KK)
  411.             ENDDO
  412. *           Condition de phase 
  413. *           On impose une vitesse nulle au 1er DDL en t = 0:
  414. *           dx1(0) = 0 ==> sum_i{i*Q1_si} = 0
  415.             RHS(NT+1) = ZERO 
  416.             DO I = 1,NT+1
  417.               Ja(NT+1,I) = ZERO
  418.             ENDDO
  419.             DO JJ = 2,2*NHBM,2
  420.               Ja(NT+1,JJ*NDDL+1) = JJ/TWO
  421. *              RHS(NT+1)  = RHS(NT+1) + (JJ/TWO)*Q1(JJ*NDDL+1)
  422.             ENDDO
  423. c           si on souhaite verifier la convergence en correction dX
  424. c             res=dnrm2(NT,RHS,1)              
  425.             Ja(NT+1,NT+1) = ZERO        
  426.             CALL COPYMAT(NT+1,Ja,J2) 
  427. *           Appel a "DGESV" (LAPACK)
  428.             CALL DGESV(NT+1,1,J2,NT+1,IPIV2,RHS,NT+1,INFO)
  429.             DO M = 1,NT
  430.                Q1(M) = Q1(M)-RHS(M)
  431.             ENDDO
  432.             OMEG = OMEG-RHS(NT+1)
  433.             CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.) 
  434. *           cas Granger-Paidoussis : rajouter raideur de couplage
  435.             DO J = 1,IPALA(/1)
  436.               IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  437.                CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)    
  438.               ENDIF
  439.             ENDDO
  440.            ENDIF
  441.  
  442. *    ----- CAS 4: CONTINUATION (systeme autonome) -----------------
  443.            IF (TYPC.EQ.'CA') THEN  
  444. *    Calcul de la correction dY, telle que: Jaa*dY = -[R;g;0]
  445. *    avec Jaa la Jacobienne augmentee: Jaa ={RX,Rw,Ra;gx,0,0;dX,dw,da}
  446.             DO II=1,NT
  447.               DO JJ = 1,NT
  448.                 Jaa(II,JJ) = RX(II,JJ)
  449.               ENDDO
  450.             ENDDO
  451. *           Calcul de la derivee: Rw = dR/dw 
  452.             CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) 
  453. *           cas Granger-Paidoussis : R,w = R,w - dFf/dw 
  454.             DO J = 1,IPALA(/1)
  455.               IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  456.                 VCTCS(4) = VV
  457.                 CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)     
  458.               ENDIF
  459.             ENDDO 
  460.             DO I = 1,NT+2
  461.               Jaa(NT+2,I) = ZERO
  462.             ENDDO
  463. *           Calcul de la derivee : R,Vmoy = dR/dVmoy
  464.             CALL HBMRV(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Ra)         
  465.             DO KK = 1,NT
  466.                 Jaa(KK,NT+1) = Rw(KK)
  467.                 Jaa(KK,NT+2) = Ra(KK)
  468.                 Jaa(NT+2,KK) = dX(KK)
  469.                 Jaa(NT+1,KK) = ZERO
  470.                 RHS2(KK) = fvec(KK)    
  471.             ENDDO
  472.             Jaa(NT+2,NT+1) = dw
  473.             Jaa(NT+2,NT+2) = dv
  474. *           Condition de phase 
  475. *           On impose une vitesse nulle au 1er DDL en t = 0:
  476. *           dx1(0) = 0 ==> g(Q) = sum_i{i*Q1_si} = 0
  477. *           dg/dw
  478.             Jaa(NT+1,NT+1) = ZERO
  479. *           dg/dv
  480.             Jaa(NT+1,NT+2) = ZERO
  481.             RHS2(NT+1) = ZERO 
  482.             DO JJ = 2,2*NHBM,2
  483.               Jaa(NT+1,JJ*NDDL+1) = JJ/2
  484.               RHS2(NT+1)  = RHS2(NT+1) + (JJ/2)*Q1(JJ*NDDL+1)
  485.             ENDDO
  486.             RHS2(NT+2) = ZERO
  487. *            RHS2(NT+2) = ZERO
  488. *            WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)
  489. *            STOP
  490.             CALL COPYMAT(NT+2,Jaa,Ja2) 
  491. *           Appel a "DGESV" (LAPACK)
  492.             CALL DGESV(NT+2,1,Ja2,NT+2,IPIV3,RHS2,NT+2,INFO)
  493. c           si on souhaite verifier la convergence en correction dX
  494. c             res=dnrm2(NT+1,RHS2,1) 
  495. c             WRITE(*,*) 'res=',res
  496. *            WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)        
  497.             DO M = 1,NT
  498.                Q1(M) = Q1(M)-RHS2(M)
  499.             ENDDO
  500.             OMEG = OMEG-RHS2(NT+1)
  501.             VV   = VV  -RHS2(NT+2)
  502.             CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  503.             DO J = 1,IPALA(/1)
  504.               IF (IPALA(J,3).EQ.101) THEN
  505.                VCTCS(4) = VV
  506.                XPALA(J,4) = VV
  507.                CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)    
  508.                ENDIF
  509.             ENDDO      
  510.            ENDIF
  511.  
  512. *    ----- CAS 5: INITIALISATION (systeme Hamiltonien) --------------
  513.            IF (TYPC.EQ.'IN') THEN
  514. *          Calcul de la correction dY, telle que: Ja*dY = -[R;g]
  515. *          avec: Ja = [Rx Rw; gx 0]  
  516.             DO II=1,NT
  517.               DO JJ = 1,NT
  518.                 Ja(II,JJ) = RX(II,JJ)
  519.               ENDDO
  520.             ENDDO
  521. *           Calcul de la derivee: Rw = dR/dw 
  522.             CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.false.)
  523.             DO KK = 1,NT
  524.               Ja(KK,NT+1) = Rw(KK)
  525.               RHS(KK) = fvec(KK)
  526.             ENDDO
  527. *           Condition d'energie
  528. *           On initialise en imposant un certain niveau d'energie,
  529. *           i.e. on fixe la norme de Q1
  530.             RHS(NT+1) = 0.5*(DDOT(NT,Q1,1,Q1,1))-XPARA
  531. cbp    <=>  RHS(NT+1) = -0.5*(DDOT(NT,Q1,1,Q1,1)) 
  532. c           car XPARA=DDOT(Q1)=sqrt(Energie_initiale) dans hbmalo
  533. cbp : incoherence avec calcul energie realise dans hbmco3
  534. c             RHS(NT+1) = DDOT(NT,Q1,1,Q1,1)-XPARA**2
  535.             DO I = 1,NT
  536.               Ja(NT+1,I) = Q1(I)
  537.             ENDDO
  538. c           si on souhaite verifier la convergence en correction dX
  539.             res=dnrm2(NT+1,RHS,1)              
  540.             Ja(NT+1,NT+1) = ZERO        
  541. *            CALL COPYMAT(NT+1,Ja,J2) 
  542. *           Appel a "DGESV" (LAPACK)
  543.             CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO)
  544.             DO M = 1,NT
  545.                Q1(M) = Q1(M)-RHS(M)
  546.             ENDDO
  547.             OMEG = OMEG-RHS(NT+1)
  548.             CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.false.) 
  549.            ENDIF
  550.  
  551. *    ----- CAS 6: CONTINUATION (systeme Hamiltonien) -----------------
  552.            IF (TYPC.EQ.'CN') THEN  
  553. *          Calcul de la correction dY, telle que: Jaa*dY = -[R;g;0]
  554. *          avec Jaa la Jacobienne augmentee: 
  555. *                       Jaa ={RX,Rw,Ra;gx,0,0;dX,dw,0}
  556. *          Ici, 'a' est un parametre artificiel qui ferme le probleme
  557. *          initialement sous-contraint.
  558.             DO II=1,NT
  559.               DO JJ = 1,NT
  560.                 Jaa(II,JJ) = RX(II,JJ)
  561.               ENDDO
  562.             ENDDO
  563. *           Calcul de la derivee: Rw = dR/dw 
  564.             CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.false.)
  565. *           Calcul de la derivee : Ra= (LoI)*X
  566.             CALL HBMDVEC(NT,NHBM,NDDL,Q1,1.D0,Ra)  
  567.             DO KK = 1,NT
  568.                 Jaa(KK,NT+1) = Rw(KK)
  569.                 Jaa(KK,NT+2) = Ra(KK)
  570.                 Jaa(NT+2,KK) = dX(KK)
  571.                 RHS2(KK) = fvec(KK)    
  572.             ENDDO
  573.             Jaa(NT+2,NT+1) = dw
  574.             Jaa(NT+2,NT+2) = ZERO
  575. *           Condition de phase 
  576. *           On minimise la variation de phase entre deux pas succesifs
  577. *           (condition integrale)
  578. *           dg/dw
  579.             Jaa(NT+1,NT+1) = ZERO
  580. *           dg/da
  581.             Jaa(NT+1,NT+2) = ZERO
  582.             CALL HBMDVEC(NT,NHBM,NDDL,QOLD,-1.D0,gx)  
  583.             DO JJ = 1,NT
  584.               Jaa(NT+1,JJ) = gx(JJ)
  585.             ENDDO
  586.             RHS2(NT+1)  = DDOT(NT,gx,1,Q1,1)
  587.             RHS2(NT+2)  = ZERO
  588. *            RHS2(NT+1) = ZERO
  589. *            RHS2(NT+2) = ZERO
  590. *            WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)
  591. *            STOP
  592.             CALL COPYMAT(NT+2,Jaa,Ja2) 
  593.             res=dnrm2(NT+2,RHS2,1) 
  594. *           Appel a "DGESV" (LAPACK)
  595.             CALL DGESV(NT+2,1,Ja2,NT+2,IPIV3,RHS2,NT+2,INFO)
  596. c           si on souhaite verifier la convergence en correction dX
  597. *            res=dnrm2(NT+2,RHS2,1) 
  598. c             WRITE(*,*) 'res=',res
  599. *            WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)        
  600.             DO M = 1,NT
  601.                Q1(M) = Q1(M)-RHS2(M)
  602.             ENDDO
  603.             OMEG = OMEG-RHS2(NT+1)
  604.             CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.false.)    
  605.            ENDIF
  606.  
  607. *    ----- EVALUATION DU RESIDU -----------------------------------
  608.            CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  609.      &           MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,JANAL2,fvec,fref,NNM) 
  610.  
  611. *          residu relatif (calcule autrement dans le cas NNM)
  612.            IF (TYPC(2:2).NE.'N') THEN
  613.              res=dnrm2(NT,fvec,1)/fref0
  614.            ENDIF
  615.  
  616.            IF (IIMPI.GE.2) THEN
  617.              Q1NRM2=dnrm2(NT,Q1,1)
  618. c              IF (IIMPI.GE.333) THEN 
  619. c              WRITE(IOIMP,111) 'Q^',its,(Q1(iou),iou=1,NT)
  620. c              WRITE(IOIMP,111) 'R^',its,(fvec(iou),iou=1,NT)
  621. c              ENDIF
  622.              WRITE(IOIMP,999) its,OMEG,Q1NRM2,res
  623.            ENDIF
  624.  
  625.  
  626. *    ----- TEST DE CONVERGENCE ------------------------------------
  627.            ITER = its
  628.            IF (res.LE.TOLMIN) THEN             
  629.              check=.false.
  630. c            Calcul de la Jacobienne pour la solution convergee
  631.              CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.not.NNM)
  632.              DO J = 1,IPALA(/1)
  633.                IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  634.                  VCTCS(4) = VV
  635.                  XPALA(J,4) = VV
  636.                  CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)
  637.                ENDIF
  638.              ENDDO             
  639.              DO J=1,NT
  640.              DO I=1,NT
  641.                 RX(I,J) = ZZ(I,J) - JAC(I,J)
  642.              ENDDO         
  643.              ENDDO                    
  644.              KTQ = MTQ
  645.              KTEMP=MTEMP
  646.              RETURN
  647.            ELSE
  648.             check = .true.
  649.            ENDIF
  650.  
  651.         ENDDO 
  652. *====== FIN DES ITERATIONS =============================================
  653.  
  654.  
  655. *=======================================================================
  656. *       Mise en forme des messages
  657. *=======================================================================  
  658.  990    FORMAT(' + newt : ',I3,'  w=',F9.5,' |Q|=',F9.5,' |R|=',E10.3,
  659.      &         ' |ref|=',E10.3,' ds=',F9.5)
  660.  999    FORMAT(' + newt : ',I3,'  w=',F9.5,' |Q|=',F9.5,' |R|=',E10.3)
  661. cdebug
  662. c  123    FORMAT('Ja(',I2,',1:8)=',6(E10.3E1,1X),'... F=',E9.3E1)
  663.  123    FORMAT('Ja(',I2,',1:6)=',6(E10.3,1X),'... F=',E10.3)
  664. c  111    FORMAT(A,I2,'=',11(1X,F10.6))
  665.  
  666.        END               
  667.  
  668.  
  669.  
  670.  
  671.  
  672.  
  673.  

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