summaryrefslogtreecommitdiff
path: root/Dragon/src/FMTSUD.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/FMTSUD.f')
-rw-r--r--Dragon/src/FMTSUD.f421
1 files changed, 421 insertions, 0 deletions
diff --git a/Dragon/src/FMTSUD.f b/Dragon/src/FMTSUD.f
new file mode 100644
index 0000000..b41b72d
--- /dev/null
+++ b/Dragon/src/FMTSUD.f
@@ -0,0 +1,421 @@
+*DECK FMTSUD
+ SUBROUTINE FMTSUD(NENTRY,IENTRY,KENTRY,IPRINT,NOPT,IOPT,IKFLU,
+ > NREG ,NGROUP,NDIM ,POLOAQ,AZMOAQ,
+ > VOLUME,XPOL ,WPOL ,XAZI ,WAZI ,
+ > FLUX ,AFLUX ,TFLUX ,WGHT ,MU ,ETA )
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* To process the angular fluxes and generate the SUS3D file.
+*
+*Copyright:
+* Copyright (C) 2008 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):
+* G. Marleau
+*
+*Parameters: input
+* NENTRY number of data structures transfered to this module.
+* IENTRY data structure type where:
+* =1 for LCM memory object;
+* =2 for XSM file;
+* =3 for sequential binary file;
+* =4 for sequential ASCII file.
+* KENTRY data structure pointer.
+* IPRINT print level.
+* NOPT number of options.
+* IOPT processing option.
+* IKFLU pointer to the FLUX data structure.
+* NREG number of regions for problem.
+* NGROUP number of groups for problem.
+* NDIM number of dimensions of problem.
+* POLOAQ polar quadrature order.
+* AZMOAQ azimuthal quadrature order.
+* VOLUME regional volumes.
+* XPOL polar quadrature points.
+* WPOL polar quadrature weights.
+* XAZI azimuthal quadrature points.
+* WAZI azimuthal quadrature weights.
+*
+*Parameters: output
+* FLUX direct and adjoint flux.
+* AFLUX angular components of the direct and adjoint flux.
+* TFLUX temporary flux vector.
+* WGHT temporary weight.
+* MU temporary mu.
+* ETA temporary eta.
+*
+*----------
+*
+ USE GANLIB
+ IMPLICIT NONE
+*----
+* Subroutine arguments
+*----
+ TYPE(C_PTR) KENTRY(NENTRY)
+ INTEGER NENTRY,IENTRY(NENTRY)
+ INTEGER IPRINT,NOPT,IOPT(NOPT),IKFLU
+ INTEGER NREG,NGROUP,NDIM,POLOAQ,AZMOAQ
+ REAL VOLUME(NREG),XPOL(POLOAQ,2),WPOL(POLOAQ),
+ > XAZI(NDIM,AZMOAQ),WAZI(AZMOAQ)
+ REAL FLUX(NREG,2,NGROUP,2)
+ DOUBLE PRECISION AFLUX(NREG,POLOAQ,AZMOAQ*2,2,NGROUP)
+ DOUBLE PRECISION TFLUX(NREG+1,POLOAQ,AZMOAQ*2)
+ REAL WGHT(POLOAQ*AZMOAQ*2),MU(POLOAQ*AZMOAQ*2),
+ > ETA(POLOAQ*AZMOAQ*2)
+*----
+* Local parameters
+*----
+ INTEGER IOUT
+ CHARACTER NAMSBR*6
+ PARAMETER (IOUT=6,NAMSBR='FMTSUD')
+ INTEGER ILCMUP,ILCMDN
+ PARAMETER (ILCMUP=1,ILCMDN=2)
+*----
+* Local variables
+*----
+ TYPE(C_PTR) IPU
+ INTEGER IFPU,IAPU,IFORM,IGROUP,IKU,IAKU,NAZ,IA,IP,IR,IQUA
+ CHARACTER*12 NAMREC
+ REAL RELERR,ERRMAX,ETATMP,FNORM
+*----
+* Initialize FLUX vectors
+*----
+ NAZ=AZMOAQ
+ DO IP=1,POLOAQ
+ XPOL(IP,2)=SQRT(1.-XPOL(IP,1)*XPOL(IP,1))
+ ENDDO
+ IF(IPRINT .GE. 10) THEN
+ WRITE(IOUT,6000) NAMSBR,NGROUP,NDIM,2*AZMOAQ,POLOAQ,NREG
+ WRITE(IOUT,6001)
+ WRITE(IOUT,6011)
+ > (XAZI(1,IA),XAZI(2,IA),WAZI(IA),IA=1,AZMOAQ)
+ WRITE(IOUT,6002)
+ WRITE(IOUT,6011)
+ > (XPOL(IP,1),XPOL(IP,2),WPOL(IP),IP=1,POLOAQ)
+ ENDIF
+ FLUX(:NREG,:2,:NGROUP,:2)=0.0
+ AFLUX(:NREG,:POLOAQ,:AZMOAQ*2,:2,:NGROUP)=0.0D0
+*----
+* Get information from FLUX data structure.
+* 1. Flux
+* 2. Angular flux
+* 3. Adjoint
+* 4. Angular adjoint
+*----
+ IPU=KENTRY(IKFLU)
+ CALL LCMSIX(IPU,'FLUXDIRECT ',ILCMUP)
+ DO IGROUP=1,NGROUP
+ WRITE(NAMREC,'(A4,I3,5X)') 'FLUX',IGROUP
+ CALL LCMGET(IPU,NAMREC,FLUX(1,1,IGROUP,1))
+ ENDDO
+ CALL LCMSIX(IPU,'ANGULAR DIR ',ILCMUP)
+ DO IGROUP=1,NGROUP
+ WRITE(NAMREC,'(A5,I3,4X)') 'AFLUX',IGROUP
+ CALL LCMGET(IPU,NAMREC,AFLUX(1,1,1,1,IGROUP))
+ ENDDO
+ CALL LCMSIX(IPU,'ANGULAR DIR ',ILCMDN)
+ CALL LCMSIX(IPU,'FLUXDIRECT ',ILCMDN)
+ CALL LCMSIX(IPU,'FLUXADJOINT ',ILCMUP)
+ DO IGROUP=1,NGROUP
+ WRITE(NAMREC,'(A4,I3,5X)') 'FLUX',IGROUP
+ CALL LCMGET(IPU,NAMREC,FLUX(1,2,IGROUP,1))
+ ENDDO
+ CALL LCMSIX(IPU,'ANGULAR ADJ ',ILCMUP)
+ DO IGROUP=1,NGROUP
+ WRITE(NAMREC,'(A5,I3,4X)') 'AFLUX',IGROUP
+ CALL LCMGET(IPU,NAMREC,AFLUX(1,1,1,2,IGROUP))
+ ENDDO
+ CALL LCMSIX(IPU,'ANGULAR ADJ ',ILCMDN)
+ CALL LCMSIX(IPU,'FLUXADJOINT ',ILCMDN)
+*----
+* Create first SUS file
+* Volume, and tracking directions
+*----
+ IFPU=FILUNIT(KENTRY(1))
+ WRITE(IFPU,1000) NREG
+ WRITE(IFPU,1001) 2*AZMOAQ*POLOAQ
+*----
+* -\Omega (\varphi+\pi,\pi-\theta)
+*----
+ IQUA=0
+ DO IA=1,AZMOAQ
+ ETATMP=XAZI(2,NAZ+1-IA)
+ IF(ETATMP .EQ. 0.0) THEN
+ ETATMP=1.0E-10
+ ENDIF
+ DO IP=1,POLOAQ
+ IQUA=IQUA+1
+ MU(IQUA)=-XPOL(IP,2)*XAZI(1,NAZ+1-IA)
+ ETA(IQUA)=-XPOL(IP,2)*ETATMP
+ WGHT(IQUA)=WPOL(IP)/WAZI(NAZ+1-IA)
+ ENDDO
+ ENDDO
+*----
+* +\Omega (\varphi,\theta)
+*----
+ DO IA=1,AZMOAQ
+ ETATMP=XAZI(2,IA)
+ IF(ETATMP .EQ. 0.0) THEN
+ ETATMP=1.0E-10
+ ENDIF
+ DO IP=1,POLOAQ
+ IQUA=IQUA+1
+ MU(IQUA)=XPOL(IP,2)*XAZI(1,IA)
+ ETA(IQUA)=XPOL(IP,2)*ETATMP
+ WGHT(IQUA)=WPOL(IP)/WAZI(IA)
+ ENDDO
+ ENDDO
+ WRITE(IFPU,1010) 'wght',
+ > (WGHT(IQUA),IQUA=1,POLOAQ*AZMOAQ*2)
+ WRITE(IFPU,1010) 'mu ',
+ > (MU(IQUA),IQUA=1,POLOAQ*AZMOAQ*2)
+ WRITE(IFPU,1010) 'eta ',
+ > (ETA(IQUA),IQUA=1,POLOAQ*AZMOAQ*2)
+ WRITE(IFPU,1010) 'vol ',
+ > (VOLUME(IR),IR=1,NREG)
+*----
+* Print angular flux
+*----
+ IFPU=FILUNIT(KENTRY(2))
+ IKU=IENTRY(2)
+ IFORM=IOPT(2)
+ ERRMAX=0.0
+ IF(IKU .EQ. 3) THEN
+ DO IGROUP=1,NGROUP
+ WRITE(IFPU) 'Group',IGROUP
+ ENDDO
+ DO IGROUP=NGROUP+1,NGROUP+6
+ WRITE(IFPU) 'Comments'
+ ENDDO
+ ELSE
+ DO IGROUP=1,NGROUP
+ WRITE(IFPU,'(A8,4X,I8)') 'Groups =',IGROUP
+ ENDDO
+ DO IGROUP=NGROUP+1,NGROUP+6
+ WRITE(IFPU,'(A8)') 'Comments'
+ ENDDO
+ ENDIF
+ DO IGROUP=1,NGROUP
+ DO IA=1,AZMOAQ*2
+ DO IP=1,POLOAQ
+ IF(IFORM .EQ. 1) THEN
+ DO IR=1,NREG
+ TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,1,IGROUP)
+ IF(IA .LE. AZMOAQ) THEN
+ FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+ ELSE
+ FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+ ENDIF
+ ENDDO
+ IF(IPRINT .GE. 100) THEN
+ WRITE(IOUT,6003) IP,IA,IGROUP
+ WRITE(IOUT,6010) (TFLUX(IR,IP,IA),IR=1,NREG)
+ ENDIF
+ ELSE
+ IR=1
+ TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,1,IGROUP)
+ DO IR=1,NREG
+ TFLUX(IR+1,IP,IA)=
+ > 2.0D0*AFLUX(IR,IP,IA,1,IGROUP)-TFLUX(IR,IP,IA)
+ IF(IA .LE. AZMOAQ) THEN
+ FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+ ELSE
+ FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+ ENDIF
+ ENDDO
+ IF(IPRINT .GE. 100) THEN
+ WRITE(IOUT,6003) IP,IA,IGROUP
+ WRITE(IOUT,6010)
+ > ((TFLUX(IR,IP,IA)+TFLUX(IR+1,IP,IA))/2.0,IR=1,NREG)
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ IF(IFORM .EQ. 1) THEN
+ IF(IKU .EQ. 3) THEN
+ WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ELSE
+ WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ENDIF
+ ELSE
+ IF(IKU .EQ. 3) THEN
+ WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ELSE
+ WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ENDIF
+ ENDIF
+ ENDDO
+ IF(IPRINT .GE. 10) THEN
+ FNORM=FLUX(1,1,1,2)/FLUX(1,1,1,1)
+ DO IGROUP=1,NGROUP
+ WRITE(6,6030) IGROUP
+ DO IR=1,NREG
+ RELERR=100.0*(FLUX(IR,1,IGROUP,2)
+ > -FNORM*FLUX(IR,1,IGROUP,1))
+ > /FLUX(IR,1,IGROUP,2)
+ ERRMAX=MAX(ERRMAX,ABS(RELERR))
+ IF(IPRINT .GE. 20) THEN
+ WRITE(6,6031) IR,FNORM*FLUX(IR,1,IGROUP,1),
+ > FLUX(IR,1,IGROUP,2),RELERR
+ ENDIF
+ ENDDO
+ ENDDO
+ WRITE(6,6020) ERRMAX
+ ENDIF
+*----
+* Print angular adjoint
+*----
+ IAPU=FILUNIT(KENTRY(3))
+ IAKU=IENTRY(3)
+ ERRMAX=0.0
+ IF(IAKU .EQ. 3) THEN
+ DO IGROUP=NGROUP,1,-1
+ WRITE(IAPU) 'Group',IGROUP
+ ENDDO
+ DO IGROUP=NGROUP+1,NGROUP+6
+ WRITE(IAPU) 'Comments'
+ ENDDO
+ ELSE
+ DO IGROUP=NGROUP,1,-1
+ WRITE(IAPU,'(A8,4X,I8)') 'Groups =',IGROUP
+ ENDDO
+ DO IGROUP=NGROUP+1,NGROUP+6
+ WRITE(IAPU,'(A8)') 'Comments'
+ ENDDO
+ ENDIF
+ DO IGROUP=NGROUP,1,-1
+ DO IA=1,AZMOAQ*2
+ DO IP=1,POLOAQ
+ IF(IFORM .EQ. 1) THEN
+ DO IR=1,NREG
+ TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,2,IGROUP)
+ IF(IA .LE. AZMOAQ) THEN
+ FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+ ELSE
+ FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+ ENDIF
+ ENDDO
+ IF(IPRINT .GE. 100) THEN
+ WRITE(IOUT,6004) IP,IA,IGROUP
+ WRITE(IOUT,6010) (TFLUX(IR,IP,IA),IR=1,NREG)
+ ENDIF
+ ELSE
+ IR=1
+ TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,2,IGROUP)
+ DO IR=1,NREG
+ TFLUX(IR+1,IP,IA)=
+ > 2.0D0*AFLUX(IR,IP,IA,2,IGROUP)-TFLUX(IR,IP,IA)
+ IF(IA .LE. AZMOAQ) THEN
+ FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+ ELSE
+ FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+ > +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+ ENDIF
+ ENDDO
+ IF(IPRINT .GE. 100) THEN
+ WRITE(IOUT,6004) IP,IA,IGROUP
+ WRITE(IOUT,6010)
+ > ((TFLUX(IR,IP,IA)+TFLUX(IR+1,IP,IA))/2.0,IR=1,NREG)
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ IF(IFORM .EQ. 1) THEN
+ IF(IAKU .EQ. 3) THEN
+ WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ELSE
+ WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ENDIF
+ ELSE
+ IF(IAKU .EQ. 3) THEN
+ WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ELSE
+ WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=1,AZMOAQ)
+ WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+ > IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+ ENDIF
+ ENDIF
+ ENDDO
+ IF(IPRINT .GE. 10) THEN
+ FNORM=FLUX(1,2,1,2)/FLUX(1,2,1,1)
+ DO IGROUP=NGROUP,1,-1
+ WRITE(6,6032) IGROUP
+ DO IR=1,NREG
+ RELERR=100.0*(FLUX(IR,2,IGROUP,2)
+ > -FNORM*FLUX(IR,2,IGROUP,1))
+ > /FLUX(IR,1,IGROUP,2)
+ ERRMAX=MAX(ERRMAX,ABS(RELERR))
+ IF(IPRINT .GE. 20) THEN
+ WRITE(6,6031) IR,FNORM*FLUX(IR,2,IGROUP,1),
+ > FLUX(IR,2,IGROUP,2),RELERR
+ ENDIF
+ ENDDO
+ ENDDO
+ WRITE(6,6020) ERRMAX
+ ENDIF
+*----
+* Processing finished, return
+*----
+ RETURN
+*----
+* Formats
+*----
+ 1000 FORMAT('Number of regions =',I5)
+ 1001 FORMAT('Nomber of angles =',I5)
+ 1002 FORMAT(1P,5E20.10)
+ 1010 FORMAT(A4,1P/(3E15.7))
+ 6000 FORMAT('Output from routine ',A6/
+ > 'Number of groups =',I5/
+ > 'Number of dimens =',I5/
+ > 'Number of azimuth =',I5/
+ > 'Number of polar =',I5/
+ > 'Number of regions =',I5)
+ 6001 FORMAT('Azimuthal quadrature')
+ 6002 FORMAT('Polar quadrature')
+ 6003 FORMAT('Polar =',I5,2X,'Azim = ',I5,2X,'Group = ',I5/
+ > 'Direct angular flux per region ')
+ 6004 FORMAT('Polar =',I5,2X,'Azim = ',I5,2X,'Group = ',I5/
+ > 'Adjoint angular flux per region')
+ 6010 FORMAT(1P,5E20.10)
+ 6011 FORMAT(1P,3E20.10)
+ 6020 FORMAT('Maximum relative on flux (%) = ',F15.7)
+ 6030 FORMAT(' Flux for group ',I5)
+ 6031 FORMAT(' Region ',I5,1P,2E20.10,5X,0P,F20.10)
+ 6032 FORMAT(' Adjoint for group ',I5)
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:33:31 +0000