Merge pull request #15087 from marcosvanella/master

FDS Source : UGLMAT HYPRE, compute only in MPI processes with nonzero unkowns.

Author: marcosvanella

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

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
 

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).
+!
+! ------------------------------------------------------------------------------------------------ 

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

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