Skip to content

Epicure-preprocess + more things #7

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

Open
loganoz wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Epicure-preprocess + more things #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5057d1e
Optimizations
Apr 20, 2025
dc8c4fe
Added HOPR Optimization
zalbanob May 14, 2025
3f62fab
Change get_rho in LESModel to be pure function such that compiles in cpu
oscarmarino May 22, 2025
afbdef1
Update AL new implementation.
oscarmarino May 22, 2025
ba0f3c2
delete prints
oscarmarino May 22, 2025
ae62452
Merge pull request #260 from zalbanob/epicure_bsc_gpu
loganoz Oct 29, 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
156 changes: 74 additions & 82 deletions Solver/src/NavierStokesSolver/SpatialDiscretization.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1354,90 +1354,82 @@ subroutine TimeDerivative_VolumetricContribution(mesh)

end subroutine TimeDerivative_VolumetricContribution

!
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
!

subroutine TimeDerivative_VolumetricContribution_Split(mesh)
use HexMeshClass
use NodalStorageClass, only: NodalStorage
use ElementClass
use DGIntegrals
use RiemannSolvers_NS
implicit none
type(HexMesh), intent (inout) :: mesh
!
! ---------------
! Local variables
! ---------------
!
integer :: i, j, k,l,eq, eID
real(kind=RP) :: Flux(1:NCONS, 1:NDIM)

!$acc parallel present(mesh) vector_length(128) num_gangs(9750) async(1)
!$acc loop gang
do eID = 1 , size(mesh % elements)

!$acc loop vector collapse(3) private(Flux)
do k = 0, mesh % elements(eID) % Nxyz(3)
do j = 0, mesh % elements(eID) % Nxyz(2)
do i = 0, mesh % elements(eID) % Nxyz(1)

call ViscousFlux_STATE( NCONS, NGRAD, mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % U_x(:,i,j,k), &
mesh % elements(eID) % storage % U_y(:,i,j,k) , mesh % elements(eID) % storage % U_z(:,i,j,k), &
mesh % elements(eID) % storage % mu_ns(1,i,j,k), 0.0_RP, &
mesh % elements(eID) % storage % mu_ns(2,i,j,k), Flux)
!$acc loop seq
do eq = 1, NCONS

mesh % elements(eID) % storage % contravariantFlux(eq,i,j,k,IX) = - Flux(eq,IX) * mesh % elements(eID) % geom % jGradXi(IX,i,j,k) &
- Flux(eq,IY) * mesh % elements(eID) % geom % jGradXi(IY,i,j,k) &
- Flux(eq,IZ) * mesh % elements(eID) % geom % jGradXi(IZ,i,j,k)

mesh % elements(eID) % storage % contravariantFlux(eq,i,j,k,IY) = - Flux(eq,IX) * mesh % elements(eID) % geom % jGradEta(IX,i,j,k) &
- Flux(eq,IY) * mesh % elements(eID) % geom % jGradEta(IY,i,j,k) &
- Flux(eq,IZ) * mesh % elements(eID) % geom % jGradEta(IZ,i,j,k)

mesh % elements(eID) % storage % contravariantFlux(eq,i,j,k,IZ) = - Flux(eq,IX) * mesh % elements(eID) % geom % jGradZeta(IX,i,j,k) &
- Flux(eq,IY) * mesh % elements(eID) % geom % jGradZeta(IY,i,j,k) &
- Flux(eq,IZ) * mesh % elements(eID) % geom % jGradZeta(IZ,i,j,k)
end do
end do ; end do ; end do

call ScalarWeakIntegrals_StdVolumeGreen( mesh % elements(eID) % Nxyz, NCONS, mesh % elements(eID) % storage % contravariantFlux, &
mesh % elements(eID) % storage % QDot)
!
! *************************************
! Compute interior contravariant fluxes
! *************************************
!
! Compute inviscid contravariant flux
! -----------------------------------
!$acc loop vector collapse(3) private(Flux)
do k = 0, mesh % elements(eID) % Nxyz(3)
do j = 0, mesh % elements(eID) % Nxyz(2)
do i = 0, mesh % elements(eID) % Nxyz(1)
!$acc loop seq
do l = 0, mesh % elements(eID) % Nxyz(1)
call TwoPointFlux_Selector(mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % Q(:,l,j,k), mesh % elements(eID) % geom % jGradXi(:,i,j,k), mesh % elements(eID) % geom % jGradXi(:,l,j,k), Flux(:,IX))
call TwoPointFlux_Selector(mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % Q(:,i,l,k), mesh % elements(eID) % geom % jGradEta(:,i,j,k), mesh % elements(eID) % geom % jGradEta(:,i,l,k), Flux(:,IY))
call TwoPointFlux_Selector(mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % Q(:,i,j,l), mesh % elements(eID) % geom % jGradZeta(:,i,j,k), mesh % elements(eID) % geom % jGradZeta(:,i,j,l), Flux(:,IZ))

!$acc loop seq
do eq = 1, NCONS
mesh % elements(eID) % storage % QDot(eq,i,j,k) = mesh % elements(eID) % storage % QDot(eq,i,j,k) &
- NodalStorage(mesh % elements(eID) % Nxyz(1)) % sharpD(i,l) * Flux(eq,IX) &
- NodalStorage(mesh % elements(eID) % Nxyz(2)) % sharpD(j,l) * Flux(eq,IY) &
- NodalStorage(mesh % elements(eID) % Nxyz(3)) % sharpD(k,l) * Flux(eq,IZ)
end do
end do

end do ; end do ; end do

enddo
!$acc end parallel loop

end subroutine TimeDerivative_VolumetricContribution_Split
#define PASTE(a) a
#define CAT(a,b) PASTE(a)b

#define SUFFIX StandardDG
#define TWO_POINT_FLUX_FUNC StandardDG_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC

#define SUFFIX Morinishi
#define TWO_POINT_FLUX_FUNC Morinishi_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC

#define SUFFIX Ducros
#define TWO_POINT_FLUX_FUNC Ducros_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC

#define SUFFIX KennedyGruber
#define TWO_POINT_FLUX_FUNC KennedyGruber_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC

#define SUFFIX Pirozzoli
#define TWO_POINT_FLUX_FUNC Pirozzoli_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC

#define SUFFIX EntropyConserving
#define TWO_POINT_FLUX_FUNC EntropyConserving_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC

#define SUFFIX Chandrasekar
#define TWO_POINT_FLUX_FUNC Chandrasekar_TwoPointFlux
#include "split_template.f90"
#undef SUFFIX
#undef TWO_POINT_FLUX_FUNC


subroutine TimeDerivative_VolumetricContribution_Split(mesh)
use HexMeshClass
use RiemannSolvers_NSKeywordsModule
use RiemannSolvers_NS, only : whichAverage
implicit none
type(HexMesh), intent(inout) :: mesh

select case (whichAverage)
case (STANDARD_AVG)
call TimeDerivative_VolumetricContribution_Split_StandardDG(mesh)
case (MORINISHI_AVG)
call TimeDerivative_VolumetricContribution_Split_Morinishi(mesh)
case (DUCROS_AVG)
call TimeDerivative_VolumetricContribution_Split_Ducros(mesh)
case (KENNEDYGRUBER_AVG)
call TimeDerivative_VolumetricContribution_Split_KennedyGruber(mesh)
case (PIROZZOLI_AVG)
call TimeDerivative_VolumetricContribution_Split_Pirozzoli(mesh)
case (ENTROPYCONS_AVG)
call TimeDerivative_VolumetricContribution_Split_EntropyConserving(mesh)
case (CHANDRASEKAR_AVG)
call TimeDerivative_VolumetricContribution_Split_Chandrasekar(mesh)
case default
print*, "Averaging not recognized."
errorMessage(STD_OUT)
error stop
end select
end subroutine
!
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
!
Expand Down
92 changes: 92 additions & 0 deletions Solver/src/NavierStokesSolver/split_template.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
CAT(subroutine TimeDerivative_VolumetricContribution_Split_,SUFFIX)(mesh)
use HexMeshClass
use NodalStorageClass, only: NodalStorage
use ElementClass
use DGIntegrals
use RiemannSolvers_NS
implicit none
type(HexMesh), intent (inout) :: mesh

! ---------------
! Local variables
! ---------------
integer :: i, j, k,l,eq, eID
real(kind=RP) :: Flux(1:NCONS, 1:NDIM), Q_acc

!$acc parallel present(mesh) vector_length(128) num_gangs(9750) async(1)
!$acc loop gang
do eID = 1 , size(mesh % elements)

!$acc loop vector collapse(3) private(Flux)
do k = 0, mesh % elements(eID) % Nxyz(3)
do j = 0, mesh % elements(eID) % Nxyz(2)
do i = 0, mesh % elements(eID) % Nxyz(1)
call ViscousFlux_STATE( NCONS, NGRAD, mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % U_x(:,i,j,k), &
mesh % elements(eID) % storage % U_y(:,i,j,k) , mesh % elements(eID) % storage % U_z(:,i,j,k), &
mesh % elements(eID) % storage % mu_ns(1,i,j,k), 0.0_RP, &
mesh % elements(eID) % storage % mu_ns(2,i,j,k), Flux)
!$acc loop seq
do eq = 1, NCONS

mesh % elements(eID) % storage % contravariantFlux(eq,i,j,k,IX) = - Flux(eq,IX) * mesh % elements(eID) % geom % jGradXi(IX,i,j,k) &
- Flux(eq,IY) * mesh % elements(eID) % geom % jGradXi(IY,i,j,k) &
- Flux(eq,IZ) * mesh % elements(eID) % geom % jGradXi(IZ,i,j,k)

mesh % elements(eID) % storage % contravariantFlux(eq,i,j,k,IY) = - Flux(eq,IX) * mesh % elements(eID) % geom % jGradEta(IX,i,j,k) &
- Flux(eq,IY) * mesh % elements(eID) % geom % jGradEta(IY,i,j,k) &
- Flux(eq,IZ) * mesh % elements(eID) % geom % jGradEta(IZ,i,j,k)

mesh % elements(eID) % storage % contravariantFlux(eq,i,j,k,IZ) = - Flux(eq,IX) * mesh % elements(eID) % geom % jGradZeta(IX,i,j,k) &
- Flux(eq,IY) * mesh % elements(eID) % geom % jGradZeta(IY,i,j,k) &
- Flux(eq,IZ) * mesh % elements(eID) % geom % jGradZeta(IZ,i,j,k)
end do
end do
end do
end do

!$acc loop vector collapse(4)
do k = 0, mesh % elements(eID) % Nxyz(3)
do j = 0, mesh % elements(eID) % Nxyz(2)
do i = 0, mesh % elements(eID) % Nxyz(1)
do eq = 1, NCONS
Q_acc = 0.0_RP
!$acc loop seq
do l = 0, mesh % elements(eID) % Nxyz(1)
Q_acc = Q_acc + NodalStorage(mesh % elements(eID) % Nxyz(1)) % hatD(i,l) * mesh % elements(eID) % storage % contravariantFlux(eq,l,j,k,IX) + &
NodalStorage(mesh % elements(eID) % Nxyz(2)) % hatD(j,l) * mesh % elements(eID) % storage % contravariantFlux(eq,i,l,k,IY) + &
NodalStorage(mesh % elements(eID) % Nxyz(3)) % hatD(k,l) * mesh % elements(eID) % storage % contravariantFlux(eq,i,j,l,IZ)
end do
mesh % elements(eID) % storage % QDot(eq,i,j,k) = Q_acc
end do
end do
end do
end do

! *************************************
! Compute interior contravariant fluxes
! *************************************
! Compute inviscid contravariant flux
! -----------------------------------
!$acc loop vector collapse(3) private(Flux)
do k = 0, mesh % elements(eID) % Nxyz(3)
do j = 0, mesh % elements(eID) % Nxyz(2)
do i = 0, mesh % elements(eID) % Nxyz(1)
!$acc loop seq
do l = 0, mesh % elements(eID) % Nxyz(1)
call TWO_POINT_FLUX_FUNC( mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % Q(:,l,j,k), mesh % elements(eID) % geom % jGradXi(:,i,j,k), mesh % elements(eID) % geom % jGradXi(:,l,j,k), Flux(:,IX))
call TWO_POINT_FLUX_FUNC( mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % Q(:,i,l,k), mesh % elements(eID) % geom % jGradEta(:,i,j,k), mesh % elements(eID) % geom % jGradEta(:,i,l,k), Flux(:,IY))
call TWO_POINT_FLUX_FUNC( mesh % elements(eID) % storage % Q(:,i,j,k), mesh % elements(eID) % storage % Q(:,i,j,l), mesh % elements(eID) % geom % jGradZeta(:,i,j,k), mesh % elements(eID) % geom % jGradZeta(:,i,j,l), Flux(:,IZ))
!$acc loop seq
do eq = 1, NCONS
mesh % elements(eID) % storage % QDot(eq,i,j,k) = mesh % elements(eID) % storage % QDot(eq,i,j,k) &
- NodalStorage(mesh % elements(eID) % Nxyz(1)) % sharpD(i,l) * Flux(eq,IX) &
- NodalStorage(mesh % elements(eID) % Nxyz(2)) % sharpD(j,l) * Flux(eq,IY) &
- NodalStorage(mesh % elements(eID) % Nxyz(3)) % sharpD(k,l) * Flux(eq,IZ)
end do
end do
end do
end do
end do
enddo
!$acc end parallel loop
CAT(end subroutine TimeDerivative_VolumetricContribution_Split_,SUFFIX)
82 changes: 57 additions & 25 deletions Solver/src/libs/mesh/HexMesh.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4089,12 +4089,15 @@ subroutine HexMesh_ComputeWallDistances(self,facesList,elementList)
else
num_of_elems = self % no_of_elements
end if

!$acc data copyin(self, Xwall)
do ii = 1, num_of_elems
if ( present(elementList) ) then
eID = elementList (ii)
else
eID = ii
end if

associate(e => self % elements(eID))
allocate(e % geom % dWall(0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3)))
if( self% IBM% active ) then
Expand All @@ -4103,20 +4106,35 @@ subroutine HexMesh_ComputeWallDistances(self,facesList,elementList)
endif

if( .not. self% IBM% active ) then
do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
xP = e % geom % x(:,i,j,k)

minimumDistance = HUGE(1.0_RP)
do fID = 1, no_of_wallDOFS
currentDistance = sum(POW2(xP - Xwall(:,fID)))
minimumDistance = min(minimumDistance, currentDistance)
if (no_of_wallDOFS .ne. 0) then
!$acc data copyin(self % elements(eID))
!$acc data copyin(self % elements(eID) % geom)
!$acc data copyin(self % elements(eID) % geom % x)
!$acc data create(self % elements(eID) % geom % dWall)
do k = 0, self % elements(eID) % Nxyz(3)
!$acc parallel loop gang collapse(2) private(xP) vector_length(512) async(k)
do j = 0, self % elements(eID) % Nxyz(2)
do i = 0, self % elements(eID) % Nxyz(1)
xP = self % elements(eID) % geom % x(:, i, j, k)
minimumDistance = HUGE(1.0_RP)
!$acc loop gang vector reduction(min:minimumDistance) private(currentDistance, minimumDistance)
do fID = 1, no_of_wallDOFS
currentDistance = sum(POW2(xP - Xwall(:, fID)))
minimumDistance = min(minimumDistance, currentDistance)
end do
self % elements(eID) % geom % dWall(i, j, k) = sqrt(minimumDistance)
end do
end do
!$acc end parallel loop
!$acc update self(self % elements(eID) % geom % dWall(:, :, k)) async(k)
end do

e % geom % dWall(i,j,k) = sqrt(minimumDistance)

end do ; end do ; end do
!$acc wait
!$acc end data
!$acc end data
!$acc end data
!$acc end data
endif
end if

end associate
end do
!
Expand All @@ -4141,22 +4159,36 @@ subroutine HexMesh_ComputeWallDistances(self,facesList,elementList)
endif

if( .not. self% IBM% active ) then
do j = 0, fe % Nf(2) ; do i = 0, fe % Nf(1)
xP = fe % geom % x(:,i,j)

minimumDistance = HUGE(1.0_RP)
do fID = 1, no_of_wallDOFS
currentDistance = sum(POW2(xP - Xwall(:,fID)))
minimumDistance = min(minimumDistance, currentDistance)
if (no_of_wallDOFS .ne. 0 ) then
!$acc data copyin(self % faces(eID))
!$acc data copyin(self % faces(eID) % geom)
!$acc data copyin(self % faces(eID) % geom % x)
!$acc data create(self % faces(eID) % geom % dWall)
do j = 0, self % faces(eID) % Nf(2)
!$acc parallel loop gang private(xP) vector_length(512) async(j)
do i = 0, self % faces(eID) % Nf(1)
xP = self % faces(eID) % geom % x(:, i, j)
minimumDistance = HUGE(1.0_RP)
!$acc loop vector reduction(min:minimumDistance) private(currentDistance, minimumDistance)
do fID = 1, no_of_wallDOFS
currentDistance = sum(POW2(xP - Xwall(:, fID)))
minimumDistance = min(minimumDistance, currentDistance)
end do
self % faces(eID) % geom % dWall(i, j) = sqrt(minimumDistance)
end do
!$acc end parallel loop
!$acc update self(self % faces(eID) % geom % dWall(:j)) async(j)
end do

fe % geom % dWall(i,j) = sqrt(minimumDistance)

end do ; end do
end if

!$acc wait
!$acc end data
!$acc end data
!$acc end data
!$acc end data
endif
endif
end associate
end do
!$acc end data

deallocate(Xwall)

Expand Down
Loading