c@a c@versb C----------------------------------------------------------------------- C CVERS Code_Saturne version 1.3 C ------------------------ C C This file is part of the Code_Saturne Kernel, element of the C Code_Saturne CFD tool. C C Copyright (C) 1998-2007 EDF S.A., France C C contact: saturne-support@edf.fr C C The Code_Saturne Kernel is free software; you can redistribute it C and/or modify it under the terms of the GNU General Public License C as published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C The Code_Saturne Kernel is distributed in the hope that it will be C useful, but WITHOUT ANY WARRANTY; without even the implied warranty C of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public License C along with the Code_Saturne Kernel; if not, write to the C Free Software Foundation, Inc., C 51 Franklin St, Fifth Floor, C Boston, MA 02110-1301 USA C C----------------------------------------------------------------------- c@verse SUBROUTINE JACOBI C ***************** C ------------------------------------------------------------- & ( CNOM , NCELET , NCEL , NFAC , & ISYM , NITMAP , IINVPE , IWARNP , NFECRA , NITERF , & EPSILP , RNORM , RESIDU , & IFACEL , IA , & DAM , XAM , SMBRP , VX , RK , RA ) C ------------------------------------------------------------- C*********************************************************************** C FONCTION : C ---------- c@foncb CFONC CFONC RESOLUTION DU SYSTEME (DAM+ XAM).VX = SMBRP PAR UNE METHODE DE CFONC JACOBI OU GAUSS SEIDEL RELAXEE CFONC CFONC VX SUPPOSE INITIALISE EN ENTREE CFONC c@fonce C----------------------------------------------------------------------- C ARGUMENTS c@argub CARGU .______________.____._____.______________________________________. CARGU ! NOM !TYPE!MODE ! ROLE ! CARGU !______________!____!_____!______________________________________! CARGU ! CNOM ! A ! -> ! NOM DE LA VARIABLE ! CARGU ! NCELET ! E ! -> ! NOMBRE D'ELEMENTS HALO COMPRIS ! CARGU ! NCEL ! E ! -> ! NOMBRE D'ELEMENTS ACTIFS ! CARGU ! NFAC ! E ! -> ! NOMBRE DE FACES INTERNES ! CARGU ! ISYM ! E ! -> ! INDICATEUR = 1 MATRICE SYM ! CARGU ! ! ! ! = 2 MATRICE NON SYM ! CARGU ! NITMAP ! E ! -> ! NOMBRE MAX D'ITER POUR RESOL ITERATIV! CARGU ! IINVPE ! E ! -> ! INDICATEUR POUR ANNULER LES INCREMENT! CARGU ! ! ! ! EN PERIODICITE DE ROTATION (=2) OU ! CARGU ! ! ! ! POUR LES ECHANGER NORMALEMENT DE ! CARGU ! ! ! ! MANIERE SCALAIRE (=1) ! CARGU ! IWARNP ! E ! -> ! NIVEAU D'IMPRESSION ! CARGU ! NFECRA ! E ! -> ! UNITE DU FICHIER SORTIE STD ! CARGU ! NITERF ! E ! <- ! NOMBRE D'ITERATIONS EFFECTUEES ! CARGU ! EPSILP ! R ! -> ! PRECISION POUR RESOL ITER ! CARGU ! RNORM ! R ! -> ! NORMALISATION DU RESIDU ! CARGU ! RESIDU ! R ! <- ! RESIDU FINAL NON NORME ! CARGU ! IFACEL(2,NFAC! TE ! -> ! No DES ELTS VOISINS D'UNE FACE INTERN! CARGU ! IA(*) ! TR ! - ! MACRO TABLEAU ENTIER ! CARGU ! DAM(NCELET ! TR ! -> ! DIAGONALE DE LA MATRICE ! CARGU ! XAM(NFAC,ISYM! TR ! -> ! EXTRADIAGONALE DE LA MATRICE ! CARGU ! SMBRP(NCELET ! TR ! -> ! SECOND MEMBRE DU SYSTEME ! CARGU ! VX (NCELET ! TR ! <-> ! SOLUTION DU SYSTEME ! CARGU ! RK (NCELET ! TR ! - ! AUXILIAIRE DE TRAVAIL ! CARGU ! RA(*) ! TR ! - ! MACRO TABLEAU REEL ! CARGU !______________!____!_____!______________________________________! c@argue C c@commb CCOMM COMMONS CCOMM .______________.____._____.______________________________________. CCOMM ! NOM !TYPE!MODE ! ROLE ! CCOMM !______________!____!_____!______________________________________! CCOMM !______________!____!_____!______________________________________! c@comme C C TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU) C L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL) C MODE : -> DONNEE, <- RESULTAT, <-> DONNEE MODIFIEE, C - TABLEAU DE TRAVAIL C*********************************************************************** C IMPLICIT NONE C C*********************************************************************** C DONNEES EN COMMON C*********************************************************************** C INCLUDE "paramx.h" INCLUDE "cstnum.h" INCLUDE "vector.h" INCLUDE "period.h" INCLUDE "parall.h" C C*********************************************************************** C C ARGUMENTS C CHARACTER*8 CNOM INTEGER NCELET , NCEL , NFAC INTEGER ISYM , NITMAP , IINVPE , IWARNP , NFECRA INTEGER NITERF INTEGER IDIMTE , ITENSO DOUBLE PRECISION EPSILP , RNORM , RESIDU C INTEGER IFACEL(2,NFAC) INTEGER IA(*) DOUBLE PRECISION DAM(NCELET ),XAM(NFAC ,ISYM),SMBRP(NCELET) DOUBLE PRECISION VX(NCELET),RK(NCELET) DOUBLE PRECISION RA(*) C C VARIABLES LOCALES C INTEGER IFAC, IEL, II, JJ, ISQRT C C*********************************************************************** C C======================================================================= C 1. CALCULS PRELIMINAIRES C======================================================================= C IF(ISYM.NE.2.AND.ISYM.NE.1) THEN WRITE(NFECRA,9000)CNOM,ISYM CALL CSEXIT (1) ENDIF C NITERF = 0 C ISQRT = 1 CALL PRODSC(NCELET,NCEL,ISQRT,SMBRP,SMBRP,RESIDU) C IF( RNORM.LE.EPZERO .OR. RESIDU.LE.EPZERO ) THEN IF(IWARNP.GE.2) THEN WRITE (NFECRA,1000) CNOM,RNORM,RESIDU ENDIF RETURN ENDIF C C C======================================================================= C 2. ITERATION COURANTE (STD OU MAILLAGE PAIR/IMPAIR) C======================================================================= C 100 CONTINUE NITERF = NITERF +1 C C DO IEL=1,NCEL RK(IEL) = VX(IEL) ENDDO C C ---> TRAITEMENT DU PARALLELISME C IF(IRANGP.GE.0) CALL PARCOM (RK) C =========== C C ---> TRAITEMENT DE LA PERIODICITE C IF(IPERIO.EQ.1) THEN C IF(IINVPE.EQ.1) THEN IDIMTE = 0 ITENSO = 0 CALL PERCOM C =========== & ( IDIMTE , ITENSO , & RK , RK , RK , & RK , RK , RK , & RK , RK , RK ) ELSEIF(IINVPE.EQ.2) THEN IDIMTE = 0 ITENSO = 11 CALL PERCOM C =========== & ( IDIMTE , ITENSO , & RK , RK , RK , & RK , RK , RK , & RK , RK , RK ) ENDIF C ENDIF C C ---> CAS STANDARD C DO IEL = 1, NCEL VX(IEL) = SMBRP(IEL) ENDDO IF(NCELET.GT.NCEL) THEN DO IEL = NCEL+1, NCELET VX(IEL) = 0.D0 ENDDO ENDIF C C C On distingue les deux cas pour permettre au compilateur C d'optimiser (ISYM = 1 et ISYM = 2). C La seule difference entre les deux blocs est l'indice de C XAM(IFAC,.) dans la ligne C VX(JJ) = VX(JJ) -XAM(IFAC,ISYM)*RK(II) C IF(ISYM.EQ.1) THEN C IF (IVECTI.EQ.1) THEN C !OCL NOVREC,VRL(16) DO IFAC = 1, NFAC II = IFACEL(1,IFAC) JJ = IFACEL(2,IFAC) VX(II) = VX(II) -XAM(IFAC,1)*RK(JJ) VX(JJ) = VX(JJ) -XAM(IFAC,1)*RK(II) ENDDO C ELSE C C VECTORISATION NON FORCEE DO IFAC = 1, NFAC II = IFACEL(1,IFAC) JJ = IFACEL(2,IFAC) VX(II) = VX(II) -XAM(IFAC,1)*RK(JJ) VX(JJ) = VX(JJ) -XAM(IFAC,1)*RK(II) ENDDO C ENDIF C ELSE C IF (IVECTI.EQ.1) THEN C !OCL NOVREC,VRL(16) DO IFAC = 1, NFAC II = IFACEL(1,IFAC) JJ = IFACEL(2,IFAC) VX(II) = VX(II) -XAM(IFAC,1)*RK(JJ) VX(JJ) = VX(JJ) -XAM(IFAC,2)*RK(II) ENDDO C ELSE C C VECTORISATION NON FORCEE DO IFAC = 1, NFAC II = IFACEL(1,IFAC) JJ = IFACEL(2,IFAC) VX(II) = VX(II) -XAM(IFAC,1)*RK(JJ) VX(JJ) = VX(JJ) -XAM(IFAC,2)*RK(II) ENDDO C ENDIF C ENDIF C C DO IEL = 1, NCEL VX(IEL) = VX(IEL)/DAM(IEL) ENDDO C C======================================================================= C 3. TEST DE CONVERGENCE ET IMPRESSIONS C======================================================================= C RESIDU = 0.D0 DO IEL=1,NCEL RESIDU = RESIDU +( DAM(IEL)*(VX(IEL)-RK(IEL)) )**2 ENDDO IF (IRANGP.GE.0) CALL PARSOM (RESIDU) C =========== RESIDU = SQRT(RESIDU) C IF( RESIDU.GT.EPSILP*RNORM ) THEN IF( NITERF.LT.NITMAP ) THEN IF(IWARNP.GE.3) THEN WRITE (NFECRA,1100)CNOM,NITERF,RESIDU,RESIDU/RNORM ENDIF GOTO 100 ELSE IF(IWARNP.GE.1) THEN WRITE (NFECRA,1100)CNOM,NITERF,RESIDU,RESIDU/RNORM WRITE (NFECRA,1200)CNOM ENDIF ENDIF ELSE IF(IWARNP.GE.2) THEN WRITE (NFECRA,1100)CNOM,NITERF,RESIDU,RESIDU/RNORM ENDIF ENDIF C C-------- C FORMATS C-------- C 1000 FORMAT(1X,A8, & ' JACOBI: SORTIE IMMEDIATE: RNORM:', E11.4,' RESIDU:',E11.4) 1100 FORMAT (1X,A8, & ' : JACOBI NITER: ',I5,' RES ABS:',E11.4,' RES NOR:',E11.4 ) 1200 FORMAT ( &'@ ',/, &'@ @@ ATTENTION : ',A8 ,' NON CONVERGENCE DE JACOBI ',/, &'@ ********* ',/, &'@ ' ) 9000 FORMAT ( &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/, &'@ @@ ATTENTION : ARRET DANS jacobi ',/, &'@ ********* ',/, &'@ APPEL DE jacobi POUR ',A8 ,' AVEC ISYM = ',I10 ,/, &'@ ',/, &'@ Le calcul ne peut pas etre execute. ',/, &'@ ',/, &'@ Contacter l''assistance. ',/, &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/) C C12345678 : JACOBI NITER: 10000 RES ABS: 1234567890 RES NOR: 124567890 C---- C FIN C---- C RETURN C END c@z