diff --git a/Build/makefile b/Build/makefile index 30c237b8bc..c82e76d026 100644 --- a/Build/makefile +++ b/Build/makefile @@ -469,7 +469,7 @@ soot.o : prec.o mesh.o cons.o func.o turb.o pois.o : prec.o cons.o geom.o : gsmv.o func.o prec.o cons.o type.o mesh.o turb.o imkl.o data.o ccib.o : gsmv.o func.o prec.o cons.o type.o mesh.o turb.o geom.o -radi.o : func.o prec.o cons.o type.o mesh.o rcal.o geom.o +radi.o : func.o prec.o cons.o type.o mesh.o rcal.o geom.o ccib.o part.o : func.o prec.o cons.o type.o mesh.o devc.o data.o geom.o soot.o ccib.o vege.o : func.o prec.o cons.o type.o mesh.o geom.o ctrl.o : prec.o cons.o type.o mesh.o func.o @@ -504,7 +504,7 @@ soot.obj : prec.obj mesh.obj cons.obj func.obj turb.obj pois.obj : prec.obj cons.obj geom.obj : gsmv.obj func.obj prec.obj cons.obj type.obj mesh.obj turb.obj imkl.obj data.obj ccib.obj : gsmv.obj func.obj prec.obj cons.obj type.obj mesh.obj turb.obj geom.obj -radi.obj : func.obj prec.obj cons.obj type.obj mesh.obj rcal.obj geom.obj +radi.obj : func.obj prec.obj cons.obj type.obj mesh.obj rcal.obj geom.obj ccib.obj part.obj : func.obj prec.obj cons.obj type.obj mesh.obj devc.obj data.obj geom.obj soot.obj ccib.obj vege.obj : func.obj prec.obj cons.obj type.obj mesh.obj geom.obj ctrl.obj : prec.obj cons.obj type.obj mesh.obj func.obj diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 794b5ad43a..a39abc4911 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -132,7 +132,8 @@ MODULE CC_SCALARS GET_CFACE_OPEN_BC_COEF,GET_FN_DIVERGENCE_CUTCELL,GET_OPENBC_TANGENTIAL_CUTFACE_VEL,& GET_CUTCELL_DDDT,GET_H_CUTFACES,GET_H_MATRIX_CC,GET_H_GUARD_CUTCELL,GET_CRTCFCC_INT_STENCILS,GET_RCFACES_H, & GET_CC_MATRIXGRAPH_H,GET_CC_IROW,GET_CC_UNKH,GET_CUTCELL_HP, GET_LINKED_FV, GET_PRES_CFACE_BCS, & - GET_PRES_CFACE, GET_PRES_CFACE_TEST, GET_UVWGAS_CFACE, GET_MUDNS_CFACE, GET_BOUNDFACE_GEOM_INFO_H, & + GET_PRES_CFACE, GET_PRES_CFACE_TEST, GET_UVWGAS_CFACE, GET_CFACE_RAD_NVEC, GET_MUDNS_CFACE, & + GET_BOUNDFACE_GEOM_INFO_H, & GET_FH_FROM_PRHS_AND_BCS,GET_LINKED_VELOCITIES, & FINISH_CC, INIT_CUTCELL_DATA,MESH_CC_EXCHANGE,NUMBER_UNKH_CUTCELLS,& ROTATED_CUBE_ANN_SOLN,ROTATED_CUBE_VELOCITY_FLUX,ROTATED_CUBE_RHS_ZZ,& @@ -154,11 +155,11 @@ MODULE CC_SCALARS SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE INTEGER :: NM,NOM,IQ,NQ,IIO,JJO,KKO,NOM2,ICFACE,SURF_B,FRONT,OD_INDEX,N_RESOLVED,N_GLOBAL,IP,K,P,TOTAL_SEND,TOTAL_RECV,IERR -REAL(EB) :: POS(3) +REAL(EB) :: POS(3),NBACK(3) INTEGER, ALLOCATABLE :: SEND_CNT(:),RECV_CNT(:),SDISPL(:),RDISPL(:),CNT3_S(:),CNT3_R(:),DSP3_S(:),DSP3_R(:),COUNTER(:) INTEGER, ALLOCATABLE :: SEND_NOM(:),RECV_NOM(:),REPLY_IDX(:),REPLY_SURF(:),REPLY_DIM(:), & REPLY_BACK(:),REPLY_BSURF(:),REPLY_BDIM(:),MAP_NM(:),MAP_FRONT(:) -REAL(EB), ALLOCATABLE :: SEND_POS(:),RECV_POS(:) +REAL(EB), ALLOCATABLE :: SEND_POS(:),RECV_POS(:),SEND_NBACK(:),RECV_NBACK(:) N_RESOLVED = 0 ALLOCATE(SEND_CNT(0:N_MPI_PROCESSES-1)) ; SEND_CNT = 0 @@ -169,10 +170,11 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE IF (NQ==0) CYCLE IF (PROCESS(NOM)==MY_RANK) THEN ! Owner is on this process: resolve directly. DO IQ=1,NQ - POS(1:3) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ) + POS(1:3) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ) + NBACK(1:3) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_NBACK(IAXIS:KAXIS,IQ) CALL SEARCH_OTHER_MESHES(POS(1),POS(2),POS(3),NOM2,IIO,JJO,KKO) IF (NOM2/=NOM) CYCLE - ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1) + ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1,NBACK) IF (ICFACE>0) THEN FRONT = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_FRONT(IQ) OD_INDEX = MESHES(NM)%CFACE(FRONT)%OD_INDEX @@ -206,10 +208,12 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE TOTAL_SEND = SDISPL(N_MPI_PROCESSES-1) + SEND_CNT(N_MPI_PROCESSES-1) TOTAL_RECV = RDISPL(N_MPI_PROCESSES-1) + RECV_CNT(N_MPI_PROCESSES-1) CNT3_S = 3*SEND_CNT ; CNT3_R = 3*RECV_CNT ; DSP3_S = 3*SDISPL ; DSP3_R = 3*RDISPL - ALLOCATE(SEND_NOM(MAX(TOTAL_SEND,1)),SEND_POS(MAX(3*TOTAL_SEND,1)),MAP_NM(MAX(TOTAL_SEND,1)), & + ALLOCATE(SEND_NOM(MAX(TOTAL_SEND,1)),SEND_POS(MAX(3*TOTAL_SEND,1)),SEND_NBACK(MAX(3*TOTAL_SEND,1)), & + MAP_NM(MAX(TOTAL_SEND,1)), & MAP_FRONT(MAX(TOTAL_SEND,1)),REPLY_BACK(MAX(TOTAL_SEND,1)),REPLY_BSURF(MAX(TOTAL_SEND,1)), & REPLY_BDIM(MAX(TOTAL_SEND,1))) - ALLOCATE(RECV_NOM(MAX(TOTAL_RECV,1)),RECV_POS(MAX(3*TOTAL_RECV,1)),REPLY_IDX(MAX(TOTAL_RECV,1)), & + ALLOCATE(RECV_NOM(MAX(TOTAL_RECV,1)),RECV_POS(MAX(3*TOTAL_RECV,1)),RECV_NBACK(MAX(3*TOTAL_RECV,1)), & + REPLY_IDX(MAX(TOTAL_RECV,1)), & REPLY_SURF(MAX(TOTAL_RECV,1)),REPLY_DIM(MAX(TOTAL_RECV,1))) ALLOCATE(COUNTER(0:N_MPI_PROCESSES-1)) ; COUNTER = SDISPL @@ -224,24 +228,28 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE IP = PROCESS(NOM) DO IQ=1,NQ P = COUNTER(IP) + 1 ; COUNTER(IP) = P - SEND_NOM(P) = NOM - SEND_POS(3*P-2:3*P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ) - MAP_NM(P) = NM - MAP_FRONT(P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_FRONT(IQ) + SEND_NOM(P) = NOM + SEND_POS(3*P-2:3*P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ) + SEND_NBACK(3*P-2:3*P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_NBACK(IAXIS:KAXIS,IQ) + MAP_NM(P) = NM + MAP_FRONT(P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_FRONT(IQ) ENDDO ENDDO ENDDO CALL MPI_ALLTOALLV(SEND_NOM,SEND_CNT,SDISPL,MPI_INTEGER,RECV_NOM,RECV_CNT,RDISPL,MPI_INTEGER,MPI_COMM_WORLD,IERR) CALL MPI_ALLTOALLV(SEND_POS,CNT3_S,DSP3_S,MPI_DOUBLE_PRECISION,RECV_POS,CNT3_R,DSP3_R,MPI_DOUBLE_PRECISION,& MPI_COMM_WORLD,IERR) + CALL MPI_ALLTOALLV(SEND_NBACK,CNT3_S,DSP3_S,MPI_DOUBLE_PRECISION,RECV_NBACK,CNT3_R,DSP3_R,MPI_DOUBLE_PRECISION,& + MPI_COMM_WORLD,IERR) ! Owner side: resolve each received query against its local cut-face data. DO K=1,TOTAL_RECV NOM = RECV_NOM(K) - POS(1:3) = RECV_POS(3*K-2:3*K) + POS(1:3) = RECV_POS(3*K-2:3*K) + NBACK(1:3) = RECV_NBACK(3*K-2:3*K) CALL SEARCH_OTHER_MESHES(POS(1),POS(2),POS(3),NOM2,IIO,JJO,KKO) ICFACE = 0 ; SURF_B = 0 - IF (NOM2==NOM) ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1) + IF (NOM2==NOM) ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1,NBACK) IF (ICFACE>0) SURF_B = MESHES(NOM)%CFACE(ICFACE)%SURF_INDEX REPLY_IDX(K) = ICFACE REPLY_SURF(K) = SURF_B @@ -304,8 +312,8 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE CALL MPI_ALLREDUCE(MPI_IN_PLACE,N_GLOBAL,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR) IF (N_GLOBAL>0) CF_BACK_EXCHANGE_ACTIVE = .TRUE. - DEALLOCATE(RECV_CNT,SDISPL,RDISPL,CNT3_S,CNT3_R,DSP3_S,DSP3_R,COUNTER,SEND_NOM,SEND_POS,MAP_NM,MAP_FRONT, & - REPLY_BACK,REPLY_BSURF,REPLY_BDIM,RECV_NOM,RECV_POS,REPLY_IDX,REPLY_SURF,REPLY_DIM) + DEALLOCATE(RECV_CNT,SDISPL,RDISPL,CNT3_S,CNT3_R,DSP3_S,DSP3_R,COUNTER,SEND_NOM,SEND_POS,SEND_NBACK,MAP_NM,MAP_FRONT, & + REPLY_BACK,REPLY_BSURF,REPLY_BDIM,RECV_NOM,RECV_POS,RECV_NBACK,REPLY_IDX,REPLY_SURF,REPLY_DIM) ENDIF CROSS_PROCESS_IF DEALLOCATE(SEND_CNT) @@ -622,7 +630,7 @@ SUBROUTINE GET_H_GUARD_CUTCELL(IPZ,HP) ELSE; D_EXT = ABS(MESHES(NOM)%ZC(KKO) - XFN) ENDIF END SELECT - ! Refinement seam only: match GET_H_MATRIX_CC sub-face selection (skip solid/no-unknown partners). + ! Refinement seam only: match GET_H_MATRIX_CC sub-face selection (skip solid/no-unknown partners). IF (WC%BOUNDARY_TYPE == INTERPOLATED_BOUNDARY) THEN IF (ALLOCATED(OM%MUNKH)) THEN; IF (OM%MUNKH(IIO,JJO,KKO) < 1) CYCLE; ENDIF IF (MESHES(NOM)%FCVAR(IIO+ISHF(IAXIS),JJO+ISHF(JAXIS),KKO+ISHF(KAXIS),CC_FGSC,X1AXIS) == CC_SOLID) CYCLE @@ -6756,6 +6764,20 @@ SUBROUTINE GET_UVWGAS_CFACE(U_CELL,V_CELL,W_CELL,IND1,IND2,UP,VP,WP,PREDFCT) RETURN END SUBROUTINE GET_UVWGAS_CFACE +! ----------------------------------- GET_CFACE_RAD_NVEC --------------------------------- +FUNCTION GET_CFACE_RAD_NVEC(CFA) RESULT(NVEC) + +TYPE(CFACE_TYPE), INTENT(IN) :: CFA +REAL(EB) :: NVEC(IAXIS:KAXIS) +INTEGER :: IBOD,IWSEL + +NVEC(IAXIS:KAXIS) = BOUNDARY_COORD(CFA%BC_INDEX)%NVEC(IAXIS:KAXIS) +IBOD = CUT_FACE(CFA%CUT_FACE_IND1)%BODTRI(1,CFA%CUT_FACE_IND2) +IWSEL= CUT_FACE(CFA%CUT_FACE_IND1)%BODTRI(2,CFA%CUT_FACE_IND2) +IF (IBOD>=1 .AND. IWSEL>=1) NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) + +END FUNCTION GET_CFACE_RAD_NVEC + ! ----------------------------------- GET_PRES_CFACE --------------------------------- SUBROUTINE GET_PRES_CFACE(PRESS,IND1,IND2,CFA) @@ -8936,7 +8958,7 @@ SUBROUTINE CCSPECIES_ADVECTION ! region and adds components to thermodynamic divergence. ! Local Variables: -REAL(EB) :: RHO_Z_PV(-2:1), VELC, FN_ZZ, ZZ_GET_N +REAL(EB) :: RHO_Z_PV(-2:1), VELC, FN_ZZ, FN_ZZ_USE, ZZ_GET_N REAL(EB), PARAMETER :: SGNFCT=1._EB INTEGER :: IOR, ICFA TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 @@ -9127,6 +9149,12 @@ SUBROUTINE CCSPECIES_ADVECTION CASE(INTERPOLATED_BOUNDARY) VELC = UVW_SAVE(IW) END SELECT + ! For a cut-cell against a mesh interface,flux-limited value built across the interface. + IF ((WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) .AND. & + RCF%CELL_LIST(1,ISIDE+2)==CC_FTYPE_CFGAS) THEN + ICC = RCF%CELL_LIST(2,ISIDE+2); JCC = RCF%CELL_LIST(3,ISIDE+2) + CALL GET_CCIFACE_RHOZZ_LIM(X1AXIS,I,J,K,ISIDE+2,ICC,JCC,N,VELC,PRFCT,FN_ZZ) + ENDIF SELECT CASE(RCF%CELL_LIST(1,ISIDE+2)) CASE(CC_FTYPE_RGGAS) SELECT CASE(X1AXIS) @@ -9342,8 +9370,12 @@ SUBROUTINE CCSPECIES_ADVECTION RHOPV(ISIDE) = PRFCT *CUT_CELL(ICC)%RHO(JCC) + (1._EB-PRFCT) *CUT_CELL(ICC)%RHOS(JCC) ZZ_GET_N = PRFCT*CUT_CELL(ICC)%ZZ(N,JCC) + (1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(N,JCC) RHO_Z_PV(ISIDE) = RHOPV(ISIDE)*ZZ_GET_N + ! At a mesh interface use the flux-limited bar{rho*Z} across the interface. + FN_ZZ_USE = FN_ZZ + IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) & + CALL GET_CCIFACE_RHOZZ_LIM(X1AXIS,I,J,K,ISIDE+2,ICC,JCC,N,VELC,PRFCT,FN_ZZ_USE) CUT_CELL(ICC)%U_DOT_DEL_RHO_Z_VOL(N,JCC) = CUT_CELL(ICC)%U_DOT_DEL_RHO_Z_VOL(N,JCC) + & - FCT*SGNFCT*(FN_ZZ-RHO_Z_PV(ISIDE))*VELC*AF + FCT*SGNFCT*(FN_ZZ_USE-RHO_Z_PV(ISIDE))*VELC*AF ENDIF ENDDO ! IFACE @@ -9473,7 +9505,7 @@ SUBROUTINE CCENTHALPY_ADVECTION ! Local Variables: -REAL(EB) :: RHO_H_S_PV(-2:1), VELC, VELC2, FN_H_S, TMP_F_GAS +REAL(EB) :: RHO_H_S_PV(-2:1), VELC, VELC2, FN_H_S, FN_H_S_USE, TMP_F_GAS INTEGER :: IOR, ICFA TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 TYPE(CC_REGFACEZ_TYPE), POINTER, DIMENSION(:) :: REGFACE_Z @@ -10070,6 +10102,12 @@ SUBROUTINE CCENTHALPY_ADVECTION ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMP_F_GAS) FN_H_S = B1%RHO_F*H_S ! bar{rho*hs} + ! At a mesh interface with a cut-cell, use the flux-limited bar{rho*h_s} built across the interface. + IF ((WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) .AND. & + RC_FACE(IFACE)%CELL_LIST(1,ISIDE+2)==CC_FTYPE_CFGAS) THEN + ICC = RC_FACE(IFACE)%CELL_LIST(2,ISIDE+2); JCC = RC_FACE(IFACE)%CELL_LIST(3,ISIDE+2) + CALL GET_CCIFACE_RHOHS_LIM(X1AXIS,I,J,K,ISIDE+2,ICC,JCC,VELC,PRFCT,FN_H_S) + ENDIF ! Finally add to Div: SELECT CASE(RC_FACE(IFACE)%CELL_LIST(1,ISIDE+2)) CASE(CC_FTYPE_RGGAS) ! Regular cell @@ -10133,7 +10171,13 @@ SUBROUTINE CCENTHALPY_ADVECTION (1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(1:N_TRACKED_SPECIES,JCC) CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMPV(ISIDE)) RHO_H_S_PV(ISIDE) = RHOPV(ISIDE)*H_S - CUT_CELL(ICC)%DVOL(JCC)=CUT_CELL(ICC)%DVOL(JCC)+FCT*(-1._EB)*(FN_H_S-RHO_H_S_PV(ISIDE))*VELC*AF ! +ve or -ve dot + ! At a mesh interface use the flux-limited bar{rho*h_s} built across the interface. + FN_H_S_USE = FN_H_S + IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) & + CALL GET_CCIFACE_RHOHS_LIM(CUT_FACE(ICF)%IJK(KAXIS+1), & + CUT_FACE(ICF)%IJK(IAXIS),CUT_FACE(ICF)%IJK(JAXIS),CUT_FACE(ICF)%IJK(KAXIS), & + ISIDE+2,ICC,JCC,VELC,PRFCT,FN_H_S_USE) + CUT_CELL(ICC)%DVOL(JCC)=CUT_CELL(ICC)%DVOL(JCC)+FCT*(-1._EB)*(FN_H_S_USE-RHO_H_S_PV(ISIDE))*VELC*AF ! +ve or -ve dot END SELECT ENDDO ! IFACE ENDDO ! ICF @@ -13449,6 +13493,112 @@ SUBROUTINE GET_RHOZZVECTOR_SCALAR_3D(N) END SUBROUTINE GET_RHOZZVECTOR_SCALAR_3D +! ----------------------------- GET_CCIFACE_RHOZZ_LIM -------------------------- +! Flux-limited bar{rho*Z_N} face value at an INTERPOLATED / PERIODIC-boundary gasphase cut-face. +SUBROUTINE GET_CCIFACE_RHOZZ_LIM(X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC,N,VELC,PRFCT,FN_ZZ) + +INTEGER, INTENT(IN) :: X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC,N +REAL(EB), INTENT(IN) :: VELC,PRFCT +REAL(EB), INTENT(OUT):: FN_ZZ +INTEGER :: DD,ILSD,ILOUT,II,JJ,KK,SHFT(3) +REAL(EB) :: RZPV(-2:1) +REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP +REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP +REAL(EB), TARGET :: ZT(0:3,1,1),UT(1,1,1),FT(1,1,1) +REAL(EB), POINTER, DIMENSION(:,:,:) :: ZTP,UTP,FTP + +IF (PRFCT > 0.5_EB) THEN; RHOP=>RHO ; ZZP=>ZZ ! PRFCT=1 -> n-level fields +ELSE ; RHOP=>RHOS; ZZP=>ZZS ! PRFCT=0 -> star fields +ENDIF + +SELECT CASE(X1AXIS) +CASE(IAXIS); SHFT=(/1,0,0/) +CASE(JAXIS); SHFT=(/0,1,0/) +CASE(KAXIS); SHFT=(/0,0,1/) +END SELECT + +! Cartesian rho*Z on the four stencil cells: cell(DD) = base_face_cell + (1+DD) along the face axis. +DO DD=-2,1 + II=IFC+(1+DD)*SHFT(1); JJ=JFC+(1+DD)*SHFT(2); KK=KFC+(1+DD)*SHFT(3) + RZPV(DD)=RHOP(II,JJ,KK)*ZZP(II,JJ,KK,N) +ENDDO + +! Override the local immediate cut-cell (side ILOC) with the detailed cut-cell rho*Z: +ILSD = ILOC-2 ! LOW_IND(1) -> -1 ; HIGH_IND(2) -> 0 +RZPV(ILSD) = (PRFCT*CUT_CELL(ICC)%RHO(JCC) +(1._EB-PRFCT)*CUT_CELL(ICC)%RHOS(JCC)) * & + (PRFCT*CUT_CELL(ICC)%ZZ(N,JCC)+(1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(N,JCC)) + +! Solid fall-back only for the local-outer stencil cell (the ghost-outer cell is remote gas): +IF (ILOC==LOW_IND) THEN; ILOUT=-2; ELSE; ILOUT=1; ENDIF +II=IFC+(1+ILOUT)*SHFT(1); JJ=JFC+(1+ILOUT)*SHFT(2); KK=KFC+(1+ILOUT)*SHFT(3) +IF (II>=1.AND.II<=IBAR .AND. JJ>=1.AND.JJ<=JBAR .AND. KK>=1.AND.KK<=KBAR) THEN + IF (CELL(CELL_INDEX(II,JJ,KK))%SOLID .OR. CCVAR(II,JJ,KK,CC_CGSC)==CC_SOLID) RZPV(ILOUT)=RZPV(ILSD) +ENDIF +ZT(0:3,1,1)=RZPV(-2:1); UT(1,1,1)=VELC +ZTP=>ZT; UTP=>UT; FTP=>FT +CALL GET_SCALAR_FACE_VALUE(UTP,ZTP,FTP,1,1,1,1,1,1,1,I_FLUX_LIMITER) +FN_ZZ=FT(1,1,1) + +RETURN +END SUBROUTINE GET_CCIFACE_RHOZZ_LIM + + +! ----------------------------- GET_CCIFACE_RHOHS_LIM -------------------------- +! Enthalpy counterpart of GET_CCIFACE_RHOZZ_LIM: flux-limited bar{rho*h_s} face value at an +! INTERPOLATED / PERIODIC-boundary gasphase cut-face. +SUBROUTINE GET_CCIFACE_RHOHS_LIM(X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC,VELC,PRFCT,FN_H_S) +USE PHYSICAL_FUNCTIONS, ONLY: GET_SENSIBLE_ENTHALPY + +INTEGER, INTENT(IN) :: X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC +REAL(EB), INTENT(IN) :: VELC,PRFCT +REAL(EB), INTENT(OUT):: FN_H_S +INTEGER :: DD,ILSD,ILOUT,II,JJ,KK,SHFT(3) +REAL(EB) :: RHS_PV(-2:1),HS1,ZZ_GET(1:N_TRACKED_SPECIES) +REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP +REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP +REAL(EB), TARGET :: ZT(0:3,1,1),UT(1,1,1),FT(1,1,1) +REAL(EB), POINTER, DIMENSION(:,:,:) :: ZTP,UTP,FTP + +IF (PRFCT > 0.5_EB) THEN; RHOP=>RHO ; ZZP=>ZZ +ELSE ; RHOP=>RHOS; ZZP=>ZZS +ENDIF + +SELECT CASE(X1AXIS) +CASE(IAXIS); SHFT=(/1,0,0/) +CASE(JAXIS); SHFT=(/0,1,0/) +CASE(KAXIS); SHFT=(/0,0,1/) +END SELECT + +! Cartesian rho*h_s on the four stencil cells: cell(DD) = base_face_cell + (1+DD) along the axis. +DO DD=-2,1 + II=IFC+(1+DD)*SHFT(1); JJ=JFC+(1+DD)*SHFT(2); KK=KFC+(1+DD)*SHFT(3) + ZZ_GET(1:N_TRACKED_SPECIES)=ZZP(II,JJ,KK,1:N_TRACKED_SPECIES) + CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,HS1,TMP(II,JJ,KK)) + RHS_PV(DD)=RHOP(II,JJ,KK)*HS1 +ENDDO + +! Override the local immediate cut-cell (side ILOC) with the detailed cut-cell rho*h_s: +ILSD = ILOC-2 +ZZ_GET(1:N_TRACKED_SPECIES)=PRFCT*CUT_CELL(ICC)%ZZ(1:N_TRACKED_SPECIES,JCC)+ & + (1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(1:N_TRACKED_SPECIES,JCC) +CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,HS1,CUT_CELL(ICC)%TMP(JCC)) +RHS_PV(ILSD)=(PRFCT*CUT_CELL(ICC)%RHO(JCC)+(1._EB-PRFCT)*CUT_CELL(ICC)%RHOS(JCC))*HS1 + +! Solid fall-back for the local-outer stencil cell only: +IF (ILOC==LOW_IND) THEN; ILOUT=-2; ELSE; ILOUT=1; ENDIF +II=IFC+(1+ILOUT)*SHFT(1); JJ=JFC+(1+ILOUT)*SHFT(2); KK=KFC+(1+ILOUT)*SHFT(3) +IF (II>=1.AND.II<=IBAR .AND. JJ>=1.AND.JJ<=JBAR .AND. KK>=1.AND.KK<=KBAR) THEN + IF (CELL(CELL_INDEX(II,JJ,KK))%SOLID .OR. CCVAR(II,JJ,KK,CC_CGSC)==CC_SOLID) RHS_PV(ILOUT)=RHS_PV(ILSD) +ENDIF +ZT(0:3,1,1)=RHS_PV(-2:1); UT(1,1,1)=VELC +ZTP=>ZT; UTP=>UT; FTP=>FT +CALL GET_SCALAR_FACE_VALUE(UTP,ZTP,FTP,1,1,1,1,1,1,1,I_FLUX_LIMITER) +FN_H_S=FT(1,1,1) + +RETURN +END SUBROUTINE GET_CCIFACE_RHOHS_LIM + + ! ------------------------- GET_ADVDIFFVECTOR_SCALAR_3D ------------------------- SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N) @@ -13460,7 +13610,7 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N) REAL(EB):: PRFCT INTEGER :: X1AXIS,IFACE,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND),ICF,IND1,IND2,IOR INTEGER :: LOCROW_1,LOCROW_2,ILOC,IROW,ICC,JCC,IW -REAL(EB):: AF,VELC,RHO_Z,FN_ZZ,FCT +REAL(EB):: AF,VELC,RHO_Z,FN_ZZ,FCT,FN_ZZ_USE REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP,UP,VP,WP REAL(EB), POINTER, DIMENSION(:,:,:,:):: ZZP TYPE(CC_REGFACEZ_TYPE), POINTER, DIMENSION(:) :: REGFACE_Z @@ -13620,6 +13770,12 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N) CASE(INTERPOLATED_BOUNDARY) VELC = UVW_SAVE(IW) END SELECT + ! For a cut-cell against a mesh interface , use the flux-limited bar{rho*Z} built across the interface. + IF ((WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) .AND. & + MESHES(NM)%RC_FACE(IFACE)%CELL_LIST(1,ILOC)==CC_FTYPE_CFGAS) THEN + ICC = MESHES(NM)%RC_FACE(IFACE)%CELL_LIST(2,ILOC); JCC = MESHES(NM)%RC_FACE(IFACE)%CELL_LIST(3,ILOC) + CALL GET_CCIFACE_RHOZZ_LIM(X1AXIS,I,J,K,ILOC,ICC,JCC,N,VELC,PRFCT,FN_ZZ) + ENDIF IROW=IND_LOC(ILOC) ! Process Local Unknown number. FCT = REAL(3-2*ILOC,EB) @@ -13652,7 +13808,15 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N) ! First (rho hs)_i,j,k: IF (CUT_FACE(ICF)%CELL_LIST(1,ILOC,IFACE) == CC_FTYPE_CFGAS) THEN IROW=IND_LOC(ILOC) ! Process Local Unknown number. - F_Z(IROW) = F_Z(IROW) + FCT*FN_ZZ*VELC*AF + FN_ZZ_USE = FN_ZZ + ! Build the flux-limited bar{rho*Z} across the interface. + IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) THEN + ICC = CUT_FACE(ICF)%CELL_LIST(2,ILOC,IFACE); JCC = CUT_FACE(ICF)%CELL_LIST(3,ILOC,IFACE) + CALL GET_CCIFACE_RHOZZ_LIM(CUT_FACE(ICF)%IJK(KAXIS+1), & + CUT_FACE(ICF)%IJK(IAXIS),CUT_FACE(ICF)%IJK(JAXIS),CUT_FACE(ICF)%IJK(KAXIS), & + ILOC,ICC,JCC,N,VELC,PRFCT,FN_ZZ_USE) + ENDIF + F_Z(IROW) = F_Z(IROW) + FCT*FN_ZZ_USE*VELC*AF ENDIF ENDDO ! IFACE diff --git a/Source/geom.f90 b/Source/geom.f90 index ba248ace61..d7697d1d19 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -5452,7 +5452,7 @@ MODULE COMPLEX_GEOMETRY_GRID INTEGER, PARAMETER :: GLOBAL_DELTA_EDGE = 3*GLOBAL_DELTA_CELL INTEGER, PARAMETER :: GLOBAL_DELTA_FACE = 3*GLOBAL_DELTA_CELL -INTEGER, PARAMETER :: LINSEARCH_LIMIT = 13 ! LINSEARCH_LIMIT-1 is the maximum size of array for linear search O(n). +INTEGER, PARAMETER :: LINSEARCH_LIMIT = 13 ! LINSEARCH_LIMIT-1 is the maximum size of array for linear search O(n). ! If Array larger -> binary search O(log(n)). ! Interpolation stencil threshold. Interpolation stencils will be defined if distance @@ -6492,9 +6492,10 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1,WC_B1 INTEGER :: IBOD, IWSEL, ICC, JCC -INTEGER :: IG, TRI, WSELEM(NOD1:NOD3), NOM, IIO, JJO, KKO, ICFACE, IOR, X1AXIS +INTEGER :: IG, TRI, WSELEM(NOD1:NOD3), NOM, IIO, JJO, KKO, ICFACE, IOR, X1AXIS, IWSEL_BACK REAL(EB):: XP(IAXIS:KAXIS),RDIR(IAXIS:KAXIS),V1(IAXIS:KAXIS),V2(IAXIS:KAXIS),V3(IAXIS:KAXIS),POS(IAXIS:KAXIS) -LOGICAL :: IS_INTERSECT=.FALSE. +REAL(EB):: POS_CUR(IAXIS:KAXIS),DIST_BACK,DIST_CUR,NBACK(IAXIS:KAXIS),DIAG_CELL +LOGICAL :: IS_INTERSECT=.FALSE.,IS_INTERSECT_CUR TYPE (SURFACE_TYPE), POINTER :: SF TYPE (WALL_TYPE), POINTER :: WC TYPE (MESH_TYPE), POINTER :: M @@ -6604,11 +6605,14 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ! Case of exposed Backing we need to find CFACE_INDEX of BACK CFACE. IF (SF%BACKING==EXPOSED .AND. SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN + NOM=0; IWSEL_BACK=0; IS_INTERSECT=.FALSE. IG = CF%BODTRI(1,IFACE) TRI = CF%BODTRI(2,IFACE) IF (VALID_GEOMETRY_FACE_DONOR(IG,TRI)) THEN XP(IAXIS:KAXIS) = (/ BC%X, BC%Y, BC%Z /) ! CFACE centroid location. RDIR(IAXIS:KAXIS)= - GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,TRI) ! Normal into the body. + ! Cast the inward ray against all other triangles and keep the CLOSEST forward intersection. + DIST_BACK = 1._EB/TWENTY_EPSILON_EB TRI_LOOP : DO IWSEL=1,GEOMETRY(IG)%N_FACES IF (IWSEL==TRI) CYCLE WSELEM(NOD1:NOD3) = GEOMETRY(IG)%FACES(NODS_WSEL*(IWSEL-1)+1:NODS_WSEL*IWSEL) @@ -6617,44 +6621,50 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, V2(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD2)-1)+1:MAX_DIM*WSELEM(NOD2)) V3(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD3)-1)+1:MAX_DIM*WSELEM(NOD3)) - ! Search for intersection point in POS(IAXIS:KAXIS): - CALL RAY_TRIANGLE_INTERSECT_PT(V1,V2,V3,XP,RDIR,IS_INTERSECT,POS) + ! Only consider EXIT faces (where the inward ray leaves the solid). + IF (DOT_PRODUCT(RDIR(IAXIS:KAXIS),GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)) <= TWENTY_EPSILON_EB) CYCLE TRI_LOOP - IF (IS_INTERSECT) EXIT TRI_LOOP + ! Search for intersection point in POS_CUR(IAXIS:KAXIS): + CALL RAY_TRIANGLE_INTERSECT_PT(V1,V2,V3,XP,RDIR,IS_INTERSECT_CUR,POS_CUR) + + IF (IS_INTERSECT_CUR) THEN + DIST_CUR = (POS_CUR(IAXIS)-XP(IAXIS))**2 + (POS_CUR(JAXIS)-XP(JAXIS))**2 + & + (POS_CUR(KAXIS)-XP(KAXIS))**2 + IF (DIST_CUR < DIST_BACK) THEN + DIST_BACK = DIST_CUR + POS = POS_CUR + IS_INTERSECT = .TRUE. + IWSEL_BACK = IWSEL ! Back GEOMETRY triangle the ray struck; its normal orients the back CFACE search. + ENDIF + ENDIF ENDDO TRI_LOOP - ELSE - IS_INTERSECT = .FALSE. ENDIF + ! Link a back CFACE only when an intersection was found and the distance in excess of the SURF THICKNESS + ! is within one cell diagonal; for longer distances (or no hit) leave the backing 'VOID'. IF (IS_INTERSECT) THEN - - ! Check that the distance in excess of the SURF THICKNESS is less than the cell diagonal size: - ! For longer distances from CFACE to BACK CFACE BC is 'VOID'. - IF(NORM2(XP-POS) - SF%THICKNESS > SQRT(DX(BC%IIG)**2 + DY(BC%JJG)**2 + DZ(BC%KKG)**2)) RETURN - - ! We Found an intersection with IWSEL in position POS(IAXIS:KAXIS): - ! Find indexes and mesh of cell containing intersection point: - CALL SEARCH_OTHER_MESHES(POS(IAXIS),POS(JAXIS),POS(KAXIS),NOM,IIO,JJO,KKO) - - ! Intersection point lies outside of the computational domain: treat backing as VOID. - IF (NOM==0) RETURN - - ! If the back CFACE lies in a different mesh than the front CFACE, defer to INITIALIZE_BACK_CFACE_EXCHANGE. - IF (NOM/=NM) THEN - IF (CFA%OD_INDEX>0) CALL ADD_BACK_CFACE_QUERY(NM,NOM,CFACE_INDEX,POS) - RETURN - ENDIF - - ! Back CFACE is in the same (local) mesh: resolve its index now, excluding the front CFACE itself. - ICFACE = GET_BACK_CFACE_INDEX(NM,IIO,JJO,KKO,POS,ICF,IFACE) - IF (ICFACE>0 .AND. CFA%OD_INDEX>0) THEN - M%BOUNDARY_ONE_D(CFA%OD_INDEX)%BACK_MESH = NM - M%BOUNDARY_ONE_D(CFA%OD_INDEX)%BACK_INDEX = ICFACE + DIAG_CELL = SQRT(DX(BC%IIG)**2 + DY(BC%JJG)**2 + DZ(BC%KKG)**2) + IF (NORM2(XP-POS) - SF%THICKNESS <= DIAG_CELL) THEN + ! We Found an intersection with IWSEL in position POS(IAXIS:KAXIS): + ! Back GEOMETRY triangle normal used to orient the back CFACE search (same-mesh and cross-mesh): + NBACK(IAXIS:KAXIS) = GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,IWSEL_BACK) + + ! Find indexes and mesh of cell containing intersection point: + CALL SEARCH_OTHER_MESHES(POS(IAXIS),POS(JAXIS),POS(KAXIS),NOM,IIO,JJO,KKO) + + IF (NOM/=0 .AND. NOM/=NM) THEN + ! Back CFACE lies in a different mesh: defer to INITIALIZE_BACK_CFACE_EXCHANGE. + IF (CFA%OD_INDEX>0) CALL ADD_BACK_CFACE_QUERY(NM,NOM,CFACE_INDEX,POS,NBACK) + ELSEIF (NOM==NM) THEN + ! Back CFACE is in the same (local) mesh: resolve its index now, excluding the front CFACE itself. + ICFACE = GET_BACK_CFACE_INDEX(NM,IIO,JJO,KKO,POS,ICF,IFACE,NBACK) + IF (ICFACE>0 .AND. CFA%OD_INDEX>0) THEN + M%BOUNDARY_ONE_D(CFA%OD_INDEX)%BACK_MESH = NM + M%BOUNDARY_ONE_D(CFA%OD_INDEX)%BACK_INDEX = ICFACE + ENDIF + ENDIF ENDIF - - ELSE ! Did not find intersection with other triangles. Leave VOID BC. - RETURN ENDIF ENDIF @@ -6743,58 +6753,76 @@ END SUBROUTINE INIT_CFACE_CELL ! ----------------------- ADD_BACK_CFACE_QUERY ----------------------------- -SUBROUTINE ADD_BACK_CFACE_QUERY(NM,NOM,FRONT_CFACE,POS) +SUBROUTINE ADD_BACK_CFACE_QUERY(NM,NOM,FRONT_CFACE,POS,NBACK) INTEGER, INTENT(IN) :: NM,NOM,FRONT_CFACE -REAL(EB), INTENT(IN) :: POS(IAXIS:KAXIS) +REAL(EB), INTENT(IN) :: POS(IAXIS:KAXIS),NBACK(IAXIS:KAXIS) INTEGER :: N,NDIM -REAL(EB), ALLOCATABLE :: XYZ_DUMMY(:,:) +REAL(EB), ALLOCATABLE :: XYZ_DUMMY(:,:),NBACK_DUMMY(:,:) INTEGER, ALLOCATABLE :: FRONT_DUMMY(:) ASSOCIATE(M3=>MESHES(NM)%OMESH(NOM)) IF (.NOT.ALLOCATED(M3%CFACE_QUERY_XYZ)) THEN M3%N_CFACE_QUERY_DIM = 64 ALLOCATE(M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,M3%N_CFACE_QUERY_DIM)) + ALLOCATE(M3%CFACE_QUERY_NBACK(IAXIS:KAXIS,M3%N_CFACE_QUERY_DIM)) ALLOCATE(M3%CFACE_QUERY_FRONT(M3%N_CFACE_QUERY_DIM)) M3%N_CFACE_QUERY = 0 ENDIF IF (M3%N_CFACE_QUERY+1 > M3%N_CFACE_QUERY_DIM) THEN ! Grow the query arrays. NDIM = M3%N_CFACE_QUERY_DIM ALLOCATE(XYZ_DUMMY(IAXIS:KAXIS,NDIM)) ; XYZ_DUMMY = M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,1:NDIM) + ALLOCATE(NBACK_DUMMY(IAXIS:KAXIS,NDIM)); NBACK_DUMMY = M3%CFACE_QUERY_NBACK(IAXIS:KAXIS,1:NDIM) ALLOCATE(FRONT_DUMMY(NDIM)) ; FRONT_DUMMY = M3%CFACE_QUERY_FRONT(1:NDIM) - DEALLOCATE(M3%CFACE_QUERY_XYZ,M3%CFACE_QUERY_FRONT) + DEALLOCATE(M3%CFACE_QUERY_XYZ,M3%CFACE_QUERY_NBACK,M3%CFACE_QUERY_FRONT) M3%N_CFACE_QUERY_DIM = 2*NDIM ALLOCATE(M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,M3%N_CFACE_QUERY_DIM)) + ALLOCATE(M3%CFACE_QUERY_NBACK(IAXIS:KAXIS,M3%N_CFACE_QUERY_DIM)) ALLOCATE(M3%CFACE_QUERY_FRONT(M3%N_CFACE_QUERY_DIM)) - M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,1:NDIM) = XYZ_DUMMY - M3%CFACE_QUERY_FRONT(1:NDIM) = FRONT_DUMMY - DEALLOCATE(XYZ_DUMMY,FRONT_DUMMY) + M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,1:NDIM) = XYZ_DUMMY + M3%CFACE_QUERY_NBACK(IAXIS:KAXIS,1:NDIM) = NBACK_DUMMY + M3%CFACE_QUERY_FRONT(1:NDIM) = FRONT_DUMMY + DEALLOCATE(XYZ_DUMMY,NBACK_DUMMY,FRONT_DUMMY) ENDIF N = M3%N_CFACE_QUERY + 1 M3%N_CFACE_QUERY = N -M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,N) = POS(IAXIS:KAXIS) -M3%CFACE_QUERY_FRONT(N) = FRONT_CFACE +M3%CFACE_QUERY_XYZ(IAXIS:KAXIS,N) = POS(IAXIS:KAXIS) +M3%CFACE_QUERY_NBACK(IAXIS:KAXIS,N) = NBACK(IAXIS:KAXIS) +M3%CFACE_QUERY_FRONT(N) = FRONT_CFACE END ASSOCIATE END SUBROUTINE ADD_BACK_CFACE_QUERY ! ----------------------- GET_BACK_CFACE_INDEX ----------------------------- -INTEGER FUNCTION GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,ICF_EXCLUDE,IFACE_EXCLUDE) +INTEGER FUNCTION GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,ICF_EXCLUDE,IFACE_EXCLUDE,NBACK) INTEGER, INTENT(IN) :: NOM,IIO,JJO,KKO,ICF_EXCLUDE,IFACE_EXCLUDE REAL(EB), INTENT(IN) :: POS(IAXIS:KAXIS) -INTEGER :: IIV(3),JJV(3),KKV(3),III,JJJ,KKK,II,JJ,KK,ICF2,JCF22,ICFF,JCFF,ICF3,JCF3 -REAL(EB) :: DIST,DIST2 -LOGICAL :: BACK_CFACE_FOUND +REAL(EB), INTENT(IN), OPTIONAL :: NBACK(IAXIS:KAXIS) ! Normal of the back GEOMETRY triangle the ray struck, used to filter orientation. +INTEGER :: IIV(3),JJV(3),KKV(3),III,JJJ,KKK,II,JJ,KK,ICF2,JCF22,ICFF,JCFF,ICF3,JCF3,ICFF_ANY,JCFF_ANY,ICAND,IWCAND +REAL(EB) :: DIST,DIST2,DIST_ANY,RD(IAXIS:KAXIS),RDN,ALIGN +LOGICAL :: BACK_CFACE_FOUND,BACK_CFACE_FOUND_ANY,USE_ORIENT +REAL(EB), PARAMETER :: ALIGN_MIN = 0.5_EB ! Accept a back CFACE only if n_cand . n_backtri >= this val. GET_BACK_CFACE_INDEX = 0 IF (NOM<1) RETURN IF (.NOT.ALLOCATED(MESHES(NOM)%CCVAR)) RETURN +! Enable the orientation filter only when the back triangle normal is provided. +USE_ORIENT = .FALSE. +IF (PRESENT(NBACK)) THEN + RDN = NORM2(NBACK) + IF (RDN > TWENTY_EPSILON_EB) THEN + RD(IAXIS:KAXIS) = NBACK(IAXIS:KAXIS)/RDN + USE_ORIENT = .TRUE. + ENDIF +ENDIF + IIV(1:3) = (/ IIO, MAX(IIO-1,1), MIN(IIO+1,MESHES(NOM)%IBAR) /) JJV(1:3) = (/ JJO, MAX(JJO-1,1), MIN(JJO+1,MESHES(NOM)%JBAR) /) KKV(1:3) = (/ KKO, MAX(KKO-1,1), MIN(KKO+1,MESHES(NOM)%KBAR) /) -DIST = 1._EB/TWENTY_EPSILON_EB; ICFF=0; JCFF=0; BACK_CFACE_FOUND=.FALSE. +DIST = 1._EB/TWENTY_EPSILON_EB; ICFF=0; JCFF=0; BACK_CFACE_FOUND=.FALSE. ! Best orientation-consistent pick. +DIST_ANY = 1._EB/TWENTY_EPSILON_EB; ICFF_ANY=0; JCFF_ANY=0; BACK_CFACE_FOUND_ANY=.FALSE. ! Nearest-centroid fallback. DO KKK=1,3 KK=KKV(KKK) DO JJJ=1,3 @@ -6809,7 +6837,24 @@ INTEGER FUNCTION GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,ICF_EXCLUDE,IFACE_EXCL DIST2 = (POS(IAXIS) - MESHES(NOM)%CUT_FACE(ICF2)%XYZCEN(IAXIS,JCF22))**2._EB + & (POS(JAXIS) - MESHES(NOM)%CUT_FACE(ICF2)%XYZCEN(JAXIS,JCF22))**2._EB + & (POS(KAXIS) - MESHES(NOM)%CUT_FACE(ICF2)%XYZCEN(KAXIS,JCF22))**2._EB - IF (DIST2=1 .AND. IWCAND>=1) THEN + IF (IWCAND<=GEOMETRY(ICAND)%N_FACES) & + ALIGN = DOT_PRODUCT(GEOMETRY(ICAND)%FACES_NORMAL(IAXIS:KAXIS,IWCAND),RD(IAXIS:KAXIS)) + ENDIF + ENDIF + ENDIF + ! Nearest-centroid candidate (used if no orientation-consistent face is found): + IF (DIST2=ALIGN_MIN) THEN DIST = DIST2; ICFF = ICF2; JCFF = JCF22; BACK_CFACE_FOUND = .TRUE. ENDIF ENDDO @@ -6817,6 +6862,10 @@ INTEGER FUNCTION GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,ICF_EXCLUDE,IFACE_EXCL ENDDO ENDDO ENDDO +! Fall back to the nearest-centroid CFACE if no orientation-consistent back face was found. +IF (.NOT.BACK_CFACE_FOUND) THEN + ICFF = ICFF_ANY; JCFF = JCFF_ANY; BACK_CFACE_FOUND = BACK_CFACE_FOUND_ANY +ENDIF ! Loop NOM CUT_FACE array to find BACKING CFACE index: IF (BACK_CFACE_FOUND) THEN DO ICF3=1,MESHES(NOM)%N_CUTFACE_MESH diff --git a/Source/radi.f90 b/Source/radi.f90 index c379b154b5..997773d1f1 100644 --- a/Source/radi.f90 +++ b/Source/radi.f90 @@ -2774,7 +2774,7 @@ MODULE RAD REAL(EB) :: KAPPA_C_FAC,LOG_KAPPA_C_FAC !< Scaling factor for soot coefficient look-up table REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: Y2RADCAL_SPECIES !< Primitive species mapping to radcal species REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: RADCAL_SPECIES2KAPPA ! Absorption coefficient look-up table -REAL(EB), ALLOCATABLE, DIMENSION(:) :: KAPPA_COND !< Array to remove condensed species from absorption coefficient look-up +REAL(EB), ALLOCATABLE, DIMENSION(:) :: KAPPA_COND !< Array to remove condensed species from absorption coefficient look-up INTEGER :: N_RADCAL_ARRAY_SIZE !< Number of radcal species present INTEGER :: RADCAL_SPECIES_INDEX(16) !< Mapping of radcal species present to radcal calling function CHARACTER(LABEL_LENGTH) :: RADCAL_SPECIES_ID(16)='NULL'!< Name of radcal species @@ -3093,7 +3093,7 @@ SUBROUTINE INIT_RADIATION END DO GET_RADCAL_SPECIES BUILD_KAPPA_ARRAY: IF (N_RADCAL_ARRAY_SIZE>0) THEN - + KAPPA_X_FAC = (KAPPA_X_MAX/KAPPA_X_MIN)**(1._EB/(N_KAPPA_X-1)) LOG_KAPPA_X_FAC = LOG(KAPPA_X_FAC) KAPPA_C_FAC = (KAPPA_C_MAX/KAPPA_C_MIN)**(1._EB/(N_KAPPA_X-1)) @@ -3106,7 +3106,7 @@ SUBROUTINE INIT_RADIATION ALLOCATE(KAPPA_COND(N_TRACKED_SPECIES),STAT=IZERO) CALL ChkMemErr('RADI','KAPPA_COND',IZERO) KAPPA_COND = 1._EB - + DO NS=1,N_TRACKED_SPECIES IF (SPECIES_MIXTURE(NS)%EVAPORATION_SMIX_INDEX>0) THEN CALL MEAN_CROSS_SECTIONS(SMIX_INDEX=NS) @@ -3339,7 +3339,7 @@ SUBROUTINE CALCULATE_FVM_ANGLES() ELSE MERI_COMP(N)=(SIN(PHIUP)-SIN(PHILOW))*F_THETA AZIM_COMP(N)=(COS(PHILOW)-COS(PHIUP))*F_THETA - AXIS_COMP(N)=0.5_EB*(PHIUP-PHILOW) * ((SIN(THETAUP))**2-(SIN(THETALOW))**2) + AXIS_COMP(N)=0.5_EB*(PHIUP-PHILOW) * ((SIN(THETAUP))**2-(SIN(THETALOW))**2) DLANG_LOCAL(1,N) = SIN(THETA)*COS(PHI) DLANG_LOCAL(2,N) = SIN(THETA)*SIN(PHI) DLANG_LOCAL(3,N) = COS(THETA) @@ -3628,7 +3628,7 @@ SUBROUTINE INTERPOLATE_IL() ENDIF ENDDO - !Interpolate neighbouring mesh cell intensities + !Interpolate neighbouring mesh cell intensities DO NNN=1,N_NEIGHBORING_MESHES NOM = NEIGHBORING_MESH(NNN) M2 => OMESH(NOM) @@ -3710,11 +3710,12 @@ SUBROUTINE RADIATION_FVM USE MATH_FUNCTIONS, ONLY : INTERPOLATE1D USE TRAN, ONLY : GET_IJK USE COMPLEX_GEOMETRY, ONLY : CC_CGSC,CC_SOLID +USE CC_SCALARS, ONLY : GET_CFACE_RAD_NVEC USE PHYSICAL_FUNCTIONS, ONLY : GET_VOLUME_FRACTION, GET_MASS_FRACTION 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,SUMILW + RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC,SUMILW,NVECL(3) INTEGER :: N,NN,IIG,JJG,KKG,I,J,K,IW,ICF,II,JJ,KK,IOR,IC,IWUP,IWDOWN, & ISTART, IEND, ISTEP, JSTART, JEND, JSTEP, & @@ -4188,9 +4189,10 @@ SUBROUTINE RADIATION_FVM CFA => CFACE(ICF) BR => BOUNDARY_RADIA(CFA%BR_INDEX) BC => BOUNDARY_COORD(CFA%BC_INDEX) + NVECL = GET_CFACE_RAD_NVEC(CFA) DO N=1,NRA DLA = (/DLX(N),DLY(N),DLZ(N)/) - DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(NVECL,DLA) ! face normal * radiation angle IF (DLF<0._EB) INRAD_F(ICF) = INRAD_F(ICF) - DLF*BR%BAND(IBND)%ILW(N) ENDDO ENDDO @@ -4294,7 +4296,7 @@ SUBROUTINE RADIATION_FVM BR => BOUNDARY_RADIA(CFA%BR_INDEX) BC => BOUNDARY_COORD(CFA%BC_INDEX) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) - DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(GET_CFACE_RAD_NVEC(CFA),DLA) ! face normal * radiation angle IF (DLF<0._EB) CYCLE CFACE_LOOP1 BR%BAND(IBND)%ILW(N) = OUTRAD_F(ICF) + RPI*(1._EB-B1%EMISSIVITY)*INRAD_F(ICF) ENDDO CFACE_LOOP1 @@ -4472,7 +4474,7 @@ SUBROUTINE RADIATION_FVM DO IFACE=1,CF%NFACE CFA => CFACE(CF%CFACE_INDEX(IFACE)) BC => BOUNDARY_COORD(CFA%BC_INDEX) - DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(GET_CFACE_RAD_NVEC(CFA),DLA) ! face normal * radiation angle IF (DLF>0._EB) THEN ! upwind BR => BOUNDARY_RADIA(CFA%BR_INDEX) AILFU = AILFU + DLF*CFA%AREA*BR%BAND(IBND)%ILW(N) @@ -4595,7 +4597,7 @@ SUBROUTINE RADIATION_FVM CFA => CFACE(ICF) IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP2 BC => BOUNDARY_COORD(CFA%BC_INDEX) - DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(GET_CFACE_RAD_NVEC(CFA),DLA) ! face normal * radiation angle IF (DLF>=0._EB) CYCLE CFACE_LOOP2 ! outgoing BR => BOUNDARY_RADIA(CFA%BR_INDEX) BC => BOUNDARY_COORD(CFA%BC_INDEX) @@ -4607,7 +4609,7 @@ SUBROUTINE RADIATION_FVM CFA => CFACE(ICF) IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP3 BC => BOUNDARY_COORD(CFA%BC_INDEX) - DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(GET_CFACE_RAD_NVEC(CFA),DLA) ! face normal * radiation angle IF (DLF>=0._EB) CYCLE CFACE_LOOP3 ! outgoing BR => BOUNDARY_RADIA(CFA%BR_INDEX) BC => BOUNDARY_COORD(CFA%BC_INDEX) diff --git a/Source/type.f90 b/Source/type.f90 index a9e102bd44..392a779432 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -1043,6 +1043,7 @@ MODULE TYPES TYPE (STORAGE_TYPE) :: CFACE_SEND_BUFFER,CFACE_RECV_BUFFER INTEGER :: N_CFACE_QUERY=0,N_CFACE_QUERY_DIM=0 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: CFACE_QUERY_XYZ !< Back-face intersection points (3,N) that NM needs NOM to resolve + REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: CFACE_QUERY_NBACK !< Back GEOMETRY triangle normal (3,N) used to orient the back CFACE search INTEGER, ALLOCATABLE, DIMENSION(:) :: CFACE_QUERY_FRONT !< Front CFACE index (in NM) associated with each query ! CC_IBM data exchange arrays: @@ -2158,7 +2159,7 @@ MODULE TYPES CHARACTER(LABEL_LENGTH) :: SMOKEVIEW_LABEL !< Smokeview label for QUANTITY CHARACTER(LABEL_LENGTH) :: SMOKEVIEW_BAR_LABEL !< Smokeview colorbar label for QUANTITY CHARACTER(LABEL_LENGTH) :: UNITS !< Units for QUANTITY - LOGICAL :: DRY !< Remove water vapor before computing a mass or volume fraction + LOGICAL :: DRY !< Remove water vapor before computing a mass or volume fraction END TYPE HVAC_QUANTITY_TYPE