c*********************************************************************** SUBROUTINE coll_qs(dt,compt,reprend,err,vare,qmax,corr) c routine private du module mod_static c résolution du système des équations de l'équilibre quasi-statique c Pour le probleme aux limites: c résolution du système d'équations différentielles non linéaires c de la structure interne par développement sur B-splines c par itération Newton avec dérivées analytiques c les équations de la SI sont écrites pour tout point dans c static_m/_r, le jacobien de la méthode NR est formé et résolu c dans resout dans l'espace des splines c pour le modèle quasi-statique, au temps t c paramètres indépendants du temps c m_qs : ordre des splines pour l'intégration par collocation c ord_qs=m_qs+r_qs=m_qs+1 : ordre de la représentation spline c ne : nombre d'inconnues (6 ou 7) c paramètres dépendants du temps c nombre de couches au temps t : n_qs c knot : nombre de points de table (nodaux), c dim_qs=knot-ord_qs : dimension de l'espace des splines c bp(ne,dim_qs) : coefficients des splines c q(n_qs), qt(knot) : abscisses, et vecteur nodal c entrées: c compt: compteur du nombre d'itérations c dt : pas temporel c sorties: c err : erreur max c corr : facteur d'amortissement Newton-Raphson c qmax : indice de la couche au max d'erreur c vare : variable du max. erreur c reprend=.TRUE. : il y a un Pb. il faut diminuer dt c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k c------------------------------------------------------------------------- USE mod_donnees, ONLY : d_lum, d_press, d_ray, d_temp, en_m23, Ipg, 1 langue, lisse, m_qs, ne, nom_qs, npt_lisse, ord_qs, pturb, r_qs USE mod_kind USE mod_nuc, ONLY : t_sup USE mod_numerique, ONLY: bsp1dn, bvald, coll, entre, gauss_band, linf, 1 lisse_sum, no_croiss USE mod_variables, ONLY: bp, dim_qs, knot, m_stat, n_qs, q, qt, r_stat, xl IMPLICIT NONE REAL (kind=dp), INTENT(in) :: dt INTEGER, INTENT(in) :: compt REAL (kind=dp), INTENT(out) :: err, corr INTEGER, INTENT(out) :: qmax, vare LOGICAL, INTENT(out) :: reprend REAL (kind=dp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: dl, ds REAL (kind=dp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: a, b, d, temp REAL (kind=dp), SAVE, ALLOCATABLE, DIMENSION(:) :: nrm REAL (kind=dp), DIMENSION(ne,ne,0:r_qs) :: ae REAL (kind=dp), DIMENSION(ne,0:r_qs) :: y REAL (kind=dp), DIMENSION(ne) :: be, dfdq, f, mini, new_max REAL (kind=dp) :: corrp, er, T_max INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: indpc, i_nrm INTEGER, DIMENSION(ne) :: loc INTEGER, SAVE :: lq=1, nbj, ncoll, nl, n_qsp INTEGER :: eq, cq, i, id, ind, indice, ip, j, ligne, spl, var LOGICAL, SAVE :: debug=.FALSE., init=.TRUE. LOGICAL :: inversible c--------------------------------------------------------------------- 2000 FORMAT(8es10.3) c----------------initialisations--------------------------- IF(init)THEN init=.FALSE. c limites des correction NR maximales, relatives pour R, L, M ALLOCATE(i_nrm(ne),nrm(ne)) ; nrm=0.15d0 nrm(1)=d_press ; nrm(2)=d_temp ; nrm(3)=d_ray ; nrm(4)=d_lum i_nrm=NINT(nrm*1.d2) c définitions et allocations diverses ncoll=(n_qs-1)*m_qs !nombre de points de collocation nl=ne*(ncoll+r_qs) !nombre de lignes du jacobien nbj=ne*ord_qs !longueur du bloc du jacobien dim_qs=(n_qs-1)*m_qs+r_qs !dimension de l'espace des splines ALLOCATE(a(nl,nbj),b(1,nl),d(0:r_qs,ord_qs),dl(ne,0:r_qs,ord_qs), 1 ds(ord_qs,0:r_qs,ord_qs),indpc(nl),temp(ne,dim_qs),xcoll(ncoll),xl(ne)) xl(1)=q(1) !r au centre xl(2)=q(1) !l au centre xl(3)=q(1) !m au centre xl(4)=q(n_qs) !Ptot à l'extérieur xl(5)=q(n_qs) !T à l'extérieur xl(6)=q(n_qs) !m à l'extérieur IF(pturb)xl(Ipg)=q(n_qs) !avec Pturb, Pgaz à l'extérieur CALL coll(q,n_qs,m_qs,xcoll) c PRINT*,'xcoll' ; WRITE(*,2000)xcoll ; PAUSE'xcoll' c les points de raccord q sont les entiers 1, 2,..n_qs, c entre deux points de raccord les abscisses des ord_qs c points de collocation ne diffèrent que par un nombre entier c les valeurs des ord_qs B-splines non id. nulles aux m_qs c points de collocation sur [i,i+1] peuvent être tabulés une fois pour toutes DO i=1,m_qs CALL linf(xcoll(i),qt,knot,lq) CALL bvald(xcoll(i),qt,ord_qs,lq,r_qs,d) ; ds(i,:,:)=d c PRINT*,i ; WRITE(*,2000)ds(i,:,:) ENDDO c PAUSE'ds' c idem pour les limites DO i=1,ne CALL linf(xl(i),qt,knot,lq) CALL bvald(xl(i),qt,ord_qs,lq,r_qs-1,d) ; dl(i,:,:)=d c PRINT*,i ; WRITE(*,2000)dl(i,:,:) ENDDO c PAUSE'dl' c d est désormais inutile DEALLOCATE(d) c écritures de quelques paramètres SELECT CASE(langue) CASE('english') WRITE(*,1001)nrm ; WRITE(2,1001)nrm 1001 FORMAT(/,'Collocation method using modified Newton-Raphson method',/, 1 'maxima of relative corrections nrm=',/, 8es10.3) CASE DEFAULT WRITE(*,1)nrm ; WRITE(2,1)nrm 1 FORMAT(/,'Algorithme de collocation avec Newton-Raphson modifié',/, 1 'corrections relatives maximales autorisées, nrm=',/,8es10.3) END SELECT c nombre de couches du modèle précédent n_qsp=n_qs ENDIF c------------------initialisations fin ------------------------------ c si le nombre de couches du modèle quasi-static a changé IF(n_qs /= n_qsp)THEN n_qsp=n_qs c redéfinitions et allocations diverses ncoll=(n_qs-1)*m_qs !nombre de points de collocation nl=ne*(ncoll+r_qs) !nombre de lignes du jacobien dim_qs=(n_qs-1)*m_qs+r_qs !DIMENSION de l'espace des splines DEALLOCATE(a,b,indpc,temp,xcoll) ALLOCATE(a(nl,nbj),b(1,nl),indpc(nl),temp(ne,dim_qs),xcoll(ncoll)) xl(4)=q(n_qs) !Ptot à l'extérieur xl(5)=q(n_qs) !T à l'extérieur xl(6)=q(n_qs) !m à l'extérieur IF(pturb)xl(Ipg)=q(n_qs) !avec Pturb, Pgas à l'extérieur c les nouveaux points de collocation CALL coll(q,n_qs,m_qs,xcoll) ENDIF c WRITE(*,2000)xcoll(1:8) ; PAUSE'coll_qs xcoll' c le jacobien quasi-statique c initialisations c ip est l'indice de la couche ou se trouvent les m_qs c points de collocation, utilisé dans static_m/r pour le c facteur de répartition fac(ip) ligne=0 ; lq=ord_qs ; a=0.d0 ; b=0.d0 ; ip=1 c pour les points de collocation c PRINT*,'ncoll,m_qs',ncoll,m_qs DO i=1,ncoll,m_qs CALL linf(xcoll(i),q,n_qs,ip) c WRITE(*,2000)REAL(i),REAL(ip),xcoll(i),q(ip),fac(ip) DO j=1,m_qs cq=i+j-1 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,xcoll(cq),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 1 dans coll_qs' y(:,0)=f ; y(:,1)=dfdq ; ind=ne*(lq-ord_qs)+1 CALL static(1,cq,0,y,be,ae,compt,dt,reprend,ip) ; IF(reprend)RETURN DO eq=1,ne ligne=ligne+1 ; b(1,ligne)=be(eq); indpc(ligne)=ind DO spl=1,ord_qs DO var=1,ne indice=ne*(spl-1)+var DO id=0,r_qs a(ligne,indice)=a(ligne,indice)+ae(eq,var,id)*ds(j,id,spl) ENDDO !id ENDDO !var ENDDO !spl c WRITE(*,2000)a(ligne,:),b(1,ligne) ; PAUSE'coll_qs ligne' ENDDO !eq ENDDO !j c WRITE(*,2000)a(ligne,:),b(1,ligne) ; PAUSE'coll_qs i' ; PRINT*,'i=',i ENDDO !i c PAUSE'coll_qs avant lim' c pour les limites, comme r_qs=1, il y a ne limites DO i=1,ne CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,xl(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb de limite dans coll_qs' y(:,0)=f c PRINT*,i ; WRITE(*,2000)xl(i) ; WRITE(*,2000)y(:,0) c WRITE(*,2000)y(:,1) ; PAUSE'lim' CALL static(2,cq,i,y,be,ae,compt,dt,reprend,ip) ; IF(reprend)RETURN ligne=ligne+1 ; b(1,ligne)=be(1) ; indpc(ligne)=ne*(lq-ord_qs)+1 DO spl=1,ord_qs DO var=1,ne indice=ne*(spl-1)+var a(ligne,indice)=a(ligne,indice)+ae(1,var,0)*dl(i,0,spl) ENDDO !var ENDDO !spl c PRINT*,ligne ; WRITE(*,2000)a(ligne,:),b(1,ligne) ; PAUSE'a' ENDDO !i c PAUSE'coll_qs, apres lim' c ----------écritures de debug--------------------- c DO i=nl-15,nl c j=26931 c DO i=j,j+5 c PRINT*,'ligne:',i ; WRITE(*,2000)a(i,:),b(1,i) ; PAUSE'coll_qs, a & b' c ENDDO c PAUSE'dernières lignes' c PRINT*,ligne,nl,nbj ; PRINT*,indpc ; PAUSE'avant gauss' c ---------------pour debug adapter les limites et décommenter c debug=.TRUE. IF(debug)THEN ligne=0 ; spl=26932-ne ; var=MIN(spl+ne,nl) DO i=1,ncoll,m_qs DO j=1,m_qs DO eq=1,ne ligne=ligne+1 IF(entre(spl,var,ligne))THEN PRINT*,'interne, b / a, ligne=',ligne,', équation:',eq WRITE(*,2000)b(1,ligne) WRITE(*,2000)a(ligne,:) ENDIF ENDDO ENDDO ENDDO DO i=1,ne ligne=ligne+1 IF(entre(spl,var,ligne))THEN PRINT*,'limites ext b, / a, ligne=',ligne WRITE(*,2000)b(1,ligne) WRITE(*,2000)a(ligne,:) ENDIF ENDDO c PAUSE'coll_qs debug' ENDIF c -----------------debug----------------------- c résolution du système linéaire c PRINT*,'coll_qs nl,nbj',nl,nbj CALL gauss_band(a,b,indpc,nl,nl,nbj,1,inversible) IF(.NOT.inversible)THEN PRINT*,'jacobien non inversible dans coll_qs' ; PAUSE'ARRET' ; STOP ENDIF c DO i=1,ligne,ne c WRITE(*,2000)(b(1,j),j=i,i+ne-1) ; PAUSE c ENDDO c PAUSE'coll_qs, apres gauss_band' c construction du vecteur temporaire des corrections temp=RESHAPE(b,SHAPE(temp)) c PRINT*,'apres RESHAPE' c la nécessité d'une limitation des corrections NR n'est effective dès qu'une c correction dépasse de nrm(var) = 15% (50% pour L) la valeur de la variable concernée c et ce, en au moins un point où la valeur de la variable dépasse mini=5% de son max mini=0.05d0*MAXVAL(bp,dim=2) c les corrections seront multipliées par corr de façon à limiter c les trop grandes variations, sauf pour psi et Pturb et au centre pour T, R, L, M corr=1.d0 ; corrp=corr DO spl=1,dim_qs B1: DO var=1,5 c le test de précision sur L est supprimé c IF(var == 4)CYCLE B1 c au centre pour R, L, M IF(spl == 1 .AND. entre(3,5,var))CYCLE B1 c si la correction est trop petite IF(ABS(bp(var,spl)) < mini(var))CYCLE B1 c amortissement des corrections B3: DO IF(corr*ABS(temp(var,spl)) < nrm(var)*ABS(bp(var,spl)))EXIT B3 corr=corr/2.d0 ENDDO B3 c écritures IF(corr < 1.d0 .AND. corrp /= corr)THEN SELECT CASE(var) CASE(1,2) WRITE(usl_static,5)nom_qs(var),spl,spl/2,mini(var),1.d0/corr,temp(var,spl), 1 nrm(var),nom_qs(var),bp(var,spl),EXP(bp(var,spl)) 5 FORMAT(/,'coefficient NR modifié sur: ',a7,', spline:', i6,', couche:', 1 i6,', mini:',es10.3,/,'1/NRm=',es10.3,', correction =',es10.3,' > ', 2 f5.2,' pour ',a7,'=',es10.3,', P ou T=',es10.3) CASE DEFAULT WRITE(usl_static,2)nom_qs(var),spl,spl/2,mini(var),1.d0/corr,temp(var,spl), 1 i_nrm(var),nom_qs(var),bp(var,spl) 2 FORMAT(/,'coefficient NR modifié sur: ',a,', spline:',i6, 1 ', couche:',i6,', mini:',es10.3,/,'1/NRm=',es10.3,', correction =' 2 es10.3,' > ',i3,'% pour ',a7,'=',es10.3) END SELECT ENDIF corrp=corr ENDDO B1 !var ENDDO !spl IF(corr == 1.d0)WRITE(usl_static,3)corr 3 FORMAT(/,'coefficient NR modifié: 1/NRm=',es10.3) c estimation de la précision, toutes les corrections sont réduites par le facteur 1/corr err=0.d0 ; er=0.d0 ; vare=0 B4: DO spl=1,dim_qs B2: DO var=1,5 c le test de précision sur L est supprimé c IF(var == 4)CYCLE B2 c au centre, on n'est pas exigeant pour T, R, L, M IF(spl == 1 .AND. entre(2,5,var))CYCLE B2 c si la correction est trop petite IF(ABS(bp(var,spl)) < mini(var))CYCLE B2 er=ABS(temp(var,spl)/bp(var,spl)) IF(er > err)THEN qmax=spl ; vare=var ; err=er ENDIF ENDDO B2 !var ENDDO B4 !spl c corrections au maximum de l'erreur WRITE(usl_static,100)qmax,(qmax-ord_qs+1)/m_qs+1 100 FORMAT(/,'erreur max. {Coeff. var. q-stat. / corr. NR} B-spline#',i6, 1 ', couche#',i6) IF(pturb)THEN WRITE(usl_static,101)nom_qs(1:5),bp(1:5,qmax),bp(7,qmax), 1 temp(1:5,qmax),temp(7,qmax) 101 FORMAT(6(1x,a7,2x),/,6es10.3,/,6es10.3) ELSE WRITE(usl_static,103)nom_qs(1:5),bp(1:5,qmax),temp(1:5,qmax) 103 FORMAT(5(1x,a7,2x),/,5es10.3,/,5es10.3) ENDIF c nouvelle solution bp=bp-temp*corr c DO spl=1,10 !dim_qs c WRITE(*,2000)bp(1:ne,spl) !; PAUSE'coll_qs bp' c ENDDO c pour éviter des erreurs d'arrondi au centre, pour R, L, M c bp(3:5,1)=0.d0 !!!!!!!!!!pas efficace à rejeter c indice de la couche ou se trouve le maximum qmax=(qmax-ord_qs+1)/m_qs+1 c T doit être inférieur à t_sup*1.01 new_max=MAXVAL(bp,dim=2) loc=MAXLOC(bp,dim=2) T_max=EXP(new_max(2)) reprend=T_max > t_sup*1.01d0 IF(reprend)THEN WRITE(*,6)loc(2),T_max,t_sup 6 FORMAT(/,'coll_qs, couche:',i4,', T=',es10.3,' > ',es10.3) RETURN ENDIF c--------------------ESSAI--------------------------- IF(lisse)CALL lisse_sum(npt_lisse,ne,bp,q,qt,ord_qs,n_qs,knot) c----------------------------------------------------- c pour l'interpolation inverse m_stat ou r_stat c avec inter, extraction de r_stat et m_stat de la solution c en_m23 = .TRUE. variables m_stat=m^2/3, r_stat=r^2 c en_m23 = .FALSE. variables m_stat=m^1/3, r_stat=r c pour lim_zc il faut les valeurs de r_stat et m_stat de la dernière itération IF(SIZE(r_stat) /= n_qs)THEN DEALLOCATE(m_stat,r_stat) ; ALLOCATE(m_stat(n_qs),r_stat(n_qs)) ENDIF DO i=1,n_qs r_stat(i)=bp(3,m_qs*(i-1)+1) ; m_stat(i)=bp(5,m_qs*(i-1)+1) c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) !test c WRITE(*,2000)f(3),r_stat(i),f(5),m_stat(i) ; PAUSE'test coll_qs' ENDDO RETURN END SUBROUTINE coll_qs