Skip to content

Modify the MPIFace process send-recv in HexMesh and MPIFace. Add test case for MLRK #271

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
hmarbona wants to merge 16 commits into from
Closed

Modify the MPIFace process send-recv in HexMesh and MPIFace. Add test case for MLRK #271

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
60f08d2
MPI partitioning for MLRK
Aug 7, 2025
99be5f2
Debug MPI partition for MLRK
Aug 7, 2025
974e24e
Debug MPI partition for MLRK
Aug 7, 2025
d5a6d7e
Debug MPI partition for MLRK
Aug 8, 2025
55ea95d
Fixed Memory Leak due to MPI Communication
Aug 12, 2025
0366c7a
Debug Memory Leak Fixes due to MPI Communications
Aug 12, 2025
f1fae21
Add option to direct partition without optimizing level for MLRK
Aug 13, 2025
b564f18
Reduce GitHub check for CI_Parallel due to wall clock limit
Aug 15, 2025
6ffc2c5
Add MLRK Cylinder test and memory monitor
Aug 24, 2025
c8c7ba1
Debug Residual Monitor
Aug 24, 2025
8d89087
Debug Cylinder_MLRK test
Aug 24, 2025
e92a585
Debug Residual Memory
Aug 24, 2025
4d93864
Minor reorder nullify
Sep 2, 2025
4d0512e
Merge branch 'loganoz:master' into master
hmarbona Sep 19, 2025
dd219ca
Update MLRK and Optimized MU Solver
Sep 19, 2025
23695ae
Debug
Sep 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
290 changes: 155 additions & 135 deletions .github/workflows/CI_parallel.yml
View file
Open in desktop

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions Solver/configure
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ TEST_CASES="./Euler/BoxAroundCircle \
./NavierStokes/Convergence_energy \
./NavierStokes/Convergence_entropy \
./NavierStokes/Cylinder \
./NavierStokes/Cylinder_MLRK \
./NavierStokes/CylinderFAS \
./NavierStokes/CylinderAnisFAS \
./NavierStokes/CylinderBR2 \
Expand Down
2 changes: 1 addition & 1 deletion Solver/src/CahnHilliardSolver/main.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,8 +168,8 @@ PROGRAM HORSES3DMainCH
call Stopwatch % destruct
CALL timeIntegrator % destruct()
CALL sem % destruct()
call Finalize_InterpolationMatrices
call DestructBoundaryConditions
call Finalize_InterpolationMatrices
call DestructGlobalNodalStorage()
CALL destructSharedBCModule

Expand Down
102 changes: 68 additions & 34 deletions Solver/src/MultiphaseSolver/SpatialDiscretization.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ end subroutine computeBoundaryFluxF
real(kind=RP), protected :: IMEX_S0 = 0.0_RP
real(kind=RP), protected :: IMEX_K0 = 1.0_RP
logical :: use_non_constant_speed_of_sound = .false.
integer :: speed_of_sound_model = 1
!
! ========
CONTAINS
Expand Down Expand Up @@ -192,6 +193,11 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh)
end select

use_non_constant_speed_of_sound = controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY)
if ( controlVariables % ContainsKey("speed of sound profile") .and. (trim(controlVariables % stringValueForKey('speed of sound profile', requestedLength = LINE_LENGTH)) == 'linear quadratic')) then
speed_of_sound_model = 2
else
speed_of_sound_model = 1
end if

call CHDiscretization % Construct(controlVariables, ELLIPTIC_CH)
call CHDiscretization % Describe
Expand All @@ -200,6 +206,12 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh)

if(use_non_constant_speed_of_sound) then
write(STD_OUT,'(A)') " Implementing artificial compressibility with a non-constant speed of sound in each phase"
if (speed_of_sound_model.eq.2) then
write(STD_OUT,'(A)') " Speed of sound profile: linear quadratic along interface"
else
write(STD_OUT,'(A)') " Speed of sound profile: linear along interface"
end if

else
write(STD_OUT,'(A)') " Implementing artificial compressibility with a constant ACM factor in each phase"
endif
Expand Down Expand Up @@ -238,7 +250,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
!
INTEGER :: k, eID, fID, i, j, ierr, locLevel, lID
real(kind=RP) :: sqrtRho, invMa2
class(Element), pointer :: e
logical :: compute_element
logical, allocatable :: face_mask(:)
real(kind=RP) :: mu_smag, delta
Expand Down Expand Up @@ -268,6 +279,24 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
!
select case (mode)
case (CTD_IGNORE_MODE,CTD_IMEX_EXPLICIT)

!
! --------------------------------------------
! Prolong Cahn-Hilliard and iNS to faces
! --------------------------------------------
!
call mesh % ProlongSolutionToFaces(NCONS, element_mask=element_mask, Level=locLevel)
!
! ----------------
! Update MPI Faces
! ----------------
!
#ifdef _HAS_MPI_
!$omp single
call mesh % UpdateMPIFacesSolution(NCONS)
!$omp end single
#endif

!$omp do schedule(runtime) private(eID)
do lID = 1, MLIter(locLevel,1)
eID = MLIter_eID(lID)
Expand All @@ -279,8 +308,24 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
endif
end do
!$omp end do
end select

#ifdef _HAS_MPI_
!$omp single
call mesh % GatherMPIFacesSolution(NCONS)
!$omp end single
#endif

! Not optimized for MLRK and MixedRK but okey
!$omp do schedule(runtime)
do fID = 1, size(mesh % faces)
associate(f => mesh % faces(fID))
f % storage(1) % c(1,:,:) = f % storage(1) % QNS(IMC,:,:)
f % storage(2) % c(1,:,:) = f % storage(2) % QNS(IMC,:,:)
end associate
end do
!$omp end do

end select
!
! -------------------------------
! Set memory to Cahn-Hilliard (C)
Expand All @@ -298,12 +343,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
end select
!$omp end single
!
! --------------------------------------------
! Prolong Cahn-Hilliard concentration to faces
! --------------------------------------------
!
call mesh % ProlongSolutionToFaces(NCOMP, element_mask=element_mask, Level=locLevel)
!
! ----------------
! Update MPI Faces
! ----------------
Expand Down Expand Up @@ -455,20 +494,13 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
call SetBoundaryConditionsEqn(NS_BC)
!$omp end single
!
! -------------------------
! Prolong solution to faces
! -------------------------
!
call mesh % ProlongSolutionToFaces(NCONS, element_mask=element_mask, Level=locLevel)
!
! ----------------
! Update MPI Faces
! ----------------
!
#ifdef _HAS_MPI_
!$omp single
call mesh % UpdateMPIFacesSolution(NCONS)
call mesh % GatherMPIFacesSolution(NCONS)
!$omp end single
#endif
!
Expand All @@ -488,8 +520,14 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
mesh % elements(eID) % storage % rho = min(max(mesh % elements(eID) % storage % rho, dimensionless % rho_min),dimensionless % rho_max)

if (use_non_constant_speed_of_sound ) then
mesh % elements(eID) % storage % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:),0.0_RP),1.0_RP)))**2
mesh % elements(eID) % storage % invMa2 = mesh % elements(eID) % storage % invMa2*mesh % elements(eID) % storage % rho
if (speed_of_sound_model.eq.1) then
! Linear profile
mesh % elements(eID) % storage % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:),0.0_RP),1.0_RP)))**2
mesh % elements(eID) % storage % invMa2 = mesh % elements(eID) % storage % invMa2*mesh % elements(eID) % storage % rho
else
! Linear quadratic profile
mesh % elements(eID) % storage % invMa2 = dimensionless % invMa2(2) + (dimensionless % invMa2(1) - dimensionless % invMa2(2)) * min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:), 0.0_RP), 1.0_RP)
end if
else
mesh % elements(eID) % storage % invMa2 = dimensionless % invMa2(1)
endif
Expand All @@ -506,7 +544,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
if (present(element_mask)) compute_element = face_mask(fID)

if (compute_element) then

mesh % faces(fID) % storage(1) % rho = dimensionless % rho(2) + (dimensionless % rho(1)-dimensionless % rho(2))* mesh % faces(fID) % storage(1) % Q(IMC,:,:)
mesh % faces(fID) % storage(2) % rho = dimensionless % rho(2) + (dimensionless % rho(1)-dimensionless % rho(2))* mesh % faces(fID) % storage(2) % Q(IMC,:,:)

Expand All @@ -515,11 +552,19 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem


if (use_non_constant_speed_of_sound ) then
mesh % faces(fID) % storage(1) % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:),0.0_RP),1.0_RP)))**2
mesh % faces(fID) % storage(2) % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:),0.0_RP),1.0_RP)))**2
if (speed_of_sound_model.eq.1) then
! Linear profile
mesh % faces(fID) % storage(1) % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:),0.0_RP),1.0_RP)))**2
mesh % faces(fID) % storage(2) % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:),0.0_RP),1.0_RP)))**2

mesh % faces(fID) % storage(1) % invMa2 = mesh % faces(fID) % storage(1) % invMa2*mesh % faces(fID) % storage(1) % rho
mesh % faces(fID) % storage(2) % invMa2 = mesh % faces(fID) % storage(2) % invMa2*mesh % faces(fID) % storage(2) % rho
mesh % faces(fID) % storage(1) % invMa2 = mesh % faces(fID) % storage(1) % invMa2*mesh % faces(fID) % storage(1) % rho
mesh % faces(fID) % storage(2) % invMa2 = mesh % faces(fID) % storage(2) % invMa2*mesh % faces(fID) % storage(2) % rho

else
! Linear quadratic profile
mesh % faces(fID) % storage(1) % invMa2 = dimensionless % invMa2(2) + (dimensionless % invMa2(1) - dimensionless % invMa2(2)) * min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:), 0.0_RP), 1.0_RP)
mesh % faces(fID) % storage(2) % invMa2 = dimensionless % invMa2(2) + (dimensionless % invMa2(1) - dimensionless % invMa2(2)) * min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:), 0.0_RP), 1.0_RP)
end if

else
mesh % faces(fID) % storage(1) % invMa2 = dimensionless % invMa2(1)
Expand All @@ -536,21 +581,11 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
!
call ViscousDiscretization % ComputeLocalGradients( NCONS, NCONS, mesh , time , mGradientVariables, Level = locLevel)
!
! --------------------
! Update MPI Gradients
! --------------------
!
#ifdef _HAS_MPI_
!$omp single
call mesh % UpdateMPIFacesGradients(NCONS)
!$omp end single
#endif
!
! -------------------------------------
! Add the Non-Conservative term to QDot
! -------------------------------------
!
!$omp do schedule(runtime) private(i,j,k,e,sqrtRho,invMa2,eID)
!$omp do schedule(runtime) private(i,j,k,sqrtRho,invMa2,eID)
do lID = 1, MLIter(locLevel,1) !
eID = MLIter_eID(lID)
compute_element = .true.
Expand Down Expand Up @@ -595,7 +630,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem

#ifdef _HAS_MPI_
!$omp single
! Not sure about the position of this w.r.t the MPI directly above
call mesh % UpdateMPIFacesGradients(NCONS)
call mesh % GatherMPIFacesGradients(NCONS)
!$omp end single
Expand Down
51 changes: 42 additions & 9 deletions Solver/src/MultiphaseSolver/main.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,9 @@ PROGRAM HORSES3DMainMU
TYPE( DGSem ) :: sem
TYPE( TimeIntegrator_t ) :: timeIntegrator
LOGICAL :: success, saveGradients, saveSource
logical :: generateMonitor = .TRUE.
logical :: optimizePartitionLevel=.FALSE.
logical :: isMLRK=.FALSE.
integer :: initial_iteration
INTEGER :: ierr
real(kind=RP) :: initial_time
Expand Down Expand Up @@ -88,23 +91,53 @@ PROGRAM HORSES3DMainMU
call GetMeshPolynomialOrders(controlVariables,Nx,Ny,Nz,Nmax)
call InitializeNodalStorage(controlVariables, Nmax)
call Initialize_InterpolationMatrices(Nmax)


! MPI Partition option for MLRK time integeration
if ((MPI_Process % doMPIAction).and.(trim(controlVariables % stringValueForKey('explicit method', requestedLength = LINE_LENGTH)) == 'multi level rk3')) then

isMLRK = .TRUE. ! The time integration is Multi-Level RK

if ( controlVariables % ContainsKey("optimized partition") ) then
optimizePartitionLevel = controlVariables % LogicalValueForKey ("optimized partition")
end if

if (optimizePartitionLevel) generateMonitor =.FALSE. ! Do not generate monitor for the first construction of sem as it will be reconstruct
end if

! Construct DGSEM library
call sem % construct ( controlVariables = controlVariables, &
Nx_ = Nx, Ny_ = Ny, Nz_ = Nz, &
success = success)

call Initialize_SpaceAndTimeMethods(controlVariables, sem % mesh)

IF(.NOT. success) error stop "Mesh reading error"
IF(.NOT. success) error stop "Boundary condition specification error"
CALL UserDefinedFinalSetup(sem % mesh, thermodynamics, dimensionless, refValues, multiphase)
Nx_ = Nx, Ny_ = Ny, Nz_ = Nz, &
success = success, generateMonitor =generateMonitor)
!
! -------------------------
! Set the initial condition
! -------------------------
!
call UserDefinedFinalSetup(sem % mesh, thermodynamics, dimensionless, refValues, multiphase)
call sem % SetInitialCondition(controlVariables, initial_iteration, initial_time)
!
! ----------------------------------------------
! Reconstruct for MLRK explicit time step method
! ----------------------------------------------
!
if (isMLRK .and. optimizePartitionLevel .and. MPI_Process % doMPIAction) then
call sem % reconstruct ( controlVariables = controlVariables, &
Nx_ = Nx, Ny_ = Ny, Nz_ = Nz, &
success = success)
call UserDefinedFinalSetup(sem % mesh, thermodynamics, dimensionless, refValues, multiphase)
call sem % SetInitialCondition(controlVariables, initial_iteration, initial_time)
end if
!
! -----------------------------
! Initialize the discretization
! -----------------------------
!
call Initialize_SpaceAndTimeMethods(controlVariables, sem % mesh)

IF(.NOT. success) error stop "Mesh reading error"
IF(.NOT. success) error stop "Boundary condition specification error"

!
! -----------------------------
! Construct the time integrator
! -----------------------------
Expand Down
51 changes: 43 additions & 8 deletions Solver/src/NavierStokesSolver/main.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ PROGRAM HORSES3DMainNS
use InterpolationMatrices , only: Initialize_InterpolationMatrices, Finalize_InterpolationMatrices
use ProblemFileFunctions
use BoundaryConditions , only: DestructBoundaryConditions
use PartitionedMeshClass
#ifdef _HAS_MPI_
use mpi
#endif
Expand All @@ -34,13 +35,17 @@ PROGRAM HORSES3DMainNS
TYPE( DGSem ) :: sem
TYPE( TimeIntegrator_t ) :: timeIntegrator
LOGICAL :: success, saveGradients, saveSensor, saveLES, saveSource
logical :: generateMonitor = .TRUE.
logical :: optimizePartitionLevel=.FALSE.
logical :: isMLRK=.FALSE.
integer :: initial_iteration
INTEGER :: ierr
real(kind=RP) :: initial_time
character(len=LINE_LENGTH) :: solutionFileName
integer, allocatable :: Nx(:), Ny(:), Nz(:)
integer :: Nmax


call SetSolver(NAVIERSTOKES_SOLVER)
!
! -----------------------------------------
Expand Down Expand Up @@ -96,21 +101,51 @@ PROGRAM HORSES3DMainNS
call InitializeNodalStorage (controlVariables ,Nmax)
call Initialize_InterpolationMatrices(Nmax)

! MPI Partition option for MLRK time integeration
if ((MPI_Process % doMPIAction).and.(trim(controlVariables % stringValueForKey('explicit method', requestedLength = LINE_LENGTH)) == 'multi level rk3')) then

isMLRK = .TRUE. ! The time integration is Multi-Level RK

if ( controlVariables % ContainsKey("optimized partition") ) then
optimizePartitionLevel = controlVariables % LogicalValueForKey ("optimized partition")
end if

if (optimizePartitionLevel) generateMonitor =.FALSE. ! Do not generate monitor for the first construction of sem as it will be reconstruct
end if

! Construct DGSEM library
call sem % construct ( controlVariables = controlVariables, &
Nx_ = Nx, Ny_ = Ny, Nz_ = Nz, &
success = success)

call Initialize_SpaceAndTimeMethods(controlVariables, sem)

IF(.NOT. success) error stop "Mesh reading error"
IF(.NOT. success) error stop "Boundary condition specification error"
CALL UserDefinedFinalSetup(sem % mesh, thermodynamics, dimensionless, refValues)
success = success, generateMonitor =generateMonitor)
!
! -------------------------
! Set the initial condition
! -------------------------
!
call sem % SetInitialCondition(controlVariables, initial_iteration, initial_time)
call UserDefinedFinalSetup(sem % mesh, thermodynamics, dimensionless, refValues)
call sem % SetInitialCondition(controlVariables, initial_iteration, initial_time)
!
! ----------------------------------------------
! Reconstruct for MLRK explicit time step method
! ----------------------------------------------
!
if (isMLRK .and. optimizePartitionLevel .and. MPI_Process % doMPIAction) then
call sem % reconstruct ( controlVariables = controlVariables, &
Nx_ = Nx, Ny_ = Ny, Nz_ = Nz, &
success = success)
call UserDefinedFinalSetup(sem % mesh, thermodynamics, dimensionless, refValues)
call sem % SetInitialCondition(controlVariables, initial_iteration, initial_time)
end if
!
! -----------------------------
! Initialize the discretization
! -----------------------------
!
call Initialize_SpaceAndTimeMethods(controlVariables, sem)

IF(.NOT. success) error stop "Mesh reading error"
IF(.NOT. success) error stop "Boundary condition specification error"

!
! -------------------
! Build the particles
Expand Down
Loading