Skip to content

Commit 4ddbbe8

Browse files
cxp484cxp484
authored
Merge pull request #16389 from cxp484/fvdom_ray_rotation
FDS Source: Bug fix for Radiation Ray rotation when particles exist
2 parents 199afe5 + 7d3c589 commit 4ddbbe8

2 files changed

Lines changed: 83 additions & 19 deletions

File tree

‎Source/cons.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -845,7 +845,6 @@ MODULE RADCONS
845845
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLANG !< Angles
846846
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLANG_OLD !< Angles in previous rotation
847847
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLANG_LOCAL !< Angles in rotating axis system
848-
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: ORIENTATION_FACTOR !< Fraction of radiation angle corresponding to a particular direction
849848
REAL(EB), ALLOCATABLE, DIMENSION(:) :: BBFRAC !< Fraction of blackbody radiation
850849
REAL(EB), ALLOCATABLE, DIMENSION(:) :: WL_LOW !< Lower wavelength limit of the spectral band
851850
REAL(EB), ALLOCATABLE, DIMENSION(:) :: WL_HIGH !< Upper wavelength limit of the spectral band
@@ -858,6 +857,7 @@ MODULE RADCONS
858857
REAL(EB), ALLOCATABLE, DIMENSION(:) :: DLB !< Mean bottom component of RAYN vector (cylindrical case)
859858
REAL(EB), ALLOCATABLE, DIMENSION(:) :: DLB_COMP !< Mean bottom component of RAYN vector (cylindrical case)
860859
REAL(EB), ALLOCATABLE, DIMENSION(:) :: RSA !< Array of solid angles
860+
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DLO_ORIENTATION !< COS of orientation vector w.r.t fvm angle vectors
861861

862862
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: DLM !< Mirroring indices
863863
INTEGER, ALLOCATABLE, DIMENSION(:) :: NRP !< Number of radiation phi angles at each theta band

‎Source/radi.f90‎

Lines changed: 82 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -2852,6 +2852,9 @@ SUBROUTINE INIT_RADIATION
28522852
CALL ChkMemErr('RADI','DLANG_OLD',IZERO)
28532853
ALLOCATE(DLANG_LOCAL(3,1:NRA),STAT=IZERO)
28542854
CALL ChkMemErr('RADI','DLANG_LOCAL',IZERO)
2855+
ALLOCATE(DLO_ORIENTATION(1:N_ORIENTATION_VECTOR,1:NRA),STAT=IZERO)
2856+
CALL ChkMemErr('RADI','DLO_ORIENTATION',IZERO)
2857+
28552858

28562859
! Set for ray rotation
28572860
ALLOW_RANDOM_RADIATION_ROTATION = RANDOMIZE_RADIATION_DIRECTIONS .AND. .NOT.CYLINDRICAL .AND. .NOT.TWO_D
@@ -3404,7 +3407,6 @@ SUBROUTINE CALCULATE_DIRECTION_COEFFICIENTS()
34043407
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
34053408

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

34103412
TNOW=CURRENT_TIME()
@@ -3507,20 +3509,19 @@ SUBROUTINE CALCULATE_DIRECTION_COEFFICIENTS()
35073509
! DLO is the integral of the orientation vector dotted with the directional solid angle of the radiation directions.
35083510
! VIEW_ANGLE_FACTOR is the reduction of the radiation due to a view angle less than 180, like a narrow field of view radiometer.
35093511
IF (SOLID_PARTICLES) THEN
3510-
IF (.NOT. ALLOCATED(COSINE_ARRAY)) ALLOCATE(COSINE_ARRAY(1:NRA))
35113512
IF (.NOT. ALLOCATED(NEAREST_RADIATION_ANGLE)) ALLOCATE(NEAREST_RADIATION_ANGLE(N_ORIENTATION_VECTOR))
35123513
IF (.NOT. ALLOCATED(VIEW_ANGLE_FACTOR)) ALLOCATE(VIEW_ANGLE_FACTOR(N_ORIENTATION_VECTOR))
35133514
VIEW_ANGLE_FACTOR = 0._EB
35143515
DO IO=1,N_ORIENTATION_VECTOR
35153516
DLO = 0._EB
35163517
DO N=1,NRA
3517-
COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
3518+
DLO_ORIENTATION(IO,N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
35183519
ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + &
35193520
ORIENTATION_VECTOR(3,IO)*DLANG(3,N)
3520-
IF (-COSINE_ARRAY(N) > COS_HALF_VIEW_ANGLE(IO)) &
3521+
IF (-DLO_ORIENTATION(IO,N) > COS_HALF_VIEW_ANGLE(IO)) &
35213522
DLO = DLO - (ORIENTATION_VECTOR(1,IO)*DLX(N) + ORIENTATION_VECTOR(2,IO)*DLY(N) + ORIENTATION_VECTOR(3,IO)*DLZ(N))
35223523
ENDDO
3523-
NEAREST_RADIATION_ANGLE(IO) = MINLOC(COSINE_ARRAY,DIM=1)
3524+
NEAREST_RADIATION_ANGLE(IO) = MINLOC(DLO_ORIENTATION(IO,1:NRA),DIM=1)
35243525
VIEW_ANGLE_FACTOR(IO) = PI/DLO
35253526
ENDDO
35263527
ENDIF
@@ -3536,7 +3537,7 @@ SUBROUTINE INTERPOLATE_IL()
35363537

35373538
INTEGER :: NRA,NM,NNN,LL,I_INTP
35383539
REAL(EB) :: TNOW,COSANG(NUMBER_RADIATION_ANGLES), DUMMY(NUMBER_RADIATION_ANGLES)
3539-
INTEGER :: NNEW,NOLD,IBND,IW,NOM,ICF
3540+
INTEGER :: NNEW,NOLD,IBND,IW,NOM,ICF,IP
35403541
INTEGER :: IDX(NUMBER_RADIATION_ANGLES,N_INTP)
35413542

35423543
REAL(EB) :: W(NUMBER_RADIATION_ANGLES,N_INTP)
@@ -3546,6 +3547,8 @@ SUBROUTINE INTERPOLATE_IL()
35463547
REAL(EB) :: ILW_OLD(NUMBER_RADIATION_ANGLES)
35473548
TYPE (OMESH_TYPE), POINTER :: M2
35483549
TYPE(CFACE_TYPE), POINTER :: CFA
3550+
TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
3551+
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
35493552

35503553
TNOW=CURRENT_TIME()
35513554

@@ -3592,6 +3595,7 @@ SUBROUTINE INTERPOLATE_IL()
35923595
ENDDO
35933596
ENDDO
35943597

3598+
! Cut-cells
35953599
DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
35963600
CFA => CFACE(ICF)
35973601
BR => BOUNDARY_RADIA(CFA%BR_INDEX)
@@ -3605,6 +3609,25 @@ SUBROUTINE INTERPOLATE_IL()
36053609
ENDDO
36063610
ENDDO
36073611

3612+
! Particles
3613+
DO IP=1,NLP
3614+
LP => LAGRANGIAN_PARTICLE(IP)
3615+
LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
3616+
IF (LPC%SOLID_PARTICLE .OR. LPC%MASSLESS_TARGET) THEN
3617+
IF (LP%ORIENTATION_INDEX>0) THEN
3618+
BR => BOUNDARY_RADIA(LP%BR_INDEX)
3619+
DO IBND=1,NUMBER_SPECTRAL_BANDS
3620+
ILW_OLD = BR%BAND(IBND)%ILW(:)
3621+
BR%BAND(IBND)%ILW(:) = 0._EB
3622+
DO I_INTP=1,N_INTP
3623+
BR%BAND(IBND)%ILW(:) = BR%BAND(IBND)%ILW(:) + W(:,I_INTP)*ILW_OLD(IDX(:,I_INTP))*RSA(IDX(:,I_INTP))
3624+
ENDDO
3625+
BR%BAND(IBND)%ILW(:) = BR%BAND(IBND)%ILW(:)/RSA(:)
3626+
ENDDO
3627+
ENDIF
3628+
ENDIF
3629+
ENDDO
3630+
36083631
!Interpolate neighbouring mesh cell intensities
36093632
DO NNN=1,N_NEIGHBORING_MESHES
36103633
NOM = NEIGHBORING_MESH(NNN)
@@ -3691,19 +3714,20 @@ SUBROUTINE RADIATION_FVM
36913714
REAL(EB) :: RAP, AX, AXU, AXD, AY, AYU, AYD, AZ, AZU, AZD, VC, RU, RD, RP, AFD, &
36923715
ILXU, ILYU, ILZU, QVAL, BBF, BBFA, NCSDROP, RSA_RAT,EFLUX,SOOT_MASS_FRACTION, &
36933716
AIU_SUM,A_SUM,VOL,VC1,AY1,AZ1,DLO,COS_DLO,AILFU, &
3694-
RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC
3717+
RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC,SUMILW
36953718

36963719
INTEGER :: N,NN,IIG,JJG,KKG,I,J,K,IW,ICF,II,JJ,KK,IOR,IC,IWUP,IWDOWN, &
36973720
ISTART, IEND, ISTEP, JSTART, JEND, JSTEP, &
36983721
KSTART, KEND, KSTEP, NSTART, NEND, NSTEP, &
36993722
I_UIID, N_UPDATES, IBND, NOM, ARRAY_INDEX,NRA, &
3700-
IMIN, JMIN, KMIN, IMAX, JMAX, KMAX, N_SLICE, M_IJK, IJK, LL
3723+
IMIN, JMIN, KMIN, IMAX, JMAX, KMAX, N_SLICE, M_IJK, IJK, LL,IO
37013724
INTEGER :: IADD,IFACE,INDCF
37023725
INTEGER, ALLOCATABLE :: IJK_SLICE(:,:)
3703-
REAL(EB) :: XID,YJD,ZKD,KAPPA_PART_SINGLE,DLF,DLA(3),TSI,TMP_EXTERIOR,TEMP_ORIENTATION(3)
3726+
REAL(EB) :: XID,YJD,ZKD,KAPPA_PART_SINGLE,DLF,DLA(3),TSI,TMP_EXTERIOR,TEMP_ORIENTATION(3),&
3727+
COS_DLO_ARR(NUMBER_RADIATION_ANGLES)
37043728
REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET
37053729
INTEGER :: IID,JJD,KKD,IP
3706-
LOGICAL :: UPDATE_INTENSITY
3730+
LOGICAL :: UPDATE_INTENSITY, IS_PARTICLE_ORIENTATION_RAMP
37073731
REAL(EB), POINTER, DIMENSION(:,:,:) :: IL,UIIOLD,KAPPA_PART,KFST4_PART,EXTCOE,SCAEFF,SCAEFF_G,IL_UP
37083732
REAL(EB), POINTER, DIMENSION(:) :: OUTRAD_W,INRAD_W,OUTRAD_F,INRAD_F,IL_F
37093733
TYPE (OMESH_TYPE), POINTER :: M2
@@ -4620,18 +4644,21 @@ SUBROUTINE RADIATION_FVM
46204644
IF (LPC%N_ORIENTATION==0) CYCLE PARTICLE_RADIATION_LOOP
46214645
BC => BOUNDARY_COORD(LP%BC_INDEX)
46224646
TEMP_ORIENTATION(1:3) = ORIENTATION_VECTOR(1:3,LP%ORIENTATION_INDEX)
4647+
IS_PARTICLE_ORIENTATION_RAMP = .FALSE.
46234648
IF (LP%INIT_INDEX > 0) THEN
46244649
IN => INITIALIZATION(LP%INIT_INDEX)
46254650
IF (ANY(IN%ORIENTATION_RAMP_INDEX > 0)) THEN
4626-
TEMP_ORIENTATION(1) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(1))
4627-
TEMP_ORIENTATION(2) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(2))
4628-
TEMP_ORIENTATION(3) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(3))
4629-
TEMP_ORIENTATION = TEMP_ORIENTATION / &
4630-
(SQRT(TEMP_ORIENTATION(1)**2+TEMP_ORIENTATION(2)**2+TEMP_ORIENTATION(3)**2) &
4631-
+TWENTY_EPSILON_EB)
4651+
CALL CALC_PARTICLE_TEMP_ORIENTATION(T,IN%ORIENTATION_RAMP_INDEX(1:3),TEMP_ORIENTATION)
4652+
IS_PARTICLE_ORIENTATION_RAMP = .TRUE.
46324653
ENDIF
46334654
ENDIF
4634-
COS_DLO = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
4655+
4656+
IF(IS_PARTICLE_ORIENTATION_RAMP) THEN
4657+
COS_DLO = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
4658+
ELSE
4659+
COS_DLO = -DLO_ORIENTATION(LP%ORIENTATION_INDEX,N)
4660+
ENDIF
4661+
46354662
IF (COS_DLO > COS_HALF_VIEW_ANGLE(LP%ORIENTATION_INDEX)) THEN
46364663
DLO = -(TEMP_ORIENTATION(1)*DLX(N) + TEMP_ORIENTATION(2)*DLY(N) + TEMP_ORIENTATION(3)*DLZ(N))
46374664
BR => BOUNDARY_RADIA(LP%BR_INDEX)
@@ -4774,8 +4801,25 @@ SUBROUTINE RADIATION_FVM
47744801
IF (LP%ORIENTATION_INDEX>0) THEN
47754802
BR => BOUNDARY_RADIA(LP%BR_INDEX)
47764803
B1%Q_RAD_IN = 0._EB
4804+
IO = LP%ORIENTATION_INDEX
4805+
IS_PARTICLE_ORIENTATION_RAMP = .FALSE.
4806+
IF (LP%INIT_INDEX > 0) THEN
4807+
IN => INITIALIZATION(LP%INIT_INDEX)
4808+
IF (ANY(IN%ORIENTATION_RAMP_INDEX > 0)) THEN
4809+
IS_PARTICLE_ORIENTATION_RAMP = .TRUE.
4810+
DO N=1,NRA
4811+
CALL CALC_PARTICLE_TEMP_ORIENTATION(T,IN%ORIENTATION_RAMP_INDEX(1:3),TEMP_ORIENTATION)
4812+
COS_DLO_ARR(N) = DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
4813+
ENDDO
4814+
ENDIF
4815+
ENDIF
47774816
DO IBND=1,NUMBER_SPECTRAL_BANDS
4778-
B1%Q_RAD_IN = B1%Q_RAD_IN + B1%EMISSIVITY * (WEIGH_CYL*SUM(BR%BAND(IBND)%ILW(1:NUMBER_RADIATION_ANGLES)) + EFLUX)
4817+
IF (IS_PARTICLE_ORIENTATION_RAMP) THEN
4818+
SUMILW = SUM(BR%BAND(IBND)%ILW,MASK = -COS_DLO_ARR(:) > COS_HALF_VIEW_ANGLE(IO))
4819+
ELSE
4820+
SUMILW = SUM(BR%BAND(IBND)%ILW,MASK = -DLO_ORIENTATION(IO,:) > COS_HALF_VIEW_ANGLE(IO))
4821+
ENDIF
4822+
B1%Q_RAD_IN = B1%Q_RAD_IN + B1%EMISSIVITY*(WEIGH_CYL*SUMILW + EFLUX)
47794823
ENDDO
47804824
ELSE
47814825
BC => BOUNDARY_COORD(LP%BC_INDEX)
@@ -4814,6 +4858,26 @@ SUBROUTINE RADIATION_FVM
48144858

48154859
END SUBROUTINE RADIATION_FVM
48164860

4861+
!> \brief Provide the temporary orientation of a particle, given a ramp.
4862+
SUBROUTINE CALC_PARTICLE_TEMP_ORIENTATION(T, ORIENTATION_RAMP_INDEX, TEMP_ORIENTATION)
4863+
4864+
REAL(EB), INTENT(IN) :: T
4865+
INTEGER, INTENT(IN) :: ORIENTATION_RAMP_INDEX(3)
4866+
REAL(EB), INTENT(OUT) :: TEMP_ORIENTATION(3)
4867+
4868+
REAL(EB) :: NORM
4869+
4870+
TEMP_ORIENTATION(1) = EVALUATE_RAMP(T, ORIENTATION_RAMP_INDEX(1))
4871+
TEMP_ORIENTATION(2) = EVALUATE_RAMP(T, ORIENTATION_RAMP_INDEX(2))
4872+
TEMP_ORIENTATION(3) = EVALUATE_RAMP(T, ORIENTATION_RAMP_INDEX(3))
4873+
4874+
NORM = SQRT(SUM(TEMP_ORIENTATION**2)) + TWENTY_EPSILON_EB
4875+
4876+
TEMP_ORIENTATION = TEMP_ORIENTATION / NORM
4877+
4878+
END SUBROUTINE CALC_PARTICLE_TEMP_ORIENTATION
4879+
4880+
48174881

48184882
!> \brief Add user-specified HRRPUV to the heat release rate term, Q
48194883
!> \param MODE Indicator of whether volumetric heat release rate, HRRPUV, is to be added to or subtracted from Q

0 commit comments

Comments
 (0)