diff options
Diffstat (limited to 'Dragon/src/SALTLS.f90')
| -rw-r--r-- | Dragon/src/SALTLS.f90 | 299 |
1 files changed, 299 insertions, 0 deletions
diff --git a/Dragon/src/SALTLS.f90 b/Dragon/src/SALTLS.f90 new file mode 100644 index 0000000..32a4500 --- /dev/null +++ b/Dragon/src/SALTLS.f90 @@ -0,0 +1,299 @@ +! +!----------------------------------------------------------------------- +! +!Purpose: +! To generate the standard tracking lines (isotropic tracking) for a +! geometry using the SALT algorithm. +! +!Copyright: +! Copyright (C) 2014 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 data structure in 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. +! NFREG number of regions. +! NBANGL number of angles. +! NQUAD number of quarter (in 2-D). +! RENO track normalisation option. A value RENO=-1 implies +! a direction dependent normalization of the tracks for +! the volume. A value reno=0, implies a global normalisation. +! NBDR number of directions for track normalization (no normalization +! when reno=1). +! 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. +! GG geometry basic information. +! +!Parameters: output +! NBTDIR number of tracks directions considered. +! MAXSGL maximum number of segments in a line. +! NTLINE total number of lines generated. +! DVNOR ratio of analytic to tracked volume. +! +!----------------------------------------------------------------------- +! +SUBROUTINE SALTLS(IFTEMP,IPRINT,IGTRK,NFREG,NBANGL,NQUAD,RENO,NBDR,IFMT,DENUSR, & + DANGLT,DDENWT,GG,NBTDIR,MAXSGL,NTLINE,DVNOR) + USE PRECISION_AND_KINDS, ONLY : PDB,SMALL,PI,TWOPI,HALFPI,INFINITY + USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC + USE SAL_TRACKING_TYPES, ONLY : NMAX2,MINLEN,NNN,ITRAC2,RTRAC2,DPIECE,CNT,CNT0,NBTRAC,LGMORE, & + LGMORE,IERR,DD0,EX0,EY0,DELX,AX,AY,DINIT,ANGTAB,ELMTAB,NB_MAX,IMPX + USE SAL_AUX_MOD, ONLY : SAL231,SAL232,SAL235 + USE SAL_TRAJECTORY_MOD, ONLY : SALTRA + IMPLICIT NONE + !---- + ! subroutine arguments + !---- + INTEGER :: IFTEMP,IPRINT,IGTRK,NFREG,NBANGL,NQUAD,RENO,NBDR,IFMT,MAXSGL,NBTDIR,NTLINE,OK + REAL(PDB) :: DENUSR,DANGLT(2,NQUAD*NBANGL),DDENWT(NQUAD,NBANGL),DVNOR(NFREG,NBDR) + TYPE(T_G_BASIC) :: GG + !---- + ! local parameters + !---- + INTEGER :: II0,IPHI,MQ,MQUAD,NPIECE,INODE,IANGL,II,IQUAD,FIRST,ICURR, & + JCURR,LASTI,NTSEG,IAVERR + LOGICAL :: LGON + REAL :: EPS0,X + REAL(PDB) :: WT,WR,DEL0,DEL2,DNEW,DOLD,ANGLE,DELM,KEEP_DELX, & + MOV_DELX,SUMM,NORM,XFACT,EPSILON_PDB,DCERR,DAVERR, & + DMVERR,DSVERR,TORIG(2) + REAL(PDB), DIMENSION(2) :: THETA0 + REAL(PDB), ALLOCATABLE, DIMENSION(:) :: VOLN,SURFN,CURRN + REAL(PDB), ALLOCATABLE, DIMENSION(:,:) :: FACNRM ! aux for normalisation + REAL, PARAMETER :: EPS3 = 1E-3 + INTEGER, PARAMETER :: FOUT =6 + ! + IMPX=IPRINT + !---- + ! recompute weights + !---- + NBTDIR=0 + SUMM=0._PDB + DO IANGL=1,NBANGL + DO IQUAD=1,NQUAD + IF(DDENWT(IQUAD,IANGL).GT.0._PDB) NBTDIR=NBTDIR+1 + SUMM=SUMM+DDENWT(IQUAD,IANGL) + ENDDO + ENDDO + NORM=0.5_PDB/SUMM + NB_MAX=1 + ALLOCATE(ANGTAB(2),ELMTAB(2)) + !---- + ! isotropic tracking loop + !---- + ! define quadrature + ! get limit intervals for radial quadrature + ! minimum radial interval to contain one trajectory + EPSILON_PDB=SQRT(EPSILON(X)) + EPS0=REAL(10.*EPSILON_PDB) + DEL0=REAL(EPS3,PDB)/DENUSR + DEL2=-1._PDB + ! start a set of tracks + ! CNT0 = address of beginning of trajectory - 1 + ! CNT = address for trajectory data - 1 + ! NBTRAC = number of trajectories in a record + NBTRAC=0 + MAXSGL=0 + IPHI=0 + ALLOCATE(VOLN(GG%NB_NODE),SURFN(GG%NB_SURF2), & + CURRN(GG%NB_SURF2),FACNRM(GG%NB_NODE,NBANGL*NQUAD),STAT=OK) + IF(OK/=0) CALL XABORT('SALTLS: NOT ENOUGH MEMORY R') + FACNRM(:GG%NB_NODE,:NBANGL*NQUAD)=0._PDB + DO IANGL=1,NBANGL + DO IQUAD=1,NQUAD + IPHI=IPHI+1 + ! KEEP IPHI + EX0=DANGLT(1,(IANGL-1)*NQUAD+IQUAD) + EY0=DANGLT(2,(IANGL-1)*NQUAD+IQUAD) + ANGLE=DACOS(EX0) + ! get theta- and theta+ from angle (to be used in SAL235 to + ! decide whether to include projections of tangents to arcs): + THETA0(1)=ANGLE+HALFPI + THETA0(2)=THETA0(1)+PI + IF(THETA0(2)>TWOPI)THETA0(2)=THETA0(2)-TWOPI + ! get projection of points onto axis orthogonal to tracking: + ! only first and last points from macro perimeter + CALL SAL235(NPIECE,THETA0,EX0,EY0,GG%IPAR,GG%RPAR,GG%PERIM_MAC2,GG%NPERIM_MAC2) + ! integrate on each piece + DOLD=DPIECE(1) + DNEW=DPIECE(2) + DELM=DNEW-DOLD + IF(DELM>DEL0)THEN + ! compute nber of intervals for step =< 1/denusr + MQUAD=1+INT(DELM*DENUSR) + DELM=DELM/MQUAD + DO MQ=1,MQUAD + DELX=DOLD+0.5_PDB*DELM + KEEP_DELX=DELX + WR=REAL(DELM) + ! initialize entering distance + DD0=-INFINITY + LGON=.TRUE. + DO WHILE (LGON) + CNT0=0 + CNT=CNT0+NNN + IERR=0 + MOV_DELX=0. + DO WHILE (IERR==0) + ! compute one trajectory: + AX=DELX*EY0 + AY=-DELX*EX0 + CALL SALTRA(DANGLT,GG%NPERIM_MAC2,GG%PERIM_MAC2,GG%ISURF2_ELEM,GG%IPAR, & + GG%RPAR,GG%PPERIM_NODE,GG%PERIM_NODE) + IF(IERR==0)THEN + ! if IERR=0,trajectory has entered into the element joint, moves + ! DELX -> DELX+EPSILON_PDB*N + MOV_DELX=MOV_DELX+EPS0 + DELX=KEEP_DELX+MOV_DELX + CNT=CNT0+NNN + ENDIF + ENDDO + ! a trajectory has been completed: store angle order nber, total weight and wr + !!! corrected by A. Hebert in ev3874 + !!! ITRAC2(CNT0+4)=IPHI + ITRAC2(CNT0+7)=IPHI + RTRAC2(CNT0+7)=DDENWT(IQUAD,IANGL)*WR*NORM + 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('SALTLS: INVALID SEGMENT LENGTH') + ENDDO + ! compute volumes + CALL SAL232(ITRAC2(II0:),RTRAC2(II0:),FACNRM,GG,SURFN,CURRN) + ! next line + IF(CNT+NNN+MINLEN>=NMAX2) CALL XABORT('SALTLS: BUFFER OVERFLOW') + NBTRAC=NBTRAC+1 + WT=RTRAC2(II0+7-1) + WR=RTRAC2(II0+8-1) + FIRST=1 + LASTI=ITRAC2(II0) + IF(GG%NB_SURF2/=0)THEN + ICURR=ITRAC2(II0+5-1) + JCURR=ITRAC2(II0+6-1) + ELSE + ICURR=0 + JCURR=0 + ENDIF + NTSEG=LASTI-FIRST+3 + MAXSGL=MAX(MAXSGL,NTSEG) + IF(IGTRK == 1) THEN + IF(IFMT == 1) THEN + TORIG(1)=AX+DANGLT(1,IPHI)*DINIT ; TORIG(2)=AY+DANGLT(2,IPHI)*DINIT ; + WRITE(IFTEMP) 1,NTSEG,WT,IPHI, & + -ICURR,(ITRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),-JCURR, & + 0.5D0,(RTRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),0.5D0, & + NBTRAC,1,MQ,1,TORIG(1),TORIG(2) + ELSE + WRITE(IFTEMP) 1,NTSEG,WT,IPHI, & + -ICURR,(ITRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),-JCURR, & + 0.5D0,(RTRAC2(II0+II-1),II=FIRST+NNN,LASTI+NNN),0.5D0 + ENDIF + IF(IPRINT>5) WRITE(FOUT,'(2X,''TRAJ# DELX '',''IERR = '',I6,3X,1P,E12.4,I5)') & + NBTRAC,DELX,IERR + ENDIF + LGON=LGMORE + ENDDO + ! end of one interval: move into beginning of next + DOLD=DOLD+DELM + ENDDO + ! end of one piece + ENDIF + ! end of trajectories for this angle + IF(IPRINT > 5) THEN + WRITE(FOUT,'('' ANGLE EX EY NTRA = '',1P,3E12.4,I8,/)') ANGLE,EX0,EY0,NBTRAC + 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 + DO IQUAD=1,NQUAD + VOLN(INODE)=VOLN(INODE)+2._PDB*FACNRM(INODE,IANGL+(IQUAD-1)*NBANGL)* & + DDENWT(IQUAD,IANGL)*NORM + ENDDO + ENDDO + ENDDO + DMVERR=0.0D0 + DSVERR=0.0D0 + DAVERR=0.0D0 + IAVERR=0 + DO INODE=1,GG%NB_NODE + DO IPHI=1,NBANGL*NQUAD + 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,NBANGL*NQUAD + WRITE(*,*) 'iphi=',IPHI + WRITE(*,*) 'facnrm : ',(FACNRM(INODE,IPHI),INODE=1,MIN(10,GG%NB_NODE)) + ENDDO + ENDIF + IF(GG%NB_NODE > NFREG) CALL XABORT('SALTLS: 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,NBANGL*NQUAD + XFACT=FACNRM(INODE,IPHI) + DVNOR(INODE,IPHI+1)=DVNOR(INODE,IPHI+1)+XFACT*GG%VOL_NODE(INODE) + 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 + DEALLOCATE(FACNRM,CURRN,SURFN,VOLN) + DEALLOCATE(ELMTAB,ANGTAB) + RETURN + 6005 FORMAT(' SALTLS: Global RMS, maximum and average errors (%) ', & + 'on region volumes :',3(2X,F10.5)) +END SUBROUTINE SALTLS |