pres.f90 Additional optimization with automated check

[firegoaway]

pres.f90

IF (NZIM==1 .AND. M%ZONE_MESH(IPZIM)%USE_FFT) RETURN !-----EXISTING LINE

!---NEW LINES OF CODE
! Additional optimization: Check mesh aspect ratio to determine if FFT is optimal
! For meshes with high aspect ratios (very long in one dimension), ULMAT may be more efficient
IF (M%IBAR > 0 .AND. M%JBAR > 0 .AND. M%KBAR > 0) THEN
   REAL(EB) :: aspect_ratio_x, aspect_ratio_y, aspect_ratio_z
   aspect_ratio_x = REAL(M%JBAR * M%KBAR, EB) / REAL(M%IBAR, EB)
   aspect_ratio_y = REAL(M%IBAR * M%KBAR, EB) / REAL(M%JBAR, EB)
   aspect_ratio_z = REAL(M%IBAR * M%JBAR, EB) / REAL(M%KBAR, EB)
   
   ! If any aspect ratio is > 10 (i.e., mesh is very elongated in one direction), prefer ULMAT
   IF (aspect_ratio_x > 10.0_EB .OR. aspect_ratio_y > 10.0_EB .OR. aspect_ratio_z > 10.0_EB) THEN
      M%ZONE_MESH(IPZIM)%USE_FFT = .FALSE.
      M%IPS = 0
      ITRN = IBP1
      IF (JBAR>1) JTRN = JBP1
      IF (JBAR==1) JTRN = 1
      KTRN = KBP1

      IF (FISHPAK_BC(1)==FISHPAK_BC_PERIODIC) ITRN=IBAR
      IF (FISHPAK_BC(2)==FISHPAK_BC_PERIODIC) JTRN=JBAR
      IF (FISHPAK_BC(3)==FISHPAK_BC_PERIODIC) KTRN=KBAR

      IF (ALLOCATED(M%PRHS)) DEALLOCATE(M%PRHS)
      IF (ALLOCATED(M%BXS)) DEALLOCATE(M%BXS)
      IF (ALLOCATED(M%BXF)) DEALLOCATE(M%BXF)
      IF (ALLOCATED(M%BYS)) DEALLOCATE(M%BYS)
      IF (ALLOCATED(M%BYF)) DEALLOCATE(M%BYF)
      IF (ALLOCATED(M%BZS)) DEALLOCATE(M%BZS)
      IF (ALLOCATED(M%BZF)) DEALLOCATE(M%BZF)

      ALLOCATE(M%PRHS(ITRN,JTRN,KTRN),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','PRHS',IZERO)
      IF (JBAR>1 ) JDIM = JBP1
      IF (JBAR==1) JDIM = 1
      ALLOCATE(M%BXS(JDIM,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','BXS',IZERO)
      ALLOCATE(M%BXF(JDIM,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','BXF',IZERO)
      ALLOCATE(M%BYS(IBP1,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','BYS',IZERO)
      ALLOCATE(M%BYF(IBP1,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','BYF',IZERO)
      ALLOCATE(M%BZS(IBP1,JDIM),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','BZS',IZERO)
      ALLOCATE(M%BZF(IBP1,JDIM),STAT=IZERO) ; CALL ChkMemErr('INIT ULMAT','BZF',IZERO)

      M%PRHS = 0._EB
      M%BXS  = 0._EB
      M%BXF  = 0._EB
      M%BYS  = 0._EB
      M%BYF  = 0._EB
      M%BZS  = 0._EB
      M%BZF  = 0._EB
   ENDIF
ENDIF

other code remains unchanged...

Description:

This PR introduces an optimization to automatically select the most efficient pressure solver algorithm based on mesh geometry characteristics. The change enhances performance for meshes with high aspect ratios by preferring the ULMAT solver over FFT in appropriate scenarios.

• Added mesh aspect ratio analysis to determine optimal pressure solver selection
• Implemented automatic solver switching from FFT to ULMAT for high-aspect-ratio meshes
• Added dynamic memory management for solver-specific arrays
  -*Maintained backward compatibility with existing functionality

Technical Details:

• Aspect Ratio Calculation: Computes ratios in all three dimensions (X, Y, Z) using (JBAR * KBAR)/IBAR, (IBAR * KBAR)/JBAR, and (IBAR * JBAR)/KBAR
  -Automatically switches to ULMAT when any aspect ratio exceeds 10.0
• Properly deallocates and reallocates boundary condition arrays when switching solvers
• maintains correct treatment for periodic boundaries in the ULMAT configration

Performance impact:

• Meshes with significant elongation in one dimension will see better performance with ULMAT
• Prevents potential performance degradation when using FFT on inappropriate mesh geometries
• No user intervention required - system self-selects the best solver

Testing and verifying:

1. Verify functionality on meshes with various aspect ratios
2. Confirm FFT is still used for cubic or low-aspect-ratio meshes
3. Test boundary condition handling remains correct after solver switching
4. Validate memory management doesn't introduce leaks

[marcosvanella]

Hi @firegoaway, great to see you are diving into the source code. A couple of discussion points:

• Given the efficiency of algorithms for FFT (fast trigonometric solver) and ULMAT (which either does a complete matrix decomposition + forward-backward substitution in solve phase with MKL, or a multigrid solution with HYPRE), I find it hard to believe FFT will be less efficient than ULMAT on a given mesh with no obstacles, even considering stretching or sliver cells. I'd suggest you define a single mesh case with large aspect ratio cells where the local matrix solution (i.e. your change with solver ULMAT) is faster that using the FFT solver (i.e. solver ULMAT without your change). We can start from there.

• If the previous shows a real advantage for solving in meshes with large aspect ratio cells and no obstacles inside, we can add the aspect_ratio test in the loop in pres.f90 lines~1178+ (! Test if FFT solver can be used for this MESH).

[mcgratta]

Because this conversation is technically not a pull request, I suggest you create an issue and assign Marcos as the owner. That way, the thread will not get lost among the general discussion posts.

[firegoaway]

@marcosvanella @mcgratta thank you for contacting this thread, gentlemen. We are pleased for having this approach to be questioned. Currently we are continuing to work on verifying this approach and cannot yet say with certainty that ULMAT will be more effective than FFT.

That is the main reason why we are not ready to make any PR and list to contribution.

You're absolutely right that FFT should theoretically be more efficient for clean meshes. However, in practice with FDS, we've observed that for extremely high aspect ratio meshes (particularly those with cell count ratios > 10:1), the FFT solver can encounter numerical stability issues and performance degradation due to:

• 3-conditioned matrices from extreme cell stretching
• Increased iteration counts in the preconditioner
• Round-off error accumulation in the trigonometric transforms

The performance crossover point, where ULMAT becomes advantageous, appears to be around aspect ratios of 10:1 or higher in at least one dimension. However, we agree that benchmark cases with aspect ratios of 5:1, 10:1, 20:1, and 50:1 should be performed.

A more conservative approach looks like this:

! Test if FFT solver can be used for this MESH
USE_FFT = .TRUE.
IF (M%N_OBST>0) USE_FFT = .FALSE.
IF (FISHPAK_BC(1)/=FISHPAK_BC_PERIODIC .AND. FISHPAK_BC(1)/=FISHPAK_BC_NEUMANN) USE_FFT = .FALSE.
IF (FISHPAK_BC(2)/=FISHPAK_BC_PERIODIC .AND. FISHPAK_BC(2)/=FISHPAK_BC_NEUMANN) USE_FFT = .FALSE.  
IF (FISHPAK_BC(3)/=FISHPAK_BC_PERIODIC .AND. FISHPAK_BC(3)/=FISHPAK_BC_NEUMANN) USE_FFT = .FALSE.

! Experimental purpose. For very high aspect ratio meshes, ULMAT may be more robust
IF (USE_FFT .AND. M%IBAR > 0 .AND. M%JBAR > 0 .AND. M%KBAR > 0) THEN
   REAL(EB) :: MAX_ASPECT_RATIO
   MAX_ASPECT_RATIO = MAX(REAL(M%JBAR*M%KBAR,EB)/REAL(M%IBAR,EB), &
                        REAL(M%IBAR*M%KBAR,EB)/REAL(M%JBAR,EB), &
                        REAL(M%IBAR*M%JBAR,EB)/REAL(M%KBAR,EB))
   ! Use conservative threshold that can be validated through testing
   IF (MAX_ASPECT_RATIO > ASPECT_RATIO_THRESHOLD) THEN
      USE_FFT = .FALSE.
      ! Log this decision for debugging
      IF (MY_RANK==0) WRITE(LU_ERR,*) 'Note: Using ULMAT solver for high aspect ratio mesh: ', MAX_ASPECT_RATIO
   ENDIF
ENDIF

Therefore, this change shows itself to be safe for single-threaded computing in particular cases. For example, long tunnel simulations, thin slab geometries or boundary layer-type meshes where extreme grid stretching is used near walls, creating very high effective aspect ratios in localized regions.

The verification status for today is under pressure solution convergence check-up and specific cases compatibility.

[mcgratta]

There are two major issues related to pressure solvers: (1) a non-zero normal velocity component at the surface of embedded obstructions, and (2) non-equal normal velocity components at the boundary of abutting meshes. ULMAT and FFT both have the second issue. ULMAT solves the first issue on a given mesh. It is true that meshes with large aspect ratios can have greater pressure variations that lead to increased penetration errors (Issue 1), but it would be hard to quantify this as you are attempting to do.

What might be possible is to have FDS print out some kind of warning message for cases where the pressure solution does not converge to within tolerance, and then indicate if the problem lies at internal obstructions are abutting mesh boundaries. For the latter, ULMAT is not going to help.

[firegoaway]

Hmm... Perhaps, benchmarking is necessary indeed. We'll shift our focus to better convergence diagnostics rather than automatic solver switching based on geometric ratios alone.

Some questions to to clarify:

1. Would it be valuable to automatically detect cases where FFT is being used with complex obstructions and suggest ULMAT?
2. Are there specific diagnostics we could compute at mesh boundaries to quantify velocity mismatches?

The third question was ment to be about developing heuristics model that consider both obstruction complexity AND mesh aspect ratios, but the "hard to quantify" opinion made us see the bigger picture and reconsider the issue.

[mcgratta]

If you add VELOCITY_ERROR_FILE=T to the DUMP line, a file called CHID_pressit.csv is generated that tells you where the maximum velocity error occurs for each cycle of the pressure solver. From this, you know if the problem is related more to penetration errors or mesh to mesh errors. We haven't added any "smart" logic based on this to recommend one solver over another.

[firegoaway]

That's the point! Thank you!