Skip to content

Commit 5f9d7c5

Browse files
andrew-plattandrew-platt
committed
TurbSim: add OMP on CalcFourierCoeffs_IEC
- Very cludgy mess. Seems to work locally though.
1 parent 1184401 commit 5f9d7c5

File tree

1 file changed

+81
-21
lines changed

1 file changed

+81
-21
lines changed

‎modules/turbsim/src/TSsubs.f90‎

Lines changed: 81 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,10 @@
1919
!**********************************************************************************************************************************
2020
MODULE TSSubs
2121

22+
#ifdef _OPENMP
23+
use omp_lib
24+
#endif
25+
2226
USE ModifiedvKrm_mod
2327

2428
use TS_Profiles
@@ -39,14 +43,14 @@ MODULE TSSubs
3943
!! real array) of the simulated velocity (wind/water speed). It returns
4044
!! values FOR ONLY the velocity components that use the IEC method for
4145
!! computing spatial coherence; i.e., for i where SCMod(i) == CohMod_IEC
42-
SUBROUTINE CalcFourierCoeffs_IEC( p, U, PhaseAngles, S, V, TRH, ErrStat, ErrMsg )
46+
SUBROUTINE CalcFourierCoeffs_IEC( p, U, PhaseAngles, S, V, TRH_in, ErrStat, ErrMsg )
4347

4448
TYPE(TurbSim_ParameterType), INTENT(IN ) :: p !< TurbSim parameters
4549
REAL(ReKi), INTENT(IN) :: U (:) !< The steady u-component wind speeds for the grid (NPoints).
4650
REAL(ReKi), INTENT(IN) :: PhaseAngles (:,:,:) !< The array that holds the random phases [number of points, number of frequencies, number of wind components=3].
4751
REAL(ReKi), INTENT(IN) :: S (:,:,:) !< The turbulence PSD array (NumFreq,NPoints,3).
4852
REAL(ReKi), INTENT(INOUT) :: V (:,:,:) !< An array containing the summations of the rows of H (NumSteps,NPoints,3).
49-
REAL(ReKi), INTENT(INOUT) :: TRH (:) !< The transfer function matrix. just used as a work array
53+
REAL(ReKi), TARGET, INTENT(INOUT) :: TRH_in (:) !< The transfer function matrix. just used as a work array
5054
INTEGER(IntKi), INTENT(OUT) :: ErrStat
5155
CHARACTER(*), INTENT(OUT) :: ErrMsg
5256

@@ -61,10 +65,17 @@ SUBROUTINE CalcFourierCoeffs_IEC( p, U, PhaseAngles, S, V, TRH, ErrStat, ErrMsg
6165
INTEGER :: IFreq
6266
INTEGER :: Indx
6367
INTEGER :: IVec ! wind component, 1=u, 2=v, 3=w
68+
integer :: OMP_threads
69+
logical :: OMPerr ! error with OMP thread requiring abort
70+
real(ReKi), pointer :: TRH(:) ! temp transfer function matrix (work array)
71+
real(ReKi), allocatable, target :: TRH_lcl(:) ! make local copy
6472

6573
INTEGER(IntKi) :: ErrStat2
6674
CHARACTER(MaxMsgLen) :: ErrMsg2
6775

76+
#ifdef _OPENMP
77+
OMP_threads = omp_get_max_threads() ! limit number of threads to reasonable number (either set externally, or by hardware limit)
78+
#endif
6879

6980

7081
ErrStat = ErrID_None
@@ -115,7 +126,32 @@ SUBROUTINE CalcFourierCoeffs_IEC( p, U, PhaseAngles, S, V, TRH, ErrStat, ErrMsg
115126
! Calculate the coherence, Veers' H matrix (CSDs), and the fourier coefficients
116127
!---------------------------------------------------------------------------------
117128

129+
130+
OMPerr = .false.
131+
!$OMP PARALLEL DO &
132+
!$OMP NUM_THREADS(OMP_threads) &
133+
!$OMP PRIVATE(IFREQ, Indx, I, J, TRH_lcl, ErrStat2, ErrMsg2, TRH) &
134+
!$OMP SHARED(p,PhaseAngles,S,V,Dist,DistU,IVec,OMPerr,TRH_in,ErrStat,ErrMsg,AbortErrLev) &
135+
!$OMP DEFAULT(None)
118136
DO IFREQ = 1,p%grid%NumFreq
137+
! if an aborting error happened in a prior thread, skip
138+
if (OMPerr) cycle
139+
140+
#ifndef _OPENMP
141+
TRH => TRH_in
142+
#else
143+
if (.not. allocated(TRH_lcl)) then
144+
call AllocAry( TRH_lcl, p%grid%NPacked, 'TRH coherence array', ErrStat2, ErrMsg2 )
145+
!$OMP CRITICAL
146+
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'CalcFourierCoeffs')
147+
!$OMP END CRITICAL
148+
if (ErrStat2 >= AbortErrLev) OMPerr = .true.
149+
if (OMPerr) cycle
150+
TRH => TRH_lcl
151+
endif
152+
#endif
153+
154+
119155
! -----------------------------------------------
120156
! Create the coherence matrix for this frequency
121157
! -----------------------------------------------
@@ -149,14 +185,28 @@ SUBROUTINE CalcFourierCoeffs_IEC( p, U, PhaseAngles, S, V, TRH, ErrStat, ErrMsg
149185
! use H matrix to calculate coefficients
150186
! -----------------------------------------------
151187

152-
CALL Coh2H( p, IVec, IFreq, TRH, S, ErrStat2, ErrMsg2 )
153-
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'CalcFourierCoeffs_IEC')
154-
IF (ErrStat >= AbortErrLev) THEN
155-
CALL Cleanup()
156-
RETURN
157-
END IF
158-
CALL H2Coeffs( IVec, IFreq, TRH, PhaseAngles, V, p%grid%NPoints )
188+
CALL Coh2H( p, IVec, IFreq, TRH, S, ErrStat2, ErrMsg2 )
189+
!$OMP CRITICAL
190+
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'CalcFourierCoeffs_IEC')
191+
!$OMP END CRITICAL
192+
if (ErrStat2 >= AbortErrLev) then
193+
OMPerr = .true.
194+
else
195+
CALL H2Coeffs( IVec, IFreq, TRH, PhaseAngles, V, p%grid%NPoints )
196+
endif
197+
198+
#ifdef _OPENMP
199+
if (allocated(TRH_lcl)) deallocate(TRH_lcl)
200+
#endif
201+
202+
159203
END DO !IFreq
204+
!$OMP END PARALLEL DO
205+
206+
if (OMPerr) then
207+
call Cleanup()
208+
return
209+
endif
160210

161211
END DO !IVec
162212

@@ -753,6 +803,11 @@ SUBROUTINE Coh2H( p, IVec, IFreq, TRH, S, ErrStat, ErrMsg )
753803
integer :: Indx, J, I, NPts
754804

755805

806+
#ifdef _OPENMP
807+
call omp_set_max_active_levels(1) ! disallow the LAPACK_pptrf to use OMP parallelization (this kills performance)
808+
#endif
809+
810+
756811
! -------------------------------------------------------------
757812
! Calculate the Cholesky factorization for the coherence matrix
758813
! -------------------------------------------------------------
@@ -2025,28 +2080,33 @@ SUBROUTINE AddMeanAndRotate(p, V, U, HWindDir, VWindDir)
20252080
REAL(ReKi) :: v3(3) ! temporary 3-component array containing velocity
20262081
INTEGER(IntKi) :: ITime ! loop counter for time step
20272082
INTEGER(IntKi) :: IPoint ! loop counter for grid points
2028-
2029-
2030-
2031-
2083+
integer(IntKi) :: OMP_threads ! number of threads for OMP (if used)
2084+
2085+
#ifdef _OPENMP
2086+
OMP_threads = omp_get_max_threads() ! limit number of threads to reasonable number (either set externally, or by hardware limit)
2087+
#endif
2088+
20322089
!..............................................................................
2033-
! Add mean wind to u' components and rotate to inertial reference
2090+
! Add mean wind to u' components and rotate to inertial reference
20342091
! frame coordinate system
2035-
!..............................................................................
2092+
!..............................................................................
2093+
!$OMP PARALLEL DO &
2094+
!$OMP NUM_THREADS(OMP_threads) &
2095+
!$OMP PRIVATE(IPoint,ITime,v3) &
2096+
!$OMP SHARED(p, U, V, HWindDir, VWindDir) &
2097+
!$OMP DEFAULT(None)
20362098
DO IPoint=1,p%grid%Npoints
20372099
DO ITime=1,p%grid%NumSteps
20382100

2039-
! Add mean wind speed to the streamwise component and
2040-
! Rotate the wind to the X-Y-Z (inertial) reference frame coordinates:
2041-
2101+
! Add mean wind speed to the streamwise component and
2102+
! Rotate the wind to the X-Y-Z (inertial) reference frame coordinates:
20422103
v3 = V(ITime,IPoint,:)
20432104
CALL CalculateWindComponents( v3, U(IPoint), HWindDir(IPoint), VWindDir(IPoint), V(ITime,IPoint,:) )
20442105

20452106
ENDDO ! ITime
2046-
20472107
ENDDO ! IPoint
2048-
2049-
2108+
!$OMP END PARALLEL DO
2109+
20502110
END SUBROUTINE AddMeanAndRotate
20512111
!=======================================================================
20522112
SUBROUTINE TS_ValidateInput(P, ErrStat, ErrMsg)

0 commit comments

Comments
 (0)