Skip to content

FDS Source: Fix a bug and reorganize code #16373

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cxp484 merged 1 commit into from
Jun 29, 2026
+76 −85
Merged

FDS Source: Fix a bug and reorganize code #16373

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Source/main.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -917,7 +917,7 @@ PROGRAM FDS

DO ITER=1,RADIATION_ITERATIONS
IF (RADIATION .AND. ALLOW_RANDOM_RADIATION_ROTATION .AND. ITER==1) THEN
START_NEW_ANGLE_CYCLE = MOD(MESHES(LOWER_MESH_INDEX)%RAD_CALL_COUNTER-1, & ! -1 because we already called it once duing initialization.
START_NEW_ANGLE_CYCLE = MOD(MAX(MESHES(LOWER_MESH_INDEX)%RAD_CALL_COUNTER-1,0), & ! -1 because we already called it once duing initialization.
ANGLE_INCREMENT*TIME_STEP_INCREMENT)==0
IF (START_NEW_ANGLE_CYCLE) THEN
CALL CALCULATE_DIRECTION_COEFFICIENTS()
Expand Down
159 changes: 75 additions & 84 deletions Source/radi.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2805,7 +2805,7 @@ SUBROUTINE INIT_RADIATION
USE RADCAL_CALC
USE WSGG_ARRAYS
REAL(EB) :: PLANCK_C2,KSI,LT,RCRHO,YY,BBF,AP0,AMEAN,RADIANCE,TRANSMISSIVITY,XX
INTEGER :: N,I,J,K,IPC,IZERO,NN,NI,II,JJ,IIM,JJM,IBND,NS,NRA,NSB,RADCAL_TEMP(16)=0,RCT_SKIP=-1
INTEGER :: N,I,J,K,IPC,IZERO,IBND,NS,NRA,NSB,RADCAL_TEMP(16)=0,RCT_SKIP=-1
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE (RAD_FILE_TYPE), POINTER :: RF
TYPE (SPECIES_TYPE), POINTER :: SS
Expand Down Expand Up @@ -2856,59 +2856,6 @@ SUBROUTINE INIT_RADIATION
! Set for ray rotation
ALLOW_RANDOM_RADIATION_ROTATION = RANDOMIZE_RADIATION_DIRECTIONS .AND. .NOT.CYLINDRICAL .AND. .NOT.TWO_D

! Determine mean direction normals and sweeping orders
! as described in the FDS Tech. Ref. Guide Vol. 1 Sec. 6.2.2.

! Calculate mirroring matrix
N = 0
DO I=1,NRT
DO J=1,NRP(I)
N = N + 1
DO K=1,3
IF (TWO_D .AND. .NOT.CYLINDRICAL) THEN
SELECT CASE(K)
CASE(1) ! X-surfaces
IIM = 1
JJM = NRP(I) - J + 1
CASE(2) ! Y-surfaces
IIM = 1
JJM = J
CASE(3) ! Z-surfaces
IIM = 1
JJM = NRP(I)/2 - J + 1
END SELECT
JJM = MODULO(JJM,NRP(I))
IF (JJM==0) JJM = NRP(I)
ELSE
SELECT CASE(K)
CASE(1) ! X-surfaces
IIM = I
JJM = NRP(I)/2 - J + 1
CASE(2) ! Y-surfaces
IIM = I
JJM = NRP(I) - J + 1
CASE(3) ! Z-surfaces
IIM = NRT - I + 1
JJM = J
END SELECT
IIM = MODULO(IIM,NRT)
JJM = MODULO(JJM,NRP(I))
IF (IIM==0) IIM = NRT
IF (JJM==0) JJM = NRP(I)
ENDIF

NN = 0
DO II = 1,IIM
DO JJ = 1,NRP(II)
NN = NN + 1
IF ((II==IIM).AND.(JJ==JJM)) NI = NN
ENDDO
ENDDO
DLM(N,K) = NI
ENDDO
ENDDO
ENDDO

!-----------------------------------------------------
!
! Radiative properties computation
Expand Down Expand Up @@ -3346,7 +3293,7 @@ END SUBROUTINE INIT_RADIATION
SUBROUTINE CALCULATE_FVM_ANGLES()

REAL(EB) :: THETAUP,THETALOW,PHIUP,PHILOW,F_THETA,THETA,PHI
INTEGER :: NRA,N,I,J
INTEGER :: NRA,N,I,J,K,IIM,JJM,NN,II,JJ,NI

NRA = NUMBER_RADIATION_ANGLES

Expand Down Expand Up @@ -3396,15 +3343,67 @@ SUBROUTINE CALCULATE_FVM_ANGLES()
ENDIF
ENDDO
ENDDO

! Calculate mirroring matrix
N = 0
DO I=1,NRT
DO J=1,NRP(I)
N = N + 1
DO K=1,3
IF (TWO_D .AND. .NOT.CYLINDRICAL) THEN
SELECT CASE(K)
CASE(1) ! X-surfaces
IIM = 1
JJM = NRP(I) - J + 1
CASE(2) ! Y-surfaces
IIM = 1
JJM = J
CASE(3) ! Z-surfaces
IIM = 1
JJM = NRP(I)/2 - J + 1
END SELECT
JJM = MODULO(JJM,NRP(I))
IF (JJM==0) JJM = NRP(I)
ELSE
SELECT CASE(K)
CASE(1) ! X-surfaces
IIM = I
JJM = NRP(I)/2 - J + 1
CASE(2) ! Y-surfaces
IIM = I
JJM = NRP(I) - J + 1
CASE(3) ! Z-surfaces
IIM = NRT - I + 1
JJM = J
END SELECT
IIM = MODULO(IIM,NRT)
JJM = MODULO(JJM,NRP(I))
IF (IIM==0) IIM = NRT
IF (JJM==0) JJM = NRP(I)
ENDIF

NN = 0
DO II = 1,IIM
DO JJ = 1,NRP(II)
NN = NN + 1
IF ((II==IIM).AND.(JJ==JJM)) NI = NN
ENDDO
ENDDO
DLM(N,K) = NI
ENDDO
ENDDO
ENDDO

END SUBROUTINE CALCULATE_FVM_ANGLES



!> \Calculate direction coefficients for with and without random rotations
!> as described in the FDS Tech. Ref. Guide Vol. 1 Sec. 6.2.2.
SUBROUTINE CALCULATE_DIRECTION_COEFFICIENTS()
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME

INTEGER :: NRA,N,I,J,IO,IERR
INTEGER :: NRA,N,IO,I,IERR
REAL(EB), ALLOCATABLE, DIMENSION(:) :: COSINE_ARRAY
REAL(EB) :: TNOW,AXIS(3),MERIDIAN(3), AZIMUTH(3),E1(3),E2(3),REF(3),PSI,DLO,MAGTMP,RNDNUM(3)

Expand Down Expand Up @@ -3462,35 +3461,27 @@ SUBROUTINE CALCULATE_DIRECTION_COEFFICIENTS()
AXIS = (/0._EB, 0._EB, 1._EB/)
ENDIF

N = 0
DO I=1,NRT
DO J=1,NRP(I)
N = N + 1
IF (CYLINDRICAL) THEN
DLX(N) = MERI_COMP(N)
DLY(N) = AZIM_COMP(N)
DLB(N) = DLB_COMP(N)
DLZ(N) = AXIS_COMP(N)
DLANG(1,N) = DLANG_LOCAL(1,N)
DLANG(2,N) = DLANG_LOCAL(2,N)
DLANG(3,N) = DLANG_LOCAL(3,N)
IF (N==1000000) WRITE(LU_ERR,'(A)') 'This line should never get executed. It is here only to prevent optimization.'
ELSEIF (TWO_D) THEN
DLX(N) = MERI_COMP(N)
DLY(N) = AZIM_COMP(N)
DLZ(N) = AXIS_COMP(N)
DLANG(1,N) = DLANG_LOCAL(1,N)
DLANG(2,N) = DLANG_LOCAL(2,N)
DLANG(3,N) = DLANG_LOCAL(3,N)
ELSE ! Right now random rotation is allowed only for 3D
DLX(N) = MERI_COMP(N)*MERIDIAN(1)+AZIM_COMP(N)*AZIMUTH(1)+AXIS_COMP(N)*AXIS(1)
DLY(N) = MERI_COMP(N)*MERIDIAN(2)+AZIM_COMP(N)*AZIMUTH(2)+AXIS_COMP(N)*AXIS(2)
DLZ(N) = MERI_COMP(N)*MERIDIAN(3)+AZIM_COMP(N)*AZIMUTH(3)+AXIS_COMP(N)*AXIS(3)
DLANG_OLD(:,N)=DLANG(:,N)
DLANG(:,N) = DLANG_LOCAL(1,N)*MERIDIAN(:) + DLANG_LOCAL(2,N)*AZIMUTH(:) + DLANG_LOCAL(3,N)*AXIS(:)
ENDIF

IF (CYLINDRICAL) THEN
DLX(1:NRA) = MERI_COMP(1:NRA)
DLY(1:NRA) = AZIM_COMP(1:NRA)
DLB(1:NRA) = DLB_COMP(1:NRA)
DLZ(1:NRA) = AXIS_COMP(1:NRA)
DLANG(1:3,1:NRA) = DLANG_LOCAL(1:3,1:NRA)
ELSEIF (TWO_D) THEN
DLX(1:NRA) = MERI_COMP(1:NRA)
DLY(1:NRA) = AZIM_COMP(1:NRA)
DLZ(1:NRA) = AXIS_COMP(1:NRA)
DLANG(1:3,1:NRA) = DLANG_LOCAL(1:3,1:NRA)
ELSE ! Right now random rotation is allowed only for 3D
DLX(1:NRA) = MERI_COMP(1:NRA)*MERIDIAN(1)+AZIM_COMP(1:NRA)*AZIMUTH(1)+AXIS_COMP(1:NRA)*AXIS(1)
DLY(1:NRA) = MERI_COMP(1:NRA)*MERIDIAN(2)+AZIM_COMP(1:NRA)*AZIMUTH(2)+AXIS_COMP(1:NRA)*AXIS(2)
DLZ(1:NRA) = MERI_COMP(1:NRA)*MERIDIAN(3)+AZIM_COMP(1:NRA)*AZIMUTH(3)+AXIS_COMP(1:NRA)*AXIS(3)
DLANG_OLD(:,1:NRA)=DLANG(:,1:NRA)
DO I =1,3
DLANG(I,1:NRA) = MERIDIAN(I)*DLANG_LOCAL(1,1:NRA) + AZIMUTH(I)*DLANG_LOCAL(2,1:NRA)+ AXIS(I)*DLANG_LOCAL(3,1:NRA)
ENDDO
ENDDO
ENDIF

! In axially symmetric case, each angle represents two symmetric angles. So weight the intensities by two.
WEIGH_CYL = 1._EB
Expand Down
Loading