From 5057d1e308bc6e9ea6e2b56bb22cb5803999b629 Mon Sep 17 00:00:00 2001 From: zalbanob Date: Sun, 20 Apr 2025 02:47:26 +0200 Subject: [PATCH 1/5] Optimizations --- .../SpatialDiscretization.f90 | 156 +++++++++--------- .../src/NavierStokesSolver/split_template.f90 | 92 +++++++++++ Solver/src/libs/mesh/HexMesh.f90 | 77 ++++++--- .../navierstokes/RiemannSolvers_NS.f90 | 1 + 4 files changed, 220 insertions(+), 106 deletions(-) create mode 100644 Solver/src/NavierStokesSolver/split_template.f90 diff --git a/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 b/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 index 93f40f593..efd87ff09 100644 --- a/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 +++ b/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 @@ -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 ! !/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! diff --git a/Solver/src/NavierStokesSolver/split_template.f90 b/Solver/src/NavierStokesSolver/split_template.f90 new file mode 100644 index 000000000..530689fc4 --- /dev/null +++ b/Solver/src/NavierStokesSolver/split_template.f90 @@ -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) \ No newline at end of file diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index 324097fa8..5bd1425a4 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -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 @@ -4103,20 +4106,33 @@ 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) + !$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 - - e % geom % dWall(i,j,k) = sqrt(minimumDistance) - - end do ; end do ; end do + !$acc end parallel loop + !$acc update self(self % elements(eID) % geom % dWall(:, :, k)) async(k) + end do + !$acc wait + !$acc end data + !$acc end data + !$acc end data + !$acc end data end if - end associate end do ! @@ -4141,22 +4157,35 @@ 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) + !$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 - - fe % geom % dWall(i,j) = sqrt(minimumDistance) - - end do ; end do + !$acc end parallel loop + !$acc update self(self % faces(eID) % geom % dWall(:j)) async(j) + end do + !$acc wait + !$acc end data + !$acc end data + !$acc end data + !$acc end data end if - end associate end do + !$acc end data deallocate(Xwall) diff --git a/Solver/src/libs/physics/navierstokes/RiemannSolvers_NS.f90 b/Solver/src/libs/physics/navierstokes/RiemannSolvers_NS.f90 index 8bbceca25..c30092445 100644 --- a/Solver/src/libs/physics/navierstokes/RiemannSolvers_NS.f90 +++ b/Solver/src/libs/physics/navierstokes/RiemannSolvers_NS.f90 @@ -66,6 +66,7 @@ module RiemannSolvers_NS public SetRiemannSolver, DescribeRiemannSolver public AveragedStates, TwoPointFlux, RiemannSolver_dFdQ, TwoPointFlux_Selector public CentralRiemannSolver_acc, RiemannSolver_Selector + public KennedyGruber_TwoPointFlux, TwoPointFluxFCN, EntropyConserving_TwoPointFlux, Pirozzoli_TwoPointFlux, Ducros_TwoPointFlux, Morinishi_TwoPointFlux, StandardDG_TwoPointFlux, chandrasekar_twopointflux abstract interface !subroutine RiemannSolverFCN(Nx, Ny, QLeft, QRight, nHat, t1, t2, flux) From dc8c4fec6d23c295273f673a6d49ed1b81a68ca6 Mon Sep 17 00:00:00 2001 From: zalbanob <171132623+zalbanob@users.noreply.github.com> Date: Wed, 14 May 2025 11:58:59 +0200 Subject: [PATCH 2/5] Added HOPR Optimization --- Solver/src/libs/mesh/HexMesh.f90 | 97 ++++++----- Solver/src/libs/mesh/MappedGeometry.f90 | 166 +++++------------- Solver/src/libs/mesh/MeshPartitioning.f90 | 182 ++++++++------------ Solver/src/libs/mesh/Read_HDF5Mesh_HOPR.f90 | 38 ++-- 4 files changed, 179 insertions(+), 304 deletions(-) diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index 5bd1425a4..ff7650ff6 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -4106,32 +4106,34 @@ subroutine HexMesh_ComputeWallDistances(self,facesList,elementList) endif if( .not. self% IBM% active ) 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) + 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 - self % elements(eID) % geom % dWall(i, j, k) = sqrt(minimumDistance) end do + !$acc end parallel loop + !$acc update self(self % elements(eID) % geom % dWall(:, :, k)) async(k) end do - !$acc end parallel loop - !$acc update self(self % elements(eID) % geom % dWall(:, :, k)) async(k) - end do - !$acc wait - !$acc end data - !$acc end data - !$acc end data - !$acc end data + !$acc wait + !$acc end data + !$acc end data + !$acc end data + !$acc end data + endif end if end associate end do @@ -4157,32 +4159,33 @@ subroutine HexMesh_ComputeWallDistances(self,facesList,elementList) endif if( .not. self% IBM% active ) 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) + 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 - - self % faces(eID) % geom % dWall(i, j) = sqrt(minimumDistance) + !$acc end parallel loop + !$acc update self(self % faces(eID) % geom % dWall(:j)) async(j) end do - !$acc end parallel loop - !$acc update self(self % faces(eID) % geom % dWall(:j)) async(j) - end do - !$acc wait - !$acc end data - !$acc end data - !$acc end data - !$acc end data - 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 diff --git a/Solver/src/libs/mesh/MappedGeometry.f90 b/Solver/src/libs/mesh/MappedGeometry.f90 index 9b8b2c710..08efbebf1 100644 --- a/Solver/src/libs/mesh/MappedGeometry.f90 +++ b/Solver/src/libs/mesh/MappedGeometry.f90 @@ -193,14 +193,16 @@ subroutine computeMetricTermsConservativeForm(self, spAxi, spAeta, spAzeta, mapp integer :: i, j, k, m, n, l real(kind=RP) :: grad_x(NDIM,NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) real(kind=RP) :: xCGL(NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) - real(kind=RP) :: auxgrad(NDIM,NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) - real(kind=RP) :: coordsProduct(NDIM,0:self % Nx,0:self % Ny,0:self % Nz) - real(kind=RP) :: Jai(NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) - real(kind=RP) :: Ja1CGL(NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) - real(kind=RP) :: Ja2CGL(NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) - real(kind=RP) :: Ja3CGL(NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) + real(kind=RP) :: auxgrad(4,NDIM,NDIM,0:self % Nx, 0:self % Ny, 0:self % Nz) + real(kind=RP) :: coordsProduct(4,NDIM,0:self % Nx,0:self % Ny,0:self % Nz) + real(kind=RP) :: Ja1CGL(4,0:self % Nx, 0:self % Ny, 0:self % Nz) + real(kind=RP) :: Ja2CGL(4,0:self % Nx, 0:self % Ny, 0:self % Nz) + real(kind=RP) :: Ja3CGL(4,0:self % Nx, 0:self % Ny, 0:self % Nz) real(kind=RP) :: JacobianCGL(0:self % Nx, 0:self % Ny, 0:self % Nz) real(kind=RP) :: x(3) + + real(kind=RP) :: jGradXi_local(NDIM), jGradEta_local(NDIM), jGradZeta_local(NDIM), jacobian_local + real(kind=RP) :: auxgrad_l1(NDIM), auxgrad_l2(NDIM), auxgrad_l3(NDIM) ! ! Compute the mapping gradient in Chebyshev-Gauss-Lobatto points ! -------------------------------------------------------------- @@ -217,125 +219,27 @@ subroutine computeMetricTermsConservativeForm(self, spAxi, spAeta, spAzeta, mapp ! Compute coordinates combination ! ------------------------------- do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - coordsProduct(:,i,j,k) = xCGL(3,i,j,k) * grad_x(2,:,i,j,k) & - - xCGL(2,i,j,k) * grad_x(3,:,i,j,k) + coordsProduct(1,:,i,j,k) = xCGL(3,i,j,k) * grad_x(2,:,i,j,k) - xCGL(2,i,j,k) * grad_x(3,:,i,j,k) + coordsProduct(2,:,i,j,k) = xCGL(1,i,j,k) * grad_x(3,:,i,j,k) - xCGL(3,i,j,k) * grad_x(1,:,i,j,k) + coordsProduct(3,:,i,j,k) = xCGL(2,i,j,k) * grad_x(1,:,i,j,k) - xCGL(1,i,j,k) * grad_x(2,:,i,j,k) end do ; end do ; end do -! -! Compute its gradient -! -------------------- - auxgrad = 0.0_RP - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - do l = 0, self % Nx - auxgrad(:,1,i,j,k) = auxgrad(:,1,i,j,k) + coordsProduct(:,l,j,k) * spAxi % DCGL(i,l) - end do - do l = 0, self % Ny - auxgrad(:,2,i,j,k) = auxgrad(:,2,i,j,k) + coordsProduct(:,i,l,k) * spAeta % DCGL(j,l) - end do - - do l = 0, self % Nz - auxgrad(:,3,i,j,k) = auxgrad(:,3,i,j,k) + coordsProduct(:,i,j,l) * spAzeta % DCGL(k,l) - end do - end do ; end do ; end do -! -! Compute the curl -! ---------------- - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - Jai(1,i,j,k) = auxgrad(3,2,i,j,k) - auxgrad(2,3,i,j,k) - Jai(2,i,j,k) = auxgrad(1,3,i,j,k) - auxgrad(3,1,i,j,k) - Jai(3,i,j,k) = auxgrad(2,1,i,j,k) - auxgrad(1,2,i,j,k) - end do ; end do ; end do -! -! Assign to the first coordinate of each metrics -! ---------------------------------------------- - Ja1CGL(1,:,:,:) = -0.5_RP * Jai(1,:,:,:) - Ja2CGL(1,:,:,:) = -0.5_RP * Jai(2,:,:,:) - Ja3CGL(1,:,:,:) = -0.5_RP * Jai(3,:,:,:) -! -! ***************************************** -! Compute the y-coordinates of the mappings -! ***************************************** -! -! Compute coordinates combination -! ------------------------------- - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - coordsProduct(:,i,j,k) = xCGL(1,i,j,k) * grad_x(3,:,i,j,k) & - - xCGL(3,i,j,k) * grad_x(1,:,i,j,k) - end do ; end do ; end do -! -! Compute its gradient -! -------------------- auxgrad = 0.0_RP - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx + do k = 0, self % Nz; do j = 0, self % Ny; do i = 0, self % Nx do l = 0, self % Nx - auxgrad(:,1,i,j,k) = auxgrad(:,1,i,j,k) + coordsProduct(:,l,j,k) * spAxi % DCGL(i,l) + auxgrad(1:4,:,1,i,j,k) = auxgrad(1:4,:,1,i,j,k) + coordsProduct(1:4,:,l,j,k) * spAxi % DCGL(i,l) end do - do l = 0, self % Ny - auxgrad(:,2,i,j,k) = auxgrad(:,2,i,j,k) + coordsProduct(:,i,l,k) * spAeta % DCGL(j,l) + auxgrad(1:4,:,2,i,j,k) = auxgrad(1:4,:,2,i,j,k) + coordsProduct(1:4,:,i,l,k) * spAeta % DCGL(j,l) end do - do l = 0, self % Nz - auxgrad(:,3,i,j,k) = auxgrad(:,3,i,j,k) + coordsProduct(:,i,j,l) * spAzeta % DCGL(k,l) - end do - end do ; end do ; end do -! -! Compute the curl -! ---------------- - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - Jai(1,i,j,k) = auxgrad(3,2,i,j,k) - auxgrad(2,3,i,j,k) - Jai(2,i,j,k) = auxgrad(1,3,i,j,k) - auxgrad(3,1,i,j,k) - Jai(3,i,j,k) = auxgrad(2,1,i,j,k) - auxgrad(1,2,i,j,k) - end do ; end do ; end do -! -! Assign to the second coordinate of each metrics -! ----------------------------------------------- - Ja1CGL(2,:,:,:) = -0.5_RP*Jai(1,:,:,:) - Ja2CGL(2,:,:,:) = -0.5_RP*Jai(2,:,:,:) - Ja3CGL(2,:,:,:) = -0.5_RP*Jai(3,:,:,:) -! -! ***************************************** -! Compute the z-coordinates of the mappings -! ***************************************** -! -! Compute coordinates combination -! ------------------------------- - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - coordsProduct(:,i,j,k) = xCGL(2,i,j,k) * grad_x(1,:,i,j,k) & - - xCGL(1,i,j,k) * grad_x(2,:,i,j,k) - end do ; end do ; end do -! -! Compute its gradient -! -------------------- - auxgrad = 0.0_RP - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - do l = 0, self % Nx - auxgrad(:,1,i,j,k) = auxgrad(:,1,i,j,k) + coordsProduct(:,l,j,k) * spAxi % DCGL(i,l) + auxgrad(1:4,:,3,i,j,k) = auxgrad(1:4,:,3,i,j,k) + coordsProduct(1:4,:,i,j,l) * spAzeta % DCGL(k,l) end do + Ja1CGL(1:4,i,j,k) = -0.5_RP * (auxgrad(1:4,3,2,i,j,k) - auxgrad(1:4,2,3,i,j,k)) + Ja2CGL(1:4,i,j,k) = -0.5_RP * (auxgrad(1:4,1,3,i,j,k) - auxgrad(1:4,3,1,i,j,k)) + Ja3CGL(1:4,i,j,k) = -0.5_RP * (auxgrad(1:4,2,1,i,j,k) - auxgrad(1:4,1,2,i,j,k)) + end do ; end do; end do - do l = 0, self % Ny - auxgrad(:,2,i,j,k) = auxgrad(:,2,i,j,k) + coordsProduct(:,i,l,k) * spAeta % DCGL(j,l) - end do - - do l = 0, self % Nz - auxgrad(:,3,i,j,k) = auxgrad(:,3,i,j,k) + coordsProduct(:,i,j,l) * spAzeta % DCGL(k,l) - end do - end do ; end do ; end do -! -! Compute the curl -! ---------------- - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - Jai(1,i,j,k) = auxgrad(3,2,i,j,k) - auxgrad(2,3,i,j,k) - Jai(2,i,j,k) = auxgrad(1,3,i,j,k) - auxgrad(3,1,i,j,k) - Jai(3,i,j,k) = auxgrad(2,1,i,j,k) - auxgrad(1,2,i,j,k) - end do ; end do ; end do -! -! Assign to the third coordinate of each metrics -! ---------------------------------------------- - Ja1CGL(3,:,:,:) = -0.5_RP * Jai(1,:,:,:) - Ja2CGL(3,:,:,:) = -0.5_RP * Jai(2,:,:,:) - Ja3CGL(3,:,:,:) = -0.5_RP * Jai(3,:,:,:) -! ! ******************** ! Compute the Jacobian ! ******************** @@ -350,37 +254,45 @@ subroutine computeMetricTermsConservativeForm(self, spAxi, spAeta, spAzeta, mapp ! Return to Gauss points ! ********************** ! - self % jGradXi = 0.0_RP - self % jGradEta = 0.0_RP + self % jGradXi = 0.0_RP + self % jGradEta = 0.0_RP self % jGradZeta = 0.0_RP - self % jacobian = 0.0_RP + self % jacobian = 0.0_RP do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx + jGradXi_local = 0 + jGradEta_local = 0 + jGradZeta_local = 0 + jacobian_local = 0 do n = 0, self % Nz ; do m = 0, self % Ny ; do l = 0, self % Nx - self % jGradXi(:,i,j,k) = self % jGradXi(:,i,j,k) + Ja1CGL(:,l,m,n) & - * spAxi % TCheb2Gauss(i,l) & + jGradXi_local(:) = jGradXi_local(:) + Ja1CGL(1:3,l,m,n) & + * spAxi % TCheb2Gauss(i,l) & * spAeta % TCheb2Gauss(j,m) & * spAzeta % TCheb2Gauss(k,n) - self % jGradEta(:,i,j,k) = self % jGradEta(:,i,j,k) + Ja2CGL(:,l,m,n) & - * spAxi % TCheb2Gauss(i,l) & + jGradEta_local(:) = jGradEta_local(:) + Ja2CGL(1:3,l,m,n) & + * spAxi % TCheb2Gauss(i,l) & * spAeta % TCheb2Gauss(j,m) & * spAzeta % TCheb2Gauss(k,n) - self % jGradZeta(:,i,j,k) = self % jGradZeta(:,i,j,k) + Ja3CGL(:,l,m,n) & + jGradZeta_local(:) = jGradZeta_local(:) + Ja3CGL(1:3,l,m,n) & * spAxi % TCheb2Gauss(i,l) & * spAeta % TCheb2Gauss(j,m) & * spAzeta % TCheb2Gauss(k,n) - self % jacobian(i,j,k) = self % jacobian(i,j,k) + JacobianCGL(l,m,n) & + jacobian_local = jacobian_local + JacobianCGL(l,m,n) & * spAxi % TCheb2Gauss(i,l) & * spAeta % TCheb2Gauss(j,m) & * spAzeta % TCheb2Gauss(k,n) end do ; end do ; end do + + self % jGradXi(:,i,j,k) = jGradXi_local(:) + self % jGradEta(:,i,j,k) = jGradEta_local(:) + self % jGradZeta(:,i,j,k) = jGradZeta_local(:) + self % jacobian(i,j,k) = jacobian_local + self % invJacobian(i,j,k) = 1._RP / jacobian_local end do ; end do ; end do - do k = 0, self % Nz ; do j = 0, self % Ny ; do i = 0, self % Nx - self % invJacobian(i,j,k) = 1._RP / self % jacobian(i,j,k) - end do ; end do ; end do + end subroutine computeMetricTermsConservativeForm ! !/////////////////////////////////////////////////////////////////////// diff --git a/Solver/src/libs/mesh/MeshPartitioning.f90 b/Solver/src/libs/mesh/MeshPartitioning.f90 index f02e974bf..f54bcef50 100644 --- a/Solver/src/libs/mesh/MeshPartitioning.f90 +++ b/Solver/src/libs/mesh/MeshPartitioning.f90 @@ -93,120 +93,78 @@ end subroutine GetElementsDomain ! !//////////////////////////////////////////////////////////////////////// ! - subroutine GetNodesPartition(mesh, no_of_domains, elementsDomain, partitions) - use Utilities, only: Qsort - implicit none - type(HexMesh), intent(in) :: mesh - integer, intent(in) :: no_of_domains - integer, intent(in) :: elementsDomain(mesh % no_of_elements) - type(PartitionedMesh_t), intent(inout) :: partitions(no_of_domains) -! -! --------------- -! Local Variables -! --------------- -! - integer :: nvertex - integer :: i - integer :: j - integer :: k - integer :: ipoint - integer :: jpoint - integer :: idomain - integer :: npoints - integer :: ielem - logical :: isnewpoint - logical :: meshIsHOPR - integer, allocatable :: points(:) - integer, allocatable :: HOPRpoints(:) - - nvertex = 8 - - meshIsHOPR = allocated (mesh % HOPRnodeIDs) - - do idomain=1,no_of_domains -! -! Get the number of elements for the partition -! -------------------------------------------- - partitions(idomain)%no_of_elements = count(elementsDomain == idomain) - allocate(partitions(idomain)%elementIDs(partitions(idomain)%no_of_elements)) -! -! This will store the partition nodes (allocated as 8 * no_of_elements) -! --------------------------------------------------------------------- - allocate(points(nvertex*partitions(idomain)%no_of_elements)) - points = 0 - if (meshIsHOPR) then - allocate(HOPRpoints(nvertex*partitions(idomain)%no_of_elements)) - HOPRpoints = 0 - end if -! -! **************************************** -! Gather each partition nodes and elements -! **************************************** -! - k = 0 - npoints = 0 - do ielem=1,mesh % no_of_elements - if (elementsDomain(ielem) == idomain) then -! -! Append a new element -! -------------------- - k = k + 1 - partitions(idomain)%elementIDs(k) = ielem -! -! Append its nodes -! ---------------- - do j=1,nvertex -! -! Get the node ID -! --------------- - jpoint = mesh % elements(ielem) % nodeIDs(j) -! -! Check if it is already stored -! ----------------------------- - isnewpoint = .true. - do i=1,npoints - ipoint = points(i) - if (jpoint == ipoint) then - isnewpoint = .false. - exit - end if - end do -! -! Store the node -! -------------- - if (isnewpoint) then - npoints = npoints + 1 - points(npoints) = jpoint - if (meshIsHOPR) HOPRpoints(npoints) = mesh % HOPRnodeIDs(jpoint) - end if - end do - end if + subroutine GetNodesPartition(mesh, no_of_domains, elementsDomain, partitions) + use Utilities, only: Qsort + implicit none + + type(HexMesh), intent(in) :: mesh + integer, intent(in) :: no_of_domains + integer, intent(in) :: elementsDomain(mesh%no_of_elements) + type(PartitionedMesh_t), intent(inout) :: partitions(no_of_domains) + + integer :: nvertex, i, j, k, ipoint, jpoint + integer :: idomain, npoints, ielem + logical :: isnewpoint, meshIsHOPR + integer, allocatable :: points(:), HOPRpoints(:) + logical, allocatable :: mask(:) ! Mask array for node uniqueness + integer :: max_nodes + + nvertex = 8 + meshIsHOPR = allocated(mesh%HOPRnodeIDs) + max_nodes = 0 + do ielem = 1, mesh%no_of_elements + max_nodes = max(maxval(mesh%elements(ielem)%nodeIDs(:)), max_nodes) end do -! -! Put the nodeIDs into the partitions structure -! --------------------------------------------- - allocate(partitions(idomain)%nodeIDs(npoints)) - partitions(idomain)%nodeIDs(:) = points(1:npoints) - - if (meshIsHOPR) then - allocate(partitions(idomain)%HOPRnodeIDs(npoints)) - partitions(idomain)%HOPRnodeIDs(:) = HOPRpoints(1:npoints) - end if -! -! Sort the nodeIDs to read the mesh file accordingly (only needed for SpecMesh) -! -------------------------------------------------- - if (.not. meshIsHOPR) call Qsort(partitions(idomain)%nodeIDs) - - partitions(idomain)%no_of_nodes = npoints -! -! **** -! Free -! **** -! - deallocate(points) - safedeallocate(HOPRpoints) - end do + allocate(mask(max_nodes)) ! One-time allocation + + do idomain = 1, no_of_domains + partitions(idomain)%no_of_elements = count(elementsDomain == idomain) + allocate(partitions(idomain)%elementIDs(partitions(idomain)%no_of_elements)) + allocate(points(nvertex * partitions(idomain)%no_of_elements)) + points = 0 + + if (meshIsHOPR) then + allocate(HOPRpoints(nvertex * partitions(idomain)%no_of_elements)) + HOPRpoints = 0 + end if + + mask = .false. ! Reset mask for new domain + k = 0 + npoints = 0 + do ielem = 1, mesh%no_of_elements + if (elementsDomain(ielem) == idomain) then + k = k + 1 + partitions(idomain)%elementIDs(k) = ielem + do j = 1, nvertex + jpoint = mesh%elements(ielem)%nodeIDs(j) + if (.not. mask(jpoint)) then + npoints = npoints + 1 + points(npoints) = jpoint + mask(jpoint) = .true. + if (meshIsHOPR) HOPRpoints(npoints) = mesh%HOPRnodeIDs(jpoint) + end if + end do + end if + end do + + allocate(partitions(idomain)%nodeIDs(npoints)) + partitions(idomain)%nodeIDs = points(1:npoints) + + if (meshIsHOPR) then + allocate(partitions(idomain)%HOPRnodeIDs(npoints)) + partitions(idomain)%HOPRnodeIDs = HOPRpoints(1:npoints) + end if + + if (.not. meshIsHOPR) call Qsort(partitions(idomain)%nodeIDs) + + partitions(idomain)%no_of_nodes = npoints + + deallocate(points) + if (allocated(HOPRpoints)) deallocate(HOPRpoints) + end do + + deallocate(mask) end subroutine GetNodesPartition ! !//////////////////////////////////////////////////////////////////////// diff --git a/Solver/src/libs/mesh/Read_HDF5Mesh_HOPR.f90 b/Solver/src/libs/mesh/Read_HDF5Mesh_HOPR.f90 index 26f75d586..f4ab29784 100644 --- a/Solver/src/libs/mesh/Read_HDF5Mesh_HOPR.f90 +++ b/Solver/src/libs/mesh/Read_HDF5Mesh_HOPR.f90 @@ -143,6 +143,7 @@ subroutine ConstructMesh_FromHDF5File_( self, fileName, nodes, Nx, Ny, Nz, MeshI integer, allocatable :: HOPRNodeMap(:) ! Map from the global node index of HORSES3D to the global node index of HOPR real(kind=RP), allocatable :: TempNodes(:,:) ! Nodes read from file to be exported to self % nodes logical :: CurveCondition + integer, allocatable :: NodeToHOPRMap(:) !--------------------------------------------------------------- ! @@ -255,10 +256,8 @@ subroutine ConstructMesh_FromHDF5File_( self, fileName, nodes, Nx, Ny, Nz, MeshI ! Read nodeIDs and add them to the self%nodes array DO k = 1, NODES_PER_ELEMENT HOPRNodeID = ElemInfo(ELEM_FirstNodeInd,l) + HCornerMap(k) - corners(:,k) = NodeCoords(:,HOPRNodeID) / Lref - - call AddToNodeMap (TempNodes , HOPRNodeMap, corners(:,k), GlobalNodeIDs(HOPRNodeID), nodeIDs(k)) + call AddToNodeMap(TempNodes , HOPRNodeMap, NodeToHOPRMap, corners(:,k), GlobalNodeIDs(HOPRNodeID), nodeIDs(k)) END DO do k = 1, FACES_PER_ELEMENT @@ -1241,9 +1240,10 @@ subroutine InitNodeMap (TempNodes , HOPRNodeMap, nUniqueNodes) allocate (TempNodes(3,nUniqueNodes)) allocate (HOPRNodeMap(nUniqueNodes)) - - TempNodes = 0._RP - HOPRNodeMap = 0 + allocate(NodeToHOPRMap(nUniqueNodes)) + NodeToHOPRMap = -1 + TempNodes = 0._RP + HOPRNodeMap = 0 end subroutine InitNodeMap @@ -1253,11 +1253,13 @@ end subroutine InitNodeMap ! ---------------------------------------------------------------- ! Add a new entry to the node map ! ---------------------------------------------------------------- - subroutine AddToNodeMap (TempNodes , HOPRNodeMap, newnode, HOPRGlobalID,nodeID) + + subroutine AddToNodeMap(TempNodes , HOPRNodeMap, NodeToHOPRMap, newnode, HOPRGlobalID,nodeID) implicit none !-------------------------------------------- real(kind=RP), intent(inout) :: TempNodes(:,:) integer , intent(inout) :: HOPRNodeMap(:) + integer , intent(inout) :: NodeToHOPRMap(:) real(kind=RP), intent(in) :: newnode(3) integer , intent(in) :: HOPRGlobalID integer , intent(out) :: nodeID ! Node ID in HORSES3D! @@ -1265,21 +1267,20 @@ subroutine AddToNodeMap (TempNodes , HOPRNodeMap, newnode, HOPRGlobalID,nodeID) integer :: i ! Counter !-------------------------------------------- - do i = 1, idx - if (HOPRGlobalID == HOPRNodeMap(i)) then - nodeID = i - return - end if - end do - + + if (NodeToHOPRMap(HOPRGlobalID) .ne. -1 ) then + nodeID = NodeToHOPRMap(HOPRGlobalID) + return + end if + idx = idx + 1 nodeID = idx - - TempNodes(:,idx) = newnode - HOPRNodeMap(idx) = HOPRGlobalID - + TempNodes(:,idx) = newnode + HOPRNodeMap(idx) = HOPRGlobalID + NodeToHOPRMap(HOPRGlobalID) = idx end subroutine AddToNodeMap + ! !/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! @@ -1318,6 +1319,7 @@ subroutine FinishNodeMap (TempNodes , HOPRNodeMap, nodes, HOPRnodeIDs) deallocate (TempNodes) deallocate (HOPRNodeMap) + deallocate(NodeToHOPRMap) end subroutine FinishNodeMap From 3f62fab4f6827d477b168331dc84849b6dfd210a Mon Sep 17 00:00:00 2001 From: Oscar Marino Date: Thu, 22 May 2025 12:41:03 +0200 Subject: [PATCH 3/5] Change get_rho in LESModel to be pure function such that compiles in cpu --- Solver/src/libs/physics/common/LESModels.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Solver/src/libs/physics/common/LESModels.f90 b/Solver/src/libs/physics/common/LESModels.f90 index 62caa159e..3416bbd28 100644 --- a/Solver/src/libs/physics/common/LESModels.f90 +++ b/Solver/src/libs/physics/common/LESModels.f90 @@ -608,7 +608,7 @@ subroutine Vreman_Describe(self) end subroutine Vreman_Describe - function get_rho(Q, dimensionless_) result(rho) + pure function get_rho(Q, dimensionless_) result(rho) !$acc routine seq implicit none real(kind=RP), intent(in) :: Q(:) From afbdef18d4c84970b7c31729889e7b619e17e04e Mon Sep 17 00:00:00 2001 From: Oscar Marino Date: Thu, 22 May 2025 12:41:38 +0200 Subject: [PATCH 4/5] Update AL new implementation. --- Solver/src/libs/sources/ActuatorLine.f90 | 158 ++++++------------ .../libs/timeintegrator/TimeIntegrator.f90 | 3 + 2 files changed, 56 insertions(+), 105 deletions(-) diff --git a/Solver/src/libs/sources/ActuatorLine.f90 b/Solver/src/libs/sources/ActuatorLine.f90 index d92e36167..b03440912 100644 --- a/Solver/src/libs/sources/ActuatorLine.f90 +++ b/Solver/src/libs/sources/ActuatorLine.f90 @@ -81,7 +81,6 @@ module ActuatorLine logical :: save_average = .false. logical :: save_instant = .false. logical :: verbose = .false. - logical :: averageSubElement = .true. character(len=LINE_LENGTH) :: file_name integer :: number_iterations integer :: save_iterations @@ -110,7 +109,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) use fluiddata use MPI_Process_Info implicit none - class(farm_t) , intent(inout) :: self + type(farm_t) , intent(inout) :: self TYPE(FTValueDictionary), intent(in) :: controlVariables real(kind=RP), intent(in) :: t0 type(HexMesh), intent(in) :: mesh @@ -119,11 +118,11 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) ! --------------- ! integer :: i, j, k, ii, fid, n_aoa, n_airfoil - CHARACTER(LEN=LINE_LENGTH) :: arg, char1 - CHARACTER(LEN=LINE_LENGTH) :: solution_file - CHARACTER(LEN=5) :: file_id + character(LEN=LINE_LENGTH) :: arg, char1 + character(LEN=LINE_LENGTH) :: solution_file + character(len=LINE_LENGTH) :: restart_name, restart_operations_name + character(LEN=5) :: file_id real(kind=RP), dimension(:), allocatable :: initial_azimutal - character(len=STRING_CONSTANT_LENGTH) :: restart_name, restart_operations_name logical :: fileExists integer :: nelem, eID, eIndex real(kind=RP) :: tolerance, r_square @@ -140,7 +139,6 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) self % save_instant = controlVariables % getValueOrDefault("actuator save instant", .false.) self % save_iterations = controlVariables % getValueOrDefault("actuator save iteration", 1) self % verbose = controlVariables % getValueOrDefault("actuator verbose", .false.) - self % averageSubElement = controlVariables % getValueOrDefault("actuator average subelement", .true.) self % tolerance_factor = controlVariables % getValueOrDefault("actuator tolerance", 0.2_RP) restart_name = controlVariables % stringValueForKey( restartFileNameKey, requestedLength = STRING_CONSTANT_LENGTH ) @@ -299,7 +297,6 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) end select write(STD_OUT,'(30X,A,A28,F10.3,F10.3)') "->", 'Tip correction constants: ', self%turbine_t(1)%blade_t(1)%tip_c1, self%turbine_t(1)%blade_t(1)%tip_c2 write(STD_OUT,'(30X,A,A28,L1)') "->", "Projection formulation: ", self % calculate_with_projection - if (.not. self%calculate_with_projection) write(STD_OUT,'(30X,A,A28,L1)') "->", "Average sub-Element: ", self % averageSubElement write(STD_OUT,'(30X,A,A28,L1)') "->", "Save blade average values: ", self % save_average if (fileExists) write(STD_OUT,'(30X,A)') 'Using restaring operations of turbines' end if @@ -526,7 +523,7 @@ end subroutine ConstructFarm !/////////////////////////////////////////////////////////////////////////////////////// subroutine DestructFarm(self) implicit none - class(Farm_t), intent(inout) :: self + type(Farm_t), intent(inout) :: self ! integer :: i, j, k @@ -573,7 +570,7 @@ subroutine UpdateFarm(self,time, mesh) use MPI_Process_Info implicit none - class(Farm_t), intent(inout) :: self + type(Farm_t), intent(inout) :: self real(kind=RP), intent(in) :: time type(HexMesh), intent(in) :: mesh @@ -660,8 +657,8 @@ subroutine UpdateFarm(self,time, mesh) x = [self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,1),self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,2),self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,3)] call FindActuatorPointElement(mesh, x, eID, xi, found) if (found) then - ! averaged state values of the cell - Qtemp = element_averageQ(mesh,eID, xi, self % averageSubElement) + ! interpolate state values in the element + Qtemp = interpolateQ(mesh,eID, xi) delta_temp = (mesh % elements(eID) % geom % Volume / product(mesh % elements(eID) % Nxyz + 1)) ** (1.0_RP / 3.0_RP) else Qtemp = 0.0_RP @@ -711,7 +708,7 @@ subroutine ForcesFarm(self, mesh, time) use fluiddata implicit none - class(Farm_t) , intent(inout) :: self + type(Farm_t) , intent(inout) :: self type(HexMesh), intent(in) :: mesh real(kind=RP),intent(in) :: time @@ -750,9 +747,6 @@ subroutine ForcesFarm(self, mesh, time) interp = GaussianInterpolation(self%epsilon_type, mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), & self%turbine_t(kk)%blade_t(jj)%chord(ii), self%gauss_epsil,self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii)) call FarmGetLocalForces(self, ii, jj, kk, mesh%elements(eID)%storage%Q(:,i,j,k), interp, local_angle, local_velocity, local_Re, local_thrust, local_rotor_force) - ! todo - ! set scalar values as private - ! reduce accumulated arrays ! minus account action-reaction effect, is the force on the fliud actuator_source(1) = actuator_source(1) - local_thrust @@ -792,7 +786,6 @@ subroutine ForcesFarm(self, mesh, time) else ! no projection - ! print *, "start loop " !$acc parallel loop gang present(self,mesh) !$omp do schedule(runtime) private(i,j,k,ii,jj,kk,actuator_source,eID,interp) do eIndex = 1, size(elementsActuated) @@ -858,7 +851,7 @@ subroutine WriteFarmForces(self,time,iter,last) use MPI_Process_Info implicit none - class(Farm_t), intent(inout) :: self + type(Farm_t), intent(inout) :: self real(kind=RP),intent(in) :: time integer, intent(in) :: iter logical, optional :: last @@ -880,10 +873,25 @@ subroutine WriteFarmForces(self,time,iter,last) isLast = .false. end if + ! for last iteration save average values only, all calculations and comunications were done before + if (isLast) then + if ( .not. self % save_average ) return + if ( .not. MPI_Process % isRoot ) return + do kk=1, self%num_turbines + write(file_id, '(I3.3)') kk + write(arg , '(A,A,A,A)') trim(self%file_name), "_Actuator_Line_average_turb_", trim(file_id) , ".dat" + open( newunit = fID , file = trim(arg) , action = "write" , access = "append" , status = "old" ) + do ii = 1, self % turbine_t(kk) % num_blade_sections + write(fid,"(6(2X,ES24.16))") self%turbine_t(kk)%blade_t(1)%r_R(ii), self%turbine_t(kk)%average_conditions(ii,:) + end do + close(fid) + end do + return + end if + save_instant = self%save_instant .and. ( mod(iter,self % save_iterations) .eq. 0 ) t = time * Lref / refValues%V - if (self%calculate_with_projection) then ! this is necessary for Gaussian weighted sum @@ -947,8 +955,7 @@ subroutine WriteFarmForces(self,time,iter,last) if ( .not. MPI_Process % isRoot ) return ! save in memory the time step forces for each element blade and the whole blades - if (.not. isLast) call FarmUpdateBladeForces(self) - ! if (.not. self%calculate_with_projection .and. .not. isLast) call self % FarmUpdateBladeForces() + call FarmUpdateBladeForces(self) !write output torque thrust to file do kk=1, self%num_turbines @@ -967,16 +974,7 @@ subroutine WriteFarmForces(self,time,iter,last) write(fid,"(10(2X,ES24.16))") t, self%turbine_t(kk)%Cp, self%turbine_t(kk)%Ct close(fid) - if (self % save_average .and. isLast) then - write(arg , '(A,A,A,A)') trim(self%file_name), "_Actuator_Line_average_turb_", trim(file_id) , ".dat" - open( newunit = fID , file = trim(arg) , action = "write" , access = "append" , status = "old" ) - do ii = 1, self % turbine_t(kk) % num_blade_sections - write(fid,"(6(2X,ES24.16))") self%turbine_t(kk)%blade_t(1)%r_R(ii), self%turbine_t(kk)%average_conditions(ii,:) - end do - close(fid) - end if - - if (save_instant .and. .not. isLast) then + if (save_instant) then do jj = 1, self%turbine_t(kk)%num_blades write(arg , '(2A,I3.3,A,I10.10,3A)') trim(self%file_name), "_Actuator_Line_instant_",jj ,"_" ,iter, "_turb_", trim(file_id), ".dat" open ( newunit = fID , file = trim(arg) , status = "unknown" , action = "write" ) @@ -1006,7 +1004,7 @@ Subroutine FarmGetLocalForces(self, ii, jj, kk, Q, interp, local_angle, local_ve use VariableConversion, only: Temperature, SutherlandsLaw #endif implicit none - class(Farm_t) :: self + type(Farm_t) :: self integer, intent(in) :: ii, jj, kk real(kind=RP), dimension(NCONS), intent(in) :: Q real(kind=RP), intent(in) :: interp @@ -1120,7 +1118,7 @@ Subroutine FarmUpdateBladeForces(self) use fluiddata Implicit None - class(Farm_t), intent(inout) :: self + type(Farm_t), intent(inout) :: self !local variables integer :: ii, jj, kk real(kind=RP), dimension(:), allocatable :: aoa @@ -1329,32 +1327,8 @@ function InterpolateAirfoilData(x1,x2,y1,y2,new_x) InterpolateAirfoilData=a*new_x+b end function -function full_element_averageQ(mesh,eID) - use HexMeshClass - use PhysicsStorage - use NodalStorageClass - implicit none - - type(HexMesh), intent(in) :: mesh - integer, intent(in) :: eID - integer :: k, j, i - - integer :: total_points - real(kind=RP), dimension(NCONS) :: full_element_averageQ, Qsum - - - Qsum(:) = 0.0_RP - total_points = 0 - do k = 0, mesh%elements(eID) % Nxyz(3) ; do j = 0, mesh%elements(eID) % Nxyz(2) ; do i = 0, mesh%elements(eID) % Nxyz(1) - Qsum(:)=Qsum(:)+mesh%elements(eID) % Storage % Q(:,i,j,k) - total_points=total_points + 1 - end do ; end do ; end do - - full_element_averageQ(:) = Qsum(:) / real(total_points,RP) - -end function full_element_averageQ - -Function semi_element_averageQ(mesh,eID,xi) result(Qe) +! high order interpolation of Q +Function interpolateQ(mesh,eID,xi) result(Qe) use HexMeshClass use PhysicsStorage use NodalStorageClass @@ -1365,59 +1339,33 @@ Function semi_element_averageQ(mesh,eID,xi) result(Qe) real(kind=RP), dimension(NDIM), intent(in) :: xi real(kind=RP), dimension(NCONS) :: Qe - real(kind=RP), dimension(NCONS) :: Qsum - - integer :: k, j, i, direction, N, ind - integer, dimension(NDIM) :: firstNodeIndex - integer :: total_points - type(NodalStorage_t), pointer :: spAxi - - ! fist get the sub element nodes index - do direction = 1, NDIM - - N = mesh % elements(eID) % Nxyz(direction) - spAxi => NodalStorage(N) + integer :: k, j, i + integer, dimension(NDIM) :: Nxyz + type(NodalStorage_t), pointer :: spAxi, spAeta, spAzeta + real(kind=RP), allocatable :: lxi(:) , leta(:), lzeta(:) !interpolants - do ind = 0, N - firstNodeIndex(direction) = ind-1 - if (xi(direction) .le. spAxi%x(ind)) exit - end do + Nxyz = mesh % elements(eID) % Nxyz - if (firstNodeIndex(direction) .eq. -1) firstNodeIndex(direction) = 0 + spAxi => NodalStorage(Nxyz(1)) + spAeta => NodalStorage(Nxyz(2)) + spAzeta => NodalStorage(Nxyz(3)) - end do + allocate( lxi(0:Nxyz(1)), leta(0:Nxyz(2)), lzeta(0:Nxyz(3)) ) - nullify(spAxi) + lxi = spAxi % lj(xi(1)) + leta = spAeta % lj(xi(2)) + lzeta = spAzeta % lj(xi(3)) - ! now average on the sub element - Qsum(:) = 0.0_RP - total_points = 0 - do k = firstNodeIndex(IZ), firstNodeIndex(IZ)+1 ; do j = firstNodeIndex(IY), firstNodeIndex(IY)+1 ; do i = firstNodeIndex(IX),firstNodeIndex(IX)+1 - Qsum(:) = Qsum(:) + mesh % elements(eID) % Storage % Q(:,i,j,k) - total_points = total_points + 1 - end do ; end do ; end do + !$acc update self(mesh%elements(eID)%Storage%Q) + Qe = 0.0_RP + do k = 0, Nxyz(3) ; do j = 0, Nxyz(2) ; do i = 0, Nxyz(1) + Qe = Qe + mesh % elements(eID) % Storage % Q(:,i,j,k) * lxi(i) * leta(j) * lzeta(k) + end do ; end do ; end do - Qe(:) = Qsum(:) / real(total_points,RP) + deallocate(lxi,leta,lzeta) + nullify(spAxi,spAeta,spAzeta) -End Function semi_element_averageQ - -Function element_averageQ(mesh,eID,xi,averageSubElement) result(Qe) - use HexMeshClass - use PhysicsStorage - Implicit None - type(HexMesh), intent(in) :: mesh - integer, intent(in) :: eID - logical, intent(in) :: averageSubElement - real(kind=RP), dimension(NDIM), intent(in) :: xi - real(kind=RP), dimension(NCONS) :: Qe - - !$acc update self(mesh%elements(eID)%Storage%Q) - if (averageSubElement) then - Qe = semi_element_averageQ(mesh, eid, xi) - else - Qe = full_element_averageQ(mesh, eid) - end if -End Function element_averageQ +END Function interpolateQ #endif end module diff --git a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 index 446031a17..134494fd2 100644 --- a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 +++ b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 @@ -800,6 +800,9 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe #if defined(NAVIERSTOKES) && (!(SPALARTALMARAS)) call sem % fwh % writeToFile( force = .TRUE. ) call sponge % writeBaseFlow(sem % mesh, k, t, last=.true.) +#endif +#if defined(NAVIERSTOKES) + if(ActuatorLineFlag) call WriteFarmForces(farm, t, k, last=.true.) #endif end if From ba0f3c2a9aea339cfeb5a49d237b9e4cf81e60b6 Mon Sep 17 00:00:00 2001 From: Oscar Marino Date: Thu, 22 May 2025 12:42:43 +0200 Subject: [PATCH 5/5] delete prints --- Solver/src/libs/sources/ActuatorLine.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Solver/src/libs/sources/ActuatorLine.f90 b/Solver/src/libs/sources/ActuatorLine.f90 index b03440912..57bec8fb1 100644 --- a/Solver/src/libs/sources/ActuatorLine.f90 +++ b/Solver/src/libs/sources/ActuatorLine.f90 @@ -728,7 +728,7 @@ subroutine ForcesFarm(self, mesh, time) #endif if (.not. self % active) return - print *, "start ForcesFarm" + ! print *, "start ForcesFarm" Non_dimensional = POW2(refValues % V) * refValues % rho / Lref t = time * Lref / refValues % V @@ -839,7 +839,7 @@ subroutine ForcesFarm(self, mesh, time) ! print *, "s: ", mesh % elements(eID) % storage % S_NS endif - print *, "finish ForcesFarm" + ! print *, "finish ForcesFarm" end subroutine ForcesFarm !