Merge remote-tracking branch 'firemodels/master'

Author: cxp484

.github/workflows/linux.yml

@@ -34,9 +34,9 @@ jobs:
     - uses: rscohn2/setup-oneapi@v0
       with:
         components: |
-          ifx
-          impi
-          mkl
+          ifx@2025.2.0
+          impi@2021.16.0
+          mkl@2025.2.0
         prune: false       
 
     - uses: actions/checkout@v5
@@ -98,7 +98,7 @@ jobs:
     steps:
     - uses: actions/checkout@v5
 
-    - uses: actions/setup-python@v5
+    - uses: actions/setup-python@v6
       with:
         python-version: '3.9'
         cache: 'pip' # caching pip dependencies

.github/workflows/windows.yml

@@ -26,9 +26,9 @@ permissions:
 env:
   # update urls for oneapi packages according to
   # https://github.com/oneapi-src/oneapi-ci/blob/master/.github/workflows/build_all.yml
-  WINDOWS_BASEKIT_URL: https://registrationcenter-download.intel.com/akdlm/IRC_NAS/09a8acaf-265f-4460-866c-a3375ed5b4ff/intel-oneapi-base-toolkit-2025.2.0.591_offline.exe
+  WINDOWS_BASEKIT_URL: https://registrationcenter-download.intel.com/akdlm/IRC_NAS/f5881e61-dcdc-40f1-9bd9-717081ac623c/intel-oneapi-base-toolkit-2025.2.1.46_offline.exe
   WINDOWS_BASEKIT_COMPONENTS: intel.oneapi.win.mkl.devel
-  WINDOWS_HPCKIT_URL: https://registrationcenter-download.intel.com/akdlm/IRC_NAS/3bbdaf75-6728-492e-a18c-be654dae9ee2/intel-oneapi-hpc-toolkit-2025.2.0.576_offline.exe
+  WINDOWS_HPCKIT_URL: https://registrationcenter-download.intel.com/akdlm/IRC_NAS/e63ac2b4-8a9a-4768-979a-399a8b6299de/intel-oneapi-hpc-toolkit-2025.2.1.46_offline.exe
   WINDOWS_HPCKIT_COMPONENTS: intel.oneapi.win.ifort-compiler:intel.oneapi.win.mpi.devel
 
 

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

Build/makefile

@@ -355,7 +355,7 @@ ompi_gnu_linux : obj = fds_ompi_gnu_linux
 ompi_gnu_linux : setup $(obj_mpi)
 	$(FCOMPL) $(FFLAGS) $(FOPENMPFLAGS) -o $(obj) $(obj_mpi) $(LFLAGSMKL)
 
-ompi_gnu_linux_db : FFLAGS = -O0 -std=f2018 -ggdb -Wall -Werror -Wunused-parameter -Wcharacter-truncation -Wno-target-lifetime -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow -frecursive -ffpe-summary=none -fall-intrinsics -fbounds-check $(GITINFOGNU) $(FFLAGSMKL) $(GFORTRAN_OPTIONS) $(FFLAGS_HYPRE) $(FFLAGS_SUNDIALS)
+ompi_gnu_linux_db : FFLAGS = -O0 -std=f2018 -ggdb -finit-real=snan -finit-integer=-999999 -Wall -Werror -Wunused-parameter -Wcharacter-truncation -Wno-target-lifetime -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow -frecursive -ffpe-summary=none -fall-intrinsics -fbounds-check $(GITINFOGNU) $(FFLAGSMKL) $(GFORTRAN_OPTIONS) $(FFLAGS_HYPRE) $(FFLAGS_SUNDIALS)
 ompi_gnu_linux_db : LFLAGSMKL = $(LFLAGSMKL_GNU_OPENMP) $(LFLAGS_HYPRE) $(LFLAGS_SUNDIALS)
 ompi_gnu_linux_db : FCOMPL = $(COMP_FC)
 ompi_gnu_linux_db : FOPENMPFLAGS = -fopenmp

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
 

Source/wall.f90

@@ -162,9 +162,7 @@ SUBROUTINE WALL_BC(T,DT,NM)
       IF (B1%NODE_INDEX/=0) CALL CALC_HVAC_BC(BC,B1,SF)
    ENDIF
 
-   IF (WC%BOUNDARY_TYPE/=OPEN_BOUNDARY .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) THEN
-      CALL CALCULATE_ZZ_F(T,DT,WALL_INDEX=IW)
-   ENDIF
+   IF (WC%BOUNDARY_TYPE/=OPEN_BOUNDARY .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CALL CALCULATE_ZZ_F(T,DT,WALL_INDEX=IW)
 
    IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CALL CALCULATE_RHO_F(BC,B1,WALL_INDEX=IW)
 
@@ -220,7 +218,6 @@ SUBROUTINE WALL_BC(T,DT,NM)
    ENDIF
 
    CALL CALCULATE_ZZ_F(T,DT,CFACE_INDEX=ICF)
-
    CALL CALCULATE_RHO_F(BC,B1,CFACE_INDEX=ICF)
 
 ENDDO CFACE_LOOP
@@ -1328,9 +1325,9 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
 ! Only set species mass fraction in the ghost cell if it is not solid
 
 IF (PRESENT(WALL_INDEX)) THEN
-   IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(ICG)%SOLID) THEN
-      ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES)    = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES)
-      ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES)
+   IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS) THEN
+      ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES)    = B1%ZZ_F(1:N_TRACKED_SPECIES)
+      ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
    ENDIF
 ENDIF
 
@@ -1496,10 +1493,14 @@ SUBROUTINE CALCULATE_RHO_F(BC,B1,WALL_INDEX,CFACE_INDEX)
    ENDIF
 ENDIF
 
-IF (PRESENT(WALL_INDEX)) THEN
-   IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN
-      RHOP(BC%II ,BC%JJ, BC%KK ) = 2._EB*B1%RHO_F - B1%RHO_G
-      RHOP(BC%II2,BC%JJ2,BC%KK2) = 2._EB*B1%RHO_F - B1%RHO_G
+! Fill exterior ghost cells with surface density and temperature
+
+IF (PRESENT(WALL_INDEX) .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN
+   IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS) THEN
+      RHOP(BC%II ,BC%JJ, BC%KK ) = B1%RHO_F
+      RHOP(BC%II2,BC%JJ2,BC%KK2) = B1%RHO_F
+      TMP(BC%II ,BC%JJ, BC%KK ) = B1%TMP_F
+      TMP(BC%II2,BC%JJ2,BC%KK2) = B1%TMP_F
    ENDIF
 ENDIF