From 9e850ca38706ac4a6bdc656cb2b8ab38c7ae5135 Mon Sep 17 00:00:00 2001 From: mcgratta Date: Thu, 9 Nov 2023 14:12:41 -0500 Subject: [PATCH] FDS Source: Make NVEC (normal vector) a component of BOUNDARY_COORD --- Source/ccib.f90 | 24 ++++++++++++++---------- Source/dump.f90 | 8 ++++---- Source/func.f90 | 6 +++--- Source/geom.f90 | 16 ++++++++-------- Source/init.f90 | 18 +++++++++++++++++- Source/part.f90 | 42 +++++++++++++++++++++--------------------- Source/radi.f90 | 15 ++++++++++----- Source/soot.f90 | 2 +- Source/type.f90 | 3 ++- Source/vege.f90 | 2 +- 10 files changed, 81 insertions(+), 55 deletions(-) diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 66672245732..c77cc546769 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -3495,10 +3495,11 @@ SUBROUTINE CC_COMPUTE_KRES(APPLY_TO_ESTIMATED_VARIABLES,NM) AF = CUT_FACE(IFC2)%AREA(IFACE2) VELN = (1._EB-PRFCT)*CUT_FACE(IFC2)%VEL( IFACE2) + PRFCT *CUT_FACE(IFC2)%VELS(IFACE2) ! - to use velocity into gasphase, projected area. + BC => MESHES(NM)%BOUNDARY_COORD(CFACE(ICFA)%BC_INDEX) DO X1AXIS=IAXIS,KAXIS - AUI = ABS(CFACE(ICFA)%NVEC(X1AXIS)*AF) + AUI = ABS(BC%NVEC(X1AXIS)*AF) ATOTV(X1AXIS) = ATOTV(X1AXIS) + AUI - UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*CFACE(ICFA)%NVEC(X1AXIS)*AUI + UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*BC%NVEC(X1AXIS)*AUI ENDDO END SELECT @@ -3541,10 +3542,11 @@ SUBROUTINE CC_COMPUTE_KRES(APPLY_TO_ESTIMATED_VARIABLES,NM) AF = CUT_FACE(IFC2)%AREA(IFACE2) VELN = (1._EB-PRFCT)*CUT_FACE(IFC2)%VEL( IFACE2) + PRFCT *CUT_FACE(IFC2)%VELS(IFACE2) ! - to use velocity into gasphase, projected area. + BC => MESHES(NM)%BOUNDARY_COORD(CFACE(ICFA)%BC_INDEX) DO X1AXIS=IAXIS,KAXIS - AUI = ABS(CFACE(ICFA)%NVEC(X1AXIS)*AF) + AUI = ABS(BC%NVEC(X1AXIS)*AF) ATOTV(X1AXIS) = ATOTV(X1AXIS) + AUI - UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*CFACE(ICFA)%NVEC(X1AXIS)*AUI + UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*BC%NVEC(X1AXIS)*AUI ENDDO END SELECT @@ -5844,12 +5846,13 @@ SUBROUTINE CC_COMPUTE_VISCOSITY(DT,NM) DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS CFA => CFACE(ICF) + BC => BOUNDARY_COORD(CFA%BC_INDEX) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) B2 => BOUNDARY_PROP2(CFA%B2_INDEX) SF => SURFACE(CFA%SURF_INDEX) ! Surface Velocity: - NVEC = CFA%NVEC + NVEC = BC%NVEC ! right now VEL_T not defined for CFACEs TVEC1=(/ 0._EB,0._EB,0._EB/) @@ -6825,7 +6828,7 @@ SUBROUTINE CFACE_THERMAL_GASVARS(ICF,B1) VVEL(VIND) = VVEL(VIND) + CUT_FACE(IND1)%INT_COEF(INPE)*CUT_FACE(IND1)%INT_FVARS(INT_VEL_IND,INPE) ENDDO ENDDO - V_TANG(IAXIS:KAXIS) = VVEL(IAXIS:KAXIS) - DOT_PRODUCT(VVEL,CFA%NVEC)*CFA%NVEC(IAXIS:KAXIS) + V_TANG(IAXIS:KAXIS) = VVEL(IAXIS:KAXIS) - DOT_PRODUCT(VVEL,BC%NVEC)*BC%NVEC(IAXIS:KAXIS) B1%U_TANG = SQRT(V_TANG(IAXIS)**2._EB+V_TANG(JAXIS)**2._EB+V_TANG(KAXIS)**2._EB) ELSE @@ -15940,11 +15943,12 @@ SUBROUTINE GET_BOUND_VEL(X1AXIS,INBFC_CFCEN,XYZ_PP,VELX1) CFA => CFACE(ICF) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) +BC => BOUNDARY_COORD(CFA%BC_INDEX) ! Velocity into Gas Region, component along X1AXIS: IF (PREDICTOR) THEN - VELX1 = -B1%U_NORMAL * CFACE(ICF)%NVEC(X1AXIS) + VELX1 = -B1%U_NORMAL * BC%NVEC(X1AXIS) ELSE - VELX1 = -B1%U_NORMAL_S* CFACE(ICF)%NVEC(X1AXIS) + VELX1 = -B1%U_NORMAL_S* BC%NVEC(X1AXIS) ENDIF RETURN @@ -16583,7 +16587,7 @@ SUBROUTINE GET_CRTCFCC_INT_STENCILS N_CVAR_START = NPE_LIST_START DO IFACE=1,CUT_FACE(ICF)%NFACE ICF1 = CUT_FACE(ICF)%CFACE_INDEX(IFACE) - DV(IAXIS:KAXIS) = CFACE(ICF1)%NVEC(IAXIS:KAXIS) + DV(IAXIS:KAXIS) = MESHES(NM)%BOUNDARY_COORD(CFACE(ICF1)%BC_INDEX)%NVEC(IAXIS:KAXIS) CALL GET_DELN(1.001_EB,DELN,DXCELL(I),DYCELL(J),DZCELL(K),NVEC=DV,CLOSE_PT=.TRUE.) ! Origin in boundary point XYZ_PP(IAXIS:KAXIS): XYZ_PP(IAXIS:KAXIS) = CUT_FACE(ICF)%XYZCEN(IAXIS:KAXIS,IFACE) @@ -16611,7 +16615,7 @@ SUBROUTINE GET_CRTCFCC_INT_STENCILS XYZ_PP(IAXIS:KAXIS) = CUT_FACE(ICF)%XYZCEN(IAXIS:KAXIS,IFACE) FOUND_INBFC(1:3) = (/ CC_FTYPE_CFINB, ICF, IFACE /) ICF1 = CUT_FACE(ICF)%CFACE_INDEX(IFACE) - DV(IAXIS:KAXIS) = CFACE(ICF1)%NVEC(IAXIS:KAXIS) + DV(IAXIS:KAXIS) = MESHES(NM)%BOUNDARY_COORD(CFACE(ICF1)%BC_INDEX)%NVEC(IAXIS:KAXIS) INT_NOUT(IAXIS:KAXIS,IFACE) = DV(IAXIS:KAXIS) INT_XN(0:INT_N_EXT_PTS) = 0._EB ! Initialize interpolation coefficients along the normal probe direction DV diff --git a/Source/dump.f90 b/Source/dump.f90 index d7742874468..9d83216490a 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -5265,7 +5265,7 @@ SUBROUTINE DUMP_CFACES_GEOM(FUNIT,STIME) IND1 = CFACE(ICF)%CUT_FACE_IND1 IND2 = CFACE(ICF)%CUT_FACE_IND2 XLOC(1:3) = (/ CFACE(ICF)%X, CFACE(ICF)%Y, CFACE(ICF)%Z /) - XLOC(4:6) = XLOC(1:3) + CUT_FACE(IND1)%INT_XN(1,IND2)*CFACE(ICF)%NVEC(IAXIS:KAXIS) + XLOC(4:6) = XLOC(1:3) + CUT_FACE(IND1)%INT_XN(1,IND2)*BOUNDARY_COORD(CFACE(ICF)%BC_INDEX)%NVEC(IAXIS:KAXIS) VERTS(6*(ICF-INTERNAL_CFACE_CELLS_LB-1)+1:6*(ICF-INTERNAL_CFACE_CELLS_LB)) = REAL(XLOC(1:6),FB) ENDDO WRITE(FUNIT) VERTS(1:6*N_INTERNAL_CFACE_CELLS) @@ -9001,7 +9001,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(NM,INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL PVEC = ( PR1 - (PR2-PR1)*Z1**2 / (Z2**2-Z1**2) ) * NVEC ! surface normal pressure force ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN - NVEC = CFA%NVEC + NVEC = BC%NVEC ! find cut-cell adjacent to CFACE IND1 = CFA%CUT_FACE_IND1 IND2 = CFA%CUT_FACE_IND2 @@ -9038,7 +9038,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(NM,INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL W_CELL = 0.5_EB*(W(IIG,JJG,KKG-1)+W(IIG,JJG,KKG)) MU_WALL = MU_DNS(IIG,JJG,KKG) ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN - NVEC = CFA%NVEC + NVEC = BC%NVEC ! right now VEL_T not defined for CFACEs TVEC1=(/ 0._EB,0._EB,0._EB/) TVEC2=(/ 0._EB,0._EB,0._EB/) @@ -9093,7 +9093,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(NM,INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL KKG = BC%KKG PVEC = RHO(IIG,JJG,KKG)*H(IIG,JJG,KKG) * NVEC ! surface normal pressure force ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN - NVEC = CFA%NVEC + NVEC = BC%NVEC ! find cut-cell adjacent to CFACE IND1 = CFA%CUT_FACE_IND1 IND2 = CFA%CUT_FACE_IND2 diff --git a/Source/func.f90 b/Source/func.f90 index af90263d0b3..282fea51a1e 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -1638,9 +1638,6 @@ SUBROUTINE PACK_CFACE(NM,OS,CFA,SURF_INDEX,RC,IC,LC,UNPACK_IT,COUNT_ONLY) ! Assign and initialize reals RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%AREA,UNPACK_IT) -RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%NVEC(1),UNPACK_IT) -RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%NVEC(2),UNPACK_IT) -RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%NVEC(3),UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%DUNDT,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%RSUM_G,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%TMP_G,UNPACK_IT) @@ -1701,6 +1698,9 @@ SUBROUTINE PACK_BOUNDARY_COORD(NM,IC,RC,OS,BC_INDEX,UNPACK_IT,COUNT_ONLY) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%Y2,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%Z1,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%Z2,UNPACK_IT) +RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%NVEC(1),UNPACK_IT) +RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%NVEC(2),UNPACK_IT) +RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%NVEC(3),UNPACK_IT) END SUBROUTINE PACK_BOUNDARY_COORD diff --git a/Source/geom.f90 b/Source/geom.f90 index 9308c737754..318a9ce4a08 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -669,7 +669,7 @@ SUBROUTINE POINT_IN_CFACE(NM,XP,YP,ZP,CFACE_INDEX,IN_CFACE) ! Normal, max normal component, define plane X2AXIS,X3AXIS to do search: VERT_CUTFACE = SIZE(MESHES(NM)%CUT_FACE(INBFC)%CFELEM, DIM=1); ALLOCATE(CFELEM(1:VERT_CUTFACE)) CFELEM(1:VERT_CUTFACE) = MESHES(NM)%CUT_FACE(INBFC)%CFELEM(1:VERT_CUTFACE,INBFC_LOC) -NVEC(IAXIS:KAXIS) => CFACE(CFACE_INDEX)%NVEC(IAXIS:KAXIS) +NVEC(IAXIS:KAXIS) => MESHES(NM)%BOUNDARY_COORD(CFACE(CFACE_INDEX)%BC_INDEX)%NVEC(IAXIS:KAXIS) ! Plane equation for INBOUNDARY cut-face plane: ! Location of first point in cf polygon is P0: @@ -8516,13 +8516,13 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ! Normal to cut-face: V2(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(2,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE) V3(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(3,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE) - CALL CROSS_PRODUCT(CFA%NVEC(IAXIS:KAXIS),V2,V3) - IF(NORM2(CFA%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN - CFA%NVEC(IAXIS:KAXIS) = CFA%NVEC(IAXIS:KAXIS)/NORM2(CFA%NVEC) + CALL CROSS_PRODUCT(BC%NVEC(IAXIS:KAXIS),V2,V3) + IF(NORM2(BC%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN + BC%NVEC(IAXIS:KAXIS) = BC%NVEC(IAXIS:KAXIS)/NORM2(BC%NVEC) ELSE IBOD =CF%BODTRI(1,IFACE) IWSEL=CF%BODTRI(2,IFACE) - CFA%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) + BC%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) ENDIF ! Boundary CFACES processed are defined of type SOLID_BOUNDARY @@ -8544,9 +8544,9 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, WC_BC => M%BOUNDARY_COORD(WC%BC_INDEX) IOR = WC_BC%IOR SELECT CASE(ABS(IOR)) - CASE(IAXIS); CFA%NVEC(IAXIS:KAXIS) = (/ REAL(SIGN(1,IOR),EB), 0._EB, 0._EB /) - CASE(JAXIS); CFA%NVEC(IAXIS:KAXIS) = (/ 0._EB, REAL(SIGN(1,IOR),EB), 0._EB /) - CASE(KAXIS); CFA%NVEC(IAXIS:KAXIS) = (/ 0._EB, 0._EB, REAL(SIGN(1,IOR),EB) /) + CASE(IAXIS); BC%NVEC(IAXIS:KAXIS) = (/ REAL(SIGN(1,IOR),EB), 0._EB, 0._EB /) + CASE(JAXIS); BC%NVEC(IAXIS:KAXIS) = (/ 0._EB, REAL(SIGN(1,IOR),EB), 0._EB /) + CASE(KAXIS); BC%NVEC(IAXIS:KAXIS) = (/ 0._EB, 0._EB, REAL(SIGN(1,IOR),EB) /) END SELECT BC%IOR = IOR diff --git a/Source/init.f90 b/Source/init.f90 index d041493add1..22baf902974 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -1140,7 +1140,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) IW = MAXLOC(M%CUT_FACE(ICF)%AREA(1:M%CUT_FACE(ICF)%NFACE),DIM=1) CFA => M%CFACE( M%CUT_FACE(ICF)%CFACE_INDEX(IW) ) BC => M%BOUNDARY_COORD(CFA%BC_INDEX) - IF (CFA%NVEC(KAXIS)>-TWO_EPSILON_EB .AND. CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + IF (BC%NVEC(KAXIS)>-TWO_EPSILON_EB .AND. CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN IF (BC%KKG > M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)) THEN M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE) = MIN( M%KBAR , M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)+BC%KKG ) ENDIF @@ -3142,6 +3142,14 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT) BC%JJG = JJG BC%KKG = KKG BC%IOR = IOR +SELECT CASE(BC%IOR) + CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/) + CASE(-1) ; BC%NVEC=(/-1._EB, 0._EB, 0._EB/) + CASE( 2) ; BC%NVEC=(/ 0._EB, 1._EB, 0._EB/) + CASE(-2) ; BC%NVEC=(/ 0._EB,-1._EB, 0._EB/) + CASE( 3) ; BC%NVEC=(/ 0._EB, 0._EB, 1._EB/) + CASE(-3) ; BC%NVEC=(/ 0._EB, 0._EB,-1._EB/) +END SELECT BC%X = XW BC%Y = YW BC%Z = ZW @@ -3510,6 +3518,14 @@ SUBROUTINE INIT_THIN_WALL_CELL(NM,I,J,K,OBST_INDEX,ITW,IOR,SURF_INDEX,IEC) BC%JJG = J BC%KKG = K BC%IOR = IOR +SELECT CASE(BC%IOR) + CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/) + CASE(-1) ; BC%NVEC=(/-1._EB, 0._EB, 0._EB/) + CASE( 2) ; BC%NVEC=(/ 0._EB, 1._EB, 0._EB/) + CASE(-2) ; BC%NVEC=(/ 0._EB,-1._EB, 0._EB/) + CASE( 3) ; BC%NVEC=(/ 0._EB, 0._EB, 1._EB/) + CASE(-3) ; BC%NVEC=(/ 0._EB, 0._EB,-1._EB/) +END SELECT BC%X = M%X(I) BC%Y = M%Y(J) BC%Z = M%Z(K) diff --git a/Source/part.f90 b/Source/part.f90 index f585d1d80c3..38e319df839 100644 --- a/Source/part.f90 +++ b/Source/part.f90 @@ -825,17 +825,17 @@ SUBROUTINE PARTICLE_FACE_INSERT(WALL_INDEX,CFACE_INDEX) ENDIF ELSEIF (PRESENT(CFACE_INDEX)) THEN CALL RANDOM_CFACE_XYZ(NM,CFA,CFA_X,CFA_Y,CFA_Z) - BC%X = CFA_X + CFA%NVEC(1)*VENT_OFFSET*DX(IIG) - BC%Y = CFA_Y + CFA%NVEC(2)*VENT_OFFSET*DY(JJG) - BC%Z = CFA_Z + CFA%NVEC(3)*VENT_OFFSET*DZ(KKG) + BC%X = CFA_X + BC%NVEC(1)*VENT_OFFSET*DX(IIG) + BC%Y = CFA_Y + BC%NVEC(2)*VENT_OFFSET*DY(JJG) + BC%Z = CFA_Z + BC%NVEC(3)*VENT_OFFSET*DZ(KKG) IF (INSERT_TYPE==1 .AND. LP%EMBER) THEN CALL RANDOM_NUMBER(RN3) BC%Z = CFA_Z + SF%EMBER_GENERATION_HEIGHT(1) + & (SF%EMBER_GENERATION_HEIGHT(2)-SF%EMBER_GENERATION_HEIGHT(1))*REAL(RN3,EB) ENDIF - LP%U = DOT_PRODUCT(CFA%NVEC,(/-B1%U_NORMAL,SF%VEL_T(1),SF%VEL_T(2)/)) - LP%V = DOT_PRODUCT(CFA%NVEC,(/SF%VEL_T(1),-B1%U_NORMAL,SF%VEL_T(2)/)) - LP%W = DOT_PRODUCT(CFA%NVEC,(/SF%VEL_T(1),SF%VEL_T(2),-B1%U_NORMAL/)) + LP%U = DOT_PRODUCT(BC%NVEC,(/-B1%U_NORMAL,SF%VEL_T(1),SF%VEL_T(2)/)) + LP%V = DOT_PRODUCT(BC%NVEC,(/SF%VEL_T(1),-B1%U_NORMAL,SF%VEL_T(2)/)) + LP%W = DOT_PRODUCT(BC%NVEC,(/SF%VEL_T(1),SF%VEL_T(2),-B1%U_NORMAL/)) ENDIF WALL_OR_CFACE_IF_2 ! Update idicies in case offset puts location in a different cell @@ -1146,7 +1146,7 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES IF (DIST TWO_EPSILON_EB) CC_VALID=.TRUE. @@ -1319,7 +1319,7 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES IF (DIST TWO_EPSILON_EB) EXIT RAND_LOCATION_LOOP @@ -1994,16 +1994,16 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM) CFA2 => CFACE(ICF) CFA2_BC => BOUNDARY_COORD(CFA2%BC_INDEX) - DIST2 = DOT_PRODUCT(CFA%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB)) + DIST2 = DOT_PRODUCT(CFA_BC%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB)) IF(DIST2<-0.99_EB) THEN ! Only for slopes less than 8 degrees. IN_CFACE=.TRUE. ! Test for case of CFACE_INDEX switching to an ICF with similar slope, if so don't do anything. - IF (.NOT.(LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA%NVEC,CFA2%NVEC))>0.99_EB) ) THEN + IF (.NOT.(LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA_BC%NVEC,CFA2_BC%NVEC))>0.99_EB) ) THEN - IF (LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA%NVEC,CFA2%NVEC))<=0.99_EB) THEN + IF (LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA_BC%NVEC,CFA2_BC%NVEC))<=0.99_EB) THEN ! Case of switching to different ICF with different slope: - DIST2 = DOT_PRODUCT(CFA2%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB)) + DIST2 = DOT_PRODUCT(CFA2_BC%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB)) ! If ICF is almost vertical pointing in the direction of velocity, set creep velocity. IF(ABS(DIST2)<0.01_EB .AND. (CFA2_BC%Z CFACE(ICF) CFA_BC => BOUNDARY_COORD(CFA%BC_INDEX) P_VECTOR = (/BC%X-CFA_BC%X, BC%Y-CFA_BC%Y, BC%Z-CFA_BC%Z/) - TEST_POS = .FALSE.; IF (LP%CFACE_INDEX == 0) TEST_POS = DOT_PRODUCT(CFA%NVEC,P_VECTOR) > TWO_EPSILON_EB + TEST_POS = .FALSE.; IF (LP%CFACE_INDEX == 0) TEST_POS = DOT_PRODUCT(CFA_BC%NVEC,P_VECTOR) > TWO_EPSILON_EB - CFACE_ATTACH : IF (DOT_PRODUCT(CFA%NVEC,GVEC)>0._EB .OR. TEST_POS) THEN + CFACE_ATTACH : IF (DOT_PRODUCT(CFA_BC%NVEC,GVEC)>0._EB .OR. TEST_POS) THEN ! Normal points down or particle in gas phase. Let particle move freely: @@ -2071,8 +2071,8 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM) IF(LPC%ADHERE_TO_SOLID) THEN - CALL CROSS_PRODUCT(VEL_VECTOR_1,CFA%NVEC,GVEC) - CALL CROSS_PRODUCT(VEL_VECTOR_2,VEL_VECTOR_1,CFA%NVEC) + CALL CROSS_PRODUCT(VEL_VECTOR_1,CFA_BC%NVEC,GVEC) + CALL CROSS_PRODUCT(VEL_VECTOR_2,VEL_VECTOR_1,CFA_BC%NVEC) CFACE_SLOPE : IF (NORM2(VEL_VECTOR_2) > TWO_EPSILON_EB .AND. ABS(LP%W) > TWO_EPSILON_EB) THEN ! The surface is tilted; particles go down slope: VEL_VECTOR_1 = VEL_VECTOR_2/NORM2(VEL_VECTOR_2) @@ -2092,7 +2092,7 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM) ELSEIF (LP%CFACE_INDEX /= 0 .AND. ABS(LP%W) TWO_EPSILON_EB) DIST2_MIN = ATAN2(LP%V,LP%U) @@ -2117,12 +2117,12 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM) PVEC_L = NORM2(P_VECTOR) IF (PVEC_L>TWO_EPSILON_EB) THEN - THETA = ACOS(DOT_PRODUCT(CFACE(ICF)%NVEC,P_VECTOR/PVEC_L)) + THETA = ACOS(DOT_PRODUCT(BOUNDARY_COORD(CFACE(ICF)%BC_INDEX)%NVEC,P_VECTOR/PVEC_L)) IF (THETA>PIO2) THEN DELTA = PVEC_L*SIN(THETA-0.5_EB*PI)+TWO_EPSILON_EB - BC%X = BC%X + DELTA*CFA%NVEC(1) - BC%Y = BC%Y + DELTA*CFA%NVEC(2) - BC%Z = BC%Z + DELTA*CFA%NVEC(3) + BC%X = BC%X + DELTA*CFA_BC%NVEC(1) + BC%Y = BC%Y + DELTA*CFA_BC%NVEC(2) + BC%Z = BC%Z + DELTA*CFA_BC%NVEC(3) ENDIF ENDIF LP%CFACE_INDEX = ICF diff --git a/Source/radi.f90 b/Source/radi.f90 index d099ee05c42..bf79e10bd25 100644 --- a/Source/radi.f90 +++ b/Source/radi.f90 @@ -3929,9 +3929,10 @@ SUBROUTINE RADIATION_FVM DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS CFA => CFACE(ICF) BR => BOUNDARY_RADIA(CFA%BR_INDEX) + BC => BOUNDARY_COORD(CFA%BC_INDEX) DO N=1,NRA DLA = (/DLX(N),DLY(N),DLZ(N)/) - DLF = DOT_PRODUCT(CFA%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle IF (DLF<0._EB) INRAD_F(ICF) = INRAD_F(ICF) - DLF*BR%BAND(IBND)%ILW(N) ENDDO ENDDO @@ -4033,8 +4034,9 @@ SUBROUTINE RADIATION_FVM CFA => CFACE(ICF) IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP1 BR => BOUNDARY_RADIA(CFA%BR_INDEX) + BC => BOUNDARY_COORD(CFA%BC_INDEX) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) - DLF = DOT_PRODUCT(CFA%NVEC,DLA) ! face normal * radiation angle + DLF = DOT_PRODUCT(BC%NVEC,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 @@ -4211,7 +4213,8 @@ SUBROUTINE RADIATION_FVM ! Loop through CFACES and assign as upwind or downwind DO IFACE=1,CF%NFACE CFA => CFACE(CF%CFACE_INDEX(IFACE)) - DLF = DOT_PRODUCT(CFA%NVEC,DLA) ! face normal * radiation angle + BC => BOUNDARY_COORD(CFA%BC_INDEX) + DLF = DOT_PRODUCT(BC%NVEC,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) @@ -4331,7 +4334,8 @@ SUBROUTINE RADIATION_FVM CFACE_LOOP2: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS CFA => CFACE(ICF) IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP2 - DLF = DOT_PRODUCT(CFA%NVEC,DLA) ! face normal * radiation angle + BC => BOUNDARY_COORD(CFA%BC_INDEX) + DLF = DOT_PRODUCT(BC%NVEC,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) @@ -4342,7 +4346,8 @@ SUBROUTINE RADIATION_FVM CFACE_LOOP3: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS CFA => CFACE(ICF) IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP3 - DLF = DOT_PRODUCT(CFA%NVEC,DLA) ! face normal * radiation angle + BC => BOUNDARY_COORD(CFA%BC_INDEX) + DLF = DOT_PRODUCT(BC%NVEC,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/soot.f90 b/Source/soot.f90 index 9b8140fcc34..8c8d22c5fba 100644 --- a/Source/soot.f90 +++ b/Source/soot.f90 @@ -257,7 +257,7 @@ SUBROUTINE CALC_DEPOSITION(WALL_INDEX,CFACE_INDEX) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) B2 => BOUNDARY_PROP2(CFA%B2_INDEX) BC => BOUNDARY_COORD(CFA%BC_INDEX) - NVEC = CFA%NVEC + NVEC = BC%NVEC TMP_G = CFA%TMP_G RHO_G = CFA%RHO_G RSUM_G= CFA%RSUM_G diff --git a/Source/type.f90 b/Source/type.f90 index 04bc0112e3c..db6d40409b2 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -193,6 +193,8 @@ MODULE TYPES REAL(EB) :: Z1 !< Lower z extent of boundary cell (m) REAL(EB) :: Z2 !< Upper z extent of boundary cell (m) + REAL(EB), DIMENSION(3) :: NVEC=0._EB !< Normal vector + END TYPE BOUNDARY_COORD_TYPE @@ -1231,7 +1233,6 @@ MODULE TYPES INTEGER :: N_INTEGERS=0 !< Number of integers to pack into restart or send/recv buffer INTEGER :: N_LOGICALS=0 !< Number of logicals to pack into restart or send/recv buffer REAL(EB) :: AREA=0._EB !< CFACE area. From CUT_FACE(CUT_FACE_IND1)%AREA(CUT_FACE_IND2) - REAL(EB) :: NVEC(3)=0._EB !< CFACE normal out unit vector. REAL(EB) :: DUNDT=0._EB REAL(EB) :: RSUM_G=0._EB !< \f$ R_0 \sum_\alpha Z_\alpha/W_\alpha \f$ in first gas phase cell REAL(EB) :: TMP_G !< Temperature (K) in adjacent gas phase cell diff --git a/Source/vege.f90 b/Source/vege.f90 index a7cf8f8c55b..a8248d212c9 100644 --- a/Source/vege.f90 +++ b/Source/vege.f90 @@ -89,7 +89,7 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_1(NM) IW = MAXLOC(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE),DIM=1) CFA => CFACE( CUT_FACE(ICF)%CFACE_INDEX(IW) ) BC => BOUNDARY_COORD(CFA%BC_INDEX) - IF (CFA%NVEC(KAXIS)>-TWO_EPSILON_EB .AND. CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + IF (BC%NVEC(KAXIS)>-TWO_EPSILON_EB .AND. CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN IF (SUM(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE))