diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 7f8f651fbbc..65c97e763c9 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -107,7 +107,7 @@ MODULE CC_SCALARS CC_CUTCELL_VELOCITY,CC_CUTFACE_VELOCITY,CC_RESTORE_UVW_UNLINKED,& CHECK_CFLVN_LINKED_CELLS,ADD_Q_DOT_CUTCELLS,CFACE_THERMAL_GASVARS,& CFACE_PREDICT_NORMAL_VELOCITY,COMPUTE_LINKED_CUTFACE_BAROCLINIC,& - COPY_CC_UNKH_TO_HS, COPY_CC_HS_TO_UNKH, COPY_UNST_DM_TO_CART, CUTFACE_VELOCITIES, & + COPY_CC_UNKH_TO_HS, COPY_CC_MUNKH_TO_UNKH, COPY_UNST_DM_TO_CART, CUTFACE_VELOCITIES, & 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, & @@ -21702,16 +21702,14 @@ SUBROUTINE NUMBER_UNKH_CUTCELLS(FLAG12,NM,IPZ,NUNKH_LC) RETURN END SUBROUTINE NUMBER_UNKH_CUTCELLS -! ------------------- COPY_CC_HS_TO_UNKH ---------------------------- +! ------------------- COPY_CC_MUNKH_TO_UNKH ---------------------------- -SUBROUTINE COPY_CC_HS_TO_UNKH(NM) - -INTEGER, INTENT(IN) :: NM +SUBROUTINE COPY_CC_MUNKH_TO_UNKH ! Local Variables: -INTEGER :: NOM,ICC,II,JJ,KK,IOR,IW,IIO,JJO,KKO +INTEGER :: NOM,ICC,IW,IIO,JJO,KKO,II,JJ,KK TYPE (OMESH_TYPE), POINTER :: OM - +TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC ! Loop over external wall cells: EXTERNAL_WALL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS @@ -21719,11 +21717,6 @@ SUBROUTINE COPY_CC_HS_TO_UNKH(NM) EWC=>EXTERNAL_WALL(IW) IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE EXTERNAL_WALL_LOOP - BC => BOUNDARY_COORD(WC%BC_INDEX) - II = BC%II - JJ = BC%JJ - KK = BC%KK - IOR = BC%IOR NOM = EWC%NOM OM => OMESH(NOM) @@ -21731,15 +21724,35 @@ SUBROUTINE COPY_CC_HS_TO_UNKH(NM) KKO=EWC%KKO_MIN JJO=EWC%JJO_MIN IIO=EWC%IIO_MIN - + BC => BOUNDARY_COORD(WC%BC_INDEX) + II = BC%II; JJ = BC%JJ; KK = BC%KK; ICC=CCVAR(II,JJ,KK,CC_IDCC) - IF (ICC > 0) THEN ! Cut-cells on this guard-cell Cartesian cell. - MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(OM%HS(IIO,JJO,KKO)) + CUT_CELL(ICC)%UNKH(1) = OM%MUNKH(IIO,JJO,KKO) ELSE - MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH) = INT(OM%HS(IIO,JJO,KKO)) + CCVAR(II,JJ,KK,CC_UNKH) = OM%MUNKH(IIO,JJO,KKO) ENDIF + ! ! Loop over all cells in mesh NOM that correspond to this boundary face (supports grid refinement): + ! DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + ! DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + ! DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + + ! ! Check if this cell in mesh NOM is a cut-cell: + ! ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC) + + ! IF (ICC > 0) THEN + ! ! Copy to cut-cell storage in mesh NOM: + ! MESHES(NOM)%CUT_CELL(ICC)%UNKH(1) = OM%MUNKH(IIO,JJO,KKO) + ! ELSE + ! ! Copy to regular cell storage in mesh NOM: + ! MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_UNKH) = OM%MUNKH(IIO,JJO,KKO) + ! ENDIF + + ! ENDDO + ! ENDDO + ! ENDDO + ENDDO EXTERNAL_WALL_LOOP ! Loop over external wall cells: @@ -21753,7 +21766,7 @@ SUBROUTINE COPY_CC_HS_TO_UNKH(NM) ENDDO EXTERNAL_WALL_LOOP2 RETURN -END SUBROUTINE COPY_CC_HS_TO_UNKH +END SUBROUTINE COPY_CC_MUNKH_TO_UNKH ! ------------------- COPY_CC_UNKH_TO_HS ---------------------------- diff --git a/Source/main.f90 b/Source/main.f90 index a577f4b5bcb..eaa3dec658c 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -1612,15 +1612,19 @@ SUBROUTINE GLOBAL_MATRIX_REASSIGN(FORCE_REASSIGN) CALL ULMAT_SOLVER_SETUP(NM) ENDDO CALL STOP_CHECK(1) - CASE (UGLMAT_FLAG,GLMAT_FLAG) - IF(ALLOCATED(ZONE_SOLVE)) CALL FINISH_GLMAT_SOLVER - CALL GLMAT_SOLVER_SETUP(1) - CALL STOP_CHECK(1) - CALL MESH_EXCHANGE(3) ! Exchange guard cell info for CCVAR(I,J,K,CGSC) -> HS. - CALL GLMAT_SOLVER_SETUP(2) - CALL MESH_EXCHANGE(3) ! Exchange guard cell info for CCVAR(I,J,K,UNKH) -> HS. - CALL GLMAT_SOLVER_SETUP(3) - CALL STOP_CHECK(1) + CASE (UGLMAT_FLAG,GLMAT_FLAG) + IF(ALLOCATED(ZONE_SOLVE)) CALL FINISH_GLMAT_SOLVER + CALL GLMAT_SOLVER_SETUP(-1) ! Initialize EWC_TYPE, copy wall types to HS + CALL STOP_CHECK(1) + CALL MESH_EXCHANGE(3) ! Exchange guard cell info for IS_WALLT -> HS + CALL GLMAT_SOLVER_SETUP(0) ! Process coarse faces, copy updated wall types to HS + CALL MESH_EXCHANGE(3) ! Re-exchange guard cell info for IS_WALLT -> HS + CALL GLMAT_SOLVER_SETUP(1) + CALL MESH_EXCHANGE(3) ! Exchange guard cell info for CCVAR(I,J,K,CGSC) -> HS + CALL GLMAT_SOLVER_SETUP(2) + CALL MESH_EXCHANGE(3) ! Exchange guard cell info for CCVAR(I,J,K,UNKH) -> HS + CALL GLMAT_SOLVER_SETUP(3) + CALL STOP_CHECK(1) END SELECT OBST_CREATED_OR_REMOVED = .FALSE. ENDIF diff --git a/Source/pres.f90 b/Source/pres.f90 index d64b1dff083..27acbb3122b 100644 --- a/Source/pres.f90 +++ b/Source/pres.f90 @@ -3016,7 +3016,7 @@ MODULE GLOBMAT_SOLVER USE COMPLEX_GEOMETRY, ONLY : CALL_FOR_GLMAT, CC_CGSC,CC_FGSC, CC_UNKH, CC_NCVARS, & NM_START,IPARM,NNZ_ROW_H,CALL_FROM_GLMAT_SETUP USE CC_SCALARS, ONLY : GET_H_CUTFACES, GET_BOUNDFACE_GEOM_INFO_H, ADD_INPLACE_NNZ_H_WHLDOM, & - COPY_CC_HS_TO_UNKH, COPY_CC_UNKH_TO_HS + COPY_CC_MUNKH_TO_UNKH, COPY_CC_UNKH_TO_HS #ifdef WITH_MKL USE MKL_CLUSTER_SPARSE_SOLVER @@ -3037,6 +3037,7 @@ MODULE GLOBMAT_SOLVER INTEGER, PARAMETER :: IS_CGSC = 1 ! Face media type: IS_GASPHASE, IS_SOLID or IS_CUTCFE. INTEGER, PARAMETER :: IS_UNKH = 2 ! H unknown number. INTEGER, PARAMETER :: IS_NCVARS = 2 ! Number of face variables in MESHES(NM)%CCVAR. +INTEGER, PARAMETER :: IS_WALLT = 100 ! Wall cell type. INTEGER, SAVE :: ILO_CELL,IHI_CELL,JLO_CELL,JHI_CELL,KLO_CELL,KHI_CELL INTEGER, SAVE :: ILO_FACE,IHI_FACE,JLO_FACE,JHI_FACE,KLO_FACE,KHI_FACE @@ -3082,6 +3083,24 @@ MODULE GLOBMAT_SOLVER CONTAINS +! --------------------------- COMPUTE_GUARD_CELL_INDEXES ----------------------- + +PURE SUBROUTINE COMPUTE_GUARD_CELL_INDEXES(IOR_IN, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) +! Compute the guard cell location in the neighboring mesh based on orientation +INTEGER, INTENT(IN) :: IOR_IN, IIO, JJO, KKO +INTEGER, INTENT(OUT) :: II_NOM, JJ_NOM, KK_NOM + +SELECT CASE(IOR_IN) +CASE( IAXIS); II_NOM = IIO + 1; JJ_NOM = JJO; KK_NOM = KKO +CASE(-IAXIS); II_NOM = IIO - 1; JJ_NOM = JJO; KK_NOM = KKO +CASE( JAXIS); II_NOM = IIO; JJ_NOM = JJO + 1; KK_NOM = KKO +CASE(-JAXIS); II_NOM = IIO; JJ_NOM = JJO - 1; KK_NOM = KKO +CASE( KAXIS); II_NOM = IIO; JJ_NOM = JJO; KK_NOM = KKO + 1 +CASE(-KAXIS); II_NOM = IIO; JJ_NOM = JJO; KK_NOM = KKO - 1 +END SELECT + +END SUBROUTINE COMPUTE_GUARD_CELL_INDEXES + ! --------------------------- GLMAT_SOLVER ------------------------------------- SUBROUTINE GLMAT_SOLVER(T,DT) @@ -3380,6 +3399,11 @@ SUBROUTINE GLMAT_SOLVER_SETUP(STAGE_FLAG) ! Local Variables: LOGICAL :: SUPPORTED_MESH=.TRUE. +LOGICAL :: FINE_INTERPOLATED_FLG, FINE_SOLID_FLG +INTEGER :: NM,IW,IIO,JJO,KKO,II_NOM,JJ_NOM,KK_NOM +TYPE(WALL_TYPE), POINTER :: WC +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC IF (FREEZE_VELOCITY) RETURN ! Fixed velocity soln. i.e. PERIODIC_TEST=102 => FREEZE_VELOCITY=.TRUE. IF (SOLID_PHASE_ONLY) RETURN @@ -3400,20 +3424,105 @@ SUBROUTINE GLMAT_SOLVER_SETUP(STAGE_FLAG) IF(PRES_ON_WHOLE_DOMAIN) N_ZONE_GLOBMAT = 0 SELECT CASE(STAGE_FLAG) -CASE(1) +CASE(-1) ! Initialization of EWC_TYPE array: + + CALL_FROM_GLMAT_SETUP = .TRUE. + + ! Factor to drop DY(J) in cylindrical coordinates. Soln assumes DTheta=1. + CYL_FCT = 0._EB; IF (CYLINDRICAL) CYL_FCT = 1._EB + + ! Check for unsupported mesh configurations: + CALL CHECK_UNSUPPORTED_MESH(SUPPORTED_MESH) + IF (.NOT.SUPPORTED_MESH) RETURN - CALL_FROM_GLMAT_SETUP = .TRUE. + ITERATE_PRESSURE = .TRUE. ! Although there is no need to do pressure iterations to drive down velocity error + ! on wall cells (i.e. the solution should give the right unique dH/dxn), leave it + ! .TRUE. to write out velocity error diagnostics. - ! Factor to drop DY(J) in cylindrical coordinates. Soln assumes DTheta=1. - CYL_FCT = 0._EB; IF (CYLINDRICAL) CYL_FCT = 1._EB + ! Copy external wall cells types to HS: + CALL COPY_CCVAR_IN_HS(IS_WALLT) - ! Check for unsupported mesh configurations: - CALL CHECK_UNSUPPORTED_MESH(SUPPORTED_MESH) - IF (.NOT.SUPPORTED_MESH) RETURN +CASE(0) + ! Copy external wall cells types from HS to OMESH%EWC_TYPE: + CALL COPY_HS_IN_CCVAR(IS_WALLT) - ITERATE_PRESSURE = .TRUE. ! Although there is no need to do pressure iterations to drive down velocity error - ! on wall cells (i.e. the solution should give the right unique dH/dxn), leave it - ! .TRUE. to write out velocity error diagnostics. + ! Here the INTERPOLATED coarse side faces note if at least one fine face is INTERPOLATED and one SOLID, set to SOLID. + MESH_LOOP_1: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + DO IW=1,N_EXTERNAL_WALL_CELLS + WC => WALL(IW) + BC => BOUNDARY_COORD(WC%BC_INDEX) + EWC => EXTERNAL_WALL(IW) + IF(WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .AND. EWC%NOM > 0 .AND. EWC%AREA_RATIO>0.9_EB) THEN + IF(ALLOCATED(OMESH(EWC%NOM)%EWC_TYPE)) THEN + ! Only coarse mesh side checked for one INTERPOLATED fine side face and one SOLID fine side face + FINE_INTERPOLATED_FLG = .FALSE. + FINE_SOLID_FLG = .FALSE. + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + ! Compute guard cell location in neighboring mesh + CALL COMPUTE_GUARD_CELL_INDEXES(BC%IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + IF(OMESH(EWC%NOM)%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) == INTERPOLATED_BOUNDARY) THEN + FINE_INTERPOLATED_FLG = .TRUE. + ELSE + FINE_SOLID_FLG = .TRUE. + ENDIF + ENDDO + ENDDO + ENDDO + IF(FINE_INTERPOLATED_FLG .AND. FINE_SOLID_FLG) WC%BOUNDARY_TYPE = SOLID_BOUNDARY + ENDIF + ENDIF + ENDDO + ENDDO MESH_LOOP_1 + + ! Copy external wall cells types to HS: + CALL COPY_CCVAR_IN_HS(IS_WALLT) + +CASE(1) + + ! Copy external wall cells types from HS to OMESH%EWC_TYPE: + CALL COPY_HS_IN_CCVAR(IS_WALLT) + + ! Here for INTERPOLATED fine side faces note if coarse side face is SOLID, if so set their type to SOLID + MESH_LOOP_2: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + DO IW=1,N_EXTERNAL_WALL_CELLS + WC => WALL(IW) + BC => BOUNDARY_COORD(WC%BC_INDEX) + EWC => EXTERNAL_WALL(IW) + IF(WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. EWC%NOM > 0) THEN + IF(ALLOCATED(OMESH(EWC%NOM)%EWC_TYPE) .AND. EWC%AREA_RATIO>0.9_EB) THEN + ! SOLID coarse face looks for INTERPOLATED fine side faces, if found sets EWC_TYPE to SOLID. + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + ! Compute guard cell location in neighboring mesh + CALL COMPUTE_GUARD_CELL_INDEXES(BC%IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + IF(OMESH(EWC%NOM)%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) == INTERPOLATED_BOUNDARY) & + OMESH(EWC%NOM)%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) = SOLID_BOUNDARY + ENDDO + ENDDO + ENDDO + ENDIF + ELSEIF(WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) THEN + IF(ALLOCATED(OMESH(EWC%NOM)%EWC_TYPE) .AND. EWC%AREA_RATIO<0.9_EB) THEN + ! Only fine mesh side checked for SOLID_BOUNDARY coarse side faces. + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + ! Compute guard cell location in neighboring mesh + CALL COMPUTE_GUARD_CELL_INDEXES(BC%IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + IF(OMESH(EWC%NOM)%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) == SOLID_BOUNDARY) & + WC%BOUNDARY_TYPE = SOLID_BOUNDARY + ENDDO + ENDDO + ENDDO + ENDIF + ENDIF + ENDDO + ENDDO MESH_LOOP_2 ! Test for CC_IBM, define CGSC and UNKH locations in CCVAR: IF (CC_IBM) THEN @@ -3492,7 +3601,6 @@ SUBROUTINE CHECK_UNSUPPORTED_MESH(SUPPORTED_MESH) USE MPI_F08 USE MESH_POINTERS USE GLOBAL_CONSTANTS, ONLY : N_MPI_PROCESSES -USE TRAN, ONLY : TRANS LOGICAL, INTENT(OUT) :: SUPPORTED_MESH @@ -3509,64 +3617,50 @@ SUBROUTINE CHECK_UNSUPPORTED_MESH(SUPPORTED_MESH) SUPPORTED_MESH = .TRUE. -! 1. Stretched grids which is untested: -GLMAT_IF_1 : IF(PRES_FLAG==GLMAT_FLAG) THEN -TRN_ME(1:2) = 0 -MESH_LOOP_TRN : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - TRN_ME(1) = TRN_ME(1) + TRANS(NM)%NOCMAX -ENDDO MESH_LOOP_TRN -TRN_ME(2)=TRN_ME(1) -IF (N_MPI_PROCESSES > 1) CALL MPI_ALLREDUCE(TRN_ME(1),TRN_ME(2),1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR) -IF (TRN_ME(2) > 0) THEN ! There is a TRNX, TRNY or TRNZ line defined for stretched grids. Not Unsupported. - IF (MY_RANK == 0) WRITE(LU_ERR,*) 'GLMAT Setup Error : Stretched grids currently unsupported.' - SUPPORTED_MESH = .FALSE. - STOP_STATUS = SETUP_STOP - RETURN -ENDIF -ENDIF GLMAT_IF_1 - IF (NMESHES == 1) RETURN -! 2. Two different cell sizes in mesh (i.e. different refinement levels): -NM = 1 -IF (MY_RANK==PROCESS(NM)) THEN - CALL POINT_TO_MESH(NM) - DX_P(IAXIS) = DX(1) - DX_P(JAXIS) = DY(1) - DX_P(KAXIS) = DZ(1) -ENDIF -IF (N_MPI_PROCESSES > 1) CALL MPI_BCAST(DX_P(1),3,MPI_DOUBLE_PRECISION,PROCESS(NM),MPI_COMM_WORLD,IERR) -! Find domain sizes to define relative epsilon: -MIN_XS(1:3) = (/ MESHES(NM)%XS, MESHES(NM)%YS, MESHES(NM)%ZS /) -MAX_XF(1:3) = (/ MESHES(NM)%XF, MESHES(NM)%YF, MESHES(NM)%ZF /) -DO NM=2,NMESHES - MIN_XS(1) = MIN(MIN_XS(1),MESHES(NM)%XS) - MIN_XS(2) = MIN(MIN_XS(2),MESHES(NM)%YS) - MIN_XS(3) = MIN(MIN_XS(3),MESHES(NM)%ZS) - MAX_XF(1) = MAX(MAX_XF(1),MESHES(NM)%XF) - MAX_XF(2) = MAX(MAX_XF(2),MESHES(NM)%YF) - MAX_XF(3) = MAX(MAX_XF(3),MESHES(NM)%ZF) -ENDDO -LX = MAX(MAX_XF(1)-MIN_XS(1),1._EB) -LY = MAX(MAX_XF(2)-MIN_XS(2),1._EB) -LZ = MAX(MAX_XF(3)-MIN_XS(3),1._EB) -TRN_ME(1:2) = 0 -MESH_LOOP_CELL : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - CALL POINT_TO_MESH(NM) - IF(ABS(DX_P(IAXIS)-DX(1)) > 10._EB*TWO_EPSILON_EB*LX) TRN_ME(1) = TRN_ME(1) + 1 - IF(ABS(DX_P(JAXIS)-DY(1)) > 10._EB*TWO_EPSILON_EB*LY) TRN_ME(1) = TRN_ME(1) + 1 - IF(ABS(DX_P(KAXIS)-DZ(1)) > 10._EB*TWO_EPSILON_EB*LZ) TRN_ME(1) = TRN_ME(1) + 1 -ENDDO MESH_LOOP_CELL -TRN_ME(2)=TRN_ME(1) -IF (N_MPI_PROCESSES > 1) CALL MPI_ALLREDUCE(TRN_ME(1),TRN_ME(2),1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR) -IF (TRN_ME(2) > 0) THEN ! Meshes at different refinement levels. Not Unsupported. - IF (MY_RANK == 0) WRITE(LU_ERR,*) 'GLMAT Setup Error : Meshes at different refinement levels unsupported.' - SUPPORTED_MESH = .FALSE. - STOP_STATUS = SETUP_STOP - RETURN +! 1. For now unsupported: CC_IBM and two different cell sizes in mesh (i.e. different refinement levels): +IF(CC_IBM) THEN + NM = 1 + IF (MY_RANK==PROCESS(NM)) THEN + CALL POINT_TO_MESH(NM) + DX_P(IAXIS) = DX(1) + DX_P(JAXIS) = DY(1) + DX_P(KAXIS) = DZ(1) + ENDIF + IF (N_MPI_PROCESSES > 1) CALL MPI_BCAST(DX_P(1),3,MPI_DOUBLE_PRECISION,PROCESS(NM),MPI_COMM_WORLD,IERR) + ! Find domain sizes to define relative epsilon: + MIN_XS(1:3) = (/ MESHES(NM)%XS, MESHES(NM)%YS, MESHES(NM)%ZS /) + MAX_XF(1:3) = (/ MESHES(NM)%XF, MESHES(NM)%YF, MESHES(NM)%ZF /) + DO NM=2,NMESHES + MIN_XS(1) = MIN(MIN_XS(1),MESHES(NM)%XS) + MIN_XS(2) = MIN(MIN_XS(2),MESHES(NM)%YS) + MIN_XS(3) = MIN(MIN_XS(3),MESHES(NM)%ZS) + MAX_XF(1) = MAX(MAX_XF(1),MESHES(NM)%XF) + MAX_XF(2) = MAX(MAX_XF(2),MESHES(NM)%YF) + MAX_XF(3) = MAX(MAX_XF(3),MESHES(NM)%ZF) + ENDDO + LX = MAX(MAX_XF(1)-MIN_XS(1),1._EB) + LY = MAX(MAX_XF(2)-MIN_XS(2),1._EB) + LZ = MAX(MAX_XF(3)-MIN_XS(3),1._EB) + TRN_ME(1:2) = 0 + MESH_LOOP_CELL : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + IF(ABS(DX_P(IAXIS)-DX(1)) > 10._EB*TWO_EPSILON_EB*LX) TRN_ME(1) = TRN_ME(1) + 1 + IF(ABS(DX_P(JAXIS)-DY(1)) > 10._EB*TWO_EPSILON_EB*LY) TRN_ME(1) = TRN_ME(1) + 1 + IF(ABS(DX_P(KAXIS)-DZ(1)) > 10._EB*TWO_EPSILON_EB*LZ) TRN_ME(1) = TRN_ME(1) + 1 + ENDDO MESH_LOOP_CELL + TRN_ME(2)=TRN_ME(1) + IF (N_MPI_PROCESSES > 1) CALL MPI_ALLREDUCE(TRN_ME(1),TRN_ME(2),1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR) + IF (TRN_ME(2) > 0) THEN ! Meshes at different refinement levels. Not Unsupported. + IF (MY_RANK == 0) WRITE(LU_ERR,*) 'GLMAT Setup Error : Meshes at different refinement levels unsupported.' + SUPPORTED_MESH = .FALSE. + STOP_STATUS = SETUP_STOP + RETURN + ENDIF ENDIF -! 3. Two (or more) disjoint domains, where at least one has all Neumann BCs and one has some Dirichlet bcs. +! 2. Two (or more) disjoint domains, where at least one has all Neumann BCs and one has some Dirichlet bcs. ! This is a topological problem that would require different Matrix types (i.e. one positive definite and one ! indefinite), which would require separate solutions. ! A possible approach to look at is to solve the whole system as indefinite, and then substract a constant in @@ -3692,12 +3786,13 @@ SUBROUTINE COPY_H_OMESH_TO_MESH USE COMP_FUNCTIONS, ONLY: CURRENT_TIME USE CC_SCALARS, ONLY: GET_H_GUARD_CUTCELL ! Local Variables: -INTEGER :: NM,NOM,II,JJ,KK,IOR,IW,IIO,JJO,KKO +INTEGER :: NM,NOM,II,JJ,KK,IOR,IW,IIO,JJO,KKO,IIG,JJG,KKG,II_NOM,JJ_NOM,KK_NOM TYPE (OMESH_TYPE), POINTER :: OM TYPE (WALL_TYPE), POINTER :: WC TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC LOGICAL :: FLG +REAL(EB) :: H_EXTERNAL_MEAN, DX_INT, DX_EXT, WEIGHT_INT, WEIGHT_EXT, NCELLS_EXT, H_FACE IF (CC_IBM) CALL_FOR_GLMAT = .FALSE. @@ -3729,34 +3824,65 @@ SUBROUTINE COPY_H_OMESH_TO_MESH IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE EXTERNAL_WALL_LOOP_1 ENDIF - II = BC%II - JJ = BC%JJ - KK = BC%KK - IOR = BC%IOR - NOM = EWC%NOM + II = BC%II; JJ = BC%JJ; KK = BC%KK + IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG + IOR = BC%IOR; NOM = EWC%NOM; IF(NOM < 1) CYCLE + ! Here if NOM==0 means it is an OBST laying on an external boundary -> CYCLE - IF(NOM < 1) CYCLE OM => OMESH(NOM) - ! This assumes all meshes at the same level of refinement: - KKO=EWC%KKO_MIN - JJO=EWC%JJO_MIN - IIO=EWC%IIO_MIN + IF (CC_IBM) THEN + ! This assumes all meshes at the same level of refinement: For now. + KKO=EWC%KKO_MIN + JJO=EWC%JJO_MIN + IIO=EWC%IIO_MIN + H(II,JJ,KK) = OM%H(IIO,JJO,KKO) + CYCLE EXTERNAL_WALL_LOOP_1 + ENDIF + ! GRID REFINEMENT: Compute mean H accounting for boundary types + H_EXTERNAL_MEAN = 0._EB + NCELLS_EXT = 0._EB + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. OM%MUNKH(IIO,JJO,KKO) <= 0) CYCLE + + ! Compute guard cell location in neighboring mesh based on IOR + CALL COMPUTE_GUARD_CELL_INDEXES(IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + + ! Check boundary type and use appropriate value + IF (OM%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) == INTERPOLATED_BOUNDARY) THEN + ! Use actual external cell value (gradient continuity) + H_EXTERNAL_MEAN = H_EXTERNAL_MEAN + OM%H(IIO,JJO,KKO) + ELSE + ! Use internal cell value (zero gradient assumption) + H_EXTERNAL_MEAN = H_EXTERNAL_MEAN + H(IIG,JJG,KKG) + ENDIF + NCELLS_EXT = NCELLS_EXT + 1._EB + ENDDO + ENDDO + ENDDO + H_EXTERNAL_MEAN = H_EXTERNAL_MEAN / NCELLS_EXT + + ! Get cell sizes for distance-weighted interpolation SELECT CASE(IOR) - CASE( IAXIS) - H(II,JJ,KK) = OM%H(IIO,JJO,KKO) - CASE(-IAXIS) - H(II,JJ,KK) = OM%H(IIO,JJO,KKO) - CASE( JAXIS) - H(II,JJ,KK) = OM%H(IIO,JJO,KKO) - CASE(-JAXIS) - H(II,JJ,KK) = OM%H(IIO,JJO,KKO) - CASE( KAXIS) - H(II,JJ,KK) = OM%H(IIO,JJO,KKO) - CASE(-KAXIS) - H(II,JJ,KK) = OM%H(IIO,JJO,KKO) + CASE( IAXIS,-IAXIS) + DX_INT = DX(IIG) + DX_EXT = MESHES(NOM)%DX(EWC%IIO_MIN) + CASE( JAXIS,-JAXIS) + DX_INT = DY(JJG) + DX_EXT = MESHES(NOM)%DY(EWC%JJO_MIN) + CASE( KAXIS,-KAXIS) + DX_INT = DZ(KKG) + DX_EXT = MESHES(NOM)%DZ(EWC%KKO_MIN) END SELECT + + ! Distance-weighted interpolation at face, then extrapolation to guard cell center + WEIGHT_EXT = DX_INT / (DX_INT + DX_EXT) + WEIGHT_INT = DX_EXT / (DX_INT + DX_EXT) + H_FACE = WEIGHT_INT * H(IIG,JJG,KKG) + WEIGHT_EXT * H_EXTERNAL_MEAN + H(II,JJ,KK) = H(IIG,JJG,KKG) + 2.0_EB * (H_FACE - H(IIG,JJG,KKG)) ENDDO EXTERNAL_WALL_LOOP_1 @@ -3788,34 +3914,64 @@ SUBROUTINE COPY_H_OMESH_TO_MESH IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE EXTERNAL_WALL_LOOP_2 ENDIF - II = BC%II - JJ = BC%JJ - KK = BC%KK - IOR = BC%IOR - NOM = EWC%NOM + II = BC%II; JJ = BC%JJ; KK = BC%KK + IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG + IOR = BC%IOR; NOM = EWC%NOM; IF(NOM < 1) CYCLE ! Here if NOM==0 means it is an OBST laying on an external boundary -> CYCLE - IF(NOM < 1) CYCLE OM => OMESH(NOM) - ! This assumes all meshes at the same level of refinement: - KKO=EWC%KKO_MIN - JJO=EWC%JJO_MIN - IIO=EWC%IIO_MIN + IF (CC_IBM) THEN + ! This assumes all meshes at the same level of refinement: For now. + KKO=EWC%KKO_MIN + JJO=EWC%JJO_MIN + IIO=EWC%IIO_MIN + HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) + CYCLE EXTERNAL_WALL_LOOP_2 + ENDIF + ! GRID REFINEMENT: Compute mean HS accounting for boundary types + H_EXTERNAL_MEAN = 0._EB + NCELLS_EXT = 0._EB + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. OM%MUNKH(IIO,JJO,KKO) <= 0) CYCLE + + ! Compute guard cell location in neighboring mesh based on IOR + CALL COMPUTE_GUARD_CELL_INDEXES(IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + + ! Check boundary type and use appropriate value + IF (OM%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) == INTERPOLATED_BOUNDARY) THEN + ! Use actual external cell value (gradient continuity) + H_EXTERNAL_MEAN = H_EXTERNAL_MEAN + OM%HS(IIO,JJO,KKO) + ELSE + ! Use internal cell value (zero gradient assumption) + H_EXTERNAL_MEAN = H_EXTERNAL_MEAN + HS(IIG,JJG,KKG) + ENDIF + NCELLS_EXT = NCELLS_EXT + 1._EB + ENDDO + ENDDO + ENDDO + H_EXTERNAL_MEAN = H_EXTERNAL_MEAN / NCELLS_EXT + + ! Get cell sizes for distance-weighted interpolation SELECT CASE(IOR) - CASE( IAXIS) - HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) - CASE(-IAXIS) - HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) - CASE( JAXIS) - HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) - CASE(-JAXIS) - HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) - CASE( KAXIS) - HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) - CASE(-KAXIS) - HS(II,JJ,KK) = OM%HS(IIO,JJO,KKO) + CASE( IAXIS,-IAXIS) + DX_INT = DX(IIG) + DX_EXT = MESHES(NOM)%DX(EWC%IIO_MIN) + CASE( JAXIS,-JAXIS) + DX_INT = DY(JJG) + DX_EXT = MESHES(NOM)%DY(EWC%JJO_MIN) + CASE( KAXIS,-KAXIS) + DX_INT = DZ(KKG) + DX_EXT = MESHES(NOM)%DZ(EWC%KKO_MIN) END SELECT + + ! Distance-weighted interpolation at face, then extrapolation to guard cell center + WEIGHT_EXT = DX_INT / (DX_INT + DX_EXT) + WEIGHT_INT = DX_EXT / (DX_INT + DX_EXT) + H_FACE = WEIGHT_INT * HS(IIG,JJG,KKG) + WEIGHT_EXT * H_EXTERNAL_MEAN + HS(II,JJ,KK) = HS(IIG,JJG,KKG) + 2.0_EB * (H_FACE - HS(IIG,JJG,KKG)) ENDDO EXTERNAL_WALL_LOOP_2 @@ -3837,70 +3993,129 @@ SUBROUTINE COPY_HS_IN_CCVAR(VAR_CC) INTEGER, INTENT(IN) :: VAR_CC ! Local Variables: -INTEGER :: NM,NOM,II,JJ,KK,IOR,IW,IIO,JJO,KKO -TYPE (OMESH_TYPE), POINTER :: OM -TYPE (WALL_TYPE), POINTER :: WC -TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC -LOGICAL :: FLG - -MESH_LOOP : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX +INTEGER :: NM,NOM +LOGICAL, PARAMETER :: WRITE_EWC_TYPE = .FALSE. +! GRID REFINEMENT: Allocate OMESH arrays for GSCH, MUNKH, and EWC_TYPE if needed +DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL POINT_TO_MESH(NM) + DO NOM=1,NMESHES + IF (.NOT.ALLOCATED(OMESH(NOM)%HS)) CYCLE + + ! Allocate GSCH if not already allocated (same bounds as HS): + IF (VAR_CC==CGSC .AND. .NOT.ALLOCATED(OMESH(NOM)%GSCH)) THEN + ALLOCATE(OMESH(NOM)%GSCH(LBOUND(OMESH(NOM)%HS,1):UBOUND(OMESH(NOM)%HS,1), & + LBOUND(OMESH(NOM)%HS,2):UBOUND(OMESH(NOM)%HS,2), & + LBOUND(OMESH(NOM)%HS,3):UBOUND(OMESH(NOM)%HS,3))) + ENDIF + + ! Allocate MUNKH if not already allocated (same bounds as HS): + IF (VAR_CC==UNKH .AND. .NOT.ALLOCATED(OMESH(NOM)%MUNKH)) THEN + ALLOCATE(OMESH(NOM)%MUNKH(LBOUND(OMESH(NOM)%HS,1):UBOUND(OMESH(NOM)%HS,1), & + LBOUND(OMESH(NOM)%HS,2):UBOUND(OMESH(NOM)%HS,2), & + LBOUND(OMESH(NOM)%HS,3):UBOUND(OMESH(NOM)%HS,3))) + ENDIF + + ! Allocate EWC_TYPE if not already allocated (same bounds as HS): + IF (VAR_CC==IS_WALLT .AND. .NOT.ALLOCATED(OMESH(NOM)%EWC_TYPE)) THEN + ALLOCATE(OMESH(NOM)%EWC_TYPE(LBOUND(OMESH(NOM)%HS,1):UBOUND(OMESH(NOM)%HS,1), & + LBOUND(OMESH(NOM)%HS,2):UBOUND(OMESH(NOM)%HS,2), & + LBOUND(OMESH(NOM)%HS,3):UBOUND(OMESH(NOM)%HS,3))) + ENDIF + + ! Copy HS data to appropriate OMESH array: + IF (VAR_CC==CGSC .AND. ALLOCATED(OMESH(NOM)%GSCH)) THEN + OMESH(NOM)%GSCH(:,:,:) = INT(OMESH(NOM)%HS(:,:,:)) + ELSEIF (VAR_CC==UNKH .AND. ALLOCATED(OMESH(NOM)%MUNKH)) THEN + OMESH(NOM)%MUNKH(:,:,:) = INT(OMESH(NOM)%HS(:,:,:)) + ELSEIF (VAR_CC==IS_WALLT .AND. ALLOCATED(OMESH(NOM)%EWC_TYPE)) THEN + OMESH(NOM)%EWC_TYPE(:,:,:) = INT(OMESH(NOM)%HS(:,:,:)) + ENDIF + ENDDO +ENDDO - ! - IF (CC_IBM .AND. VAR_CC==UNKH) THEN - - CALL COPY_CC_HS_TO_UNKH(NM) - - ELSE - - ! Loop over external wall cells: - EXTERNAL_WALL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS +! GRID REFINEMENT: Handle CC_IBM case if needed +IF (CC_IBM .AND. VAR_CC==UNKH) THEN + DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + CALL COPY_CC_MUNKH_TO_UNKH + ENDDO +ENDIF - WC=>WALL(IW) - BC=>BOUNDARY_COORD(WC%BC_INDEX) - EWC=>EXTERNAL_WALL(IW) - IF (PRES_ON_WHOLE_DOMAIN) THEN - ! Matrix connectivities kept in cases of BOUNDARY_TYPE=INTERPOLATED_BOUNDARY, NULL_BOUNDARY, or - ! SOLID_BOUNDARY where there is an OMESH, case of OBSTS right in the boundary of meshes. - FLG = WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==NULL_BOUNDARY - FLG = FLG .OR. (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. EWC%NOM > 0) - IF (.NOT.FLG) CYCLE EXTERNAL_WALL_LOOP - ELSE - ! Case of solving for H only on the gas phase, connectivity kept only for INTERPOLATED_BOUNDARY. - IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE EXTERNAL_WALL_LOOP - ENDIF +! Zero out HS after using it +IF (VAR_CC == UNKH) THEN + DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + HS(:,:,:) = 0._EB + ENDDO +ENDIF - II = BC%II - JJ = BC%JJ - KK = BC%KK - IOR = BC%IOR - NOM = EWC%NOM - ! Here if NOM==0 means it is an OBST laying on an external boundary -> CYCLE - IF(NOM < 1) CYCLE EXTERNAL_WALL_LOOP - OM => OMESH(NOM) +! DIAGNOSTIC: Write out OMESH%EWC_TYPE for external wall cells +IF (WRITE_EWC_TYPE .AND. VAR_CC == IS_WALLT) CALL WRITE_EWC_TYPE_DIAGNOSTIC - ! This assumes all meshes at the same level of refinement: - KKO=EWC%KKO_MIN - JJO=EWC%JJO_MIN - IIO=EWC%IIO_MIN +RETURN +END SUBROUTINE COPY_HS_IN_CCVAR - CCVAR(II,JJ,KK,VAR_CC) = INT(OM%HS(IIO,JJO,KKO)) +! --------------------------- WRITE_EWC_TYPE_DIAGNOSTIC ----------------------------- - ENDDO EXTERNAL_WALL_LOOP +SUBROUTINE WRITE_EWC_TYPE_DIAGNOSTIC - ENDIF +USE MESH_POINTERS - IF (VAR_CC == UNKH) THEN - HS(:,:,:) = 0._EB - ENDIF +! Local Variables: +INTEGER :: NM,NOM,IW,IIG,JJG,KKG,IOR,IIO,JJO,KKO,BNDRY_TYPE,II_NOM,JJ_NOM,KK_NOM +TYPE(WALL_TYPE), POINTER :: WC +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC -ENDDO MESH_LOOP +DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + + EXTERNAL_WALL_LOOP_DIAG: DO IW=1,N_EXTERNAL_WALL_CELLS + WC => WALL(IW) + BC => BOUNDARY_COORD(WC%BC_INDEX) + EWC => EXTERNAL_WALL(IW) + + ! Skip if no neighboring mesh + NOM = EWC%NOM; IF (NOM < 1) CYCLE + + ! Get internal cell indices + IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG + IOR = BC%IOR + + ! Check if OMESH%EWC_TYPE is allocated + IF (.NOT.ALLOCATED(OMESH(NOM)%EWC_TYPE)) THEN + WRITE(0,*) 'WARNING: MESH',NM,'IW=',IW,'OMESH(',NOM,')%EWC_TYPE NOT ALLOCATED' + CYCLE + ENDIF + + ! Loop over all external cells in neighboring mesh + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + + ! Compute guard cell location in neighboring mesh based on IOR + CALL COMPUTE_GUARD_CELL_INDEXES(IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + + BNDRY_TYPE = OMESH(NOM)%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) + + ! Only print if boundary type is defined (not IS_UNDEFINED) + IF (BNDRY_TYPE /= IS_UNDEFINED) THEN + WRITE(0,'(A,I3,A,I5,A,3I4,A,I2,A,I3,A,3I4,A,3I4,A,I3)') & + 'MESH',NM,' IW=',IW,' INT(',IIG,JJG,KKG,') IOR=',IOR, & + ' -> MESH',NOM,' EXT(',IIO,JJO,KKO,') GUARD(',II_NOM,JJ_NOM,KK_NOM,') TYPE=',BNDRY_TYPE + ENDIF + + ENDDO + ENDDO + ENDDO + + ENDDO EXTERNAL_WALL_LOOP_DIAG + +ENDDO RETURN -END SUBROUTINE COPY_HS_IN_CCVAR - +END SUBROUTINE WRITE_EWC_TYPE_DIAGNOSTIC ! ------------------------------- COPY_CCVAR_IN_HS ---------------------------------- @@ -3910,14 +4125,29 @@ SUBROUTINE COPY_CCVAR_IN_HS(VAR_CC) INTEGER, INTENT(IN) :: VAR_CC ! Local Variables: -INTEGER :: NM +INTEGER :: NM,IW,II,JJ,KK +TYPE(WALL_TYPE), POINTER :: WC +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL POINT_TO_MESH(NM) - HS(0:IBP1,0:JBP1,0:KBP1) = REAL(CCVAR(0:IBP1,0:JBP1,0:KBP1,VAR_CC),EB) - - ! Now cut-cells add their single Cartesian UNKH value in HS: - IF(CC_IBM .AND. VAR_CC==UNKH) CALL COPY_CC_UNKH_TO_HS(NM) + + ! Special case for IS_WALLT: populate HS from WALL array instead of CCVAR + IF (VAR_CC == IS_WALLT) THEN + HS(:,:,:) = REAL(IS_UNDEFINED,EB) ! Initialize + DO IW=1,N_EXTERNAL_WALL_CELLS + WC => WALL(IW) + BC => BOUNDARY_COORD(WC%BC_INDEX) + II = BC%II; JJ = BC%JJ; KK = BC%KK + HS(II,JJ,KK) = REAL(WC%BOUNDARY_TYPE,EB) + ENDDO + ELSE + ! Standard case: copy from CCVAR to HS + HS(0:IBP1,0:JBP1,0:KBP1) = REAL(CCVAR(0:IBP1,0:JBP1,0:KBP1,VAR_CC),EB) + + ! Now cut-cells add their single Cartesian UNKH value in HS: + IF(CC_IBM .AND. VAR_CC==UNKH) CALL COPY_CC_UNKH_TO_HS(NM) + ENDIF ENDDO @@ -4304,7 +4534,7 @@ SUBROUTINE GET_BCS_H_MATRIX USE CC_SCALARS, ONLY : GET_CC_UNKH, GET_CFACE_OPEN_BC_COEF ! Local Variables: -INTEGER :: NM,NM1,JLOC,JCOL,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND),IERR,IIG,JJG,KKG,II,JJ,KK,IW,ILH,JLH,KLH,IRC +INTEGER :: NM,NM1,JLOC,JCOL,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND),IERR,IIG,JJG,KKG,IW,ILH,JLH,KLH,IRC REAL(EB):: AF,IDX,BIJ TYPE(WALL_TYPE), POINTER :: WC TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC @@ -4328,7 +4558,6 @@ SUBROUTINE GET_BCS_H_MATRIX IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG IF((.NOT.CC_IBM .AND. CCVAR(IIG,JJG,KKG,UNKH)<1) & .OR. ZONE_SOLVE(PRESSURE_ZONE(IIG,JJG,KKG))%CONNECTED_ZONE_PARENT/=IPZ) CYCLE - II = BC%II; JJ = BC%JJ; KK = BC%KK ! Unknowns on related cells: IF(CCVAR(IIG,JJG,KKG,CGSC)==IS_GASPHASE .OR. PRES_ON_WHOLE_DOMAIN) THEN IND(LOW_IND) = CCVAR(IIG,JJG,KKG,UNKH) ! internal cell. @@ -4404,13 +4633,18 @@ SUBROUTINE GET_H_MATRIX REAL(EB), POINTER, DIMENSION(:) :: DX1,DX2,DX3 INTEGER :: I,J,K,I1,I2,I3,IIP,JJP,KKP,IIM,JJM,KKM INTEGER :: ILOC,JLOC,IROW,JCOL,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND) -REAL(EB):: AF,IDX,BIJ,KFACE(1:2,1:2) +REAL(EB):: AF,IDX,BIJ,KFACE(1:2,1:2) ! KFACE still used for regular internal faces TYPE(CC_REGFACE_TYPE), POINTER, DIMENSION(:) :: RGF TYPE(WALL_TYPE), POINTER :: WC TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC INTEGER :: IIG,JJG,KKG,II,JJ,KK,IOR,LOCROW,IW -INTEGER :: WC_JD(1:2,1:2) +INTEGER :: IIO,JJO,KKO,NOM ! Grid refinement: loop indices and neighboring mesh number +INTEGER :: II_NOM,JJ_NOM,KK_NOM ! Grid refinement: guard cell indices in neighboring mesh +INTEGER :: IUNK_INT,IUNK_EXT,IROW_INT ! Grid refinement: unknown numbers for internal and external cells +REAL(EB):: AF_INT,AF_EXT,DX_INT,DX_EXT ! Grid refinement: areas and distances for both cells +TYPE(MESH_TYPE), POINTER :: M2 ! Grid refinement: pointer to neighboring mesh +TYPE(OMESH_TYPE), POINTER :: OM ! Grid refinement: pointer to OMESH data LOGICAL :: FLG IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT @@ -4505,61 +4739,153 @@ SUBROUTINE GET_H_MATRIX ! Next, Wall faces of type INTERPOLATED_BOUNDARY or PERIODIC_BOUNDARY: ! Here We have to do something about WALL cells that are also cut-faces, who wins? Make cut-faces take precedence. + ! + ! MODIFICATION FOR GRID REFINEMENT: + ! Loop over ALL cells in neighboring mesh and compute flux coupling coefficients for each. + ! Column indices are searched on-the-fly (no pre-storage needed). + ! LOCROW = LOW_IND WALL_LOOP_1 : DO IW=1,N_EXTERNAL_WALL_CELLS WC => WALL(IW) BC => BOUNDARY_COORD(WC%BC_INDEX) EWC=>EXTERNAL_WALL(IW) + FLG = WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY .OR. WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY IF(PRES_ON_WHOLE_DOMAIN) & FLG = FLG .OR. WC%BOUNDARY_TYPE==NULL_BOUNDARY .OR. (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. EWC%NOM > 0) IF ( .NOT.FLG .OR. EWC%NOM<1) CYCLE ! Here if NOM==0 means it is an OBST laying on an external boundary -> CYCLE - IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG + + IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG; II = BC%II; JJ = BC%JJ; KK = BC%KK; IF(ZONE_SOLVE(PRESSURE_ZONE(IIG,JJG,KKG))%CONNECTED_ZONE_PARENT/=IPZ) CYCLE - - WC_JD(1,1) = WC%JD11_INDEX - WC_JD(1,2) = WC%JD12_INDEX - WC_JD(2,1) = WC%JD21_INDEX - WC_JD(2,2) = WC%JD22_INDEX - - II = BC%II; JJ = BC%JJ; KK = BC%KK; IOR = BC%IOR - ! Check if CC_IBM -> If either IIG,JJG,KKG or II,JJ,KK cell is type IS_CUTCFE or IS_SOLID cycle: + + IOR = BC%IOR + + ! Check if CC_IBM -> If IIG,JJG,KKG cell is type IS_CUTCFE or IS_SOLID cycle: IF ( .NOT.PRES_ON_WHOLE_DOMAIN .AND. CC_IBM ) THEN - IF(CCVAR(II ,JJ ,KK ,CC_CGSC) /= IS_GASPHASE) CYCLE IF(CCVAR(IIG,JJG,KKG,CC_CGSC) /= IS_GASPHASE) CYCLE ENDIF + + ! IUNK_INT and IROW_INT: + IUNK_INT = CCVAR(IIG,JJG,KKG,UNKH) ! internal. + IROW_INT = IUNK_INT - ZSL%UNKH_IND(NM1) - ! Unknowns on related cells: - IND(LOW_IND) = CCVAR(IIG,JJG,KKG,UNKH) ! internal. - IND(HIGH_IND) = CCVAR(II,JJ,KK,UNKH) ! guard-cell. - IND_LOC(LOW_IND) = IND(LOW_IND) - ZSL%UNKH_IND(NM1) ! All row indexes must refer to ind_loc. - IND_LOC(HIGH_IND)= IND(HIGH_IND)- ZSL%UNKH_IND(NM1) + ! Get neighboring mesh information: + NOM = EWC%NOM + M2 => MESHES(NOM) + OM => OMESH(NOM) + + ! Distance from internal cell center to boundary face (half cell width): SELECT CASE(IOR) CASE( IAXIS) - AF = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*R(IIG-1)) * DZ(KKG); IDX= 1._EB/DXN(IIG-1) + DX_INT = 0.5_EB * DXN(IIG-1) CASE(-IAXIS) - AF = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*R(IIG )) * DZ(KKG); IDX= 1._EB/DXN(IIG) + DX_INT = 0.5_EB * DXN(IIG) CASE( JAXIS) - AF = DX(IIG)*DZ(KKG); IDX= 1._EB/DYN(JJG-1) + DX_INT = 0.5_EB * DYN(JJG-1) CASE(-JAXIS) - AF = DX(IIG)*DZ(KKG); IDX= 1._EB/DYN(JJG) + DX_INT = 0.5_EB * DYN(JJG) CASE( KAXIS) - AF = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*RC(IIG ))* DX(IIG); IDX= 1._EB/DZN(KKG-1) + DX_INT = 0.5_EB * DZN(KKG-1) CASE(-KAXIS) - AF = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*RC(IIG ))* DX(IIG); IDX= 1._EB/DZN(KKG) + DX_INT = 0.5_EB * DZN(KKG) END SELECT - ! Now add to Adiff corresponding coeff: - BIJ = IDX*AF - ! Cols ind(1) ind(2) - KFACE(1,1)= BIJ; KFACE(1,2)=-BIJ ! Row ind(1) - KFACE(2,1)=-BIJ; KFACE(2,2)= BIJ ! Row ind(2) - ILOC = LOCROW ! Local row number in Kface, only for cell IIG,JJG,KKG. - DO JLOC = LOW_IND,HIGH_IND ! Local col number in Kface, JD - IROW = IND_LOC(ILOC) ! Unknown number. - JCOL = WC_JD(ILOC,JLOC) - ZSL%ROW_H(IROW)%D(JCOL) = ZSL%ROW_H(IROW)%D(JCOL) + KFACE(ILOC,JLOC) - ENDDO + + ! Loop over ALL cells in the neighboring mesh that share this boundary face: + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + + ! Get unknown number for external cell in neighboring mesh: + IUNK_EXT = OM%MUNKH(IIO,JJO,KKO) + + ! For UGLMAT, skip if external cell is solid: + IF ( .NOT.PRES_ON_WHOLE_DOMAIN .AND.IUNK_EXT <= 0) CYCLE + + ! GRID REFINEMENT: Compute guard cell location in neighboring mesh based on IOR + CALL COMPUTE_GUARD_CELL_INDEXES(IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + + ! Check if external cell has INTERPOLATED_BOUNDARY type + ! If not INTERPOLATED_BOUNDARY, skip (zero coefficient, no flux coupling) + IF (ALLOCATED(OM%EWC_TYPE)) THEN + IF (OM%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) /= INTERPOLATED_BOUNDARY) CYCLE + ENDIF + + ! Area of INTERNAL cell face (current mesh): + SELECT CASE(IOR) + CASE( IAXIS, -IAXIS) + IF (CYLINDRICAL) THEN + AF_INT = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*R(IIG)) * DZ(KKG) + ELSE + AF_INT = DY(JJG) * DZ(KKG) + ENDIF + CASE( JAXIS, -JAXIS) + AF_INT = DX(IIG) * DZ(KKG) + CASE( KAXIS, -KAXIS) + IF (CYLINDRICAL) THEN + AF_INT = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*RC(IIG)) * DX(IIG) + ELSE + AF_INT = DX(IIG) * DY(JJG) + ENDIF + END SELECT + + ! Area of EXTERNAL cell face (other mesh): + SELECT CASE(IOR) + CASE( IAXIS, -IAXIS) + IF (CYLINDRICAL) THEN + AF_EXT = ((1._EB-CYL_FCT)*M2%DY(JJO) + CYL_FCT*M2%R(IIO)) * M2%DZ(KKO) + ELSE + AF_EXT = M2%DY(JJO) * M2%DZ(KKO) + ENDIF + CASE( JAXIS, -JAXIS) + AF_EXT = M2%DX(IIO) * M2%DZ(KKO) + CASE( KAXIS, -KAXIS) + IF (CYLINDRICAL) THEN + AF_EXT = ((1._EB-CYL_FCT)*M2%DY(JJO) + CYL_FCT*M2%RC(IIO)) * M2%DX(IIO) + ELSE + AF_EXT = M2%DX(IIO) * M2%DY(JJO) + ENDIF + END SELECT + + ! Use the MINIMUM area (actual contact/intersection area): + AF = MIN(AF_INT, AF_EXT) + + ! Distance from boundary face to external cell center (half cell width): + SELECT CASE(IOR) + CASE( IAXIS, -IAXIS) + DX_EXT = 0.5_EB * M2%DXN(IIO) + CASE( JAXIS, -JAXIS) + DX_EXT = 0.5_EB * M2%DYN(JJO) + CASE( KAXIS, -KAXIS) + DX_EXT = 0.5_EB * M2%DZN(KKO) + END SELECT + + ! Inverse of total distance between cell centers: + IDX = 1._EB / (DX_INT + DX_EXT) + + ! Flux coupling coefficient: + BIJ = IDX * AF + + ! (1) Add to DIAGONAL of internal cell row: + DO JLOC = 1, ZSL%ROW_H(IROW_INT)%NNZ + IF (IUNK_INT == ZSL%ROW_H(IROW_INT)%JD(JLOC)) THEN + ZSL%ROW_H(IROW_INT)%D(JLOC) = ZSL%ROW_H(IROW_INT)%D(JLOC) + BIJ + EXIT + ENDIF + ENDDO + + ! (2) Add to OFF-DIAGONAL coupling to external cell: + DO JLOC = 1, ZSL%ROW_H(IROW_INT)%NNZ + IF (IUNK_EXT == ZSL%ROW_H(IROW_INT)%JD(JLOC)) THEN + ZSL%ROW_H(IROW_INT)%D(JLOC) = ZSL%ROW_H(IROW_INT)%D(JLOC) - BIJ + EXIT + ENDIF + ENDDO + + ENDDO ! IIO loop + ENDDO ! JJO loop + ENDDO ! KKO loop + ENDDO WALL_LOOP_1 ! Contribution to Laplacian matrix from RC and cut-faces: @@ -4583,13 +4909,14 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM INTEGER :: NM INTEGER :: X1AXIS,IFACE,I,I1,J,K,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND) INTEGER :: LOCROW,IIND,NII,ILOC -INTEGER :: NREG,IIM,JJM,KKM,IIP,JJP,KKP,LOW_FACE,HIGH_FACE,IW,II,JJ,KK,IIG,JJG,KKG,ICF +INTEGER :: NREG,IIM,JJM,KKM,IIP,JJP,KKP,LOW_FACE,HIGH_FACE,IW,IIG,JJG,KKG,II,JJ,KK,ICF TYPE(CC_REGFACE_TYPE), POINTER, DIMENSION(:) :: RGF TYPE(CC_RCFACE_TYPE), POINTER :: RCF TYPE(WALL_TYPE), POINTER :: WC TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC -INTEGER :: WC_JD(1:2,1:2) +INTEGER :: IIO,JJO,KKO,NOM ! Grid refinement: loop indices and neighboring mesh number +INTEGER :: II_NOM,JJ_NOM,KK_NOM ! Grid refinement: guard cell indices in neighboring mesh LOGICAL :: FLG ! Write number of pressure unknowns to output: @@ -4695,6 +5022,11 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM ! Next, Wall faces of type INTERPOLATED_BOUNDARY or PERIODIC_BOUNDARY: ! Here We have to do something about WALL cells that are also cut-faces, who wins? Make cut-faces take precedence. + ! + ! MODIFICATION FOR GRID REFINEMENT: + ! Loop over ALL cells in the neighboring mesh (EWC%IIO_MIN:IIO_MAX, etc.) to handle cases where + ! one coarse cell in this mesh connects to multiple fine cells in the neighboring mesh. + ! LOCROW = LOW_IND WALL_LOOP_1 : DO IW=1,N_EXTERNAL_WALL_CELLS WC => WALL(IW) @@ -4704,19 +5036,57 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM IF(PRES_ON_WHOLE_DOMAIN) & FLG = FLG .OR. WC%BOUNDARY_TYPE==NULL_BOUNDARY .OR. (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. EWC%NOM > 0) IF ( .NOT.FLG .OR. EWC%NOM<1) CYCLE ! Here if NOM==0 means it is an OBST laying on an external boundary -> CYCLE - IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG; II = BC%II; JJ = BC%JJ; KK = BC%KK + + IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG; II = BC%II; JJ = BC%JJ; KK = BC%KK; IF(ZONE_SOLVE(PRESSURE_ZONE(IIG,JJG,KKG))%CONNECTED_ZONE_PARENT/=IPZ) CYCLE - ! Check if CC_IBM -> If either IIG,JJG,KKG or II,JJ,KK cell is type IS_CUTCFE or IS_SOLID cycle: + + ! Check if CC_IBM -> If IIG,JJG,KKG cell is type IS_CUTCFE or IS_SOLID cycle: + ! for cut-faces, RC_faces this will be taken care of in GET_CC_MATRIXGRAPH_H IF (CC_IBM) THEN - IF(CCVAR(II ,JJ ,KK ,CC_CGSC) /= IS_GASPHASE) CYCLE - IF(CCVAR(IIG,JJG,KKG,CC_CGSC) /= IS_GASPHASE) CYCLE + IF(CCVAR(IIG,JJG,KKG,CC_CGSC) /= IS_GASPHASE) CYCLE ENDIF - ! Unknowns on related cells: - IND(LOW_IND) = CCVAR(IIG,JJG,KKG,UNKH) ! internal. - IND(HIGH_IND) = CCVAR(II,JJ,KK,UNKH) ! guard-cell. - IND_LOC(LOW_IND) = IND(LOW_IND) - ZSL%UNKH_IND(NM_START) ! All row indexes must refer to ind_loc. - IND_LOC(HIGH_IND)= IND(HIGH_IND)- ZSL%UNKH_IND(NM_START) - CALL ADD_INPLACE_NNZ_H_WHLDOM(LOCROW,LOCROW,IND,IND_LOC,IPZ) + + ! Get the neighboring mesh number: + NOM = EWC%NOM + + ! Loop over ALL cells in the neighboring mesh that share this boundary face. + ! For same refinement: IIO_MIN==IIO_MAX (single cell) + ! For 2:1 refinement: IIO_MIN to IIO_MAX spans 2 cells in one direction (4 cells for face, 8 for edge) + ! For 4:1 refinement: spans 4 cells in one direction (16 cells for face) + ! + DO KKO = EWC%KKO_MIN, EWC%KKO_MAX + DO JJO = EWC%JJO_MIN, EWC%JJO_MAX + DO IIO = EWC%IIO_MIN, EWC%IIO_MAX + + ! For UGLMAT, check if the external cell is gas phase: + IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. OMESH(NOM)%MUNKH(IIO,JJO,KKO) <= 0) CYCLE + + ! GRID REFINEMENT: Compute guard cell location in neighboring mesh based on IOR + CALL COMPUTE_GUARD_CELL_INDEXES(BC%IOR, IIO, JJO, KKO, II_NOM, JJ_NOM, KK_NOM) + + ! Check if external cell has INTERPOLATED_BOUNDARY type + ! If not INTERPOLATED_BOUNDARY, skip (no connectivity, zero gradient) + IF (ALLOCATED(OMESH(NOM)%EWC_TYPE)) THEN + IF (OMESH(NOM)%EWC_TYPE(II_NOM,JJ_NOM,KK_NOM) /= INTERPOLATED_BOUNDARY) CYCLE + ENDIF + + ! Unknown numbers on related cells: + IND(LOW_IND) = CCVAR(IIG,JJG,KKG,UNKH) ! internal cell in this mesh + IND(HIGH_IND) = OMESH(NOM)%MUNKH(IIO,JJO,KKO) ! cell in neighboring mesh + + ! Convert to local row numbering: + IND_LOC(LOW_IND) = IND(LOW_IND) - ZSL%UNKH_IND(NM_START) + IND_LOC(HIGH_IND) = IND(HIGH_IND) - ZSL%UNKH_IND(NM_START) + + ! Add this connectivity to the matrix graph. + ! This function adds IND(HIGH_IND) to the column list (JD array) of row IND(LOW_IND), + ! and vice versa for symmetric matrices. + CALL ADD_INPLACE_NNZ_H_WHLDOM(LOCROW,LOCROW,IND,IND_LOC,IPZ) + + ENDDO ! IIO loop + ENDDO ! JJO loop + ENDDO ! KKO loop + ENDDO WALL_LOOP_1 ! Finally Add nonzeros corresponding to RC_FACE, CUT_FACE @@ -4784,44 +5154,6 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM NULLIFY(RGF) ENDDO AXIS_LOOP_2 - ! Now Wall faces column locations: - LOCROW = LOW_IND - WALL_LOOP_2 : DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC => WALL(IW) - BC => BOUNDARY_COORD(WC%BC_INDEX) - WC_JD(1:2,1:2) = IS_UNDEFINED - IF (.NOT.PRES_ON_WHOLE_DOMAIN .AND. WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE - IF(IW <= N_EXTERNAL_WALL_CELLS) THEN - ! Here if NOM==0 means it is an OBST laying on an external boundary -> CYCLE - EWC=>EXTERNAL_WALL(IW); IF(EWC%NOM < 1) CYCLE - ENDIF - IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG; II = BC%II; JJ = BC%JJ; KK = BC%KK - IF((.NOT.CC_IBM .AND. CCVAR(IIG,JJG,KKG,UNKH)<1) & - .OR. ZONE_SOLVE(PRESSURE_ZONE(IIG,JJG,KKG))%CONNECTED_ZONE_PARENT/=IPZ) CYCLE - ! Check if CC_IBM -> If either IIG,JJG,KKG or II,JJ,KK cell is type IS_CUTCFE or IS_SOLID cycle: - IF (CC_IBM) THEN - IF(CCVAR(II ,JJ ,KK ,CC_CGSC) /= IS_GASPHASE) CYCLE - IF(CCVAR(IIG,JJG,KKG,CC_CGSC) /= IS_GASPHASE) CYCLE - ENDIF - ! Unknowns on related cells: - IND(LOW_IND) = CCVAR(IIG,JJG,KKG,UNKH) ! internal. - IND(HIGH_IND) = CCVAR(II,JJ,KK,UNKH) ! guard-cell. - IND_LOC(LOW_IND) = IND(LOW_IND) - ZSL%UNKH_IND(NM_START) ! All row indexes must refer to ind_loc. - IND_LOC(HIGH_IND)= IND(HIGH_IND)- ZSL%UNKH_IND(NM_START) - DO IIND=LOW_IND,HIGH_IND - NII = ZSL%ROW_H(IND_LOC(LOCROW))%NNZ - DO ILOC=1,NII - IF ( IND(IIND) == ZSL%ROW_H(IND_LOC(LOCROW))%JD(ILOC) ) THEN - WC_JD(LOCROW,IIND) = ILOC - EXIT - ENDIF - ENDDO - ENDDO - WC%JD11_INDEX = WC_JD(1,1) - WC%JD12_INDEX = WC_JD(1,2) - WC%JD21_INDEX = WC_JD(2,1) - WC%JD22_INDEX = WC_JD(2,2) - ENDDO WALL_LOOP_2 ! Finally cut-face: IF (CC_IBM) CALL GET_CC_MATRIXGRAPH_H(NM,NM_START,IPZ,.FALSE.) ENDDO MESH_LOOP_2 @@ -5362,7 +5694,7 @@ SUBROUTINE FINISH_GLMAT_SOLVER #endif ! Local variables: -INTEGER :: MAXFCT, MNUM, PHASE, NRHS, ERROR, MSGLVL +INTEGER :: MAXFCT, MNUM, PHASE, NRHS, ERROR, MSGLVL, NOM #ifdef WITH_MKL INTEGER :: PERM(1) #endif @@ -5397,6 +5729,13 @@ SUBROUTINE FINISH_GLMAT_SOLVER #endif ENDIF +! Deallocate OMESH arrays used for grid refinement: +DO NOM=1,NMESHES + IF (ALLOCATED(OMESH(NOM)%GSCH)) DEALLOCATE(OMESH(NOM)%GSCH) + IF (ALLOCATED(OMESH(NOM)%MUNKH)) DEALLOCATE(OMESH(NOM)%MUNKH) + IF (ALLOCATED(OMESH(NOM)%EWC_TYPE)) DEALLOCATE(OMESH(NOM)%EWC_TYPE) +ENDDO + DEALLOCATE(ZONE_SOLVE) RETURN diff --git a/Source/type.f90 b/Source/type.f90 index a13c1211546..767a6276588 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -1047,7 +1047,7 @@ MODULE TYPES INTEGER, ALLOCATABLE, DIMENSION(:) :: ICC_UNKZ_CT_S, ICC_UNKZ_CC_S, ICC_UNKZ_CT_R, ICC_UNKZ_CC_R,& ICF_UFFB_CF_S, ICF_UFFB_CF_R INTEGER, ALLOCATABLE, DIMENSION(:) :: UNKZ_CT_S, UNKZ_CC_S, UNKZ_CT_R, UNKZ_CC_R - INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: MUNKH,GSCH + INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: MUNKH,GSCH,EWC_TYPE ! Face variables data (velocities): INTEGER, ALLOCATABLE, DIMENSION(:) :: IIO_FC_R,JJO_FC_R,KKO_FC_R,AXS_FC_R,IIO_FC_S,JJO_FC_S,KKO_FC_S,AXS_FC_S