FDS Source: Bug fix for Radiation Ray rotation when particles exist

Source/cons.f90

@@ -845,7 +845,6 @@ MODULE RADCONS
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLANG              !< Angles
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLANG_OLD          !< Angles in previous rotation
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLANG_LOCAL        !< Angles in rotating axis system
-REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: ORIENTATION_FACTOR !< Fraction of radiation angle corresponding to a particular direction
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: BBFRAC             !< Fraction of blackbody radiation
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: WL_LOW             !< Lower wavelength limit of the spectral band
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: WL_HIGH            !< Upper wavelength limit of the spectral band
@@ -858,6 +857,7 @@ MODULE RADCONS
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: DLB                !< Mean bottom component of RAYN vector (cylindrical case)
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: DLB_COMP           !< Mean bottom component of RAYN vector (cylindrical case)
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: RSA                !< Array of solid angles
+REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLO_ORIENTATION    !< COS of orientation vector w.r.t fvm angle vectors
 
 INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: DLM                !< Mirroring indices
 INTEGER, ALLOCATABLE, DIMENSION(:)    :: NRP                !< Number of radiation phi angles at each theta band

Source/radi.f90

@@ -2852,6 +2852,9 @@ SUBROUTINE INIT_RADIATION
 CALL ChkMemErr('RADI','DLANG_OLD',IZERO)
 ALLOCATE(DLANG_LOCAL(3,1:NRA),STAT=IZERO)
 CALL ChkMemErr('RADI','DLANG_LOCAL',IZERO)
+ALLOCATE(DLO_ORIENTATION(1:N_ORIENTATION_VECTOR,1:NRA),STAT=IZERO)
+CALL ChkMemErr('RADI','DLO_ORIENTATION',IZERO)
+
 
 ! Set for ray rotation
 ALLOW_RANDOM_RADIATION_ROTATION = RANDOMIZE_RADIATION_DIRECTIONS .AND. .NOT.CYLINDRICAL .AND. .NOT.TWO_D
@@ -3404,7 +3407,6 @@ SUBROUTINE CALCULATE_DIRECTION_COEFFICIENTS()
 USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
 
 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)
 
 TNOW=CURRENT_TIME()
@@ -3507,20 +3509,19 @@ SUBROUTINE CALCULATE_DIRECTION_COEFFICIENTS()
 ! DLO is the integral of the orientation vector dotted with the directional solid angle of the radiation directions.
 ! VIEW_ANGLE_FACTOR is the reduction of the radiation due to a view angle less than 180, like a narrow field of view radiometer.
 IF (SOLID_PARTICLES) THEN
-   IF (.NOT. ALLOCATED(COSINE_ARRAY)) ALLOCATE(COSINE_ARRAY(1:NRA))
    IF (.NOT. ALLOCATED(NEAREST_RADIATION_ANGLE)) ALLOCATE(NEAREST_RADIATION_ANGLE(N_ORIENTATION_VECTOR))
    IF (.NOT. ALLOCATED(VIEW_ANGLE_FACTOR)) ALLOCATE(VIEW_ANGLE_FACTOR(N_ORIENTATION_VECTOR))
    VIEW_ANGLE_FACTOR = 0._EB
    DO IO=1,N_ORIENTATION_VECTOR
       DLO = 0._EB
       DO N=1,NRA
-         COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
+         DLO_ORIENTATION(IO,N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
                            ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + &
                            ORIENTATION_VECTOR(3,IO)*DLANG(3,N)
-         IF (-COSINE_ARRAY(N) > COS_HALF_VIEW_ANGLE(IO)) &
+         IF (-DLO_ORIENTATION(IO,N) > COS_HALF_VIEW_ANGLE(IO)) &
             DLO = DLO - (ORIENTATION_VECTOR(1,IO)*DLX(N) + ORIENTATION_VECTOR(2,IO)*DLY(N) + ORIENTATION_VECTOR(3,IO)*DLZ(N))
       ENDDO
-      NEAREST_RADIATION_ANGLE(IO) = MINLOC(COSINE_ARRAY,DIM=1)
+      NEAREST_RADIATION_ANGLE(IO) = MINLOC(DLO_ORIENTATION(IO,1:NRA),DIM=1)
       VIEW_ANGLE_FACTOR(IO) = PI/DLO
    ENDDO
 ENDIF
@@ -3536,7 +3537,7 @@ SUBROUTINE INTERPOLATE_IL()
 
 INTEGER :: NRA,NM,NNN,LL,I_INTP
 REAL(EB) :: TNOW,COSANG(NUMBER_RADIATION_ANGLES), DUMMY(NUMBER_RADIATION_ANGLES)
-INTEGER :: NNEW,NOLD,IBND,IW,NOM,ICF
+INTEGER :: NNEW,NOLD,IBND,IW,NOM,ICF,IP
 INTEGER :: IDX(NUMBER_RADIATION_ANGLES,N_INTP)
 
 REAL(EB) :: W(NUMBER_RADIATION_ANGLES,N_INTP)
@@ -3546,6 +3547,8 @@ SUBROUTINE INTERPOLATE_IL()
 REAL(EB) :: ILW_OLD(NUMBER_RADIATION_ANGLES)
 TYPE (OMESH_TYPE), POINTER :: M2
 TYPE(CFACE_TYPE), POINTER :: CFA
+TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
+TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
 
 TNOW=CURRENT_TIME()
 
@@ -3592,6 +3595,7 @@ SUBROUTINE INTERPOLATE_IL()
       ENDDO
    ENDDO
 
+   ! Cut-cells
    DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
       CFA => CFACE(ICF)
       BR  => BOUNDARY_RADIA(CFA%BR_INDEX)
@@ -3605,6 +3609,25 @@ SUBROUTINE INTERPOLATE_IL()
       ENDDO
    ENDDO
 
+   ! Particles
+   DO IP=1,NLP
+      LP => LAGRANGIAN_PARTICLE(IP)
+      LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
+      IF (LPC%SOLID_PARTICLE .OR. LPC%MASSLESS_TARGET) THEN
+         IF (LP%ORIENTATION_INDEX>0) THEN
+            BR => BOUNDARY_RADIA(LP%BR_INDEX)
+            DO IBND=1,NUMBER_SPECTRAL_BANDS
+               ILW_OLD = BR%BAND(IBND)%ILW(:)
+               BR%BAND(IBND)%ILW(:) = 0._EB
+               DO I_INTP=1,N_INTP
+                  BR%BAND(IBND)%ILW(:) = BR%BAND(IBND)%ILW(:) + W(:,I_INTP)*ILW_OLD(IDX(:,I_INTP))*RSA(IDX(:,I_INTP))
+               ENDDO
+               BR%BAND(IBND)%ILW(:) = BR%BAND(IBND)%ILW(:)/RSA(:)
+            ENDDO
+         ENDIF
+      ENDIF
+   ENDDO
+
    !Interpolate neighbouring mesh cell intensities 
    DO NNN=1,N_NEIGHBORING_MESHES
       NOM = NEIGHBORING_MESH(NNN)
@@ -3691,19 +3714,20 @@ SUBROUTINE RADIATION_FVM
 REAL(EB) :: RAP, AX, AXU, AXD, AY, AYU, AYD, AZ, AZU, AZD, VC, RU, RD, RP, AFD, &
             ILXU, ILYU, ILZU, QVAL, BBF, BBFA, NCSDROP, RSA_RAT,EFLUX,SOOT_MASS_FRACTION, &
             AIU_SUM,A_SUM,VOL,VC1,AY1,AZ1,DLO,COS_DLO,AILFU, &
-            RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC
+            RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC,SUMILW
 
 INTEGER  :: N,NN,IIG,JJG,KKG,I,J,K,IW,ICF,II,JJ,KK,IOR,IC,IWUP,IWDOWN, &
             ISTART, IEND, ISTEP, JSTART, JEND, JSTEP, &
             KSTART, KEND, KSTEP, NSTART, NEND, NSTEP, &
             I_UIID, N_UPDATES, IBND, NOM, ARRAY_INDEX,NRA, &
-            IMIN, JMIN, KMIN, IMAX, JMAX, KMAX, N_SLICE, M_IJK, IJK, LL
+            IMIN, JMIN, KMIN, IMAX, JMAX, KMAX, N_SLICE, M_IJK, IJK, LL,IO
 INTEGER  :: IADD,IFACE,INDCF
 INTEGER, ALLOCATABLE :: IJK_SLICE(:,:)
-REAL(EB) :: XID,YJD,ZKD,KAPPA_PART_SINGLE,DLF,DLA(3),TSI,TMP_EXTERIOR,TEMP_ORIENTATION(3)
+REAL(EB) :: XID,YJD,ZKD,KAPPA_PART_SINGLE,DLF,DLA(3),TSI,TMP_EXTERIOR,TEMP_ORIENTATION(3),&
+            COS_DLO_ARR(NUMBER_RADIATION_ANGLES)
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET
 INTEGER :: IID,JJD,KKD,IP
-LOGICAL :: UPDATE_INTENSITY
+LOGICAL :: UPDATE_INTENSITY, IS_PARTICLE_ORIENTATION_RAMP
 REAL(EB), POINTER, DIMENSION(:,:,:) :: IL,UIIOLD,KAPPA_PART,KFST4_PART,EXTCOE,SCAEFF,SCAEFF_G,IL_UP
 REAL(EB), POINTER, DIMENSION(:)     :: OUTRAD_W,INRAD_W,OUTRAD_F,INRAD_F,IL_F
 TYPE (OMESH_TYPE), POINTER :: M2
@@ -4620,18 +4644,21 @@ SUBROUTINE RADIATION_FVM
                   IF (LPC%N_ORIENTATION==0) CYCLE PARTICLE_RADIATION_LOOP
                   BC => BOUNDARY_COORD(LP%BC_INDEX)
                   TEMP_ORIENTATION(1:3) = ORIENTATION_VECTOR(1:3,LP%ORIENTATION_INDEX)
+                  IS_PARTICLE_ORIENTATION_RAMP = .FALSE.
                   IF (LP%INIT_INDEX > 0) THEN
                      IN => INITIALIZATION(LP%INIT_INDEX)
                      IF (ANY(IN%ORIENTATION_RAMP_INDEX > 0)) THEN
-                        TEMP_ORIENTATION(1) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(1))
-                        TEMP_ORIENTATION(2) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(2))
-                        TEMP_ORIENTATION(3) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(3))
-                        TEMP_ORIENTATION = TEMP_ORIENTATION / &
-                                           (SQRT(TEMP_ORIENTATION(1)**2+TEMP_ORIENTATION(2)**2+TEMP_ORIENTATION(3)**2) &
-                                           +TWENTY_EPSILON_EB)
+                        CALL CALC_PARTICLE_TEMP_ORIENTATION(T,IN%ORIENTATION_RAMP_INDEX(1:3),TEMP_ORIENTATION)
+                        IS_PARTICLE_ORIENTATION_RAMP = .TRUE.
                      ENDIF
                   ENDIF
-                  COS_DLO = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
+
+                  IF(IS_PARTICLE_ORIENTATION_RAMP) THEN
+                     COS_DLO = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
+                  ELSE
+                     COS_DLO = -DLO_ORIENTATION(LP%ORIENTATION_INDEX,N)
+                  ENDIF
+
                   IF (COS_DLO > COS_HALF_VIEW_ANGLE(LP%ORIENTATION_INDEX)) THEN
                      DLO = -(TEMP_ORIENTATION(1)*DLX(N) + TEMP_ORIENTATION(2)*DLY(N) + TEMP_ORIENTATION(3)*DLZ(N))
                      BR => BOUNDARY_RADIA(LP%BR_INDEX)
@@ -4774,8 +4801,25 @@ SUBROUTINE RADIATION_FVM
          IF (LP%ORIENTATION_INDEX>0) THEN
             BR => BOUNDARY_RADIA(LP%BR_INDEX)
             B1%Q_RAD_IN = 0._EB
+            IO = LP%ORIENTATION_INDEX
+            IS_PARTICLE_ORIENTATION_RAMP = .FALSE.
+            IF (LP%INIT_INDEX > 0) THEN
+               IN => INITIALIZATION(LP%INIT_INDEX)
+               IF (ANY(IN%ORIENTATION_RAMP_INDEX > 0)) THEN
+                  IS_PARTICLE_ORIENTATION_RAMP = .TRUE.
+                  DO N=1,NRA
+                     CALL CALC_PARTICLE_TEMP_ORIENTATION(T,IN%ORIENTATION_RAMP_INDEX(1:3),TEMP_ORIENTATION)
+                     COS_DLO_ARR(N) = DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
+                  ENDDO
+               ENDIF
+            ENDIF
             DO IBND=1,NUMBER_SPECTRAL_BANDS
-               B1%Q_RAD_IN = B1%Q_RAD_IN + B1%EMISSIVITY * (WEIGH_CYL*SUM(BR%BAND(IBND)%ILW(1:NUMBER_RADIATION_ANGLES)) + EFLUX)
+               IF (IS_PARTICLE_ORIENTATION_RAMP) THEN
+                  SUMILW = SUM(BR%BAND(IBND)%ILW,MASK = -COS_DLO_ARR(:) > COS_HALF_VIEW_ANGLE(IO))
+               ELSE
+                  SUMILW = SUM(BR%BAND(IBND)%ILW,MASK = -DLO_ORIENTATION(IO,:) > COS_HALF_VIEW_ANGLE(IO))
+               ENDIF
+               B1%Q_RAD_IN = B1%Q_RAD_IN + B1%EMISSIVITY*(WEIGH_CYL*SUMILW + EFLUX)
             ENDDO
          ELSE
             BC => BOUNDARY_COORD(LP%BC_INDEX)
@@ -4814,6 +4858,26 @@ SUBROUTINE RADIATION_FVM
 
 END SUBROUTINE RADIATION_FVM
 
+!> \brief Provide the temporary orientation of a particle, given a ramp.
+SUBROUTINE CALC_PARTICLE_TEMP_ORIENTATION(T, ORIENTATION_RAMP_INDEX, TEMP_ORIENTATION)
+
+REAL(EB), INTENT(IN)  :: T
+INTEGER,  INTENT(IN)  :: ORIENTATION_RAMP_INDEX(3)
+REAL(EB), INTENT(OUT) :: TEMP_ORIENTATION(3)
+
+REAL(EB) :: NORM
+
+TEMP_ORIENTATION(1) = EVALUATE_RAMP(T, ORIENTATION_RAMP_INDEX(1))
+TEMP_ORIENTATION(2) = EVALUATE_RAMP(T, ORIENTATION_RAMP_INDEX(2))
+TEMP_ORIENTATION(3) = EVALUATE_RAMP(T, ORIENTATION_RAMP_INDEX(3))
+
+NORM = SQRT(SUM(TEMP_ORIENTATION**2)) + TWENTY_EPSILON_EB
+
+TEMP_ORIENTATION = TEMP_ORIENTATION / NORM
+
+END SUBROUTINE CALC_PARTICLE_TEMP_ORIENTATION
+
+
 
 !> \brief Add user-specified HRRPUV to the heat release rate term, Q
 !> \param MODE Indicator of whether volumetric heat release rate, HRRPUV, is to be added to or subtracted from Q