diff options
Diffstat (limited to 'Dragon/src/SALTLC.f90')
| -rw-r--r-- | Dragon/src/SALTLC.f90 | 360 |
1 files changed, 360 insertions, 0 deletions
diff --git a/Dragon/src/SALTLC.f90 b/Dragon/src/SALTLC.f90 new file mode 100644 index 0000000..1e6023e --- /dev/null +++ b/Dragon/src/SALTLC.f90 @@ -0,0 +1,360 @@ +! +!----------------------------------------------------------------------- +! +!Purpose: +! To generate the cyclic tracking lines (specular tracking) for a +! geometry using the SALT algorithm. +! +!Copyright: +! Copyright (C) 2015 Ecole Polytechnique de Montreal. +! This library is free software; you can redistribute it and/or +! modify it under the terms of the GNU Lesser General Public +! License as published by the Free Software Foundation; either +! version 2.1 of the License, or (at your option) any later version +! +!Author(s): +! A. Hebert +! +!Parameters: input +! IFTEMP pointer to a temporary TRACKING file in update or creation +! mode. +! IPRINT print level. +! IGTRK flag to generate the tracking file. In the case where +! IGTRK=1, the tracking is performed and +! used to evaluate the track normalisation factor and the +! tracking file is generated. When IGTRK=0, the tracking is +! still performed and used to evaluate the +! track normalisation factor but the tracking file is not +! generated. +! NDIM problem dimensions. +! NFREG number of regions. +! NBANGL number of angles. +! RENO track normalisation option. A value RENO=-1 implies +! a direction dependent normalization of the tracks +! for the volume while a value RENO=0, implies +! a global normalisation. +! NBDR number of directions for track normalization. +! IFMT tracking file format: +! IFMT=0 for short file; +! IFMT=1 long file required by TLM: module. +! DENUSR user defined track density. +! DANGLT angle cosines. +! DDENWT angular density for each angle. +! NBSANG number of subtracks for each angles. +! GG geometry basic information. +! +!Parameters: output +! MAXSUB maximum number of subtracks in a line. +! MAXSGL maximum number of segments in a line. +! NTLINE total number of lines generated. +! DVNOR ratio of analytic to tracked volume. +! +!----------------------------------------------------------------------- +! +SUBROUTINE SALTLC(IFTEMP,IPRINT,IGTRK,NDIM,NFREG,NBANGL,RENO,NBDR,IFMT,DENUSR,DANGLT, & + DDENWT,NBSANG,GG,MAXSUB,MAXSGL,NTLINE,DVNOR) + USE PRECISION_AND_KINDS, ONLY : PDB,INFINITY,SMALL + USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC,ANGGEO,TYPGEO + USE SAL_TRACKING_TYPES, ONLY : NMAX2,MINLEN,NNN,ITRAC2,RTRAC2,DELR,CNT,CNT0,NBTRAC,IERR,DD0,EX0, & + EY0,DELX,DINIT,EX,ANGTAB,ELMTAB,TORIG,DNEW,N_AXIS,NB_TOT,NB_MAX, & + N_AXIS_KEEP,IMPX + USE SAL_AUX_MOD, ONLY : SAL231,SAL232,SAL237,SAL220_1 + USE SAL_TRAJECTORY_MOD, ONLY : SALTRA + IMPLICIT NONE + !---- + ! Subroutine arguments + !---- + INTEGER :: IFTEMP,IPRINT,IGTRK,NDIM,NFREG,NBANGL,MAXSUB,MAXSGL,NTLINE,RENO,NBDR,IFMT, & + NBSANG(5,NBANGL) + REAL(PDB) :: DENUSR,DANGLT(NDIM,4*NBANGL),DDENWT(4*NBANGL),DVNOR(NFREG,NBDR) + TYPE(T_G_BASIC) :: GG + !---- + INTEGER :: AXIS(2),PIECE,NPIECE,IANGL,II0,II,KEEP_NAXIS,NN,MM,IPHI,P1,P2,INODE,OK,ISURF, & + JCURR,IJK1,ITR,ITRS,JPHI,NA,NAOLD,NEST,NSEG,NTRACK,NTSEG,IAVERR,NMAX3, & + ICYCL,NCYCLE + REAL(PDB) :: WR,WT,PROJTAB(6),KEEP_DELX,MOV_DELX,XFACT,ANGLE,COSSURF,DCERR,DAVERR,DMVERR, & + DSVERR,EPS0 + REAL :: X + REAL(PDB), ALLOCATABLE, DIMENSION(:) :: VOLN + REAL(PDB), ALLOCATABLE, DIMENSION(:,:) :: FACNRM ! aux for normalisation + INTEGER, PARAMETER :: FOUT =6 + INTEGER, ALLOCATABLE, DIMENSION(:) :: ITRACK_TMP + REAL(PDB), ALLOCATABLE, DIMENSION(:) :: RTRACK_TMP + INTEGER, POINTER, DIMENSION(:) :: ITRAC3 + REAL(PDB), POINTER, DIMENSION(:) :: RTRAC3 + ! + IMPX=IPRINT + CALL SAL220_1(ANGGEO) + ! + DELR=1.0D0/DENUSR + NB_MAX=2*MAXVAL(NBSANG(1,:)+NBSANG(2,:)) + IF((TYPGEO.EQ.7).OR.(TYPGEO.EQ.12)) NB_MAX=2*NB_MAX + ALLOCATE(ANGTAB(2*NB_MAX),ELMTAB(2*NB_MAX),TORIG(2,NB_MAX),N_AXIS_KEEP(NB_MAX),STAT=OK) + IF(OK.NE.0) CALL XABORT('SALTLC: Not enough memory ird') + MAXSUB=0 + MAXSGL=0 + NBTRAC=0 + CNT0=0 + CNT=CNT0+NNN + ALLOCATE(VOLN(GG%NB_NODE),FACNRM(GG%NB_NODE,2*NBANGL),STAT=OK) + IF(OK/=0) CALL XABORT('SALTLC: NOT ENOUGH MEMORY R') + FACNRM(:GG%NB_NODE,:2*NBANGL)=0._PDB + ! initialize entering distance + DD0=-INFINITY + EPS0=REAL(10.*SQRT(EPSILON(X))) + NCYCLE=1 + IF((TYPGEO == 5).OR.((TYPGEO == 9))) NCYCLE=2 + DO ICYCL=1,NCYCLE + DO IANGL=1,NBANGL + ! keep IANGL + NN=NBSANG(2,IANGL) ; MM=NBSANG(1,IANGL) + IF(ICYCL == 1) THEN + EX0=DANGLT(1,IANGL) ; EY0=DANGLT(2,IANGL) + ELSE IF(ICYCL == 2) THEN + EX0=-DANGLT(1,IANGL) ; EY0=DANGLT(2,IANGL) + ELSE IF(ICYCL == 3) THEN + EX0=-DANGLT(1,IANGL) ; EY0=-DANGLT(2,IANGL) + ELSE IF(ICYCL == 4) THEN + EX0=DANGLT(1,IANGL) ; EY0=-DANGLT(2,IANGL) + ENDIF + IF(EX0 < 0.0) MM=-MM + IF(EY0 < 0.0) NN=-NN + ! projections of geometry outline onto the two rotation or symmetry axis + P1=GG%PPERIM_MAC2(3); P2=GG%PPERIM_MAC2(4)-1 + CALL SAL237(EX0,EY0,MM,NN,PROJTAB,AXIS) + ! track in direction IANGL and its relative directions + NPIECE=AXIS(1) + ! tracking vector + DO PIECE=1,NPIECE + ! axis nber, position on the axis + N_AXIS=AXIS(2) + DELX=PROJTAB(3)+PROJTAB(5)*(REAL(PIECE)-0.5) + KEEP_NAXIS=N_AXIS + KEEP_DELX=DELX + ! radial weight of the track + WR=PROJTAB(6) + IERR=0 + MOV_DELX=0. + DO WHILE (IERR==0) + ! compute one trajectory: + ! (1) compute one trajectory + ! (2) keep entering points if we have more trajectories + CALL SALTRA(DANGLT(:,:2*NBANGL),GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%ISURF2_ELEM,GG%IPAR, & + GG%RPAR,GG%PPERIM_NODE,GG%PERIM_NODE,GG%IBC2_ELEM,GG%IDATA_BC2,GG%BCDATA, & + GG%PPERIM_MAC2,GG%DIST_AXIS) + IF(IERR==1) DINIT=DNEW + IF(IERR==0) THEN + ! if ierr=0,trajectory has entered into the seam of element joint, + ! moves delx -> delx+epsilon_pdb*n + MOV_DELX=MOV_DELX+EPS0 + DELX=KEEP_DELX+MOV_DELX + N_AXIS=KEEP_NAXIS + CNT=CNT0+NNN + ENDIF + ENDDO + ! if ierr=-1,tracking have not re-entering point=no trajectory + IF(IERR/=-1) THEN + ! a trajectory has been completed: store angle order nber, total weight and WR + ANGLE=DACOS(EX) + ITRAC2(CNT0+4)=IANGL + RTRAC2(CNT0+7)=0.5_PDB*WR/DDENWT(IANGL) + RTRAC2(CNT0+8)=WR + II0=CNT0+1 + IF(IPRINT > 3) CALL SAL231(RTRAC2(II0:),ITRAC2(II0:),DELX,EX0,EY0,ANGLE) + DO II=1,ITRAC2(II0) + IF(RTRAC2(II0+II+NNN-1) <= 0.0) CALL XABORT('SALTLC: INVALID SEGMENT LENGTH') + ENDDO + ! compute volumes + CALL SAL232(ITRAC2(II0:),RTRAC2(II0:),FACNRM,GG) + ! next line + IF(CNT+NNN+MINLEN>=NMAX2) THEN + NMAX3=CNT+NNN+MINLEN+1000 + ALLOCATE(ITRAC3(2*NMAX3),RTRAC3(NMAX3),STAT=OK) + IF(OK/=0) CALL XABORT('SALTLC: NMAX2 overflow.') + RTRAC3(:NMAX2)=RTRAC2(:NMAX2) + ITRAC3(:2*NMAX2)=ITRAC2(:2*NMAX2) + DEALLOCATE(RTRAC2,ITRAC2) + RTRAC2=>RTRAC3 + ITRAC2=>ITRAC3 + NMAX2=CNT+NNN+MINLEN+1000 + ENDIF + NBTRAC=NBTRAC+1 + ! + ! total weight and space weight + COSSURF=RTRAC2(II0+1) + WT=RTRAC2(II0-1+7) + NTRACK=ITRAC2(II0) + NB_TOT=ITRAC2(II0+1) + IF(NB_TOT > NB_MAX) CALL XABORT('SALTLC: NB_MAX overflow.') + ! + ! identify entering end leaving surfaces + JCURR=0 + IJK1=NNN+NTRACK + NTSEG=NTRACK+2*NB_TOT + ! + MAXSGL=MAX(MAXSGL,NTSEG) + MAXSUB=MAX(MAXSUB,NB_TOT) + ALLOCATE(ITRACK_TMP(NTSEG),RTRACK_TMP(NTSEG)) + ! + ! loop over sub-trajectories + NSEG=0 + NAOLD=0 + ITR=NNN + DO ITRS=1,NB_TOT + IJK1=IJK1+1 + NEST=ITRAC2(II0-1+IJK1) + IJK1=IJK1+1 + JPHI=ITRAC2(II0-1+IJK1) + IF(IPRINT > 5) WRITE(FOUT,'(I6,"*",6X,I6,5X,I6)') ITRS,NEST,JPHI + IF(TYPGEO <= 7) NA=(JPHI-1)/NBANGL+1 ! Cartesian + ISURF=N_AXIS_KEEP(ITRS) + IF(ITRS == 1) THEN + ITRACK_TMP(1)=-ISURF + RTRACK_TMP(1)=0.5 + IF(IPRINT > 5) WRITE(FOUT,*) ' -> surface',ISURF + ENDIF + IF(ISURF == 0) THEN + WRITE(FOUT,*) ' JPHI=',JPHI,' -> surface',ISURF + CALL XABORT('SALTLC: symmetry not implemented.') + ENDIF + IF(IPRINT > 5) THEN + IF(TYPGEO <= 7) THEN + WRITE(FOUT,*) ' JPHI=',JPHI,' NAOLD=',NAOLD,' NA=',NA,' -> surface',ISURF + ELSE + WRITE(FOUT,*) ' JPHI=',JPHI,' -> surface',ISURF + ENDIF + ENDIF + NSEG=NSEG+1 + ITRACK_TMP(NSEG)=-ISURF + RTRACK_TMP(NSEG)=0.5 + IF(ITRS > 1) THEN + NSEG=NSEG+1 + ITRACK_TMP(NSEG)=-ISURF + RTRACK_TMP(NSEG)=0.5 + ENDIF + DO II=1,NEST + ITR=ITR+1 + IF(IPRINT > 5) THEN + WRITE(FOUT,'(5X,I6,"*",3(1P,I6,1P,E10.2,I7))') II, & + ITRAC2(II0-1+ITR),RTRAC2(II0-1+ITR),ITRAC2(II0-1+ITR+NMAX2) + ENDIF + NSEG=NSEG+1 + ITRACK_TMP(NSEG)=ITRAC2(II0-1+ITR) + RTRACK_TMP(NSEG)=RTRAC2(II0-1+ITR) + ENDDO + NAOLD=NA + ENDDO + ! case of a geometry with specular reflective condition which is not a + ! rectangular -> anisotropy treatment not supported + IF(JCURR == 0) JCURR=1 + NSEG=NSEG+1 + IF(NSEG /= NTSEG) CALL XABORT('SALTLC: NTSEG inconsistency') + ITRACK_TMP(NSEG)=-JCURR + RTRACK_TMP(NSEG)=0.5 + IF(IPRINT > 5) THEN + WRITE(FOUT,*) 'SALTLC: EXCELT entry with',NSEG,'segments:' + WRITE(FOUT,*) (ITRACK_TMP(II),II=1,NSEG) + WRITE(FOUT,*) (RTRACK_TMP(II),II=1,NSEG) + ENDIF + IF(IGTRK == 1) THEN + IF(IFMT == 1) THEN + WRITE(IFTEMP) NB_TOT,NTSEG,WT,(ANGTAB(II),II=2,2*NB_TOT,2),(ITRACK_TMP(II),II=1,NTSEG), & + (RTRACK_TMP(II),II=1,NTSEG),NBTRAC,1,1,1,((TORIG(II0,II),II0=1,NDIM),II=1,NB_TOT) + ELSE + WRITE(IFTEMP) NB_TOT,NTSEG,WT,(ANGTAB(II),II=2,2*NB_TOT,2),(ITRACK_TMP(II),II=1,NTSEG), & + (RTRACK_TMP(II),II=1,NTSEG) + ENDIF + ENDIF + ! + DEALLOCATE(RTRACK_TMP,ITRACK_TMP) + CNT0=CNT + IF(CNT0+NNN+MINLEN >= NMAX2) THEN + CALL XABORT('SALTLC: NMAX2 overflow(2)') + ELSE + CNT=CNT0+NNN + ENDIF + ENDIF + ENDDO + ! end of trajectories for this angle + IF(IPRINT > 5) THEN + WRITE(FOUT,'('' SALTLC: ANGLE EX EY = '',1P,3E12.4/)') ANGLE,EX0,EY0 + ENDIF + ENDDO + ENDDO + NTLINE=NBTRAC + !---- + ! Compute merged normalization factors + !---- + IF(RENO/=1) THEN + DO INODE=1,GG%NB_NODE + VOLN(INODE)=0._PDB + DO IANGL=1,NBANGL + VOLN(INODE)=VOLN(INODE)+(FACNRM(INODE,IANGL)+FACNRM(INODE,NBANGL+IANGL)) & + /DDENWT(IANGL) + ENDDO + ENDDO + DMVERR=0.0D0 + DSVERR=0.0D0 + DAVERR=0.0D0 + IAVERR=0 + DO INODE=1,GG%NB_NODE + DO IPHI=1,2*NBANGL + IF(RENO==0) THEN + IF(ABS(VOLN(INODE))>SMALL) THEN + FACNRM(INODE,IPHI)=GG%VOL_NODE(INODE)/VOLN(INODE) + ENDIF + ELSE IF(RENO==-1) THEN + IF(ABS(FACNRM(INODE,IPHI))>SMALL) THEN + FACNRM(INODE,IPHI)=GG%VOL_NODE(INODE)/FACNRM(INODE,IPHI) + ENDIF + ENDIF + ENDDO + IF(ABS(VOLN(INODE))>SMALL) THEN + IAVERR=IAVERR+1 + DCERR=100.0D0*(1.0D0-GG%VOL_NODE(INODE)/VOLN(INODE)) + DMVERR=MAX(DMVERR,ABS(DCERR)) + DSVERR=DSVERR+DCERR*DCERR + DAVERR=DAVERR+DCERR + ENDIF + ENDDO + DSVERR=SQRT(DSVERR/DBLE(IAVERR)) + DAVERR=DAVERR/DBLE(IAVERR) + IF(IPRINT > 0) WRITE(FOUT,6005) DSVERR,DMVERR,DAVERR + IF(IPRINT > 5) THEN + DO IPHI=1,2*NBANGL + WRITE(*,*) 'iphi=',IPHI + WRITE(*,*) 'facnrm : ',(FACNRM(INODE,IPHI),INODE=1,MIN(10,GG%NB_NODE)) + ENDDO + ENDIF + IF(GG%NB_NODE > NFREG) CALL XABORT('SALTLC: bug.') + DVNOR(:,:)=0._PDB + DO INODE=1,GG%NB_NODE + IF(RENO==0) THEN + DVNOR(INODE,1)=DVNOR(INODE,1)+FACNRM(INODE,1)*GG%VOL_NODE(INODE) + ELSE IF(RENO==-1) THEN + DO IPHI=1,2*NBANGL + XFACT=FACNRM(INODE,IPHI) + DVNOR(INODE,IPHI+1)=DVNOR(INODE,IPHI+1)+XFACT*GG%VOL_NODE(INODE) + DVNOR(INODE,2*NBANGL+IPHI+1)=DVNOR(INODE,IPHI+1) + ENDDO + ENDIF + ENDDO + DO IPHI=1,NBDR + DVNOR(:,IPHI)=DVNOR(:,IPHI)/GG%VOL_NODE(:) + ENDDO + IF(IPRINT > 4) THEN + DO IPHI=1,NBDR + WRITE(*,*) 'iphi=',IPHI + WRITE(*,*) 'dvnor : ',(DVNOR(INODE,IPHI),INODE=1,MIN(10,NFREG)) + ENDDO + ENDIF + ENDIF + !---- + ! Scratch storage deallocation + !---- + DEALLOCATE(FACNRM,VOLN) + DEALLOCATE(N_AXIS_KEEP,TORIG,ELMTAB,ANGTAB) + RETURN + 6005 FORMAT(' SALTLC: Global RMS, maximum and average errors (%) ', & + 'on region volumes :',3(2X,F10.5)) +END SUBROUTINE SALTLC |