diff --git a/Build/Scripts/set_compilers.sh b/Build/Scripts/set_compilers.sh index 5137d1bb7e3..11bc912ae40 100755 --- a/Build/Scripts/set_compilers.sh +++ b/Build/Scripts/set_compilers.sh @@ -54,7 +54,11 @@ set_compiler_from_env_var COMP_CXX FIREMODELS_CXX set_COMP_CXX set_compiler_from_env_var COMP_FC FIREMODELS_FC set_COMP_FC # Determine compiler list based on build target -if [[ "$FDS_BUILD_TARGET" == *"osx"* ]]; then +if [[ "$FDS_BUILD_TARGET" == "ompi_intel"* ]]; then + select_compiler_from_system COMP_CC mpicc icx + select_compiler_from_system COMP_CXX mpicxx icpx + select_compiler_from_system COMP_FC mpifort +elif [[ "$FDS_BUILD_TARGET" == *"osx"* ]]; then select_compiler_from_system COMP_CC mpicc clang gcc select_compiler_from_system COMP_CXX mpicxx clang++ g++ select_compiler_from_system COMP_FC mpifort diff --git a/Source/pres.f90 b/Source/pres.f90 index 901bff5bdae..e9ae0fe90fd 100644 --- a/Source/pres.f90 +++ b/Source/pres.f90 @@ -3059,6 +3059,14 @@ MODULE GLOBMAT_SOLVER ! Handle for ZONE_SOLVER array entries: TYPE(ZONE_SOLVE_TYPE), POINTER :: ZSL +! Communicator for each zone: +#ifdef WITH_HYPRE +TYPE ZSL_COMM_TYPE + TYPE(MPI_COMM) :: COMM +END TYPE ZSL_COMM_TYPE +TYPE(ZSL_COMM_TYPE), ALLOCATABLE, DIMENSION(:) :: ZSL_COMM +#endif + ! Matrix types: INTEGER, PARAMETER :: SYMM_INDEFINITE =-2 INTEGER, PARAMETER :: SYMM_POSITIVE_DEFINITE= 2 @@ -3209,15 +3217,17 @@ SUBROUTINE GLMAT_SOLVER(T,DT) #ifdef WITH_HYPRE IF (ZSL%MTYPE==SYMM_INDEFINITE .AND. ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) ZSL%F_H(MAX(ZSL%NUNKH_LOCAL,1)) = 0._EB ! Solve the system using HYPRE: - CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%F_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%F_H, HYPRE_IERR) - CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSOLVE(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%PARCSR_A_H, ZSL%HYPRE_ZSL%PAR_F_H, & - ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR) - IF (CHECK_POISSON .AND. HYPRE_SOLVER_SETPRINTLEVEL>0) THEN - CALL HYPRE_PARCSRPCGGETNUMITERATIONS(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%NUM_ITERATIONS, HYPRE_IERR) - CALL HYPRE_PARCSRPCGGETFINALRELATIVE(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%FINAL_RES_NORM, HYPRE_IERR) + IF(ZSL%NUNKH_LOCAL > 0) THEN + CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%F_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%F_H, HYPRE_IERR) + CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSOLVE(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%PARCSR_A_H, ZSL%HYPRE_ZSL%PAR_F_H, & + ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR) + IF (CHECK_POISSON .AND. HYPRE_SOLVER_SETPRINTLEVEL>0) THEN + CALL HYPRE_PARCSRPCGGETNUMITERATIONS(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%NUM_ITERATIONS, HYPRE_IERR) + CALL HYPRE_PARCSRPCGGETFINALRELATIVE(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%FINAL_RES_NORM, HYPRE_IERR) + ENDIF + CALL HYPRE_IJVECTORGETVALUES(ZSL%HYPRE_ZSL%X_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%X_H, HYPRE_IERR) ENDIF - CALL HYPRE_IJVECTORGETVALUES(ZSL%HYPRE_ZSL%X_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%X_H, HYPRE_IERR) #endif END SELECT LIBRARY_SELECT @@ -3931,8 +3941,7 @@ SUBROUTINE GET_H_MATRIX_LUDCMP INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MB_FACTOR #endif #ifdef WITH_HYPRE -INTEGER :: ONEV(1), END_ROW -REAL(EB) :: ZEROV(1) +INTEGER ::END_ROW, COLOR CHARACTER(FN_LENGTH) :: FN_PARCSRPCG_MATRIX #endif !.. All other variables @@ -3963,6 +3972,7 @@ SUBROUTINE GET_H_MATRIX_LUDCMP STOP_STATUS = SETUP_STOP RETURN ENDIF + ALLOCATE(ZSL_COMM(0:N_ZONE_GLOBMAT)) #endif END SELECT @@ -4104,97 +4114,100 @@ SUBROUTINE GET_H_MATRIX_LUDCMP IF (ALLOCATED(ZSL%HYPRE_ZSL%INDICES)) DEALLOCATE(ZSL%HYPRE_ZSL%INDICES) ALLOCATE( ZSL%HYPRE_ZSL%INDICES(MAX(1,ZSL%NUNKH_LOCAL)) ) - CALL HYPRE_IJMATRIXCREATE(MPI_COMM_WORLD,ZSL%LOWER_ROW-1,ZSL%UPPER_ROW-1,& - ZSL%LOWER_ROW-1,ZSL%UPPER_ROW-1,ZSL%HYPRE_ZSL%A_H,HYPRE_IERR) - CALL HYPRE_IJMATRIXSETOBJECTTYPE(ZSL%HYPRE_ZSL%A_H,HYPRE_PARCSR,HYPRE_IERR) - CALL HYPRE_IJMATRIXINITIALIZE(ZSL%HYPRE_ZSL%A_H,HYPRE_IERR) - IF (ZSL%MTYPE==SYMM_INDEFINITE) THEN - IF(ZSL%NUNKH_TOTAL==1) THEN ! Single unknown, zero coefficient matrix (1,1). Set coefficient to 1. - IF(ZSL%NUNKH_LOCAL==1) THEN - ZSL%NNZ_D_MAT_H(1) = 1 - DEALLOCATE(ZSL%JD_MAT_H); ALLOCATE(ZSL%JD_MAT_H(1,1)); ZSL%JD_MAT_H(1,1) = 1 - DEALLOCATE(ZSL%D_MAT_H ); ALLOCATE(ZSL%D_MAT_H(1,1) ); ZSL%D_MAT_H(1,1) = 1._EB - ENDIF - ELSE ! More than one unknown - IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN - END_ROW = ZSL%NUNKH_LOCAL-1 - ELSE - END_ROW = ZSL%NUNKH_LOCAL - ENDIF - ! Rows 1 to ZM%NUNKH-1, last column, set all to zero: - DO IROW=1,END_ROW - DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW) - IF ( ZSL%JD_MAT_H(JCOL,IROW) /= ZSL%NUNKH_TOTAL ) CYCLE ! Make zero matrix entries in last column. - ZSL%D_MAT_H(JCOL,IROW) = 0._EB - ENDDO - ENDDO - ! Last row, all zeros except the diagonal that keeps diagonal number: - IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN - IROW = ZSL%NUNKH_LOCAL - DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW) - IF ( ZSL%JD_MAT_H(JCOL,IROW) /= ZSL%NUNKH_TOTAL ) ZSL%D_MAT_H(JCOL,IROW) = 0._EB + COLOR = MPI_UNDEFINED; IF(ZSL%NUNKH_LOCAL > 0) COLOR = 1 + CALL MPI_COMM_SPLIT(MPI_COMM_WORLD,COLOR,0,ZSL_COMM(IPZ)%COMM,HYPRE_IERR) + + IF(ZSL%NUNKH_LOCAL > 0) THEN + CALL HYPRE_IJMATRIXCREATE(ZSL_COMM(IPZ)%COMM,ZSL%LOWER_ROW-1,ZSL%UPPER_ROW-1,& + ZSL%LOWER_ROW-1,ZSL%UPPER_ROW-1,ZSL%HYPRE_ZSL%A_H,HYPRE_IERR) + + CALL HYPRE_IJMATRIXSETOBJECTTYPE(ZSL%HYPRE_ZSL%A_H,HYPRE_PARCSR,HYPRE_IERR) + CALL HYPRE_IJMATRIXINITIALIZE(ZSL%HYPRE_ZSL%A_H,HYPRE_IERR) + + IF (ZSL%MTYPE==SYMM_INDEFINITE) THEN + IF(ZSL%NUNKH_TOTAL==1) THEN ! Single unknown, zero coefficient matrix (1,1). Set coefficient to 1. + IF(ZSL%NUNKH_LOCAL==1) THEN + ZSL%NNZ_D_MAT_H(1) = 1 + DEALLOCATE(ZSL%JD_MAT_H); ALLOCATE(ZSL%JD_MAT_H(1,1)); ZSL%JD_MAT_H(1,1) = 1 + DEALLOCATE(ZSL%D_MAT_H ); ALLOCATE(ZSL%D_MAT_H(1,1) ); ZSL%D_MAT_H(1,1) = 1._EB + ENDIF + ELSE ! More than one unknown + IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN + END_ROW = ZSL%NUNKH_LOCAL-1 + ELSE + END_ROW = ZSL%NUNKH_LOCAL + ENDIF + ! Rows 1 to ZM%NUNKH-1, last column, set all to zero: + DO IROW=1,END_ROW + DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW) + IF ( ZSL%JD_MAT_H(JCOL,IROW) /= ZSL%NUNKH_TOTAL ) CYCLE ! Make zero matrix entries in last column. + ZSL%D_MAT_H(JCOL,IROW) = 0._EB + ENDDO ENDDO + ! Last row, all zeros except the diagonal that keeps diagonal number: + IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN + IROW = ZSL%NUNKH_LOCAL + DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW) + IF ( ZSL%JD_MAT_H(JCOL,IROW) /= ZSL%NUNKH_TOTAL ) ZSL%D_MAT_H(JCOL,IROW) = 0._EB + ENDDO + ENDIF ENDIF ENDIF - ENDIF - IF(ZSL%NUNKH_LOCAL > 0) THEN + ZSL%JD_MAT_H = ZSL%JD_MAT_H - 1 DO IROW=1,ZSL%NUNKH_LOCAL ZSL%HYPRE_ZSL%INDICES(IROW)=ZSL%LOWER_ROW+IROW-2 CALL HYPRE_IJMATRIXSETVALUES(ZSL%HYPRE_ZSL%A_H, 1, ZSL%NNZ_D_MAT_H(IROW), ZSL%HYPRE_ZSL%INDICES(IROW), & ZSL%JD_MAT_H(1:ZSL%NNZ_D_MAT_H(IROW),IROW), ZSL%D_MAT_H(1:ZSL%NNZ_D_MAT_H(IROW),IROW), HYPRE_IERR) ENDDO - ELSE - ZSL%HYPRE_ZSL%INDICES(1)=ZSL%LOWER_ROW-1 - ONEV(1) = 1; ZEROV(1) = 0._EB - ! Define a zero coefficient in A(ZSL%LOWER_ROW-1,ZSL%LOWER_ROW-1): - CALL HYPRE_IJMATRIXSETVALUES(ZSL%HYPRE_ZSL%A_H, 1, ONEV(1), ZSL%HYPRE_ZSL%INDICES(1), & - ZSL%HYPRE_ZSL%INDICES(1), ZEROV(1), HYPRE_IERR) - ENDIF - CALL HYPRE_IJMATRIXASSEMBLE(ZSL%HYPRE_ZSL%A_H, HYPRE_IERR) - CALL HYPRE_IJMATRIXGETOBJECT(ZSL%HYPRE_ZSL%A_H, ZSL%HYPRE_ZSL%PARCSR_A_H, HYPRE_IERR) - ! Create right hand side vector - CALL HYPRE_IJVECTORCREATE(MPI_COMM_WORLD, ZSL%LOWER_ROW-1, ZSL%UPPER_ROW-1, ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) - CALL HYPRE_IJVECTORSETOBJECTTYPE(ZSL%HYPRE_ZSL%F_H, HYPRE_PARCSR, HYPRE_IERR) - CALL HYPRE_IJVECTORINITIALIZE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) - ! Create solution vector - CALL HYPRE_IJVECTORCREATE(MPI_COMM_WORLD, ZSL%LOWER_ROW-1, ZSL%UPPER_ROW-1, ZSL%HYPRE_ZSL%X_H, HYPRE_IERR) - CALL HYPRE_IJVECTORSETOBJECTTYPE(ZSL%HYPRE_ZSL%X_H, HYPRE_PARCSR, HYPRE_IERR) - CALL HYPRE_IJVECTORINITIALIZE(ZSL%HYPRE_ZSL%X_H, HYPRE_IERR) - ! Set values - CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%F_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%F_H, HYPRE_IERR) - CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%X_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%X_H, HYPRE_IERR) - ! Assemble vectors - CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) - CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%X_H, HYPRE_IERR) - ! Get rhs and soln objects - CALL HYPRE_IJVECTORGETOBJECT(ZSL%HYPRE_ZSL%F_H, ZSL%HYPRE_ZSL%PAR_F_H, HYPRE_IERR) - CALL HYPRE_IJVECTORGETOBJECT(ZSL%HYPRE_ZSL%X_H, ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR) - ! Create solver (Parallel Compressed Sparse Row Preconditioned Conjugate Gradient) - CALL HYPRE_PARCSRPCGCREATE(MPI_COMM_WORLD, ZSL%HYPRE_ZSL%SOLVER, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSETMAXITER(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_MAXIT, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSETTOL(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_TOL, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSETTWONORM(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETTWONORM, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSETPRINTLEVEL(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETPRINTLEVEL, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSETLOGGING(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETLOGGING, HYPRE_IERR) - - ! Set up the Algebraic Multi-Grid (AMG) preconditioner and specify any parameters - CALL HYPRE_BOOMERAMGCREATE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_IERR) - CALL HYPRE_BOOMERAMGSETPRINTLEVEL(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_SETPRINTLEVEL, HYPRE_IERR) - CALL HYPRE_BOOMERAMGSETCOARSENTYPE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_COARSENINGTYPE, HYPRE_IERR) - CALL HYPRE_BOOMERAMGSETRELAXTYPE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_SETRELAXTYPE, HYPRE_IERR) - CALL HYPRE_BOOMERAMGSETNUMSWEEPS(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_NUMSWEEPS, HYPRE_IERR) - CALL HYPRE_BOOMERAMGSETTOL(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_TOL, HYPRE_IERR) - CALL HYPRE_BOOMERAMGSETMAXITER(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_MAXITER, HYPRE_IERR) - CALL HYPRE_PARCSRPCGSETPRECOND(ZSL%HYPRE_ZSL%SOLVER, HYPRE_PRECOND_ID, ZSL%HYPRE_ZSL%PRECOND, HYPRE_IERR) - ! Solver setup - IF(WRITE_PARCSRPCG_MATRIX) THEN - WRITE(FN_PARCSRPCG_MATRIX,'(A,A)') TRIM(CHID),'_hypre_matrix.txt' - CALL HYPRE_PARCSRMATRIXPRINT(ZSL%HYPRE_ZSL%PARCSR_A_H, TRIM(FN_PARCSRPCG_MATRIX), LEN(TRIM(FN_PARCSRPCG_MATRIX)), HYPRE_IERR) + CALL HYPRE_IJMATRIXASSEMBLE(ZSL%HYPRE_ZSL%A_H, HYPRE_IERR) + CALL HYPRE_IJMATRIXGETOBJECT(ZSL%HYPRE_ZSL%A_H, ZSL%HYPRE_ZSL%PARCSR_A_H, HYPRE_IERR) + + ! Create right hand side vector + CALL HYPRE_IJVECTORCREATE(ZSL_COMM(IPZ)%COMM, ZSL%LOWER_ROW-1, ZSL%UPPER_ROW-1, ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) + CALL HYPRE_IJVECTORSETOBJECTTYPE(ZSL%HYPRE_ZSL%F_H, HYPRE_PARCSR, HYPRE_IERR) + CALL HYPRE_IJVECTORINITIALIZE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) + ! Create solution vector + CALL HYPRE_IJVECTORCREATE(ZSL_COMM(IPZ)%COMM, ZSL%LOWER_ROW-1, ZSL%UPPER_ROW-1, ZSL%HYPRE_ZSL%X_H, HYPRE_IERR) + CALL HYPRE_IJVECTORSETOBJECTTYPE(ZSL%HYPRE_ZSL%X_H, HYPRE_PARCSR, HYPRE_IERR) + CALL HYPRE_IJVECTORINITIALIZE(ZSL%HYPRE_ZSL%X_H, HYPRE_IERR) + ! Set values + CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%F_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%F_H, HYPRE_IERR) + CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%X_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%X_H, HYPRE_IERR) + ! Assemble vectors + CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR) + CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%X_H, HYPRE_IERR) + ! Get rhs and soln objects + CALL HYPRE_IJVECTORGETOBJECT(ZSL%HYPRE_ZSL%F_H, ZSL%HYPRE_ZSL%PAR_F_H, HYPRE_IERR) + CALL HYPRE_IJVECTORGETOBJECT(ZSL%HYPRE_ZSL%X_H, ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR) + + ! Create solver (Parallel Compressed Sparse Row Preconditioned Conjugate Gradient) + CALL HYPRE_PARCSRPCGCREATE(ZSL_COMM(IPZ)%COMM, ZSL%HYPRE_ZSL%SOLVER, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSETMAXITER(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_MAXIT, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSETTOL(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_TOL, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSETTWONORM(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETTWONORM, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSETPRINTLEVEL(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETPRINTLEVEL, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSETLOGGING(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETLOGGING, HYPRE_IERR) + + ! Set up the Algebraic Multi-Grid (AMG) preconditioner and specify any parameters + CALL HYPRE_BOOMERAMGCREATE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_IERR) + CALL HYPRE_BOOMERAMGSETPRINTLEVEL(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_SETPRINTLEVEL, HYPRE_IERR) + CALL HYPRE_BOOMERAMGSETCOARSENTYPE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_COARSENINGTYPE, HYPRE_IERR) + CALL HYPRE_BOOMERAMGSETRELAXTYPE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_SETRELAXTYPE, HYPRE_IERR) + CALL HYPRE_BOOMERAMGSETNUMSWEEPS(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_NUMSWEEPS, HYPRE_IERR) + CALL HYPRE_BOOMERAMGSETTOL(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_TOL, HYPRE_IERR) + CALL HYPRE_BOOMERAMGSETMAXITER(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_MAXITER, HYPRE_IERR) + CALL HYPRE_PARCSRPCGSETPRECOND(ZSL%HYPRE_ZSL%SOLVER, HYPRE_PRECOND_ID, ZSL%HYPRE_ZSL%PRECOND, HYPRE_IERR) + ! Solver setup + IF(WRITE_PARCSRPCG_MATRIX) THEN + WRITE(FN_PARCSRPCG_MATRIX,'(A,A)') TRIM(CHID),'_hypre_matrix.txt' + CALL HYPRE_PARCSRMATRIXPRINT(ZSL%HYPRE_ZSL%PARCSR_A_H, TRIM(FN_PARCSRPCG_MATRIX), & + LEN(TRIM(FN_PARCSRPCG_MATRIX)), HYPRE_IERR) + ENDIF + CALL HYPRE_PARCSRPCGSETUP(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%PARCSR_A_H,& + ZSL%HYPRE_ZSL%PAR_F_H, ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR) ENDIF - CALL HYPRE_PARCSRPCGSETUP(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%PARCSR_A_H,& - ZSL%HYPRE_ZSL%PAR_F_H, ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR) #endif END SELECT LIBRARY_SELECT @@ -5377,6 +5390,7 @@ SUBROUTINE FINISH_GLMAT_SOLVER IF (UGLMAT_SOLVER_LIBRARY==HYPRE_FLAG) THEN #ifdef WITH_HYPRE CALL HYPRE_FINALIZE(HYPRE_IERR) + DEALLOCATE(ZSL_COMM) #endif ENDIF diff --git a/Utilities/Misc/geometry/colocated_schemes/README b/Utilities/Misc/geometry/colocated_schemes/README new file mode 100644 index 00000000000..2db81cb8dc5 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/README @@ -0,0 +1,9 @@ +! Routines to test colocated schemes for Low Mach model time integration. +! Test cases: +! 1. taylor_green_vortex.m : Test overall accuracy on structured grids with Nx x Ny cells. +! Same parameters as ns2d set of verification cases in FDS_Verification_Guide.pdf. +! Performs 2 tests : a) L2 error norm of momentum at T_end over the whole domain. +! b) RMS error of U-Velocity in center point over all time integration steps, +! (same as done in the FDS verification guide). +! +! ------------------------------------------------------------------------------------------------ diff --git a/Utilities/Misc/geometry/colocated_schemes/build_FB_cell_face.m b/Utilities/Misc/geometry/colocated_schemes/build_FB_cell_face.m new file mode 100644 index 00000000000..b7b6e49a0e9 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/build_FB_cell_face.m @@ -0,0 +1,48 @@ +function [FB_cell, FB_face_n] = build_FB_cell_face(p_cell, rho_cell, cells, faces, Lx, Ly) +% p_cell : Nc x 1 pressure (or Bernoulli-pressure part) at cells +% rho_cell : Nc x 1 density at cells +% Outputs: +% FB_cell : Nc x 2 vector -p * ∇(1/ρ) at cells +% FB_face_n : Nf x 1 scalar face-normal avg: 0.5(FB_P+FB_N)·n_f + +Nc = numel(cells); Nf = numel(faces); + +invR = 1 ./ rho_cell; % 1/ρ at cells +gradInvR = ls_grad_scalar_2d(invR, cells, faces, Lx, Ly); % ∇(1/ρ) at cells + +% Cell-centered F_B = - p * ∇(1/ρ) +FB_cell = - p_cell(:) .* gradInvR; % implicit expansion: (Nc×1) .* (Nc×2) → (Nc×2) + +% Face-normal average like F_A +FB_face_n = zeros(Nf,1); +for f = 1:Nf + P = faces(f).owner; + N = faces(f).neigh; % periodic mesh ⇒ N>0 + nf = faces(f).nf(:); + if N>0 + FBav = 0.5*(FB_cell(P,:) + FB_cell(N,:)); + else + FBav = FB_cell(P,:); % (kept for completeness) + end + FB_face_n(f) = FBav * nf; +end +end + +% function FB_face_n = build_FB_face_n_compact(p, rho, cells, faces, Lx, Ly) +% Nf = numel(faces); +% FB_face_n = zeros(Nf,1); +% for f = 1:Nf +% P = faces(f).owner; +% N = faces(f).neigh; % periodic ⇒ N>0 +% nf = faces(f).nf(:); +% +% dvec = periodic_delta_2d(cells(N).xc - cells(P).xc, Lx, Ly); +% dn = dot(dvec, nf); % compact face distance along nf +% invrP = 1.0 / rho(P); +% invrN = 1.0 / rho(N); +% d_invrho_dn = (invrN - invrP) / dn; +% +% pf = 0.5*(p(P) + p(N)); % centered face pressure +% FB_face_n(f) = - pf * d_invrho_dn; +% end +% end diff --git a/Utilities/Misc/geometry/colocated_schemes/build_adv_body_forces_2d.m b/Utilities/Misc/geometry/colocated_schemes/build_adv_body_forces_2d.m new file mode 100644 index 00000000000..0e87709e026 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/build_adv_body_forces_2d.m @@ -0,0 +1,77 @@ +function [FA_cell, FA_face_n] = build_adv_body_forces_2d(u_in, rho, rho0, gvec, U_f, cells, faces, stage) +% F_A = div(u⊗u) - grad(0.5|u|^2) - (div u) u - (1/rho)(rho-rho0) g +% u : Nc x 2 (cell-centered) +% rho : Nc x 1 +% U_f : Nf x 1 (face-normal velocity along faces(f).nf, owner->neigh) +% returns: +% FA_cell : Nc x 2 (cell-centered vector) +% FA_face_n : Nf x 1 (face-normal avg of F_A for Poisson RHS) + +global H_projection + +%if stage==2 +% u = reconstruct_u_from_faces_w(U_f, cells, faces); +%else + u = u_in; +%end +Nc = size(u,1); Nf = numel(faces); +Kc = 0.5*sum(u.^2,2); % 0.5|u|^2 at cells + + + +sum_divuV = zeros(Nc,1); % Σ s_pf U_f A_f +sum_divUU = zeros(Nc,2); % Σ s_pf (u_f^c U_f) A_f +sum_gradKV= zeros(Nc,2); % Σ s_pf (K_f^c n_f) A_f + +for f = 1:Nf + P = faces(f).owner; N = faces(f).neigh; + nf = faces(f).nf(:); Af = faces(f).Af; + Uf = U_f(f); + + % centered face values (periodic mesh ⇒ N>0) + if N>0 + uf = 0.5*(u(P,:)+u(N,:)); + Kf = 0.5*(Kc(P)+Kc(N)); + else + uf = u(P,:); Kf = Kc(P); + end + + adv_flux = (uf .* Uf) * Af; % vector + gradK_flux = (Kf * nf.') * Af; % vector + + % owner (+), neighbor (-) + sum_divUU(P,:) = sum_divUU(P,:) + adv_flux; + sum_gradKV(P,:)= sum_gradKV(P,:) + gradK_flux; + sum_divuV(P) = sum_divuV(P) + Uf*Af; + + if N>0 + sum_divUU(N,:) = sum_divUU(N,:) - adv_flux; + sum_gradKV(N,:)= sum_gradKV(N,:) - gradK_flux; + sum_divuV(N) = sum_divuV(N) - Uf*Af; + end +end + +FA_cell = zeros(Nc,2); +g = gvec(:).'; +for p = 1:Nc + Vp = cells(p).V; + divUU = sum_divUU(p,:)/Vp; % ∇·(u⊗u) + if(H_projection) + gradK = sum_gradKV(p,:)/Vp; % ∇(0.5|u|^2) + divu_u = (sum_divuV(p)/Vp)*u(p,:); % (∇·u) u + else + gradK = zeros(1,2); + divu_u = zeros(1,2); + end + body = - ((rho(p)-rho0)/rho(p)) * g; + FA_cell(p,:) = divUU - gradK - divu_u + body; +end + +% face-normal average of F_A (for Poisson RHS) +FA_face_n = zeros(Nf,1); +for f = 1:Nf + P = faces(f).owner; N = faces(f).neigh; nf = faces(f).nf(:); + Fav = (N>0) * 0.5*(FA_cell(P,:)+FA_cell(N,:)) + (N==0)*FA_cell(P,:); + FA_face_n(f) = Fav * nf; +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/build_poisson_matrix_compact.m b/Utilities/Misc/geometry/colocated_schemes/build_poisson_matrix_compact.m new file mode 100644 index 00000000000..9e3b92314b0 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/build_poisson_matrix_compact.m @@ -0,0 +1,87 @@ +function A = build_poisson_matrix_compact(cells, faces, Lx, Ly) +% Returns sparse A (Nc x Nc) for: sum_f ((H_N - H_P)/dn)*Af = RHS + +do_A_tests = 0; % If 1, do sanity tests on matrix. + +Nc = numel(cells); +Nf = numel(faces); + +% Compact (two-point) Poisson matrix in flux form +% A*phi produces per-cell SUM_f (Af/|dPN|) * (phi_P - phi_N), i.e., flux balance. +I = zeros(4*Nf,1); J = I; V = I; m = 0; + +for f = 1:Nf + P = faces(f).owner; + N = faces(f).neigh; % assume periodic/internal ⇒ N>0 + if N<=0, continue; end + Af = faces(f).Af; + nf = faces(f).nf(:); % unit normal, owner→neighbor + + % owner→neighbor minimal-image displacement (periodic safe) + dvec = periodic_delta_2d(cells(N).xc - cells(P).xc, Lx, Ly); + dn = dot(dvec, nf); % can be ± + %dPN = max(eps, abs(dn)); % ALWAYS use absolute distance + T = Af / dn; %PN; % transmissibility (≥0) + + % symmetric 2x2 face stencil: + % row P: +T at P, -T at N + m=m+1; I(m)=P; J(m)=P; V(m)= T; + m=m+1; I(m)=P; J(m)=N; V(m)=-T; + % row N: -T at P, +T at N + m=m+1; I(m)=N; J(m)=N; V(m)= T; + m=m+1; I(m)=N; J(m)=P; V(m)=-T; +end + +A = -sparse(I(1:m), J(1:m), V(1:m), Nc, Nc); + + +% Tests: +if(do_A_tests) +% symmetry +fprintf('||A - A''|| = %.3e\n', norm(A - A','fro')); + +% nullspace (periodic): one zero eigenvalue (the constant) +lam = eigs(A,4,'smallestreal'); +disp(sort(real(lam)).'); % expect ~[0, +, +, +] + +z = A*ones(numel(cells),1); +fprintf('||A*1||_2 = %.3e\n', norm(z,2)); + + +% divergence-of-gradient pairing (with compact normals) +H = rand(Nc,1); H = H - mean(H); +[dHdn_f, ~] = grad_from_compact_normals(H, cells, faces, Lx, Ly); +LHS = A*H; % face flux balance per cell +RHS = zeros(Nc,1); +for p=1:Nc + acc=0; for f=cells(p).faces + s = +1; if faces(f).neigh==p, s=-1; end + acc = acc + s*dHdn_f(f) * faces(f).Af; + end + RHS(p) = acc; +end +fprintf('pairing ||A*H - sum_f (dHdn_f Af s_pf)||_2 = %.3e\n', norm(LHS-RHS,2)); + +bad = 0; +for f=1:Nf + P=faces(f).owner; N=faces(f).neigh; + dvec = periodic_delta_2d(cells(N).xc - cells(P).xc, Lx, Ly); + if dot(dvec, faces(f).nf(:)) <= 0, bad=bad+1; end +end +fprintf('faces with nf·(xN-xP) <= 0: %d\n', bad); + +% 1) normals point owner->neighbor on average (not required to be perfectly aligned) +bad = 0; +for f = 1:numel(faces) + P=faces(f).owner; N=faces(f).neigh; + if N<=0, continue; end + dvec = periodic_delta_2d(cells(N).xc - cells(P).xc, Lx, Ly); + if dot(dvec, faces(f).nf(:)) < 0, bad = bad + 1; end +end +fprintf('faces with (xN-xP)·nf < 0: %d\n', bad); +% It’s OK if some are negative (skew), but we rely on ABS(dPN) in both A and q_out. + +% 2) row-sum (nullspace) check: should be 0 for unpinned periodic A +fprintf('||A*1||_2 = %.3e\n', norm(A*ones(numel(cells),1),2)); +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/build_poisson_rhs_H.m b/Utilities/Misc/geometry/colocated_schemes/build_poisson_rhs_H.m new file mode 100644 index 00000000000..82584ae1bfb --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/build_poisson_rhs_H.m @@ -0,0 +1,40 @@ +function rhs = build_poisson_rhs_H(U_f, divU_th, dt, cells, faces, FA_face_n, FB_face_n, Fextra_face_n) +% U_f : Nf x 1 face-normal velocity (owner→neighbor) +% divU_th : Nc x 1 target thermodynamic divergence per cell +% dt : scalar +% FA_face_n : Nf x 1 face-normal advective+body (scalar). [] → zeros +% FB_face_n : Nf x 1 face-normal -p ∂n(1/ρ). [] → zeros +% Fextra_face_n: Nf x 1 optional extra face term (e.g., explicit viscous). [] → zeros +% Output: +% rhs(p) = - [ (div_th*V - Σ U_f A_f)/dt + Σ FA_face_n A_f + Σ FB_face_n A_f + Σ Fextra A_f ] + +Nc = numel(cells); Nf = numel(faces); +if nargin<6 || isempty(FA_face_n), FA_face_n = zeros(Nf,1); end +if nargin<7 || isempty(FB_face_n), FB_face_n = zeros(Nf,1); end +if nargin<8 || isempty(Fextra_face_n),Fextra_face_n= zeros(Nf,1); end + +rhs = zeros(Nc,1); + +for p = 1:Nc + Vp = cells(p).V; + S_U = 0.0; % Σ s_pf U_f A_f + S_FA = 0.0; % Σ s_pf F_A,f A_f + S_FB = 0.0; % Σ s_pf F_B,f A_f + S_FX = 0.0; % Σ s_pf Fextra,f A_f + + fids = cells(p).faces; + for k = 1:numel(fids) + f = fids(k); + Af = faces(f).Af; + + s = +1; if faces(f).neigh==p, s = -1; end % outward sign for cell p + + S_U = S_U + s * U_f(f) * Af; + S_FA = S_FA + s * FA_face_n(f) * Af; + S_FB = S_FB + s * FB_face_n(f) * Af; + S_FX = S_FX + s * Fextra_face_n(f) * Af; + end + + rhs(p) = - ( (divU_th(p)*Vp - S_U)/dt + S_FA + S_FB + S_FX ); +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/build_square_mesh_periodic_2d.m b/Utilities/Misc/geometry/colocated_schemes/build_square_mesh_periodic_2d.m new file mode 100644 index 00000000000..4a60042c3b3 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/build_square_mesh_periodic_2d.m @@ -0,0 +1,89 @@ +function [cells, faces] = build_square_mesh_periodic_2d(Nx, Ny, Lx, Ly) + dx = Lx/Nx; dy = Ly/Ny; + Nc = Nx*Ny; + xw = -Lx/2; + ys = -Ly/2; + + % --- Cells --- + cells(Nc) = struct('xc',[],'V',[],'faces',[]); + id = @(i,j) (j-1)*Nx + i; % 1-based linear index + for j = 1:Ny + yc = ys + (j-0.5)*dy; + for i = 1:Nx + xc_ij = xw + (i-0.5)*dx; + p = id(i,j); + cells(p).xc = [xc_ij, yc]; + cells(p).V = dx*dy; + cells(p).faces = []; % will fill later + end + end + + % --- Faces (right and top only; periodic wraps add neighbors) --- + faces = struct('owner',{},'neigh',{},'nf',{},'Af',{},'cf',{}); + fcount = 0; + + % 1) Vertical faces: between (i,j) and (i+1,j), wrap at i=Nx → 1 + for j = 1:Ny + yc = ys + (j-0.5)*dy; + for i = 1:Nx + P = id(i,j); + ir = i+1; if ir>Nx, ir = 1; end + N = id(ir,j); + + cf_x = xw + i*dx; if cf_x>Lx, cf_x = cf_x - Lx; end + cf = [cf_x, yc]; + nf = [1, 0]; % normal from P to N + Af = dy; + + fcount = fcount + 1; + faces(fcount).owner = P; + faces(fcount).neigh = N; + faces(fcount).nf = nf; + faces(fcount).Af = Af; + faces(fcount).cf = cf; + + cells(P).faces(end+1) = fcount; + cells(N).faces(end+1) = fcount; + end + end + + % 2) Horizontal faces: between (i,j) and (i,j+1), wrap at j=Ny → 1 + for j = 1:Ny + jr = j+1; if jr>Ny, jr = 1; end + y_face = ys + j*dy; if y_face>Ly, y_face = y_face - Ly; end + for i = 1:Nx + P = id(i,j); + Np = id(i,jr); + + xc = xw + (i-0.5)*dx; + cf = [xc, y_face]; + nf = [0, 1]; % normal from P to N + Af = dx; + + fcount = fcount + 1; + faces(fcount).owner = P; + faces(fcount).neigh = Np; + faces(fcount).nf = nf; + faces(fcount).Af = Af; + faces(fcount).cf = cf; + + cells(P).faces(end+1) = fcount; + cells(Np).faces(end+1) = fcount; + end + end + + % Build dpn: + for f = 1:numel(faces) + P = faces(f).owner; N = faces(f).neigh; + nf = faces(f).nf(:); + if N > 0 + dr_raw = cells(N).xc(:) - cells(P).xc(:); + dr = periodic_delta_2d(dr_raw, Lx, Ly); + faces(f).dPN = abs( dr' * nf ); + else + faces(f).dPN = []; % or a BC-specific value + end +end + +end + diff --git a/Utilities/Misc/geometry/colocated_schemes/build_viscous_forces_2d.m b/Utilities/Misc/geometry/colocated_schemes/build_viscous_forces_2d.m new file mode 100644 index 00000000000..1713cef27cb --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/build_viscous_forces_2d.m @@ -0,0 +1,102 @@ +function [divTau_cell, FA_visc_cell, FA_visc_face_n,traction_f] = build_viscous_forces_2d(u, rho, mu, lambda, cells, faces, Lx, Ly) +% Compute (∇·tau) at cells using face tractions of the Newtonian viscous stress. +% u: Nc x 2 velocities (u,v) at cells +% rho: Nc x 1 density +% mu: Nc x 1 dynamic viscosity at cells (can vary) +% lambda: Nc x 1 bulk viscosity at cells (if [], uses Stokes: lambda = -2/3 mu) +% cells: struct array with .xc(1x2), .V, .faces(1xnf) +% faces: struct array with .owner, .neigh (0 if boundary), .nf(1x2 unit), .Af, .cf(1x2) +% +% returns: +% divTau_cell: Nc x 2 (∇·tau) at cells +% FA_visc_face_n : Nf x 1 face-normal scalar to add into FA_face_n +% traction_f: Nf x 2 tau_f * n_f (vector without area) + +Nc = size(u,1); +Nf = numel(faces); + +if isempty(lambda) + lambda = -2/3 * mu; % Stokes' hypothesis +end + +% 1) cell-centered velocity gradients by LS +[gradU, gradV] = ls_grad_velocity_2d(u, cells, faces, Lx, Ly); +% gradU(p,:) = [du/dx, du/dy]; gradV(p,:) = [dv/dx, dv/dy] + +% 2) face traction tau_f * n_f +traction_f = zeros(Nf,2); + +for f = 1:Nf + P = faces(f).owner; + N = faces(f).neigh; + nf = faces(f).nf(:); % unit normal (P -> N) + % build face gradient by centered average (robust, symmetric) + if N>0 + Gp = [gradU(P,1), gradU(P,2); % ∇u at P (2x2) + gradV(P,1), gradV(P,2)]; % rows: du/dx du/dy ; dv/dx dv/dy + Gn = [gradU(N,1), gradU(N,2); + gradV(N,1), gradV(N,2)]; + Gf = 0.5*(Gp + Gn); + + muf = 0.5*(mu(P) + mu(N)); + lamf= 0.5*(lambda(P) + lambda(N)); + else + % boundary face: simple one-sided fallback using cell gradient + Gf = [gradU(P,1), gradU(P,2); + gradV(P,1), gradV(P,2)]; + muf = mu(P); + lamf= lambda(P); + % NOTE: For wall/inlet/outlet you can/should replace this + % with BC-consistent gradient (ghost-cell or analytic). + end + + % divergence at face (trace of gradient) + divUf = Gf(1,1) + Gf(2,2); + + % stress tensor at face: tau = mu*(Grad u + Grad u^T) + lambda*(div u) I + tau = muf*(Gf + Gf.') + lamf*divUf*eye(2); + + % traction: t_f = tau * n_f (vector 2x1) + traction_f(f,:) = (tau * nf).'; +end + +% 3) assemble (∇·tau)_cell by divergence theorem +divTau_cell = zeros(Nc,2); +for p = 1:Nc + sumF = [0,0]; + fids = cells(p).faces; + for k = 1:numel(fids) + f = fids(k); + Af = faces(f).Af; + nf = faces(f).nf(:); % orientation P->N + + if faces(f).owner == p + % outward normal of cell p is +nf + sumF = sumF + (traction_f(f,:)*Af); + elseif faces(f).neigh == p + % outward normal of cell p is -nf + sumF = sumF - (traction_f(f,:)*Af); + else + % should not happen + end + end + divTau_cell(p,:) = sumF / cells(p).V; +end + +% Cell viscous part of F_A (vector) +FA_visc_cell = - divTau_cell ./ rho; % Nc x 2, implicit expansion + +% Face-normal average +FA_visc_face_n = zeros(Nf,1); +for f = 1:Nf + P = faces(f).owner; + N = faces(f).neigh; + nf = faces(f).nf(:); + if N>0 + Fav = 0.5 * (FA_visc_cell(P,:) + FA_visc_cell(N,:)); + else + Fav = FA_visc_cell(P,:); + end + FA_visc_face_n(f) = Fav * nf; % dot product +end +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/cell_divergence_from_faces.m b/Utilities/Misc/geometry/colocated_schemes/cell_divergence_from_faces.m new file mode 100644 index 00000000000..90b1b71e07a --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/cell_divergence_from_faces.m @@ -0,0 +1,15 @@ +function divU = cell_divergence_from_faces(Uface, cells, faces) +% Returns (∇·u)_p from updated face-normal velocities +Nc = numel(cells); +divU = zeros(Nc,1); +for p = 1:Nc + S = 0.0; + for k = 1:numel(cells(p).faces) + f = cells(p).faces(k); + Af = faces(f).Af; + s = +1; if faces(f).neigh==p, s = -1; end + S = S + s * Uface(f) * Af; + end + divU(p) = S / cells(p).V; +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/check_H_divU.m b/Utilities/Misc/geometry/colocated_schemes/check_H_divU.m new file mode 100644 index 00000000000..9c6ecd3aac5 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/check_H_divU.m @@ -0,0 +1,51 @@ +% function check_H_divU(H_cell, Uface, cells, faces) +% Nc = numel(cells); V = reshape([cells.V],[],1); +% +% % cell divergence from projected face speeds +% divU = zeros(Nc,1); +% for p=1:Nc +% S=0; for k=1:numel(cells(p).faces) +% f=cells(p).faces(k); s=+1; if faces(f).neigh==p, s=-1; end +% S = S + s*Uface(f)*faces(f).Af; +% end +% divU(p) = S / V(p); +% end +% +% % face-averaged H +% Nf = numel(faces); Hf = zeros(Nf,1); +% for f=1:Nf +% P=faces(f).owner; N=faces(f).neigh; +% Hf = (N > 0) * 0.5*(H_cell(P) + H_cell(N)) + (N == 0) * H_cell(P); +% end +% +% LHS = sum( V .* (H_cell(:).*divU) ); +% RHS = sum( Hf .* Uface .* [faces.Af].' ); +% fprintf('(H,divU) face pairing: LHS=%.3e RHS=%.3e diff=%.3e\n', LHS, RHS, LHS-RHS); +% end + +function check_H_divU(H_cell, Uface, cells, faces) +Nc = numel(cells); V = reshape([cells.V],[],1); + +% cell divergence from faces (owner→neigh, s_pf = +1 for owner, -1 for neigh) +divU = zeros(Nc,1); +for p=1:Nc + acc = 0.0; + flist = cells(p).faces(:); + for k=1:numel(flist) + f = flist(k); + s = +1; if faces(f).neigh==p, s = -1; end + acc = acc + s * Uface(f) * faces(f).Af; + end + divU(p) = acc / V(p); +end + +LHS = sum( V .* (H_cell(:) .* divU) ); + +% RHS = sum_f (H_owner - H_neigh) * U_f * A_f (internal faces only) +owners = [faces.owner]'; neighs = [faces.neigh]'; Af = [faces.Af]'; +mask = (neighs > 0); +dH = H_cell(owners(mask)) - H_cell(neighs(mask)); +RHS = sum( dH .* Uface(mask) .* Af(mask) ); + +fprintf('(H,divU) pure pairing: LHS=%.6e RHS=%.6e diff=%.3e\n', LHS, RHS, LHS-RHS); +end diff --git a/Utilities/Misc/geometry/colocated_schemes/check_product_rule_cellwise.m b/Utilities/Misc/geometry/colocated_schemes/check_product_rule_cellwise.m new file mode 100644 index 00000000000..2d6be3ea39a --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/check_product_rule_cellwise.m @@ -0,0 +1,60 @@ +function check_product_rule_cellwise(u_cell, H_cell, cells, faces) +Nc = numel(cells); Nf = numel(faces); +V = reshape([cells.V],[],1); + +% Green–Gauss grad H and Hf +sum_flux = zeros(Nc,2); +Hf = zeros(Nf,1); +for f = 1:Nf + P = faces(f).owner; + N = faces(f).neigh; + nf = faces(f).nf(:); + Af = faces(f).Af; + + if N > 0 + Hf(f) = 0.5*(H_cell(P) + H_cell(N)); + else + Hf(f) = H_cell(P); + end + contrib = (Hf(f) * nf.').*Af; % 1x2 + sum_flux(P,:) = sum_flux(P,:) + contrib; + if N > 0 + sum_flux(N,:) = sum_flux(N,:) - contrib; + end +end +gradH = sum_flux ./ V; % Nc x 2 + +% RHS and divU using cell-interpolated face-normal speed +rhs_cell = zeros(Nc,1); +divU = zeros(Nc,1); +for p = 1:Nc + S_rhs = 0; S_div = 0; + fids = cells(p).faces; + for k = 1:numel(fids) + f = fids(k); + nf = faces(f).nf(:); + Af = faces(f).Af; + s = +1; if faces(f).neigh==p, s = -1; end + + P = faces(f).owner; + N = faces(f).neigh; + if N > 0 + uf = 0.5*(u_cell(P,:) + u_cell(N,:)); + else + uf = u_cell(P,:); + end + ufn = uf * nf; % scalar + + S_rhs = S_rhs + s * Hf(f) * ufn * Af; + S_div = S_div + s * ufn * Af; + end + rhs_cell(p) = S_rhs; + divU(p) = S_div / V(p); +end + +lhs_cell = V .* ( sum(u_cell .* gradH, 2) + H_cell(:).*divU ); +res_cell = rhs_cell - lhs_cell; + +fprintf('Cell product rule: L2=%.3e Linf=%.3e sum(RHS)=%.3e sum(LHS)=%.3e diff=%.3e\n', ... + norm(res_cell,2), norm(res_cell,inf), sum(rhs_cell), sum(lhs_cell), sum(rhs_cell)-sum(lhs_cell)); +end diff --git a/Utilities/Misc/geometry/colocated_schemes/compute_dt_explicit.m b/Utilities/Misc/geometry/colocated_schemes/compute_dt_explicit.m new file mode 100644 index 00000000000..4a5086aa6d4 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/compute_dt_explicit.m @@ -0,0 +1,101 @@ +function out = compute_dt_explicit(cells, faces, Uface, mu, rho, CFL, VN, Lx, Ly) +%COMPUTE_DT_EXPLICIT Stable dt for Explicit Euler on unstructured FV mesh. +% +% dt_adv_p = CFL * Vp / sum_f |U_f| Af +% dt_dif_p = VN * Vp / sum_f ( 2 * nu_f * Af / dn ) +% +% Inputs: +% cells : struct array with fields .V, .xc(1x2), .faces(1:nfp) +% faces : struct array with fields .owner, .neigh, .nf(1x2 unit), .Af, .cf(1x2) +% Uface : Nf x 1 face-normal velocities (along faces(f).nf, owner->neigh) +% mu : Nc x 1 dynamic viscosity at cells +% rho : Nc x 1 density at cells +% CFL : scalar convective Courant number (e.g. 0.5 … 0.9) +% VN : scalar Von Neumann number for diffusion (e.g. 0.5 is safe) +% Lx,Ly : domain lengths (for periodic minimum-image distances) +% +% Output 'out' struct: +% out.dt_adv_min : min convective dt over cells +% out.dt_dif_min : min diffusive dt over cells +% out.dt : min(out.dt_adv_min, out.dt_dif_min) +% out.dt_adv_cell : Nc x 1 per-cell convective dt +% out.dt_dif_cell : Nc x 1 per-cell diffusive dt +% out.limiter_cell : Nc x 1: -1 advective-limited, +1 diffusive-limited, 0 neither +% +% Notes: +% - Assumes every face has owner & neigh (periodic). For boundary faces w/ neigh=0, +% treat dn and face-averages accordingly (not needed for fully periodic meshes). +% - Make sure Uface uses the same face-normal orientation as faces(f).nf. + +Nc = numel(cells); +Nf = numel(faces); + +% Precompute per-cell accumulators +sum_absU_A = zeros(Nc,1); % Σ |U_f| Af +sum_dif = zeros(Nc,1); % Σ 2 * nu_f * Af / dn + +% Handy inline for periodic minimum-image delta +wrap = @(dr,L)(dr - L.*(dr > 0.5*L) + L.*(dr <= -0.5*L)); + +for f = 1:Nf + P = faces(f).owner; + N = faces(f).neigh; % periodic mesh => N>0 + nf = faces(f).nf(:); + Af = faces(f).Af; + + % Convective accumulator (same for both adjacent cells, with opposite sign handled later): + absU_A = abs(Uface(f)) * Af; + + % Diffusion: face nu and compact dn + if N > 0 + % arithmetic face average (ok for ν) + nu_f = 0.5*(mu(P)/rho(P) + mu(N)/rho(N)); + dvec = [cells(N).xc(1) - cells(P).xc(1), cells(N).xc(2) - cells(P).xc(2)]; + dvec = [wrap(dvec(1),Lx), wrap(dvec(2),Ly)]; + dn = max( eps, dot(dvec, nf) ); % ensure positive + else + % boundary: one-sided estimate (not used for periodic) + nu_f = (mu(P)/rho(P)); + % estimate dn from owner-centroid to face-centroid projected on nf + dvec = [faces(f).cf(1) - cells(P).xc(1), faces(f).cf(2) - cells(P).xc(2)]; + dn = max( eps, dot(dvec, nf) ); + end + + dif_f = 2.0 * nu_f * (Af / dn); + + % Add to owner and neighbor cells + sum_absU_A(P) = sum_absU_A(P) + absU_A; + sum_dif(P) = sum_dif(P) + dif_f; + + if N > 0 + sum_absU_A(N) = sum_absU_A(N) + absU_A; + sum_dif(N) = sum_dif(N) + dif_f; + end +end + +Vp = reshape([cells.V],[],1); + +% Per-cell dt’s (protect against division by zero) +dt_adv_cell = inf(Nc,1); +dt_dif_cell = inf(Nc,1); + +mask_adv = sum_absU_A > 0; +dt_adv_cell(mask_adv) = CFL * Vp(mask_adv) ./ sum_absU_A(mask_adv); + +mask_dif = sum_dif > 0; +dt_dif_cell(mask_dif) = VN * Vp(mask_dif) ./ sum_dif(mask_dif); + +% Outputs +out.dt_adv_cell = dt_adv_cell; +out.dt_dif_cell = dt_dif_cell; +out.dt_adv_min = min(dt_adv_cell); +out.dt_dif_min = min(dt_dif_cell); +out.dt = min(out.dt_adv_min, out.dt_dif_min); + +% Limiter label per cell +out.limiter_cell = zeros(Nc,1); +adv_better = dt_adv_cell < dt_dif_cell; +dif_better = dt_dif_cell < dt_adv_cell; +out.limiter_cell(adv_better) = -1; % advective-limited +out.limiter_cell(dif_better) = +1; % diffusive-limited +end diff --git a/Utilities/Misc/geometry/colocated_schemes/do_projection_stage.m b/Utilities/Misc/geometry/colocated_schemes/do_projection_stage.m new file mode 100644 index 00000000000..0d4f65e6150 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/do_projection_stage.m @@ -0,0 +1,59 @@ +function [u_out, Uface_out, H_cell, p_out, dHdn_f, gradH] = ... + do_projection_stage(u_in, u_in1, Uface_in, Uface_in1, p_in, rho, mu, lambda, ... + cells, faces, Lx, Ly, dt, rho0, gvec, divU_th, ... + A, poisson_opts, add_meanK, stage) + +global H_projection + +% ---------- build viscous ---------- +[~, FA_visc_cell, FA_visc_face_n,~] = ... + build_viscous_forces_2d(u_in1, rho, mu, lambda, cells, faces, Lx, Ly); + +% ---------- advective + body (uses GIVEN face U at this stage) ---------- +[FAdv_cell, FAdv_face_n] = ... + build_adv_body_forces_2d(u_in1, rho, rho0, gvec, Uface_in, cells, faces, stage); + +% ---------- baroclinic/B-term in cells then faces ---------- +[FB_cell, FB_face_n] = build_FB_cell_face(p_in, rho, cells, faces, Lx, Ly); + +% Combine face & cell contributions for RHS and cell update +FA_face_n = FAdv_face_n + FA_visc_face_n; % faces +FA_cell = FAdv_cell + FA_visc_cell; % cells + +% ---------- Poisson RHS (flux form) ---------- +rhs = build_poisson_rhs_H(Uface_in, divU_th, dt, cells, faces, FA_face_n, FB_face_n); + +% ---------- Solve Poisson for H (your sign convention) ---------- +[H_cell, dHdn_f, gradH] = poisson_solve(A, rhs, cells, faces, Lx, Ly, poisson_opts); + +% ---------- (optional) gauge tweak to help p match analytic TG) ---------- +if(H_projection) + if add_meanK + K_cell = 0.5*sum(u_in.^2,2); + H_gauge = H_cell + mean(K_cell); + else + H_gauge = H_cell; + end +end + +% ---------- Project faces (mass consistency) ---------- +Uface_out = project_faces(Uface_in, Uface_in1, dt, FA_face_n, FB_face_n, dHdn_f, stage); + +% ---------- Project cells (momentum update) ---------- +if(stage==1) + u_out = reconstruct_u_from_faces_w(Uface_out, cells, faces); +% u_out2 = project_cells(u_in, dt, FA_cell, divTau_cell, rho, FB_cell, gradH); +% u_out = 0.5*(u_out1 + u_out2); +else + u_out = project_cells(u_in, u_in1, dt, FA_cell, FB_cell, gradH, stage); +end + +if(H_projection) + % ---------- Pressure from Bernoulli for NEXT stage inputs ---------- + K_out = 0.5*sum(u_out.^2,2); + p_out = rho .* (H_gauge - K_out); +else + p_out = H_cell; +end + +end diff --git a/Utilities/Misc/geometry/colocated_schemes/find_cell_by_center_centered.m b/Utilities/Misc/geometry/colocated_schemes/find_cell_by_center_centered.m new file mode 100644 index 00000000000..333dd54b09b --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/find_cell_by_center_centered.m @@ -0,0 +1,24 @@ +function p = find_cell_by_center_centered(cells, Nx, Ny, Lx, Ly, i, j) +% Find the cell index p whose centroid is closest to logical indices (i,j) +% on a uniform grid spanning [-Lx/2, Lx/2] × [-Ly/2, Ly/2] with periodic wrap. + +dx = Lx / Nx; +dy = Ly / Ny; + +% target cell-center coordinates for (i,j) +xt = -Lx/2 + (i - 0.5) * dx; +yt = -Ly/2 + (j - 0.5) * dy; + +% collect centroids +Nc = numel(cells); +Xc = zeros(Nc,2); +for k = 1:Nc + Xc(k,:) = cells(k).xc(:).'; +end + +% minimal-image periodic distances to target (works for centered domain) +dxv = Xc(:,1) - xt; dxv = dxv - round(dxv / Lx) * Lx; +dyv = Xc(:,2) - yt; dyv = dyv - round(dyv / Ly) * Ly; + +[~, p] = min(dxv.^2 + dyv.^2); +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/grad_from_compact_normals_w.m b/Utilities/Misc/geometry/colocated_schemes/grad_from_compact_normals_w.m new file mode 100644 index 00000000000..a15723908ff --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/grad_from_compact_normals_w.m @@ -0,0 +1,57 @@ +function [dphidn_f, gradphi_cell] = grad_from_compact_normals_w(phi, cells, faces, Lx, Ly, wtype, beta) +if nargin<6 || isempty(wtype), wtype = 'trans'; end % 'area' | 'trans' | 'diamond' | 'blend' +if nargin<7 || isempty(beta), beta = 0.3; end + +Nf = numel(faces); Nc = numel(cells); +dphidn_f = zeros(Nf,1); + + +% face-normal derivatives from two-point differences +for f = 1:Nf + P = faces(f).owner; N = faces(f).neigh; nf = faces(f).nf(:); + + if N>0 + dr_raw = cells(N).xc(:) - cells(P).xc(:); + dr = periodic_delta_2d(dr_raw, Lx, Ly); + dphidn_f(f) = (phi(N) - phi(P)) / dot(dr,nf); + else + disp('N==0') + dphidn_f(f) = 0.0; % set per BC if needed + end +end + +% LS gradient per cell with weights +gradphi_cell = zeros(Nc,2); +for p = 1:Nc + M = zeros(2,2); b = zeros(2,1); + flist = cells(p).faces(:); + for k = 1:numel(flist) + f = flist(k); + nf = faces(f).nf(:); + Af = faces(f).Af; + % outward normal for cell p + s_pf = +1; if faces(f).neigh==p, s_pf = -1; end + nout = s_pf*nf; + + % choose weight + dr_raw = cells(N).xc(:) - cells(P).xc(:); + dr = periodic_delta_2d(dr_raw, Lx, Ly); + dPN = abs(dot(dr,nf)); + switch lower(wtype) + case 'area' , wf = Af; + case 'trans' , wf = Af/max(dPN,eps); + case 'diamond' , wf = Af*dPN; + case 'blend' , wf = (1-beta)*Af + beta*Af/max(dPN,eps); + otherwise , wf = Af; + end + + qf = dphidn_f(f); % owner->neigh normal derivative + % Use owner->neigh normal in RHS (no extra sign): + M = M + wf * (nout * nout.'); + b = b + wf * qf * nf; + end + %trM = trace(M); if trM<=0, trM=1; end + %M = M + (1e-12*trM)*eye(2); + gradphi_cell(p,:) = (M\b).'; +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/l2_error_ke.m b/Utilities/Misc/geometry/colocated_schemes/l2_error_ke.m new file mode 100644 index 00000000000..f313e475cc5 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/l2_error_ke.m @@ -0,0 +1,22 @@ +function [L2_abs, L2_rel, K_num_tot, K_ex_tot, K_ex_formula] = l2_error_ke(u_num, rho, cells, t, tg) + +global ux0 uy0 + +Nc = numel(cells); +xc = zeros(Nc,1); yc = zeros(Nc,1); V = zeros(Nc,1); +for p = 1:Nc, xc(p)=cells(p).xc(1); yc(p)=cells(p).xc(2); V(p)=cells(p).V; end +u_ex = taylor_green_uv_2d(xc, yc, t, tg.U0, tg.kx, tg.ky, tg.nu); +k_num = 0.5*rho(:).*sum(u_num.^2,2); +k_ex = 0.5*rho(:).*sum(u_ex.^2,2); +num_sq = V.*(k_num-k_ex).^2; +den_sq = V.*(k_ex).^2; +L2_abs = sqrt(sum(num_sq)/sum(V)); +L2_rel = sqrt(sum(num_sq)/max(1e-300,sum(den_sq))); +K_num_tot = sum(V.*k_num); +K_ex_tot = sum(V.*k_ex); +% Correct analytic total KE (mean k = ρ U0^2 / 4 e^{-2 ν k^2 t}) +k2 = tg.kx*tg.kx + tg.ky*tg.ky; +Vol = sum(V); +rho_bar = sum(V.*rho(:))/Vol; +K_ex_formula = rho_bar*(tg.U0^2/4)*exp(-2*tg.nu*k2*t)*Vol+(ux0^2+uy0^2)/2*Vol; +end diff --git a/Utilities/Misc/geometry/colocated_schemes/l2_error_momentum.m b/Utilities/Misc/geometry/colocated_schemes/l2_error_momentum.m new file mode 100644 index 00000000000..da413baaf4b --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/l2_error_momentum.m @@ -0,0 +1,14 @@ +function [L2_abs, L2_rel, M_num, M_ex] = l2_error_momentum(u_num, rho, cells, t, tg) +% tg has fields: U0, kx, ky, nu +Nc = numel(cells); +xc = zeros(Nc,1); yc = zeros(Nc,1); V = zeros(Nc,1); +for p = 1:Nc, xc(p)=cells(p).xc(1); yc(p)=cells(p).xc(2); V(p)=cells(p).V; end +u_ex = taylor_green_uv_2d(xc, yc, t, tg.U0, tg.kx, tg.ky, tg.nu); +m_num = rho(:).*u_num; m_ex = rho(:).*u_ex; +num_sq = V.*sum((m_num-m_ex).^2,2); +den_sq = V.*sum(m_ex.^2,2); +L2_abs = sqrt(sum(num_sq)/sum(V)); +L2_rel = sqrt(sum(num_sq)/max(1e-300,sum(den_sq))); +M_num = [sum(V.*m_num(:,1)), sum(V.*m_num(:,2))]; +M_ex = [sum(V.*m_ex(:,1)), sum(V.*m_ex(:,2)) ]; +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/ls_grad_scalar_2d.m b/Utilities/Misc/geometry/colocated_schemes/ls_grad_scalar_2d.m new file mode 100644 index 00000000000..58c59bba4be --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/ls_grad_scalar_2d.m @@ -0,0 +1,25 @@ +function gradPhi = ls_grad_scalar_2d(phi, cells, faces, Lx, Ly) +% phi: Nc x 1 cell scalar +Nc = numel(cells); +gradPhi = zeros(Nc,2); +for p = 1:Nc + xp = cells(p).xc; + fids = cells(p).faces; + M = zeros(2,2); b = [0;0]; + for k = 1:numel(fids) + f = fids(k); + P = faces(f).owner; N = faces(f).neigh; + if P==p && N>0, q = N; + elseif N==p && P>0, q = P; + else, continue; end % (no boundary faces in periodic mesh) + rq = periodic_delta_2d(cells(q).xc - xp, Lx, Ly); + Af = faces(f).Af; + w = Af / (dot(rq,rq) + 1e-14); % robust weight + M = M + w * (rq(:)*rq(:).'); + b = b + w * (phi(q) - phi(p)) * rq(:); + end + trM = trace(M); if trM<=0, trM = 1; end + M = M + (1e-12*trM)*eye(2); % tiny regularization + gradPhi(p,:) = (M\b).'; +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/ls_grad_velocity_2d.m b/Utilities/Misc/geometry/colocated_schemes/ls_grad_velocity_2d.m new file mode 100644 index 00000000000..1e99e0a78a7 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/ls_grad_velocity_2d.m @@ -0,0 +1,34 @@ +function [gradU, gradV] = ls_grad_velocity_2d(u, cells, faces, Lx, Ly) +Nc = size(u,1); +gradU = zeros(Nc,2); gradV = zeros(Nc,2); + +for p = 1:Nc + xp = cells(p).xc; + fids = cells(p).faces; + M = zeros(2,2); bu = [0;0]; bv = [0;0]; + + for k = 1:numel(fids) + f = fids(k); + P = faces(f).owner; N = faces(f).neigh; + if P==p && N>0, q = N; + elseif N==p && P>0, q = P; + else, continue + end + rq_raw = cells(q).xc - xp; + rq = periodic_delta_2d(rq_raw, Lx, Ly); % <-- wrap here + Af = faces(f).Af; + w = Af / dot(rq,rq); + + M = M + w * (rq(:)*rq(:).'); + du = u(q,1) - u(p,1); + dv = u(q,2) - u(p,2); + bu = bu + w * du * rq(:); + bv = bv + w * dv * rq(:); + end + + trM = trace(M); if trM<=0, trM = 1; end + M = M + (1e-12*trM)*eye(2); + gradU(p,:) = (M\bu).'; + gradV(p,:) = (M\bv).'; +end +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/periodic_delta_2d.m b/Utilities/Misc/geometry/colocated_schemes/periodic_delta_2d.m new file mode 100644 index 00000000000..30b981a1d14 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/periodic_delta_2d.m @@ -0,0 +1,10 @@ +function dr = periodic_delta_2d(dx_raw, Lx, Ly) +% Wrap dx_raw = [dx, dy] into (-L/2, L/2] +dr = dx_raw; +% x +if dr(1) > 0.5*Lx, dr(1) = dr(1) - Lx; +elseif dr(1) <= -0.5*Lx, dr(1) = dr(1) + Lx; end +% y +if dr(2) > 0.5*Ly, dr(2) = dr(2) - Ly; +elseif dr(2) <= -0.5*Ly, dr(2) = dr(2) + Ly; end +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/plot_L2_errors.m b/Utilities/Misc/geometry/colocated_schemes/plot_L2_errors.m new file mode 100644 index 00000000000..a2f448b4182 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/plot_L2_errors.m @@ -0,0 +1,27 @@ +function plot_L2_errors(hist) +%PLOT_L2_ERRORS Plot L2 absolute errors for momentum and kinetic energy. +% Expects fields: hist.t, hist.L2_m_abs, hist.L2_k_abs (vectors) + +t = hist.t(:); +eMm = hist.L2_m_abs(:); +eKk = hist.L2_k_abs(:); + +% Guard against zeros for log scale +eMm_plot = max(eMm, eps); +eKk_plot = max(eKk, eps); + +figure('Name','L2 Errors','Color','w'); + +subplot(1,2,1) +semilogy(t, eMm_plot, '-o','LineWidth',1.8,'MarkerSize',5); grid on +xlabel('t','FontSize',16); ylabel('L2(M)','FontSize',16) +title('L2 momentum error','FontSize',18) +set(gca,'FontSize',14) + +subplot(1,2,2) +semilogy(t, eKk_plot, '-s','LineWidth',1.8,'MarkerSize',5); grid on +xlabel('t','FontSize',16); ylabel('L2(KE)','FontSize',16) +title('L2 kinetic energy error','FontSize',18) +set(gca,'FontSize',14) + +end diff --git a/Utilities/Misc/geometry/colocated_schemes/plot_TG_field.m b/Utilities/Misc/geometry/colocated_schemes/plot_TG_field.m new file mode 100644 index 00000000000..d5952b84271 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/plot_TG_field.m @@ -0,0 +1,18 @@ +function [fid]=plot_TG_field(t,Nx,Ny,Xc,Yc,u,p) + +Ux = reshape(u(:,1), [Nx,Ny])'; % Ny x Nx +Uy = reshape(u(:,2), [Nx,Ny])'; +Pmat = reshape(p, [Nx,Ny])'; + +fid=0; clf + +nexttile; contourf(Xc, Yc, Ux, 20, 'LineStyle','none'); axis equal tight; +colorbar; title(['u_x t=' num2str(t)]); xlabel('x'); ylabel('y'); + +nexttile; contourf(Xc, Yc, Uy, 20, 'LineStyle','none'); axis equal tight; +colorbar; title('u_y'); xlabel('x'); ylabel('y'); + +nexttile; contourf(Xc, Yc, Pmat, 20, 'LineStyle','none'); axis equal tight; +colorbar; title('pressure p'); xlabel('x'); ylabel('y'); + +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/plot_mass_mom_ke.m b/Utilities/Misc/geometry/colocated_schemes/plot_mass_mom_ke.m new file mode 100644 index 00000000000..16fef56d5cc --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/plot_mass_mom_ke.m @@ -0,0 +1,80 @@ +function plot_mass_mom_ke(hist, cells) +%PLOT_MASS_MOM_KE Plot mass residual, , and . +% Inputs: +% hist : struct from run_ssprk2 (fields used: t, L2_mass, M_num, M_ex, +% K_num_tot, K_ex_formula) +% cells : mesh cells (used to compute domain volume) + +Vol = sum([cells.V]); + +% Averages +Mx_num = hist.M_num(:,1) / Vol; +My_num = hist.M_num(:,2) / Vol; +Mmag_num = sqrt(Mx_num.^2 + My_num.^2); + +hasAnalyticM = isfield(hist,'M_ex') && all(size(hist.M_ex)==size(hist.M_num)) ... + && all(~isnan(hist.M_ex(:))); +if hasAnalyticM + Mx_ex = hist.M_ex(:,1) / Vol; My_ex = hist.M_ex(:,2) / Vol; + Mmag_ex = sqrt(Mx_ex.^2 + My_ex.^2); +else + Mx_ex = []; My_ex = []; Mmag_ex = []; +end + +KE_num_avg = hist.K_num_tot / Vol; +hasAnalyticK = isfield(hist,'K_ex_formula') && all(~isnan(hist.K_ex_formula(:))); +if hasAnalyticK + KE_ex_avg = hist.K_ex_tot / Vol; %hist.K_ex_formula +else + KE_ex_avg = []; +end + +% ---- Plots ---- +figure('Name','Mass/Momentum/KE','Color','w','Position',[100 100 1200 400]); + +fs_axis = 12; % axis/tick font size +fs_label = 14; % x/y label size +fs_title = 16; % title size +fs_leg = 14; % legend size + +% (1) Mass conservation residual +subplot(1,3,1) +semilogy(hist.t, hist.L2_mass, '-o','LineWidth',1.5,'MarkerSize',6); grid on +xlabel('t','FontSize',fs_label,'FontWeight','bold'); +ylabel('||mass||_2','FontSize',fs_label,'FontWeight','bold'); +title('Mass conservation residual','FontSize',fs_title,'FontWeight','bold'); +set(gca,'FontSize',fs_axis); + +% (2) Average momentum +subplot(1,3,2) +plot(hist.t, Mmag_num, '-o', 'LineWidth',1.5,'MarkerSize',6,'DisplayName','|| num'); hold on +plot(hist.t, Mx_num, '-s', 'LineWidth',1.5,'MarkerSize',6,'DisplayName',' num'); +plot(hist.t, My_num, '-d', 'LineWidth',1.5,'MarkerSize',6,'DisplayName',' num'); +if ~isempty(Mmag_ex) + plot(hist.t, Mmag_ex, '--','LineWidth',1.5,'DisplayName','|| ex'); + plot(hist.t, Mx_ex, '--','LineWidth',1.5,'DisplayName',' ex'); + plot(hist.t, My_ex, '--','LineWidth',1.5,'DisplayName',' ex'); +end +hold off; grid on +xlabel('t','FontSize',fs_label,'FontWeight','bold'); +ylabel('Momentum / Volume','FontSize',fs_label,'FontWeight','bold'); +title('Domain-averaged momentum','FontSize',fs_title,'FontWeight','bold'); +legend('Location','best','FontSize',fs_leg); +set(gca,'FontSize',fs_axis); + +% (3) Average kinetic energy +subplot(1,3,3) +plot(hist.t, KE_num_avg, '-o','LineWidth',1.5,'MarkerSize',6,'DisplayName',' num'); hold on +if ~isempty(KE_ex_avg) + plot(hist.t, KE_ex_avg, '--','LineWidth',1.5,'DisplayName',' ex'); +end +hold off; grid on +xlabel('t','FontSize',fs_label,'FontWeight','bold'); +ylabel('KE / Volume','FontSize',fs_label,'FontWeight','bold'); +title('Domain-averaged kinetic energy','FontSize',fs_title,'FontWeight','bold'); +legend('Location','best','FontSize',fs_leg); +set(gca,'FontSize',fs_axis); + +sgtitle('Diagnostics: Volume-averaged quantities','FontSize',fs_title+2,'FontWeight','bold'); + +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/poisson_solve.m b/Utilities/Misc/geometry/colocated_schemes/poisson_solve.m new file mode 100644 index 00000000000..8bd50d46944 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/poisson_solve.m @@ -0,0 +1,177 @@ +function [H_cell, dHdn_f, gradH, info] = poisson_solve(A, rhs, cells, faces, Lx, Ly, opts) +%POISSON_SOLVE Solve compact Poisson for H and return zero-mean H, face normals and LS gradients. +% +% Inputs +% A : Nc x Nc sparse Poisson matrix (compact Laplacian) +% rhs : Nc x 1 right-hand-side (flux form) +% cells : struct array with fields .V (scalar), .xc(1x2), .faces(1:nf) +% faces : struct array with fields .owner, .neigh, .nf(1x2 unit owner->neigh), .Af +% Lx,Ly : domain lengths (for periodic minimum-image deltas) +% opts : struct with optional fields: +% .fix = 'pin' (default) or 'mean' % nullspace fix +% .pin_idx = 1 (default) % index to pin if fix='pin' +% .tol = 1e-12 % PCG tolerance +% .maxit = 500 % PCG max iterations +% .use_ilu = true % try ILU preconditioner +% .ilu_droptol = 1e-6 % ILU drop tolerance +% +% Outputs +% H_cell : Nc x 1 solution with volume-weighted zero mean +% dHdn_f : Nf x 1 compact face-normal derivative ( (H_N-H_P)/((xN-xP)·nf) ) +% gradH : Nc x 2 LS reconstruction of ∇H from dHdn_f +% info : struct with .flag, .relres, .iter, .time_solve +% +% Notes +% - Domain assumed fully periodic (pure Neumann); we remove the nullspace. +% - For non-periodic cases, extend the matrix/RHS assembly accordingly. +global compact_laplacian + +% Flag to set the gradient in cell centers either to divergence gradient (1) +% or Least Squares(0): +do_divergence_gradient = 1; + + +Nc = numel(cells); Nf = numel(faces); + +% ----- defaults +if nargin < 7, opts = struct(); end +opts = set_default(opts, 'fix', 'pin'); +opts = set_default(opts, 'pin_idx', 1); +opts = set_default(opts, 'tol', 1e-12); +opts = set_default(opts, 'maxit', 5000); +opts = set_default(opts, 'use_ilu', true); +opts = set_default(opts, 'ilu_droptol', 1e-6); + +% ----- nullspace fix and change signs to make eigs of A > 0: +switch lower(opts.fix) + case 'pin' + [Afix, rhsfix] = pin_one_dof(-A, -rhs, opts.pin_idx); + take = 1:Nc; + case 'mean' + [Afix, rhsfix] = mean_zero_constraint(-A, -rhs, cells); + take = 1:Nc; % first Nc entries are H; last one is Lagrange multiplier + otherwise + error('opts.fix must be ''pin'' or ''mean''.'); +end + +% ----- solve +t0 = tic; +L = []; U = []; +if opts.use_ilu + try + setup.type = 'ilutp'; + setup.droptol = opts.ilu_droptol; + [L,U] = ilu(Afix, setup); + catch + L = []; U = []; + end +end + +try + [Hsol, flag, relres, iter] = pcg(Afix, rhsfix, opts.tol, opts.maxit, L, U); +catch + [Hsol, flag, relres, iter] = pcg(Afix, rhsfix, opts.tol, opts.maxit); +end +info.flag = flag; info.relres = relres; info.iter = iter; info.time_solve = toc(t0); +if flag ~= 0 + warning('poisson_solve: PCG did not fully converge (flag=%d, relres=%.2e, iter=%d).', flag, relres, iter); +end + +% Extract H and enforce zero-mean +H_cell = Hsol(take); +V = reshape([cells.V], [], 1); +Hbar = sum(V .* H_cell) / sum(V); +H_cell = H_cell - Hbar; + +if(compact_laplacian) + +[dHdn_f, gradH] = grad_from_compact_normals_w(H_cell, cells, faces, Lx, Ly,'trans'); +%[dHdn_f, gradH] = grad_from_compact_normals(H_cell, cells, faces, Lx, Ly); +% ----- face-normal compact derivative dH/dn +% dHdn_f = zeros(Nf,1); +% for f = 1:Nf +% P = faces(f).owner; +% N = faces(f).neigh; +% nf = faces(f).nf(:); +% dvec = periodic_delta_2d(cells(N).xc - cells(P).xc, Lx, Ly); +% dn = dot(dvec, nf); +% dHdn_f(f) = (H_cell(N) - H_cell(P)) / dn; +% end +% +% % ----- ∇H at cells ----- +% gradH = zeros(Nc,2); +% +% if do_divergence_gradient == 1 +% % === Green–Gauss (divergence theorem) using H_f = 0.5(H_P+H_N) === +% sum_flux = zeros(Nc,2); % Σ H_f n_f A_f +% for f = 1:Nf +% P = faces(f).owner; +% N = faces(f).neigh; +% nf = faces(f).nf(:); +% Af = faces(f).Af; +% +% % face-averaged H +% if N > 0 +% Hf = 0.5*(H_cell(P) + H_cell(N)); +% else +% % (boundary case; replace with BC value if needed) +% Hf = H_cell(P); +% end +% +% contrib = (Hf * nf.').*Af; % 1x2 vector +% +% % owner (+), neighbor (−) +% sum_flux(P,:) = sum_flux(P,:) + contrib; +% if N > 0 +% sum_flux(N,:) = sum_flux(N,:) - contrib; +% end +% end +% gradH = sum_flux ./ V; % V is Nc×1 of cell volumes +% +% else +% % === Least–Squares from compact normals: M_p ∇H_p = Σ Af (∂H/∂n)_f n_f === +% for p = 1:Nc +% M = zeros(2,2); b = [0;0]; +% fids = cells(p).faces; +% for k = 1:numel(fids) +% f = fids(k); +% nf = faces(f).nf(:); +% Af = faces(f).Af; +% qf = dHdn_f(f); % (H_N - H_P)/dn +% M = M + Af * (nf*nf.'); +% b = b + Af * qf * nf; +% end +% trM = trace(M); if trM<=0, trM=1; end +% M = M + (1e-12*trM)*eye(2); % tiny Tikhonov for robustness +% gradH(p,:) = (M\b).'; +% end +% end + +else + + [gradH, ~, dHdn_f] = gradH_sparse_GG(H, cells, faces); + +end + +end + +% ======= helpers ======= + +function [Afix, rhsfix] = pin_one_dof(A, rhs, pind) +Afix = A; rhsfix = rhs; +Afix(pind,:) = 0; Afix(:,pind) = 0; +Afix(pind,pind) = 1; +rhsfix(pind) = 0; +end + +function [Aaug, rhsaug] = mean_zero_constraint(A, rhs, cells) +Nc = numel(cells); +V = reshape([cells.V], [], 1); +W = V / sum(V); % volume-weighted mean +Aaug = [A, W; W.', 0]; +rhsaug= [rhs; 0]; +end + +function s = set_default(s, name, val) +if ~isfield(s,name) || isempty(s.(name)), s.(name) = val; end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/project_cells.m b/Utilities/Misc/geometry/colocated_schemes/project_cells.m new file mode 100644 index 00000000000..c77a174fa5c --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/project_cells.m @@ -0,0 +1,13 @@ +function u_new = project_cells(u, u1, dt, FA_cell, FB_cell, gradH, stage) +% FA_cell : Nc x 2 (advective + body) cell vector +% divTau_cell : Nc x 2 (∇·τ) at cells +% rho : Nc x 1 +% FB_cell : Nc x 2 (-p ∇(1/ρ)) at cells +% gradH : Nc x 2 ∇H at cells +% Output: +% u_new : Nc x 2 +if(stage==1) + u_new = u - dt*( FA_cell + FB_cell + gradH ); +else + u_new = 0.5*(u+u1) - dt/2*( FA_cell + FB_cell + gradH ); +end diff --git a/Utilities/Misc/geometry/colocated_schemes/project_faces.m b/Utilities/Misc/geometry/colocated_schemes/project_faces.m new file mode 100644 index 00000000000..9c34c2e7e02 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/project_faces.m @@ -0,0 +1,10 @@ +function Uface_new = project_faces(Uface, Uface1, dt, FA_face_n, FB_face_n, dHdn_f, stage) +% Inputs: all Nf x 1, owner→neigh sign convention +% Output: Uface_new (Nf x 1) +if(stage==1) + % Predictor: + Uface_new = Uface - dt*( FA_face_n + FB_face_n + dHdn_f ); +else + % Corrector: + Uface_new = 0.5*(Uface+Uface1) - dt/2*( FA_face_n + FB_face_n + dHdn_f ); +end diff --git a/Utilities/Misc/geometry/colocated_schemes/reconstruct_u_from_faces_w.m b/Utilities/Misc/geometry/colocated_schemes/reconstruct_u_from_faces_w.m new file mode 100644 index 00000000000..6395e8ea294 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/reconstruct_u_from_faces_w.m @@ -0,0 +1,48 @@ +function u_cell = reconstruct_u_from_faces_w(Uface, cells, faces, wtype, beta) +if nargin<4, wtype='trans'; end +if nargin<5, beta=0.3; end + +Nc = numel(cells); dim = numel(cells(1).xc); +u_cell = zeros(Nc, dim); + +for p = 1:Nc + M = zeros(dim,dim); b = zeros(dim,1); + flist = cells(p).faces(:); + for k = 1:numel(flist) + f = flist(k); + nf = faces(f).nf(:); % owner→neighbor + Af = faces(f).Af; + + % outward normal for this cell: + s_pf = +1; if faces(f).neigh == p, s_pf = -1; end + nout = s_pf * nf; + + % choose weight + if isfield(faces,'dPN') && ~isempty(faces(f).dPN) + dPN = faces(f).dPN; + else + % fallback estimate (only used for weighting) + oth = faces(f).owner + faces(f).neigh - p; + dvec = cells(oth).xc(:) - cells(p).xc(:); + dPN = abs(dvec.'*nout); + end + switch lower(wtype) + case 'area' , wf = Af; + case 'trans' , wf = Af/max(dPN,eps); + case 'diamond' , wf = Af*dPN; + case 'blend' , wf = (1-beta)*Af + beta*Af/max(dPN,eps); + otherwise , wf = Af; + end + + % --- Correct RHS build (key fix) --- + % EITHER: b = b + wf * (s_pf*Uface(f)) * nout; + % OR (equivalent, safer): + b = b + wf * Uface(f) * nf; % uses nf, not nout + + M = M + wf * (nout * nout.'); + end + trM = trace(M); if trM<=0, trM=1; end + M = M + (1e-12*trM)*eye(dim); + u_cell(p,:) = (M\b).'; +end +end diff --git a/Utilities/Misc/geometry/colocated_schemes/run_ssprk2.m b/Utilities/Misc/geometry/colocated_schemes/run_ssprk2.m new file mode 100644 index 00000000000..e7ec3fad783 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/run_ssprk2.m @@ -0,0 +1,198 @@ +function [u, Uface, p, H_cell, t, hist] = run_ssprk2( ... + u, Uface, p, rho, mu, lambda, ... + cells, faces, Lx, Ly, params) +global compact_laplacian Nx Ny Xc Yc dstep_plot do_ssprk2 H_projection +% ----------------- defaults ----------------- +Nc = numel(cells); +defaults = struct( ... + 'rho0', 1.0, ... + 'gvec', [0,0], ... + 'CFL', 0.8, ... + 'VN', 0.5, ... + 'divU_th', zeros(Nc,1), ... + 't_begin', 0.0, ... + 't_end', 1.0, ... + 'dt', 0.01, ... + 'variable_time_step', 1, ... + 'max_steps', 10^6, ... + 'add_meanK', true, ... + 'poisson_opts', struct('fix','pin','pin_idx',1,'use_ilu',true), ... + 'tg_params', [] ... % struct with fields U0,kx,ky,nu (optional) + ); +params = set_defaults(params, defaults); + +% Prebuild compact Poisson matrix +A = build_poisson_matrix_compact(cells, faces, Lx, Ly); + +% History +hist = struct('step', [], 't',[], 'dt',[], 'L2_mass',[], ... + 'L2_m_abs',[], 'L2_m_rel',[], 'M_num',[], 'M_ex',[], ... + 'L2_k_abs',[], 'L2_k_rel',[], 'K_num_tot',[], ... + 'K_ex_tot',[], 'K_ex_formula',[], 'u_num_mid', [], ... + 'u_ex_mid', []); + +% Initial time parameters: +t = params.t_begin; +dt= params.dt; +step = 0; + +% Initial Mass, Momentum and TKE: +divU_after = cell_divergence_from_faces(Uface, cells, faces); +res_mass = divU_after - params.divU_th; +L2_mass = norm(res_mass,2); +tg = params.tg_params; % struct U0,kx,ky,nu +% Momentum +[L2m_abs, L2m_rel, M_num, M_ex] = l2_error_momentum(u, rho, cells, t, tg); +% KE +[L2k_abs, L2k_rel, K_num_tot, K_ex_tot, K_ex_formula] = ... + l2_error_ke(u, rho, cells, t, tg); + +% ---- Save history ---- +hist.step(end+1,1) = step; +hist.t(end+1,1) = t; +hist.dt(end+1,1) = dt; +hist.L2_mass(end+1,1) = L2_mass; + +hist.L2_m_abs(end+1,1)= L2m_abs; +hist.L2_m_rel(end+1,1)= L2m_rel; +hist.M_num(end+1,1:2) = M_num; +hist.M_ex(end+1,1:2) = M_ex; + +hist.L2_k_abs(end+1,1)= L2k_abs; +hist.L2_k_rel(end+1,1)= L2k_rel; +hist.K_num_tot(end+1,1) = K_num_tot; +hist.K_ex_tot(end+1,1) = K_ex_tot; +hist.K_ex_formula(end+1,1)= K_ex_formula; + +hist.u_num_mid(end+1,1:2)=u(params.p_mid,1:2); + +% Optional short log +fprintf(['SSPRK2 %5d t=%.6e dt=%.3e ||mass||_2=%.3e ', ... + 'L2m=%.2e L2k=%.2e Knum=%.3e Kex=%.3e\n'], ... + step, t, dt, L2_mass, L2m_abs, L2k_abs, K_num_tot/(Lx*Ly), K_ex_tot/(Lx*Ly)); + +[fid]=plot_TG_field(t,Nx,Ny,Xc,Yc,u,p); +figure(gcf) +pause(0.2) + +% Time loop +while (t < params.t_end) && (step < params.max_steps) + step = step + 1; + + % ---- dt from CFL/VN ---- + if (params.variable_time_step==1) + dts = compute_dt_explicit(cells, faces, Uface, mu, rho, params.CFL, params.VN, Lx, Ly); + else + dts.dt = dt; + end + dt = min(dts.dt, params.t_end - t); + + % ---- Stage 1 ---- + [u1, Uface1, H1, p1] = do_projection_stage(u,u,Uface,Uface,p, rho, mu, lambda, ... + cells, faces, Lx, Ly, dt, ... + params.rho0, params.gvec, params.divU_th, ... + A, params.poisson_opts, params.add_meanK,1); + + if(do_ssprk2==1) % SSPRK2 + % ---- Stage 2 ---- + [u2, Uface2, H2, p2] = do_projection_stage(u,u1, ... + Uface,Uface1,p1, rho, mu, lambda, ... + cells, faces, Lx, Ly, dt, ... + params.rho0, params.gvec, params.divU_th, ... + A, params.poisson_opts, params.add_meanK,2); + else % Forward Euler. + u2=u1; H2=H1; Uface2=Uface1; + end + + % ---- SSPRK2 combo ---- + u_next = u2; + Uface_next = Uface2; + H_cell = H2; + + % Pressure from Bernoulli (with optional gauge tweak) + if(H_projection) + K_next = 0.5*sum(u_next.^2,2); + if params.add_meanK + H_eff = H_cell + mean(K_next); + else + H_eff = H_cell; + end + p_next = rho .* (H_eff - K_next); + else + H_eff = H_cell; + p_next = H_cell; + end + + % ---- Mass residual ---- + divU_after = cell_divergence_from_faces(Uface_next, cells, faces); + res_mass = divU_after - params.divU_th; + L2_mass = norm(res_mass,2); + + % ---- Momentum / KE checks (if TG params available) ---- + if ~isempty(params.tg_params) + tg = params.tg_params; % struct U0,kx,ky,nu + + % Momentum + [L2m_abs, L2m_rel, M_num, M_ex] = l2_error_momentum(u_next, rho, cells, t+dt, tg); + + % KE + [L2k_abs, L2k_rel, K_num_tot, K_ex_tot, K_ex_formula] = ... + l2_error_ke(u_next, rho, cells, t+dt, tg); + else + L2m_abs=NaN; L2m_rel=NaN; M_num=[NaN NaN]; M_ex=[NaN NaN]; + L2k_abs=NaN; L2k_rel=NaN; K_num_tot=NaN; K_ex_tot=NaN; K_ex_formula=NaN; + end + + % ---- Commit & time advance ---- + u = u_next; + Uface = Uface_next; + p = p_next; + t = t + dt; + + % ---- check pressure work ------ + % check_product_rule_cellwise(u, H_eff, cells, faces) + % check_H_divU(H_cell, Uface, cells, faces) + + % ---- Save history ---- + hist.step(end+1,1) = step; + hist.t(end+1,1) = t; + hist.dt(end+1,1) = dt; + hist.L2_mass(end+1,1) = L2_mass; + + hist.L2_m_abs(end+1,1)= L2m_abs; + hist.L2_m_rel(end+1,1)= L2m_rel; + hist.M_num(end+1,1:2) = M_num; + hist.M_ex(end+1,1:2) = M_ex; + + hist.L2_k_abs(end+1,1)= L2k_abs; + hist.L2_k_rel(end+1,1)= L2k_rel; + hist.K_num_tot(end+1,1) = K_num_tot; + hist.K_ex_tot(end+1,1) = K_ex_tot; + hist.K_ex_formula(end+1,1)= K_ex_formula; + + hist.u_num_mid(end+1,1:2)=u(params.p_mid,1:2); + + % Optional short log + fprintf(['SSPRK2 %5d t=%.6e dt=%.3e ||mass||_2=%.3e ', ... + 'L2m=%.2e L2k=%.2e Knum=%.3e Kex=%.3e\n'], ... + step, t, dt, L2_mass, L2m_abs, L2k_abs, K_num_tot/(Lx*Ly), K_ex_tot/(Lx*Ly)); + + % Plot tg field: + conditional = mod(step,dstep_plot)<1e-3; + if mod(step,dstep_plot)<1e-3 + [fid]=plot_TG_field(t,Nx,Ny,Xc,Yc,u,p); + figure(gcf) + pause(0.2) + end +end +end + +function s = set_defaults(s, d) +fn = fieldnames(d); +for k=1:numel(fn) + if ~isfield(s,fn{k}) || isempty(s.(fn{k})) + s.(fn{k}) = d.(fn{k}); + end +end +end + diff --git a/Utilities/Misc/geometry/colocated_schemes/taylor_green.m b/Utilities/Misc/geometry/colocated_schemes/taylor_green.m new file mode 100644 index 00000000000..c7b300e385b --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/taylor_green.m @@ -0,0 +1,137 @@ +function [params,hist,cells,faces]=taylor_green(Nx_in,Ny_in,nu,dstep_plot_in) + +% Nx; Ny % cells in x,y +% nu % kinematic viscosity +% Rest of parameters as in FDS Verification Guide Taylor-Green vortex +% problem (ns2d). +% Taylor green Vortex Problem, colocated grid, non-conservative Low Mach +% approximation: +% -------------------------------------------------------------- +global compact_laplacian ux0 uy0 Nx Ny Xc Yc dstep_plot do_ssprk2 ... + H_projection +%close all +%clc + +compact_laplacian=1; +do_ssprk2 =1; +H_projection=1; + +%% ------------------------ FDS params ------------------------ +Nx = Nx_in; Ny = Ny_in; % Number of cells inx and y. +Lx = 2*pi; Ly = 2*pi; % domain size +rho0 = 1.0; % density +mu0 = rho0*nu; % dynamic viscosity +U0 = 2.0; % TG amplitude +kx = 1; ky = 1; % TG wavenumbers +ux0 = 1.; uy0 = 1.; % TG constant advection velocities in x,y. +gvec = [0, 0]; % No gravity. +CFL = 0.25; % CFL constraint. +VN = 0.5; % Diffusion constraint. +t_begin = 0; % Time domain parameters.. +t_end=2.*pi; +dt = 0.001; +variable_time_step = 1; +dstep_plot=dstep_plot_in; + +%% ------------------------ Build mesh ------------------------- +% cells, faces, Uface (Nf×1), mu (Nc×1), rho (Nc×1), Lx, Ly +[cells, faces] = build_square_mesh_periodic_2d(Nx, Ny, Lx, Ly); +Nc = numel(cells); Nf = numel(faces); +% Convenience vectors +xc = zeros(Nc,2); +V = zeros(Nc,1); +for p = 1:Nc + xc(p,:) = cells(p).xc; + V(p) = cells(p).V; +end + +%% ------------------------ Fields init ------------------------ +% Taylor–Green vortex (2D, periodic box) +u = zeros(Nc,2); x = xc(:,1); y = xc(:,2); + +% Coordinates 2d: +% Our linear index is id(i,j) = (j-1)*Nx + i, so reshape(...,[Nx,Ny])'. +Xc = reshape(xc(:,1), [Nx,Ny])'; +Yc = reshape(xc(:,2), [Nx,Ny])'; + +rho = rho0 * ones(Nc,1); +mu = mu0 * ones(Nc,1); lambda = -2/3 * mu; % Stokes' hypothesis. + +% Initial Fields: +% Face-normal velocities: +Uface = zeros(Nf,1); +% Do faces first then reconstruct_u_from_faces: +for f = 1:Nf + P = faces(f).owner; N = faces(f).neigh; nf = faces(f).nf(:); + xf= faces(f).cf; + ubar(1) = -U0 * cos(kx*xf(1)) .* sin(ky*xf(2)) + ux0; + ubar(2) = U0 * sin(kx*xf(1)) .* cos(ky*xf(2)) + uy0; + Uface(f) = ubar * nf; % dot product +end + +% Cell centered velocities: +% u = -U0 cos(kx x) sin(ky y), v = U0 sin(kx x) cos(ky y) +u(:,1) = -U0 * cos(kx*x) .* sin(ky*y) + ux0; +u(:,2) = U0 * sin(kx*x) .* cos(ky*y) + uy0; + +% Cell centered pressure (not used in P0 fractional step): +p = -(rho0*U0^2/4) * (cos(2*kx*x) + cos(2*ky*y)); + +%% Try SSPRK2 : +params = struct; +params.Lx = Lx; +params.Ly = Ly; +params.rho0 = rho0; +params.gvec = gvec; % TG: no body force +params.CFL = CFL; +params.VN = VN; +params.divU_th = zeros(numel(cells),1); % incompressible +params.t_begin = 0.0; +params.t_end = t_end; +params.dt = dt; +params.variable_time_step = variable_time_step; +params.max_steps = 1e5; +params.add_meanK = true; +params.poisson_opts = struct('fix','pin','pin_idx',1,'use_ilu',true); +params.tg_params=struct('U0',U0,'kx',kx,'ky',ky,'nu',nu,'ux0',ux0,'uy0',uy0); + +% Cell at Nx/2, Ny/2: +% Choose the logical center +if mod(Nx,2)==0, i_mid = Nx/2; else, i_mid = (Nx+1)/2; end +if mod(Ny,2)==0, j_mid = Ny/2; else, j_mid = (Ny+1)/2; end +params.p_mid=find_cell_by_center_centered(cells,Nx,Ny,Lx,Ly,i_mid,j_mid); +fprintf('Center cell (i=%d,j=%d)->index p= %d\n',i_mid,j_mid,params.p_mid); + +figure +[u,Uface,p,H_cell,t_end,hist]= run_ssprk2(u, Uface, p, rho, mu, lambda, ... + cells, faces, Lx, Ly, params); + + + +return +end + +% Some code snippets: + +% Reconstruct cell centered u from face normal velocities (LS): +% u = reconstruct_u_from_faces_w(Uface, cells, faces); + +% % Rescale velocities to match mean analytical KE: +% Vol = sum(V); +% rho_mean = sum(V .* rho(:)) / Vol; +% K_tgt = rho_mean * (U0^2/4)*Vol + (ux0^2+uy0^2)/2*Vol; +% [u_scaled, Uface_scaled, alpha, K_curr, K_tgt] = ... +% rescale_to_target_ke(u, Uface, rho, cells, K_tgt); +% u = u_scaled; +% Uface = Uface_scaled; + +% face normal velocities as averages of cell centered velocities: +% for f = 1:Nf +% P = faces(f).owner; N = faces(f).neigh; nf = faces(f).nf(:); +% if N > 0 +% ubar = 0.5*(u(P,:)+u(N,:)); +% else +% ubar = zeros(1,2); +% end +% Uface(f) = ubar * nf; +% end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/taylor_green_uv_2d.m b/Utilities/Misc/geometry/colocated_schemes/taylor_green_uv_2d.m new file mode 100644 index 00000000000..3aa3451af24 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/taylor_green_uv_2d.m @@ -0,0 +1,11 @@ +function UV = taylor_green_uv_2d(x, y, t, U0, kx, ky, nu) + +global ux0 uy0 + +% x,y can be vectors or arrays (same size). Returns [u,v] stacked along 2nd dim. +k2 = kx*kx + ky*ky; +decay = exp(-nu*k2*t); +u = -U0 * cos(kx*(x-ux0*t)) .* sin(ky*(y-uy0*t)) * decay + ux0; +v = U0 * sin(kx*(x-ux0*t)) .* cos(ky*(y-uy0*t)) * decay + uy0; +UV = [u(:), v(:)]; +end \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/taylor_green_uv_2d_point.m b/Utilities/Misc/geometry/colocated_schemes/taylor_green_uv_2d_point.m new file mode 100644 index 00000000000..0b60babbce6 --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/taylor_green_uv_2d_point.m @@ -0,0 +1,11 @@ +function [up, vp] = taylor_green_uv_2d_point(xpi, ypi, t, U0, kx, ky, nu, ux0, uy0) +% TAYLOR_GREEN_UV_2D_POINT TG velocities at a single (xp,yp). +if nargin < 9, uy0 = 0; end +if nargin < 8, ux0 = 0; end +xp = xpi-ux0*t; +yp = ypi-uy0*t; +k2 = kx*kx + ky*ky; +decay = exp(-nu * k2 * t); +up = -U0 * cos(kx*xp) * sin(ky*yp) * decay + ux0; +vp = U0 * sin(kx*xp) * cos(ky*yp) * decay + uy0; +end diff --git a/Utilities/Misc/geometry/colocated_schemes/taylor_green_vortex.asv b/Utilities/Misc/geometry/colocated_schemes/taylor_green_vortex.asv new file mode 100644 index 00000000000..18ad92e4efe --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/taylor_green_vortex.asv @@ -0,0 +1,136 @@ +% taylor_green_vortex: +% +% Routine that calls Taylor-Green time integration for different number +% of cells in Nx and Ny, as well as nu, and plots conservation, error norms +% of colocated SSPRK2 solver. +% ------------------------------------------------------------------------- +close all +clear all +clc + +Nx=[8 16 32 64]; Ny=Nx; +nu=[0. 0.1]; +dstep_plot = 100; + +nNu = numel(nu); nNx = numel(Nx); +ErrM = nan(nNu,nNx); % L2 momentum abs error at T_end +RMSu = nan(nNu,nNx); % RMS u(mid) error +Nxv = Nx(:)'; % row +for inu=1:length(nu) + for ix=1:length(Nx) + + [params,hist,cells,faces]=taylor_green(Nx(ix),Ny(ix),nu(inu),dstep_plot); + + % Plot Mass (divergence), momentum and Kinetic Energy vs time: + plot_mass_mom_ke(hist, cells); + + % Plot L2 errors of momentum and KE vs time: + plot_L2_errors(hist); + + % --- grab last entry --- + idx = numel(hist.t); + tend = hist.t(idx); + nsteps= hist.step(idx); + + L2_Mom_err = hist.L2_m_abs(idx); + L2_KE_err = hist.L2_k_abs(idx); + + % --- print nicely --- + fprintf(' nu, Nx, Ny = %.6g, %s %.6g, %s %.6g, Final steps, time t = %5d, %s %.6g\n',nu(inu),' ',Nx(ix),' ',Ny(ix),nsteps,' ',tend); + fprintf(' L2 momentum abs error : %.3e\n', L2_Mom_err); + fprintf(' L2 kinetic energy abs : %.3e\n', L2_KE_err); + + % Now do L2 (RMS error) on the time series for velocities at p_mid: + U0 = params.tg_params.U0; + kx = params.tg_params.kx; + ky = params.tg_params.ky; + ux0= params.tg_params.ux0; + uy0= params.tg_params.uy0; + xp=cells(params.p_mid).xc(1); yp=cells(params.p_mid).xc(2); + for i=1:idx + t=hist.t(i); + [up, vp] = taylor_green_uv_2d_point(xp,yp,t,U0,kx,ky,nu(inu),ux0,uy0); + hist.u_ex_mid(i,1:2) = [up vp]; + end + + RMS_u_mid = sqrt(1/idx*sum((hist.u_num_mid(:,1)-hist.u_ex_mid(:,1)).^2)); + RMS_v_mid = sqrt(1/idx*sum((hist.u_num_mid(:,2)-hist.u_ex_mid(:,2)).^2)); + + fprintf(' RMS u error mid : %.3e\n', RMS_u_mid); + fprintf(' RMS v error mid : %.3e\n', RMS_v_mid); + + TG_Stats{inu,ix}.params = params; + TG_Stats{inu,ix}.hist = hist; + TG_Stats{inu,ix}.nsteps = nsteps; + TG_Stats{inu,ix}.L2_Mom_err = L2_Mom_err; + TG_Stats{inu,ix}.L2_KE_err = L2_KE_err; + TG_Stats{inu,ix}.RMS_u_mid = RMS_u_mid; + TG_Stats{inu,ix}.RMS_v_mid = RMS_v_mid; + + % --- Plot u(midpoint) vs time: numerical vs exact --- + t = hist.t(:); + u_num = hist.u_num_mid(:,1); % numerical u at midpoint + u_ex = hist.u_ex_mid(:,1); % exact u at midpoint + + figure('Name','Midpoint u vs time','Color','w'); + plot(t, u_num, '-o', 'LineWidth',1.8, 'MarkerSize',5, 'DisplayName','u_{num}(mid)'); + hold on + plot(t, u_ex, '--', 'LineWidth',2.0, 'DisplayName','u_{exact}(mid)'); + grid on + xlabel('t','FontSize',14); + ylabel('u at midpoint','FontSize',14); + title(['Midpoint u(t): numerical vs exact ' num2str(Nx(ix)) '^2 cells'],'FontSize',16); + legend('Location','best'); + set(gca,'FontSize',12); + + end + + % Plot L2 momentum error norm at T_end, and RMS error of u-velocity for + % mid - point: + h = TG_Stats{1,1}.params.Lx ./ Nxv; + for ix = 1:nNx + if ~isempty(TG_Stats{inu,ix}) + ErrM(inu,ix) = TG_Stats{inu,ix}.L2_Mom_err; + RMSu(inu,ix) = TG_Stats{inu,ix}.RMS_u_mid; + end + end + + % ==== Plot L2 momentum error vs grid ==== + figure + tiledlayout(1,1,'Padding','compact','TileSpacing','compact'); + + % L2 err momentum vs h = L_x / N_x + nexttile; + hold on + plot(h, ErrM(inu,:), '-s','LineWidth',1.8,'MarkerSize',6, ... + 'DisplayName', sprintf('\\nu=%.3g', nu(inu))); + plot(h,3*h,'--k','DisplayName', sprintf('O(\\Delta x)')) + plot(h,0.8*h.^2,'-k','DisplayName', sprintf('O(\\Delta x^2)')) + set(gca,'YScale','log','XScale','log','FontSize',14); grid on + xlabel('\Delta x','FontSize',16); ylabel('||M||_2 at T_{end}','FontSize',16); + title(['L2 momentum error vs h, nu=' num2str(nu(inu))],'FontSize',18); + legend('Location','best','FontSize',16); + axis([0.08 1. 10^-3 3]) + box on + + % ==== Plot RMS u(mid) error vs grid ==== + figure + tiledlayout(1,1,'Padding','compact','TileSpacing','compact'); + + % RMS error mid point u-velocity vs h + nexttile; + hold on + plot(h, RMSu(inu,:), '-s','LineWidth',1.8,'MarkerSize',6, ... + 'DisplayName', sprintf('Numerical \\nu=%.3g', nu(inu))); + plot(h,2*h,'--k','DisplayName', sprintf('O(\\Delta x)')) + plot(h,.3*h.^2,'-k','DisplayName', sprintf('O(\\Delta x^2)')) + set(gca,'YScale','log','XScale','log','FontSize',14); grid on + xlabel('\Delta x','FontSize',16); ylabel('U-velocity RMS error','FontSize',16); + title(['U-velocity RMS error vs h, nu=' num2str(nu(inu))],'FontSize',18); + legend('Location','best','FontSize',16); + axis([0.08 1 10^-3 2]) + box on + +end + +return \ No newline at end of file diff --git a/Utilities/Misc/geometry/colocated_schemes/taylor_green_vortex.m b/Utilities/Misc/geometry/colocated_schemes/taylor_green_vortex.m new file mode 100644 index 00000000000..18ad92e4efe --- /dev/null +++ b/Utilities/Misc/geometry/colocated_schemes/taylor_green_vortex.m @@ -0,0 +1,136 @@ +% taylor_green_vortex: +% +% Routine that calls Taylor-Green time integration for different number +% of cells in Nx and Ny, as well as nu, and plots conservation, error norms +% of colocated SSPRK2 solver. +% ------------------------------------------------------------------------- +close all +clear all +clc + +Nx=[8 16 32 64]; Ny=Nx; +nu=[0. 0.1]; +dstep_plot = 100; + +nNu = numel(nu); nNx = numel(Nx); +ErrM = nan(nNu,nNx); % L2 momentum abs error at T_end +RMSu = nan(nNu,nNx); % RMS u(mid) error +Nxv = Nx(:)'; % row +for inu=1:length(nu) + for ix=1:length(Nx) + + [params,hist,cells,faces]=taylor_green(Nx(ix),Ny(ix),nu(inu),dstep_plot); + + % Plot Mass (divergence), momentum and Kinetic Energy vs time: + plot_mass_mom_ke(hist, cells); + + % Plot L2 errors of momentum and KE vs time: + plot_L2_errors(hist); + + % --- grab last entry --- + idx = numel(hist.t); + tend = hist.t(idx); + nsteps= hist.step(idx); + + L2_Mom_err = hist.L2_m_abs(idx); + L2_KE_err = hist.L2_k_abs(idx); + + % --- print nicely --- + fprintf(' nu, Nx, Ny = %.6g, %s %.6g, %s %.6g, Final steps, time t = %5d, %s %.6g\n',nu(inu),' ',Nx(ix),' ',Ny(ix),nsteps,' ',tend); + fprintf(' L2 momentum abs error : %.3e\n', L2_Mom_err); + fprintf(' L2 kinetic energy abs : %.3e\n', L2_KE_err); + + % Now do L2 (RMS error) on the time series for velocities at p_mid: + U0 = params.tg_params.U0; + kx = params.tg_params.kx; + ky = params.tg_params.ky; + ux0= params.tg_params.ux0; + uy0= params.tg_params.uy0; + xp=cells(params.p_mid).xc(1); yp=cells(params.p_mid).xc(2); + for i=1:idx + t=hist.t(i); + [up, vp] = taylor_green_uv_2d_point(xp,yp,t,U0,kx,ky,nu(inu),ux0,uy0); + hist.u_ex_mid(i,1:2) = [up vp]; + end + + RMS_u_mid = sqrt(1/idx*sum((hist.u_num_mid(:,1)-hist.u_ex_mid(:,1)).^2)); + RMS_v_mid = sqrt(1/idx*sum((hist.u_num_mid(:,2)-hist.u_ex_mid(:,2)).^2)); + + fprintf(' RMS u error mid : %.3e\n', RMS_u_mid); + fprintf(' RMS v error mid : %.3e\n', RMS_v_mid); + + TG_Stats{inu,ix}.params = params; + TG_Stats{inu,ix}.hist = hist; + TG_Stats{inu,ix}.nsteps = nsteps; + TG_Stats{inu,ix}.L2_Mom_err = L2_Mom_err; + TG_Stats{inu,ix}.L2_KE_err = L2_KE_err; + TG_Stats{inu,ix}.RMS_u_mid = RMS_u_mid; + TG_Stats{inu,ix}.RMS_v_mid = RMS_v_mid; + + % --- Plot u(midpoint) vs time: numerical vs exact --- + t = hist.t(:); + u_num = hist.u_num_mid(:,1); % numerical u at midpoint + u_ex = hist.u_ex_mid(:,1); % exact u at midpoint + + figure('Name','Midpoint u vs time','Color','w'); + plot(t, u_num, '-o', 'LineWidth',1.8, 'MarkerSize',5, 'DisplayName','u_{num}(mid)'); + hold on + plot(t, u_ex, '--', 'LineWidth',2.0, 'DisplayName','u_{exact}(mid)'); + grid on + xlabel('t','FontSize',14); + ylabel('u at midpoint','FontSize',14); + title(['Midpoint u(t): numerical vs exact ' num2str(Nx(ix)) '^2 cells'],'FontSize',16); + legend('Location','best'); + set(gca,'FontSize',12); + + end + + % Plot L2 momentum error norm at T_end, and RMS error of u-velocity for + % mid - point: + h = TG_Stats{1,1}.params.Lx ./ Nxv; + for ix = 1:nNx + if ~isempty(TG_Stats{inu,ix}) + ErrM(inu,ix) = TG_Stats{inu,ix}.L2_Mom_err; + RMSu(inu,ix) = TG_Stats{inu,ix}.RMS_u_mid; + end + end + + % ==== Plot L2 momentum error vs grid ==== + figure + tiledlayout(1,1,'Padding','compact','TileSpacing','compact'); + + % L2 err momentum vs h = L_x / N_x + nexttile; + hold on + plot(h, ErrM(inu,:), '-s','LineWidth',1.8,'MarkerSize',6, ... + 'DisplayName', sprintf('\\nu=%.3g', nu(inu))); + plot(h,3*h,'--k','DisplayName', sprintf('O(\\Delta x)')) + plot(h,0.8*h.^2,'-k','DisplayName', sprintf('O(\\Delta x^2)')) + set(gca,'YScale','log','XScale','log','FontSize',14); grid on + xlabel('\Delta x','FontSize',16); ylabel('||M||_2 at T_{end}','FontSize',16); + title(['L2 momentum error vs h, nu=' num2str(nu(inu))],'FontSize',18); + legend('Location','best','FontSize',16); + axis([0.08 1. 10^-3 3]) + box on + + % ==== Plot RMS u(mid) error vs grid ==== + figure + tiledlayout(1,1,'Padding','compact','TileSpacing','compact'); + + % RMS error mid point u-velocity vs h + nexttile; + hold on + plot(h, RMSu(inu,:), '-s','LineWidth',1.8,'MarkerSize',6, ... + 'DisplayName', sprintf('Numerical \\nu=%.3g', nu(inu))); + plot(h,2*h,'--k','DisplayName', sprintf('O(\\Delta x)')) + plot(h,.3*h.^2,'-k','DisplayName', sprintf('O(\\Delta x^2)')) + set(gca,'YScale','log','XScale','log','FontSize',14); grid on + xlabel('\Delta x','FontSize',16); ylabel('U-velocity RMS error','FontSize',16); + title(['U-velocity RMS error vs h, nu=' num2str(nu(inu))],'FontSize',18); + legend('Location','best','FontSize',16); + axis([0.08 1 10^-3 2]) + box on + +end + +return \ No newline at end of file