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 LAGNWC C ***************** C ------------------------------------------------------------- & ( IDBIA0 , IDBRA0 , & NDIM , NCELET , NCEL , NFAC , NFABOR , NFML , NPRFML , & NNOD , LNDNOD , LNDFAC , LNDFBR , NCELBR , & NBPMAX , NVP , NVP1 , NVEP , NIVEP , & NPT , NPTNEW , NEW , & IFACEL , IFABOR , IFMFBR , IFMCEL , IPRFML , & IPNFAC , NODFAC , IPNFBR , NODFBR , & ITYCEL , ICOCEL , & IFRLAG , ISORTI , IWORKP , IA , & XYZCEN , SURFAC , SURFBO , CDGFAC , CDGFBO , XYZNOD , VOLUME , & SURFBN , ETTP , RA ) C ------------------------------------------------------------- C C*********************************************************************** C FONCTION : C ---------- c@foncb CFONC CFONC SOUS-PROGRAMME DU MODULE LAGRANGIEN : CFONC ------------------------------------- CFONC CFONC INJECTION EN CONTINUE DES PARTICULES CFONC c@fonce C----------------------------------------------------------------------- C ARGUMENTS c@argub CARGU .______________.____._____.______________________________________. CARGU ! NOM !TYPE!MODE ! ROLE ! CARGU !______________!____!_____!______________________________________! CARGU ! ! ! ! ! CARGU ! IDBIA0 ! E ! -> ! NUMERO DE LA 1ERE CASE LIBRE DANS IA ! CARGU ! IDBRA0 ! E ! -> ! NUMERO DE LA 1ERE CASE LIBRE DANS RA ! CARGU ! NDIM ! E ! -> ! DIMENSION DE L'ESPACE ! CARGU ! NCELET ! E ! -> ! NOMBRE D'ELEMENTS HALO COMPRIS ! CARGU ! NCEL ! E ! -> ! NOMBRE D'ELEMENTS ACTIFS ! CARGU ! NFAC ! E ! -> ! NOMBRE DE FACES INTERNES ! CARGU ! NFABOR ! E ! -> ! NOMBRE DE FACES DE BORD ! CARGU ! NFML ! E ! -> ! NOMBRE DE FAMILLES D ENTITES ! CARGU ! NPRFML ! E ! -> ! NOMBRE DE PROPRIETESE DES FAMILLES ! CARGU ! NNOD ! E ! -> ! NOMBRE DE SOMMETS ! CARGU ! LNDNOD ! E ! -> ! DIM. CONNECTIVITE CELLULES->FACES ! CARGU ! LNDFAC ! E ! -> ! LONGUEUR DU TABLEAU NODFAC ! CARGU ! LNDFBR ! E ! -> ! LONGUEUR DU TABLEAU NODFBR ! CARGU ! NCELBR ! E ! -> ! NOMBRE D'ELEMENTS AYANT AU MOINS UNE ! CARGU ! ! ! ! FACE DE BORD ! CARGU ! NBPMAX ! E ! -> ! NOMBRE MAX DE PARTICULIES AUTORISE ! CARGU ! NVP ! E ! -> ! NOMBRE DE VARIABLES PARTICULAIRES ! CARGU ! NVP1 ! E ! -> ! NVP SANS POSITION, VFLUIDE, VPART ! CARGU ! NVEP ! E ! -> ! NOMBRE INFO PARTICULAIRES (REELS) ! CARGU ! NIVEP ! E ! -> ! NOMBRE INFO PARTICULAIRES (ENTIERS) ! CARGU ! NPT ! E ! <- ! NOMBRE COURANT DE PARTICULES ! CARGU ! NPTNEW ! E ! -> ! NOMBRE TOTAL DE NOUVELLES PARTICULES ! CARGU ! ! ! ! POUR TOUTES LES ZONES D'INJECTION ! CARGU ! NEW ! E ! -> ! NOMBRE DE NOUVELLES PART A INJECTER ! CARGU ! ! ! ! POUR LA ZONE D'INJECTION COURANTE ! CARGU ! IZONE ! E ! -> ! COULEUR DE LA ZONE D'INJECTION ! CARGU ! IFACEL ! TE ! -> ! ELEMENTS VOISINS D'UNE FACE INTERNE ! CARGU ! (2, NFAC) ! ! ! ! CARGU ! IFABOR ! TE ! -> ! ELEMENT VOISIN D'UNE FACE DE BORD ! CARGU ! (NFABOR) ! ! ! ! CARGU ! IFMFBR ! TE ! -> ! NUMERO DE FAMILLE D'UNE FACE DE BORD ! CARGU ! (NFABOR) ! ! ! ! CARGU ! IFMCEL ! TE ! -> ! NUMERO DE FAMILLE D'UNE CELLULE ! CARGU ! (NCELET) ! ! ! ! CARGU ! IPRFML ! TE ! -> ! PROPRIETES D'UNE FAMILLE ! CARGU ! (NFML,NPRFML! ! ! ! CARGU ! IPNFAC ! TE ! -> ! POSITION DU PREMIER NOEUD DE CHAQUE ! CARGU ! (LNDFAC) ! ! ! FACE INTERNE DANS NODFAC ! CARGU ! NODFAC ! TE ! -> ! CONNECTIVITE FACES INTERNES/NOEUDS ! CARGU ! (NFAC+1) ! ! ! ! CARGU ! IPNFBR ! TE ! -> ! POSITION DU PREMIER NOEUD DE CHAQUE ! CARGU ! (LNDFBR) ! ! ! FACE DE BORD DANS NODFBR ! CARGU ! NODFBR ! TE ! -> ! CONNECTIVITE FACES DE BORD/NOEUDS ! CARGU ! (NFABOR+1) ! ! ! ! CARGU ! ICOCEL ! TE ! -> ! CONNECTIVITE CELLULES -> FACES ! CARGU ! (LNDNOD) ! ! ! FACE DE BORD SI NUMERO NEGATIF ! CARGU ! ITYCEL ! TE ! -> ! CONNECTIVITE CELLULES -> FACES ! CARGU ! (NCELET+1) ! ! ! POINTEUR DU TABLEAU ICOCEL ! CARGU ! ISORTI ! TE ! <-> ! POUR CHAQUE PARTICULE : ! CARGU ! (NBPMAX) ! ! ! * NUMERO DE SA CELLULE ! CARGU ! ! ! ! * 0 SI SORTIE DU DOMAINE ! CARGU ! IWORKP(NPTNEW! TE ! -> ! NUMERO DE LA FACE D'INJECTION ! CARGU ! IA(*) ! TR ! - ! MACRO TABLEAU ENTIER ! CARGU ! XYZCEN ! TR ! -> ! POINT ASSOCIES AUX VOLUMES DE CONTROL! CARGU ! (NDIM,NCELET ! ! ! ! CARGU ! SURFAC ! TR ! -> ! VECTEUR SURFACE DES FACES INTERNES ! CARGU ! (NDIM,NFAC) ! ! ! ! CARGU ! SURFBO ! TR ! -> ! VECTEUR SURFACE DES FACES DE BORD ! CARGU ! (NDIM,NFABOR)! ! ! ! CARGU ! CDGFAC ! TR ! -> ! CENTRE DE GRAVITE DES FACES INTERNES ! CARGU ! (NDIM,NFAC) ! ! ! ! CARGU ! CDGFBO ! TR ! -> ! CENTRE DE GRAVITE DES FACES DE BORD ! CARGU ! (NDIM,NFABOR)! ! ! ! CARGU ! XYZNOD ! TR ! -> ! COORDONNES DES NOEUDS ! CARGU ! (NDIM,NNOD) ! ! ! ! CARGU ! VOLUME(NCELET! TR ! -> ! VOLUME D'UN DES NCELET ELEMENTS ! CARGU ! SURFBN(NFABOR! TR ! -> ! SURFACE DES FACES DE BORD ! CARGU ! ETTP ! TR ! -> ! TABLEAUX DES VARIABLES LIEES ! CARGU ! (NBPMAX,NVP)! ! ! AUX PARTICULES ETAPE COURANTE ! 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*********************************************************************** C IMPLICIT NONE C C********************************************************************** C DONNEES EN COMMON C********************************************************************** C INCLUDE "paramx.h" INCLUDE "numvar.h" INCLUDE "optcal.h" INCLUDE "entsor.h" INCLUDE "period.h" INCLUDE "lagpar.h" INCLUDE "lagran.h" C C********************************************************************** C C ARGUMENTS C INTEGER IDBIA0 , IDBRA0 INTEGER NDIM , NCELET , NCEL , NFAC , NFABOR INTEGER NFML , NPRFML INTEGER NNOD , LNDNOD , LNDFAC , LNDFBR , NCELBR INTEGER NBPMAX , NVP , NVP1 , NVEP , NIVEP INTEGER NPT , NPTNEW , NEW INTEGER IFACEL(2,NFAC) , IFABOR(NFABOR) INTEGER IFMFBR(NFABOR) , IFMCEL(NCELET) INTEGER IPRFML(NFML,NPRFML) INTEGER IPNFAC(NFAC+1) , NODFAC(LNDFAC) INTEGER IPNFBR(NFABOR+1) , NODFBR(LNDFBR) INTEGER ICOCEL(LNDNOD) , ITYCEL(NCELET+1) INTEGER ISORTI(NBPMAX) INTEGER IFRLAG(NFABOR) , IWORKP(NPTNEW) INTEGER IA(*) C DOUBLE PRECISION XYZCEN(NDIM,NCELET) DOUBLE PRECISION SURFAC(NDIM,NFAC), SURFBO(NDIM,NFABOR) DOUBLE PRECISION CDGFAC(NDIM,NFAC), CDGFBO(NDIM,NFABOR) DOUBLE PRECISION XYZNOD(NDIM,NNOD), VOLUME(NCELET) DOUBLE PRECISION SURFBN(NFABOR) DOUBLE PRECISION ETTP(NBPMAX,NVP) DOUBLE PRECISION RA(*) C C VARIABLES LOCALES C INTEGER IDEBIA, IDEBRA, IFINIA, IFINRA C INTEGER NP , IEL , N1 , IFAC , KFAC , NBP INTEGER II , JJ , IN , ISORT INTEGER INDIAN , IFAOLD , IFANEW INTEGER ITYPFO , ICONFO(100) INTEGER IDEHOR , IERRIE , ICECPT INTEGER ICELCR, IPERCR ,ITEPAS, ISENS, IPER C DOUBLE PRECISION RD(1) DOUBLE PRECISION XF, YF, ZF DOUBLE PRECISION UP, VP, WP, UF, VF, WF DOUBLE PRECISION PTA(3), PTB(3), VECT(3), VECTN(3) C C======================================================================= C 0. GESTION MEMOIRE C======================================================================= C IDEBIA = IDBIA0 IDEBRA = IDBRA0 C C======================================================================= C 1. INITIALISATION C======================================================================= C C Traitement de la periodicite C IF (IPERIO.EQ.1) THEN C ICELCR = IDEBIA IPERCR = ICELCR + NCELET-NCEL IFINIA = IPERCR + NCELET-NCEL CALL IASIZE('LAGNWC', IFINIA) C =========== C DO IEL = 1,NCELET-NCEL IA(ICELCR+IEL-1) = 0 IA(IPERCR+IEL-1) = 0 ENDDO CALL PERLOC(NCELET, NCEL, IA(ICELCR), IA(IPERCR)) C =========== C ENDIF C C======================================================================= C 2. Injection des particules C======================================================================= C C BOUCLE SUR LES NOUVELLES PARTICULES : C DO NP = 1,NEW C C RAZ du compteur de cellules traversees C ICECPT = 0 C C-->IDEHOR = 1 si exterieur du domaine de calcul C IDEHOR = 0 C C Incrementation du pointeur sur les particules C NPT = NPT + 1 C C On sauvegarde la maille de depart C ISORT = ISORTI(NPT) C C Tirage aleatoire C N1 = 1 CALL ZUFALL(N1,RD) C C Nouvelle position : on le fait dans la C direction de la vitesse mais possibilite C de le faire par rapport a la normale C XF = ETTP(NPT,JXP) + RD(1) *DTP *ETTP(NPT,JUP) YF = ETTP(NPT,JYP) + RD(1) *DTP *ETTP(NPT,JVP) ZF = ETTP(NPT,JZP) + RD(1) *DTP *ETTP(NPT,JWP) C UP = ETTP(NPT,JUP) VP = ETTP(NPT,JVP) WP = ETTP(NPT,JWP) C UF = ETTP(NPT,JUF) VF = ETTP(NPT,JVF) WF = ETTP(NPT,JWF) C IF ( XF.EQ. ETTP(NPT,JXP) & .AND. YF.EQ. ETTP(NPT,JYP) & .AND. ZF.EQ. ETTP(NPT,JZP) ) GOTO 300 C C Est-ce que le point est dans le domaine (Clone de LAGCEL) C C Numero de face d 'injection dans l'element C IFANEW = IWORKP(NPT) C 100 CONTINUE C IEL = ISORT IFAOLD = IFANEW INDIAN = 0 ICECPT = ICECPT + 1 C C ---> Elimination des particules qui posent problemes C la particule reste au niveau de la face d'entree C (boucles infinies) C IF (ICECPT.GT.30) THEN IDEHOR = 1 GOTO 200 ENDIF C C --> balayage des KFAC faces entourant la cellule IEL C Elles sont stockees entre ITYCEL(IEL) et ITYCEL(IEL+1)-1 C (donc KFAC ne peut pas valoir ITYCEL(IEL+1)...) C KFAC = ITYCEL(IEL)-1 C DO WHILE (INDIAN.EQ.0) C KFAC = KFAC + 1 C C --> Erreur : la particule reste au niveau de la face C d'entree C IF (KFAC.EQ.ITYCEL(IEL+1)) THEN IDEHOR = 1 GOTO 200 ENDIF C IFAC = ICOCEL(KFAC) C C-->boucle sur les faces internes C C-->si la face interne a deja ete traitee dans la cellule precedente C son numero est dans IFAOLD et on ne la retraite pas une seconde fois C IF (IFAC.GT.0 .AND. IFAC.NE.IFAOLD) THEN C IN = 0 DO NBP = IPNFAC(IFAC),IPNFAC(IFAC+1)-1 IN = IN + 1 ICONFO(IN) = NODFAC(NBP) ENDDO ITYPFO = IPNFAC(IFAC+1) - IPNFAC(IFAC) + 1 ICONFO(ITYPFO) = ICONFO(1) C CALL OUESTU C =========== & ( NFECRA , NDIM , NNOD , & IERRIE , & ETTP(NPT,JXP) , ETTP(NPT,JYP) , ETTP(NPT,JZP) , & XF , YF , ZF , & CDGFAC(1,IFAC) , CDGFAC(2,IFAC) , CDGFAC(3,IFAC) , & XYZCEN(1,IEL) , XYZCEN(2,IEL) , XYZCEN(3,IEL) , & ITYPFO , ICONFO , XYZNOD , & INDIAN ) C IF (IERRIE.EQ.1) THEN IDEHOR = 1 GOTO 200 C C-->si la particule passe dans la cellule voisine C ELSE IF (INDIAN.EQ.1) THEN C C-->si la particule NPT est dans la cellule II alors le voisin ne C peut etre que la cellule JJ, et vice versa, et inversement, et C ainsi de suite. C IFANEW = IFAC C II = IFACEL(1,IFAC) JJ = IFACEL(2,IFAC) C IF (IEL.EQ.II) THEN ISORT = JJ ELSE IF (IEL.EQ.JJ) THEN ISORT = II ENDIF C C Traitement de la periodicite C IF (ISORT.GT.NCEL) THEN C ITEPAS = ISORT ISORT = IA(ICELCR+ITEPAS-NCEL-1) C C On recupere les informations sur la periodicite C IPER = IA(IPERCR+ITEPAS-NCEL-1)/2-1 IF (MOD(IA(IPERCR+ITEPAS-NCEL-1),2).EQ.0) THEN ISENS = 1 ELSE ISENS = 0 ENDIF C C MODIFICATION DE LA POSITION C PTA(1)= XF PTA(2)= YF PTA(3)= ZF C IEL = ISORT-NCEL C CALL LAGPER(NCELET, NCEL, ISENS, IPER, IEL, PTA, PTB) C =========== C XF = PTB(1) YF = PTB(2) ZF = PTB(3) C C MODIFICATION DE LA VITESSE C VECT(1) = UP VECT(2) = VP VECT(3) = WP C CALL LAGVEC(ISENS, IPER, VECT, VECTN) C =========== C UP = VECTN(1) VP = VECTN(2) WP = VECTN(3) C C MODIFICATION DE LA VITESSE FLUIDE VUE C VECT(1) = UF VECT(2) = VF VECT(3) = WF C CALL LAGVEC(ISENS, IPER, VECT, VECTN) C =========== C UF = VECTN(1) VF = VECTN(2) WF = VECTN(3) C IFANEW = 0 C ENDIF C C--> Retour pour balayage des face de la cellule suivante C GOTO 100 C ENDIF C C--> Balayage des faces de bord (reperees par leur valeur negative C dans ICOCEL) C C resultat : INDIAN = 0 le rayon PQ ne sort pas de la cellule par C ~~~~~~~~ cette face C INDIAN = -1 meme cellule C INDIAN = 1 interaction avec la frontiere C ELSE IF (IFAC.LT.0 .AND. IFAC.NE.IFAOLD) THEN C IFAC = -IFAC C IN = 0 DO NBP = IPNFBR(IFAC),IPNFBR(IFAC+1)-1 IN = IN + 1 ICONFO(IN) = NODFBR(NBP) ENDDO ITYPFO = IPNFBR(IFAC+1) - IPNFBR(IFAC) + 1 ICONFO(ITYPFO) = ICONFO(1) C CALL OUESTU C =========== & ( NFECRA , NDIM , NNOD , & IERRIE , & ETTP(NPT,JXP) , ETTP(NPT,JYP) , ETTP(NPT,JZP) , & XF , YF , ZF , & CDGFBO(1,IFAC) , CDGFBO(2,IFAC) , CDGFBO(3,IFAC) , & XYZCEN(1,IEL) , XYZCEN(2,IEL) , XYZCEN(3,IEL) , & ITYPFO , ICONFO , XYZNOD , & INDIAN ) C C-->si la trajectoire de la particule traverse la face de bord C alors on laisse la particule au niveau de la face d'entree C IF (IERRIE.EQ.1 .OR. INDIAN.EQ.1) THEN IDEHOR = 1 GOTO 200 ENDIF C ENDIF C C fin de DO WHILE ENDDO C C-->Fin de la boucle principale sur les particules C 200 CONTINUE C C Si le point est dans le domaine, alors on injecte la C particule du point sinon on laisse la particule au niveau C de l'entree C IF (IDEHOR.EQ.0) THEN ETTP(NPT,JXP) = XF ETTP(NPT,JYP) = YF ETTP(NPT,JZP) = ZF ISORTI(NPT) = ISORT ENDIF C 300 CONTINUE C ENDDO C C********************************************************************** C C---- C FIN C---- C RETURN C END c@z